## Hyperaccurate currents in stochastic thermodynamics

Daniel Maria Busiello

Ecole Polytechnique F´ ed´rale e de Lausanne (EPFL), Institute of Physics Laboratory of Statistical Biophysics, 1015 Lausanne, Switzerland.

## Simone Pigolotti ∗

Biological Complexity Unit, Okinawa Institute of Science and Technology Graduate University, Onna, Okinawa 904-0495, Japan.

(Dated: December 13, 2019)

Thermodynamic observables of mesoscopic systems can be expressed as integrated empirical currents. Their fluctuations are bound by thermodynamic uncertainty relations. We introduce the hyperaccurate current as the integrated empirical current with the least fluctuations in a given non-equilibrium system. For steady-state systems described by overdamped Langevin equations, we derive an equation for the hyperaccurate current by means of a variational principle. We show that the hyperaccurate current coincides with the entropy production if and only if the latter saturates the thermodynamic uncertainty relation, and it can be substantially more precise otherwise. The hyperaccurate current can be used to improve estimates of entropy production from experimental data.

Stochastic thermodynamics is a theory describing the non-equilibrium behavior of mesoscopic physical systems, from colloidal particles [1-3] to molecular motors [46]. In these systems, thermodynamic observables are stochastic quantities. A vast class of these observables can be expressed as linear functionals of the increments of a stochastic trajectory. Such observables are called integrated empirical currents . For continuous systems whose state is specified by a vector glyph[vector] x , an integrated empirical current R t ( ) (from now on simply 'current') evolves according to the dynamics [7]

equations also satisfy the same bound [11]. Interestingly, the bound of Eq. (2) does not hold for discrete-time processes [12] and looser bounds have been derived for this case [13, 14]. Thermodynamic uncertainty relations have been generalized to periodically driven systems out of steady state [15, 16] and to observables other than currents [17]. These results have been recently unified with a geometrical interpretation in the space of observables [18].

$$\frac { d R } { d t } = \vec { c } \circ \frac { d \vec { x } } { d t } \quad \quad \ \ ( 1 ) \quad \begin{array} { c } \alpha \colon \\ \text{pr} \end{array}$$

where glyph[vector] c = ( ) is a vector field that determines the curglyph[vector] c glyph[vector] x rent, and ◦ indicates the Stratonovich prescription. The total entropy production at steady state and the heat released into a thermal reservoir are examples of thermodynamic observables that can be expressed as currents.

It has been recently observed that, at steady state, all currents must satisfy the so-called thermodynamic uncertainty relation [5, 8, 9]

$$\frac { \sigma _ { R } ^ { 2 } } { \langle R \rangle ^ { 2 } } \geq \frac { 2 } { \langle S \rangle }. \quad \quad \ \ ( 2 ) \quad \ s c _ { d \epsilon }$$

The left-hand side of Eq. (2) is the coefficient of variation squared (CV ) of an arbitrary current 2 R , observed at steady state during a time t . In the right-hand side, 〈 S 〉 is the total entropy produced on average in the same time interval. Equation (2) was originally demonstrated for discrete-state systems described by master equations, first in the long-time limit [9] and later for finite times [10]. Continuous-state systems described by Langevin

∗ simone.pigolotti@oist.jp

Conceptually, the importance of Eq. (2) is that it sets a universal minimum amount of dissipation necessary to achieve currents of a given precision. Equation (2) is also of more practical interest: by seeking for currents approaching the bound, one can estimate the entropy production in a more accurate way than with other methods [19]. To this aim, it is important to know which current R approximates the bound best and how close to saturation it is. It was shown that the only current that can saturate the bound is the entropy production itself [20]. However, it is still unclear what happens when the entropy production does not saturate the bound.

In this Rapid Communication, we introduce the hyperaccurate current as the current with the lowest CV 2 in a given stochastic system. For continuous systems described by a set of overdamped Langevin equations, we derive the Euler-Lagrange equations that must be satisfied by the hyperaccurate current, and solve them in concrete examples.

We consider mesoscopic physical systems that can be described by N slow degrees of freedom glyph[vector] x = glyph[vector] x t ( ) = x 1 ( ) t , x 2 ( ) t , . . . , x N ( ). t Such degrees of freedom evolve according to a set of overdamped Langevin equations

$$\frac { d } { d t } \vec { x } = \hat { \mu } \cdot \vec { F } + \vec { \nabla } \cdot \hat { D } + \sqrt { 2 } \hat { \sigma } \cdot \vec { \xi } \text{ \quad \ \ } \text{(3)}$$

where glyph[vector] ξ = glyph[vector] ξ t ( ) = ξ 1 ( ) t , . . . , ξ N ( ) t is a Gaussian white noise with mean 〈 ξ i ( ) t 〉 = 0 and autocorrelation 〈 ξ i ( ) t ξ j ( t ′ ) 〉 = δ t ( -t ′ ) δ ij . Here the noise is interpreted

in the Ito sense. The symmetric matrix ˆ = ˆ( µ µ glyph[vector] x ) is the motility tensor and the vector glyph[vector] F = glyph[vector] F glyph[vector] ( x ) is the force acting on the system. The matrix ˆ = ˆ( σ σ glyph[vector] x ) is related to the symmetric diffusion matrix ˆ D = ˆ ( D glyph[vector] x ) by the relation ˆ σ T ˆ = σ ˆ . D We assume the Einstein relation ˆ D = k B Tµ ˆ to hold, where k B is the Boltzmann constant and T the temperature. We further assume that the matrices ˆ, σ ˆ , D and ˆ are non-degenerate. µ We associate to Eqs. (3) the Fokker-Planck equation

$$\partial _ { t } P ( \vec { x } ; t ) = \vec { \nabla } \cdot \left [ - \hat { \mu } \cdot \vec { F } P ( \vec { x } ; t ) + \hat { D } \cdot \vec { \nabla } P ( \vec { x } ; t ) \right ]. \ \ ( 4 ) \quad \text{by m}.$$

We call P st = P st ( glyph[vector] x ) the stationary solution of Eq. (4), P glyph[vector] ( x t glyph[vector] ; | y t ; ′ ) the propagator, glyph[vector] J = glyph[vector] J glyph[vector] , ( x t ) = ˆ µ · glyph[vector] FP glyph[vector] , ( x t ) -ˆ D glyph[vector] ·∇ P glyph[vector] , ( x t ) the flux, and glyph[vector] J st = glyph[vector] J st ( glyph[vector] x ) = ˆ µ FP · glyph[vector] st ( glyph[vector] x ) -ˆ D · glyph[vector] ∇ P st ( glyph[vector] x ) the stationary flux. We substitute Eqs. (3) and (4) into Eq. (1), finding an explicit evolution equation for a generic current

$$\frac { d R } { d t } = \frac { \vec { c } \cdot \vec { J } + \vec { \nabla } \cdot ( \hat { D } \cdot \vec { c } \, P ) } { P } + \sqrt { 2 } \vec { c } \cdot \hat { \sigma } \cdot \vec { \xi }. \quad ( 5 ) \quad \text{by a}$$

Equation (5) is interpreted in the Ito sense. Important examples of currents are the heat released in the thermal bath Q , with glyph[vector] c = glyph[vector] F [21], and the total entropy production S at steady state, with glyph[vector] c = ˆ D -1 · glyph[vector] J st /P st . Substituting this latter choice into Eq. (5) directly yields the evolution equation for the entropy production derived in [22].

We consider the evolution a current at steady state and use Eq. (5) to derive the uncertainty bound of Eq. (2) in a straightforward way. We introduce the bound term

$$\frac { d R _ { \text{bound} } } { d t } = \sqrt { 2 } \frac { \langle R \rangle } { \langle S \rangle } \frac { \vec { J } ^ { \text{st} } } { P ^ { \text{st} } } \cdot ( \hat { \sigma } ^ { T } ) ^ { - 1 } \cdot \vec { \xi }. \quad \ \ ( 6 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The bound term is defined so that its variance over the mean of the current squared saturates the uncertainty bound of Eq. (2), i.e.,

$$\frac { \sigma _ { R _ { \text{bound} } } ^ { 2 } } { \langle R \rangle ^ { 2 } } = \frac { 2 } { \langle S \rangle }. \quad \quad \quad ( 7 )$$

We now decompose an arbitrary current R t ( ) into the sum of the bound term and a deviation term

$$R _ { \text{dev} } ( t ) = R ( t ) - R _ { \text{bound} } ( t ). \quad \ \ ( 8 ) \quad - \,,$$

In terms of this decomposition, the left-hand side of the uncertainty bound reads

$$\frac { \sigma _ { R } ^ { 2 } } { \langle R \rangle ^ { 2 } } = \frac { \sigma _ { R _ { \text{bound} } } ^ { 2 } } { \langle R \rangle ^ { 2 } } + \frac { \sigma _ { R _ { \text{dev} } } ^ { 2 } } { \langle R \rangle ^ { 2 } } + 2 \frac { \sigma _ { R _ { \text{bound} }, R _ { \text{dev} } } ^ { 2 } } { \langle R \rangle ^ { 2 } }. \quad ( 9 ) \quad \overset { \text{perav} } { \int d \vec { y } }.$$

An explicit computation shows that the covariance σ 2 R bound ,R dev always vanishes, see SI. This implies

$$\frac { \sigma _ { R } ^ { 2 } } { \langle R \rangle ^ { 2 } } = \frac { 2 } { \langle S \rangle } + \frac { \sigma _ { R _ { \text{dev} } } ^ { 2 } } { \langle R \rangle ^ { 2 } } \geq \frac { 2 } { \langle S \rangle }. \quad \ \ ( 1 0 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

Equation (10) means that the variance of R dev is responsible for the deviation from the bound.

This calculation constitutes a short and direct demonstration of the thermodynamic uncertainty relation for a system governed by Langevin equations [11]. An advantage of this approach is to provide an explicit expression for the deviation from the bound. In particular, a current R saturates the uncertainty bound only when σ 2 R dev = 0. A necessary condition for this to hold is that the noise amplitude of R dev must vanish. Imposing this condition by means of Eqs. (5), (6), and (8) yields

$$\stackrel { \bullet \, } { \stackrel { \bullet } { \stackrel { \bullet } { \stackrel { \bullet } { \stackrel { \bullet } { \, } } } } } ( 4 ), \quad \sqrt { 2 } \left ( \vec { c } - \frac { \vec { J } ^ { \text{st} } } { P ^ { \text{st} } } \cdot \hat { D } ^ { - 1 } \frac { \langle R \rangle } { \langle S \rangle } \right ) \cdot \hat { \sigma } = 0 \quad \Leftrightarrow \quad \vec { c } \in \frac { \vec { J } ^ { \text{st} } } { P ^ { \text{st} } } \cdot \hat { D } ^ { - 1 }. \\ \stackrel { \bullet } { \, } - \hat { D } \cdot \quad \text{When $\vec{when}$ $\vec{or}$ certificate of the condition in $E_{\sigma}$ } \text{$(11)} \, \text{$\text{than}$ } R \, \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$\text$}$}$$

When glyph[vector] c satisfies the condition in Eq. (11), then R ∝ S . This means that only the entropy production, or a current proportional to it, can saturate the uncertainty bound [20]. As a corollary, if the entropy production does not saturate the bound, the bound can not be saturated by any current.

To understand such cases, we define the hyperaccurate current R h as the current with the minimum CV , among 2 all possible choices of glyph[vector] c glyph[vector] ( x ). Since σ 2 R / R 〈 〉 2 = 〈 R 2 〉 / R 〈 〉 2 -1, we seek for the hyperaccurate current by minimizing 〈 R 2 〉 / R 〈 〉 2 with respect to the function glyph[vector] c glyph[vector] ( x ).

The average value of R reads

$$\langle R \rangle = t \left \langle \frac { d R } { d t } \right \rangle = t \int d \vec { x } \ \vec { c } ( \vec { x } ) \cdot \vec { J } ^ { s t } ( \vec { x } ) \quad \ ( 1 2 )$$

where in the last equality we used Eq. (5). Similarly, we express the second moment as 〈 R 2 〉 = [ 〈 ∫ t 0 dt ′ ( dR/dt ′ )] 2 〉 . We use these expressions to evaluate the first variation of 〈 R 2 〉 / R 〈 〉 2 with respect to glyph[vector] c glyph[vector] ( x ) and impose that it must vanish (see [ ? ]). This procedure results in the EulerLagrange equation

$$\begin{array} { c } & \hat { D } ^ { - 1 } ( \vec { x } ) \cdot \vec { J } ^ { \text{st} } ( \vec { x } ) \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left \langle \frac { \vec { J } ^ { \text{st} } ( \vec { y } ) \cdot \vec { c } _ { \text{h} } ( \vec { y } ) } { P ^ { \text{st} } ( \vec { x } ) P ^ { \text{st} } ( \vec { y } ) } \right \rangle _ { y } + \\ \text{to the } & & + P ^ { \text{st} } ( \vec { x } ) \vec { \nabla } _ { \vec { x } } \left \{ \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left \langle \frac { \vec { \nabla } _ { \vec { y } } \cdot \left [ P ^ { \text{st} } ( \vec { y } ) \hat { D } ( \vec { y } ) \cdot \vec { c } _ { \text{h} } ( \vec { y } ) \right ] } { P ^ { \text{st} } ( \vec { x } ) P ^ { \text{st} } ( \vec { y } ) } \right \rangle _ { y } \right \rangle _ { y } \end{array} \right \rangle _ { y } = \\ \text{of the } & & = t \hat { D } ^ { - 1 } ( \vec { x } ) \cdot \vec { J } ^ { \text{st} } ( \vec { x } ) \frac { \langle R _ { \text{h} } ^ { 2 } \rangle } { 2 \langle R _ { \text{h} } \rangle } - t P ^ { \text{st} } ( \vec { x } ) \vec { c } _ { \text{h} } ( \vec { x } ) \\ & & \text{where } \, \vec { c } _ { \text{h} } ( \vec { x } ) \text{ is the vector field associated to the $hv$-}$$

where glyph[vector] c h ( glyph[vector] x ) is the vector field associated to the hyperaccurate current, and we denoted with 〈 . . . 〉 y = ∫ dglyph[vector] P y ( glyph[vector] x t glyph[vector] ; | y t ; ′′ ) P st ( glyph[vector] y ) the average over the initial state. In principle, also the Fano factor 〈 R 2 h 〉 / (2 〈 R h 〉 ) on the right-hand side of Eq. (13) implicitly depends on glyph[vector] c h ( glyph[vector] y ). However, we can exploit the fact that rescaling glyph[vector] c h ( glyph[vector] y ) by an arbitrary multiplicative factor does not change its CV . The solution of Eq. (13) is therefore defined up to 2 an arbitrary multiplicative constant. From now on, we shall fix this constant by setting σ 2 R h / 2 〈 R h 〉 = 1.

In the long time limit, Eq. (13) reduces to the simpler form

$$\int d \vec { y } \ \hat { K } ( \vec { x }, \vec { y } ) \cdot \vec { c } _ { h } ( \vec { y } ) = \vec { J } ^ { s t } ( \vec { x } ) \quad \quad ( 1 4 ) \quad \underset { \Phi \cdot } { \text{mol} }$$

see SI, where we defined the integral kernel

$$\hat { K } ( \vec { x }, \vec { y } ) = \frac { \vec { J } ^ { s t } ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \phi ( \vec { x }, \vec { y } ) \vec { J } ^ { s t } ( \vec { y } ) + P ^ { s t } ( \vec { x } ) \hat { D } ( \vec { x } ) \delta ( \vec { x } - \vec { y } )$$

$$- P ^ { \text{st} } ( \vec { x } ) \left [ \hat { D } ( \vec { x } ) \cdot \vec { \nabla } _ { \vec { x } } \right ] \cdot \vec { \nabla } _ { \vec { y } } \left [ \frac { \phi ( \vec { x }, \vec { y } ) } { P ^ { \text{st} } ( \vec { x } ) } \right ] \cdot \hat { D } ( \vec { y } ) P ^ { \text{st} } ( \vec { y } ) \cdot \ ( 1 5 ) \quad \ \. \ \. \.$$

and the function

$$\phi ( \vec { x }, \vec { y } ) = \int _ { 0 } ^ { + \infty } d t \, \left [ P ( \vec { x } ; t | \vec { y } ; 0 ) - P ^ { s t } ( \vec { x } ) \right ]. \quad ( 1 6 ) \, \begin{array} { c } \underset { \text{the $F$} } { \text{syste} } \\ \text{beco} \\ \text{the $F$} \end{array}$$

If the kernel ˆ ( K glyph[vector] , x glyph[vector] y ) can be inverted, then glyph[vector] c h ( glyph[vector] x ) can be expressed as

$$\vec { c } _ { h } ( \vec { x } ) = \int d \vec { y } \ \hat { K } ^ { - 1 } ( \vec { x }, \vec { y } ) \cdot \vec { J } ^ { \text{st} } ( \vec { y } ) \quad \quad ( 1 7 )$$

where ∫ dz K ˆ -1 ( glyph[vector] x, glyph[vector] z ) · ˆ ( K glyph[vector], z glyph[vector] y ) = δ glyph[vector] ( x -glyph[vector] y ).

We are now in the position to study whether the entropy production can still be hyperaccurate when it does not saturate the bound. To this aim, we assume R h ∝ S , i.e., glyph[vector] c h ∝ ˆ D -1 · glyph[vector] J st /P st and substitute this choice into Eq. (14), obtaining

$$\int d \vec { y } \, \phi ( \vec { x }, \vec { y } ) \frac { \vec { J } ^ { \text{st} } ( \vec { y } ) \cdot \hat { D } ^ { - 1 } ( \vec { y } ) \cdot \vec { J } ^ { \text{st} } ( \vec { y } ) } { P ^ { \text{st} } ( \vec { y } ) } \in P ^ { \text{st} } ( \vec { x } ). \quad ( 1 8 ) \quad \\.$$

We interpret the left hand side of Eq. (18) as the integral operator ∫ dglyph[vector] y φ glyph[vector] , ( x glyph[vector] y ) acting on the function g glyph[vector] ( y ) = glyph[vector] J st ( glyph[vector] y ) · ˆ D -1 ( glyph[vector] y ) · glyph[vector] J st ( glyph[vector] y /P ) st ( glyph[vector] y ). Such integral operator shares the same eigenfunctions of the Fokker-Planck equation (4). In particular, the stationary solution in the right-hand side of Eq. (18) is a right eigenfunction associated to a non degenerate eigenvalue equal to zero. Therefore, Eq. (18) can be satisfied only if g y ( ) ∝ P st ( glyph[vector] y ), i.e., if the quantity glyph[vector] J st · ˆ D -1 · glyph[vector] J st / P ( st ) 2 is constant. But this is precisely the condition for the entropy production to saturate the uncertainty bound [22]. We therefore conclude that, when the entropy production does not saturate the bound, it cannot be identified as the hyperaccurate current.

By definition, the CV 2 of the hyperaccurate current provides the tightest possible bound on the CV 2 of a current, the hyperaccurate bound B h . Since we set σ 2 R h / 2 〈 R h 〉 = 1, B h depends solely on the average of R h

$$\frac { \sigma _ { R } ^ { 2 } } { \langle R \rangle ^ { 2 } } \geq \mathcal { B } _ { \text{h} } = \frac { \sigma _ { R _ { \text{h} } } ^ { 2 } } { \langle R _ { \text{h} } \rangle ^ { 2 } } = \frac { 2 } { \langle R _ { \text{h} } \rangle }. \quad \quad ( 1 9 ) \quad \text{I}$$

By using Eqs. (12) and (17) to express the average of the hyperaccurate current, we obtain

$$\mathcal { B } _ { \text{h} } = \frac { 2 } { t$$

We now study the hyperaccurate current in two concrete models, where we take ˆ µ = ˆ D = ˆ for simplicI ity, with ˆ the identity I matrix. Our first example is a molecular motor in a one-dimensional periodic potential U x ( ) = sin(2 πx ) subject to a constant non-conservative force f . The system is described by the Langevin equation

$$\frac { d x } { d t } = f - \frac { d U ( x ) }$$

In this case, Eq. (15) is one dimensional. We numerically solve it by discretizing the interval [0 , 1] with a mesh ∆, so that the integral in Eq. (14) becomes a linear system of equations and the integral kernel in Eq. (15) becomes a matrix. We estimate this matrix by solving the Fokker-Planck equation numerically with the same spatial mesh ∆ (see [ ? ] for details).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00738v2-figure-eae354fa-7ebe-4865-bbc0-81a942fe5844.jpg)

glyph[arrowtp]

glyph[arrowtp]

FIG. 1. Hyperaccurate current of a molecular motor model, Eq. (21). (a) CV 2 of the hyperaccurate current and the entropy production as a function of the force f . The continuous line is the uncertainty bound of Eq. (2). Inset: Ratio between the CV 2 of the hyperaccurate current and that of the entropy production as a function of f . (b) Comparison of c x ( ) for the hyperaccurate current in red (lighter gray) and for the entropy production in blue (darker gray) for two different values of the force f , shown in the figures.

In this model, both R h and S are quite close to the bound, Fig. 1(a), with appreciable differences only for intermediate values of f (see also [22]). The CV 2 of R h is lower than that of S as predicted, although their difference is rather small [less than 1% in the range of f we considered; inset of Fig. 1(a)]. Inspecting c h ( x ), we find that it is rather similar to the one characterizing

the entropy production for low values of the force and substantially different at larger values of the force, Fig. 1(b).

As a second example, we consider the two dimensional Langevin dynamics on a torus [0 1] , × [0 , 1]:

$$\begin{array} {$$

with the non-conservative force F x ( 2 ) = f cos(2 πx 2 ). The stationary probability distribution is homogeneous, P st ( x , x 1 2 ) = 1 and the steady state flux is glyph[vector] J st ( x , x 1 2 ) = F x ( 2 ) , 0.

Since the dynamics is invariant under translations along the x 1 axis, then glyph[vector] c h ( x , x 1 2 ) = c h 1 , ( x , x 1 2 ) , c h 2 , ( x , x 1 2 ) cannot depend on x 1 . Writing Eq. (14) by components, we find that c h 2 , ( x 2 ) = 0 (see SI). Consequently, Eq. (14) reduces to the onedimensional equation in the unknown c h 1 , ( x 2 )

$$\int d y _ { 2 } \ K ( x _ { 2 }, y _ { 2 }$$

where the kernel is

$$K ( x _ { 2 }, y _ { 2 } ) & = f ^ { 2 } \$$

Since the coordinate x 2 evolves according to a simple diffusion process with periodic boundary conditions, the function φ x , y ( 2 2 ) can be explicitly expressed as

$$\phi ( x _ { 2 }, y _ { 2 } ) = \sum _ { n = 0 } ^ { + \infty } \frac { 1 } { 2 \pi ^ { 2 } n ^ { 2 } } \cos [ 2 \pi n ( x _ { 2 } - y _ { 2 } ) ] \quad ( 2 5 ) \quad \tan ^ { \cdot } { \ t a i n t } _ { \text{type} }$$

Expanding the solution c h 1 , ( x 2 ) in a Fourier basis and substituting into Eq. (23), the Fourier coefficients can be analytically calculated at any order (see [ ? ]).

In this case the CV 2 of the hyperaccurate current is much lower than that of the entropy production far from equilibrium, i.e., when f glyph[greatermuch] 1 [see Fig. 2(a)]. The hyperaccurate current converges to the entropy production when the system is near equilibrium and the bound tends to be saturated. Farther from equilibrium, the hyperaccurate current is markedly different from the entropy production, Fig. 2(b).

In this Rapid Communication, we introduced the hyperaccurate current for systems described by overdamped Langevin equations. We have shown with examples that the hyperaccurate current can be substantially more accurate than the entropy production, in cases where the latter significantly departs from the uncertainty bound. By its definition, the hyperaccurate current provides the tightest possible uncertainty bound to the CV 2 of an arbitrary current. Our theory can be extended to discretestate or discrete-time systems and possibly employed to study non integrated currents or non stationary dynamics. We leave these investigations for future work.

(

(

glyph[arrowtp]

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00738v2-figure-d37e56ee-927c-45d0-be09-e4551f414edc.jpg)



FIG. 2. Hyperaccurate current in a two-dimensional model, Eq. (22). (a) The black line is the thermodynamic uncertainty bound. The blue (top gray) line is the CV 2 of the hyperaccurate current. The red (middle gray) line is the CV 2 of the entropy production. All curves are plotted as a function of the nonconservative force f . The points represent random currents generated by adding to the coefficients κ h,n Gaussian random variables with mean zero and variance equal to f (dark-gray points) and 4 f (light-gray points). (b) Red (dark-gray) lines represent c h, 1 ( x 2 ) for different values of f . The blue line (dark-gray) represents c 1 ( x 2 ) = F x ( 2 ), whose associated current is the entropy production.

It is worthwhile discussing how the results presented here can help in estimating entropy production in experiments. Naive estimators of entropy production often require very large sample size and/or observation times to provide accurate results. Reference [19] proposes to use Eq. (2) as a tool to estimate entropy production, or at least bound it. This strategy relies on the fact that, empirically, the CV 2 of a current is much easier to estimate than the entropy production. One crucial ingredient of this strategy is to identify a current whose CV 2 is sufficiently close to the bound. Reference [19] tackles this problem by means of a Monte Carlo scheme. This approach is relatively simple to implement, but has the disadvantages of being computationally costly

and prone to overfitting, especially in high-dimensional systems. These difficulties are circumvented by the theory developed in this Rapid Communication. One possible strategy is therefore to build an approximate model of the physical system at hand, evaluate its hyperaccurate current using the theory developed in this Rapid Communication, and then measure the CV 2 of the hyperaccurate current in experiments. To pursue this strategy, it will be key to develop efficient numerical schemes [24] to solve the integral equation (14) and therefore compute the hyperaccurate current in systems more complex than the simple examples considered in

- [1] U. Seifert, Reports on progress in physics 75 , 126001 (2012).
- [2] V. Blickle, T. Speck, L. Helden, U. Seifert, and C. Bechinger, Physical review letters 96 , 070603 (2006).
- [3] I. A. Martinez, E. Rold´n, L. Dinis, ´ a and R. A. Rica, Soft matter 13 , 22 (2017).
- [4] F. J¨licher, A. Ajdari, u and J. Prost, Reviews of Modern Physics 69 , 1269 (1997).
- [5] P. Pietzonka, A. C. Barato, and U. Seifert, Physical Review E 93 , 052145 (2016).
- [6] D. M. Busiello, C. Jarzynski, and O. Raz, New Journal of Physics 20 , 093015 (2018).
- [7] R. Chetrite and H. Touchette, in Annales Henri Poincar´ e , Vol. 16 (Springer, 2015) pp. 2005-2057.
- [8] A. C. Barato and U. Seifert, Physical review letters 114 , 158101 (2015).
- [9] T. R. Gingrich, J. M. Horowitz, N. Perunov, and J. L. England, Physical review letters 116 , 120601 (2016).
- [10] J. M. Horowitz and T. R. Gingrich, Physical Review E 96 , 020103 (2017).
- [11] A. Dechant and S.-i. Sasa, Journal of Statistical Mechanics: Theory and Experiment 2018 , 063209 (2018).
- [12] N. Shiraishi, arXiv preprint arXiv:1706.00892 (2017).
- [13] K. Proesmans and C. Van den Broeck, EPL (Europhysics Letters) 119 , 20001 (2017).
- [14] D. Chiuchi` u and S. Pigolotti, Physical Review E 97 , 032109 (2018).
- [15] A. Dechant, Journal of Physics A: Mathematical and Theoretical 52 , 035001 (2018).
- [16] T. Koyuk, U. Seifert, and P. Pietzonka, Journal of Physics A: Mathematical and Theoretical 52 , 02LT02 (2018).
- [17] Y. Hasegawa and T. Van Vu, Physical Review Letters 123 , 110602 (2019).
- [18] G. Falasco, M. Esposito, and J.-C. Delvenne, arXiv preprint arXiv:1906.11360 (2019).
- [19] J. Li, J. M. Horowitz, T. R. Gingrich, and N. Fakhri, Nature communications 10 , 1666 (2019).
- [20] Y. Hasegawa and T. Van Vu, Physical Review E 99 , 062126 (2019).
- [21] K. Sekimoto, Stochastic energetics , Vol. 799 (Springer, 2010).
- [22] S. Pigolotti, I. Neri, E. Rold´n, ´ a and F. J¨licher, Physical u review letters 119 , 140604 (2017).
- [23] See Supplemental Material for additional calculations and mathematical details.

this Rapid Communication. The results of Fig. 2(a) show that, even perturbing the hyperaccurate current, one can obtain currents that are substantially more accurate than the entropy production. This supports the idea that the hyperaccurate current computed in an approximate model of a physical system can be sufficiently close to the bound to provide a reliable estimate of entropy production, if measured in an experiment.

We acknowledge D. Chiuchi´, u E. Fried, F. J¨licher, u A. Maritan, I. Neri, L. Peliti, and E Rold´ an for many ´ discussions.

[24] L. M. Delves and J. Mohamed, Computational methods for integral equations (CUP Archive, 1988).

## SUPPLEMENTARY INFORMATION

In this document, we provide additional calculations and mathematical details complementing the manuscript 'Hyperaccurate currents in stochastic thermodynamics' (from now on 'Main Text'). The document is organized as follows. In Section I, we show how to evaluate averages of stochastic integrals that we often use in the following. In Section II we show that σ 2 R bound ,R dev = 0. In Section III we derive the Euler-Lagrange equation for the vector field glyph[vector] c h ( glyph[vector] x ) defining the hyperaccurate current. The long-time limit is presented in Section IV. In Section V we discretize the integral kernel of the Euler-Lagrange equation for solving it numerically on a one-dimensional grid. Finally, in Section VI we analytically compute glyph[vector] c h ( glyph[vector] x ) for the two-dimensional model presented in the Main Text.

## I. EVALUATING TWO-TIME AVERAGES WITH THE DOOB TRANSFORM

In this Section, we show how to evaluate averages of stochastic integrals of the form

$$I = \left \langle \int \, d t ^ { \prime } \int d t ^ { \prime \prime } \, \vec { f } ( \vec { x } ( t ^ { \prime } ), t ^ { \prime }, \vec { y } ( t ^ { \prime \prime } ), t ^ { \prime \prime } ) \cdot \vec { \xi } ( t ^ { \prime \prime } ) \, \right \rangle$$

where glyph[vector] f glyph[vector] ( x t ( ′ ) , t ′ , glyph[vector] y t ( ′′ ) , t ′′ ) is an arbitrary function of an Ito process at two different times (see also [22]). For t ′′ ≥ t ′ , the average in Eq. (26) always vanishes due to the non-anticipating properties of the Wiener process in the Ito calculus. The case t ′′ &lt; t ′ requires more care. To evaluate the average in this case, we introduce the Doob transform of the process.

Doob transform maps a stochastic process conditioned on a future event to an unconditioned stochastic process with an additional drift term. In our case, we consider the Langevin equation:

$$\frac { d } { d t } \vec { x } = \hat { \mu } \cdot \vec { F } + \vec { \nabla } \cdot \hat { D } + \sqrt { 2 } \hat { \sigma } \cdot \vec { \xi }$$

and impose a future condition glyph[vector] x t ( ′ ) = glyph[vector] x ′ , with t ′ &gt; t . Doob showed that the ensemble of trajectories { glyph[vector] x } of Eq. ,(27) conditioned on the future event glyph[vector] x t ( ′ ) = glyph[vector] x ′ is equal to the unconditioned ensemble of trajectories { glyph[vector] z } generated by the Langevin equation

$$\frac { d } { d t } \vec { z } = \hat { \mu } \cdot \vec { F } + \vec { \nabla } \cdot \hat { D } + \sqrt { 2 } \hat { \sigma } \cdot \vec { \eta } ( t ) + 2 \hat { D } \cdot \vec { \nabla } _ { \vec { y } } \log P ( \vec { x } ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } )$$

where glyph[vector] η t ( ′′ ) is unbiased white noise. Comparing Eq. (27) and Eq. (28), we obtain that the conditioned averages can be transformed into unconditioned ones by substituting the noise term

$$\vec { \xi } ( t ) = \vec { \eta } ( t ) + \sqrt { 2 } \hat { \sigma } ^ { T } \vec { \nabla } _ { \vec { y } } \log P ( \vec { x } ; t ^ { \prime } | \vec { y } ; t ).$$

The second term in the right hand side represents an additional drift. Substituting this expression in Eq. (26) and using that the average of any function multiplied by glyph[vector] η vanishes, we obtain

$$I = \sqrt { 2 } \left \langle \int \, d t ^ { \prime } \int d t ^ { \prime \prime } \, \vec { f } ( \vec { x } ( t ^ { \prime } ), t ^ { \prime }, \vec { y } ( t ^ { \prime \prime } ), t ^ { \prime \prime } ) \cdot \hat { \sigma } ^ { T } \cdot \vec { \nabla } _ { \vec { y } } \log P ( \vec { x } ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) \cdot \right \rangle$$

$$\text{II.} \text{ } \text{PROOF THAT } \sigma _ { R _ { b o u n d }, R _ { d e v } } ^ { 2 } = 0$$

We now demonstrate that the covariance term σ 2 R bound ,R dev = 〈 R bound R dev 〉-〈 R bound 〉〈 R dev 〉 vanishes, where R bound and R dev are defined in Eqs. (6) and (8) of the Main Text. Since 〈 R bound 〉 = 0, we express the covariance as

$$\Gamma _ { \Lambda \Lambda \Lambda } \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambla \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambma \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambd \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lamb Da \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \LambDA \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lamboda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lambda \Lamb ] + \\ \sigma _ { R _ { \Lambda \Lambda \Lambda \Lambda \Lambda } } } ^ { 2 } _ { R _ { \Lambda \Lambda \Lambda } } } + 2 \frac { \langle R \rangle } { \langle \Sigma \rangle } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } d \vec { x } \, d \vec { y } \, \frac { \psi ( \vec { y } ) } { P ^ { s t } ( \vec { y } ) } \frac { \vec { y } ^ { s t } ( \vec { y } ) } { P ^ { s t } ( \vec { y } ) } \frac { \vec { y } ^ { s t } ( \vec { y } ) } { P ^ { s t } ( \vec { y } ) } \left [ \vec { \nabla _ { \vec { y } } } _ { \vec { y } } P ( \vec { x } t ^ { \prime } | \vec { y } t ^ { \prime \prime } ) ] P ^ { s t } ( \vec { y } )$$

where in the last term we expressed the average over the noise ξ t ( ′′ ) conditioned at future time using the Doob transform, see Section I. The quantity ψ is defined as

$$\psi ( \vec { x } ) = \vec { c } ( \vec { x } ) \cdot \vec { J } ^ { s t } ( \vec { x } ) + \vec { \nabla } _ { \vec { x } } \cdot [ \hat { D } ( \vec { x } ) \cdot \vec { c } ( \vec { x } ) P ^ { s t } ( \vec { x } ) ].$$

We now rewrite the last expression in Eq. (31) using that 〈 glyph[vector] J st · glyph[vector]/P c st 〉 = ∫ dglyph[vector] x glyph[vector] J st ( glyph[vector] x ) · glyph[vector] c glyph[vector] ( x ) = 〈 R 〉 (see Eq. (12) in the Main Text) and integrating by parts over glyph[vector] y the last term. We obtain

$$\sigma _ { R _ { b o u d }, R _ { d e v } } ^ { 2 } = - 2 \frac { \langle R \rangle ^ { 2 } } { \langle \Sigma \rangle } + 2 \frac { \langle R \rangle ^ { 2 } } { \langle \Sigma \rangle } + 2 \frac { \langle R \rangle } { \langle \Sigma \rangle } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t } d \vec { t } \, d \vec { j } \, \frac { \psi ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } P ( \vec { x } ; t ^ { \prime } | \vec { j } ; t ^ { \prime \prime } ) \vec { \nabla } _ { \vec { j } } \cdot \vec { J } ^ { s t } ( \vec { j } ) = 0$$

since glyph[vector] ∇ · glyph[vector] y glyph[vector] J st ( glyph[vector] y ) = -∂ P t st ( glyph[vector] x ) = 0 at steady state. Here and in the following, when integrating by parts we assume that the boundary term always vanish due to appropriate boundary conditions on P glyph[vector] ( x ). This result directly leads to Eq. (10) in the Main Text.

## III. DERIVATION OF THE EULER-LAGRANGE EQUATIONS

In this Section we derive the Euler-Lagrange equations for the hyperaccurate current. We want to minimize the quantity 〈 R 2 〉 / R 〈 〉 2 . Our first step is to derive tractable expressions for the first two moments of the current. The average reads

$$\langle R \rangle = t \left \langle \frac { \vec { c } \cdot \vec { J } ^ { \text{st} } } { P ^ { \text{st} } } \right \rangle$$

The second moment can be expressed as

$$\langle R ^ { 2 } \rangle = 2 t \langle \vec { c } \cdot \hat { D } \cdot \vec { c } \rangle + 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left [ \left \langle \frac { \psi ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \frac { \psi ( \vec { y } ) } { P ^ { s t } ( \vec { y } ) } \right \rangle + \sqrt { 2 } \left \langle \frac { \psi ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \vec { c } ( \vec { y } ) \cdot \hat { \sigma } ^ { T } ( \vec { y } ) \cdot \vec { \xi } ( t ^ { \prime \prime } ) \right \rangle \right ].$$

To evaluate the last term on the right hand side of Eq. (35) we again apply the Doob transform, obtaining

$$\langle R ^ { 2 } \rangle & = 2 t \langle \vec { c } \cdot \hat { D } \cdot \vec { c } \rangle + \\ & + 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d \vec { x } \, d \vec { y } \, P ( \vec { x }, t ^ { \prime } | \vec { y }, t ^ { \prime \prime } ) P ^ { s t } ( \vec { y }, t ^ { \prime \prime } ) \left [ \frac { \psi ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \frac { \psi ( \vec { y } ) } { P ^ { s t } ( \vec { y } ) } + 2 \frac { \psi ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \vec { c } ( \vec { y } ) \cdot \hat { D } ( \vec { y } ) \cdot \frac { \nabla _ { \vec { y } } P ( \vec { x }, t ^ { \prime } | \vec { y }, t ^ { \prime \prime } ) } { P ( \vec { x }, t ^ { \prime } | \$$

Integrating by parts the last term over glyph[vector] y and substituting the definition of ψ we obtain

$$\langle R ^ { 2 } \rangle & = 2 t \langle \vec { c } \cdot \hat { D } \cdot \vec { c } \rangle + \\ & + 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left ( \frac { \vec { c } ( \vec { x } ) \cdot \vec { J } ^ { s t } ( \vec { x } ) + \vec { \nabla } _ { \vec { x } } \cdot \left [ \hat { D } ( \vec { x } ) \cdot \vec { c } ( \vec { x } ) P ^ { s t } ( \vec { x } ) \right ] } { P ^ { s t } ( \vec { x } ) } \frac { \vec { c } ( \vec { y } ) \cdot \vec { J } ^ { s t } ( \vec { y } ) - \vec { \nabla } _ { \vec { y } } \cdot [ \hat { D } ( \vec { y } ) \cdot \vec { c } ( \vec { y } ) P ^ { s t } ( \vec { y } ) ] } { P ^ { s t } ( \vec { y } ) } \right > \quad \quad \quad$$

(38)

We introduce a new variable

$$\vec { \Gamma } ( \vec { x } ) = P ^ { \text{st} } ( \vec { x } ) \hat { D } ( \vec { x } ) \cdot \vec { c } ( \vec { x } ).$$

With this definition, the current is equal to the entropy production if Γ( glyph[vector] glyph[vector] x ) = glyph[vector] J st ( glyph[vector] x ). In terms of the new variable, the average of R reads

$$\langle R \rangle = t \int d \vec { x } \ \vec { \Gamma } ( \vec { x } ) \cdot \hat { D } ^ { - 1 } \cdot \frac { \vec { J } ^ { s t } ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) }$$

and, using Eq. (37), the second moment can be expressed as

$$\langle R ^ { 2 } \rangle = \int d \vec { x } \, d \vec { y } \, \left [ \vec { \Gamma } ( \vec { x } ) \cdot \hat { K } ^ { ( 1 ) } \cdot \vec { \Gamma } ( \vec { y } ) + \vec { \nabla } _ { \vec { x } } \cdot \vec { \Gamma } ( \vec { x } ) K ^ { ( 2 ) } \vec { \nabla } _ { \vec { y } } \cdot \vec { \Gamma } ( \vec { y } ) + \vec { \nabla } _ { \vec { x } } \cdot \vec { \Gamma } ( \vec { x } ) \vec { K } ^ { ( 3 ) } \cdot \vec { \Gamma } ( \vec { y } ) + \vec { \Gamma } ( \vec { x } ) \cdot \vec { K } ^ { ( 4 ) } \vec { \nabla } _ { \vec { y } } \cdot \vec { \Gamma } ( \vec { y } ) \right ] \, ( 4 1 )$$

with

$$\hat { K } ^ { ( 1 ) } & = \frac { 2 \hat { D } ^ { - 1 } ( \vec { x } ) } { P ^ { \text{st} } ( x ) } t \delta ( \vec { x } - \vec { y } ) \\ & + 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } P ( x ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) P ^ { \text{st} } ( \vec { y } ) \frac { \hat { D } ^ { - 1 } ( \vec { x } ) \cdot \vec { J } ^ { \text{st} } ( \vec { x } ) } { P ^ { \text{st} } ( \vec { y } ) ^ { 2 } } \frac { \vec { J } ^ { \text{st} } ( \vec { y } ) \cdot \hat { D } ^ { - 1 } ( \vec { y } ) } { P ^ { \text{st} } ( \vec { x } ) ^ { 2 } } \\ K ^ { ( 2 ) } & = - 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } P ( x ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) P ^ { \text{st} } ( \vec { y } ) \frac { 1 } { P ^ { \text{st} } ( \vec { x } ) P ^ { \text{st} } ( \vec { y } ) } \\ \vec { K } ^ { ( 3 ) } & = 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } P ( x ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) P ^ { \text{st} } ( \vec { y } ) \frac { 1 } { P ^ { \text{st} } ( \vec { x } ) } \frac { \vec { J } ^ { \text{st} } ( \vec { y } ) \cdot \hat { D } ^ { - 1 } ( \vec { y } ) } { P ^ { \text{st} } ( \vec { x } ) } \\ \vec { K } ^ { ( 4 ) } & = - 2 \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } P ( x ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) P ^ { \text{st} } ( \vec { y } ) \frac { \vec { J } ^ { \text{st} } ( \vec { x } ) \cdot \hat { D } ^ { - 1 } ( \vec { x } ) } { P ^ { \text{st} } ( \vec { x } ) ^ { 2 } } \frac { 1 } { P ^ { \text{st} } ( \vec { y } ) } \\ \text{that } \hat { K } ^ { ( 1 ) } & \text{ is a tensor, } K ^ { ( 2 ) } \text{ is a scalar, and } \vec { K } ^ { ( 3 ) }, \, \vec { K } ^ { ( 4 ) } \text{ are vectors. } We now compute the first variation respect
\hat { K } ^ { ( 1 ) } \quad \text{$\cdot$} \quad \text{$\cdot$}$$

Note that ˆ K (1) is a tensor, K (2) is a scalar, and glyph[vector] K (3) , glyph[vector] K (4) are vectors. We now compute the first variation respect to the function Γ( glyph[vector] glyph[vector] x )

$$\delta \left ( \frac { \langle R ^ { 2 } \rangle } { \langle R \rangle ^ { 2 } } \right ) = \frac { 1 } { \langle R \rangle ^ { 2 } } \frac { \delta \langle R ^ { 2 } \rangle } { \delta \vec { \Gamma } } \delta \vec { \eta } - \frac { 2 \langle R ^ { 2 } \rangle } { \langle R \rangle ^ { 3 } } \frac { \delta \langle R \rangle } { \delta \vec { \Gamma } } \delta \vec { \eta }$$

where while

$$\frac { \delta \langle R ^ { 2 } \rangle } { \delta \vec { \Gamma } } & = \int d \vec { x } d \vec { y } \left [ \vec { \eta } ( \vec { x } ) \cdot \hat { K } ^ { ( 1 ) } \cdot \vec { \Gamma } ( \vec { y } ) + \vec { \Gamma } ( \vec { x } ) \cdot \hat { K } ^ { ( 1 ) } \cdot \vec { \eta } ( \vec { y } ) \\ & - \vec { \eta } ( \vec { x } ) \cdot \left ( \vec { \nabla } _ { \vec { x } } K ^ { ( 2 ) } \right ) \vec { \nabla } _ { \vec { y } } \cdot \vec { \Gamma } ( \vec { y } ) - \vec { \nabla } _ { \vec { x } } \cdot \vec { \Gamma } ( \vec { x } ) \left ( \vec { \nabla } _ { \vec { y } } K ^ { ( 2 ) } \right ) \cdot \vec { \eta } ( \vec { y } ) \\ & - \vec { \eta } ( \vec { x } ) \cdot \vec { \Gamma } ( \vec { y } ) \left ( \vec { \nabla } _ { \vec { x } } \cdot \vec { K } ^ { ( 3 ) } \right ) + \vec { \nabla } _ { \vec { x } } \cdot \vec { \Gamma } ( \vec { x } ) \vec { K } ^ { ( 3 ) } \cdot \vec { \eta } ( \vec { y } ) \\ & + \vec { \eta } ( \vec { x } ) \cdot \vec { K } ^ { ( 4 ) } \vec { \nabla } _ { \vec { y } } \cdot \vec { \Gamma } ( \vec { y } ) - \vec { \Gamma } ( \vec { x } ) \cdot \vec { \eta } ( \vec { y } ) \left ( \vec { \nabla } _ { \vec { y } } \cdot \vec { K } ^ { ( 4 ) } \right ) \right ]$$

where we already integrated by parts. We reorganize this expression by swapping the variables x and y in the appropriate terms and noting that K (1) and K (2) are self-adjoint, whereas K (3) T = -K (4) . Here the superscript ' T ' denote the transposed operator, i.e., the operator obtained by swapping x and y . This results in

$$\frac { \delta ( R ^ { 2 } ) } { \delta \vec { \Gamma } } \cdot \delta \vec { \eta } = \int d \vec { x } \ \vec { \eta } ( \vec { x } ) \cdot \vec { G } ( \vec { x }, t )$$

with

$$\vec { G } ( \vec { x }, t ) = 2 \, \int d \vec { y } \, \hat { K } ^ { ( 1 ) } \cdot \vec { \Gamma } ( \vec { y } ) - \left ( \vec { \nabla } _ { \vec { x } } K ^ { ( 2 ) } \right ) \vec { \nabla } _ { \vec { y } } \cdot \vec { \Gamma } ( \vec { y } ).$$

Imposing from Eqs. (43) that the first variation vanishes for any choice of η and using Eqs. (44), and (46) leads to the condition

$$\vec { G } ( \vec { x }, t ) = 2 t \hat { D } ^ { - 1 } ( \vec { x } ) \cdot \frac { \vec { J } ^ { s t } ( \vec { x } ) } { P ^ { s t } ( \vec { x } ) } \frac { \langle R ^ { 2 } \rangle } { \langle R \rangle }$$

Substituting the expression of glyph[vector] G glyph[vector] , ( x t ) and expressing the equation in terms of glyph[vector] c h ( glyph[vector] x ) we finally obtain Eq. (13) in the Main Text.

$$\frac { \delta \langle R \rangle } { \delta \vec { \Gamma } } \delta \vec { \eta } = t \int d \vec { x } \ \vec { \eta } ( \vec { x } ) \cdot \hat { D } ^ { - 1 } \cdot \frac { \vec { J } ^ { \text{st} } ( \vec { x } ) } { P ^ { \text{st} } ( \vec { x } ) }$$

## IV. LONG-TIME LIMIT

We now derive the long-time limit of Eq. (13) in the Main Text, that we rewrite as

$$\Lambda \, \Phi ( \vec { x } ) \cdot \vec { c } _ { \mathfrak { h } } ( \vec { x } ) \cdot \vec { c } _ { \mathfrak { h } } ( \vec { x } ) + \frac { P ^ { \mathfrak { s } t } ( \vec { x } ) } { t } \hat { D } ( \vec { x } ) \cdot \vec { \nabla } _ { \vec { z } } \left \{ \int d \vec { y } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \, \frac { P ( \vec { x } ; t ^ { \prime } | j j ; t ^ { \prime \prime } ) } { P ^ { \mathfrak { s } t } ( \vec { x } ) } \vec { \nabla } _ { \vec { y } } \cdot \left [ P ^ { \mathfrak { s } t } ( \vec { y } ) \hat { D } ( \vec { y } ) \cdot \vec { c } _ { \mathfrak { h } } ( \vec { y } ) \right ] \right \} + \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \  \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \. \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \ \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \: \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ } \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \  \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \quad \, \ \quad \\ \, \quad \, \, \,$$

In the limit t → ∞ , the first two term on the left-hand side of Eq. (49) converge to a finite value, whereas the last term on the left-hand side and the right-hand side diverge as t . To avoid this divergence, we subtract from both sides the contribution glyph[vector] J st ( glyph[vector] x ) 〈 R h 〉 / 2, that we rewrite as

$$\frac { \vec { J } ^ { s t } ( \vec { x } ) ( R _ { h } ) } { 2 } = \frac { \vec { J } ^ { s t } ( \vec { x } ) } { t } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left \langle \frac { \vec { c } _ { h } \cdot \vec { J } ^ { s t } } { P ^ { s t } } \right \rangle = \frac { \vec { J } ^ { s t } ( \vec { x } ) } { t } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \int d \vec { i } \, P ^ { s t } ( \vec { x } ) \vec { c } _ { h } ( \vec { y } ) \cdot \vec { J } ^ { s t } ( \vec { y } ) \quad ( 5 0 )$$

Equation (49) then becomes

$$\cdots \quad & \cdot \cdot \cdot \cdot \\ P ^ { s t } ( \vec { x } ) \hat { D } ( \vec { x } ) \cdot \hat { c } _ { h } ( \vec { x } ) + \frac { P ^ { s t } ( \vec { x } ) } { t } \hat { D } ( \vec { x } ) \cdot \vec { \nabla } _ { \vec { z } } \left \{ \int d \vec { y } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \, \frac { P ( \vec { x } ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) } { P ^ { s t } ( \vec { x } ) } \vec { \nabla } _ { \vec { y } } \cdot \left [ P ^ { s t } ( \vec { y } ) \hat { D } ( \vec { y } ) \cdot \vec { c } _ { h } ( \vec { y } ) \right ] \right \} + \quad ( 5 1 ) \\ + \frac { \vec { J } ^ { s t } ( \vec { x } ) } { t P ^ { s t } ( \vec { x } ) } \int d \vec { y } \int _ { 0 } ^ { t } d t ^ { \prime } \int _ { 0 } ^ { t ^ { \prime } } d t ^ { \prime \prime } \left [ P ( \vec { x } ; t ^ { \prime } | \vec { y } ; t ^ { \prime \prime } ) - P ^ { s t } ( \vec { x } ) \right ] \vec { J } ^ { s t } ( \vec { y } ) \cdot \vec { c } _ { h } ( \vec { y } ) = \vec { J } ^ { s t } ( \vec { x } ) \frac { \sigma _ { R _ { h } } ^ { 2 } } { 2 ( R _ { h } ) }. \\ \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \,.$$

At this point it is safe to take the limit t →∞ . Performing the integral over time and using that at steady state the propagator depends only on the time difference t ′ -t ′′ we obtain

$$P ^ { s t } ( \vec { x } ) \hat { D } ( \vec { x } ) \cdot \vec { c } _ { \mathfrak { h } } ( \vec { x } ) + P ^ { s t } ( \vec { x } ) \left \{ \int d \vec { j } \, \hat { D } ( \vec { x } ) \cdot \vec { \nabla } _ { \vec { x } } \left [ \frac { \phi ( \vec { x }, \vec { y } ) } { P ^ { s t } ( \vec { x } ) } \right ] \vec { \nabla } _ { \vec { j } } \cdot \left [ P ^ { s t } ( \vec { y } ) \hat { D } ( \vec { y } ) \cdot \vec { c } _ { \mathfrak { h } } ( \vec { y } ) \right ] \right \} +$$

$$+ \, \frac { \vec { J } ^ { \text{st} } ( \vec { x } ) } { P ^ { \text{st} } ( \vec { x } ) } \int d \vec { y } \, \phi ( \vec { x }, \vec { y } ) \vec { J } ^ { \text{st} } ( \vec { y } ) \cdot \vec { c } _ { \text{h} } ( \vec { y } ) = \vec { J } ^ { \text{st} } ( \vec { x } ) \frac { \sigma _ { R _ { \text{h} } } ^ { 2 } } { 2 \langle R _ { \text{h} } \rangle }$$

where the function φ glyph[vector] , ( x glyph[vector] y ) is defined in Eq. (16) of the Main Text. Integrating the second term on the right hand side of Eq. (53) by parts and rewriting the whole expression as an integral equation leads to Eqs. (14) and (15) of the Main Text.

## V. DISCRETIZATION OF THE INTEGRAL KERNEL IN ONE DIMENSION

In this Section we discretize the integral kernel of the one dimensional model in the Main Text and show how to solve the Euler-Lagrange equation numerically in this case. We start by writing the expression of kernel in Eq. (17) in one dimension

$$K ( x, y ) = \frac { ( J ^ { s t } ) ^ { 2 } } { P ^ { s t } ( x ) } \phi ( x, y ) + P ^ { s t } ( x ) D ( x ) \delta ( x - y ) - D ( x ) P ^ { s t } ( x ) \nabla _ { x y } \left [ \frac { \phi ( x, y ) } { P ^ { s t } ( x ) } \right ] P ^ { s t } ( y ) D ( y )$$

where we used that the the stationary flux is constant in one dimension.

Writing explicitly the spatial derivatives results in

$$K ( x, y ) = \frac { ( J ^ { s t } ) ^ { 2 } } { P ^ { s t } ( x ) } \phi ( x, y ) + P ^ { s t } ( x ) D ( x ) \delta ( x - y ) + \left [ - \nabla _ { x y } \phi ( x, y ) + \frac { \nabla _ { x } P ^ { s t } ( x ) } { P ^ { s t } ( x ) } \nabla _ { y } \phi ( x, y ) \right ] D ( x ) P ^ { s t } ( y ) D ( y ) \ \ ( 5 5 )$$

The system is periodic in the interval [0 , 1]. We discretize this interval by introducing a mesh ∆. In this way, a function of x becomes a function of a discrete index i = glyph[floorleft] x/ ∆ , where glyph[floorright] glyph[floorleft] . . . glyph[floorright] denotes the integer part. Similarly, functions of x and y become matrices with indices i and j . In particular, the integral Eq. (14) in the Main Text becomes a linear system

$$\sum _ { j } M _ { i j } ( c _ { \mathrm h } ) _ { j } = J ^ { \mathrm s t }$$

where we call M ij the matrix obtained by discretizing the kernel of the integral equation. Such matrix reads

$$M _ { i j } = \frac { ( J ^ { s t } ) ^ { 2 } } { P _ { i } ^ { s t } } ( \phi ) _ { i j } + P _ { i } ^ { s t } D _ { i } \delta _ { i j } + \left ( - ( \nabla _ { x y } \phi ) _ { i j } + \frac { ( \nabla _ { x } P ^ { s t } ) _ { i } } { P _ { i } ^ { s t } } ( \nabla _ { y } \phi ) _ { i j } \right ) D _ { i } P _ { j } ^ { s t } D _ { j }$$

where δ ij is the Kronecker delta. We use the notation ( ∇ x P st ) , ( i ∇ y φ ) ij , and ( ∇ xy φ ) ij for the discretized derivatives, that are defined as

$$\text{at are gemneq as} & & ( \nabla _ { x } P ^ { s t } ) _ { i } = \frac { - P ^ { s t } _ { i + 2 } + 8 P ^ { s t } _ { i + 1 } - 8 P ^ { s t } _ { i - 1 } + P ^ { s t } _ { i - 2 } } { 1 2 \Delta } + \mathcal { O } ( \Delta ^ { 4 } ) \\ & ( \nabla _ { y } \phi ) _ { i j } = \frac { - \phi _ { i, j + 2 } + 8 \phi _ { i, j + 1 } - 8 \phi _ { i, j - 1 } + \phi _ { i, j - 2 } } { 1 2 \Delta } + \mathcal { O } ( \Delta ^ { 4 } ) \\ & ( \nabla _ { x y } \phi ) _ { i j } = \frac { - ( \nabla _ { y } \phi ) _ { i + 2, j } + 8 ( \nabla _ { y } \phi ) _ { i + 1, j } - 8 ( \nabla _ { y } \phi ) _ { i - 1, j } + ( \nabla _ { y } \phi ) _ { i - 2, j } } { 1 2 \Delta } + \mathcal { O } ( \Delta ^ { 4 } ) \\ \text{the 1 D example nested in the Main Text } \text{$\text{$\text{$\text{$\text{$\text{$\text{$D}(\text{$r}$} } } \, = $ 1$} \quad T$\text{ estimate the non-associate $P(\text{$r}$} \, $d(\text{$r}$} \, = $ 1$}$$

In the 1 D example presented in the Main Text, we set D x ( ) = 1. To estimate the propagator P x t y ( ; | ; 0), we numerically solve the Fokker-Planck equation associated to Eq. (21) of the Main Text

$$\partial _ { t } P ( x ; t | y ; 0 ) = - \nabla _ { x } \left [ ( \nabla _ { x } ( f x - U ( x ) ) P ( x ; t | y ; 0 ) \right ] + \nabla _ { x } ^ { 2 } P ( x ; t | y ; 0 ).$$

For the numerical integration, we approximate the initial condition P x ( ; 0 | y ; 0) = δ x ( -y ) with a Gaussian distribution with mean y and variance 3 × 10 -6 . We use the built-in solver in Mathematica with a spatial mesh ∆ = 0 002, a . time step ∆ t = 10 -5 and an 'accuracy goal' equal to half the Machine Precision (53 bits). We reach stationarity (with an error on the order of 10 -5 ) after about 10 5 time steps, i.e., at a final time t f = 1. We compute φ x, y ( ) from the propagator by integrating over time using the trapezoidal rule. Stationary flux and probability distribution are computed from the stationary solution P x t ( ; f | y ; 0) = P st ( glyph[vector] x ) + O (10 -5 ).

## VI. HYPERACCURATE CURRENT FOR THE TWO DIMENSIONAL MODEL

In this Section we derive an explicit expression for the hyperaccurate current for the two-dimensional model in Eq. (24) of the Main Text. First of all, glyph[vector] c h can not depend on x 1 because of invariance under translations along the x 1 axis. The Euler-Lagrange equations for this model then read

$$c _ { h, 1 } ( x _ { 2 } ) + f ^ { 2 } \int d y _ { 2 } \cos ( 2 \pi x _ { 2 } ) \phi ( x _ { 2 }, y _ { 2 } ) \cos ( 2 \pi y _ { 2 } ) c _ { h, 1 } ( y _ { 2 } ) = f \cos ( 2 \pi x _ { 2 } )$$

$$c _ { \mathfrak { h }, 2 } ( x _ { 2 } ) = - \int d y _ { 2 } \, \partial _ { x _ { 2 }, y _ { 2 } } \phi ( x _ { 2 }, y _ { 2 } ) \, c _ { \mathfrak { h }, 2 } ( y _ { 2 } )$$

where φ x , y ( 2 2 ) can be computed explicitly from the diffusion equation

$$\phi ( x _ { 2 }, y _ { 2 } ) = \sum _ { n = 0 } ^ { + \infty } \frac { 1 } { 2 \pi ^ { 2 } n ^ { 2 } } \cos ( 2 \pi n ( x _ { 2 } - y _ { 2 } ) )$$

Eq. (61) is solved by c h 2 , ( x 2 ) = 0, so that the hyperaccurate vector field glyph[vector] c h ( x 2 ) = ( c h 1 , ( x 2 ) , 0) is governed by the one-dimensional Eq. (60). We now expand the solution in a Fourier basis

$$c _ { h, 1 } ( y ) = \kappa _ { h, 0 } + \chi _ { h, 0 } + \sum _ { n = 1 } ^ { + \infty } \kappa _ { h, n } \cos ( 2 \pi n y ) + \sum _ { n = 1 } ^ { + \infty } \chi _ { h, n } \sin ( 2 \pi n y )$$

Inserting Eq. (63) in Eq. (60), we get

$$\kappa _ { h, 0 } + \chi _ { h, 0 } + \sum _ { n = 1 } ^ { + \infty } \kappa _ { h, n } \cos ( 2 \pi n y _ { 2 } ) + f ^ { 2 } \cos ( 2 \pi y _ { 2 } ) \sum _ { n = 2 } ^ { + \infty } \frac { 1 } { 2 \pi ^ { 2 } n ^ { 2 } } \frac { \kappa _ { h, n + 1 } + \kappa _ { h, n - 1 } } { 4 } \, \cos ( 2 \pi n y _ { 2 } ) + \frac { f ^ { 2 } } { 2 \pi ^ { 2 } } \cos ( 2 \pi y _ { 2 } ) ^ { 2 } \left ( \frac { \kappa _ { h, 2 } } { 4 } + \frac { \kappa _ { h, 0 } } { 2 } \right ) +$$

$$+ \sum _ { n = 1 } ^ { + \infty } \chi _ { h, n } \sin ( 2 \pi n y _ { 2 } ) + f ^ { 2 } \cos ( 2 \pi y _ { 2 } ) \sum _ { n = 2 } ^ { + \infty } \frac { 1 } { 2 \pi ^ { 2 } n ^ { 2 } } \frac { \chi _ { h, n + 1 } + \chi _ { h, n - 1 } } { 4 } \sin ( 2 \pi n y _ { 2 } ) + \frac { f ^ { 2 } } { 2 \pi ^ { 2 } } \frac { \chi _ { h, 2 } } { 4 } \cos ( 2 \pi y _ { 2 } ) \sin ( 2 \pi y _ { 2 } ) = f \cos ( 2 \pi y _ { 2 } ) \\ \slash e \, \in \, \real \,.$$

(64)

To make progress, we use the properties of trigonometric functions

$$\int _ { 0 } ^ { 1 } d y _ { 2 } \cos ( 2 \pi n ( x _ { 2 } - y _ { 2 } ) ) \cos ( 2 \pi y _ { 2 } ) \cos ( 2 \pi m y _ { 2 } ) = \frac { 1 } { 4 } \cos ( 2 \pi n y ) \delta _ { m, n \pm 1 } + \frac { 1 } { 4 } \cos ( 2 \pi y ) \delta _ { n, 1 } \delta _ { m, 0 }$$

$$\overset { \cdot } { \int } _ { 0 } ^ { 1 } d y _ { 2 } \cos ( 2 \pi n ( x _ { 2 } - y _ { 2 } ) ) \cos ( 2 \pi y _ { 2 } ) \sin ( 2 \pi m y _ { 2 } ) = \frac { 1 } { 4 } \sin ( 2 \pi n y ) \delta _ { m, n \pm 1 } - \frac { 1 } { 4 } \sin ( 2 \pi y ) \delta _ { n, 1 } \delta _ { m, 0 }$$

$$\int _ { 0 } ^ { 1 } d y _ { 2 } \cos ( 2 \pi n ( x _ { 2 } - y _ { 2 } ) ) \cos ( 2 \pi y _ { 2 } ) \sin ( 2 \pi m y _ { 2 } ) & = \frac { 1 } { 4 } \sin ( 2 \pi n y ) \delta _ { m, n \pm 1 } - \frac { 1 } { 4 } \sin ( 2 \pi y ) \delta _ { n, 1 } \delta _ { m, 0 } \\ \cos ( 2 \pi n y _ { 2 } ) \cos ( 2 \pi y _ { 2 } ) & = \frac { 1 } { 2 } \left ( \cos ( 2 \pi ( n + 1 ) y _ { 2 } ) + \cos ( 2 \pi ( n - 1 ) y _ { 2 } ) \right ) \\ \sin ( 2 \pi n y _ { 2 } ) \cos ( 2 \pi y _ { 2 } ) & = \frac { 1 } { 2 } \left ( \sin ( 2 \pi ( n + 1 ) y _ { 2 } ) + \sin ( 2 \pi ( n - 1 ) y _ { 2 } ) \right ) \\ \nolimits _ { \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \nolimits \nolimits.$$

Expressing the products of trigonometric functions using these relations we obtain

$$\text{Expressing the products of trigonometric functions using these relations we obtain} \\ \kappa _ { h, 0 } + \chi _ { h, 0 } + \sum _ { n = 1 } ^ { + \infty } \kappa _ { h, n } \cos ( 2 \pi n y _ { 2 } ) + f ^ { 2 } \sum _ { n = 2 } ^ { + \infty } \frac { \kappa _ { h, n + 1 } + \kappa _ { h, n - 1 } } { 1 6 \pi ^ { 2 } n ^ { 2 } } \left ( \cos ( 2 \pi ( n + 1 ) y _ { 2 } ) + \cos ( 2 \pi ( n - 1 ) y _ { 2 } ) \right ) + \\ + \frac { f ^ { 2 } } { 4 \pi ^ { 2 } } \left ( \cos ( 4 \pi y _ { 2 } ) + 1 \right ) \left ( \frac { \kappa _ { h, 2 } } { 4 } + \frac { \kappa _ { h, 0 } } { 2 } \right ) + \\ + \sum _ { n = 1 } ^ { + \infty } \chi _ { h, n } \sin ( 2 \pi n y _ { 2 } ) + f ^ { 2 } \sum _ { n = 2 } ^ { + \infty } \frac { \kappa _ { h, n + 1 } + \chi _ { h, n - 1 } } { 1 6 \pi ^ { 2 } n ^ { 2 } } \left ( \sin ( 2 \pi ( n + 1 ) y _ { 2 } ) + \sin ( 2 \pi ( n - 1 ) y _ { 2 } ) \right ) + \\ + \frac { f ^ { 2 } } { 4 \pi ^ { 2 } } \frac { \kappa _ { h, 2 } } { 4 } \left ( \sin ( 4 \pi y _ { 2 } ) + 1 \right ) = f \cos ( 2 \pi y _ { 2 } )$$

The coefficients associated to the sine vanish, χ m = 0, since on the r.h.s we have the stationary flux, which is a cosine. The non-vanishing coefficients can be written at any order

$$\text{shing coefficients can be written at any order} \\ \kappa _ { h, 0 } + \frac { f ^ { 2 } } { 4 \pi ^ { 2 } } \left ( \frac { \kappa _ { h, 2 } } { 4 } + \frac { \kappa _ { h, 0 } } { 2 } \right ) = 0 \\ \kappa _ { h, 1 } + \frac { f ^ { 2 } } { 8 \pi ^ { 2 } } \frac { \kappa _ { h, 3 } + \kappa _ { h, 1 } } { 8 } = f \\ \kappa _ { h, 2 } + \frac { f ^ { 2 } } { 1 8 \pi ^ { 2 } } \frac { \kappa _ { h, 4 } + \kappa _ { h, 2 } } { 8 } + \frac { f ^ { 2 } } { 4 \pi ^ { 2 } } \left ( \frac { \kappa _ { h, 2 } } { 4 } + \frac { \kappa _ { h, 0 } } { 2 } \right ) = 0 \\ \kappa _ { h, n } + f ^ { 2 } \frac { \kappa _ { h, n } + \kappa _ { h, n - 2 } } { 1 6 \pi ^ { 2 } ( n - 1 ) ^ { 2 } } + f ^ { 2 } \frac { \kappa _ { h, n + 2 } + \kappa _ { h, n } } { 1 6 \pi ^ { 2 } ( n + 1 ) ^ { 2 } } = 0 \quad \text{for $n>2$} \\ \text{it the second summation and changed the indices in $m=n+1$ and $m=n-1$.  The curves in Fig.2b of}$$

where we split the second summation and changed the indices in m = n +1 and m = n -1. The curves in Fig.2b of the Main Text were obtained by truncating the expansion to the 30th order. Higher order coefficients were smaller than 10 -9 for all values of the force in the explored range f ∈ [0 1 . , 100].
Daniel Maria Busiello Simone Pigolotti
Hyperaccurate currents in stochastic thermodynamics
Thermodynamic observables of mesoscopic systems can be expressed as
integrated empirical currents. Their fluctuations are bound by thermodynamic
uncertainty relations. We introduce the hyperaccurate current as the integrated
empirical current with the least fluctuations in a given non-equilibrium
system. For steady-state systems described by overdamped Langevin equations, we
derive an equation for the hyperaccurate current by means of a variational
principle. We show that the hyperaccurate current coincides with the entropy
production if and only if the latter saturates the thermodynamic uncertainty
relation, and it can be substantially more precise otherwise. The hyperaccurate
current can be used to improve estimates of entropy production from
experimental data.