# Dynamic correlations in the conserved Manna sandpile

Anirban Mukherjee  Punyabrata Pradhan  Department of Physics of Complex
Systems, S. N. Bose National Centre for Basic Sciences, Block-JD, Sector-III,
Salt Lake, Kolkata 700106, India.

###### Abstract

We study dynamic correlations for current and mass, as well as the associated
power spectra, in the one-dimensional conserved Manna sandpile. We show that,
in the thermodynamic limit, the variance of cumulative bond current up to time
$T$ grows subdiffusively as $T^{1/2-\mu}$ with the exponent $\mu\geq 0$
depending on the density regimes considered and, likewise, the power spectra
of current and mass at low frequency $f$ varies as $f^{1/2+\mu}$ and
$f^{-3/2+\mu}$, respectively; our theory predicts that, far from criticality,
$\mu=0$ and, near criticality, $\mu=(\beta+1)/2\nu_{\perp}z>0$ with $\beta$,
$\nu_{\perp}$ and $z$ being the order-parameter, correlation-length and
dynamic exponents, respectively. The anomalous suppression of fluctuations
near criticality signifies a “dynamic hyperuniformity”, characterized by a set
of fluctuation relations, in which current, mass and tagged-particle
displacement fluctuations are shown to have a precise quantitative
relationship with the density-dependent activity (or, its derivative). In
particular, the relation, ${\cal D}_{s}(\bar{\rho})=a(\bar{\rho})/\bar{\rho}$,
between the self-diffusion coefficient ${\cal D}_{s}(\bar{\rho})$, activity
$a(\bar{\rho})$ and density $\bar{\rho}$ explains a previous simulation
observation [Eur. Phys. J. B 72, 441 (2009)] that, near criticality, the self-
diffusion coefficient in the Manna sandpile has the same scaling behavior as
the activity.

##  I Introduction

Long-ranged temporal correlations are ubiquitous in nature Mandelbrot (1983).
They usually manifest as the so-called “$1/f$” or the flicker noise, having a
characteristic low-frequency power spectra with a power-law form $1/f^{\psi}$,
where $0<\psi<2$, over a wide range of frequency $f$. In fact, the $1/f$ noise
has been seen in a variety of seemingly unrelated systems, such as solar
flares Frontera and Fuligni (1979); Aschwanden and et al (2014), forest fires
Palmieri and Jensen (2020), electrical activities in brain Novikov et al.
(1997); Hesse and Gross (2014), stock-market fluctuations Tebaldi (2021),
water flow in rivers Mandelbrot (1983) and resistance fluctuations in
conductors Voss and Clarke (1976), among others; for details, see reviews
Jensen (1998); Bak (1997). However there is no general theory explaining the
relative abundance of $1/f$ noise in nature.

Bak, Tang, and Wisenfeld (BTW) proposed sandpiles as paradigmatic models of
”self-organized criticality” (SOC) in order to provide a generic mechanism for
long-ranged correlations in natural systems, and $1/f$ noise in particular.
Bak et al. (1987, 1988). Sandpiles are spatially extended and threshold-
activated systems, in which dynamical activities spread through cascades of
toppling events (initialed when a local threshold is crossed), resulting in
“avalanche”-like dynamical activities and long-ranged correlations in the
systems. They were envisaged as model-systems driven by slow addition of
“energy”, or grains, with local energy conservation in the bulk and
dissipation at the boundary. Due to an intriguing interplay between drive and
dissipation, the system evolves, apparently without fine-tuning of any
parameters, towards a nonequilibrium steady state characterized through
avalanches at all scales, i.e., a scale-invariant critical state with power-
law distributions concerning various observables. Later, several variants,
known as the conserved or “fixed-energy” sandpiles, were proposed, where there
is no dissipation, but the total number of grains (or, total energy) remains
conserved, thus allowing the critical state to be reached by tuning the global
density.

There have been numerous studies of sandpiles, and significant progress has
been made in characterizing the static and dynamic properties of both critical
and off-critical states of the systems. However, apart from some exact Dhar
and Ramaswamy (1989); Dhar (1990); Priezzhev et al. (1996); Majumdar and Dhar
(1992) and mathematically rigorous Fey et al. (2010); Podder and Rolla (2021);
Cipriani et al. (2017) results, the majority of the studies have been carried
out using simulations Dickman et al. (1998, 2001); Vespignani et al. (2000)
and phenomenological field-theoretical descriptions Vespignani et al. (1998);
Le Doussal and Wiese (2015); Ramasco et al. (2004). This is primarily due to
the fact that the steady-state measure of a driven interacting-particle system
such as sandpiles is in most cases a-priori unknown and, as a result, analytic
calculations, beginning with a microscopic dynamical description, prove to be
quite challenging Dickman (2002); Le Doussal and Wiese (2015); da Cunha et al.
(2009). Perhaps not surprisingly, despite the fact that an explanation of
$1/f$ noise was the main motivation for BTW’s introduction of sandpile models,
a good theoretical understanding of their time-dependent properties,
particularly the exact hydrodynamics and the related transport coefficients
governing the large-scale relaxations, remains lacking Dhar (2006); Dickman et
al. (2000); Redig (2006); Pradhan (2021); Ali (1995). Another fascinating
aspect of sandpiles has only recently been discovered. That is, the critical
state of sandpiles can be hyperuniform Hexner and Levine (2015, 2017), meaning
that the static sub-system mass fluctuations scale with sub-system size in an
anomalously slow manner Torquato and Stillinger (2003); Torquato (2016).
However, we currently have a limited knowledge of hyperuniform states of
matter, and one can only speculate as to how such a state emerges dynamically
in the first place Nizam et al. (2021). In fact, except for a few exact
results concerning static properties of the hyperuniform state Dandekar and
Dhar (2013); Dandekar (2020), there has been little theoretical progress in
this direction and a systematic approach for identifying the precise
microscopic dynamical origin of the anomalous fluctuations remains elusive. In
this scenario, a closer examination of the underlying dynamical mechanism that
results in such a state is desirable.

Here we address the above issues in the context of conserved stochastic
sandpiles and specifically focus on a continuous-time variant of the
celebrated Manna sandpile Manna (1991); Dickman et al. (2001), which has drawn
a lot of attention in the past Dickman et al. (2000). The conserved Manna
sandpile is a paradigm for systems, exhibiting a nonequilibrium absorbing
phase transition from a dynamically active state to an absorbing one having no
activities, upon tuning the global density. In fact, through several
simulation studies in the recent past, it is known that the critical state of
the conserved Manna sandpile is indeed hyperuniform Hexner and Levine (2015).
In another recent studies, it has been shown that the (near-)critical state is
characterized by the singular transport coefficients, leading to anomalous
relaxation and particle-transport in the system Chatterjee et al. (2018);
Tapader et al. (2021).

Historically, time-dependent properties of sandpiles have been studied in
terms of power spectra of dynamical activity such as instantaneous toppling
events in the systems Kertesz and Kiss (1990); Laurson et al. (2005). In the
original slowly-driven version of sandpiles, BTW had reported $1/f^{\psi}$
power-law behavior of the power spectrum for the activity, with the exponent
$\psi<2$ Bak et al. (1987), although their claim was refuted as several
simulation studies later found the exponent $\psi=2$ Jensen et al. (1989);
Kertesz and Kiss (1990). Subsequently, a more careful scaling analysis of
simulation data however revealed a nontrivial power spectra with the exponent
$\psi<2$ Laurson et al. (2005). On the theoretical front, a dynamic
renormalization-group analysis of phenomenological field-theoretical equations
describing a “running” sandpile (driven with nonzero rate of grain addition)
allowed an analytical calculation of the exponents $\psi$, involving activity
as well as output current, where $1/f$-type noise was observed, with
$\psi=-1$; quite interestingly, the temporal correlations in the long-time
regime were found to be anti-correlated, with $\psi=-1<0$. This could well be
the earliest signature of “dynamic hyperuniformity” in sandpiles. Recently,
this particular aspect of dynamic hyperuniformity, i.e., hyperuniformity in
the temporal domain, was also analyzed in a variant of the slowly-driven
sandpiles, called the Oslo ricepile Garcia-Millan et al. (2018). Similar low-
frequency behavior of the activity power spectrum with $\psi=-1$ (i.e., anti-
correlated) was observed in Ref. Yadav et al. (2012) for a directed
deterministic sandpile on a ladder with a finite driving rate; in the slow-
driving limit though, Maslov et al. Maslov et al. (1999) previously showed the
model to exhibit $1/f^{\psi}$ power spectrum, with $\psi=1$, for the total
mass fluctuation. In a slightly different study Baiesi and Maes (2006) of a
driven sandpile, albeit on a periodic domain, the power-spectrum of activity
had been found to be $1/f^{\psi}$ with $\psi=1$. In a conserved deterministic
lattice gas in two dimensions Giometto and Jensen (2012), subsystem mass
fluctuation was found to exhibit the power spectrum $S_{M}(f)\sim
f^{-\psi_{M}}$, where $\psi=1.5$ away from criticality and $\psi\approx 1.8$
near criticality.

In this paper, we theoretically investigate the time-dependent correlations
for current and mass in the (quasi-)steady state of the one-dimensional
conserved Manna sandpile. We establish a direct quantitative relationship
between various static and dynamic fluctuation properties in terms of the
density-dependent activity - the system’s “order parameter”, and its
derivative. We show that, in the thermodynamic limit, with system size
$L\rightarrow\infty$ and density $\bar{\rho}$ fixed, the variance of the local
(bond) current ${\mathcal{Q}(T)}$ up to time $T$ grows subdiffusively with
time. That is, we have $\left\langle\mathcal{Q}^{2}(T)\right\rangle\sim
T^{\alpha}$, where, away from criticality (in the time regime $T\ll L^{2}$),
the exponent $\alpha=1/2$ and, near criticality (in the regime $T\ll L^{z}$),
the current fluctuation is further suppressed with the exponent
$\alpha=1/2-\mu$, where $\mu=(\beta+1)/2\nu_{\perp}z>0$ and $\beta$,
$\nu_{\perp}$ and $z$ are the activity, correlation-length and dynamic
exponents, respectively; thus the anomalous suppression of the current
fluctuation near criticality serves as the dynamic precursor to the
hyperuniform state formed at the critical point. Indeed we find that the time-
dependent (two-point) correlation function for the current is long-ranged
(power-law) and negative, resulting in the low-frequency behavior of the
corresponding power spectrum $S_{\cal J}(f)\sim f^{\psi_{\mathcal{J}}}$, which
vanishes at low frequency, where $\psi_{\mathcal{J}}=1/2$ away from
criticality (strictly speaking, in the time regime $1/L^{2}\ll f\ll 1$ for
finite $L$) and $\psi_{\mathcal{J}}=1/2+\mu$ near criticality (in the time
regime $1/L^{z}\ll f\ll 1$). We show that the power spectrum $S_{M}(f)$ for
subsystem-mass fluctuation on the other hand diverges $S_{M}(f)\sim
f^{-\psi_{M}}$ at low frequency, where $\psi_{M}=3/2$ away from criticality
($1/L^{2}\ll f\ll 1$) and $\psi_{M}=3/2-\mu$ near criticality ($1/L^{z}\ll
f\ll 1$). These two exponents are not independent though and they are
connected by a scaling relation $\psi_{M}=2-\psi_{\mathcal{J}}$. In the
opposite limit of $T\rightarrow\infty$ (more specifically, $T\gg L^{2}$), with
system size $L$ finite (but, still large) and density $\bar{\rho}$ fixed, the
scaled current fluctuation in the steady state is shown to be proportional to
the activity $a(\bar{\rho})$ as
$\lim_{T\rightarrow\infty}\left\langle\mathcal{Q}^{2}(T)\right\rangle/T=2a(\bar{\rho})/L$.
On the other hand, the steady-state fluctuation of the subsystem current
$\bar{Q}(l,T)$, i.e., the cumulative (summed over bonds) current in a
subsystem of size $l$, in the thermodynamic limit ($L\rightarrow\infty$)
interestingly depends on the order of the limits taken. When the infinite-
subsystem-size limit is taken first and then the infinite-time limit, the
scaled subsystem current fluctuation
$\sigma_{Q}^{2}(\bar{\rho})\equiv\lim_{l\rightarrow\infty}\lim_{T\rightarrow\infty}\langle\bar{Q}^{2}(l,T)\rangle/lT=2a(\bar{\rho})$
converges to twice the activity. We derive a nonequilibrium fluctuation
relation, $\sigma^{2}(\bar{\rho})=\sigma_{Q}^{2}(\bar{\rho})/2D(\bar{\rho})$,
which connects the (scaled) subsystem mass fluctuation
$\sigma^{2}(\bar{\rho})=\langle M_{l}^{2}\rangle-\langle M_{l}\rangle^{2}$,
with $M_{l}$ being mass in a subsystem of size $l$ and $D(\bar{\rho})$ the
density-dependent bulk-diffusion coefficient, to the (scaled) subsystem
current fluctuation. The relation explains why the mass fluctuation should
vanish upon approaching criticality and helps characterize the dynamical
origin of hyperuniformity in the system.

We also study the steady-state mean-square displacements of tagged-particles,
which are characterized by the self-diffusion coefficient of the individual
particle. We theoretically show that the self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$ is identically equal to the ratio,
$a(\bar{\rho})/\bar{\rho}$, of the activity to the global number density of
the system, i.e., $\mathcal{D}_{s}(\bar{\rho})=a(\bar{\rho})/\bar{\rho}$, a
fluctuation relation, which connects the (scaled) tagged-particle displacement
fluctuation to the density-dependent activity. This relation immediately
explains a previous simulation observation of Ref. da Cunha et al. (2009)
that, upon approaching criticality, the self-diffusion coefficient in the
conserved Manna sandpile vanishes in the same fashion as the activity. Our
theoretical results are in a reasonably good agreement with simulations.

The plan of the paper is as follows. In sec. II.1, we define the conserved
Manna sandpile and the various quantities of interest. In sec. II.2, we
present our calculation method where we introduce an approximate truncation
scheme, helping us to calculate the dynamic correlations. Then in sec. II.3,
we study the dynamic properties of bond current fluctuations in the system and
provide a scaling argument to explain its behavior near criticality, followed
by the calculation of the corresponding power spectrum in sec. II.4.1. We then
proceed to calculate the variance of the cumulative subsystem (i.e., space-
time integrated) current and elucidate its relationship with the particle
mobility in sec. II.4.3. We study the self-diffusion coefficient of tagged
particles and the power spectrum for subsystem mass fluctuation in sec. II.5
and sec. II.6, respectively. Finally, in sec. III, we summarize with some
concluding remarks.

##  II Dynamic correlations in the steady state

###  II.1 Model and definitions

We consider the continuous-time variant Dickman et al. (2001) of the conserved
(“fixed energy”) Manna sandpile Manna (1991) on a ring of $L$ sites. Any site
$i$, with $i=0$, $1$, $\ldots$, $L-1$, can have $m_{i}\geq 0$ number of
particles, with $m_{i}=0$, $1$, $\ldots$, $N$; the total number of particles

$\displaystyle N=\sum_{i=0}^{L-1}m_{i},$ (1)

remains conserved; in this paper we throughout denote the global density as
$\bar{\rho}=N/L$. The dynamical rules are as following: An active site - a
site with $m_{i}>1$ \- topples with rate $1$ by randomly, and independently,
transferring each of the two particles to one of its nearest neighbours.

The sytem violates detailed balance and eventually reaches a nonequilibrium
(quasi-)steady state, which is not described by the familiar equilibrium
Boltzmann-Gibbs distribution and whose probability measure is apriori unknown.
The steady state of the system is usually charactarized through a global order
parameter, called the activity $a(\bar{\rho})$, defined as the density of
active sites,

$\displaystyle a(\bar{\rho})=\frac{\left\langle N_{a}\right\rangle}{L},$ (2)

where $N_{a}$ is the total number of active sites in the system and
$\langle.\rangle$ denotes the steady-state average. Interestingly the system
has a nontrivial spatio-temporal structure and, upon tuning the global density
$\bar{\rho}$, undergoes an absorbing phase transition. That is, above a
critical density $\rho_{c}$, there are dynamical activities in the system,
but, below the critical density, the dynamical activities in the steady state
cease and consequently there are no movements of particles in the system. The
absorbing phase transition in the conserved Manna sandpile has been intesively
studied in the past and can be charactarized by the following critical
exponents Dickman et al. (2001) \- the order-parameter exponent $\beta$, the
correlation-length exponent $\nu_{\perp}$ and the dynamic exponent $z$: Up on
approaching criticality from above, we have the following scaling behavior of
activity $a(\Delta)\sim\Delta^{\beta}$, correlation length
$\xi\sim\Delta^{-\nu_{\perp}}$, and the relaxation time $\tau_{r}\sim L^{z}$
where relative density $\Delta=\bar{\rho}-\rho_{c}>0$ and the critical density
$\rho_{c}\approx 0.94885$.

One can write update rules in an infinitesimal time-interval between time $t$
and $t+\differential{t}$ as given below,

$\displaystyle
m_{i}(t+\differential{t})=\begin{cases}m_{i}(t)+1&\frac{1}{2}\hat{a}_{i+1}\differential{t}\\\
m_{i}(t)+1&\frac{1}{2}\hat{a}_{i-1}\differential{t}\\\
m_{i}(t)+2&\frac{1}{4}\hat{a}_{i+1}\differential{t}\\\
m_{i}(t)+2&\frac{1}{4}\hat{a}_{i-1}\differential{t}\\\
m_{i}(t)-2&\hat{a}_{i}\differential{t}\\\
m_{i}(t)&[1-\Sigma\differential{t}],\end{cases}$ (3)

where $\Sigma=(3/4)(\hat{a}_{i+1}+\hat{a}_{i-1})+\hat{a}_{i}$; here
$\hat{a}_{i}$ is an indicator function with $\hat{a}_{i}=1$ if the site is
active and $\hat{a}_{i}=0$ otherwise. Using the above update rules eq.(3), we
can write the time evolution equation of the first moment of local mass as

$\displaystyle\derivative{t}\left\langle
m_{i}(t)\right\rangle=\left[\left\langle\hat{a}_{i-1}(t)-2\hat{a}_{i}(t)+\hat{a}_{i+1}(t)\right\rangle\right].$
(4)

Denoting the local density $\rho_{i}(t)=\left\langle m_{i}(t)\right\rangle$,
we can alternatively write the above equation as

$\displaystyle\derivative{t}\rho_{i}(t)=\sum_{k}\Delta_{i,k}a_{k}(t)$ (5)

where $\Delta_{i,k}$ is the discrete Laplacian and
$a_{k}(t)=\left\langle\hat{a}_{k}\right\rangle(t)$ is the average
instantaneous activity. On the large spatio-temporal scales and by taking the
diffusive scaling limit $i\rightarrow x=i/L$ and $t\rightarrow\tau=t/L^{2}$,
we can write the hydrodynamic time-evolution equation for the local density
field $\rho(x,\tau)$ as in eq.(5) Chatterjee et al. (2018); Tapader et al.
(2021),

$\displaystyle\frac{\partial\rho(x,\tau)}{\partial\tau}=\frac{\partial^{2}a(\rho)}{\partial
x^{2}}\equiv\partialderivative{x}\quantity(D(\rho)\partialderivative{\rho}{x}),$
(6)

where $D(\rho)$ is the density-dependent bulk-diffusion coefficient. It has
been previously demonstrated in Ref. Chatterjee et al. (2018); Tapader et al.
(2021), the bulk-diffusion coefficient can be written in terms of the
derivative of the activity $a(\rho)$ wrt density $\rho$,

$\displaystyle D(\rho)=\frac{da(\rho)}{d\rho}\equiv a^{\prime}(\rho).$ (7)

Indeed the relaxation processes occuring on a large (coarse-grained) scale are
primarily governed by the bulk-diffusion coefficient - the fact that we later
use to introduce a truncation scheme [see eq.(23)] for analytically
calculating various time-dependendent correlation functions, which would not
have been possible otherwise in a system with nontrivial correlations as in
sandpiles. Notably, one can recast the density evolution equation (5) in a
microscopic form of the continuity equation,

$\displaystyle\derivative{t}\rho_{i}(t)=\left\langle\mathcal{J}_{i}(t)-\mathcal{J}_{i+1}(t)\right\rangle,$
(8)

where the microscopic instantaneous current $\mathcal{J}_{i}(t)$ is defined as
the number of particles crossing a bond $(i,i+1)$ in an infinitesimal time
interval $\quantity(t,t+\differential{t})$. It is useful to define a related
observable - the cumulative, or time-integrated, bond current
$\mathcal{Q}_{i}(t)$ upto time $t$, which is used to calculate various other
correlation functions, such as that involving mass and activity, and is easily
measured in simulations. At the microscopic level, the time-integrated current
$\mathcal{Q}_{i}(t)$ is defined as the total number of particles transferred
accross $i^{th}$ bond, connecting the nearest-neighbour pair of sites
$\quantity(i,i+1)$, during a time interval $\quantity[0,t]$. That is, the
time-integrated current accross the $i^{th}$ bond during an infinitesimal time
interval $\quantity[t,t+\differential{t}]$ is simply
$\mathcal{J}_{i}(t)\differential{t}$ with

$\displaystyle\mathcal{J}_{i}(t)=\lim_{\Delta t\rightarrow
0}\frac{\quantity(\mathcal{Q}_{i}(t+\Delta t)-\mathcal{Q}_{i}(t))}{\Delta
t}\equiv\frac{d\mathcal{Q}_{i}(t)}{\differential{t}},$ (9)

and

$\displaystyle\mathcal{\mathcal{Q}}_{i}(T)=\int^{T}_{0}\differential{t}\mathcal{J}_{i}(t).$
(10)

On the average level, we therefore have

$\displaystyle\left\langle\mathcal{J}_{i}(t)\right\rangle=\left\langle\frac{d\mathcal{Q}_{i}(t)}{\differential{t}}\right\rangle=\frac{d\left\langle\mathcal{Q}_{i}(t)\right\rangle}{\differential{t}}.$
(11)

We now decompose the instantaneous current into two parts \- a diffusive
component $\mathcal{J}_{i}^{(d)}(t)$ and a fluctuating component
$\mathcal{J}_{i}^{(fl)}(t)$,

$\displaystyle\mathcal{J}_{i}(t)=\mathcal{J}_{i}^{(d)}(t)+\mathcal{J}_{i}^{(fl)}(t),$
(12)

where, motivated by eq.(4), we identify the diffusive current as

$\mathcal{J}_{i}^{(d)}(t)=\hat{a}_{i}(t)-\hat{a}_{i+1}(t).$ (13)

Indeed, as we see later, in that case only the diffusive-current component
possesses long-ranged temporal correlations, varying slowly in time as a power
law, and the fluctuating component on the other hand is simply a delta-
correlated one, explaining the motivation behind the above decomposition. Note
that, due to the fact
$\left\langle\mathcal{J}_{i}\right\rangle=\left\langle\mathcal{J}_{i}^{(d)}\right\rangle=\left\langle\hat{a}_{i}\right\rangle-\left\langle\hat{a}_{i+1}\right\rangle$,
we must have $\left\langle\mathcal{J}_{i}^{(fl)}\right\rangle=0$. Indeed the
fluctuating current component can be related to the (conserve) noise term in
an appropriately coarse-grained fluctuating hydrodynamic theory, which can be
then used to study the large-scale fluctuation properties of the system
Bertini et al. (2015); Sadhu and Derrida (2016).

In the following sections, we study the fluctuation properties of various
compenents of currents, instantaneous and the fluctuating one as decomposed in
eq.13. Throughout the paper, we use the following notation for correlation
function $C_{r}^{AB}(t,t^{\prime})$ involving any two local observables
$A_{i}(t)$ and $B_{j}(t^{\prime})$ with $t\geq t^{\prime}$,

$\displaystyle C_{r=|i-j|}^{AB}(t,t^{\prime})=\left\langle
A_{i}(t)B_{j}(t^{\prime})\right\rangle-\left\langle
A_{i}(t)\right\rangle\left\langle B_{j}(t^{\prime})\right\rangle,$ (14)

where $r=|j-i|$ is the relative distance. Further, we denote the spatial
Fourier transform of the correlation function $C_{r}^{AB}(t,t^{\prime})$ as

$\displaystyle\tilde{C}_{q}^{AB}(t,t^{\prime})=\sum\limits_{r=0}^{L-1}C_{r}^{AB}(t,t^{\prime})e^{{\bf\it
i}qr},$ (15)

where $q=2\pi k/L$ and $k=0$, $1$, $\ldots$, $L-1$ and the inverse Fourier
transform as

$\displaystyle
C_{r}^{AB}(t,t^{\prime})=\frac{1}{L}\sum\limits_{q}\tilde{C}_{q}^{AB}(t,t^{\prime})e^{-{\bf\it
i}qr}.$ (16)

By introducing a truncation scheme as discussed below, we can theoretically
compute the statistics of different combinations of various local currents
$\mathcal{J}_{i}$, $\mathcal{J}_{i}^{(d)}$, $\mathcal{J}_{i}^{(fl)}$ as well
as mass $m_{i}$, essentially in terms of the following two correlation
functions -
$\left\langle\mathcal{Q}_{i}(t)\mathcal{Q}_{j}(t^{\prime})\right\rangle$ and
$\left\langle m_{i}(t)\mathcal{Q}_{j}(t^{\prime})\right\rangle$.

###  II.2 Theory

For the conserved Manna sandpile, we write the stochastic update equation of
the integrated current $\mathcal{Q}_{i}(t)$ in an infinitesimal time
(continuous) interval $\quantity[t,t+\differential{t}]$,

$\displaystyle\mathcal{Q}_{i}(t+\differential{t})=\begin{cases}\textbf{{events}}&\textbf{{probabilities}}\\\
\mathcal{Q}_{i}(t)+1&\frac{1}{2}\hat{a}_{i}(t)\differential{t}\\\
\mathcal{Q}_{i}(t)+2&\frac{1}{4}\hat{a}_{i}(t)\differential{t}\\\
\mathcal{Q}_{i}(t)-1&\frac{1}{2}\hat{a}_{i+1}(t)\differential{t}\\\
\mathcal{Q}_{i}(t)-2&\frac{1}{4}\hat{a}_{i+1}(t)\differential{t}\\\
\mathcal{Q}_{i}(t)&1-\Sigma\differential{t}\end{cases}$ (17)

where $\Sigma=\quantity(\hat{a}_{i}(t)-\hat{a}_{i+1}(t))$. Using the above
update rules, the time-evolution equation for the first moment of the time-
integrated current $\mathcal{Q}_{i}(t)$ can be written as

$\displaystyle\derivative{t}\left\langle\mathcal{Q}_{i}(t)\right\rangle=\left\langle\hat{a}_{i}(t)\right\rangle-\left\langle\hat{a}_{i+1}(t)\right\rangle.$
(18)

Similarly, using eq.(17), we find the time-evolution equation for the second
moment
$\left\langle\mathcal{Q}_{i}(t)\mathcal{Q}_{i+r}(t^{\prime})\right\rangle=C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})$
of the integrated current at two different times $t$ and $t^{\prime}$, for
$t>t^{\prime}$, as given below (see appendix C.1 for details),

$\displaystyle\derivative{t}C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})=\quantity(C^{\hat{a}\mathcal{Q}}_{r}(t,t^{\prime})-C^{\hat{a}\mathcal{Q}}_{r-1}(t,t^{\prime})).$
(19)

The above equation, which is central to our study, is however difficult to
solve exactly due to an infinite hierarchy of correlation functions involved
and so we employ below an approximation scheme.

We note that the evolution of space and time dependendent activity-current
correlation function $C^{\hat{a}\mathcal{Q}}_{r}(t,t^{\prime})$ appearing in
eq.(19) contains higher-order correlation functions, involving activity,
current and some other observables. The calculations of the higher-order
correlations would eventually result in a rapidly growing complexity in the
hierarchy of correlation functions, which do not constitute a closed set of
equations. More specifically, one can start with the infinitesimal time update
equation for $\hat{a}_{i}$ itself,

$\displaystyle\hat{a}_{i}(t+\differential{t})=\begin{cases}\textbf{{events}}&\textbf{{probabilities}}\\\
\hat{a}_{i}(t)+1&\frac{1}{2}\hat{a}_{i+1}(t)\hat{p}_{i}(t)\delta_{m_{i},1}\differential{t}\\\
\hat{a}_{i}(t)+1&\frac{1}{2}\hat{a}_{i-1}(t)\hat{p}_{i}(t)\delta_{m_{i},1}\differential{t}\\\
\hat{a}_{i}(t)+1&\frac{1}{4}\hat{a}_{i+1}(t)\hat{p}_{i}(t)\differential{t}\\\
\hat{a}_{i}(t)+1&\frac{1}{4}\hat{a}_{i-1}(t)\hat{p}_{i}(t)\differential{t}\\\
\hat{a}_{i}(t)-1&\hat{a}_{i}(t)\quantity(\delta_{m_{i},2}+\delta_{m_{i},3})\differential{t}\\\
\hat{a}_{i}(t)&1-\Sigma\differential{t},\end{cases}$ (20)

where $\hat{p}_{i}(t)=\quantity(1-\hat{a}_{i}(t))$ and
$\Sigma=\quantity(\hat{p}_{i}(t)/2)(\delta_{m_{i},1}+1/2)\quantity[\hat{a}_{i+1}+\hat{a}_{i-1}]+\hat{a}_{i}\quantity(\delta_{m_{i},2}+\delta_{m_{i},3})$.
From this update equation, for $t>t^{\prime}$, we can write the evolution
equation for $C^{a\mathcal{Q}}_{r}(t,t^{\prime})$ as

$\displaystyle\derivative{t}C^{a\mathcal{Q}}_{r}(t,t^{\prime})\equiv\derivative{t}\left\langle\hat{a}_{0}(t)\mathcal{Q}_{r}(t^{\prime})\right\rangle_{c}=\left\langle\quantity{\derivative{t}\hat{a}_{0}(t)}\mathcal{Q}_{r}(t^{\prime})\right\rangle_{c},$
(21)

where the observable inside the curly bracket evolves according to the
equation,

$\displaystyle\derivative{t}\hat{a}_{0}(t)=-\hat{a}_{0}(t)\quantity(\delta_{m_{0},2}+\delta_{m_{0},3})+$
$\displaystyle\quantity[\hat{p}_{0}(t)\quantity(\hat{a}_{1}(t)+\hat{a}_{L-1}(t))]\quantity(\frac{1}{2}\delta_{m_{0},1}+\frac{1}{4}).$
(22)

Clearly, to solve for $C_{r}^{a\mathcal{Q}}(t,t^{\prime})$, one needs to
calculate the correlation functions
$\left\langle\delta_{m_{0},2}\hat{a}_{0}(t)\mathcal{Q}_{r}(t^{\prime})\right\rangle$
and
$\left\langle\delta_{m_{0},1}\hat{a}_{0}(t)\hat{a}_{1}(t)\mathcal{Q}_{r}(t^{\prime})\right\rangle$,
which would in turn involve another set of even higher-order correlation
functions; of course, in this way, one generates an infinite hierarchy of
equations, which is difficult to handle.

We bypass the difficulty by employing the following approximation scheme,
which in fact truncates the otherwise infinite hierarchy. To this end, we
approximate the local diffusive current, which is the gradient of the
instantaneous local activity and written as

$\displaystyle\mathcal{J}_{i}^{(d)}(t)=\quantity[\hat{a}_{i}(t)-\hat{a}_{i+1}(t)]\equiv
D(\bar{\rho})\quantity[m_{i}(t)-m_{i+1}(t)],~{}~{}~{}~{}$ (23)

where the bulk-diffusion coefficent is given by
$D(\bar{\rho})=a^{\prime}(\bar{\rho})$ as in eq.(7) and can be treated as a
constant. Essentially, we have assumed in the above truncation scheme that
fluctuations around the global density $\bar{\rho}$ are small and the local
current (the gradient of activity) relaxes diffusively via the gradient of
local mass. As demonstrated later, this approximation goes beyond a simple
mean-field theory and captures relevant correlations quite well on the large
(hydrodynamic) time scales. More specifically, to compute correlations between
any observable $A(t)$ (e.g.,, current or mass) and the current observable
$\mathcal{J}_{i}^{(d)}(t)$, we replace the diffusive-current observable
$\mathcal{J}_{i}^{(d)}(t)$ by its truncated form as given in the rhs of
eq(23), to obtain the following equality (approximate),

$\langle A(t)\mathcal{J}_{i}^{(d)}(t^{\prime})\rangle\simeq
a^{\prime}(\bar{\rho})\quantity[A(t)\\{m_{i}(t^{\prime})-m_{i+1}(t^{\prime})\\}].$
(24)

Now we can proceed further by first rewriting eq.(19) as

$\displaystyle\derivative{t}C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})=a^{\prime}(\bar{\rho})\quantity(C^{m\mathcal{Q}}_{r}(t,t^{\prime})-C^{m\mathcal{Q}}_{r-1}(t,t^{\prime})),$
(25)

and then expressing the time-evolution equation for the mass-current
correlations $C^{m\mathcal{Q}}_{r}(t,t^{\prime})=\left\langle
m_{i}(t)\mathcal{Q}_{i+r}(t^{\prime})\right\rangle-\langle
m_{i}(t)\rangle\langle\mathcal{Q}_{i+r}(t^{\prime})\rangle$ as

$\displaystyle\derivative{t}C^{m\mathcal{Q}}_{r}(t,t^{\prime})\simeq
a^{\prime}(\bar{\rho})\sum\limits_{k}\Delta_{r,k}C^{m\mathcal{Q}}_{k}(t,t^{\prime}),$
(26)

where have used eq.(23) in the intermediate steps (see appendix sec. C.2). It
is worth noting that, in eq.(23) or in eqs.(25) and (26), the activity appears
simply as a global density-dependent constant prefactor $a(\bar{\rho})$ and
thus we obtain a closed set of equations, involving only mass and current
correlations; equations (25) and (26) can be solved, albeit in terms of the
activity $a(\bar{\rho})$, which however remains undetermined in our theory.
Interestingly, it was previously possible to exactly calculate the dynamic
correlation functions in simple exclusion processes Sadhu and Derrida (2016)
because, in that case, one already gets a closed set of equations for mass and
current correlations and, furthermore, because the steady-state measure is a
product one, allowing one to calculate various static (time-independent)
quantities, which enter into the expression of the current and mass
fluctuations. Without an explicit knowledge of the steady state measure Sadhu
and Dhar (2009), the Manna sandpile, on the other hand, is nontrivial due to
the nonzero spatial correlations present in the system, making explicit
calculations of the static quantities, such as the density-dependent activity,
quite difficult. Nevertheless, as shown below, by extending the formalism
developed in the context of simple exclusion processes Sadhu and Derrida
(2016), one can calculate various dynamic correlations in terms of activity
and obtain fluctuation relations, which precisely quantify the underlying
relationship between dynamic and static fluctuations in the system.

At this stage, it is useful to introduce the Fourier representation of the
respective correlation functions as given in eq.(15), and we can then write
eqs.(25) and (26) in the respective Fourier modes,

$\displaystyle\derivative{t}\tilde{C}_{q}^{\mathcal{Q}\mathcal{Q}}(t,t^{\prime})=a^{\prime}(\bar{\rho})\tilde{C}_{q}^{m\mathcal{Q}}(t,t^{\prime})\quantity[1-e^{{\bf\it
i}q}],$ (27)

and

$\displaystyle\derivative{t}\tilde{C}_{q}^{m\mathcal{Q}}(t,t^{\prime})=-a^{\prime}(\bar{\rho})\lambda_{q}\tilde{C}_{q}^{m\mathcal{Q}}(t,t^{\prime}),$
(28)

where

$\lambda_{q}=2\quantity[1-\cos q].$ (29)

Now, eqs.(27) and (28), can be integrated to have

$\displaystyle\tilde{C}_{q}^{\mathcal{Q}\mathcal{Q}}(t,t^{\prime})=$
$\displaystyle\int\limits_{t^{\prime}}^{t}\differential{t^{\prime\prime}}a^{\prime}(\bar{\rho})\tilde{C}^{m\mathcal{Q}}_{q}(t^{\prime\prime},t^{\prime})\quantity[1-e^{{\bf\it
i}q}]$ $\displaystyle+$
$\displaystyle\tilde{C}_{q}^{\mathcal{Q}\mathcal{Q}}(t^{\prime},t^{\prime}),$
(30)

and

$\displaystyle\tilde{C}^{m\mathcal{Q}}_{q}(t,t^{\prime})=e^{-a^{\prime}(\bar{\rho})\lambda_{q}(t-t^{\prime})}\tilde{C}^{m\mathcal{Q}}_{q}(t^{\prime},t^{\prime}).$
(31)

respectively. However, to fully solve for the unequal-time correlation
functions $\tilde{C}_{q}^{\mathcal{Q}\mathcal{Q}}(t,t^{\prime})$ and
$\tilde{C}^{m\mathcal{Q}}_{q}(t,t^{\prime})$, we need to calculate their
respective equal-time counterparts as well. First we obtain the evolution
equation for the equal-time mass-current correlation function
$C^{m\mathcal{Q}}_{r}(t^{\prime},t^{\prime})$ and, then writing
$C^{m\mathcal{Q}}_{r}(t^{\prime},t^{\prime})$ in the Fourier space (using
eq.(148); see appendix C.4 for details), we get

$\displaystyle\derivative{t^{\prime}}\tilde{C}^{m\mathcal{Q}}_{q}(t^{\prime},t^{\prime})=-a^{\prime}(\bar{\rho})\lambda_{q}\tilde{C}_{q}^{m\mathcal{Q}}(t^{\prime},t^{\prime})+\tilde{f}_{q}(t^{\prime}),$
(32)

where the Fourier transform of the source term $\tilde{f}_{q}(t^{\prime})$ in
the steady state is given by

$\displaystyle\tilde{f}_{q}=\tilde{C}^{m\hat{a}}_{q}\quantity(1-e^{-{\bf\it
i}q})-2a(\bar{\rho})\quantity(1-e^{-{\bf\it
i}q})\quantity[1+\frac{\lambda_{q}}{4}].$ (33)

Now eq.(32) can be directly integrated to obtain

$\displaystyle\tilde{C}^{m\mathcal{Q}}_{q}(t^{\prime},t^{\prime})=\int\limits_{0}^{t^{\prime}}\differential{t^{\prime\prime}}e^{-a^{\prime}(\bar{\rho})\lambda_{q}(t^{\prime}-t^{\prime\prime})}\tilde{f}_{q}(t^{\prime\prime}),$
(34)

substituting which into eq.(31), we get the unequal-time mass-current
correlation function in terms of $\tilde{f}_{q}$,

$\displaystyle\tilde{C}^{m\mathcal{Q}}_{q}(t,t^{\prime})=\int\limits_{0}^{t^{\prime}}\differential{t^{\prime\prime}}e^{-a^{\prime}(\bar{\rho})\lambda_{q}(t-t^{\prime\prime})}\tilde{f}_{q}(t^{\prime\prime}).$
(35)

To calculate the above correlation, we need to calculate the activity-mass
correlation as in eq.(33). Importantly, as shown below, we can calculate the
static (time-independent) activity-mass correlation function exactly in the
steady state (see appendix C.5). Using the steady-state condition
$dC^{mm}_{r}(t,t)/dt=0$, we obtain

$\displaystyle\derivative{t}C_{r}^{mm}(t,t)=$
$\displaystyle\sum\limits_{k}2\left\langle
m_{0}\Delta_{rk}\hat{a}_{k}\right\rangle_{c}+B_{r}=0,$ (36)

where $B_{r}$ is the source term having the form,

$\displaystyle B_{r}(\bar{\rho})=$ $\displaystyle
7a(\bar{\rho})\delta_{0,r}-4a(\bar{\rho})\quantity(\delta_{0,r+1}+\delta_{0,r-1})+$
$\displaystyle\frac{a(\bar{\rho})}{2}\quantity(\delta_{0,r+2}+\delta_{0,r-2}).$
(37)

Equation (36) can be solved by employing a generating function,

$\displaystyle G\quantity(z)=\sum\limits_{r=0}^{\infty}C^{m\hat{a}}_{r}z^{r},$
(38)

for the equal-time mass-activity correlation; see appendix.(C.5). Here we
directly provide the solution of the generating function in terms of the
static density-dependent activity,

$\displaystyle
G\quantity(z)=\frac{3a(\bar{\rho})}{2}-\frac{a(\bar{\rho})}{4}z,$ (39)

implying the mass-activity static correlation to be

$\displaystyle
C^{m\hat{a}}_{r}=\begin{cases}\frac{3a(\bar{\rho})}{2}&\text{for}\hskip
4.0ptr=0,\\\ -\frac{a(\bar{\rho})}{4}&\text{for}\hskip
4.0pt\absolutevalue{r}=1,\\\ 0&\text{otherwise},\end{cases}$ (40)

which is in fact exact. Then, by writing Fourier transform of the above
equation, we have

$\displaystyle\tilde{C}_{q}^{m\hat{a}}=a(\bar{\rho})+\frac{a(\bar{\rho})}{4}\lambda_{q},$
(41)

and substituting the above in eq.(33), we straightforwardly obtain

$\displaystyle\tilde{f}_{q}=-a(\bar{\rho})\quantity(1-e^{-{\bf\it
i}q})\quantity(1+\frac{\lambda_{q}}{4}).$ (42)

###  II.3 Time-integrated current fluctuation

In this section, we calculate the time-integrated bond-current fluctuation, by
using the theory developed in the previous section. To this end, we substitute
eq.(42) into eq.(35) and get an explicit solution the first term on the rhs of
eq.(II.2). Similarly, we can calculate the second term on the rhs of eq.(II.2)
as given below,

$\displaystyle
C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t)\simeq\int\limits_{0}^{t}\differential{t^{\prime}}\Gamma_{r}(t^{\prime})+$
$\displaystyle\frac{1}{L}\int\limits_{0}^{t}\differential{t^{\prime}}a^{\prime}(\bar{\rho})\sum\limits_{q=0}^{L-1}\tilde{C}^{m\mathcal{Q}}_{q}(t^{\prime},t^{\prime})\quantity[1-e^{{\bf\it
i}q}]\quantity(2-\lambda_{qr}),$ (43)

where $\lambda_{qr}=2[1-\cos(qr)]$; see appendix sec C.3 for details. Here the
quantity $\Gamma_{r}(t)$ is the strength of the steady-state correlation
function for the fluctuating current
$\left\langle\mathcal{J}_{0}^{(fl)}(t)\mathcal{J}_{r}^{(fl)}(t^{\prime})\right\rangle=\Gamma_{r}(t)\delta(t-t^{\prime})$
as derived later in eq.(59). Since we are interested only in the steady-state
properties, the strength $\Gamma_{r}$ is replaced by its steady-state value,

$\displaystyle\Gamma_{r}(\bar{\rho})=3a(\bar{\rho})\delta_{0,r}-\frac{a(\bar{\rho})}{2}\quantity(\delta_{0,r+1}+\delta_{0,r-1}),$
(44)

which depends on global density $\bar{\rho}$, through the density-dependent
activity $a(\bar{\rho})$; for the detail calculation of the strength
$\Gamma_{r}$, see section II.4.1. Now the Fourier transform of eq.(II.3) leads
to the second term of rhs of eq.(II.2). By using the inverse Fourier transform
of eq.(II.2), we finally obtain the desired space- and time-dependent current
correlation function in the steady-state,

$\displaystyle C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})=$ $\displaystyle
t^{\prime}\Gamma_{r}-a^{\prime}(\bar{\rho})a(\bar{\rho})\frac{1}{L}\sum\limits_{q}\int\limits_{0}^{t^{\prime}}\differential{t^{\prime\prime}}\int\limits_{0}^{t^{\prime\prime}}\differential{t^{\prime\prime\prime}}e^{-a^{\prime}(\bar{\rho})\lambda_{q}(t^{\prime\prime}-t^{\prime\prime\prime})}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})(2-\lambda_{qr})$
$\displaystyle-a^{\prime}(\bar{\rho})a(\bar{\rho})\frac{1}{L}\sum\limits_{q}\int\limits_{t^{\prime}}^{t}\differential{t^{\prime\prime}}\int\limits_{0}^{t^{\prime}}\differential{t^{\prime\prime\prime}}e^{-a^{\prime}(\bar{\rho})\lambda_{q}(t^{\prime\prime}-t^{\prime\prime\prime})}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})e^{-{\bf\it
i}qr}.$ (45)

The asymptotic behavior of the above equation can be straightforwardly
obtained as given below (see appendix B.1 for details),

$\displaystyle\left\langle\mathcal{Q}^{2}(T)\right\rangle\simeq\begin{cases}\frac{2a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}T^{\frac{1}{2}}&\text{for}\hskip 4.0pt1\ll T\ll
L^{2},\\\ \frac{2a(\bar{\rho})}{L}T&\text{for}\hskip 4.0ptT\gg
L^{2}.\end{cases}$ (46)

In simulation, we verify a special case of eq.(II.3) by putting $r=0$ and
$t=t^{\prime}\equiv T$, i.e., the time-integrated bond current fluctuation
$\left\langle\mathcal{Q}^{2}(T)\right\rangle\equiv
C_{0}^{\mathcal{Q}\mathcal{Q}}(T,T)$ (here average current
$\left\langle\mathcal{Q}(T)\right\rangle=0$ for the steady-state measurement).
In fig.1, we plot $\left\langle Q^{2}(T)\right\rangle$, obtained from
simulation (plotted in solid lines), for various densities $\bar{\rho}=2.0$
(red line), $\bar{\rho}=1.5$ (yellow line), $\bar{\rho}=1.2$ (blue line),
$\bar{\rho}=1.0$ (green line), $\bar{\rho}=0.97$ (violet line) as a function
of $T$.

Figure 1: Scaled fluctuations of cumulative (time-integrated) bond-current up
to time $T$, obtained from simulations (solid lines), is plotted as a function
of scaled time $T/L^{2}$ for different densities $\bar{\rho}=0.97$ (violet),
$\bar{\rho}=1.0$ (green), $\bar{\rho}=1.2$ (blue), $\bar{\rho}=1.5$ (yellow),
$\bar{\rho}=2.0$ (red), and for system size $L=1000$; theory as in eq.(II.3)
with $r=0$, $t=t^{\prime}=T$ (black dotted line) is in an excellent agreement
with simulation for $\bar{\rho}=2.0$. Three guiding lines signify the initial-
time subdiffusive growth $\left\langle Q^{2}(T)\right\rangle\sim T^{1/2}$ [as
in the first part of eq.(46)] and the late-time diffusive growth $\left\langle
Q^{2}(T)\right\rangle\sim T$ [as in the second part of eq.(46)] away from
criticality, and the initial-time anomalously suppressed subdiffusive growth
$\left\langle Q^{2}(T)\right\rangle\sim T^{1/2-\mu}$ [as in eq.(48)] near
criticality.

Indeed the dynamical behaviors as predicted by the asymptotics in eq.(46) are
different in two different time regimes: On smaller initial time scales $1\ll
T\ll L^{2}$, the time-integrated current grows sub-diffusively as $T^{1/2}$
and, on larger (hydrodynamic) time scales $T\gg L^{2}$, grows linearly as $T$.
We compare the simulation result with that obtained from our theory eq.(II.3)
with $r=0$ and $t=t^{\prime}=T$, for $\bar{\rho}=2.0$ (black dashed line) and
for system size $L=1000$; one can see an excellent agreement between
simulation and theory.

As mentioned previously, our theory is expected to be valid at hydrodynamic
times ($T\gg L^{z}$) and the small-time ($T\ll L^{z}$) behavior of current
fluctuation, especially near criticality, is not quite well captured by
eq.II.3). However the small-time behavior near criticality can still be
obtained qualitatively by using the following standard scaling analysis.
Indeed, first resorting to a simple dimensional analysis, we can see that the
activity scales as $a(\Delta)\sim\Delta^{\beta}\sim T^{-\beta/\nu_{\perp}z}$,
where the relative density $\Delta=\bar{\rho}-\rho_{c}\ll 1$ and we use the
following scaling relations: correlation length $\xi\sim T^{1/z}$ the relative
density $\Delta\sim\xi^{-1/\nu_{\perp}}\sim T^{-1/\nu_{\perp}z}$, with $z$
being the dynamic exponent. Thus, by writing
$a(\Delta)/\sqrt{a^{\prime}(\Delta)}\sim T^{-\mu}$, we straightforwardly have,
in the initial-time regime $1\ll T\ll L^{z}$, the scaling behavior of the
current fluctuation
$\left\langle\mathcal{Q}^{2}(T)\right\rangle\sim(a/\sqrt{a^{\prime}})T^{1/2}\sim
T^{\alpha}$, where the exponents

$\displaystyle\alpha=\frac{1}{2}-\mu$ (47)

and

$\displaystyle\mu=\frac{\beta+1}{2\nu_{\perp}z}.$ (48)

More precisely, near criticality we expect the following scaling form for the
time-integrated bond-current fluctuation to hold,

$\displaystyle\left\langle\mathcal{Q}^{2}(T)\right\rangle\simeq L^{\alpha
z}\mathcal{G}\quantity(\Delta
L^{1/\nu_{\perp}},\frac{T}{L^{z}})=T^{\alpha}\mathcal{F}\quantity(\Delta
L^{1/\nu_{\perp}},\frac{T}{L^{z}}),$ (49)

where $\mathcal{G}$ and $\mathcal{F}$ are two scaling functions. To determine
the exponent $\alpha$ in eq.(49) from simulations, we take density value such
that $\Delta L^{1/\nu_{\perp}}\rightarrow 0$ and then we plot in Fig. 2 the
scaled variance of time-integrated bond current
$T^{-\alpha}\left\langle\mathcal{Q}^{2}(T)\right\rangle$ as a function of the
scaled time $T/L^{z}$ for rather quite smaller system sizes $L=100$ (violet),
$L=200$ (green) and $L=500$ (blue), respectively and for $\bar{\rho}=0.95$. We
get a reasonably good scaling collapse of simulation data, with the exponent
estimated to be $\alpha\approx 0.297$; our theoretical prediction of the
exponent $\alpha\approx 0.26$, computed from scaling relations eqs.(47) and
(48) by using $\beta\approx 0.42$, $z\approx 1.66$ and $\nu_{\perp}\approx
1.81$ Dickman et al. (2001), slightly underestimates that obtained from
simulation though.

Figure 2: The scaled variance
$T^{-\alpha}\left\langle\mathcal{Q}^{2}(T)\right\rangle$ of time-integrated
bond current $\mathcal{Q}(T)$ up to time $T$, obtained from simulations, is
plotted as a function of scaled time $T/L^{z}$ for system size $L=100$
(magenta), $L=200$ (green) and $L=500$ (sky-blue) and for (near-critical)
density $\bar{\rho}=0.95$, where, to achieve the scaling collapse, we use
$\alpha\approx 0.297$ and $z\approx 1.66$. The value of $\alpha\approx 0.297$
obtained from simulations is not far from $\alpha\approx 0.26$ obtained from
theory as in eqs.(47) and (48) with $\beta\approx 0.42$, $\nu_{\perp}\approx
1.81$ Dickman et al. (2001).

The variance of time-integrated bond current is known to grow subdiffusively
in the initial-time regime in diffusive systems, such as symmetric simple
exclusion processes De Masi and Ferrari (2002); Sadhu and Derrida (2016). Now,
away from criticality, the relaxation processes in sandpiles are diffusive
Chatterjee et al. (2018); Tapader et al. (2021) and, therefore not
surprisingly, the current fluctuation exhibits a subdiffusive growth, as
derived in the first part of eq.(46). Near-critical relaxation processes in
the Manna sandpile, on the other hand, are anomalous and we observe nontrivial
scaling behavior. When compared to a normal diffusive system, the Manna
sandpile near criticality exhibits strong suppression of current fluctuations
and much slower subdiffusive growth of temporal fluctuations due to the lack
of local activity, as described in eq.(49); we call it a “dynamic
hyperuniformity”, which is quite analogous to hyperuniformity studied in the
spatial domain Torquato and Stillinger (2003); Hexner and Levine (2015) and is
reminiscent of that identified in the context of temporal statistics of
avalanches in a particular version sandpile, called the Oslo ricepile Garcia-
Millan et al. (2018). Of course, the dynamic hyperuniformity, or the
anomalously subdiffusive growth of temporal fluctuations, can be equivalently
characterized in terms of the current and mass power spectra, or the
respective dynamic correlation functions, as discussed in the following
sections.

###  II.4 Current fluctuation and its power spectrum

####  II.4.1 Instantaneous current

In this section, we calculate in the steady-state the time-dependent (unequal-
time) current-current correlation function
$C_{r}^{\mathcal{J}\mathcal{J}}(t)\equiv
C_{r}^{\mathcal{J}\mathcal{J}}(t,t^{\prime}=0)$ of the instantaneous bond
current by taking time derivative of time-integrated bond current correlation
as given below,

$\displaystyle
C^{\mathcal{J}\mathcal{J}}_{r}(t)=\quantity[\derivative{t}\derivative{t^{\prime}}C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})]_{t^{\prime}=0},$
(50)

where $t\geq t^{\prime}$. Now, after differentiating eq.(II.3), we can write
the time-dependent current correlation as

$\displaystyle C^{\mathcal{J}\mathcal{J}}_{r}(t)=\Gamma_{r}\delta(t)-$
$\displaystyle
a^{\prime}(\bar{\rho})a(\bar{\rho})\quantity[\frac{1}{L}\sum\limits_{q}e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})e^{-{\bf\it
i}qr}],$ (51)

where $\Gamma_{r}$ is the strength of the fluctuating current
$\mathcal{J}^{fl}(t)$ and calculated below (see eq.(59)). We note that, as
$\lambda_{q}\geq 0$ for any $q$, the current correlation function is negative
$C^{\mathcal{J}\mathcal{J}}_{r}(t)<0$ for any $t>0$. Moreover, provided that
we first take the infinite-system-size limit $L\rightarrow\infty$, the time-
integrated bond current correlation function
$C^{\mathcal{J}\mathcal{J}}_{0}(t)$ over a large time interval $[-T,T]$ decays
as a function of time $T$ as given below,

$\displaystyle\int\limits_{-T}^{T}C^{\mathcal{J}\mathcal{J}}_{0}(t)\differential{t}\simeq\frac{a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}T^{-\frac{1}{2}}.$ (52)

Finally, by taking the limit $T\rightarrow\infty$, we obtain the following
identity,

$\displaystyle\int\limits_{-\infty}^{\infty}C^{\mathcal{J}\mathcal{J}}_{0}(t)\differential{t}=0;$
(53)

see appendix B.2 for details. Indeed, the above result is a direct consequence
of the negative current correlation present in the system and explains why the
time-integrated bond current fluctuation, as derived in eq.(46), grows
subdiffusively in the initial-time regime, $1\ll t\ll L^{2}$. The asymptotic
form of the time-dependent instantaneous current correlation function, for
$t>0$, in the thermodynamic limit can be written as

$\displaystyle
C^{\mathcal{J}\mathcal{J}}_{0}(t)\simeq-\frac{a(\bar{\rho})}{4\sqrt{\pi
a^{\prime}(\bar{\rho})}}t^{-\frac{3}{2}},$ (54)

where the density dependendent term
$a(\bar{\rho})/\sqrt{a^{\prime}(\bar{\rho})}$ in the prefactor is the same as
that in eq.(46) (see appendix B.2). Again, by employing the dimensional
scaling argument used in the previous section to derive eq.(48), we obtain a
modified power-law decay of the instantaneous current correlation near
criticality,

$\displaystyle
C^{\mathcal{J}\mathcal{J}}_{0}(t)\sim-t^{-\quantity(\frac{3}{2}+\mu)}.$ (55)

Clearly the decay is faster than that away from criticality. As discussed
previously, the faster decay of the near-critical current correlation function
is due to the fact that the activity is very small in the vicinity of
criticality, thus resulting in the anomalous suppression of fluctuations.
Indeed the suppressed fluctuation is charactarized by the exponent $\mu>0$,
whereas $\mu=0$ signifies the subdiffusive growth of the time-integrated
current, expected in a normal diffusive systems.

####  II.4.2 Fluctuating current

Now we discuss the dynamic properties of the fluctuating part
$\mathcal{J}^{(fl)}_{i}(t)$ in the instantaneous bond current, which has
already been defined in eq.(12) and whose strength appears in the actual
current correlation functions (e.g., see eqs.(52) and (50)). Here we derive
the the general space and time dependence of the correlation function
$C^{\mathcal{J}^{(fl)}\mathcal{J}^{(fl)}}_{r}(t,t^{\prime})$ of the
fluctuating current $\mathcal{J}^{(fl)}_{i}(t)$, by using the relation
[obtained from the definition in eq.(12)],

$\displaystyle C^{\mathcal{J}^{(fl)}\mathcal{J}^{(fl)}}_{r}(t,0)=$
$\displaystyle
C^{\mathcal{J}\mathcal{J}}_{r}(t,0)-C^{\mathcal{J}\mathcal{J}^{(d)}}_{r}(t,0)-C^{\mathcal{J}^{(d)}\mathcal{J}}_{r}(t,0)$
$\displaystyle+C^{\mathcal{J}^{(d)}\mathcal{J}^{(d)}}_{r}(t,0),$ (56)

and a second relation

$\displaystyle
C^{\mathcal{J}\mathcal{J}^{(d)}}_{r}(t,0)\equiv\derivative{t}C^{\mathcal{Q}\mathcal{J}^{(d)}}_{r}(t,0)=C^{\mathcal{J}^{(d)}\mathcal{J}^{(d)}}_{r}(t,0).$
(57)

We see that the second and the fourth terms of eq.(II.4.2) cancel each other.
Again, by using the following relation, for $t>t^{\prime}$,

$\displaystyle
C^{\mathcal{J}^{(d)}\mathcal{J}}_{r}(t,t^{\prime})=\derivative{t^{\prime}}C^{\mathcal{J}^{(d)}\mathcal{Q}}_{r}(t,t^{\prime})=\derivative{t^{\prime}}\derivative{t}C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime}),$
(58)

in eq.(II.4.2) along with $t^{\prime}=0$, we finally have the time-
dependendent correlation fuction for the fluctuating current,

$\displaystyle
C^{\mathcal{J}^{\quantity(fl)}\mathcal{J}^{\quantity(fl)}}_{r}(t,t^{\prime}=0)\equiv
C^{\mathcal{J}^{\quantity(fl)}\mathcal{J}^{\quantity(fl)}}_{r}(t)=\delta(t)\Gamma_{r}(\bar{\rho}),$
(59)

where $\Gamma_{r}(\bar{\rho})$ is the density-dependendent strength of the
fluctuating current $\mathcal{J}^{(fl)}$.

Figure 3: Top panel: The scaled variance
$\langle[\mathcal{Q}^{(fl)}_{l}(T)]^{2}\rangle/lT$ of space-time-integrated
fluctuating current $\mathcal{Q}^{(fl)}_{l}(T)$ up to time $T$ [as defined in
eq.(62)], obtained from simulations, is plotted as a function of relative
density $\Delta$ for different subsystem sizes $l=1$ (solid red line), $l=2$
(solid green line), $l=5$ (solid blue line), $l=10$ (solid yellow line); we
have used system size $L=1000$ and final time $T=1000$. The theoretical
prediction as in eq.(63) with $l=1$ (black dashed line) is in an excellent
agreement with the corresponding simulation; also, one can see that the
variance for large $l$ converges quite rapidly to the theoretically predicted
value $2a(\bar{\rho})$ [i.e., eq.(63) for $l=\infty$] (black dot-dashed line).
Bottom panel: A scaling collapse of the scaled variance of space-time-
integrated fluctuating current minus the asymptotic value $2a(\bar{\rho})$ for
different subsystem sizes is observed, when plotted as a function of relative
density $\Delta$, and it is in excellent agreement with theory in eq.(63).

The analytical expression of the strength $\Gamma_{r}$ as given in eq.(44) has
some interesting properties, which are due to the two-particle transfer rule
in the Manna sandpile and are noticeably different from that in the variant of
sandpile with one-particle transfer Jain (2005) and symmetric simple exclusion
processes studied in Ref. Sadhu and Derrida (2016). As in the simple exclusion
processes, the strength $\Gamma_{r}$ for the sandpile with one-particle
transfer rule can be shown to be simply delta correlated in space, i.e.,
$\Gamma_{r}=2a(\bar{\rho})\delta_{r,0}$ ani ; this is because both the models
have a steady-state with a product measure and therefore do not have any
spatial correlations. But, in the case of the Manna sandpile, there are
nonzero spatial correlations, leading to the spatially correlated fluctuating
current, i.e., $\Gamma_{r}\neq 0$ for $r\neq 0$, as shown in eq.(44). Indeed,
as our calculation shows (for details, see appendix sec. C.3), eq.(44) is
exact in the case of Manna sandpile and we have

$\Gamma_{0}=3a(\bar{\rho}),$ (60)

being the strength of the fluctuating current $\mathcal{J}^{(fl)}_{i}$ accross
a single bond $(i,i+1)$. Moreover, we find that there exsists a sum-rule

$\sum\limits_{r}\Gamma_{r}=2a(\bar{\rho}),$ (61)

which, as shown later in eq.(69), is directly related to the scaled space-time
integrated current fluctuations and therefore related to another transport
coefficient, called the mobility, or equivalently, the conductivity, defined
as the ratio between average current and an externally applied small biasing
force Bertini et al. (2015); Chatterjee et al. (2018). It has been derived in
Ref.Chatterjee et al. (2018) that the conductivity in the Manna sandpile is
nothing but the density-dependent activity $a(\bar{\rho})$ itself. Remarkably,
as shown in the next section, here we show that one can indeed relate the
conductivity directly to the current fluctuation in the system.

To verify eq.(59) in simulation, let us first define a cumulative (space-time
integrated) fluctuating current across a subsystem of size $l$ and up to time
$T$,

$\displaystyle\mathcal{Q}_{l}^{(fl)}(T)=\int\limits_{0}^{T}\differential{t}\sum\limits_{i=0}^{l-1}\mathcal{J}_{i}^{(fl)}(t).$
(62)

Then, using eqs. (44) and(59) and after some algebraic manipulations, we
obtain, for $l<L$, a fluctuation relation, which immediately connects the
scaled current fluctuation and the density-dependent activity,

$\displaystyle\frac{1}{lT}\left\langle\quantity(\mathcal{Q}_{l}^{\quantity(fl)}(T))^{2}\right\rangle=2a(\bar{\rho})\quantity(1+\frac{1}{2l}).$
(63)

In top panel of fig.(3), we plot the lhs of eq.(63) as a function of the
relative density $\Delta=\bar{\rho}-\rho_{c}$ for different subsystem sizes
$l=1$ (red line), $l=2$ (green line), $l=5$ (blue line), $l=10$ (yellow line);
we take system size $L=1000$ and final time $T=100$. Note that the variance of
subsystem fluctuating current for subsystem size $l=1$ is actually the
strength $\Gamma_{0}$ of the fluctuating bond current; the corresponding
analytical result $\Gamma_{0}=3a(\bar{\rho})$ (dashed black line) as in
eq.(60) shows an excellent agreement with simulations. For comparison, in the
same fig.3, we also plot $2a(\Delta)$ as a function of $\Delta$ (the black
dash-dotted line), to demonstrate that, as subsystem size $l$ increases, the
scaled variance (lhs of eq.(63)) indeed converges towards $2a(\bar{\rho})$, as
predicted in eq.(63). To show this convergence more quantitatively, in the
bottom panel we plot the scaled quantity
$2l\quantity[\left\langle\quantity(\mathcal{Q}_{l}^{(fl)})^{2}\right\rangle-2a]$
for various subsystem sizes $l=1$ (red line), $l=2$ (green line), $l=5$ (blue
line) and $l=10$ (yellow line); we see that all the curves collapse
excellently onto each other and the collasped master curve match excellently
with the analytically predicted value $2a(\Delta)$ derived in eq.(63).

####  II.4.3 Space-time integrated current

In this section we calculate the steady-state variance
$\left\langle\bar{Q}^{2}(l,T)\right\rangle-\left\langle\bar{Q}(l,T)\right\rangle^{2}$
of the cumulative (space-time integrated) actual particle current
$\bar{Q}(l,T)=\sum_{i=0}^{l-1}\mathcal{Q}_{i}(T)$ across a subsystem of size
$l$ and up to time $T$, which can be written as

$\displaystyle\left\langle\bar{Q}^{2}(l,T)\right\rangle-\left\langle\bar{Q}(l,T)\right\rangle^{2}=\left\langle\bar{Q}^{2}(l,T)\right\rangle$
$\displaystyle=lC^{\mathcal{Q}\mathcal{Q}}_{0}(t,0)+\sum\limits_{r=1}^{l-1}2(l-r)C^{\mathcal{Q}\mathcal{Q}}_{r}(t,0),$
(64)

where we have used the fact that the average steady-state current is zero,
i.e., $\left\langle\bar{Q}(l,T)\right\rangle=0$. Now, by using the following
identity,

$\displaystyle\sum\limits_{r=1}^{l-1}2(l-r)(2-\lambda_{rn})=2\quantity(\frac{\lambda_{ln}-l\lambda_{n}}{\lambda_{n}}),$
(65)

we can rewrite eq.(II.4.3) as

$\displaystyle\left\langle\bar{Q}(l,T)^{2}\right\rangle=2a(\bar{\rho})lT+a(\bar{\rho})T\quantity(1-\delta_{l,L})-2a(\bar{\rho})\frac{a^{\prime}(\bar{\rho})}{L}\sum\limits_{q}\frac{a^{\prime}(\bar{\rho})\lambda_{q}T-1+\exp\quantity(-\lambda_{q}a^{\prime}(\bar{\rho})T)}{\quantity(\lambda_{q}a^{\prime}(\bar{\rho}))^{2}}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})\frac{\lambda_{ql}}{\lambda_{q}}.$
(66)

In fig.4, we plot the subsystem current fluctuation
$\left\langle\bar{Q}^{2}(l,T)\right\rangle$ obtained from simulations as a
function of relative density $\Delta=\rho-\rho_{c}$ for various subsystem
sizes $l$ and final times $T$: $l=2500$, $T=100$ (sky-blue line) and $l=100$,
$T=10^{5}$ (violet line). In the same figure we also compare the simulation
results with theory eq.(66): $l=2500$, $T=100$ (dot-dashed line) and $l=100$,
$T=10^{5}$ (dashed line); we observe excellent agreement between simulations
and theory. Here we note that the results for the larger subsystem size $l$
and smaller $T$ almost coincide with twice of local activity, $2a(\Delta)$, as
a function of $\Delta$ (red dashed-dotted line).

Importantly, the asymptotic expression of the variance of cumulative subsystem
(space-time integrated) current as in eq.(66) depends on the order of limits
of the two variables $T\gg 1$ and $l\gg 1$, i.e.,

$\displaystyle\frac{\left\langle\bar{Q}^{2}(l,T)\right\rangle}{lT}\simeq\begin{cases}\frac{2a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}\frac{l}{\sqrt{T}}&\text{for}\hskip 4.0ptT\gg 1,l\gg
1,\\\
2a(\bar{\rho})-\frac{8a(\bar{\rho})\sqrt{a^{\prime}(\bar{\rho})}}{3\sqrt{\pi}}\frac{\sqrt{T}}{l}&\text{for}\hskip
4.0ptl\gg 1,T\gg 1;\end{cases}$ (67)

see appendix sec. B.3 for details. The first expression in the above equation
have been obtained by taking the limit in the following order, first $T\gg 1$
and then the limit $l\gg 1$. In this particular order of limits, the scaled
fluctuation $\left\langle\bar{Q}^{2}(l,T)\right\rangle/lT$ decreases as
$1/\sqrt{T}$ and eventually vanishes in the limit of $T\rightarrow\infty$. On
the other hand, if we take the limit in the opposite order, $l\gg 1$ first and
then $T\gg 1$, we obtain the second asymptotic expression in eq.(67). That is,
in the limit $l\rightarrow\infty$, the scaled subsystem-current fluctuation
$\left\langle\bar{Q}^{2}(l,T)\right\rangle/lT$ tends to $2a(\bar{\rho})$ as
one increases $T$,

$\displaystyle\sigma_{Q}^{2}(\bar{\rho})\equiv\lim\limits_{l\rightarrow\infty}\lim\limits_{T\rightarrow\infty}\frac{\left\langle\bar{Q}^{2}(l,T)\right\rangle}{lT}=2a(\bar{\rho});$
(68)

here the infinite-subsystem-size limit is taken first, and then the infinite-
time limit. Note that, in all the above cases, we have taken the large system
size limit $L/l\gg 1$ at the very beginning. In fact, one can immediately
identify the rhs of eq.(68) as the mobility (equivalently, the conductivity)
for the Manna sandpile, as calculated in Ref.Chatterjee et al. (2018). Indeed,
eq. (68) can be thought of as a nonequilibrium version of the celebrated
Green-Kubo relations well known for equilibrium systems MARCONI et al. (2008).

Figure 4: Scaled space-time-integrated current fluctuations as a function of
relative density. The simulation data for subsystem size $l=2500$ and $T=100$
is plotted in solid blue line and data for $l=100$ and $T=10^{5}$ is plotted
in solid violet line. We compared the analytical result eq.(66) (corresponding
dashed black line) with the simulation, which is in excellent agreement.
Simulation data have been taken for $L=5000$ in both cases. In the inset we
compared the scaled total current fluctuation and twice of activity
$2a(\bar{\rho})$ as a function of $\Delta$.

Interestingly, if we take $l=L\gg 1$, which corresponds to the bond currents
summed over the whole system, we have the following identity

$\displaystyle\lim\limits_{L\rightarrow\infty}\frac{\left\langle\bar{Q}^{2}(L,T)\right\rangle}{LT}=2a(\rho)=\sum\limits_{r}\Gamma_{r}.$
(69)

Notably the above equality is valid for any finite time $T$. This is because
the sum of the diffusive currents over the full system,
$\sum_{i=1}^{L}\mathcal{J}_{i}^{(d)}$, is zero by definition (see eq.(13)).
Consequently the rhs of eq.(69) is equal to the space-time inetgral of the
fluctuating-current correlation function
$\sum_{r=-\infty}^{\infty}\int_{-\infty}^{\infty}\differential{t}C^{\mathcal{J}^{\quantity(fl)}\mathcal{J}^{\quantity(fl)}}_{r}(t,0)=2a(\bar{\rho})$,
obtained using eqs.(44), (59) and (61). In the inset of fig.4, the scaled
variance $\left\langle\bar{Q}^{2}(L,T)\right\rangle/LT$ (solid red line) and
twice the local activity $2a(\Delta)$ (black dashed line) are plotted as a
function of $\Delta$, which is in excellent agreement with eq.(69). Clearly,
below the critical point $\Delta<0$, the system goes into an absorbing state
and, as a result, the current fluctuation is identically zero. Usually
activity is considered to be the order parameter in the sandpiles. Indeed, as
the identity eq.(69) suggests, the space-time integrated current fluctuation
can serve as an order parameter and thus characterizes the dynamical state of
the system. Later we show that the self-diffusion coefficient of tagged
particles can be expressed in terms of the activity and, as previously noted
in da Cunha et al. (2009), it can be considered an alternative description of
the system’s order parameter.

####  II.4.4 Power spectrum

The two-point time-dependent correlation function for instantaneous bond
current can be characterized also through the power spectrum analysis, which
we perform in this section. From the Wiener-Khinchin theorem MacDonald (2006),
the power spectrum for the instantaneous bond current $\mathcal{J}_{i}(t)$ is
expressed in terms of the Fourier transform of the time-correlation functions,

$\displaystyle
S_{\mathcal{J}}(f)=\int\limits_{-\infty}^{\infty}\differential{t}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)e^{2\pi{\bf\it
i}ft}.$ (70)

Setting $r=0$ in eq.(II.4.1), we perform the integration in the rhs of the
above equation, leading to the following expression,

$S_{\mathcal{J}}(f)=\frac{2a(\bar{\rho})}{L}+\frac{2a(\bar{\rho})}{L}\sum_{q}\quantity(1+\frac{\lambda_{q}}{4})\frac{4\pi^{2}f^{2}}{\lambda_{q}^{2}a^{\prime}(\bar{\rho})^{2}+4\pi^{2}f^{2}}.$
(71)

Now, by subtracting the $f=0$ mode,

$\displaystyle
S_{\mathcal{J}}(0)=\lim\limits_{T\rightarrow\infty}\frac{\left\langle\mathcal{Q}^{2}(T)\right\rangle}{T}=\frac{2a(\bar{\rho})}{L},$
(72)

from the lhs of eq.(71), we rewrite eq. (71) in terms of the modified power
spectrum $\tilde{S}_{\mathcal{J}}(f)=S_{\mathcal{J}}(f)-S_{\mathcal{J}}(0)$,

$\displaystyle\tilde{S}_{\mathcal{J}}(f)=\frac{2a(\bar{\rho})}{L}\sum_{q}\quantity(1+\frac{\lambda_{q}}{4})\frac{4\pi^{2}f^{2}}{\lambda_{q}^{2}a^{\prime}(\bar{\rho})^{2}+4\pi^{2}f^{2}}.$
(73)

We can now straightforwardly obtain the asymptotic form of eq.(73) for small
frequency $1/L^{2}\ll f\ll 1$. To do this, we first replace the sum in eq.(73)
as an integral over the variable $x=q/2\pi$,

$\displaystyle\tilde{S}_{\mathcal{J}}(f)\simeq
4a(\bar{\rho})\int\limits_{1/L}^{1/2}\differential{x}\frac{1+\frac{\lambda(x)}{4}}{1+\frac{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})}{4\pi^{2}f^{2}}},$
(74)

where $\lambda(x)\simeq 4\pi^{2}x^{2}$. Then, by performing the variable
transformation

$\displaystyle
y=\frac{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})}{4\pi^{2}f^{2}},$ (75)

and doing some algebraic manipulations, we immediately obtain the modified
power spectrum of current,

$\displaystyle\tilde{S}_{\mathcal{J}}(f)\simeq
a(\bar{\rho})\sqrt{\frac{f}{2\pi
a^{\prime}(\bar{\rho})}}\int\limits_{0}^{\infty}\differential{y}\frac{y^{-\frac{3}{4}}}{(1+y)}=\frac{\sqrt{\pi}a(\bar{\rho})}{\sqrt{a^{\prime}(\bar{\rho})}}f^{\frac{1}{2}};$
(76)

see appendix sec. A for details. Again, by using the previous dimensional
scaling argument where $a/\sqrt{a^{\prime}}\sim f^{\mu}$ with $\mu$ given in
eq.(48) (see sec.II.3), we obtain the desired scaling behavior of the
subtracted power spectrum near criticality,

$\displaystyle\tilde{S}_{\mathcal{J}}(f)\sim f^{\psi_{\mathcal{J}}},$ (77)

where $\psi_{\mathcal{J}}=1/2+\mu$. Since $\mu>0$, with decreasing frequency,
the near-critical power spectrum in the above equation (77) decays faster than
that away from criticality (given by eq.(76)). In simulations, we calculate
the power spectrum by discretizing time over a small interval $\delta t$ and
calculate the discrete Fourier transform

$\displaystyle\tilde{\mathcal{J}}_{n;T}=\delta
t\sum\limits_{k=0}^{T-1}e^{{\bf\it i}2\pi f_{n}k}\mathcal{J}_{i}(k),$ (78)

where $f_{n}=n/T$ with $T$ being large. Then we define the power spectrum of
the bond current as

$\displaystyle
S_{n}=\lim\limits_{T\rightarrow\infty}\frac{1}{T}\left\langle\absolutevalue{\tilde{\mathcal{J}}_{n;T}}^{2}\right\rangle.$
(79)

In fig.(5), we plot the substracted power spectrum
$\tilde{S}_{\mathcal{J}}(f)$, obtained from simulations, for various densities
$\rho=0.97$ (violet line), $\rho=1.0$ (green line), $\rho=1.2$ (blue line),
$\rho=1.5$ (yellow line), $\rho=2.0$ (red line). For $\rho=2.0$, we also plot
$\tilde{S}_{\mathcal{J}}(f)$ obtained from theory eq.(73) (black dashed line),
which shows an excellent agreement with simulation; top-most guiding line -
$f^{1/2}$ [behavior away from criticality as in eq.(76)] and the bottom-most
guiding line - $f^{1/2+\mu}$ [behavior near criticality as in eq.(77)].

Figure 5:  Power spectrum of instantaneous current is plotted as a function of
frequency $f$ for various densities $\bar{\rho}=0.97$ (magenta), $1.0$
(green), $1.2$ (sky-blue), $1.5$ (yellow), and $2.0$ (red) and for system size
$L=1000$. The top and bottom dashed guiding lines represent the asymptotic
behaviors, the $f^{1/2}$ (far from criticality) and $f^{1/2+\mu}$ (near
criticality) scaling, respectively, with $\mu\approx 0.24$, obtained by using
$\beta\approx 0.42$, $\nu_{\perp}\approx 1.81$ and $z\approx 1.66$ in eq.
(48). The dashed black line represents the theoretical result eq.(71) for
$\rho=2.0$ and is in excellent agreement with the correponding simulation (red
line).

###  II.5 Tagged-particle displacement fluctuation

In this section, we study the fluctuations in tagged particle displacements as
a function of time. We can relate the sum of all individual time-integrated
tagged particle displacement $\sum_{\alpha=1}^{N}X_{\alpha}(T)$, where the net
displacement $X_{\alpha}(T)$ of $\alpha^{th}$ particle in a time interval
$[0,T]$, to the space-time integrated current by the following relation,

$\sum\limits_{\alpha=1}^{N}X_{\alpha}(T)=\sum\limits_{i=0}^{L-1}\mathcal{Q}_{i}(T)=\bar{\mathcal{Q}}(L,T).$
(80)

In the limit of large $T\gg L^{2}$, the self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$ can be defined through the mean-square tagged
particle displacement of the $\alpha^{th}$ particle as given below,

$\left\langle X^{2}_{\alpha}(T)\right\rangle\simeq
2\mathcal{D}_{s}(\bar{\rho})T.$ (81)

To compute the lhs of the above equation, we write the variance of the sum
$\sum_{\alpha=1}^{N}X_{\alpha}(T)$ as

$\displaystyle\left\langle\quantity[\sum\limits_{\alpha=1}^{N}X_{\alpha}(T)]^{2}\right\rangle=$
$\displaystyle\sum\limits_{\alpha,n}\sum\limits_{\alpha^{\prime},{n^{\prime}}}\left\langle\delta
X_{\alpha}(t_{n})\delta X_{\alpha^{\prime}}(t_{n^{\prime}})\right\rangle$
$\displaystyle=$ $\displaystyle 2a(\bar{\rho})LT,$ (82)

where $\delta X_{\alpha}(t)$ is the microscopic displacement of the
$\alpha^{th}$ particle in a small time interval $(t_{n},t_{n}+\delta t)$

$X_{\alpha}(T)=\sum\limits_{n}\delta X_{\alpha}(t_{n}),$ (83)

and we have used eq.(69) in the last line of eq.(II.5). Now using
$\left\langle\delta X_{\alpha}(t)\delta
X_{\alpha^{\prime}}(t^{\prime})\right\rangle\simeq 0$ for $t\neq t^{\prime}$
and therefore $\left\langle
X_{\alpha}^{2}(T)\right\rangle=\sum\limits_{n}\left\langle\delta
X_{\alpha}^{2}(t_{n})\right\rangle$, we get

$\displaystyle\left\langle\quantity[\sum\limits_{\alpha=1}^{N}X_{\alpha}(T)]^{2}\right\rangle\simeq\sum\limits_{\alpha=1}^{N}\left\langle
X_{\alpha}^{2}(T)\right\rangle=N\left\langle X^{2}_{\alpha}(T)\right\rangle.$
(84)

Comparing eqs.(II.5), (81) and (84), we obtain the following relations,

$\mathcal{D}_{s}(\bar{\rho})=\frac{a(\bar{\rho})}{\bar{\rho}}=\frac{1}{\bar{\rho}}\quantity[\lim\limits_{L,T\rightarrow\infty}\frac{\left\langle\bar{\mathcal{Q}}^{2}(L,T)\right\rangle}{2LT}],$
(85)

which connects the self-diffusion coefficient, activity and the space-time
integrated current fluctuation. Alternatively, one can show the above relation
using a slightly different argument as follows. First we note that
$\left\langle X_{\alpha}^{2}(T)\right\rangle=\left\langle
N_{\alpha}^{(h)}(T)\right\rangle$, where $N_{\alpha}^{(h)}(T)$ is the total
number of hops, performed by the $\alpha^{th}$ particle up to time $T$ da
Cunha et al. (2009). Summing over all partciles we obtain

$\displaystyle\sum\limits_{\alpha}\left\langle
X_{\alpha}^{2}(T)\right\rangle=\sum\limits_{\alpha}\left\langle
N_{\alpha}^{(h)}(T)\right\rangle=2\left\langle N^{(tp)}(T)\right\rangle,$ (86)

where $N^{(tp)}(T)$ is the total number of toppling in the whole system up to
time $T$ and we have used the fact that
$\sum\limits_{\alpha}N_{\alpha}^{(h)}(T)=2N^{(tp)}(T)$. Also in the limit of
large $T$, we have in the leading order of $T$,

$\displaystyle\left\langle N^{(tp)}(T)\right\rangle\simeq a(\bar{\rho})TL,$
(87)

where $a(\bar{\rho})$ is the activity at density $\bar{\rho}$. By summing
eq.(81) over all particles,

$\displaystyle\sum_{\alpha}\left\langle X_{\alpha}^{2}(T)\right\rangle\simeq
2\mathcal{D}_{s}(\bar{\rho})TN,$ (88)

and, then by using eqs. (87) and (88) in eq.(86), we obtain the relation as
given in eq.(85).

Figure 6:  In the top panel, we display the typical space-time trajectories
for three tagged particles corresponding to the densities $\rho=0.97$
(violet), $1.5$ (green), $2.0$ (blue) respectively with the scaled time axis.
In the bottom panel, we plot the mean-square fluctuation of tagged particle
displacement up to time $T$ (solid red line) as a function of relative density
$\Delta$, where the double angular braces $\left\langle\left\langle
X^{2}(T)\right\rangle\right\rangle=\sum_{\alpha}\left\langle
X_{\alpha}^{2}(T)\right\rangle/N$ denote the average over trajectories as well
as particles. Simulations (solid red line) show excellent agreement with the
theoretically obtained self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$ (dashed black line) as in eq.(85). In the same
panel, we also plot the bulk-diffusion coefficient
$D(\bar{\rho})=a^{\prime}(\bar{\rho})$ as a function of
$\Delta=\bar{\rho}-\rho_{c}$ (solid blue line), using eq.(7) from simulation,
which has a contrasting behavior as compared to the self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$.

In the top panel of fig. 6, we plot typical trajectories of a particular
tagged particle for different densities, $\rho=0.97$ (violet line), $\rho=1.5$
(green line), $\rho=2.0$ (blue line), as a function of scaled time. In the
bottom panel, we plot the mean square fluctuation of tagged particle
displacement (solid red line) up to time $T$, $\left\langle\left\langle
X^{2}(T)\right\rangle\right\rangle/2T$, as a function of
$\Delta=\rho-\rho_{c}$, where the double angular braces
$\left\langle\left\langle
X^{2}(T)\right\rangle\right\rangle=\sum_{\alpha}\left\langle
X_{\alpha}^{2}(T)\right\rangle/N$ denote the average over trajectories as well
as particles; in simulations, during the particle transfer at any site $i$,
two particles are chosen randomly from a particular stack. Our theoretical
expression of the density dependent self-diffusion coefficient
$\mathcal{D}_{s}(\rho)$, as given in the equality eq.(85), is also plotted
(dashed black line); one can see an excellent agreement between the simulation
and the theoretical prediction. For comparison, we also plot the bulk-
diffusion coefficient $D(\bar{\rho})$ (solid blue line), defined in eq.(7).
Here one should note that the self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$ and the bulk-diffusion coefficient
$D(\bar{\rho})$ are, in principle, two different quantities and strikingly
they have quite contrasting behaviors, especially near criticality. Indeed,
upon approaching criticality where activity decays as
$a(\bar{\rho})\sim(\bar{\rho}-\rho_{c})^{\beta}$ with $\beta<1$, the self-
diffusion coefficient being the ratio of activity to density (see eq.(85))
vanishes as $\mathcal{D}_{s}(\bar{\rho})\sim(\bar{\rho}-\rho_{c})^{\beta}$ \-
exactly in the same manner as the activity behaves near criticality, but the
bulk-diffusion coefficient being derivative of activity wrt the density [see
eq.(7)] diverges as $D(\bar{\rho})\sim 1/(\bar{\rho}-\rho_{c})^{1-\beta}$
Chatterjee et al. (2018); Tapader et al. (2021). Moreover, far from
criticality and in the limit of large density $\bar{\rho}\gg 1$, though the
self-diffusion coefficient and the bulk-diffusion coefficient both vanish,
however they do so in different manners. In that case, as the activity is
expected to behave as $a(\bar{\rho})\simeq 1-\text{const.}/\bar{\rho}$, the
self-diffusivity decays as $\mathcal{D}_{s}(\bar{\rho})\sim 1/\bar{\rho}$, but
the bulk-diffusivity decays much faster, $D(\bar{\rho})\sim 1/\bar{\rho}^{2}$.
Lastly, in the active phase, where $\bar{\rho}>\rho_{c}$, while the bulk-
diffusion coefficient $D(\bar{\rho})$ is a monotonically decreasing function
of density $\bar{\rho}$ shown in fig. 6, the self-diffusion coefficient
$\mathcal{D}_{s}(\bar{\rho})$ is however a non-monotonic function of
$\bar{\rho}$.

Importantly, unlike in the symmetric simple exclusion process where both the
time-integrated bond current and the tagged particle displacement fluctuations
grow sub-diffusively as $T^{1/2}$ Sadhu and Derrida (2016); De Masi and
Ferrari (2002), in the conserved Manna sandpile only the current fluctuation
grows sub-diffusively, whereas the tagged particle displacement fluctuation
always grows linearly with time. This is perhaps not surprising, given the
fact that, in the Manna sandpile, there are no restrictions in the particle
crossings, which are otherwise not allowed in the symmetric exclusion process.

###  II.6 Mass fluctuation and power spectrum

In the previous sections, we studied various properties of current
fluctuations in detail. Similarly, in this section, starting from the
microscopic update rules combined with the previously introduced truncation
scheme, we shall derive various dynamic properties of mass fluctuations. The
basic quantity is the two-point dynamic correlation function
$C^{mm}_{r}(t,t^{\prime})=\left\langle
m_{0}(t)m_{r}(t^{\prime})\right\rangle-\left\langle
m_{0}(t)\right\rangle\left\langle m_{r}(t^{\prime})\right\rangle$. By using
the microscopic update rules, we write the time evolution equation for
$C_{r}^{mm}(t,0)\equiv C_{r}^{mm}(t)$ as

$\displaystyle\derivative{t}C_{r}^{mm}(t)=\sum\limits_{k}\Delta_{0,k}\left\langle\hat{a}_{k}(t)m_{r}(0)\right\rangle.$
(89)

Using the earlier truncation approximation eq.(23), we write the above
equation as

$\displaystyle\derivative{t}C_{r}^{mm}(t)\simeq
a^{\prime}(\bar{\rho})\sum\limits_{k}\Delta_{r,k}C^{mm}_{k}(t).$ (90)

The solution of eq.(90) can be written, by using the Fourier representation,
as

$\displaystyle\tilde{C}_{q}^{mm}(t)\simeq
e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\tilde{C}_{q}^{mm}(0),$ (91)

where $\tilde{C}_{q}^{mm}$ is the Fourier transform of $C_{r}^{mm}$. The
equal-time mass correlation can be solved by using the approximation eq.(23)
in eq.(36) and we can write the time evolution of $C_{r}^{mm}(t,t)$ in the
steady state as

$\displaystyle\derivative{t}C_{r}^{mm}(t,t)\simeq
2a^{\prime}(\bar{\rho})\sum\limits_{k}\Delta_{0,k}\left\langle
m_{k}m_{r}\right\rangle+B_{r}=0.$ (92)

Similar to what was done earlier to solve eq.(36), the above equation can be
solved exactly using a generating function,

$G(z)=\frac{1}{a^{\prime}(\bar{\rho})}\quantity(\frac{3a(\bar{\rho})}{2}-\frac{a(\bar{\rho})}{4}z).$
(93)

According to the above generating function, we have the steady-state
correlations $C_{0}^{mm}=\left\langle
m_{0}^{2}\right\rangle-\bar{\rho}^{2}=3a/2a^{\prime}$,
$C_{1}^{mm}=\left\langle
m_{0}m_{1}\right\rangle-\bar{\rho}^{2}=-a/4a^{\prime}$ and all other
correlations being zero. Thus we immediately arrive at a relation between the
scaled subsystem-mass fluctuation and the activity,

$\sigma^{2}(\rho)\equiv\lim\limits_{l\rightarrow\infty}\frac{\left\langle(\Delta
M_{l})^{2}\right\rangle}{l}=\sum\limits_{r=-\infty}^{r=\infty}C^{mm}_{r}=\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})},$
(94)

where $\Delta M_{l}=M_{l}-\left\langle M_{l}\right\rangle$. Interestingly, by
using eqs. (7) and (68), the above identity can be recast into a
nonequilibrium version of the Green-Kubo-like relation Chatterjee et al.
(2018),

$\sigma^{2}(\bar{\rho})=\frac{\sigma^{2}_{Q}(\bar{\rho})}{2D(\bar{\rho})},$
(95)

connecting the (scaled) subsystem-mass fluctuation $\sigma^{2}(\bar{\rho})$,
the (scaled) subsysten-current fluctuation $\sigma^{2}_{Q}(\bar{\rho})$ and
the bulk-diffusion coefficient $D(\bar{\rho})$ (a slightly different form of
the above relation is usually referred as the Einstein relation in the
literature Bertini et al. (2015); Chatterjee et al. (2018)).

Next we write the solution of eq.(91) using the generating function in eq.(93)
as

$\displaystyle\tilde{C}_{q}^{mm}(t)\simeq
e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\quantity(1+\frac{\lambda_{q}}{4}).$
(96)

Finally, using the inverse Fourier transformation, we get,

$C_{r}^{mm}(t)\simeq\frac{1}{L}\sum\limits_{q}e^{-{\bf\it
i}qr}e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\quantity(1+\frac{\lambda_{q}}{4}).$
(97)

We now consider subsystem mass $M_{l}(t)=\sum_{r=0}^{l-1}m_{r}(t)$ for $l<L$
and calculate the equal-time correlation function for mass
$C^{M_{l}M_{l}}(t,0)\equiv C^{M_{l}M_{l}}(t)$ by using the following
expression,

$\displaystyle
C^{M_{l}M_{l}}(t)=lC_{0}^{mm}(t)+\sum\limits_{r=1}^{l-1}(l-r)\quantity(C_{r}^{mm}(t)+C_{-r}^{mm}(t)).$
(98)

Then by substituting eq.(97) in eq.(98), we get the equal-time correlation for
subsystem mass,

$\displaystyle
C^{M_{l}M_{l}}(t)\simeq\frac{1}{L}\sum\limits_{q}e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\quantity(1+\frac{\lambda_{q}}{4})\frac{\lambda_{lq}}{\lambda_{q}}.$
(99)

For $t=0$, the correlation function $C^{M_{l}M_{l}}(0)$ is nothing but the
equal-time subsystem mass fluctuation, which can be written in the large
system size $L\rightarrow\infty$ limit as given below

$\displaystyle C^{M_{l}M_{l}}(0)=\left\langle(M_{l}-\left\langle
M_{l}\right\rangle)^{2}\right\rangle=\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}l\quantity[1+\frac{1}{2l}].$
(100)

Then, by taking the large subsystem size $l\rightarrow\infty$ limit, where
$1\ll l\ll L$, we recover the Einstein relation, already derived in eq.(94).
Moreover the asymptotic form of eq.(99) can be written as

$\displaystyle
C^{M_{l}M_{l}}(t)-C^{M_{l}M_{l}}(0)\simeq-\frac{2a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}t^{\frac{1}{2}}$ (101)

for large time $1\ll t\ll L^{2}$; see appendix sec. B.4 for details.. Using
the Fourier transform of eq.(99), we write the power spectrum of the subsystem
mass fluctuation as

$\displaystyle
S_{M}(f)=\lim\limits_{T\rightarrow\infty}\int\limits_{-T}^{T}\differential{t}C^{M_{l}M_{l}}(t)e^{2\pi{\bf\it
i}ft},$ (102)

which can be written, by using eq.(99), as

$\displaystyle
S_{M}(f)=\frac{1}{L}\sum\limits_{q}\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\quantity(1+\frac{\lambda_{q}}{4})\frac{2\lambda_{q}a^{\prime}(\bar{\rho})}{\lambda_{q}^{2}a^{\prime}(\bar{\rho})^{2}+4\pi^{2}f^{2}}\frac{\lambda_{lq}}{\lambda_{q}}.$
(103)

In fig. 7, we plot the power spectrum of the subsystem mass fluctuation,
obtained form simulation for $L=1000$ and $l=500$, for various densities. We
compare the analytical expression eq.(103) (dotted black line) with simulation
result for $\bar{\rho}=2.0$ (solid red line), which is in excellent agreement
with theory. The asymptotic expression of the power spectrum in eq.(103) can
be obtained by simplifying the integral as given below

$\displaystyle S_{M}(f)=$ $\displaystyle
4a(\bar{\rho})\int\limits_{1/L}^{1/2}\differential{x}\frac{\lambda(lx)\quantity(1+\lambda(x)/4)}{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})+4\pi^{2}f^{2}}$
$\displaystyle\simeq$
$\displaystyle\frac{a(\bar{\rho})}{2\sqrt{\pi^{3}a^{\prime}(\bar{\rho})}}f^{-\frac{3}{2}}.$
(104)

Here, in the first step, we have replaced the sum in the rhs of eq.(103) as an
integral $(1/L)\sum_{q}\rightarrow\int_{0}^{2\pi}\differential{q}$ in the
limit $L\rightarrow\infty$ and we used $q=2\pi x$, $\lambda(x)=4\pi^{2}x^{2}$
and eq.(75); see appendix sec. B.5 for details. The above asymptotic form of
the power spectrum can be used to calculate the behavior near criticality by
using the dimensional scaling argument as performed before in eq.(47) where we
write $a/\sqrt{a^{\prime}}\sim f^{\mu}$. In other words, near criticality, the
decay of the power spectrum $S_{M}(f)$ as a function of frequency $f$ becomes
slower and is given by

$\displaystyle S_{M}(f)\sim f^{-\psi_{M}},$ (105)

where $\psi_{M}=3/2-\mu$. In fig. 7, we plot $S_{M}(f)$ as a function of
frequency for densities $\bar{\rho}=2.0$ (red), $1.5$ (orange), $1.2$ (blue),
$1.0$ (green), $0.97$ (violet); we observe that, as one approaches
criticality, the decay of the power spectrum indeed becomes slightly slower,
in accordance with our theoretical prediction in eq.(105). The slower decay of
the power spectrum up on approaching criticality can be physically understood
from the current power spectrum as follows. Due to the slower temporal growth
of the time-integrated bond current fluctuation (see eq.(47)), the near-
critical subsystem mass correlation also decays slower as a function of time,
i. e.,

$\displaystyle C^{M_{l}M_{l}}(t)-C^{M_{l}M_{l}}(0)\sim-t^{1/2-\mu},$ (106)

which is due to the fact that the time-integrated current grows slower with
time and consequently the subsystem tends to retain a particular amount of
mass for a much longer period. Indeed, this phenomenon can be thought of as
the hyperuniformity of mass fluctuations in the temporal domain - a dynamic
hyperuniformity of mass fluctuation, analogous to that of current fluctuations
as described previously in eq.(47).

Figure 7:  Power spectrum of subsystem mass fluctuations are plotted for
$L=1000$ and $l=500$. From top to bottom the points represent the data for the
densities $\bar{\rho}=2.0$ (red), $1.5$ (orange), $1.2$ (blue), $1.0$ (green),
$0.97$ (violet) respectively. The upper guiding line represents the $f^{-3/2}$
behavior [eq. (II.6)], whereas the lower guiding line represents
$f^{-\Psi_{M}}$ behavior [eq.(105)] where $\Psi_{M}=3/2-\mu\approx 1.26$. The
dashed black line represents the theoretical result eq.(103) for $\rho=2.0$.

We note that the two exponents $\psi_{\mathcal{J}}$ and $\psi_{M}$ for the
current and mass power spectra defined in eqs.(77) and (105), respectively,
are in fact related, due to the mass conservation as expressed in the
continuity equation (8). By using the Fourier transform of an observable
$A_{r}(t)=\int\limits_{-\infty}^{\infty}\differential{f}\sum\limits_{q}e^{-2\pi{\bf\it
i}ft}e^{-{\bf\it i}qr}\tilde{\mathcal{A}}_{q}(f)$, we can write eq.(8) as
$-2\pi{\bf\it
i}f\tilde{\mathcal{M}}_{q}(f)=\tilde{\mathcal{J}}_{q}(f)\quantity(e^{{\bf\it
i}q}-1)$. On a large scale $q\rightarrow 0$, we have

$S_{M}(f)\sim f^{-2}S_{\mathcal{J}}(f),$ (107)

and therefore, from eqs.(77) and (105), we obtain the scaling relation

$\psi_{\mathcal{J}}=2-\psi_{M}.$ (108)

##  III Summary and conclusions

In this paper, we study the steady-state dynamical properties of current and
mass in the active phase of the one-dimensional conserved Manna sandpile, and
we establish a direct quantitative relationship between the system’s static
and dynamic properties. First starting with a microscopic dynamical
description, we introduce a truncation scheme, that is approximate and is
expected to be valid only for long (hydrodynamic) times, but allows us to
theoretically investigate the time-dependent (two-point, unequal-time)
correlation functions for current and mass, as well as the associated power
spectra. In particular, we find that, in the thermodynamic limit, the two-
point time-dependent correlation function for the (bond) current has a delta
peak at time $t=0$ and, for time $t>0$, the correlation is negative and a
long-ranged one, decaying as $t^{-(3/2+\mu)}$. Far from criticality, we show
that the exponent $\mu=0$, resulting in a subdiffusive, $T^{1/2}$, growth of
the variance of the cumulative (time-integrated) current up to time $T$. This
type of subdiffusive growth of temporal fluctuation, which has previously been
obtained in symmetric simple exclusion processes Sadhu and Derrida (2016), is
somewhat expected for diffusive systems with normal fluctuation properties,
such as sandpiles away from criticality Tapader et al. (2021). However the
scenario changes near the critical point. Indeed, as one approaches
criticality, the activity in the system vanishes, contributing to an anomalous
suppression of the temporal current fluctuations and thus a positive value of
the exponent $\mu>0$, which has been expressed in terms of the standard static
exponents [see eq.(48)]; likewise, near criticality, the power spectrum of
current at low frequency $f$ varies as $f^{1/2+\mu}$. A similar argument can
be made for the temporal subsystem-mass fluctuation, which is induced by the
boundary currents and is also suppressed near criticality because the current
fluctuation is suppressed. The anomalously reduced mass fluctuation is
manifested in the corresponding power spectrum, which, at low frequency and
near criticality, varies as $f^{-3/2+\mu}$ with $\mu>0$; on the other hand,
far from criticality, the exponent $\mu=0$, implying $f^{-3/2}$ power
spectrum, expected in a normal diffusive system. We also derive, within our
theory, a nonequilibrium version of the Green-Kubo-like fluctuation-response
relation [see eq.(95)], or the Einstein relation Bertini et al. (2015);
Chatterjee et al. (2018), which connects dynamic and static fluctuations in
the system.

We finally investigate the mean-square displacement of tagged-particles and
show that the self-diffusion coefficient for an individual tagged particle is
identically equal to the ratio of the activity to density [see first equality
in eq.(85)]. The identity readily explains a previous simulation observation
of Ref. da Cunha et al. (2009) that the self-diffusion coefficient in the
Manna sandpile vanishes in precisely the same fashion as the activity does
upon approaching criticality. Notably, the near-critical behavior of the self-
diffusion coefficient differs markedly from that of the bulk-diffusion
coefficient, which was previously identified in Refs. Chatterjee et al.
(2018); Tapader et al. (2021) as the derivative $a^{\prime}(\bar{\rho})$ of
the activity wrt density $\bar{\rho}$ and clearly diverges near criticality.
Furthermore, while the bulk-diffusion coefficient is a monotonic function of
density, the self-diffusion coefficient is a non-monotonic one. Interestingly,
the self-diffusion coefficient can be related to the current fluctuation in
the system as expressed in the second equality in eq.(85).

Our findings are reminiscent of somewhat similar observations of dynamic
hyperuniformity, where the existence of anti-correlations in the temporal
fluctuations in sandpiles were pointed out Hwa and Kardar (1989); Garcia-
Millan et al. (2018). However, until recently Chatterjee et al. (2018), the
precise relationship between dynamic and static fluctuations, such as that
between mass and current fluctuations, was unknown, and is encoded in the
Green-Kubo-like fluctuation relation as derived here in eq.(95). The
relationship demonstrates that there are indeed two mechanisms responsible for
the vanishing of mass fluctuation near criticality: Physically, the
anomalously suppressed current fluctuation, combined with the diverging bulk-
diffusion coefficient, near criticality result in the vanishing, or
hyperuniform, density fluctuation observed recently in the conserved Manna
sandpiles Hexner and Levine (2015).

Our results are in fact a consequence of the mass conservation, as reflected
in the scaling relation (108), and are thus expected to be applicable in a
broad class of conserved sandpiles. As our analysis suggests, the anomalous
suppression of current fluctuations near criticality could be a generic
feature of hyperuniform state of matter and should serve as the dynamical
signature of such states, which have been observed in similar other systems in
the recent past Corte et al. (2008); Hexner et al. (2017). In particular, our
findings could help in determining the precise dynamical nature of the off-
and near-critical states in sandpiles by shedding light on the microscopic
dynamical origin of long-ranged temporal correlations in these systems.

###### Acknowledgements.

We thank Deepak Dhar for useful discussions and comments on the manuscript.
A.M. acknowledges financial support from the the Department of Science and
Technology, India [Fellowship No. DST/INSPIRE Fellowship/2017/IF170275] for
part of the work carried out under his senior research fellowship. We
acknowledge the Thematic Unit of Excellence on Computational Materials
Science, funded by the Department of Science and Technology, India for
computational facility.

## References

  * Mandelbrot (1983) B. B. Mandelbrot, _The fractal geometry of nature_ (W. H. Freeman and Comp., 1983), 3rd ed.
  * Frontera and Fuligni (1979) F. Frontera and F. Fuligni, The Astrophysical Journal 232, 590 (1979), ISSN 1538-4357, URL http://dx.doi.org/10.1086/157317.
  * Aschwanden and et al (2014) M. J. Aschwanden and et al, Space Science Reviews 198, 47–166 (2014), ISSN 1572-9672, URL http://dx.doi.org/10.1007/s11214-014-0054-6.
  * Palmieri and Jensen (2020) L. Palmieri and H. J. Jensen, Frontiers in Physics 8 (2020), ISSN 2296-424X, URL http://dx.doi.org/10.3389/fphy.2020.00257.
  * Novikov et al. (1997) E. Novikov, A. Novikov, D. Shannahoff-Khalsa, B. Schwartz, and J. Wright, Physical Review E 56, R2387–R2389 (1997), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.56.R2387.
  * Hesse and Gross (2014) J. Hesse and T. Gross, Frontiers in Systems Neuroscience 8 (2014), ISSN 1662-5137, URL http://dx.doi.org/10.3389/fnsys.2014.00166.
  * Tebaldi (2021) C. Tebaldi, Frontiers in Physics 8 (2021), ISSN 2296-424X, URL http://dx.doi.org/10.3389/fphy.2020.616408.
  * Voss and Clarke (1976) R. F. Voss and J. Clarke, Physical Review B 13, 556–573 (1976), ISSN 0556-2805, URL http://dx.doi.org/10.1103/PhysRevB.13.556.
  * Jensen (1998) H. J. Jensen, _Self-Organized Criticality_ (Cambridge University Press, 1998), ISBN 9780511622717, URL http://dx.doi.org/10.1017/CBO9780511622717.
  * Bak (1997) P. Bak, _How nature works: the science of self-organized criticality_ (Oxford University Press, 1997).
  * Bak et al. (1987) P. Bak, C. Tang, and K. Wiesenfeld, Physical Review Letters 59, 381–384 (1987), ISSN 0031-9007, URL http://dx.doi.org/10.1103/PhysRevLett.59.381.
  * Bak et al. (1988) P. Bak, C. Tang, and K. Wiesenfeld, Physical Review A 38, 364–374 (1988), ISSN 0556-2791, URL http://dx.doi.org/10.1103/PhysRevA.38.364.
  * Dhar and Ramaswamy (1989) D. Dhar and R. Ramaswamy, Physical Review Letters 63, 1659–1662 (1989), ISSN 0031-9007, URL http://dx.doi.org/10.1103/PhysRevLett.63.1659.
  * Dhar (1990) D. Dhar, Physical Review Letters 64, 1613–1616 (1990), ISSN 0031-9007, URL http://dx.doi.org/10.1103/PhysRevLett.64.1613.
  * Priezzhev et al. (1996) V. B. Priezzhev, D. Dhar, A. Dhar, and S. Krishnamurthy, Physical Review Letters 77, 5079–5082 (1996), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.77.5079.
  * Majumdar and Dhar (1992) S. Majumdar and D. Dhar, Physica A: Statistical Mechanics and its Applications 185, 129–145 (1992), ISSN 0378-4371, URL http://dx.doi.org/10.1016/0378-4371(92)90447-X.
  * Fey et al. (2010) A. Fey, L. Levine, and D. B. Wilson, Physical Review Letters 104 (2010), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.104.145703.
  * Podder and Rolla (2021) M. Podder and L. T. Rolla, Journal of Statistical Physics 182 (2021), ISSN 1572-9613, URL http://dx.doi.org/10.1007/s10955-021-02731-3.
  * Cipriani et al. (2017) A. Cipriani, R. S. Hazra, and W. M. Ruszel, Probability Theory and Related Fields 172, 829–868 (2017), ISSN 1432-2064, URL http://dx.doi.org/10.1007/s00440-017-0821-x.
  * Dickman et al. (1998) R. Dickman, A. Vespignani, and S. Zapperi, Physical Review E 57, 5095–5105 (1998), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.57.5095.
  * Dickman et al. (2001) R. Dickman, M. Alava, M. A. Muñoz, J. Peltola, A. Vespignani, and S. Zapperi, Physical Review E 64 (2001), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.64.056104.
  * Vespignani et al. (2000) A. Vespignani, R. Dickman, M. A. Muñoz, and S. Zapperi, Physical Review E 62, 4564–4582 (2000), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.62.4564.
  * Vespignani et al. (1998) A. Vespignani, R. Dickman, M. A. Muñoz, and S. Zapperi, Physical Review Letters 81, 5676–5679 (1998), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.81.5676.
  * Le Doussal and Wiese (2015) P. Le Doussal and K. J. Wiese, Physical Review Letters 114 (2015), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.114.110601.
  * Ramasco et al. (2004) J. J. Ramasco, M. A. Muñoz, and C. A. da Silva Santos, Physical Review E 69 (2004), ISSN 1550-2376, URL http://dx.doi.org/10.1103/PhysRevE.69.045105.
  * Dickman (2002) R. Dickman, Physical Review E 66 (2002), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.66.036122.
  * da Cunha et al. (2009) S. D. da Cunha, R. R. Vidigal, L. R. da Silva, and R. Dickman, The European Physical Journal B 72, 441–449 (2009), ISSN 1434-6036, URL http://dx.doi.org/10.1140/epjb/e2009-00367-0.
  * Dhar (2006) D. Dhar, Physica A: Statistical Mechanics and its Applications 369, 29–70 (2006), ISSN 0378-4371, URL http://dx.doi.org/10.1016/j.physa.2006.04.004.
  * Dickman et al. (2000) R. Dickman, M. A. Muñoz, A. Vespignani, and S. Zapperi, Brazilian Journal of Physics 30, 27–41 (2000), ISSN 0103-9733, URL http://dx.doi.org/10.1590/S0103-97332000000100004.
  * Redig (2006) F. Redig, Les Houches, Session LXXXIII 2005 pp. 657–729 (2006).
  * Pradhan (2021) P. Pradhan, Frontiers in Physics 9 (2021), ISSN 2296-424X, URL http://dx.doi.org/10.3389/fphy.2021.641233.
  * Ali (1995) A. A. Ali, Physical Review E 52, R4595–R4598 (1995), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.52.R4595.
  * Hexner and Levine (2015) D. Hexner and D. Levine, Physical Review Letters 114 (2015), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.114.110602.
  * Hexner and Levine (2017) D. Hexner and D. Levine, Physical Review Letters 118 (2017), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.118.020601.
  * Torquato and Stillinger (2003) S. Torquato and F. H. Stillinger, Physical Review E 68 (2003), ISSN 1095-3787, URL http://dx.doi.org/10.1103/PhysRevE.68.041113.
  * Torquato (2016) S. Torquato, Physical Review E 94 (2016), ISSN 2470-0053, URL http://dx.doi.org/10.1103/PhysRevE.94.022122.
  * Nizam et al. (2021) U. S. Nizam, G. Makey, M. Barbier, S. S. Kahraman, E. Demir, E. E. Shafigh, S. Galioglu, D. Vahabli, S. Husnugil, M. H. Guneş, et al., Journal of Physics: Condensed Matter 33, 304002 (2021), ISSN 1361-648X, URL http://dx.doi.org/10.1088/1361-648x/abf9b8.
  * Dandekar and Dhar (2013) R. Dandekar and D. Dhar, EPL (Europhysics Letters) 104, 26003 (2013), ISSN 1286-4854, URL http://dx.doi.org/10.1209/0295-5075/104/26003.
  * Dandekar (2020) R. Dandekar, Europhysics Letters 132, 10008 (2020), ISSN 1286-4854, URL http://dx.doi.org/10.1209/0295-5075/132/10008.
  * Manna (1991) S. S. Manna, Journal of Physics A: Mathematical and General 24, L363–L369 (1991), ISSN 1361-6447, URL http://dx.doi.org/10.1088/0305-4470/24/7/009.
  * Chatterjee et al. (2018) S. Chatterjee, A. Das, and P. Pradhan, Physical Review E 97 (2018), ISSN 2470-0053, URL http://dx.doi.org/10.1103/PhysRevE.97.062142.
  * Tapader et al. (2021) D. Tapader, P. Pradhan, and D. Dhar, Physical Review E 103 (2021), ISSN 2470-0053, URL http://dx.doi.org/10.1103/PhysRevE.103.032122.
  * Kertesz and Kiss (1990) J. Kertesz and L. B. Kiss, Journal of Physics A: Mathematical and General 23, L433–L440 (1990), ISSN 1361-6447, URL http://dx.doi.org/10.1088/0305-4470/23/9/006.
  * Laurson et al. (2005) L. Laurson, M. J. Alava, and S. Zapperi, Journal of Statistical Mechanics: Theory and Experiment 2005, L11001–L11001 (2005), ISSN 1742-5468, URL http://dx.doi.org/10.1088/1742-5468/2005/11/L11001.
  * Jensen et al. (1989) H. J. Jensen, K. Christensen, and H. C. Fogedby, Physical Review B 40, 7425–7427 (1989), ISSN 0163-1829, URL http://dx.doi.org/10.1103/PhysRevB.40.7425.
  * Garcia-Millan et al. (2018) R. Garcia-Millan, G. Pruessner, L. Pickering, and K. Christensen, EPL (Europhysics Letters) 122, 50003 (2018), ISSN 1286-4854, URL http://dx.doi.org/10.1209/0295-5075/122/50003.
  * Yadav et al. (2012) A. C. Yadav, R. Ramaswamy, and D. Dhar, Physical Review E 85 (2012), ISSN 1550-2376, URL http://dx.doi.org/10.1103/PhysRevE.85.061114.
  * Maslov et al. (1999) S. Maslov, C. Tang, and Y.-C. Zhang, Physical Review Letters 83, 2449–2452 (1999), ISSN 1079-7114, URL http://dx.doi.org/10.1103/PhysRevLett.83.2449.
  * Baiesi and Maes (2006) M. Baiesi and C. Maes, Europhysics Letters (EPL) 75, 413–419 (2006), ISSN 1286-4854, URL http://dx.doi.org/10.1209/epl/i2006-10137-2.
  * Giometto and Jensen (2012) A. Giometto and H. J. Jensen, Physical Review E 85 (2012), ISSN 1550-2376, URL http://dx.doi.org/10.1103/PhysRevE.85.011128.
  * Bertini et al. (2015) L. Bertini, A. De Sole, D. Gabrielli, G. Jona-Lasinio, and C. Landim, Reviews of Modern Physics 87, 593–636 (2015), ISSN 1539-0756, URL http://dx.doi.org/10.1103/RevModPhys.87.593.
  * Sadhu and Derrida (2016) T. Sadhu and B. Derrida, Journal of Statistical Mechanics: Theory and Experiment 2016, 113202 (2016), ISSN 1742-5468, URL http://dx.doi.org/10.1088/1742-5468/2016/11/113202.
  * Sadhu and Dhar (2009) T. Sadhu and D. Dhar, Journal of Statistical Physics 134, 427–441 (2009), ISSN 1572-9613, URL http://dx.doi.org/10.1007/s10955-009-9683-7.
  * De Masi and Ferrari (2002) A. De Masi and P. A. Ferrari, Journal of Statistical Physics 107, 677–683 (2002), ISSN 0022-4715, URL http://dx.doi.org/10.1023/A:1014577928229.
  * Jain (2005) K. Jain, Physical Review E 72 (2005), ISSN 1550-2376, URL http://dx.doi.org/10.1103/PhysRevE.72.017105.
  * (56) The details of the results will be presented elsewhere (????).
  * MARCONI et al. (2008) U. MARCONI, A. PUGLISI, L. RONDONI, and A. VULPIANI, Physics Reports 461, 111–195 (2008), ISSN 0370-1573, URL http://dx.doi.org/10.1016/j.physrep.2008.02.002.
  * MacDonald (2006) D. K. C. MacDonald, _Noise and Fluctuations_ , Dover Books on Physics (Dover Publications, Mineola, NY, 2006).
  * Hwa and Kardar (1989) T. Hwa and M. Kardar, Physical Review Letters 62, 1813–1816 (1989), ISSN 0031-9007, URL http://dx.doi.org/10.1103/PhysRevLett.62.1813.
  * Corte et al. (2008) L. Corte, P. M. Chaikin, J. P. Gollub, and D. J. Pine, Nature Physics 4, 420–424 (2008), ISSN 1745-2481, URL http://dx.doi.org/10.1038/nphys891.
  * Hexner et al. (2017) D. Hexner, P. M. Chaikin, and D. Levine, Proceedings of the National Academy of Sciences 114, 4294–4299 (2017), ISSN 1091-6490, URL http://dx.doi.org/10.1073/pnas.1619260114.
  * Gasper and Rahman (2004) G. Gasper and M. Rahman, _Basic Hypergeometric Series_ (Cambridge University Press, 2004), ISBN 9780511526251, URL http://dx.doi.org/10.1017/CBO9780511526251.

##  Appendix A Some algebraic identities and special integrals

We can deduce several algebric properties of $\lambda_{n}$, which are the
following,

$\displaystyle\sum\limits_{n=1}^{L-1}\frac{\lambda_{nl}}{\lambda_{n}^{2}}=\frac{1}{12}l(l-L)\left(l^{2}-lL-2\right);$
(109)

$\displaystyle\sum\limits_{n=1}^{L-1}\frac{1}{\lambda_{n}}=\frac{L^{2}-1}{12};$
(110)

$\displaystyle\sum\limits_{n=1}^{L-1}\frac{\lambda_{nr}}{\lambda_{n}}=r(L-r);$
(111)

$\displaystyle\sum\limits_{r=1}^{l-1}2(l-r)(2-\lambda_{rn})=2\quantity(\frac{\lambda_{ln}-l\lambda_{n}}{\lambda_{n}}).$
(112)

$\displaystyle\sum\limits_{n=1}^{L-1}\lambda_{nl}=2L,\hskip
4.0pt\text{for}\hskip 4.0ptl=1,2,\ldots.$ (113)

Eq.(110) is a special case of eq.(109) for $l=1$.

The integrals appeared in the context of asymptotic analysis, i.e., in eq.(76)
and later in eq.(139) have very generic solutions in terms of hypergeometric
functions. Generically, we can write those integrals in the following form,

$\displaystyle
I(y)=\int\limits_{0}^{y}\differential{z}\frac{z^{-k}}{1+z}=\frac{y^{1-k}\,_{2}F_{1}(1,1-k;2-k;-y)}{1-k},$
(114)

where ${}_{2}F_{1}$ is the hypergeometric function Gasper and Rahman (2004),
defined as

$\displaystyle_{2}F_{1}(a,b;c;z)=\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int\limits_{0}^{1}\differential{t}t^{b-1}(1-t)^{c-b-1}(1-tz)^{-a}.$
(115)

In the limit of ${y\rightarrow\infty}$, we have,

$\displaystyle\lim\limits_{y\rightarrow\infty}I(y)=\sqrt{2}\pi\hskip
5.0pt\text{for}\hskip 5.0ptk=\frac{1}{4},\frac{3}{4}.$ (116)

We have used the above result in eqs. (76), (139) to obtain the asymptotic for
of the power spectrum of current and subsystem mass respectively.

##  Appendix B Asymptotic analysis

In this section, we provide the calculation details of the results eq.(46),
(52), (54), (67), (101), presented in the main text.

###  B.1 Time-integrated bond current fluctuation

We now derive the asymptotic approximation of the time-integrated bond current
correlation using eq.(46). We can write the unequal-time time-integrated bond
current correlation as

$\displaystyle C^{\mathcal{Q}\mathcal{Q}}_{0}(t,t^{\prime})=$
$\displaystyle\frac{2a(\bar{\rho})}{L}t^{\prime}+2a(\bar{\rho})a^{\prime}(\bar{\rho})\frac{1}{L}\sum\limits_{q}\frac{1-e^{-\lambda_{q}a^{\prime}(\bar{\rho})t^{\prime}}}{\lambda_{q}^{2}{a^{\prime}}^{2}}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})-$
$\displaystyle
a(\bar{\rho})a^{\prime}(\bar{\rho})\frac{1}{L}\sum\limits_{q}\frac{1-e^{-\lambda_{q}a^{\prime}(\bar{\rho})t^{\prime}}+e^{-\lambda_{q}a^{\prime}(\bar{\rho})t}-e^{-\lambda_{q}a^{\prime}(\bar{\rho})(t-t^{\prime})}}{\lambda_{q}^{2}{a^{\prime}}^{2}}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4})$
(117)

In the infinite system size limit $L\rightarrow\infty$, we write the above sum
in the following integral form,

$\displaystyle C^{\mathcal{Q}\mathcal{Q}}_{0}(t,t^{\prime})\simeq$
$\displaystyle
4a^{\prime}(\bar{\rho})a(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}\frac{1-e^{-\lambda(x)a^{\prime}(\bar{\rho})t^{\prime}}}{\lambda(x)^{2}a^{\prime}(\bar{\rho})^{2}}\lambda(x)\quantity(1+\frac{\lambda(x)}{4})$
$\displaystyle-2a^{\prime}(\bar{\rho})a(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}\frac{1-e^{-\lambda(x)a^{\prime}(\bar{\rho})t^{\prime}}+e^{-\lambda(x)a^{\prime}(\bar{\rho})t}-e^{-\lambda(x)a^{\prime}(\bar{\rho})(t-t^{\prime})}}{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})}\lambda(x)\quantity(1+\frac{\lambda(x)}{4}).$
(118)

Note that the integral in the above equation can be expressed in terms of an
integral of the form as given below,

$\displaystyle
a(\bar{\rho})a^{\prime}(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}\frac{1-e^{-\lambda(x)a^{\prime}(\bar{\rho})t^{\prime}}}{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})}\lambda(x)\quantity(1+\frac{\lambda(x)}{4})\simeq\frac{a(\bar{\rho})\sqrt{t^{\prime}}}{2\sqrt{\pi
a^{\prime}(\bar{\rho})}},$ (119)

where, for $t\gg 1$, we have defined
$x=\quantity({y}/{4\pi^{2}a^{\prime}(\bar{\rho})t^{\prime}})^{1/2}$ and used
$\lambda(x)\simeq 4\pi^{2}x^{2}$ and
$\int\limits_{0}^{\infty}\differential{y}y^{-3/2}(1-e^{-y})=2\sqrt{\pi}$, to
explicitly calculate the integral. Using eq.(119) in each of the relevant
terms of the rhs in eq.(B.1) and then after some straightforward algebraic
manipulations, we obtain the following asymptotic form of the time-dependent
integrated bond-current correlation,

$\displaystyle
C^{\mathcal{Q}\mathcal{Q}}_{0}(t,t^{\prime})\simeq\frac{a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}\quantity(\sqrt{t}+\sqrt{t^{\prime}}-\sqrt{\absolutevalue{t-t^{\prime}}}).$
(120)

Now, by putting $t^{\prime}=t\equiv T$, the above asymptotic leads to the
first part (i.e., corresponding to the limit $1\ll T\ll L^{2}$) of eq.(46) in
the main text.

###  B.2 Time-dependent instantaneous current correlation

The steady state unequal-time correlation of instantaneous bond current
$C^{\mathcal{J}\mathcal{J}}_{0}(t,0)=\left\langle\mathcal{J}_{0}(0)\mathcal{J}_{r}(t)\right\rangle-\left\langle\mathcal{J}_{0}(0)\right\rangle\left\langle\mathcal{J}_{r}(0)\right\rangle$,
for $t\geq 0$, is given by the following expression as derived in main text
(see eq.(II.4.1))

$\displaystyle
C^{\mathcal{J}\mathcal{J}}_{0}(t,0)=\delta(t)3a(\bar{\rho})-a^{\prime}(\bar{\rho})a(\bar{\rho})\frac{1}{L}\sum\limits_{q}e^{-a^{\prime}(\bar{\rho})\lambda_{q}t}\lambda_{q}\quantity(1+\frac{\lambda_{q}}{4}).$
(121)

First, we perform the time integral in a finite time domain $[-T,T]$ as given
below,

$\displaystyle\int\limits_{-T}^{T}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)\differential{t}=\Gamma_{0}(\bar{\rho})-\frac{2a(\bar{\rho})}{L}\sum\limits_{q}\quantity(1+\frac{\lambda_{q}}{4})+2a(\bar{\rho})\quantity[\frac{1}{L}\sum\limits_{q}e^{-\lambda_{q}a^{\prime}(\bar{\rho})T}\quantity(1+\frac{\lambda_{q}}{4})],$
(122)

which, using the relations eq.(113) and
$\Gamma_{0}(\bar{\rho})=3a(\bar{\rho})$ as in the main text in eq.(60), we
simplify the above sum as

$\displaystyle\int\limits_{-T}^{T}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)\differential{t}=\frac{2a(\bar{\rho})}{L}+2a(\bar{\rho})\quantity[\frac{1}{L}\sum\limits_{q}e^{-\lambda_{q}a^{\prime}(\bar{\rho})T}\quantity(1+\frac{\lambda_{q}}{4})].$
(123)

Now, first taking the infinite-system size limit, i.e., the limit
$L\rightarrow\infty$, we can further write the sum as an integral,

$\displaystyle\lim\limits_{L\rightarrow\infty}\int\limits_{-T}^{T}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)\differential{t}=4a(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}e^{-\lambda(x)a^{\prime}(\bar{\rho})T}\quantity(1+\frac{\lambda(x)}{4}),$
(124)

where $q=2\pi x$. Now using $\lambda(x)\simeq 4\pi^{2}x^{2}$,
$\int_{0}^{\infty}\differential{y}e^{-y}y^{-1/2}=\sqrt{\pi}$ and a variable
transformation $x=y^{1/2}/2\pi\sqrt{a^{\prime}(\bar{\rho})T}$, we can
explicitly calculate the integral as in eq.(124) as

$\displaystyle\lim\limits_{L\rightarrow\infty}\int\limits_{-T}^{T}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)\differential{t}\simeq\frac{a(\bar{\rho})}{\pi\sqrt{a^{\prime}(\bar{\rho})T}}\int\limits_{0}^{\infty}\differential{y}e^{-y}y^{-1/2}=\frac{a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}T^{-1/2},$ (125)

which is the result in the main text in eq.(52). Finally, by taking the limit
$T\rightarrow\infty$, we get

$\displaystyle\int\limits_{-\infty}^{\infty}C^{\mathcal{J}\mathcal{J}}_{0}(t,0)\differential{t}=0,$
(126)

which is the result presented in main text in eq.(53).

Similarly, we can find the asymptotic form of
$c^{\mathcal{J}\mathcal{J}}_{0}(t,0)$ presented in main text in eq.(54). In
the limit $L\rightarrow\infty$, the temporal current correlation
$c^{\mathcal{J}\mathcal{J}}_{0}(t,0)$ for $t>0$ can be written as an integral,
as given below

$\displaystyle
c^{\mathcal{J}\mathcal{J}}_{0}(t,0)\simeq-2a^{\prime}(\bar{\rho})a(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}e^{-a^{\prime}(\bar{\rho})\lambda(x)t}\lambda(x)\quantity(1+\frac{\lambda(x)}{4}).$
(127)

Now, again using $\lambda(x)\simeq 4\pi^{2}x^{2}$ and a variable
transformation $x=y^{1/2}/2\pi\sqrt{a^{\prime}(\bar{\rho})T}$, we calculate
the above integral as

$\displaystyle
c^{\mathcal{J}\mathcal{J}}_{0}(t,0)\simeq-\frac{a(\bar{\rho})}{\sqrt{4\pi^{2}a^{\prime}(\bar{\rho})}}t^{-3/2}\int\limits_{0}^{\infty}\differential{y}\sqrt{y}e^{-y},$
(128)

where we have ignored the subleading term $\order{t^{-5/2}}$. Finally, using
$\int\limits_{0}^{\infty}\differential{y}\sqrt{y}e^{-y}=\sqrt{\pi}/2$, we get
the result presented in main text in eq.(54),

$\displaystyle
c^{\mathcal{J}\mathcal{J}}_{0}(t,0)\simeq-\frac{a(\bar{\rho})}{4\sqrt{\pi
a^{\prime}(\bar{\rho})}}t^{-3/2}.$ (129)

###  B.3 Spacetime-integrted current fluctuation

Here we derive the asymptotic dependence of eq.(66) on subsystem size $l$ and
time $T$ (see eq.(67)); first by taking the limit $T\gg 1$ and $l\gg 1$ and
then followd by the reverse order of limit $l\gg 1$ and $T\gg 1$. In both
cases, we take infinite system size limit $L\rightarrow\infty$, therefore
$l/L\rightarrow 0$ is always satisfy.

####  B.3.1 Case I: $T\gg 1$, $l\gg 1$

In this case, we write eq.(66) in the following simplified form,

$\displaystyle\left\langle\bar{Q}^{2}(l,T)\right\rangle=\frac{2a(\bar{\rho})Tl^{2}}{L}+2a(\bar{\rho})a^{\prime}(\bar{\rho})\frac{1}{L}\sum\limits_{q}\frac{1-e^{-a^{\prime}(\bar{\rho})\lambda_{q}T}}{\lambda^{2}_{q}{a^{\prime}}^{2}}\quantity(1+\frac{\lambda_{q}}{4})\lambda_{ql}.$
(130)

In the limit $L\rightarrow\infty$, the sum in the above equation can be
converted in to the following integral,

$\displaystyle\left\langle\bar{Q}^{2}(l,T)\right\rangle\simeq
4a(\bar{\rho})a^{\prime}(\bar{\rho})\int\limits_{0}^{1/2}\differential{x}\frac{1-e^{-a^{\prime}(\bar{\rho})\lambda(x)T}}{\lambda^{2}(x){a^{\prime}}^{2}}\quantity(1+\frac{\lambda(x)}{4})\lambda(lx).$
(131)

Using the approximation $\lambda(lx)\simeq 4\pi^{2}l^{2}x^{2}$ for finite
subsystem size $l$ and a variable transformation
$x=y^{1/2}/2\pi\sqrt{a^{\prime}(\bar{\rho})T}$, we get

$\displaystyle\frac{1}{lT}\left\langle\bar{Q}^{2}(l,T)\right\rangle\simeq\frac{a(\bar{\rho})}{\pi\sqrt{a^{\prime}(\bar{\rho})}}\frac{l}{\sqrt{T}}\int\limits_{0}^{\infty}\differential{y}(1-e^{-y})y^{-3/2}=\frac{2a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}\frac{l}{\sqrt{T}},$ (132)

where we use
$\int_{0}^{\infty}\differential{y}(1-e^{-y})y^{-3/2}=2\sqrt{\pi}$. This result
appeared in eq.(67) of the main-text.

####  B.3.2 Case II: $l\gg 1$, $T\gg 1$

To compute the asymptotic form in this limit, we use the approximation
$\lambda(lx)\simeq 2$ to write eq.(66) in the following form,

$\displaystyle\frac{1}{lT}\left\langle\bar{Q}^{2}(l,T)\right\rangle_{c}\simeq
2a(\bar{\rho})+\frac{a(\bar{\rho})}{l}-\frac{8a(\bar{\rho})a^{\prime}(\bar{\rho})}{lT}\int\limits_{0}^{1/2}\frac{a^{\prime}(\bar{\rho})\lambda(x)T-1+\exp\quantity(-\lambda(x)a^{\prime}(\bar{\rho})T)}{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})}\quantity(1+\frac{\lambda(x)}{4}).$
(133)

Again using the variable transform
$x=y^{1/2}/2\pi\sqrt{a^{\prime}(\bar{\rho})T}$ and
$\int_{0}^{\infty}\differential{y}\quantity(y-1+e^{-y})y^{-5/2}=4\sqrt{\pi}/3$,
we get

$\displaystyle\frac{1}{lT}\left\langle\bar{Q}^{2}(l,T)\right\rangle_{c}\simeq$
$\displaystyle
2a(\bar{\rho})+\frac{a(\bar{\rho})}{l}-\frac{2a(\bar{\rho})a^{\prime}(\bar{\rho})}{lT}\int\limits_{0}^{\infty}\differential{y}\quantity(y-1+e^{-y})y^{-5/2}\frac{T^{3/2}}{\pi\sqrt{a^{\prime}(\bar{\rho})}}$
$\displaystyle=$ $\displaystyle
2a(\bar{\rho})+\frac{a(\bar{\rho})}{l}-\frac{8a(\bar{\rho})}{3}\sqrt{\frac{a^{\prime}(\bar{\rho})}{\pi}}\frac{\sqrt{T}}{l},$
(134)

which also appeared in eq.(67).

###  B.4 Temporal-correlation of subsystem mass

The asymptotic form of subsystem mass temporal correlation
$C^{M_{l}M_{l}}(t,0)$, appeared in eq.(101) is derived here. At $t=0$,
$C^{M_{l}M_{l}}(t,0)$ is maximum and after that, it decays as a function of
time $t$. So, in order to extract the temporal dependence of
$C^{M_{l}M_{l}}(t,0)$, we write eq.(99) as

$\displaystyle
C^{M_{l}M_{l}}(0,0)-C^{M_{l}M_{l}}(t,0)=\frac{1}{L}\sum\limits_{q}\quantity(1-e^{-\lambda_{q}a^{\prime}(\bar{\rho})t})\quantity(1+\frac{\lambda_{q}}{4})\frac{\lambda_{lq}}{\lambda_{q}},$
(135)

which, in the limit $L\rightarrow\infty$, $l\gg 1$, can be written as

$\displaystyle C^{M_{l}M_{l}}(0,0)-C^{M_{l}M_{l}}(t,0)\simeq
4\int\limits_{0}^{1/2}\differential{x}\frac{a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\quantity(1-e^{-\lambda(x)a^{\prime}t})\quantity(1+\frac{\lambda(x)}{4})\frac{1}{\lambda(x)}.$
(136)

The above equation can be further simplified using the approximation and exact
results of the previous sec. B.3.1 and we get

$\displaystyle
C^{M_{l}M_{l}}(0,0)-C^{M_{l}M_{l}}(t,0)=\frac{a(\bar{\rho})}{\sqrt{a^{\prime}(\bar{\rho})\pi^{2}}}t^{\frac{1}{2}}\int_{0}^{\infty}\differential{y}y^{-\frac{3}{2}}\quantity(1-e^{-y})=\frac{2a(\bar{\rho})}{\sqrt{\pi
a^{\prime}(\bar{\rho})}}t^{\frac{1}{2}}.$ (137)

Thus we derived eq.(101) of the main-text.

###  B.5 Power spectrum of subsystem mass fluctuation

In order to find the asymptotic form of the power spectrum of subsystem mass
fluctuation $S_{M_{l}}(f)$, in the limit $L\rightarrow\infty$ and $l\gg 1$, we
write eq.(103) as

$\displaystyle
S_{M_{l}}(f)\simeq\frac{8a(\bar{\rho})}{a^{\prime}(\bar{\rho})}\int\limits_{0}^{1/2}\frac{\lambda(x)a^{\prime}(\bar{\rho})}{\lambda^{2}(x){a^{\prime}}^{2}(\bar{\rho})+4\pi^{2}f^{2}}\quantity(1+\frac{\lambda(x)}{4})\frac{1}{\lambda(x)}.$
(138)

Using the variable transform defined in eq.(75) and
$\lambda(x)=4\pi^{2}x^{2}$, the above equation can be written in the following
simplified form,

$\displaystyle\tilde{S}_{M_{l}}(f)$
$\displaystyle\simeq\frac{a}{2\pi^{2}\sqrt{2\pi
a^{\prime}(\bar{\rho})}}f^{-3/2}\times\int\limits_{0}^{\infty}\frac{y^{-3/4}}{1+y}\differential{y}=\frac{a}{2\sqrt{\pi^{3}a^{\prime}(\bar{\rho})}}f^{-3/2},$
(139)

where, we used eq.(116) and ignored the term of $\order{f^{1/2}}$. Thus we
derive eq.(103), appeared in the main-text.

##  Appendix C Evolution equations of correlation functions

###  C.1 Different time current-current correlation

The stochastic update rules for the different-time and different-space product
function of currents $\mathcal{Q}_{i}(t)\mathcal{Q}_{j}(t^{\prime})$, can be
written for $t>t^{\prime}$ as

$\displaystyle\mathcal{Q}_{i}(t+\differential{t})\mathcal{Q}_{j}(t^{\prime})=\begin{cases}\mathbf{events}&\mathbf{probabilities}\\\
\quantity(\mathcal{Q}_{i}(t)+1)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{2}\hat{a}_{i}(t)\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+2)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{4}\hat{a}_{i}(t)\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-1)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{2}\hat{a}_{i+1}(t)\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-2)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{4}\hat{a}_{i+1}(t)\differential{t}\\\
\mathcal{Q}_{i}(t)\mathcal{Q}_{j}(t^{\prime})&1-\Sigma\differential{t},\end{cases}$
(140)

where $\Sigma=3\quantity(\hat{a}_{i}(t)+\hat{a}_{i+1}(t))/4$. Using these
rules, we write the evolution equation of the two point current-current
correlation function as

$\displaystyle\partialderivative{t}C^{\mathcal{Q}\mathcal{Q}}_{r}(t,t^{\prime})=\left\langle\quantity[\hat{a}_{0}(t)-\hat{a}_{1}(t)]\mathcal{Q}_{r}(t^{\prime})\right\rangle=\quantity[C_{r}^{\hat{a}\mathcal{Q}}(t,t^{\prime})-C_{r-1}^{\hat{a}\mathcal{Q}}(t,t^{\prime})],$
(141)

which has appeared in eq.(19) in the main-text.

###  C.2 Different time mass-current correlation

The stochastic update rules for the different-time and different-space product
function of mass and current $m_{i}(t)\mathcal{Q}_{j}(t^{\prime})$, can be
written for $t>t^{\prime}$ as

$\displaystyle
m_{i}(t+\differential{t})\mathcal{Q}_{j}(t^{\prime})=\begin{cases}\quantity(m_{i}(t)+1)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{2}\hat{a}_{i+1}(t)dt\\\
\quantity(m_{i}(t)+1)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{2}\hat{a}_{i-1}(t)dt\\\
\quantity(m_{i}(t)+1)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{4}\hat{a}_{i+1}(t)dt\\\
\quantity(m_{i}(t)+2)\mathcal{Q}_{j}(t^{\prime})&\frac{1}{4}\hat{a}_{i-1}(t)dt\\\
\quantity(m_{i}(t)-2)\mathcal{Q}_{j}(t^{\prime})&\hat{a}_{i}(t)dt\\\
m_{i}(t)\mathcal{Q}_{j}(t^{\prime})&[1-\Sigma dt],\end{cases}$ (142)

where $\Sigma=3\quantity(\hat{a}_{i}(t)+\hat{a}_{i+1}(t))/4+\hat{a}_{i}(t)$.
Using these rules, we write the evolution equation of the two point mass-
current correlation function as

$\displaystyle\partialderivative{t}C^{m\mathcal{Q}}_{r}(t,t^{\prime})=\left\langle\quantity[\hat{a}_{L-1}(t)-2\hat{a}_{0}(t)+\hat{a}_{1}(t)]\mathcal{Q}_{r}(t^{\prime})\right\rangle\simeq
a^{\prime}(\bar{\rho})\sum\limits_{k}\Delta_{r,k}C^{m\mathcal{Q}}_{k}(t,t^{\prime}),$
(143)

which has appeared in the eq.(26).

###  C.3 Equal time current-current correlation

In the Manna sandpile, during each toppling two particles can hop to each
neighbouring site intependently and it may simulataneously create current at
two neighbouring bonds. In the following, we write the update equation of the
two point product function of integrated current,

$\mathcal{Q}_{i}(t+\differential{t})\mathcal{Q}_{j}(t+\differential{t})=\\\
\begin{cases}&\mathit{probabilities}\\\
\quantity(\mathcal{Q}_{i}(t)-1)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i+1,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-1)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{i}\delta_{i-1,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-1)\mathcal{Q}_{j}(t)&\frac{1}{2}\hat{a}_{i+1}(1-\delta_{i+1,j}-\delta_{i,j})\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+1)\mathcal{Q}_{j}(t)&\frac{1}{2}\hat{a}_{i}(1-\delta_{i-1,j}-\delta_{i,j})\differential{t}\\\
\mathcal{Q}_{i}(t)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{j+1}(1-\delta_{i-1,j}-\delta_{i,j})\differential{t}\\\
\mathcal{Q}_{i}(t)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{j}(1-\delta_{i+1,j}-\delta_{i,j})\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+2)\quantity(\mathcal{Q}_{j}(t)+2)&\frac{1}{4}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-2)\quantity(\mathcal{Q}_{j}(t)-2)&\frac{1}{4}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)+2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i}\quantity(1-\delta_{i,j})\differential{t}\\\
\mathcal{Q}_{i}(t)\quantity(\mathcal{Q}_{j}(t)+2)&\frac{1}{4}\hat{a}_{j}\quantity(1-\delta_{i,j})\differential{t}\\\
\quantity(\mathcal{Q}_{i}(t)-2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i+1}\quantity(1-\delta_{i,j})\differential{t}\\\
\mathcal{Q}_{i}(t)\quantity(\mathcal{Q}_{j}(t)-2)&\frac{1}{4}\hat{a}_{j+1}\quantity(1-\delta_{i,j})\differential{t}\\\
\mathcal{Q}_{i}(t)\mathcal{Q}_{j}(t)&\quantity(1-\Sigma\differential{t}),\end{cases}$
(144)

where $\Sigma\differential{t}$ is the probability of happening all the events,
mentioned in the update rules in the time-interval $t$ and
$t+\differential{t}$. Using this update rules and eq.(23), we can write the
evolution equation of the equal-time correlation function of current as

$\displaystyle\derivative{t}C_{r}^{\mathcal{Q}\mathcal{Q}}(t,t)=$
$\displaystyle\Gamma_{r}(t)+\left\langle\quantity[\hat{a}_{0}(t)-\hat{a}_{1}(t)]\mathcal{Q}_{r}(t)\right\rangle_{c}+\left\langle\mathcal{Q}_{0}\quantity[\hat{a}_{r}-\hat{a}_{r+1}]\right\rangle_{c}$
$\displaystyle\simeq$
$\displaystyle\Gamma_{r}(t)+a^{\prime}\quantity[c^{m\mathcal{Q}}_{r}(t,t)-c^{m\mathcal{Q}}_{r-1}(t,t)]+a^{\prime}\quantity[c^{m\mathcal{Q}}_{L-r}(t,t)-c^{m\mathcal{Q}}_{L-r-1}(t,t)].$
(145)

where,

$\displaystyle\Gamma_{r}\equiv\Gamma_{i,j}=\frac{3}{2}\delta_{i,j}\left\langle\hat{a}_{i+1}+\hat{a}_{i}\right\rangle-\frac{1}{2}\delta_{i+1,j}\left\langle\hat{a}_{i+1}\right\rangle-\frac{1}{2}\delta_{i-1,j}\left\langle\hat{a}_{i}\right\rangle.$
(146)

The solution of eq.(C.3) appeared in eq.(II.3) in the main-text.

###  C.4 Mass and integrated current correlation

The update rules of the temporal evolution equation of the product function of
mass and integrated current, are given as

$\displaystyle
m_{i}(t+\differential{t})\mathcal{Q}_{j}(t+\differential{t})=\begin{cases}\quantity(m_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i+1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{i}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{i-1}\delta_{i-2,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{i-1}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(\mathcal{Q}_{j}(t)+2)&\frac{1}{4}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+2)\quantity(\mathcal{Q}_{j}(t)+2)&\frac{1}{4}\hat{a}_{i-1}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)+2)\quantity(\mathcal{Q}_{j}(t)-2)&\frac{1}{4}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(\mathcal{Q}_{j}(t)-2)&\frac{1}{4}\hat{a}_{i}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\mathcal{Q}_{j}(t)&\frac{1}{2}\hat{a}_{i+1}\quantity(1-\delta_{i,j}-\delta_{i+1,j})\differential{t}\\\
\quantity(m_{i}(t)+1)\mathcal{Q}_{j}(t)&\frac{1}{2}\hat{a}_{i-1}\quantity(1-\delta_{i-1,j}-\delta_{i-2,j})\differential{t}\\\
\quantity(m_{i}(t)-2)\mathcal{Q}_{j}(t)&\frac{1}{2}\hat{a}_{i}\quantity(1-\delta_{i,j}-\delta_{i-1,j})\differential{t}\\\
\quantity(m_{i}(t)+2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i+1}\quantity(1-\delta_{i,j})\differential{t}\\\
\quantity(m_{i}(t)-2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i}\quantity(1-\delta_{i-1,j})\differential{t}\\\
\quantity(m_{i}(t)-2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i}\quantity(1-\delta_{i,j})\differential{t}\\\
\quantity(m_{i}(t)+2)\mathcal{Q}_{j}(t)&\frac{1}{4}\hat{a}_{i-1}\quantity(1-\delta_{i-1,j})\differential{t}\\\
m_{i}(t)\quantity(\mathcal{Q}_{j}(t)+1)&\frac{1}{2}\hat{a}_{j}\quantity(1-\delta_{i,j}-\delta_{i,j-1}-\delta_{i,j+1})\differential{t}\\\
m_{i}(t)\quantity(\mathcal{Q}_{j}(t)-1)&\frac{1}{2}\hat{a}_{j+1}\quantity(1-\delta_{i,j}-\delta_{i,j+1}-\delta_{i,j+2})\differential{t}\\\
m_{i}(t)\quantity(\mathcal{Q}_{j}(t)+2)&\frac{1}{4}\hat{a}_{j}\quantity(1-\delta_{i,j}-\delta_{i,j+1})\differential{t}\\\
m_{i}(t)\quantity(\mathcal{Q}_{j}(t)-2)&\frac{1}{4}\hat{a}_{j+1}\quantity(1-\delta_{i,j}-\delta_{i,j+1})\differential{t}.\end{cases}$
(147)

Using these rules, we write the evolution equation of the two point mass and
integrated current correlation as

$\displaystyle\derivative{t}C_{r}^{m\mathcal{Q}}(t,t)=f_{r}(t)+\sum\limits_{k}\Delta_{r,k}C_{k}^{\hat{a}\mathcal{Q}}(t,t),$
(148)

where $f_{r}(t)$, the source term of the equal time mass-integrated current
correlation, which, in the steady state can be written as

$\displaystyle
f_{r}(t)=\quantity[C^{m\hat{a}}_{r}(t,t)-C^{m\hat{a}}_{r+1}(t,t)]+\frac{7a(\bar{\rho})}{2}\delta_{0,r+1}-\frac{7a(\bar{\rho})}{2}\delta_{0,r}+\frac{a(\bar{\rho})}{2}\quantity(\delta_{0,r-1}-\delta_{0,r+2}).$
(149)

The Fourier transform of eq. (148) is used to get eq.(32) in the main-text.

###  C.5 Mass-mass correlation function

The equal time mass-mass correlation function is important to calculate the
mass-activity correlation function and the power spectrum of subsystem mass.
The update rules are the following,

$\displaystyle
m_{i}(t+\differential{t})m_{j}(t+\differential{t})=\begin{cases}\quantity(m_{i}(t)+1)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(m_{j}(t)-2)&\frac{1}{2}\hat{a}_{i+1}\delta_{i+1,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i+1}\delta_{i+2,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(m_{j}(t)-2)&\frac{1}{2}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i}\delta_{i+1,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i-1}\delta_{i-2,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(m_{j}(t)-2)&\frac{1}{2}\hat{a}_{i-1}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)+1)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{i-1}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+1)m_{j}(t)&\frac{1}{2}\hat{a}_{i+1}\quantity(1-\delta_{i,j}-\delta_{i+1,j}-\delta_{i+2,j})\differential{t}\\\
\quantity(m_{i}(t)-2)m_{j}(t)&\frac{1}{2}\hat{a}_{i}\quantity(1-\delta_{i-1,j}-\delta_{i,j}-\delta_{i+1,j})\differential{t}\\\
\quantity(m_{i}(t)+1)m_{j}(t)&\frac{1}{2}\hat{a}_{i-1}\quantity(1-\delta_{i-2,j}-\delta_{i-1,j}-\delta_{i,j})\differential{t}\\\
m_{i}(t)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{j+1}\quantity(1-\delta_{i,j}-\delta_{i,j+1}-\delta_{i,j+2})\differential{t}\\\
m_{i}(t)\quantity(m_{j}(t)-2)&\frac{1}{2}\hat{a}_{j}\quantity(1-\delta_{i,j-1}-\delta_{i,j}-\delta_{i,j+1})\differential{t}\\\
m_{i}(t)\quantity(m_{j}(t)+1)&\frac{1}{2}\hat{a}_{j-1}\quantity(1-\delta_{i,j-2}-\delta_{i,j-1}-\delta_{i,j})\differential{t}\\\
\quantity(m_{i}(t)+2)\quantity(m_{j}(t)+2)&\frac{1}{4}\hat{a}_{i+1}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+2)\quantity(m_{j}(t)-2)&\frac{1}{4}\hat{a}_{i+1}\delta_{i+1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(m_{j}(t)+2)&\frac{1}{4}\hat{a}_{i}\delta_{i-1,j}\differential{t}\\\
\quantity(m_{i}(t)-2)\quantity(m_{j}(t)-2)&\frac{1}{4}\hat{a}_{i}\delta_{i,j}\differential{t}\\\
\quantity(m_{i}(t)+2)m_{j}(t)&\frac{1}{4}\hat{a}_{i+1}\quantity(1-\delta_{i,j}-\delta_{i+1,j})\differential{t}\\\
\quantity(m_{i}(t)-2)m_{j}(t)&\frac{1}{4}\hat{a}_{i}\quantity(1-\delta_{i-1,j}-\delta_{i,j})\differential{t}\\\
m_{i}(t)\quantity(m_{j}(t)+2)&\frac{1}{4}\hat{a}_{j+1}\quantity(1-\delta_{i,j}-\delta_{i,j+1})\differential{t}\\\
m_{i}(t)\quantity(m_{j}(t)-2)&\frac{1}{4}\hat{a}_{j}\quantity(1-\delta_{i,j-1}-\delta_{i,j})\differential{t}\\\