\section{Introduction}
\label{sec:intro}
\noindent 

The field of active matter has developed over the last decades to provide a physical description of classical many-body systems operating far from thermodynamic equilibrium~\cite{Fodorreview,Jülicherreview,Marchettireview}. A prominent class of such active matter are living systems: Schools of fishes~\cite{Filella2018}, flocks of birds~\cite{Ballerini2008}, and colonies of bacteria~\cite{Damton2010,Thutupalli2015a} can all exhibit collective dynamics that are manifestly out of equilibrium. However, the non-equilibrium activity of biological assemblies at smaller subcelluar scales is not always straightforward to discern~\cite{mackintosh2010,Gnesottoreview}. Examples include the stochastic fluctuations of biological assemblies such as chromosomes~\cite{Weber}, the cytoskeleton~\cite{Lau2003,Guo2014,Fakhri2014,Jülicherreview,Brangwynne2008}, and cellular membranes~\cite{Betz2009,Turlier2016,Ben-Isaac2011}. Indeed, while these fluctuations can at first sight appear indistinguishable from thermal Brownian motion, they are in many cases driven by energy-consuming processes at molecular scales~\cite{Gnesottoreview,Fodorreview,Needleman2017b, Suzuki2015,mackintosh2010,Jülicherreview}.
This molecular scale activity can propagate to mesoscopic scales, giving rise to non-equilibrium dynamics that breaks detailed balance~\cite{Battle2016,Gladrow2016,Gladrow2017,Mura2018,
Gnesottoreview,Gnesotto2018,Li2018} or that violates the fluctuation dissipation theorem~\cite{Turlier2016,Betz2009,Mizuno2007,Martin2001,Fodor2015,Guo2014}. Soft driven assemblies can also be realized in synthetic systems, including chemical fueled synthetic fibers~\cite{Boekhoven2015} and crystals of active colloidal particles~\cite{Palacci2013}.
Numerous experimental studies showed how molecular non-equilibrium processes affect the mesoscopic mechanical properties of \textit{in vivo} biological assemblies~\cite{Lau2003,Guo2014,Fakhri2014,Jülicherreview},  \textit{in vitro} reconstituted cytoskeletal networks~\cite{Liu2006,Koenderink2009,Alvarado2017} and synthetic materials~\cite{Bertrand2012}.  
 It still remains unclear, however, how to characterize the non-equilibrium fluctuations of soft driven assemblies.

 To make further progress on characterizing active systems, various candidates for a reliable and informative non-equilibrium measure have been proposed. A natural and commonly used measure of the time-irreversibility of a process is the entropy production rate. In some cases, this measure is related to  the energy dissipation in a system~\cite{Seifertreview}. Recent studies made significant progress in inferring the entropy production rate from the observed trajectories~\cite{Li2018,Frishman2018,Seara2018}. In general, for complex systems it is unclear how to interpret measures of the partial entropy production rate or how to relate the measured quantities to the real entropy production rate of the full system. It is possible, however, to set a lower bound to the total entropy production from the observation of a few mesoscopic degrees of freedom~\cite{Roldan2018,Mura2018, Bisker2017,Esposito2012,Polettini2017}. An alternative approach of using area enclosing rates (AER) of stochastic trajectories in phase space as a metric for the breaking of detailed balance was presented in~\cite{Ghanta2017} and~\cite{Gonzalez2018}. A closely related concept---the cycling frequencies of the stochastic trajectory---was used to analyse the non-equilibrium behavior emerging in the mode-space trajectories of a probe filament in an active gel~\cite{Gladrow2016,Gladrow2017} or the dynamics of  a driven disordered elastic 
network near its isostatic point~\cite{Gnesotto2018}.

Despite the multitude of available non-equilibrium measures, it is still unclear how to use them to extract useful information about the nature of active driving in a system. The cycling frequencies can be used to investigate the  non-equilibrium dynamics emerging on different lengthscales in driven elastic networks~\cite{Mura2018}. In particular, the cycling frequencies measured from trajectories of two probe particles in an internally driven elastic network display a power-law behavior as function of the distance between the particles, with an exponent that depends on the dimensionality of the system; the prefactor of this scaling law depends on the statistical properties of the internal driving. Thus these experimentally accessible cycling frequencies and their associated scaling behavior provide a promising candidate for a non-equilibrium measure that may give provide access to properties of the internal driving.

In this work we present a detailed derivation of the scaling behavior of cycling frequencies for $d$-dimensional elastic networks with internal driving. Thus, we derive a theoretical framework that allows us to relate the cycling frequency to the lengthscale of the observation and to the properties of the network and of the active noise. Furthermore, we clarify the relation between the cycling frequencies and other non-equilibrium measures such as the area enclosing rate and the entropy production rate. 
\begin{figure}
\centering
\includegraphics[width=8cm]{intropanel5}
\caption{
A) 1-dimensional elastic chain of beads at temperature $T$ with spatially varying white-noise driving intensity $\alpha_i$, and the corresponding diffusion matrix $D$. The tracked pair of beads is indicated in green. B) 2-dimensional and C) 3-dimensional elastic networks investigated in Sec.~\ref{sec:Cubic lattice}.}
\label{fig:intropanel}
\end{figure}
\section{Driven elastic networks}
\label{sec1}
We use overdamped networks of elastically coupled beads suspended in a viscous fluid as a simple model for soft subcellular  assemblies~\cite{Yucht2013,Broedersz2014,Osmanovic2017,Mao2018}. For simplicity we choose units in which the elastic spring constant, the damping coefficient of the beads and the Boltzmann constant $k_B$ equal 1. The fluid is assumed to be at thermal equilibrium and the resulting thermal fluctuations in the system are thus modelled as Gaussian white noise processes acting independently on all the beads with the same amplitude $T$. Additional active force fluctuations drive the system out of thermal equilibrium, and are implemented as independent Gaussian white noise processes with spatially varying amplitudes $\alpha_i$. 

By modelling the active forces as ``white", we essentially restrict our model to systems in which the correlation times of the active driving are shorter than the intrinsic relaxation times of the network. This model is mathematically equivalent to embedding the beads in local thermal baths at temperatures $T+\alpha_i$~\cite{Falasco2015}.

This simplified description allows us to study the dynamics of the system using a Fokker-Planck equation:
\begin{equation}
\partial_t p(\mathbf{x},t)=-\nabla \cdot (\mathbf{A}\mathbf{x}-\mathbf{D}\nabla)p(\mathbf{x},t)\coloneqq-\nabla \cdot \mathbf{j}(\mathbf{x},t), 
\label{eq:fokker}
\end{equation}
Here $\mathbf{x}$ represents the displacements of the beads relative to their equilibrium positions, $p(\mathbf{x},t)$ is the probability distribution of $\mathbf{x}$ at time $t$. We also assume that the forces are linear in $\mathbf{x}$, i.e. $\mathbf{f}(\mathbf{x})=\mathbf{A}\mathbf{x}$ with a symmetric matrix $\mathbf{A}$, and $\mathbf{D}=\mathop{\rm diag}\{T+\alpha_1,\ldots, T+\alpha_N\}$ is the diffusion matrix. The right hand side of the Eq.~\eqref{eq:fokker} can be interpreted as the divergence of the probability current density $\mathbf{j}(\mathbf{x},t)=(\mathbf{A}\mathbf{x}-\mathbf{D}\nabla)p(\mathbf{x},t)$. At  steady-state, the non-vanishing dissipative probability currents constitute a measure for non-equilibrium in a system and thus play the key role in our approach.

\section{Cycling frequencies \& phase space torque}
\label{sec2}
The steady-state probability currents, $\mathbf{j}(\mathbf{x})$, are mathematical objects that capture the presence of non-equilibrium activity by revealing time-irreversibility of the dynamics at the level of the Fokker-Planck equation. This time-irreversibility manifests through the emergence of a mean velocity field $\mathbf{v}(\mathbf{x})$ in the coordinate space, which is related to the probability current through $\mathbf{j}(\mathbf{x})=\mathbf{v}(\mathbf{x})p(\mathbf{x})$~\cite{Risken}. Therefore, from an experimental perspective, an ideal way to quantify the non-equilibrium dynamics of a system would be to measure such a velocity field $\mathbf{v}(\mathbf{x})$. However, inferring the full $\mathbf{v}(\mathbf{x})$ field is a challenge on its own. The most straightforward approaches require a discretisation of the phase space. Such a measurement would require tracking many degrees of freedom for long time periods, which is difficult in practice.

Instead of inferring $\mathbf{v}(\mathbf{x})$ in full detail, one can alternatively measure some coarse-grained quantities related to $\mathbf{v}(\mathbf{x})$, which still retain key information about the non-equilibrium dynamics of the system. For instance, we could track a pair of degrees of freedom $\mathbf{x}_{\rm r}=\{x_i,x_j\}$ and measure the average angular velocity $\langle \dot{\beta}_{ij}\rangle$, or equivalently, the rate at which the trajectory revolves around the origin in this reduced 2-dimensional subspace (Fig.~\ref{fig:torque}). This simple measurement does not require any discretisation of phase space or inference of the force field. We shall refer to $\langle \dot{\beta}_{ij}\rangle$ as the \textit{cycling frequency}. 

In general, $\langle \dot{\beta}_{ij}\rangle$ may contain only limited information about $\mathbf{v}(\mathbf{x})$. For linear systems, however, the mean phase space velocity can be written as~\cite{Weiss2003}
\begin{equation}
\mathbf{v}(\mathbf{x})=\mathbf{\Omega}\mathbf{x},\quad \mathbf{\Omega}=\mathbf{A}+\mathbf{D}\mathbf{C}^{-1}
\label{eq:omega}
\end{equation}
where $\mathbf{C}=\langle \mathbf{x}\otimes\mathbf{x}\rangle$ is the steady-state covariance matrix obeying the Lyapunov equation
\begin{equation}
\mathbf{A}\mathbf{C}+\mathbf{C}\mathbf{A}^T=-2\mathbf{D}.
\label{eq:Lyapunov}
\end{equation}
Equation~\eqref{eq:omega} sets strong constraints on the structure of $\mathbf{v}(\mathbf{x})$: for dynamics projected on any 2-dimensional subspace $\{x_i,x_j\}$ the probability currents have an elliptical structure. The remaining information about the amplitude of the currents is set by the cycling frequencies $\langle \dot{\beta}_{ij} \rangle$.

To show this, we denote the velocity of the system in the reduced $ij$-subspace by $\mathbf{v}_{\rm r}(\mathbf{x})=\{ v_i(\mathbf{x}),v_j(\mathbf{x})\}$. Note, what we observe while looking at the $\mathbf{x}_r$-subspace only is a conditional mean $\langle \mathbf{v}_{\rm r}(\mathbf{x})|\mathbf{x}_r \rangle$. Similarly to Eq.~\eqref{eq:omega}, we find that $\langle\mathbf{v}_{\rm r}(\mathbf{x})|\mathbf{x}_{\rm r}\rangle=\mathbf{\Omega}_{\rm r}\mathbf{x}_{\rm r}$, with $\mathbf{\Omega}_{\rm r}=\mathbf{A}_\textup{eff}+\mathbf{D}_{[{\rm r},{\rm r}]}\mathbf{C}_{[{\rm r},{\rm r}]}^{-1}$~\cite{Mura2018}. Here $\mathbf{A}_\textup{eff}$ is a matrix such that $\langle \mathbf{f}_{\rm r}(\mathbf{x})|\mathbf{x}_{\rm r}\rangle=\langle \{f_i(\mathbf{x}),f_j(\mathbf{x})\}|\mathbf{x}_{\rm r}\rangle=\mathbf{A}_\textup{eff}\mathbf{x}_{\rm r}$; note, $\mathbf{C}_{[{\rm r},{\rm r}]}$ and $\mathbf{D}_{[{\rm r},{\rm r}]}$ are matrices of size $[2 \times 2]$, given by $\mathbf{C}_{[{\rm r},{\rm r}]}=\{\{c_{ii},c_{ij}\},\{c_{ji},c_{jj}\}\}$ and $\mathbf{D}_{[{\rm r},{\rm r}]}=\{\{d_{ii},d_{ij}\},\{d_{ji},d_{jj}\}\}.$

Next, we show that the eigenvalues of $\mathbf{\Omega}_{\rm r}$  coincide with the cycling frequencies $\langle \dot{\beta}_{ij} \rangle$. 
First, note that $\langle \dot{\beta}_{ij} \rangle$ is invariant under orientation preserving linear transformations of the reduced subspace. We can therefore work in \textit{covariance identity coordinates}, $\mathbf{y}_{\rm r}$, such that $\mathbf{C}_{[{\rm r},{\rm r}]}=\mathbf{I}$. In this basis, $\mathbf{\Omega}_{\rm r}$ takes a particularly simple form~\cite{Weiss2003}
\begin{equation}
\mathbf{\Omega}_{\rm r}=\left(\begin{array}{cc}
0 & \omega_{ij} \\
-\omega_{ij} & 0 \\
\end{array}\right),
\label{eq:omegaCIC}
\end{equation}
with the imaginary parts of its eigenvalues on the antidiagonal. This form of $\mathbf{\Omega}_{\rm r}$ implies that for $\mathbf{C}_{[{\rm r},{\rm r}]}=\mathbf{I}$, the probability current field has a circular structure. Using Eq.~\eqref{eq:omegaCIC}, we find that
\begin{align}
\label{eq:indepomeg}
\langle\dot{\beta}_{ij}|\mathbf{y}_r\rangle&=\mean{\left.\frac{\mathbf{y}_{\rm r}\times\dot{\mathbf{y}}_{\rm r}}{|\mathbf{y}_{\rm r}|^2}\right\vert\mathbf{y}_{\rm r}}=\frac{y_2\mean{\dot{y}_1|\mathbf{y}_r}-y_1\mean{\dot{y}_2|\mathbf{y}_r}}{y_1^2+y_2^2}\nonumber\\
&=\frac{y_2\mean{v_1(\mathbf{y})|\mathbf{y}_r}-y_1\mean{v_2(\mathbf{y})|\mathbf{y}_r}}{y_1^2+y_2^2}=\omega_{ij}.
\end{align}
This means that the conditional average of the angular velocity $\langle\dot{\beta}_{ij}|\mathbf{y}_r\rangle$ is $\mathbf{y}_r$-independent and equals $\omega_{ij}$ at all points in the reduced phase space. Hence averaging over $\mathbf{y}_{\rm r}$ leads to $\langle\dot{\beta}_{ij}\rangle=\omega_{ij}$, with $\omega_{ij}^2=\det\mathbf{\Omega}_{\rm r}$. 

Recently, new approaches have been developed to infer current fields in nonlinear systems by considering an expansion of the inferred force field~\cite{Frishman2018}. Up to first order these methods reduce to calculating area enclosing rates, which are indeed closely related to the cycling frequencies, as we discuss further below.

\begin{figure}
\centering
\includegraphics[width=0.8\linewidth]{torque4}
\caption{ Schematic trajectory in the coordinate space $[x_i,x_j]$ of two tracer beads. The light blue area enclosed in the triangle represents $\mathbf{x}_{\rm r}\times\mathbf{f}_{\rm r}(\mathbf{x})/2$ appearing in Eq.~\eqref{eq:torque}, which gives the area enclosing rate upon averaging over phase space.}
\label{fig:torque}
\end{figure}

One of our central objectives is to derive a relation between the observed currents and the properties of the system and the active driving. 
Given that the cycling frequencies are set by $\omega_{ij}$---the imaginary parts of the eigenvalues of $\mathbf{\Omega}_{\rm r}$---we can make further progress by showing that for a general linear system~\cite{Mura2018}
\begin{equation}
\omega_{ij}=\frac{1}{2}\frac{\mean{\mathbf{x}_{\rm r}\times\mathbf{f}_{\rm r}(\mathbf{x})}}{\sqrt{\det\mathbf{C}_r}},
\label{eq:torque}
\end{equation}
where $\mean{\mathbf{x}_{\rm r}\times\mathbf{f}_{\rm r}(\mathbf{x})}=\langle x_if_j(\mathbf{x})-x_jf_i(\mathbf{x})\rangle=(\mathbf{C}\mathbf{A}^T-\mathbf{A}\mathbf{C})_{ij}$ is the mean phase space torque (Fig.~\ref{fig:torque}). Intuitively, Eq.~\eqref{eq:torque} implies that for an overdamped linear system the mean phase space angular velocity is proportional to the mean phase space torque. A detailed derivation is presented in ~\cite{Mura2018}. Moreover, in covariance identity coordinates Eq.~\eqref{eq:torque} reduces to $\omega_{ij}=\frac{1}{2}(A^T-A)_{ij}=\Omega_{ij}$, in accord with previous studies~\cite{Weiss2003, Gladrow2016,Gladrow2017}. One can equivalently identify $(\mathbf{C}\mathbf{A}^T-\mathbf{A}\mathbf{C})_{ij}$ with mean area enclosing rates in the $ij$-subspace, as considered in~\cite{Ghanta2017}. Here, we focus on the cycling frequencies, since they are more directly related to the probability currents, which constitute the basis of our work. In some instances, however, the area enclosing rates turn out to be particularly advantageous to work with. In these cases we  shall briefly discuss how switching to the area enclosing rates simplifies the analysis (See Sec.~\ref{sec3}).

Since the cycling frequencies contain information about the amplitudes of the phase space probability currents, they can be related to the entropy production rate. For linear systems, in covariance identity coordinates, the full entropy production rate can be expressed as a weighted sum of the cycling frequencies squared: $\Pi=\sum_{i-\mathrm{odd}}^n \omega_{i,i+1}^2[(\mathbf{D}^{-1})_{i,i}+(\mathbf{D}^{-1})_{i+1,i+1}]$. However, in many experimental contexts it is typically impossible to track all the degrees of freedom or to resolve all steps of a process. This practical limitation motivated the introduction of various measures of reduced/apparent entropy production rate~\cite{Roldan2018,Mura2018, Bisker2017,Esposito2012,Polettini2017}. For the case of a two-point measurement, as the ones discussed in this paper, one can consider a reduced entropy production rate directly related to the cycling frequency: $\Pi_{\rm r}^{(2)}~=~\omega_{ij}^2\mathrm{Tr} (\mathbf{C}_{[{\rm r},{\rm r}]}\mathbf{D}_{[{\rm r},{\rm r}]}^{-1})$, which gives a lower bound on the full contribution to the total entropy production rate from the observed pair of degrees of freedom~\cite{Mura2018}.

\section{One-dimensional chain and diffusion equation}
\label{sec3}
While Eq.~\eqref{eq:torque} sets a relation between the cycling frequencies observable in an experiment and the properties of both the network and the active noise distribution, it is not straightforward to explicitly derive these  properties from the cycling frequencies. In the following sections, we use Eq.~\eqref{eq:torque} to build a framework for extracting specific information about the system from the behavior of the cycling frequencies $\omega_{ij}$.

We first consider the simplest case of a 1-dimensional chain of $2N-1$ beads coupled by harmonic springs. This example will help us build intuition for more complex lattices. To obtain insight into how non-equilibrium behavior manifests at different length scales, we consider a two-point non-equilibrium measure. Specifically, we study how the cycling frequency in the subspace of displacements of two chosen beads $\{ x_i,x_{i+r} \}$ depends on the distance $r$ between these beads. We do this for two scenarios:
\begin{enumerate}[i)]
\item Active noise present only at single site
\item Random spatial distribution of activities $\{\alpha_i\}$
\end{enumerate}
In case i) we plot  the cycling frequency, $\omega_{\rm single}(r)$, between the active bead and another bead at distance $r$. In case ii) we consider an ensemble of active noise distributions $\{\alpha_i\}$, in which the amplitude of the active noise at each site is drawn randomly from a probability distribution $p_\alpha$, with mean $\bar{\alpha}$ and variance $\sigma_\alpha^2$. The amplitudes $\alpha_i$ are spatially uncorrelated. We then calculate $\langle \omega^2(r) \rangle$ - the squared cycling frequency between two beads at distance $r$ averaged over the activity distributions $\alpha_i$. 
\begin{figure}
\centering
\includegraphics[width=8cm]{omega1D_4new}
\caption{A) Scaling behavior of the cycling frequencies as a fuction of distance between the beads, obtained for a 1D chain and different types of activity distributions, as indicated in the legend. B) Comparison between the cycling frequencies and the Area Enclosing Rates (AER) above the weak noise limit ($\bar{\alpha}/T=4$). All data points correspond to results obtained by numerically solving Lyapunov equation.}
\label{fig:omega1D}
\end{figure}
The first observation is that in both scenarios the cycling frequencies follow a power-law as a function of distance, as shown in Fig.~\ref{fig:omega1D}. The exponent in the random distribution scenario is independent of the probability distribution $p_\alpha$ of the intensities, but different from the single activity case. 

To understand the origin of the power-law behavior and to calculate the exponents, we use Eq.~\eqref{eq:torque} to derive analytical expressions for $\omega_{\rm single}(r)$ and $\langle\omega(r)^2\rangle$.
In the case of a 1-dimensional chain with spatially uncorrelated noise, the expression for the cycling frequency  (Eq.~\eqref{eq:torque}) reduces to:
\begin{equation}
\omega_{ij}=\frac{\widetilde{\partial}_{2}^2 c_{ij}}{\sqrt{\det\mathbf{C}_{[{\rm r},{\rm r}]}}},
\label{eq:omchain}
\end{equation}
where $c_{ij}$ indicates the elements of the covariance matrix, and $\widetilde{\partial}_{2}^2 c_{ij}$ denotes the discrete second derivative across rows: $\widetilde{\partial}_{2}^2 c_{ij}=c_{i,j+1}-2c_{i,j}+c_{i,j-1}$. Thus, this result reduces the problem of calculating $\omega_{ij}$ to finding the covariance matrix $\mathbf{C}$.

Motivated by the structure of $\mathbf{D}$, we decompose $\mathbf{C}=T\overline{\mathbf{C}}+\bar{\alpha}\mathbf{C}^{*}$ into equilibrium $(\overline{\mathbf{C}})$ and non-equilibrium $(\mathbf{C}^{*})$ parts (Fig.\ref{fig:numerC}A). For the 1-dimensional chain, the Lyapunov equation (see Eq.~\eqref{eq:Lyapunov}) is equivalent to:
\begin{align}
\widetilde{\partial}_{1}^2 \overline{c}_{ij}+\widetilde{\partial}_{2}^2 \overline{c}_{ij}&=-2\delta_{ij} \label{eq:diffec} \\
\widetilde{\partial}_{1}^2 c^{*}_{ij}+\widetilde{\partial}_{2}^2 c^{*}_{ij}&=-2\delta_{ij}\frac{\alpha_i}{\bar{\alpha}} \label{eq:diffnc},
\end{align}

At equilibrium detailed balance is preserved, which implies $\widetilde{\partial}_{2}^2 \overline{c}_{ij}=0\ \forall_{i\neq j}$ (see Eq.~\eqref{eq:omchain}). We can therefore replace $\widetilde{\partial}_{2}^2 c_{ij}$ with $\widetilde{\partial}_{2}^2 c^{*}_{ij}$ in Eq.~\eqref{eq:omchain}. Then, expanding the expression for the cycling frequency in powers of $(\bar{\alpha}/T)$, we get:
\begin{equation}
\omega_{ij}=\frac{\bar{\alpha}}{T}\frac{\pdyc^{*}_{ij}}{\sqrt{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}}}+\mathcal{O}\left(\frac{\bar{\alpha}^2}{T^2}\right)
\label{eq:omegalin}
\end{equation}
Up to linear order in $(\bar{\alpha}/T)$ the contributions from $\overline{\mathbf{C}}$ and $\mathbf{C}^{*}$ separate; consequently $\omega_{ij}$ becomes linear in $\{\alpha_i\}$. This linearity appearing in the limit of weak activity will later allow us to calculate cycling frequencies averaged over different realizations of the activity $\{\alpha_i\}$. Note, if instead of $\omega_{ij}$, we consider the area enclosing rates ($\mathrm{AER}=\bar{\alpha}\pdyc^{*}_{ij}$), the factor $\sqrt{\det \mathbf{C}_{[{\rm r},{\rm r}]}}$ does not enter, implying that the expression for the area enclosing rates is linear in $\{\alpha_i\}$ irrespective of the magnitude of $\bar{\alpha}$.


\subsection{Activity at a single site}
To obtain insight into what determines the cycling frequencies in a concrete example, we first find the solution to Eq.~\eqref{eq:diffnc} for the case of activity appearing only at a single site in the center of the chain. Later, we will use this solution to construct $\mathbf{C}^{*}$ for a more general case. For now, let us assume that $\alpha_i=\alpha\delta_{iN}$ and $\bar{\alpha}=\alpha$. 

We can think of Eq.~\eqref{eq:diffnc} as a discretized stationary diffusion equation with a single source with a divergence of 2 in the middle of the $\mathbf{C}^{*}$ matrix, and with absorbing boundary conditions at the edges. The absorbing boundaries in the diffusion equation reflect the fixed boundary conditions for the elastic chain. We denote by $r=\sqrt{(i-N)^2+(j-N)^2}$ the distance from the center of $\mathbf{C}^{*}$. If we consider a continuous analogue of our discrete problem and neglect the boundary conditions, we can assume a rotational symmetry of the solution $c^{*}(r)$. The corresponding continuous diffusion equation takes the form: $\frac{1}{r}\partial_r r\partial_r c^{*}(r)=0$. Consequently, $r\partial_rc^{*}(r)=-a$ and $c^{*}(r)=-a\ln(r)+b$. 

One could also argue for $\partial_rc^{*}(r)\sim 1/r$ scaling of the ``covariance current", by demanding that the total ``covariance flux" through a circle of radius $r$, centred at $0$ is independent of $r$ and equals 2---the divergence of the source. This also allows us to identify $a=1/\pi$. The integration constant $b$ is system-size dependent and has to be set such that the covariance vanishes at the edge of the covariance matrix.


\begin{figure}
\centering
\includegraphics[width=8cm]{newC2}
\caption{A) Profiles of the matrices $\overline{\mathbf{C}}$ and $\mathbf{C}^{*}$ in the single activity case. For visual purposes, the discrete data points have been interpolated to a 2D-surface. $\overline{c}_{ij}$ is linear in both indices, resulting in $\omega_{ij}=0$ for $\alpha_i=0$. B) Values of $c^{*}_{ij}$ versus the distance from the center of the matrix $r=\sqrt{(i-N)^2+(j-N)^2}$. Orange points correspond to the entries at positions (i,N). The inset depicts the same plot in log-linear scale. }
\label{fig:numerC}
\end{figure}

The functional form that we obtain from this approximate analysis accurately describes the actual numerically obtained values of $c^{*}_{ij}$ far from the boundaries (see Fig.~\ref{fig:numerC}B). Note, the deviations that appear close to the boundaries are due to neglecting the absorbing boundary conditions and not due to the discrete nature of the problem. We can also consider a continuous limit of the problem, in which the chain is replaced by a string, by taking the limit $N\to\infty$, $k\to \infty$, while keeping $k/N=const$. In this limit the discrete diffusion equation is replaced by a continuous one, but the boundary effects still play the same role. We shall return to the continuous limit in Sec.~\ref{subsec4}, where we  discuss more complex networks.

We can use the approximate form of $c^{*}(r)$ together with  Eq.\eqref{eq:omegalin} to calculate $\omega_{\rm single}(r)$. For large $N$, we can replace $\pdyc^{*}_{N,N+r}$ with $\partial_r^2 c^{*}(r)$, to arrive at:
\begin{equation}
\omega_{\rm single}(r)=\frac{\alpha}{T}\frac{1}{\pi r^2}\frac{1}{\sqrt{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)}}+\mathcal{O}\left(\frac{\alpha^2}{T^2}\right)
\label{eq:omegasin}
\end{equation}
Interestingly, it turns out that $b$, which is in general unknown, does not enter the equation for $\omega(r)$. 

To find the equilibrium part $\overline{\mathbf{C}}$, we make the following observation. At equilibrium all cycling frequencies $\omega_{ij}$ must vanish, which combined with Eq.~\eqref{eq:omchain} gives: $\widetilde{\partial}_{2}^2\overline{c}_{ij}=0\ \forall_{i\neq j}$. Using Eq.~\eqref{eq:diffec} at point $(i,i)$ and the symmetry of $\overline{\mathbf{C}}$ we find $\widetilde{\partial}_{2}^2\overline{c}_{ii}=-1$. In general, we can therefore write the equation for $\overline{\mathbf{C}}$ as $\widetilde{\partial}_{2}^2\overline{c}_{ij}=-\delta_{ij}$. Note that this condition is equivalent to a discrete stationary diffusion equation in 1 dimension, with a single source at site $i$ and with absorbing boundary conditions. This implies that $\overline{c}_{ij}$ is linear in both indices and one can easily verify that $\overline{c}_{ij}=\min(i,j)-ij/(2N)$ satisfies Eq.~\eqref{eq:diffec}. Using this solution, we find that $(\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r))^{-\frac{1}{2}}=(\frac{1}{2}r(N-r))^{-\frac{1}{2}}\sim r^{-\frac{1}{2}}$ for $r\ll N$. Therefore, $\omega_{\rm single}(r)\sim r^{-\frac{5}{2}}$, as shown together with the numerical results in Fig.~\ref{fig:1hot}.
\begin{figure}
\centering
\includegraphics[width=8cm]{Fig_single_source3}
\caption{ A) Scaling behavior of $\omega_{\rm single}(r)$ as a function of the distance between the beads and comparison with the analytical prediction (dashed line) in Eq.~\eqref{eq:omegasin}. Below, the behaviour of the B) equilibrium and C) non-equilibrium contributions to $\omega_{\rm single}(r)$. The contribution presented in C) coincides with the area enclosing rates. All data points correspond to results obtained by numerically solving Lyapunov equation.
}
\label{fig:1hot}
\end{figure}
The good agreement between the numerical and analytical results allows us to conclude that the scaling exponent is determined by the $\ln(r)$-like profile of $\mathbf{C}^{*}$, which in turn is set by the dimensionality of the system. One could in principle find an analytical solution for $c^{*}(r)$ that accounts for the boundary conditions, but the $\ln(r)$ scaling captures the essential features.

\subsection{Spatially varying activity}
Equipped with the results from the previous section, we now consider a system with spatially varying activity. First, we will further clarify the connection between calculating cycling frequencies and solving a discretized steady-state diffusion equation for the covariance function. To this end, we consider a generic activity distribution $\{ \alpha_i \}$ and plot the corresponding active part of the covariance matrix $\mathbf{C}^{*}$ obtained by solving Eq.~\eqref{eq:diffnc} (see Fig.~\ref{fig:nowypiekny}). From the form of Eq.~\eqref{eq:omchain}, we see that the cycling frequencies $\{\omega_{ij}\}_{j=1,\ldots,2N-1}$ are proportional to the curvature of the line $\{(i,j,c^{*}_{ij})\}_{j=1,\ldots,2N-1}$, as illustrated by the plots in Fig.~\ref{fig:nowypiekny}.
\begin{figure}
\centering
\includegraphics[width=8cm]{nowypiekny2}
\caption{Profile of $\mathbf{C}^{*}$ in the case of spatially correlated amplitudes of the active noise, together with the lines $\{(i,j,c^{*}_{ij})\}_{j=i,\ldots,2N-1}$ for $i=10, 17, 23$. The color of a line indicates the sign of its curvature. The red (blue) points in the plots of $|\omega_{i,j}|$ correspond to positive (negative) cycling frequencies. Here, the profile of $\mathbf{C}^{*}$ is presented for activities $\{\alpha_i-\bar{\alpha}\}$. This transformation of activities is justified in App.~\ref{appendix1}.}
\label{fig:nowypiekny}
\end{figure}
The connection with the steady-state diffusion equation (see Eq.~\eqref{eq:diffnc}) allows us to understand how a given distribution of activities translates to a particular profile of the $\mathbf{C}^{*}$ matrix and how this in turn determines the behavior of the cycling frequencies. 

In general the amplitudes of the active noise may be spatially correlated. Here, however, we restrict ourselves exclusively to the case of spatially uncorrelated activities, which is valid in the limit of distances larger than the correlation length of the amplitudes. To the ``$i^{th}$" bead we assign a randomly sampled amplitude $\alpha_i$. We assume all $\alpha_i$ to be pairwise independent and identically distributed with distribution $p(\alpha)$. For simplicity we index the beads so that the bead in the center of the system has index $0$. To calculate $\omega(2r)=\omega_{-r,r}$, we need to determine $\widetilde{\partial}_{2}^2 c^{*}_{-r,r}$. For a given activity distribution $\{ \alpha_i\}$ we can exploit the linearity of Eq.~\eqref{eq:diffnc} to obtain the corresponding $\mathbf{C}^{*}(\{ \alpha_i\} )$ as a superposition of single-source solutions. Thus, we can write:
\begin{equation}
\widetilde{\partial}_{2}^2 c^{*}_{-r,r}(\{\alpha_i\} )=\sum_z \widetilde{\partial}_{2}^2 c^{*}_{-r,r}(\{ \alpha_i\delta_{iz} \} )
\label{eq:superposition}
\end{equation}
For beads far enough from the boundary, we approximate $\mathbf{C}^{*} (\{ \alpha_i\delta_{ij}\} )$ by a logarithmic decay centred at $(j,j)$  to obtain:
\begin{align}
\widetilde{\partial}_{2}^2 c^{*}_{-r,r}(\{ \alpha_i\delta_{iz} \} )&=\frac{1}{\pi}\frac{\alpha_z}{\bar{\alpha}} \frac{(r+z)^2-(r-z)^2}{((r+z)^2+(r-z)^2)^2}\nonumber \\
&=\frac{1}{\pi}\frac{\alpha_z}{\bar{\alpha}}\frac{rz}{(r^2+z^2)^2}
\label{eq:partfromz}
\end{align}
Combining Eqs.~\eqref{eq:superposition} and \eqref{eq:partfromz}, we calculate $\langle(\widetilde{\partial}_{2}^2 c^{*}_{-r,r}(\{\alpha_i\} ))^2\rangle$, which is the main factor in the expression for $\langle\omega^2(2r)\rangle$.
\begin{widetext}
\begin{align}
(\pi\bar{\alpha})^2\langle(\widetilde{\partial}_{2}^2 c^{*}_{-r,r}(\{\alpha_i\} ))^2\rangle&=\mean{\left(\sum_{z}\frac{rz\alpha_z}{(r^2+z^2)^2}\right)\left(\sum_{z'}\frac{rz'\alpha_{z'}}{(r^2+z'^2)^2}\right)}=
\sum_{z}\frac{r^2z^2\langle\alpha_z^2\rangle}{(r^2+z^2)^4}+\sum_{z}\frac{rz\mean{\alpha_z}}{(r^2+z^2)^2}\sum_{z'\neq z}\frac{rz'\mean{\alpha_{z'}} }{(r^2+z'^2)^2}\nonumber\\
&=\sum_{z}\frac{r^2z^2(\langle\alpha_z^2\rangle-\mean{\alpha_z}^2)}{(r^2+z^2)^4} \approx \sigma_\alpha^2\int_{-\infty}^{\infty} \frac{r^2z^2dz}{(r^2+z^2)^4}=\frac{\pi\sigma_\alpha^2}{16r^3}
\label{eq:calculate1D}
\end{align}
\end{widetext}
In the second line of this result, we used that $\sum_{z'\neq z}\frac{rz'}{(r^2+z'^2)^2}=-\frac{rz}{(r^2+z^2)^2}$ and approximated the sum by an integral. Evaluating the integral and rescaling $2r\to r$, we arrive at the final result:
\begin{equation}
\langle\omega^2(r)\rangle_\alpha=\frac{\sigma_\alpha^2}{T^2}\frac{1}{2\pi r^3}\frac{1}{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)}.
\label{eq:scaling}
\end{equation}
Given the asymptotic behavior $1/\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)\sim r^{-1}$ for $r\ll N$ we conclude that in the limit of weak activity $\sqrt{\langle\omega^2(r)\rangle}\sim r^{-2}$ (See Fig.~\ref{fig:omega1D}). Importantly, apart from reproducing the observed exponent of the power-law, our result gives a correct prediction for the prefactor, which contains information about the variance of the active forces.

\section{$d$-dimensional lattices}
\label{sec4}
\begin{figure}
\centering
\includegraphics[width=8cm]{d_dimension3}
\caption{Scaling behavior of the cycling frequencies as function of the distance between the beads. The results are obtained for different lattices and a folded Gaussian activity distribution~\cite{Mura2018}. Triangular and square markers represent triangular and square or cubic lattices, respectively. Light (dark) blue represent triangular networks with zero (finite) rest length springs. In all cases  we used $\frac{\bar{\alpha}}{T}=0.15$. For computational convenience we determined the ensemble average by performing a spatial average. All data points correspond to results obtained by numerically solving Lyapunov equation.
}
\label{fig:omega23D}
\end{figure}
\subsection{Cubic lattice}
\label{sec:Cubic lattice}
To explain the origin of the scaling behavior of $\omega(r)$ for multidimensional networks we now focus on the simplest possible case of a $d$-dimensional cubic lattice. Importantly, the calculation presented for this case also provides us with intuition for more complex lattices. 
Let us denote the bead indices corresponding to $d$ independent directions with $n_1,\ldots ,n_d$. We will denote the elements of the covariance matrix $\mathbf{C}$ as $c_{n_1,\ldots,n_d;\bar{n}_1,\ldots,\bar{n}_d}\coloneqq c_{\mathbf{n},\bar{\mathbf{n}}}$. We assume zero-restlength for the springs, so that the degrees of freedom corresponding to different directions decouple. Therefore, by $\mathbf{C}$ we actually mean the covariance matrix of only these degrees of freedom that correspond to a single chosen direction, for instance the one corresponding to the index $n_1$.

For this particular network, the Lyapunov equation is equivalent to:
\begin{equation}
\label{eq:diff}
\left(\sum_{i=1}^{d}\tilde{\partial}_{n_i}^2+\sum_{i=1}^{d}\tilde{\partial}_{\bar{n}_i}^2\right)c_{\mathbf{n},\bar{\mathbf{n}}}=-2d_{\mathbf{n},\bar{\mathbf{n}}}
\end{equation}
Similarly to the 1-dimensional case, here we recognise a discretized stationary diffusion equation in $2d$ dimensions, with the divergence of the sources given by the elements of $\mathbf{D}$.
For convenience we index the beads such that the one in the middle of the lattice is $(0,\ldots,0)$. 

Our goal is to calculate the cycling frequency $\omega_{r,\ldots,r;-r,\ldots,-r}=\omega(2\sqrt{d}r)$. Here we consider a particular case of relative position of the beads with respect to the principal directions of the lattice. It turns out, however, that the final result depends only on the distance between the beads. From Eq.~\eqref{eq:torque}, we find that:
\begin{equation}
\label{eq:omegalind}
\omega_{r,\ldots,r;-r,\ldots,-r}\approx \frac{(\bar{\alpha}/T)}{\sqrt{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}}}\sum_{i=1}^{d}\left. \tilde{\partial}_{n_i}^2 c^{*}_{\mathbf{n},\bar{\mathbf{n}}}\right|_{r,\ldots,r;-r,\ldots,-r}
\end{equation}
Following the procedure used for the 1-dimensional chain, we begin with finding the solution to a single source problem with one active bead at site $(0,\ldots,0)$. As before, we will then use this solution as a Green's function for our diffusion problem with a generic activity distribution.

Taking a continuous limit of the diffusion equation and neglecting the boundary conditions, we expect that $\partial_r c(r)\sim 1/r^{2d-1}$ and consequently $c(r)\sim 1/r^{2d-2}$, where $r$ is the distance from the center of the $2d$-dimensional covariance matrix. Therefore, for a single active bead at site $(0,\ldots,0)$, we obtain:
\begin{equation}
c^{*}_{\mathbf{n},\bar{\mathbf{n}}}=a_d\left( \sum_{i=1}^{d} n_i^2 +\sum_{i=1}^{d} \bar{n}_i^2 \right)^{-(d-1)}\mkern-18mu =a_d\left( \mathbf{n}^2 +\bar{\mathbf{n}}^2 \right)^{-(d-1)}
\end{equation}
The constant $a_d=(d-2)!/(2\pi^d)$ can be obtained from the divergence theorem, as we did in the 1-dimensional case. 

The contribution to $\omega_{\mathbf{n},\bar{\mathbf{n}}}$ from a single activity at site $(0,\ldots,0)$ is then given by (App.~\ref{appendix1}): 
\begin{equation}
\label{eq:partials}
\sum_{i=1}^{d}\partial_{n_i}^2c^{*}_{\mathbf{n},\bar{\mathbf{n}}}=2 d(d-1)a_d\frac{\mathbf{n}^2-\bar{\mathbf{n}}^2}{(\mathbf{n}^2+\bar{\mathbf{n}}^2)^{d+1}}
\end{equation}
Performing calculations analogous to those for the 1~-~dimensional chain, we arrive at (App.~\ref{appendix1})
\begin{align}
\mean{\omega_{d=2}^2(r)}_\alpha&=\frac{\sigma_\alpha^2}{T^2}\frac{8}{5\pi^3r^6}\frac{1}{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}}\\
\mean{\omega_{d=3}^2(r)}_\alpha&=\frac{\sigma_\alpha^2}{T^2}\frac{27}{8\pi^4r^9}\frac{1}{\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}}
\end{align}
Importantly, we obtain exactly the same results when considering different directions across the lattice, such as $\omega_{(r,0,\ldots,0),(-r,0,\ldots,0)}$. In general, for a $d$-dimensional lattice we expect:
\begin{align}
\partial_r^2c(r)&\sim r^{-2d}\\
\omega^2_{{\rm single},d}(r)&\sim r^{-4d}/\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)\\
\langle\omega^2_d(r)\rangle_\alpha&\sim r^{-3d}/\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)
\end{align}

For completeness we investigate the behavior of $\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)$ for different dimensions. At equilibrium all cycling frequencies vanish, leading to:
\begin{equation}
\sum_{i=1}^{d}\tilde{\partial}_{n_i}^2 \overline{c}_{\mathbf{n},\bar{\mathbf{n}}}=0\quad \forall_{\mathbf{n}\neq \bar{\mathbf{n}}}
\label{eq:cond1}
\end{equation}
For all points on the diagonal of the covariance matrix $(\mathbf{n},\mathbf{n})=(n_1,\ldots,n_d;n_1,\ldots,n_d)$ the diffusion equation (see Eq.~\eqref{eq:diff}) reads
\begin{equation}
\sum_{i=1}^{d}\tilde{\partial}_{n_i}^2\overline{c}_{\mathbf{n},\mathbf{n}}+\sum_{i=1}^{d}\tilde{\partial}_{\bar{n}_i}^2\overline{c}_{\mathbf{n},\mathbf{n}}=-2
\label{eq:cond15}
\end{equation}
Using the symmetry of the system, we conclude that the two sums in Eq.~\eqref{eq:cond15} are equal, which together with Eq.~\eqref{eq:cond1} imply that
\begin{equation}
(\sum_{i=1}^{d} \tilde{\partial}_{n_i}^2)\overline{c}_{\mathbf{n},\bar{\mathbf{n}}}=-\delta_{\mathbf{n},\bar{\mathbf{n}}}
\label{eq:cond2}
\end{equation}
for all points $(\mathbf{n},\bar{\mathbf{n}})$. This result can be interpreted in the following way: for a given $(\bar{n}_1,\ldots,\bar{n}_d)$, $\overline{c}_{\mathbf{n},\bar{\mathbf{n}}}$ as a function of $(n_1,\ldots,n_d)$ is a solution to a $d$-dimensional discretised stationary diffusion equation with a single source at position $(\bar{n}_1,\ldots,\bar{n}_d)$, and with absorbing boundary conditions. Note, there is an interesting symmetry of the diffusion equation implied by the symmetry $\overline{c}_{\mathbf{n},\bar{\mathbf{n}}}=\overline{c}_{\bar{\mathbf{n}},\mathbf{n}}$: the solution at point $\mathbf{n}$ from a source at point $\bar{\mathbf{n}}$ is equal to the solution at point $\bar{\mathbf{n}}$ from a source at point $\mathbf{n}$. While this property of the diffusion equation would be obvious in an infinite space, it surprisingly holds also in the presence of absorbing boundaries. 

It can further be shown that
\begin{align}
\overline{c}_{\mathbf{n},\bar{\mathbf{n}}}\sim\ln[(n_1-\bar{n}_1)^2+(n_2-\bar{n}_2)^2]\quad \mathrm{for}\ d=2\nonumber\\
\overline{c}_{\mathbf{n},\bar{\mathbf{n}}}\sim \left( \sum_{i=1}^{d} (n_i-\bar{n}_i)^2 \right)^{-\left( \frac{d-2}{2} \right)}\ \mathrm{for}\ d> 2
\label{eq:eCscaling}
\end{align}
This result can also be understood using a simple dimensionality argument: a diffusion problem in  $d$ dimensions with a source forming a $d_s$-dimensional plane can be mapped to a $(d-d_s)$-dimensional diffusion problem with a point source. In our case we are dealing with a diffusion problem in a $2d$-dimensional space, with a $d$-dimensional source. Reducing the $2d$-dimensional problem to a point source problem in $d$ dimensions, we arrive exactly at Eq.~\eqref{eq:eCscaling}.
From this equation we conclude that for dimensions $d\geq2$, the diagonal terms of $\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)$ strongly dominate over the off-diagonal ones. In fact one can verify that for dimensions $d\geq 2$ and for systems large enough $\det\overline{\mathbf{C}}_{[{\rm r},{\rm r}]}(r)$ depends on $r$ only weakly and does not influence the scaling behavior of $\omega^2(r)$ anymore (See Fig.~\ref{fig:omega23D}). This is a consequence of the shorter range of elastic interactions in higher dimensions.

It is important to note here that, as discussed after introducing Eq.~\eqref{eq:omegalin}, the area enclosing rates do not depend on $\det\mathbf{C}_{[{\rm r},{\rm r}]}$. This allows us to perform calculations analogous to the ones presented in this section, without assuming the limit of weak activities. As a result we predict a scaling $\sigma_\alpha r^{-3d/2}$ for the area enclosing rates even for high amplitudes of the active noise.


\subsection{Generic lattices}
\label{subsec4}
In Sec.~\ref{sec:Cubic lattice} we investigated the simplest possible case of a $d$-dimensional zero-rest length cubic lattice $\mathcal{L}^d$. For such systems the Lyapunov equation for the covariance matrix could be viewed as a discretized steady state diffusion equation defined on a space $\mathcal{L}^d~\times ~\mathcal{L}^d~\sim ~\mathcal{L}^{2d}$. For instance, for a 2-dimensional square lattice we had to solve a diffusion equation in a 4-dimensional cube. A natural question is how general the connection is between the Lyapunov equation and diffusion equations. It turns out that for many zero-rest length lattices there is simple procedure for translating a particular lattice structure to a corresponding diffusion equation for the covariance. The condition which allows us to identify the terms appearing in the Lyapunov equation with second derivatives, as in Eq.~\eqref{eq:diffnc}, is that at all sites of the lattice a spring pointing in one direction is accompanied by a spring pointing in the opposite direction. If this is the case we can directly read out the diffusion equation from the structure of the lattice, as illustrated in Fig.~\ref{fig:lattice2}. Each such pair of springs gives rise to diffusive terms in the corresponding directions, with diffusion constant proportional to the spring constant.

\begin{figure}
\centering
\includegraphics[width=8cm]{lattice6}
\caption{An exemplary lattice for which the Lyaunov equation can be interpreted as a diffusion equation with non-isotropic diffusion. The presented lattice is equivalent to the triangular lattice, if we consider them in the zero-rest length case.}
\label{fig:lattice2}
\end{figure}

In the case of finite rest length elastic networks with linearized forces, the same condition  allows us to write the Lyapunov equation as a discretized second order partial differential equation for the covariance. Importantly, the displacements in $x$- and $y$-directions are no longer decoupled and one has to solve a differential equation for three different covariances: $c^{xx}$, $c^{xy}$, $c^{yy}$. An example of such an equation for a 2-dimensional triangular lattice is included in App.~\ref{appendix2}. Importantly, the structure of the network determines not only the equation for the covariance matrix but also the cycling frequencies according to Eq.~\eqref{eq:torque}.

For a wide range of networks, including randomly diluted networks \cite{Gnesotto2018}, the condition given above is not satisfied and there is no straightforward way of translating the Lyapunov equation to a continuous diffusion-like equation for the covariance. Nevertheless, for a given network $\mathcal{G}$, which can be thought of as a graph, we can still interpret the Lyapunov equation as a Poisson equation on a graph $\mathcal{G}\times\mathcal{G}$ and relate the cycling frequency between a pair of degrees of freedom to the covariance flux through a corresponding vertex of $\mathcal{G}\times\mathcal{G}$. The theory of graph Laplacians, introduced by Kirchhoff in his study of the properties of resistor networks, has found applications in elasticity theory, graph theory, and computer science~\cite{Mohar1991,Schaeffer2007,Fisher1966,Kirchhoff1847}.

Our numerical calculations reveal that the exponents of the power-laws for the cycling frequencies observed for various lattices are set by the dimensionality of the system and are independent of the detailed structure of the lattice~\cite{Mura2018}. Heuristically, this can be understood as follows: any $d$-dimensional rigid network, with a given average coordination number can be seen as an approximation to a continuous $d$-dimensional elastic medium. For such an elastic medium, a continuous diffusion equation, as the one we used to study $d$-dimensional cubic lattices, would be an exact equation for the covariance field, and we conjecture that the cycling frequencies for a continuous medium can be obtained by taking appropriate limits of our results for a discrete system. Note that a study of the cycling frequencies directly at the level of a continuous system would require introducing the Fokker-Planck equation for fields and make the analysis considerably more difficult. Since our results presented for the $d$-dimensional cubic lattice should coincide with the results for a $d$-dimensional continuous medium, we argue that our analytical calculation captures the essential origins of the power-law phenomenon for all lattices that approximate a continuous medium well.

\section{Conclusions}
Here we considered a simple model of an internally driven elastic assembly. Using this model, we investigated the properties of cycling frequencies - a two-point measure of non-equilibrium dynamics, which can be used in experimental and theoretical studies of active systems. We discussed how to relate the cycling frequencies to other commonly used non-equilibrium measures, such as the area enclosing rates or the reduced entropy production rate. Furthermore, based on our analytical approach, reinforced with numerical results, we predict that  the cycling frequencies follow a power-law as a function of distance between the two probes in an elastic network. The exponent of the power-law depends on the dimensionality of the system, but not on the detailed structure of the network. In the case of a random spatial distribution of activities, we showed that the mean cycling frequencies $\sqrt{\langle\omega^2(r)\rangle}$ are proportional to the standard deviation of the intensities of the active noise $\sigma_\alpha$. Interestingly, the case of a single-site activity gives a different exponent than the one with randomly distributed activities.


For more complex activity distributions, the connection between the Lyapunov equation and a diffusion equation, that we presented in sections~\ref{sec3} and~\ref{sec4}, provides some intuition for how the cycling frequencies in a system depend on the structure of the active noise. Since the diffusive terms in the Lyapunov equation originate solely from the structure of the lattice, we expect that a similar framework can be used to study the behavior of the cycling frequencies for more complex distributions of the active noise, which include spatial correlations~\cite{Mura2019}.

The analytical approach we developed aims at finding a mapping between the properties of the active noise and two-point non-equilibrium measures. Based on the results we obtained and their robustness to the detailed structure of a network, we argue that the cycling frequencies and the area enclosing rates are promising tools for studying the nature of the stochastic driving in an active elastic assembly. Our approach can be tested on reconstituted acto-myosin networks~\cite{Liu2006,Koenderink2009,Alvarado2017} and other noise-driven biological~\cite{Weber,Betz2009,Turlier2016} or synthetic systems~\cite{Boekhoven2015,Palacci2013,Bertrand2012}, which can be well approximated by an elastic assembly at steady-state. Such  experiments could be performed in chromosomes, membranes or tissues, using embedded colloidal particles, or fluorescently tagged cellular organelles.


\begin{acknowledgments}
We thank F.~Gnesotto, B.~Remlein, and P.~Ronceray for many stimulating discussions.
This work was supported by the German Excellence Initiative via
 the program NanoSystems Initiative Munich (NIM), the Graduate School of Quantitative Biosciences Munich (QBM),  and was funded by the Deutsche Forschungsgemeinshaft (DFG, German Research Foundation) - 418389167.
 \end{acknowledgments}


