Scaling behavior of non-equilibrium measures in internally driven elastic assemblies

TL;DR

This paper introduces a non-invasive method to quantify non-equilibrium activity in elastic assemblies by analyzing cycling frequencies, revealing their power-law behavior and dependence on system dimensionality and noise amplitude.

Contribution

It establishes a theoretical and numerical framework linking cycling frequencies to non-equilibrium measures in elastic networks with internal activity.

Findings

Cycling frequencies follow a power-law decay with distance.

Behavior of cycling frequencies encodes system dimensionality.

Analytical mapping predicts non-equilibrium measure behavior.

Abstract

Detecting and quantifying non-equilibrium activity is essential for studying internally driven assemblies, including synthetic active matter and complex living systems such as cells or tissue. We discuss a non-invasive approach of measuring non-equilibrium behavior based on the breaking of detailed balance. We focus on "cycling frequencies" - the average frequency with which the trajectories of pairs of degrees of freedom revolve in phase space, and explain their connection with other non-equilibrium measures, including the area enclosing rate and the entropy production rate. We test our approach on simple toy-models comprised of elastic networks immersed in a viscous fluid with site-dependent internal driving. We prove both numerically and analytically that the cycling frequencies obey a power-law as a function of distance between the tracked degrees of freedom. Importantly, the behavior of the cycling frequencies contains information about the dimensionality of the system and the amplitude of active noise. The mapping we use in our analytical approach thus offers a convenient framework for predicting the behavior of two-point non-equilibrium measures for a given activity distribution in the network.

Equations

$${\partial }_{t}p(\mathbf{x},t)=-\nabla \cdot (\mathbf{Ax}-\mathbf{D}\nabla )p(\mathbf{x},t):=-\nabla \cdot \mathbf{j}(\mathbf{x},t),$$

$$\mathbf{v}(\mathbf{x})=\mathbf{\Omega x},\mathbf{\Omega }=\mathbf{A}+\mathbf{D}{\mathbf{C}}^{-1}$$

$$\mathbf{AC}+\mathbf{C}{\mathbf{A}}^T=-2\mathbf{D}.$$

$$\mathbf{\Omega}_{\rm r}=\left(\begin{array}[]{cc}0&\omega_{ij}\\ -\omega_{ij}&0\\ \end{array}\right),$$

$$⟨{\dot{\beta }}_{ij}|{\mathbf{y}}_r⟩$$

$$=\frac{y_2⟨v_1(\mathbf{y})|{\mathbf{y}}_r⟩-y_1⟨v_2(\mathbf{y})|{\mathbf{y}}_r⟩}{y_1^2+y_2^2}={\omega }_{ij}.$$

$${\omega }_{ij}=\frac{1}{2}\frac{⟨{\mathbf{x}}_{\mathrm{r}}\times {\mathbf{f}}_{\mathrm{r}}(\mathbf{x})⟩}{\sqrt{\det {\mathbf{C}}_r}},$$

$${\omega }_{ij}=\frac{{\tilde{\partial }}_2^2{c}_{ij}}{\sqrt{\det {\mathbf{C}}_{[\mathrm{r},\mathrm{r}]}}},$$

$${\tilde{\partial }}_1^2{\overline{c}}_{ij}+{\tilde{\partial }}_2^2{\overline{c}}_{ij}$$

$${\tilde{\partial }}_1^2{c}_{ij}^{\ast }+{\tilde{\partial }}_2^2{c}_{ij}^{\ast }$$

$${\omega }_{ij}=\frac{\bar{\alpha }}{T}\frac{{\tilde{\partial }}_2^2{c}_{ij}^{\ast }}{\sqrt{\det {\overline{\mathbf{C}}}_{[\mathrm{r},\mathrm{r}]}}}+\mathcal{O}\left(\frac{{\bar{\alpha }}^2}{T^2}\right)$$

$${\omega }_{\mathrm{single}}(r)=\frac{\alpha }{T}\frac{1}{\pi r^2}\frac{1}{\sqrt{\det {\overline{\mathbf{C}}}_{[\mathrm{r},\mathrm{r}]}(r)}}+\mathcal{O}\left(\frac{{\alpha }^2}{T^2}\right)$$

$${\tilde{\partial }}_2^2{c}_{-r,r}^{\ast }(\{{\alpha }_i\})=\sum _z{\tilde{\partial }}_2^2{c}_{-r,r}^{\ast }(\{{\alpha }_i{\delta }_{iz}\})$$

$${\tilde{\partial }}_2^2{c}_{-r,r}^{\ast }(\{{\alpha }_i{\delta }_{iz}\})$$

$$=\frac{1}{\pi }\frac{{\alpha }_z}{\bar{\alpha }}\frac{rz}{(r^2+z^2{)}^2}$$

$$(\pi \bar{\alpha }{)}^2⟨({\tilde{\partial }}_2^2{c}_{-r,r}^{\ast }(\{{\alpha }_i\}){)}^2⟩$$

$$=\sum _z\frac{r^2{z}^2(⟨{\alpha }_z^2⟩-{⟨{\alpha }_z⟩}^2)}{(r^2+z^2{)}^4}\approx {\sigma }_{\alpha }^2{\int }_{-\infty }^{\infty }\frac{r^2{z}^{2}dz}{(r^2+z^2{)}^4}=\frac{\pi {\sigma }_{\alpha }^2}{16r^3}$$

$$⟨{\omega }^2(r){⟩}_{\alpha }=\frac{{\sigma }_{\alpha }^2}{T^2}\frac{1}{2\pi r^3}\frac{1}{\det {\overline{\mathbf{C}}}_{[\mathrm{r},\mathrm{r}]}(r)}.$$

$$\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}}}$$

$${\omega }_{r,\dots ,r;-r,\dots ,-r}\approx \frac{(\bar{\alpha }/T)}{\sqrt{\det {\overline{\mathbf{C}}}_{[\mathrm{r},\mathrm{r}]}}}\sum _{i=1}^d{{\tilde{\partial }}_{n_i}^2{c}_{\mathbf{n},\bar{\mathbf{n}}}^{\ast }|}_{r,\dots ,r;-r,\dots ,-r}$$

$$c_{\mathbf{n},\bar{\mathbf{n}}}^{\ast }=a_d{\left(\sum _{i=1}^d{n}_i^2+\sum _{i=1}^d{\bar{n}}_i^2\right)}^{-(d-1)}=a_d{\left({\mathbf{n}}^2+{\bar{\mathbf{n}}}^2\right)}^{-(d-1)}$$

$$\sum _{i=1}^d{\partial }_{n_i}^2{c}_{\mathbf{n},\bar{\mathbf{n}}}^{\ast }=2d(d-1)a_d\frac{{\mathbf{n}}^2-{\bar{\mathbf{n}}}^2}{({\mathbf{n}}^2+{\bar{\mathbf{n}}}^2{)}^{d+1}}$$

$${⟨{\omega }_{d=2}^2(r)⟩}_{\alpha }$$

$${⟨{\omega }_{d=3}^2(r)⟩}_{\alpha }$$

$${\partial }_r^{2}c(r)$$

$${\omega }_{\mathrm{single},d}^2(r)$$

$$⟨{\omega }_d^2(r){⟩}_{\alpha }$$

$$\sum _{i=1}^d{\tilde{\partial }}_{n_i}^2{\overline{c}}_{\mathbf{n},\bar{\mathbf{n}}}=0{\forall }_{\mathbf{n}\ne \bar{\mathbf{n}}}$$

$$\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$$

$$(\sum _{i=1}^d{\tilde{\partial }}_{n_i}^2){\overline{c}}_{\mathbf{n},\bar{\mathbf{n}}}=-{\delta }_{\mathbf{n},\bar{\mathbf{n}}}$$

$${\overline{c}}_{\mathbf{n},\bar{\mathbf{n}}}\sim \ln [(n_1-{\bar{n}}_1{)}^2+(n_2-{\bar{n}}_2{)}^2]\mathrm{for}d=2$$

$${\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$$

$$c_{\mathbf{n},\bar{\mathbf{n}}}^{\ast }=a_d{\left(\sum _{i=1}^d{n}_i^2+\sum _{i=1}^d{\bar{n}}_i^2\right)}^{-(d-1)}\mathrm{with}{a}_d=\frac{(d-2)!}{2{\pi }^d}$$

$${\partial }_{n_i}{c}_{\mathbf{n},\bar{\mathbf{n}}}^{\ast }=-(d-1)a_d\frac{2n_i}{(n_1^2+\dots +{\bar{n}}_d^2{)}^d}$$

Full text

I Introduction

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 Fodor and Cristina Marchetti (2018); Jülicher et al. (2007); Marchetti et al. (2013). A prominent class of such active matter are living systems: Schools of fishes Filella et al. (2018), flocks of birds Ballerini et al. (2008), and colonies of bacteria Damton et al. (2010); Thutupalli et al. (2015) 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 MacKintosh and Schmidt (2010); Gnesotto et al. (2018a). Examples include the stochastic fluctuations of biological assemblies such as chromosomes Weber et al. (2012), the cytoskeleton Lau et al. (2003); Guo et al. (2014); Fakhri et al. (2014); Jülicher et al. (2007); Brangwynne et al. (2008), and cellular membranes Betz et al. (2009); Turlier et al. (2016); Ben-Isaac et al. (2011). 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 Gnesotto et al. (2018a); Fodor and Cristina Marchetti (2018); Needleman and Dogic (2017); Suzuki et al. (2015); MacKintosh and Schmidt (2010); Jülicher et al. (2007). This molecular scale activity can propagate to mesoscopic scales, giving rise to non-equilibrium dynamics that breaks detailed balance Battle et al. (2016); Gladrow et al. (2016, 2017); Mura et al. (2018); Gnesotto et al. (2018a, b); Li et al. (2018) or that violates the fluctuation dissipation theorem Turlier et al. (2016); Betz et al. (2009); Mizuno et al. (2007); Martin et al. (2001); Fodor et al. (2015); Guo et al. (2014). Soft driven assemblies can also be realized in synthetic systems, including chemical fueled synthetic fibers Job Boekhoven, Wouter E. Hendriksen, Ger J. M. Koper, Rienk Eelkema (2015) and crystals of active colloidal particles Palacci et al. (2013). Numerous experimental studies showed how molecular non-equilibrium processes affect the mesoscopic mechanical properties of in vivo biological assemblies Lau et al. (2003); Guo et al. (2014); Fakhri et al. (2014); Jülicher et al. (2007), in vitro reconstituted cytoskeletal networks Liu et al. (2006); Koenderink et al. (2009); Alvarado et al. (2017) and synthetic materials Bertrand et al. (2012). 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 Seifert (2012). Recent studies made significant progress in inferring the entropy production rate from the observed trajectories Li et al. (2018); Frishman and Ronceray (2018); Seara et al. (2018). 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 Roldán et al. (2018); Mura et al. (2018); Bisker et al. (2017); Esposito (2012); Polettini and Esposito (2017). 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 Ghanta et al. (2017) and Gonzalez et al. (2018). 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 Gladrow et al. (2016, 2017) or the dynamics of a driven disordered elastic network near its isostatic point Gnesotto et al. (2018b).

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 Mura et al. (2018). 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.

II Driven elastic networks

We use overdamped networks of elastically coupled beads suspended in a viscous fluid as a simple model for soft subcellular assemblies Yucht et al. (2013); Broedersz and Mackintosh (2014); Osmanović and Rabin (2017); Mao and Lubensky (2018). 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}$ Falasco et al. (2015).

This simplified description allows us to study the dynamics of the system using a Fokker-Planck equation:

[TABLE]

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. (1) 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.

III Cycling frequencies & phase space torque

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})$ Risken (1996). 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. 2). 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 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 Weiss (2003)

[TABLE]

where $\mathbf{C}=\langle\mathbf{x}\otimes\mathbf{x}\rangle$ is the steady-state covariance matrix obeying the Lyapunov equation

[TABLE]

Equation (2) 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. (2), 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}$ Mura et al. (2018). 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 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 Weiss (2003)

[TABLE]

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. (4), we find that

[TABLE]

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 Frishman and Ronceray (2018). 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.

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 Mura et al. (2018)

[TABLE]

where $\left\langle\mathbf{x}_{\rm r}\times\mathbf{f}_{\rm r}(\mathbf{x})\right\rangle=\langle x_{i}f_{j}(\mathbf{x})-x_{j}f_{i}(\mathbf{x})\rangle=(\mathbf{C}\mathbf{A}^{T}-\mathbf{A}\mathbf{C})_{ij}$ is the mean phase space torque (Fig. 2). Intuitively, Eq. (6) 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  Mura et al. (2018). Moreover, in covariance identity coordinates Eq. (6) reduces to $\omega_{ij}=\frac{1}{2}(A^{T}-A)_{ij}=\Omega_{ij}$, in accord with previous studies Weiss (2003); Gladrow et al. (2016, 2017). 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 Ghanta et al. (2017). 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. IV).

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 Roldán et al. (2018); Mura et al. (2018); Bisker et al. (2017); Esposito (2012); Polettini and Esposito (2017). 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 Mura et al. (2018).

IV One-dimensional chain and diffusion equation

While Eq. (6) 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. (6) 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:

1. i)

Active noise present only at single site 2. ii)

Random spatial distribution of activities $\{\alpha_{i}\}$

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}$.

The first observation is that in both scenarios the cycling frequencies follow a power-law as a function of distance, as shown in Fig. 3. 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. (6) 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. (6)) reduces to:

[TABLE]

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.4A). For the 1-dimensional chain, the Lyapunov equation (see Eq. (3)) is equivalent to:

[TABLE]

At equilibrium detailed balance is preserved, which implies $\widetilde{\partial}_{2}^{2}\overline{c}_{ij}=0\ \forall_{i\neq j}$ (see Eq. (7)). We can therefore replace $\widetilde{\partial}_{2}^{2}c_{ij}$ with $\widetilde{\partial}_{2}^{2}c^{*}_{ij}$ in Eq. (7). Then, expanding the expression for the cycling frequency in powers of $(\bar{\alpha}/T)$, we get:

[TABLE]

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}\widetilde{\partial}_{2}^{2}c^{*}_{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}$.

IV.1 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. (9) 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. (9) 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_{r}c^{*}(r)=-a$ and $c^{*}(r)=-a\ln(r)+b$.

One could also argue for $\partial_{r}c^{*}(r)\sim 1/r$ scaling of the “covariance current”, by demanding that the total “covariance flux” through a circle of radius $r$, centred at [math] 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.

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. 4B). 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. V.2, where we discuss more complex networks.

We can use the approximate form of $c^{*}(r)$ together with Eq.(10) to calculate $\omega_{\rm single}(r)$. For large $N$, we can replace $\widetilde{\partial}_{2}^{2}c^{*}_{N,N+r}$ with $\partial_{r}^{2}c^{*}(r)$, to arrive at:

[TABLE]

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. (7) gives: $\widetilde{\partial}_{2}^{2}\overline{c}_{ij}=0\ \forall_{i\neq j}$. Using Eq. (8) 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. (8). 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. 5.

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.

IV.2 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. (9) (see Fig. 6). From the form of Eq. (7), 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. 6.

The connection with the steady-state diffusion equation (see Eq. (9)) 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 [math]. 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. (9) to obtain the corresponding $\mathbf{C}^{*}(\{\alpha_{i}\})$ as a superposition of single-source solutions. Thus, we can write:

[TABLE]

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:

[TABLE]

Combining Eqs. (12) and (13), 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$.

[TABLE]

In the second line of this result, we used that $\sum_{z^{\prime}\neq z}\frac{rz^{\prime}}{(r^{2}+z^{\prime 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:

[TABLE]

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. 3). 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.

V $d$-dimensional lattices

V.1 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:

[TABLE]

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. (6), we find that:

[TABLE]

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:

[TABLE]

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. A):

[TABLE]

Performing calculations analogous to those for the 1 - dimensional chain, we arrive at (App. A)

[TABLE]

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:

[TABLE]

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:

[TABLE]

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. (16)) reads

[TABLE]

Using the symmetry of the system, we conclude that the two sums in Eq. (26) are equal, which together with Eq. (25) imply that

[TABLE]

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

[TABLE]

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. (28). From this equation we conclude that for dimensions $d\geq 2$, 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. 7). 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. (10), 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.

V.2 Generic lattices

In Sec. V.1 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. (9), 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. 8. Each such pair of springs gives rise to diffusive terms in the corresponding directions, with diffusion constant proportional to the spring constant.

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. B. Importantly, the structure of the network determines not only the equation for the covariance matrix but also the cycling frequencies according to Eq. (6).

For a wide range of networks, including randomly diluted networks Gnesotto et al. (2018b), 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 Mohar, B and Alavi, Y and Chartrand, G and Oellermann, Ortrud and Schwenk (1991); Schaeffer (2007); Fisher (1966); G. Kirchhoff (1847).

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 Mura et al. (2018). 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.

VI 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 IV and V, 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 Mura et al. (2019).

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 Liu et al. (2006); Koenderink et al. (2009); Alvarado et al. (2017) and other noise-driven biological Weber et al. (2012); Betz et al. (2009); Turlier et al. (2016) or synthetic systems Job Boekhoven, Wouter E. Hendriksen, Ger J. M. Koper, Rienk Eelkema (2015); Palacci et al. (2013); Bertrand et al. (2012), 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.

Acknowledgements.

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.

Appendix A Calculation for a d-dimensional zero-restlength cubic lattice

First we find the profile of the $\mathbf{C}^{*}$ matrix for the case of a single activity at site $(0,\ldots,0)$. Neglecting the boundary conditions and assuming a rotational symmetry of the solution we find:

[TABLE]

[TABLE]

[TABLE]

Adding contributions from the all the second derivatives appearing in Eq. (17) we get:

[TABLE]

In order to get an expression for $\sum_{i=1}^{d}\partial_{x_{i}}^{2}c$ in the case of one active bead at site $(z_{1},\ldots,z_{d})$, one simply has to substitute $n_{i}\to(n_{i}-z_{i})$, $\bar{n}_{i}\to(\bar{n}_{i}-z_{i})$ in Eq. (33). Therefore, the contribution to $\omega$ from an activity $\alpha_{z_{1},\ldots,z_{d}}$ at site $(z_{1},\ldots,z_{d})$ reads:

[TABLE]

Finally, we can proceed to calculating $\left\langle\omega^{2}(2\sqrt{d}\ r)\right\rangle$.

[TABLE]

Step $(1)$ follows from $\left\langle\alpha_{z_{1},\ldots,z_{d}}\alpha_{\tilde{z}_{1},\ldots,\tilde{z}_{d}}\right\rangle=\sigma_{\alpha}^{2}\delta_{z_{1},\tilde{z}_{1}}\cdots\delta_{z_{d},\tilde{z}_{d}}$. Thereby we have assumed that $\bar{\alpha}=0$. This can be achieved by replacing the noise amplitudes $\alpha_{i}$ with $\alpha_{i}-\bar{\alpha}$. This transformation is justified, because any shift of the active noise amplitudes by a constant value does not affect $\widetilde{\partial}^{2}c^{*}_{ij}$ (compare with $\widetilde{\partial}^{2}\overline{c}_{ij}=0\ \forall_{i\neq j}$). Note that $\alpha_{i}-\bar{\alpha}$ are introduced just for convenience and one should not think of them as of any noise amplitudes. Step $(2)$ results from the fact that the terms odd in $z_{i}$ sum up to [math]. In the last step we approximated the sum by an integral.

Appendix B Calculation for finite restlength triangular lattice

Here we derive equation for the covariance matrix for the case of a finite restlength triangular lattice. We index the beads in the lattice as shown in Fig. 9

Let us denote by $f^{x}_{ij}$ and $f^{y}_{ij}$ the $x$- and $y$-components of the force acting on bead $(i,j)$, and by $x_{ij}$, $y_{ij}$ the $x$- and $y$-displacements of bead $(i,j)$. Expanding the force up to linear order in displacements we find:

[TABLE]

with $\alpha=1/4$, $\beta=\sqrt{3}/4$, $\gamma=3/4$, $\delta=\sqrt{3}/4$. It is convenient to rewrite the Lyapunov equation in the following way:

[TABLE]

Let us denote the elements of the covariance matrix by: $\left\langle x_{ij}x_{kl}\right\rangle=c^{xx}_{ij;kl}$, $\left\langle x_{ij}y_{kl}\right\rangle=c^{xy}_{ij;kl}$, $\left\langle y_{ij}y_{kl}\right\rangle=c^{yy}_{ij;kl}$ and introduce the discrete derivative operators:

[TABLE]

Then, the Lyapunov equation translates to:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

If we want to move to a continuous picture, we replace $\tilde{\partial}_{1}^{2}\to\partial_{1}^{2}$, $\tilde{\partial}_{2}^{2}\to\partial_{2}^{2}$, $\tilde{\partial}_{3}^{2}\to(\partial_{1}+\partial_{2})^{2}$. In this picture $c^{xx}$, $c^{xy}$, $c^{yx}$, $c^{yy}$ should be seen as functions on a 4-dimensional cube.

One can also write down equations for the cycling frequencies using:

[TABLE]

In the last step we used the Lyapunov equation together with the fact that $\mathbf{D}$ is a diagonal matrix.