Hyperaccurate currents in stochastic thermodynamics

TL;DR

This paper introduces the concept of hyperaccurate currents in stochastic thermodynamics, which are the integrated empirical currents with minimal fluctuations, providing improved estimates of entropy production in non-equilibrium systems.

Contribution

The authors derive a variational principle for hyperaccurate currents in overdamped Langevin systems and analyze their relation to entropy production and thermodynamic uncertainty relations.

Findings

Hyperaccurate current coincides with entropy production when the uncertainty relation is saturated.

Hyperaccurate currents can significantly reduce fluctuations compared to standard currents.

The approach improves entropy production estimates from experimental data.

Abstract

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.

Equations

$$\frac{dR}{dt}=\vec{c}\circ \frac{d\vec{x}}{dt}$$

$$\frac{{\sigma }_R^2}{⟨R{⟩}^2}\ge \frac{2}{⟨S⟩}.$$

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

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

$$\frac{dR}{dt}$$

$$\frac{d{R}_{\mathrm{bound}}}{dt}=\sqrt{2}\frac{⟨R⟩}{⟨S⟩}\frac{{\vec{J}}^{\mathrm{st}}}{P^{\mathrm{st}}}\cdot ({\hat{\sigma }}^T{)}^{-1}\cdot \vec{\xi }.$$

$$\frac{{\sigma }_{R_{\mathrm{bound}}}^2}{⟨R{⟩}^2}=\frac{2}{⟨S⟩}.$$

$$R_{\mathrm{dev}}(t)=R(t)-R_{\mathrm{bound}}(t).$$

$$\frac{{\sigma }_R^2}{⟨R{⟩}^2}=\frac{{\sigma }_{R_{\mathrm{bound}}}^2}{⟨R{⟩}^2}+\frac{{\sigma }_{R_{\mathrm{dev}}}^2}{⟨R{⟩}^2}+2\frac{{\sigma }_{R_{\mathrm{bound}},R_{\mathrm{dev}}}^2}{⟨R{⟩}^2}.$$

$$\frac{{\sigma }_R^2}{⟨R{⟩}^2}=\frac{2}{⟨S⟩}+\frac{{\sigma }_{R_{\mathrm{dev}}}^2}{⟨R{⟩}^2}\ge \frac{2}{⟨S⟩}.$$

$$\sqrt{2}\left(\vec{c}-\frac{{\vec{J}}^{\mathrm{st}}}{P^{\mathrm{st}}}\cdot {\hat{D}}^{-1}\frac{⟨R⟩}{⟨S⟩}\right)\cdot \hat{\sigma }=0\iff \vec{c}\propto \frac{{\vec{J}}^{\mathrm{st}}}{P^{\mathrm{st}}}\cdot {\hat{D}}^{-1}.$$

$$⟨R⟩=t⟨\frac{dR}{dt}⟩=t\int d\vec{x}\vec{c}(\vec{x})\cdot {\vec{J}}^{\mathrm{st}}(\vec{x})$$

$${\hat{D}}^{-1}(\vec{x})\cdot {\vec{J}}^{\mathrm{st}}(\vec{x}){\int }_0^{t}d{t}^{\prime }{\int }_0^{t^{\prime }}d{t}^{\prime \prime }{⟨\frac{{\vec{J}}^{\mathrm{st}}(\vec{y})\cdot {\vec{c}}_{\mathrm{h}}(\vec{y})}{P^{\mathrm{st}}(\vec{x})P^{\mathrm{st}}(\vec{y})}⟩}_y+$$

$$\int d\vec{y}\hat{K}(\vec{x},\vec{y})\cdot {\vec{c}}_{\mathrm{h}}(\vec{y})={\vec{J}}^{\mathrm{st}}(\vec{x})$$

$$\hat{K}(\vec{x},\vec{y})=\frac{{\vec{J}}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})}\varphi (\vec{x},\vec{y}){\vec{J}}^{\mathrm{st}}(\vec{y})+P^{\mathrm{st}}(\vec{x})\hat{D}(\vec{x})\delta (\vec{x}-\vec{y})$$

$$-P^{\mathrm{st}}(\vec{x})\left[\hat{D}(\vec{x})\cdot {\vec{\nabla }}_{\vec{x}}\right]\cdot {\vec{\nabla }}_{\vec{y}}\left[\frac{\varphi (\vec{x},\vec{y})}{P^{\mathrm{st}}(\vec{x})}\right]\cdot \hat{D}(\vec{y})P^{\mathrm{st}}(\vec{y}).$$

$$\varphi (\vec{x},\vec{y})={\int }_0^{+\infty }dt\left[P(\vec{x};t|\vec{y};0)-P^{\mathrm{st}}(\vec{x})\right].$$

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

$$\int d\vec{y}\varphi (\vec{x},\vec{y})\frac{{\vec{J}}^{\mathrm{st}}(\vec{y})\cdot {\hat{D}}^{-1}(\vec{y})\cdot {\vec{J}}^{\mathrm{st}}(\vec{y})}{P^{\mathrm{st}}(\vec{y})}\propto P^{\mathrm{st}}(\vec{x}).$$

$$\frac{{\sigma }_R^2}{⟨R{⟩}^2}\ge {\mathcal{B}}_{\mathrm{h}}=\frac{{\sigma }_{R_{\mathrm{h}}}^2}{⟨R_{\mathrm{h}}{⟩}^2}=\frac{2}{⟨R_{\mathrm{h}}⟩}.$$

$${\mathcal{B}}_{\mathrm{h}}=\frac{2}{t}{\left(\int d\vec{x}d\vec{y}{\vec{J}}^{\mathrm{st}}(\vec{y})\cdot {\hat{K}}^{-1}(\vec{y},\vec{x})\cdot {\vec{J}}^{\mathrm{st}}(\vec{x})\right)}^{-1}.$$

$$\frac{dx}{dt}=f-\frac{dU(x)}{dx}+\sqrt{2}\xi .$$

$$\frac{d{x}_1}{dt}$$

$$\frac{d{x}_2}{dt}$$

$$\int d{y}_{2}K(x_2,y_2)c_{\mathrm{h},1}(y_2)=f\cos (2\pi x_2)$$

$$K(x_2,y_2)=f^2\cos (2\pi x_2)\varphi (x_2,y_2)\cos (2\pi y_2)+\delta (x_2-y_2).$$

$$\varphi (x_2,y_2)=\sum _{n=0}^{+\infty }\frac{1}{2{\pi }^2{n}^2}\cos [2\pi n(x_2-y_2)]$$

$$I=⟨\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 })⟩$$

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

$$\frac{d}{dt}\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 })$$

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

Full text

††preprint: APS/123-QED

Stochastic thermodynamics is a theory describing the non-equilibrium behavior of mesoscopic physical systems, from colloidal particles Seifert (2012); Blickle et al. (2006); Martinez et al. (2017) to molecular motors Jülicher et al. (1997); Pietzonka et al. (2016); Busiello et al. (2018). 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 $\vec{x}$, an integrated empirical current $R(t)$ (from now on simply “current”) evolves according to the dynamics Chetrite and Touchette (2015)

[TABLE]

where $\vec{c}=\vec{c}(\vec{x})$ is a vector field that determines the current, and $\circ$ 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 Barato and Seifert (2015); Pietzonka et al. (2016); Gingrich et al. (2016)

[TABLE]

The left-hand side of Eq. (2) is the coefficient of variation squared ($\mathrm{CV}^{2}$) of an arbitrary current $R$, observed at steady state during a time $t$. In the right-hand side, $\langle S\rangle$ 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 Gingrich et al. (2016) and later for finite times Horowitz and Gingrich (2017). Continuous-state systems described by Langevin equations also satisfy the same bound Dechant and Sasa (2018). Interestingly, the bound of Eq. (2) does not hold for discrete-time processes Shiraishi (2017) and looser bounds have been derived for this case Proesmans and Van den Broeck (2017); Chiuchiù and Pigolotti (2018). Thermodynamic uncertainty relations have been generalized to periodically driven systems out of steady state Dechant (2018); Koyuk et al. (2018) and to observables other than currents Hasegawa and Van Vu (2019). These results have been recently unified with a geometrical interpretation in the space of observables Falasco et al. (2019).

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 Li et al. (2019). 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 Hasegawa and Van Vu (2019b). 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 $\mathrm{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 $\vec{x}=\vec{x}(t)=x_{1}(t),x_{2}(t),\dots,x_{N}(t)$. Such degrees of freedom evolve according to a set of overdamped Langevin equations

[TABLE]

where $\vec{\xi}=\vec{\xi}(t)=\xi_{1}(t),\dots,\xi_{N}(t)$ is a Gaussian white noise with mean $\langle\xi_{i}(t)\rangle=0$ and autocorrelation $\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=\delta(t-t^{\prime})\delta_{ij}$. Here the noise is interpreted in the Ito sense. The symmetric matrix $\hat{\mu}=\hat{\mu}(\vec{x})$ is the motility tensor and the vector $\vec{F}=\vec{F}(\vec{x})$ is the force acting on the system. The matrix $\hat{\sigma}=\hat{\sigma}(\vec{x})$ is related to the symmetric diffusion matrix $\hat{D}=\hat{D}(\vec{x})$ by the relation $\hat{\sigma}^{T}\hat{\sigma}=\hat{D}$. We assume the Einstein relation $\hat{D}=k_{B}T\hat{\mu}$ to hold, where $k_{B}$ is the Boltzmann constant and $T$ the temperature. We further assume that the matrices $\hat{\sigma}$, $\hat{D}$, and $\hat{\mu}$ are non-degenerate. We associate to Eqs. (3) the Fokker-Planck equation

[TABLE]

We call $P^{\mathrm{st}}=P^{\mathrm{st}}(\vec{x})$ the stationary solution of Eq. (4), $P(\vec{x};t|\vec{y};t^{\prime})$ the propagator, $\vec{J}=\vec{J}(\vec{x},t)=\hat{\mu}\cdot\vec{F}P(\vec{x},t)-\hat{D}\cdot\vec{\nabla}P(\vec{x},t)$ the flux, and $\vec{J}^{\mathrm{st}}=\vec{J}^{\mathrm{st}}(\vec{x})=\hat{\mu}\cdot\vec{F}P^{\mathrm{st}}(\vec{x})-\hat{D}\cdot\vec{\nabla}P^{\mathrm{st}}(\vec{x})$ the stationary flux. We substitute Eqs. (3) and (4) into Eq. (1), finding an explicit evolution equation for a generic current

[TABLE]

Equation (5) is interpreted in the Ito sense. Important examples of currents are the heat released in the thermal bath $Q$, with $\vec{c}=\vec{F}$ Sekimoto (2010), and the total entropy production $S$ at steady state, with $\vec{c}=\hat{D}^{-1}\cdot\vec{J}^{\mathrm{st}}/P^{\mathrm{st}}$. Substituting this latter choice into Eq. (5) directly yields the evolution equation for the entropy production derived in Pigolotti et al. (2017).

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

[TABLE]

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.,

[TABLE]

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

[TABLE]

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

[TABLE]

An explicit computation shows that the covariance $\sigma^{2}_{R_{\mathrm{bound}},R_{\mathrm{dev}}}$ always vanishes, see SI. This implies

[TABLE]

Equation (10) means that the variance of $R_{\mathrm{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 Dechant and Sasa (2018). 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 $\sigma^{2}_{R_{\mathrm{dev}}}=0$. A necessary condition for this to hold is that the noise amplitude of $R_{\mathrm{dev}}$ must vanish. Imposing this condition by means of Eqs. (5), (6), and (8) yields

[TABLE]

When $\vec{c}$ satisfies the condition in Eq. (11), then $R\propto S$. This means that only the entropy production, or a current proportional to it, can saturate the uncertainty bound Hasegawa and Van Vu (2019b). 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_{\mathrm{h}}$ as the current with the minimum $\mathrm{CV}^{2}$, among all possible choices of $\vec{c}(\vec{x})$. Since $\sigma^{2}_{R}/\langle R\rangle^{2}=\langle R^{2}\rangle/\langle R\rangle^{2}-1$, we seek for the hyperaccurate current by minimizing $\langle R^{2}\rangle/\langle R\rangle^{2}$ with respect to the function $\vec{c}(\vec{x})$.

The average value of $R$ reads

[TABLE]

where in the last equality we used Eq. (5). Similarly, we express the second moment as $\langle R^{2}\rangle=\langle[\int_{0}^{t}dt^{\prime}(dR/dt^{\prime})]^{2}\rangle$. We use these expressions to evaluate the first variation of $\langle R^{2}\rangle/\langle R\rangle^{2}$ with respect to $\vec{c}(\vec{x})$ and impose that it must vanish (see SI). This procedure results in the Euler-Lagrange equation

[TABLE]

where $\vec{c}_{\mathrm{h}}(\vec{x})$ is the vector field associated to the hyperaccurate current, and we denoted with $\langle\dots\rangle_{y}=\int d\vec{y}P(\vec{x};t|\vec{y};t^{\prime\prime})P^{\mathrm{st}}(\vec{y})$ the average over the initial state. In principle, also the Fano factor $\langle R_{\mathrm{h}}^{2}\rangle/(2\langle R_{\mathrm{h}}\rangle)$ on the right-hand side of Eq. (13) implicitly depends on $\vec{c}_{\mathrm{h}}(\vec{y})$. However, we can exploit the fact that rescaling $\vec{c}_{\mathrm{h}}(\vec{y})$ by an arbitrary multiplicative factor does not change its $\mathrm{CV}^{2}$. The solution of Eq. (13) is therefore defined up to an arbitrary multiplicative constant. From now on, we shall fix this constant by setting $\sigma^{2}_{R_{\mathrm{h}}}/2\langle R_{\mathrm{h}}\rangle=1$.

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

[TABLE]

see SI, where we defined the integral kernel

[TABLE]

and the function

[TABLE]

If the kernel $\hat{K}(\vec{x},\vec{y})$ can be inverted, then $\vec{c}_{\mathrm{h}}(\vec{x})$ can be expressed as

[TABLE]

where $\int dz~{}\hat{K}^{-1}(\vec{x},\vec{z})\cdot\hat{K}(\vec{z},\vec{y})=\delta(\vec{x}-\vec{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_{\mathrm{h}}\propto S$, i.e., $\vec{c}_{\mathrm{h}}\propto\hat{D}^{-1}\cdot\vec{J}^{\mathrm{st}}/P^{\rm st}$ and substitute this choice into Eq. (14), obtaining

[TABLE]

We interpret the left hand side of Eq. (18) as the integral operator $\int d\vec{y}~{}\phi(\vec{x},\vec{y})$ acting on the function $g(\vec{y})=\vec{J}^{\mathrm{st}}(\vec{y})\cdot\hat{D}^{-1}(\vec{y})\cdot\vec{J}^{\mathrm{st}}(\vec{y})/P^{\rm st}(\vec{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)\propto P^{\mathrm{st}}(\vec{y})$, i.e., if the quantity $\vec{J}^{\mathrm{st}}\cdot\hat{D}^{-1}\cdot\vec{J}^{\mathrm{st}}/(P^{\rm st})^{2}$ is constant. But this is precisely the condition for the entropy production to saturate the uncertainty bound Pigolotti et al. (2017). We therefore conclude that, when the entropy production does not saturate the bound, it cannot be identified as the hyperaccurate current.

By definition, the $\mathrm{CV}^{2}$ of the hyperaccurate current provides the tightest possible bound on the $\mathrm{CV}^{2}$ of a current, the hyperaccurate bound $\mathcal{B}_{\mathrm{h}}$. Since we set $\sigma^{2}_{R_{\mathrm{h}}}/2\langle R_{\mathrm{h}}\rangle=1$, $\mathcal{B}_{\mathrm{h}}$ depends solely on the average of $R_{\mathrm{h}}$

[TABLE]

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

[TABLE]

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

[TABLE]

In this case, Eq. (15) is one dimensional. We numerically solve it by discretizing the interval $[0,1]$ with a mesh $\Delta$, 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 $\Delta$ (see SI for details).

In this model, both $R_{\mathrm{h}}$ and $S$ are quite close to the bound, Fig. 1(a), with appreciable differences only for intermediate values of $f$ (see also Pigolotti et al. (2017)). The $\mathrm{CV}^{2}$ of $R_{\mathrm{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_{\mathrm{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]\times[0,1]$:

[TABLE]

with the non-conservative force $F(x_{2})=f\cos(2\pi x_{2})$. The stationary probability distribution is homogeneous, $P^{\mathrm{st}}(x_{1},x_{2})=1$ and the steady state flux is $\vec{J}^{\mathrm{st}}(x_{1},x_{2})=F(x_{2}),0$.

Since the dynamics is invariant under translations along the $x_{1}$ axis, then $\vec{c}_{\mathrm{h}}(x_{1},x_{2})=c_{\mathrm{h},1}(x_{1},x_{2}),c_{\mathrm{h},2}(x_{1},x_{2})$ cannot depend on $x_{1}$. Writing Eq. (14) by components, we find that $c_{\mathrm{h},2}(x_{2})=0$ (see SI). Consequently, Eq. (14) reduces to the one-dimensional equation in the unknown $c_{\mathrm{h},1}(x_{2})$

[TABLE]

where the kernel is

[TABLE]

Since the coordinate $x_{2}$ evolves according to a simple diffusion process with periodic boundary conditions, the function $\phi(x_{2},y_{2})$ can be explicitly expressed as

[TABLE]

Expanding the solution $c_{\mathrm{h},1}(x_{2})$ in a Fourier basis and substituting into Eq. (23), the Fourier coefficients can be analytically calculated at any order (see SI).

In this case the $\mathrm{CV}^{2}$ of the hyperaccurate current is much lower than that of the entropy production far from equilibrium, i.e., when $f\gg 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 $\mathrm{CV}^{2}$ of an arbitrary current. Our theory can be extended to discrete-state or discrete-time systems and possibly employed to study non integrated currents or non stationary dynamics. We leave these investigations for future work.

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 Li et al. (2019) 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 $\mathrm{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 $\mathrm{CV}^{2}$ is sufficiently close to the bound. Reference Li et al. (2019) 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 $\mathrm{CV}^{2}$ of the hyperaccurate current in experiments. To pursue this strategy, it will be key to develop efficient numerical schemes Delves and Mohamed (1988) to solve the integral equation (14) and therefore compute the hyperaccurate current in systems more complex than the simple examples considered in 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ú, E. Fried, F. Jülicher, A. Maritan, I. Neri, L. Peliti, and É Roldán for many discussions.

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 $\sigma^{2}_{R_{\mathrm{bound}},R_{\mathrm{dev}}}=0$. In Section III we derive the Euler-Lagrange equation for the vector field $\vec{c}_{\mathrm{h}}(\vec{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 $\vec{c}_{\mathrm{h}}(\vec{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

[TABLE]

where $\vec{f}(\vec{x}(t^{\prime}),t^{\prime},\vec{y}(t^{\prime\prime}),t^{\prime\prime})$ is an arbitrary function of an Ito process at two different times (see also Pigolotti et al. (2017)).

For $t^{\prime\prime}\geq t^{\prime}$, the average in Eq. (26) always vanishes due to the non-anticipating properties of the Wiener process in the Ito calculus. The case $t^{\prime\prime}<t^{\prime}$ 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:

[TABLE]

and impose a future condition $\vec{x}(t^{\prime})=\vec{x}^{\prime}$, with $t^{\prime}>t$. Doob showed that the ensemble of trajectories $\{\vec{x}\}$ of Eq. ,(27) conditioned on the future event $\vec{x}(t^{\prime})=\vec{x}^{\prime}$ is equal to the unconditioned ensemble of trajectories $\{\vec{z}\}$ generated by the Langevin equation

[TABLE]

where $\vec{\eta}(t^{\prime\prime})$ 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

[TABLE]

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 $\vec{\eta}$ vanishes, we obtain

[TABLE]

II Proof that $\sigma^{2}_{R_{\mathrm{bound}},R_{\mathrm{dev}}}=0$

We now demonstrate that the covariance term $\sigma^{2}_{R_{\mathrm{bound}},R_{\mathrm{dev}}}=\langle R_{\mathrm{bound}}R_{\mathrm{dev}}\rangle-\langle R_{\mathrm{bound}}\rangle\langle R_{\mathrm{dev}}\rangle$ vanishes, where $R_{\mathrm{bound}}$ and $R_{\mathrm{dev}}$ are defined in Eqs. (6) and (8) of the Main Text. Since $\langle R_{\mathrm{bound}}\rangle=0$, we express the covariance as

[TABLE]

where in the last term we expressed the average over the noise $\xi(t^{\prime\prime})$ conditioned at future time using the Doob transform, see Section I. The quantity $\psi$ is defined as

[TABLE]

We now rewrite the last expression in Eq. (31) using that $\langle\vec{J}^{\mathrm{st}}\cdot\vec{c}/P^{\mathrm{st}}\rangle=\int d\vec{x}~{}\vec{J}^{\mathrm{st}}(\vec{x})\cdot\vec{c}(\vec{x})=\langle R\rangle$ (see Eq. (12) in the Main Text) and integrating by parts over $\vec{y}$ the last term. We obtain

[TABLE]

since $\vec{\nabla}_{\vec{y}}\cdot\vec{J}^{\mathrm{st}}(\vec{y})=-\partial_{t}P^{\mathrm{st}}(\vec{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(\vec{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 $\langle R^{2}\rangle/\langle R\rangle^{2}$. Our first step is to derive tractable expressions for the first two moments of the current. The average reads

[TABLE]

The second moment can be expressed as

[TABLE]

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

[TABLE]

Integrating by parts the last term over $\vec{y}$ and substituting the definition of $\psi$ we obtain

[TABLE]

We introduce a new variable

[TABLE]

With this definition, the current is equal to the entropy production if $\vec{\Gamma}(\vec{x})=\vec{J}^{\mathrm{st}}(\vec{x})$. In terms of the new variable, the average of $R$ reads

[TABLE]

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

[TABLE]

with

[TABLE]

Note that $\hat{K}^{(1)}$ is a tensor, $K^{(2)}$ is a scalar, and $\vec{K}^{(3)}$, $\vec{K}^{(4)}$ are vectors. We now compute the first variation respect to the function $\vec{\Gamma}(\vec{x})$

[TABLE]

where

[TABLE]

while

[TABLE]

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

[TABLE]

with

[TABLE]

Imposing from Eqs. (43) that the first variation vanishes for any choice of $\eta$ and using Eqs. (44), and (III) leads to the condition

[TABLE]

Substituting the expression of $\vec{G}(\vec{x},t)$ and expressing the equation in terms of $\vec{c}_{\mathrm{h}}(\vec{x})$ we finally obtain Eq. (13) in the Main Text.

IV Long-time limit

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

[TABLE]

In the limit $t\rightarrow\infty$, the first two term on the left-hand side of Eq. (IV) 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 $\vec{J}^{\mathrm{st}}(\vec{x})\langle R_{\mathrm{h}}\rangle/2$, that we rewrite as

[TABLE]

Equation (IV) then becomes

[TABLE]

At this point it is safe to take the limit $t\rightarrow\infty$. Performing the integral over time and using that at steady state the propagator depends only on the time difference $t^{\prime}-t^{\prime\prime}$ we obtain

[TABLE]

where the function $\phi(\vec{x},\vec{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

[TABLE]

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

Writing explicitly the spatial derivatives results in

[TABLE]

The system is periodic in the interval $[0,1]$. We discretize this interval by introducing a mesh $\Delta$. In this way, a function of $x$ becomes a function of a discrete index $i=\lfloor x/\Delta\rfloor$, where $\lfloor\dots\rfloor$ 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

[TABLE]

where we call $M_{ij}$ the matrix obtained by discretizing the kernel of the integral equation. Such matrix reads

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. We use the notation $(\nabla_{x}P^{\mathrm{st}})_{i}$, $(\nabla_{y}\phi)_{ij}$, and $(\nabla_{xy}\phi)_{ij}$ for the discretized derivatives, that are defined as

[TABLE]

In the $1D$ 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

[TABLE]

For the numerical integration, we approximate the initial condition $P(x;0|y;0)=\delta(x-y)$ with a Gaussian distribution with mean $y$ and variance $3\times 10^{-6}$. We use the built-in solver in Mathematica with a spatial mesh $\Delta=0.002$, a time step $\Delta 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 $\phi(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^{\mathrm{st}}(\vec{x})+\mathcal{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, $\vec{c}_{\mathrm{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

[TABLE]

where $\phi(x_{2},y_{2})$ can be computed explicitly from the diffusion equation

[TABLE]

Eq. (61) is solved by $c_{\mathrm{h},2}(x_{2})=0$, so that the hyperaccurate vector field $\vec{c}_{\mathrm{h}}(x_{2})=(c_{\mathrm{h},1}(x_{2}),0)$ is governed by the one-dimensional Eq. (60). We now expand the solution in a Fourier basis

[TABLE]

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

[TABLE]

To make progress, we use the properties of trigonometric functions

[TABLE]

Expressing the products of trigonometric functions using these relations we obtain

[TABLE]

The coefficients associated to the sine vanish, $\chi_{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

[TABLE]

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 $30$th order. Higher order coefficients were smaller than $10^{-9}$ for all values of the force in the explored range $f\in[0.1,100]$.