{'fraction_non_alphanumeric': 0.04823096398102249, 'fraction_numerical': 0.014573537859285578, 'mean_word_length': 4.354140562706458, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scaling of Entropy Production under Coarse Graining in Active Disordered Media

Entropy production plays a fundamental role in the study of non-equilibrium systems by offering a quantitative handle on the degree of time-reversal symmetry breaking. It depends crucially on the degree of freedom considered as well as on the scale of description. It was hitherto unknown how the entropy production at one resolution of the degrees of freedom is related to the entropy production at another resolution. This relationship is of particular relevance to coarse grained and continuum descriptions of a given phenomenon. In this work, we derive the scaling of the entropy production under iterative coarse graining on the basis of the correlations of the underlying microscopic transition rates. Our approach unveils a natural criterion to distinguish equilibrium-like and genuinely non-equilibrium macroscopic phenomena based on the sign of the scaling exponent of the entropy production per mesostate.

Entropy production plays a fundamental role in the study of non-equilibrium systems by offering a quantitative handle on the degree of time-reversal symmetry breaking. It depends crucially on the degree of freedom considered as well as on the scale of description. It was hitherto unknown how the entropy production at one resolution of the degrees of freedom is related to the entropy production at another resolution. This relationship is of particular relevance to coarse grained and continuum descriptions of a given phenomenon. In this work, we derive the scaling of the entropy production under iterative coarse graining on the basis of the correlations of the underlying microscopic transition rates for non-interacting particles in active disordered media. Our approach unveils a natural criterion to distinguish equilibrium-like and genuinely non-equilibrium macroscopic phenomena based on the sign of the scaling exponent of the entropy production per mesostate.
Introduction. Under the umbrella term of "active matter", the study of systems driven by injection and dissipation of energy at the single-agent level has played a prominent role in the development of non-equilibrium physics and expanded its interface with biology. One of the key challenges that arise when developing models of biological matter is to quantify their degree of "nonequilibriumness", i.e. the extent to which their phenomenology differs from that of a collection of passive particles driven by a bath. The rate of entropy production [1,2] allows for such differentiation by capturing time-reversal symmetry breaking at the microscopic scale [3]. While the injection of energy ensures a strong departure from equilibrium at the single-agent level, these systems do not necessarily exhibit nonequilibrium features at larger spatio-temporal scales [4][5][6][7][8]. The question of whether equilibrium is effectively restored at this mesoscopic level requires new methods to quantify how entropy production varies under spatial coarse-graining [9][10][11][12][13][14][15]. We offer a novel perspective on this issue by studying a single-particle driven-diffusion process obtained by perturbing homogeneous diffusion with a non-conservative quenched random forcing. A similar model was recently studied as an effective description for the collective motion of active matter in a random potential [16] and can be seen as a minimal description of a molecular motor self-propelling on a disordered network of cytoskeletal filaments [17]. These models are particular examples of active disordered media [18][19][20][21], which we discuss for the first time from a thermodynamic perspective. On a more abstract level, our model may be seen as a many-particle system randomly exploring a complex phase space. From this perspective, our work determines the scaling behaviour of entropy production in a wider class of systems, including biochemical reaction networks [15].
We first analyse our model on a one-dimensional ring, where it is exactly solvable [22], and we identify a trivial scaling behaviour of the mesoscopic entropy production under block coarse-graining. We then move to higher dimensional lattices, where we show how the mesoscopic entropy production decays algebraically with block size under block coarse-graining. In order to characterise the non-trivial scaling exponents, we draw on a novel field theoretic formalism based on the Martin-Siggia-Rose construction [23,24] and demonstrate that the scaling of the entropy production can be related to the small wavenumber behaviour of the probability current's spectral density by arguments reminiscent of those employed in the treatment of hyperuniform fluctuations [25]. Our main result, Eq. (18), provides a natural criterion to distinguish between equilibrium-like and genuinely nonequilibrium macroscopic phenomena based on the sign of the scaling exponent for the entropy production per mesostate.
Entropy production and coarse graining. The starting point of our analysis is a Markovian jump process satisfying the master equatioṅ P n (t) = m (P m (t)w m,n − P n (t)w n,m ) (1) for the probability P n (t), with w n,m the non-negative transition rate from state n to state m = n. The average rate of internal entropy production at steady-state is defined as [2]Ṡ i = 1 2 n,m J n,m ln π n w n,m π m w m,n where π n = lim t→∞ P n (t) is the steady-state probability mass function and J n,m = π n w n,m − π m w m,n is the net probability current from state n to state m. The entropy productionṠ i is non-negative and vanishes for systems that satisfy detailed balance. ComputingṠ i from Eq. (2) requires complete knowledge of the set of microscopic probability currents j nm = π n w n,m , which renders this observable sensitive to time-reversal symmetry breaking induced by energy injection at the microscopic scale. However, this 'fully resolved' entropy production might be of little interest in the discussion of effective descriptions at the mesoscopic scale. In recent years, various works have addressed the issue of coarse-graining, which amounts to a partial loss of information about the microscopic currents j nm [8-10, 13, 26-28]. The perennial difficulty is that the resulting mesoscopic description is in general non-Markovian [29,30]. Esposito [9] has identified a decomposition of the full entropy production under phase-space partitioning into three non-negative contributions. Assuming a separation of timescales between intra-mesostate and inter-mesostate transitions, it was also shown that the mathematical form of the entropy production, Eq. (2), is recovered at the mesoscopic level. Previous work has focused on a single coarse-graining step, partly due to constraints of Markovianity. However, if the state space is sufficiently large, it is natural to ask whether such coarse-graining procedure could be performed iteratively, in a spirit similar to Kadanoff's "block spin" renormalisation [31]. Characterising the scale dependence of suitable observables such as the entropy production per mesostatė with N (meso) the number of mesostates at a given coarsegraining level, will then convey important information regarding the degree of activity at different scales. To carry out this programme we first denote the steady-state probability current from mesostate α to mesostate β by with L the characteristic coarse-graining length scale, such that j αβ (1) = π n w nm for α = {n} and β = {m}. We then sidestep the issue of Markovianity by postulating an effective mesoscopic entropy production of the forṁ as would be measured by any observer faithfully applying Eq.
(2) to a process observed at a given resolution L. The observableṠ (meso) i (L) should be thought of as the entropy production of a different (Markovian) process constrained to displaying the coarse-grained currents of the original process. A similar approach has been recently discussed to characterise the scaling of energy dissipation in nonequilibrium reaction networks [15].
The model. We will now introduce a minimal drivendiffusion model [32] on a regular lattice, such that the corresponding non-equilibrium steady-state is amenable to iterative block coarse-graining. Diffusion in a stable potential is the prototypical equilibrium phenomenon but there are many ways to modify the familiar diffusive dynamics into a non-equilibrium process, for example by allowing for an unstable potential [33]. An alternative modification, which guarantees a steady-state, is to impose periodic boundary conditions in such a way that a global potential function cannot be defined. To see how this is done, we recall [34] that for a Markov jump process, Eq. (1), to have an equilibrium steady state, the affinity A({w}; n 1 , . . . , n M ) = ln w n1,n2 w n2,n3 ... w n M ,n1 w n1,n M ... w n3,n2 w n2,n1 . (6) of every cycle {n 1 , n 2 , ..., n M , n 1 } has to be exactly zero. Henceforth, we will use the convention that ln(0/0) = 0 as some of the w n,m may vanish. Whenever A = 0 for any cycle, the system is intrinsically out of equilibrium and we should expect steady-state probability currents [33]. A straightforward way of inducing non-equilibrium behaviour is therefore to allow for quenched disorder in the transition rates, which can be interpreted as a nonconservative random forcing [35,36]. The resulting disordered steady-state is then non-equilibrium with exceedingly high probability. In the following we therefore consider a homogeneous diffusion process on a periodic lattice in d dimensions and allow for a quenched perturbation to the nearest-neighbour hopping rates such that with ζ n,m > −h a set of zero-mean random variables, Fig. 1. Henceforth, · will denote averages over this random variable. Ring topology (d = 1). The one-dimensional version of this model has previously been considered in the context of random walks in random environments [35,36]. Starting from an exact result by Derrida [22] for the net current, J = J n,n+1 , the entropy production, Eq. (2), can be calculated asṠ i = JA({w}; {1, 2, . . . , N }), on the basis of constant J = π n w n,n+1 − π n+1 w n+1,n and the affinity A, Eq. (6), taken for the cycle passing through all sites of the ring once [34], Suppl. Sec. ??. After substituting Eq. (7) into Eqs. (6) we eventually obtaiṅ where O(ζ k ) stands for any term proportional to ζ i1±1,i1 ...ζ i k ±1,i k with any indices i 1 , ..., i k . As can be seen by setting ζ i,j = 0 for all i, j, the entropy production vanishes at zeroth order in the perturbation, as expected. Supp. Sec. I explores the weak disorder limit of Eq. (8) in more detail. We now apply the coarse-graining procedure based on Eqs. (4) and (5). In one dimension, the interface between distinct mesostates consists of a single edge (Fig. 1) and the net current j αβ (L) − j βα (L) = J is independent of the block size L. The mesoscopic entropy productioṅ S (meso) i , Eq. (5), is thus given by a sum over a subset of the contributions to the microscopic entropy productioṅ S i , Eq. (2). By invoking translational invariance and linearity of expectation, the total entropy production rate in a system with originally N states coarse-grained into FIG. 1. Coarse-graining procedure for a diffusion process on a periodic square lattice perturbed by a spatially quenched, non-conservative disorder. In one dimension (top), the transition rates locally define a random walk with increments νn = ζn,n+1 − ζn+1,n and the affinity A({w}; {1, .., N }) of the closed cycle across all sites is to leading order proportional to the displacement of the random walk after N steps, Suppl. Sec. ??. In higher dimensions (bottom), the random potential picture breaks down locally and the coarse-graining becomes non-trivial due to the current no longer being uniform. mesostate blocks of size L is = N L J ln π n w n,n+1 π n+1 w n+1,n , with arbitrary state index n. For uncorrelated noise in the weak disorder limit and using π n /π m = r n /r m [22], Ṡ (meso) i = 2λ/(N hL)+O(λ 2 ), where λ denotes the variance of ζ n,m . Irrespective of the noise strength, the dependence of the total entropy production on the block size L is solely due to the absence of terms from currents within a microstate block. In one dimension, the entropy production per mesostate ṡ , is thus independent of L. No current averaging at interfaces between blocks takes place. The situation is qualitatively different in d > 1, as we will demonstrate now.
Periodic lattices with d > 1. In higher dimensions, the equilibrium condition A({w}; {n 1 , ..., n M }) = 0 is generally broken at the local rather than global scale. No analytical expression for the steady-state currents is available and we resort to a perturbation theory in weak disorder, based on a Martin-Siggia-Rose field theory [23,24], which allows us to extract various static correlation functions in arbitrary dimensions and for a wide class of disorders, see Suppl. Sec. ??. There, we introduce the net microscopic probability current J (x) = (J (1) (x), ..., J (d) (x)) as the continuum limit of J n,m on a hypercubic lattice, together with its Fourier transform J (k) = (J (1) (k), ..., J (d) (k)). We follow the convention with k = 2πn/(N ) (n ∈ Z d ), assuming a hypercubic system, and V = (N ) d the phase space volume, with the dimensionful lattice spacing. We assume a disorder characterised by the covariance in Fourier space for |k| → 0, corresponding, for η = 0, to indicates the disorder affecting edges parallel to the a-th dimension of the lattice (see Suppl. Sec. ?? for details of the specification of backward rates). The spectral density tensor of the probability current evaluated at tree level then reads (Suppl. Sec. ??), for k = 0, Eq. (13) matches the general form of the spectral density of a divergence-free, isotropic vector field, which is well known from the theory of turbulence of incompressible fluids [37,38]. For η < 0, the vanishing of the spectral density as k → 0 indicates that the probability current is hyperuniform [25], i.e. exhibits an anomalous suppression of fluctuations at large wavelengths. The case η = 0 corresponds to uncorrelated (white) noise, We can now carry out the coarse-graining procedure for the entropy production. First, we note that the mesoscopic entropy production, Eq. (5), is given by a sum over contributions from neighbouring mesostates, α = β in (5). By linearity of expectation, the disorder average Ṡ (meso) i is therefore the expected contribution from a single interface multiplied by the number of interfaces. It follows that, for a hypercubic lattice with coordination number 2d, where {α, β} is any pair of neighbouring mesostates. Asymptotically in large L, Eq. (14) can be approximated where denotes the variance of the probability current integrated across an interface of linear dimension L. This relation is derived in Suppl. Sec. ??. In the continuum limit, the asymptotic scaling of the entropy production per mesostate is therefore controlled by the asymptotic variance of the current integrated across the interface between states, which is in turn determined by the small wavenumber behaviour of the current spectral density tensor introduced in Eq. (13). In fact, the relationship between Eqs. (13) and (16) is exactly the type of problem addressed in the study of hyperuniform systems [25].
Using results from these studies, one obtains A more thorough derivation of these scaling laws is provided in Suppl. Sec. ??. We can think of σ 2 J (L) ∼ L d−2 as scaling with the perimeter of the interface. In this sense, it is instructive to draw a comparison with the case of a non-solenoidal random vector field with spectral density J (a) (k)J (b) (k ) = (λ/V )δ ab δ k+k ,0 , in which case the variance of the integrated current instead scales with the area of the interface, σ 2 J (L) ∼ L d−1 . The requirement that the steady-state is divergence-free thus plays an important role by imposing long-range correlations in the currents, even when these are not present in the substrate, i.e. for η = 0. Combining Eqs. (15) and (17) we eventually arrive at (Suppl. Sec. ??) which constitutes our key result.
The scaling exponent changes sign at η = 1, corresponding to , suggesting a quantitative distinction between steady-states that are increasingly "equilibrium-like" at larger scales and genuinely non-equilibrium states where dissipation occurs at all scales. Numerical simulations for η = 0 in d = 2, 3, as well as the full range −1 < η < d in d = 2 show excellent agreement with our analytical prediction (see Figure  2). Investigating entropy production numerically indicates that this scaling behaviour of the entropy production is unchanged in the strong disorder regime ( Figure  3).
Concluding remarks and outlook. Based on the effective entropy production introduced in Eq. (5), we have studied the mesoscopic scaling of the entropy production per mesostate under iterative block coarse graining for a family of non-equilibrium disordered steady-states with generic substrate correlations. We demonstrate that coarse-graining of degrees of freedom results in nontrivial scaling of the entropy production, Eqs. (18), so that its value on one scale is related to its value on another scale in a non-trivial fashion that ultimately draws on the correlations of the transition rates, Eq. (12). To characterise our model we have developed a static fieldtheoretic framework, which complements recent work [39] by considering problems involving multiplicative noise. Our main result, Eq. (18), offers conditions for which active disordered media appear equilibrium-like, η ≤ 1, or genuinely non-equilibrium, η > 1, at the large scale and in a statistical sense, i.e. when the behaviour is averaged over many realisations.

Experiments.
A natural application of the theory above of active disordered media are active particles [40] on irregular surfaces, which phenomenologically behave as randomly driven passive particles [16]. This is due to ratchet effects, i.e. local asymmetries in the potential driving a net current [41]. In this case, correlations in the medium can be induced by controlling the ruggedness of the substrate [42]. Further, in vitro experiments involv-ing active transport by molecular motors in a network of cytoskeletal filaments [43,44] often involve cell extracts where these filaments are uniformly disordered (η = 0). Non-trivial, long-range correlations in this type of systems could be induced e.g. by coupling tagged filaments to an external magnetic field, as done in [45] with actin.

ACKNOWLEDGMENTS
LC would like to thank Rosalba Garcia-Millan, Salvatore Torquato, Nikolai Antonov and Polina Kakin for fruitful discussion at various stages of this work. LC and GS acknowledge support from the Francis Crick Institute, which receives its core funding from Cancer Research UK (FC001317), the UK Medical Research Council (FC001317), and the Wellcome Trust (FC001317).
for a lattice with N sites, where The entropy production, Eq. (2), can be calculated asṠ i = JA({w}; {1, 2, . . . , N }), on the basis of constant J = π n w n,n+1 − π n+1 w n+1,n and the affinity A, Eq. (6), taken for the cycle passing through all sites of the ring once [2]. Eq. (7) can be substituted into Eqs. (6) and (S1) to obtain where O(ζ k ) stands for any term proportional to ζ i1±1,i1 ...ζ i k ±1,i k with any indices i 1 , ..., i k . This is the expression we refer to in the main text, Eq. (8). In this regime, the current is dominated by the average of the approximate potential gradients, (ζ n,n+1 − ζ n+1,n )/h. The physical picture provided by Eq. (S3) shows that the weak disorder expansion, |ζ n,m | h, is equivalent to a linear response theory, J ∝ A, with the cycle affinity playing the role of the thermodynamic force [3]. The behaviour away from the weak disorder limit is more subtle. In particular, the current J and thereforeṠ i are no longer a function only of the affinity but depend on the whole random potential profile. This is due to the presence of traps, i.e. local minima of the random potential, at which probability tends to concentrate [4].
So far, we have not made any assumption about the autocorrelation of the noise ζ. However, by relating Eq. (S4) to the square of the distance travelled by a random walker with zero-mean increments ν n = ζ n,n+1 − ζ n+1,n , Fig. 1, we conclude that the system-size scaling of the disorder expectation Ṡ i in the large N limit is controlled by the Hurst exponent H ∈ (0, 1) [5] of the random walk. In particular, where H = 1/2 corresponds to uncorrelated (or short-range correlated) disorder. By the central limit theorem for sufficiently uncorrelated ν n , and more generally when N n=1 ν n is a Gaussian random variable, the affinity is also Gaussian. AsṠ i , Eq. (S4), is essentially the square of the affinity, its distribution is of the chi-squared type. These results are indeed confirmed by simulations, Fig. S1. arXiv:2109.05288v2 [cond-mat.stat-mech] 18 Mar 2022 FIG. S1. In one dimension and for weak disorder, the scaling of the entropy production with system size is controlled by the Hurst exponent according to Eq. (S5). Top: the exponent δ obtained by fitting numerical estimates of Ṡ i against N δ for different Hurst exponents H ∈ (0, 1), shown with error bar in blue, is in perfect agreement with the analytical prediction, δ = 2H − 2, Eq. (S5), shown as a black dashed line. Bottom: the distribution of the entropy productionṠi as a random variable is well described by the chi-squared distribution, shown as a black dashed line, as predicted by the theory.

II. A STATIC PATH INTEGRAL APPROACH TO DIFFUSION WITH QUENCHED NOISE
In the following we develop the field theoretic formalism used to study the correlation function of the steady-state probability currents for the model of diffusion on a lattice with spatially quenched, non-conservative noise in the hopping rates. This formalism is based on the well-known response-field construction by Martin-Siggia-Rose (MSR) and Janssen-De Dominicis [6,7] but extends it in two ways: first, our approach deals with multiplicative noise in a static framework, complementing recent work by Antonov and co-workers [8], which focused on the case of additive noise. Second, we show how the MSR 'trick' of imposing physically relevant relations between the noise field and physical observables by means of suitable resolutions of identity can be exploited to directly probe observables that depend explicitly on the noise (in this case, the probability current). The formalism is first developed in one dimension, d = 1, for the case of white noise and it is subsequently generalised to higher integer dimensions and correlated noise. It turns out that we ultimately do not need the full field theory nor the renormalisation group, because the upper critical dimension d c = 2 coincides with the lowest dimension we draw on the field theory. Nevertheless, the framework provides us with a firm basis to reliably determine the upper critical dimension and will form the well-tested foundation for the analysis of more advanced problems.
The analytical expressions derived in the following are in agreement with numerical results, as shown in Fig. S1 for the one-dimensional problem and in Figs. 2 and S2 for its generalisation to higher dimensional lattices. In one dimension, the steady state probability mass function π n for a given noise realisation {ζ m,m±1 } is solved exactly by means of Kirchhoff's theorem [2] and the homogeneous current is obtained straightforwardly as J = π n w n,n+1 − π n+1 w n1,n for an arbitrary n. The entropy production is then calculated asṠ i = JA on the basis of Eq. (6). In higher dimensions, a numerical approximation to the steady state probability is instead obtained by evolving a homogeneous initial condition, P n (t = 0) = N −d , according to the master equation, Eq. (1).
A. One dimensional case, uncorrelated noise.
We consider a one-dimensional diffusion process on a ring of size N characterised by a homogeneous hopping rate h and a quenched perturbation ζ i,i±1 to the hopping rate from site i to site i ± 1 satisfying ζ i,i±1 + h > 0. We denote by φ i the steady-state probability mass function at site i, with i = 1, 2, · · · , N . By normalisation, i φ i = 1. In the following we apply periodic boundary conditions to all indices i, so that i = 0 is equivalent to i = N and i = N + 1 is equivalent to i = 1. The steady-state probability mass function is determined implicitly by the master equation where for the second equality we have redefined the perturbation to the hopping rates according to The right hand side of Eq. (S6) can be written more suggestively in continuum notation as a drift-diffusion equation, where the symmetric part of the quenched noise, denoted by δ i , appears both as an inhomogeneous perturbation of the diffusion constant h and, in the drift term, as a random forcing. The antisymmetric part of the noise, ζ i , on the other hand, only affects the drift term. For simplicity, we will henceforth assume that the perturbation to the transition rate matrix satisfies The lattice master equation (S6) thus reduces to Numerical investigation indicates that this assumption does not modify the scaling behaviour of the mesoscopic entropy production, Figure S2. For weak disorder, this observation can be rationalised by arguing that the diffusive behaviour is dominated by its homogeneous component, while the drift term is dominated at large spatial scales by the term of lowest order in the spatial derivatives. Importantly, for i.i.d. rates, the steady-state remains generically non-equilibrium, which can be verified by computing the cycle affinity, Eq. (6), and exhibits non-zero steady-state currents. We further define the probability current field as the net current from site i to site i + 1, such that and introduce a symmetrised, local, average current In one dimension at stationarity, J i,i+1 = J i−1,i = J i and there is therefore no difference between the local averaged currents and the more microscopic ones. The non-negativity of the transition rates means that ζ i + h ≥ 0. However, in order to be able to cast Eqs. (S10) and (S11) into a path integral form, we will assume that ζ i is a weak Gaussian noise with variance λ and probability This assumption is justified a posteriori by comparing the analytical predictions in the weak disorder limit with numerical simulations where the positivity of the transition rates is confirmed at initialisation. Introducing the short-hand notation for the functional integral measure . Moreover, while out analytical approach relies on a weak disorder limit to allow for rate perturbation to be Gaussian without running into unphysical negative rates, numerical experiments suggest that the scaling laws we obtain remain valid in the strong disorder regime (independent of whether (S9) is imposed), which is explored by considering a noise on the transition rates that is homogeneously distributed in the range |ζ| ∈ (1 − , 1] for various choices of the variance (right panel for d = 2 and η = 0). Theoretical predictions for the asymptotic scaling behaviour are plotted for reference (black curves in both panels).

the expectation of a general noise-dependent observable O[ζ] can be written as the path integral
To impose the steady-state Eq. (S10) at every site, we define and introduce the following resolution of identity written in terms of Dirac delta functions fixing the functional relation between the probability mass function φ i and the noise ζ i To go from Eq. (S19) to Eq. (S20) we have used the Fourier representation of the Dirac delta function, whereφ i is a purely imaginary auxiliary field, sometimes referred to as the response field. While M is independent of φ for a general master equation and can thus be taken outside of the φ integral, the Jacobian retains a dependence on the random variables ζ 1 , ..., ζ N due to the multiplicative nature of the latter. While the sign of the auxiliary field determining that of the exponent in Eq. (S20) is arbitrary, our choice gives the conventional sign for the response propagator [7]. The determinant M[ζ] can be computed by means of Faddeev-Popov ghosts [6,7], a standard procedure which requires the introducing two new Grassmann fields, denoted ξ i andξ i , satisfying the anti-commutation relations {ξ i , ξ j } = {ξ i ,ξ j } = {ξ i , ξ j } = 0. Using the Grassmann path integral representation of the determinant [7], This path integral has the same structure as that in Eq. (S20) because differentiating f i with respect to φ j effectively returns the Markov matrix of the process, now acting to the right on the Grassmann field ξ j . The Jacobian (S23) couples the Faddeev-Popov ghosts ξ i andξ i to the noise field ζ i , which is in turn coupled to the probability mass function φ i and its responseφ i via Eq. (S20). Integrating out the noise will thus result in a coupling between the ghosts and the MSR fields, as we will demonstrate after introducing the current field J i and its responseJ i .
In principle, the current J i is fully defined in Eq. (S13). However, this definition ceases to be useful once the noise field ζ i has been integrated out, rendering J i intractable. To avoid this, we introduce a resolution of the identity similar to Eq. (S20) to enforce the definition Eq. (S13) of the probability current J i . It reads where the trivial Jacobian M ≡ 1 has been omitted. Inserting Eqs. (S20), (S23) and (S25) into the Gaussian path integral (S16) we obtain Since the noise ζ enters the exponent only linearly or quadratically, the functional integral over ζ is Gaussian and can be performed in closed form. In particular, isolating all terms on the right-hand side of Eq. (S26) that depend on ζ, In the first equality we have re-indexed the summand, applying periodic boundary conditions, so that the noise field always appears with the same index and the completion of squares producing the last line can be carried out straightforwardly. We can now rewrite Eq. (S26) as where we have split the action into its bilinear part, denoted A 0 , and an interacting part, denoted A int , which includes all higher order terms in products of fields. Denoting by the discrete first and second spatial derivatives, respectively, these contributions are given explicitly by and where h.o.t. denotes higher order terms in the spatial derivatives. In the following we are interested in the large scale (small wavenumber) features of the density and current correlation function, and we will thus ignore contributions We proceed by introducing Fourier-summed fields in the non-linear action. A typical interaction term in the Fourier representation of the action reads The field φ i represents a probability mass function and it is therefore natural to study fluctuations around the normalised homogeneous steady-state φ i = N −1 , so that Φ k=0 = 1. Instead of implementing Φ k=0 = 1 through the action along the lines of Eqs. (S10) and (S13), we perform a change of variable φ i → N −1 + φ i at the level of Eqs. (S33) and (S34), which leaves the bilinear action A 0 unchanged but generates a number of new terms in A int . For example, with γ = λ, γ = 2N −1 λ and γ = N −2 λ at bare level. The full shifted action reads, again at bare level, The various coupling are more easily identifiable in the diagrammatic notation that we shall introduce shortly.

B. Propagators and diagrammatics
The bare propagators of the theory can then be extracted directly from the bilinear action, Eq. (S33). They read Φ kΦk = 1 4h sin 2 (k/2) + r N δ k+k ,0 = (S40) where we have allowed for a mass r for the purpose of infrared regularisation, which would enter the bilinear action A 0 of Eq. (S33) via terms of the form rφ i φ i and rξ i ξ i . All physical observables are evaluated in the limit of vanishing mass r since the theory is massless as a matter of conservation of probability. The diagrammatic notation for the amputated interaction vertices is with dashed propagators denoting spatial derivatives. At bare level, γ = κ = 2λ and σ = π = χ = 4λ. Vertices appearing in the same column of the list (S44) are generated from the same interaction term of Eq. (S34) upon performing the shift φ i → φ i + N −1 . Their coupling are thus related via γ /γ = κ /κ = σ /σ = N −2 , γ /γ = κ /κ = σ /σ = 2N −1 and π /π = χ /χ = N −1 .

C. Correlation functions in one dimension
We are now ready to calculate the density and current correlation functions for the one-dimensional model. These are expressed as a power series in the small disorder strength λ, which we truncate to first order (tree level) to allow for direct comparison with the result obtained by expanding the analytical results by Derrida [1], Eq. (S3). This approximation is controlled in d = 1 because we are interested in system of finite size N and will be justified on relevance grounds in higher dimensions, where we will work in the continuum limit. The correlation function for the probability mass function at tree level is given by for k = 0. For the current we have with Kronecker delta δ k ,0 reflecting the constraint that J i = J j = J is homogeneous. For k = 0, the three diagrams contributing to the right-hand side of Eq. (S46) only differ in their symmetry factors and overall sign, which produces the desired cancellation. When k = 0, only the first diagram, which does not involve spatial derivatives, contributes and the cancellation does not occur. Since J k=0 = N J, we thus have that J 2 = N −2 J 0 J 0 = 4λN −3 , recovering the leading order contribution obtained by expanding Derrida's [1] exact result, Eq. (S3), with the antisymmetric condition ζ i,i+1 − ζ i,i−1 = 2ζ i ,
where J (a) i is once again the symmetric local average at site i of the probability current along the ath dimension of the lattice. The functional relationship between each component of the current vector field and the other fields, Eq. (S51), is then imposed by means of d resolutions of the identity of the form (S25). In the case of uncorrelated noise, the derivation of the action follows precisely the same lines as above and generates d copies of the type of terms we have already seen, each with spatial derivatives taken with respect to a different dimension. Explicitly, cf. Eq. (S33) and cf. Eq. (S34). Since vertices arising from Eq. (S53) involve vector fields, the intepretation of diagrams needs to be modified slightly. For example, the mixed propagators in Fourier representation becomes Similarly, Using the field theory in d > 1 to characterise correlations is most covenient in the continuum limit. We thus take the limit N → ∞ at fixed volume V = (N ) d , with the dimensionful lattice spacing. Based on Eqs. (S31) and (S32) we define the differential operators in the continuum ∇ d = ∇ and ∆ d = 2 ∆, together with the fields φ andφ now being defined for all x ∈ (0, N ] d , such that where the dimension of φ(x) is pre-determined by the shift performed above, which in the continuum reads φ(x) → V −1 + φ(x), while the exponent β is undetermined for the time being. Converting sums over lattice sites into integrals over space according to and demanding that the action of Eqs. (S52) and (S53) is dimensionless term by term, defines the remaining fields and the dimensionful couplings of the continuum theory D = h 2−β andλ = d+2−2β λ. The exponent β is fixed by imposing that the diffusion constant D remains finite in the limit → 0, whence β = 2. Since nowλ = λ d−2 , we identify d c = 2 as the upper critical dimension of the theory. For d > 2 all interactions become irrelevant and we expect the large scale behaviour of the theory to be well described by the tree-level diagrammatics, thus justifying the truncation of higher order terms. Precisely at d = 2, renormalisation group theory [6] predicts logarithmic corrections to tree-level scaling but these were not observed conclusively in our numerical investigation.
In the continuum we thus have the action A 0 + A int with Absorbing the shift V −1 of φ into new couplings as done using dashed variables in Eq. (S44) means that these differ now from each other in their engineering dimension, for example Maintaining a finite volume V is important for the validity of the present theory, but it also means that Fourier transforming it by mean of Eq. (11) results in sums over suitable modes and δ-functions of the form We may occasionally approximate sums over momenta by integrals. To make the notation more suggestive of continuous k ∈ R d we write the Fourier transformed fields as Φ(k), J (k), ξ(k) etc. The continuum propagators are structurally identical to their discrete counterpart and read Generalising the formalism to allow for correlations in the transition rates is relatively straightforward. What changes is the form of the probability functional for the noise, originally Eq. (S14), which should now include a nontrivial dependence on the power spectrum. The power spectrum Q(k) characterises the noise correlation function for the different components of the disorder in Fourier space as The case of uncorrelated noise covered above corresponds to Q(k) =λ, i.e. to the case where Q(k) is independent of k. The noise probability functional corresponding to Eq. (S66) reads where the sample space of ζ(k) is constrained to ζ(k) = ζ * (−k) as to maintain ζ(x), the inverse Fourier transform of ζ(k), being real. The role of Q(k) is to penalise spatial fluctuations of ζ(k) depending on its wavenumber. Since our setup is rotationally symmetric we assume the general algebraic form of the correlator to be Q(k) ∼λk −η for small k, where k = |k|. For η > 0, the right hand side of Eq. (S67) is regularised by excluding the homogeneous mode, k = 0, from the summation. The derivation of the field theory with non-trivial correlator proceeds as outline above up until Eq. (S27), which has now a slightly different form, where we have introduced the short-hand F (k) for the Fourier transform of the real expression We thus find that the generalisation to correlated noise characterised by a power spectrum Q(k) =λk −η amounts to augmenting the interaction terms (in Fourier space) of the original theory by a factor k −η . The Fourier representation, originally Eq. (S37), of the interaction vertices to be used in our subsequent calculation is upgraded to the its continuum form for general η according to with γ = κ = 2λ/V 2 and σ = 4λ/V 2 . The effect of this modification on the spectral density of the currents becomes apparent already at tree level, where we now obtain, for k = 0, This is the result we invoke in the main text, Eq. (13). The special case of uncorrelated rates is recovered at η = 0. The momentum-dependent factor appearing in brackets in Eq. (S76) is the fingerprint of a solenoidal vector field [9,10]. For k = 0 and η ≤ 0, the second and third diagrams vanish and we instead find consistently with Eq. (S46).

III. LARGE L ASYMPTOTIC FORM OF EQ. (18)
We start from Eq. (14) for the expectation of the entropy production per mesostate with linear block dimensions L, where {α, β} is any pair of neighbouring mesostates. Without loss of generality we assume that the interface is a (d − 1)-dimensional hypercube of constant x coordinate and denote {i α } the set of "boundary" states of α connected to β by a single edge. Denoting e = (1, 0, · · · , 0) the d-dimensional unit displacement vector along the x axis, let {i α + e} correspond to the set of "boundary" states of β connected to α by a single edge. Inspecting the logarithmic factor on the right-hand side of Eq. (14), where we have expanded the steady-state probability about its mean, π i =π + δπ i , withπ = N −d . Each logarithm in the expression above contains a sum over L d−1 zero-mean random variables hδπ i + π i ζ i,i+e and hδπ i+e + π i+e ζ i+e,i , respectively. While the expected value of each sum vanishes by linearity, the typical magnitude can be estimated by computing its standard deviation. We thus write down the variance of the sum The first contribution on the right hand side of Eq. (S80), associated with the sum of variances, scales like L −(d−1) due to statistical homogeneity and can thus be taken to be small compared to the factor 1 appearing in the logarithm Eq. (S79) in the large L regime. The second contribution is associated with cross correlations among the summands. It involves L d−1 (L d−1 − 1) terms so it is at most of order L 0 . However, this scaling only applies to infinite-range correlations and the more relevant decay of density as well as noise correlations over large distances will in general produce a scaling L ν with ν < 0. The second contributions is thus also asymptotically small compared to the factor 1 appearing in the logarithm. For sufficiently large L we can therefore expand Eq. (S79) to leading order and obtain ln j αβ (L) j βα (L) (S82)

IV. ASYMPTOTIC SCALING OF INTEGRATED CURRENTS FROM CORRELATION FUNCTION
Here we outline the calculation required to compute the asymptotic variance of the integrated current across a mesostate interface, Eq. (16), starting from the spectral density tensor of the probability current vector field, Eq. (13). These arguments closely follow the standard treatment of hyperuniform fluctuations, reviewed in [11]. Without loss of generality we assume that the interface, denoted Ω henceforth, is a d − 1 dimensional hypersurface of constant x = 0 coordinate embedded in d dimensional space (see Figure S3). We introducer as the d − 1 dimensional vector satisfying r = (r x ,r). The variance of the integrated current across the interface Ω is thus given by is the overlap surface fraction, with I Ω the indicator function over the interface, and v Ω is the interface surface. What we require to proceed further is therefore an expression for the component C xx of the current correlation function for displacements on the interface. To obtain this, we start from the corresponding component of the spectral density, Eq. (S76), and introducek as the d − 1 dimensional vector satisfying k = (k x ,k) to write The correlation function C xx is thus obtained from Eq. (S87) by Fourier back-transforming according to the convention Eq. (11), where we have approximated sums over momenta, which involve an infrared cutoff at a typical inverse length-scale V −1/d , as integrals in the limit of large V on the basis of the infrared convergence of the integral Eq. (S91) for η < d. Following the usual convention [6], dashed differentials are a short-hand for d d k = d d k/(2π) d . It is interesting to note that, while the full spectral density Eq. (13) is not hyperuniform for η = 0, the relevant component of the corresponding correlation function is hyperuniform for η < 1 with hyperuniformity exponent 1−η when we concentrate on the lower-dimensional object that is the mesostate interface. For the correlation function C xx (0,r) origination from Eq. (S91) to be well-defined in the range η < d, the Fourier back-transform needs to be regularised in the ultraviolet.
Here, we do so by introducing a soft cutoff, suggestive of the introduction of a microscopic lattice spacing of typical length scale β = . Denoting C R xx (0,r; β) the regularised correlation function, we write for η < d and d > 2, where 2 F 1 (a, b, c; z) denotes the hypergeometric function [12]. We note for later use that 2 F 1 (a, b, c; z) has the asymptotic expansion at large |z| unless a − b is an integer [12]. To go from Eq. (S92) to (S93) we first perform a change of variable to spherical coordinates using for n = d − 1 and n > 1, withk = |k| the wavenumber and θ = arccos[(k/|k|) · (r/|r|)] the angle between the wavevectork and the displacementr. Since the integrand of Eq. (S92) is a function only ofk and θ, the remaining angular coordinates can be integrated out. To go from Eq. (S93) to (S94) we perform a second change of variable u(θ) = cos(θ) and expand the ensuing trigonometric functions as power series inku|r| before carrying out the integral over u. In the following, we shall ignore the numerical prefactor in Eq. (S93), since is independent of both θ and k and thus irrelevant for the window size scaling. The case d = 2 needs to be treated separately and produces for η < d: which is characterised by the same leading order asymptotic scaling with |r| at |r|/β 1 as Eq. (S95) for d → 2. Taking the interface as a hypersphere of radius L = L , whence v Ω ∼ L d−1 , it was shown [11] that, for |r| < 2L, (S102) The term in Eq. (S102) originating from the zeroth order in the series expansion of γ approaches, in the asymptotic limit of L → ∞, the value of the spectral density, Eq. (S91), atk = 0, lim L→∞ Ω d d−1r C xx (0,r) = lim k→0 2λΓ 1+η 2 V 2 √ πΓ 1 + η 2 |k| 1−η (S103) and thus vanishes when the latter is hyperuniform, namely for η < 1. Based on the expansion Eq. (S96), the behaviour of the correlation function, Eq. (S95), for |r|/β 1 is C R xx (0,r; β) ∼ |r| η−d therefore each term in the integrand in the right-hand side of Eq. (S102) behaves like |r| η−2+m with m ∈ {0, 1, 3, ...}. The scaling of the variance with window size L is thus controlled by the infrared divergence as L → ∞. For η < 1 − m, the integrals are infrared convergent, which allows us to let Ω → R d−1 without modifying the leading order asymptotic scaling with L. Proceeding case by case: • 0 < η < d: For η > 1, the relevant integrals are infrared divergent. For 0 < η < 1 the term originating from the zeroth order in the expansion of γ is removed by the cancellation Eq. (S103). For η = 1, based on Eq. (S91), the correlation function C xx takes the form of a Dirac delta in real space so that the term originating from the zeroth order in the expansion of γ is finite but L independent and does not contribute to the scaling. Overall, the scaling is dominated by the infrared divergent behaviour at large |r|. Substituting Eq. (S95) together with the asymptotic expansion (S96) into Eq. (S102) and performing a change of variabler →r/L to isolate the L dependence, we arrive at Var(J; Ω) ∼ L d−2+η for 0 < η < d.
• η = 0: The case of η = 0 needs to be treated separately because the integral of the first order term in the expansion of γ(r; Ω) is logarithmically divergent. Substituting Eq. (S95) for the regularised correlation function into Eq. (S102) and setting η = 0 we are left with integrals of the form (S104) for m ≥ 1. For m > 1 and based on the expansion Eq. (S96), ther integral is dominated by the behaviour of the integrand in the large |r| regime and we recover the L d−2 scaling obtained upon setting η = 0 in the scaling law for the case 0 < η < d. For m = 1 we exploit the spherical symmetry of the integral to write and invoke the asymptotic expansion Eq. (S96) to obtain the leading order behaviour which dominates the asymptotic scaling with L. The variance of the integrated current thus scales like Var(J; Ω) ∼ L d−2 log(L/β) for η = 0.