\section*{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~\ref{I}, we show how to evaluate averages of stochastic integrals that we often use in the following. In Section~\ref{II} we show that $\sigma^2_{R_\mathrm{bound},R_\mathrm{dev}}=0$. In Section~\ref{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~\ref{IV}.  In Section~\ref{V} we discretize the integral kernel of the Euler-Lagrange equation for solving it numerically on a one-dimensional grid. Finally, in Section~\ref{VI} we analytically compute $\vec{c}_\mathrm{h}(\vec{x})$  for the two-dimensional model presented in the Main Text.

\section{Evaluating two-time averages with the Doob transform}\label{sec:doob}
\label{I}

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

\begin{equation}\label{eq:doobtrick}
I=\left\langle \int~d t' \int d t''~\vec{f}(\vec{x}(t'),t',\vec{y}(t''),t'') \cdot\vec{\xi}(t'') \right \rangle
\end{equation}
where $\vec{f}(\vec{x}(t'),t',\vec{y}(t''),t'')$ is an arbitrary function of
an Ito process at two different times (see also \cite{pigolotti2017generic}).

For $t''\ge t'$, the average in Eq. \eqref{eq:doobtrick} always vanishes due to the non-anticipating properties of the Wiener process in the Ito calculus. The case $t''<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:
\begin{equation}
\frac{d}{d t} \vec{x} = \hat{\mu} \cdot \vec{F} + \vec{\nabla} \cdot \hat{D} + \sqrt{2} \hat{\sigma} \cdot \vec{\xi}
\label{Langevin}
\end{equation}
and impose a future condition $\vec{x}(t') = \vec{x}'$, with $t' > t$.  Doob showed that the ensemble of trajectories $\{\vec{x}\}$ of Eq. ,\eqref{Langevin} conditioned on the future event $\vec{x}(t') = \vec{x}'$ is equal to the unconditioned ensemble of trajectories $\{\vec{z}\}$ generated by the Langevin equation
\begin{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'|\vec{y}; t'')
\label{LangevinDoob}
\end{equation}
where $\vec{\eta}(t'')$ is unbiased white noise. Comparing Eq. \eqref{Langevin} and Eq. \eqref{LangevinDoob}, we obtain that the conditioned averages can be transformed into unconditioned ones by substituting the noise term
\begin{equation}
\vec{\xi}(t) = \vec{\eta}(t) + \sqrt{2} \hat{\sigma}^T \vec{\nabla}_{\vec{y}} \log P(\vec{x}; t'|\vec{y}; t).
\end{equation}
The second term in the right hand side represents an additional drift. Substituting
this expression in Eq. \eqref{eq:doobtrick} and using that the average of any function multiplied by $\vec{\eta}$ vanishes, we obtain

\begin{equation}\label{eq:doobtrick2}
I=\sqrt{2}\left\langle \int~d t' \int d t''~ \vec{f}(\vec{x}(t'),t',\vec{y}(t''),t'')\cdot  \hat{\sigma}^T \cdot\vec{\nabla}_{\vec{y}} \log P(\vec{x}; t'|\vec{y}; t'')\cdot \right \rangle
\end{equation}

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

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
\begin{eqnarray}\label{eq:covar}
\sigma^2_{R_\mathrm{bound},R_\mathrm{dev}}&=&\langle R_\mathrm{bound}~(R-R_\mathrm{bound}) \rangle=\nonumber\\
&=& - \sigma^2_{R_{\rm bound}} + 2 \frac{\langle R \rangle}{\langle \Sigma \rangle} \left\langle \frac{\vec{J}^{\mathrm{st}}}{P^{\mathrm{st}}} \cdot (\hat{\sigma}^T)^{-1} \hat{\sigma}^T \cdot \vec{c} \right\rangle + \nonumber \\
&+& 2 \frac{\langle R \rangle}{\langle \Sigma \rangle} \int_0^t d t' \int_0^{t'} d t'' \int d \vec{x}~ d \vec{y}~ \frac{\psi(\vec{y})}{P^{\mathrm{st}}(\vec{y})} \frac{\vec{J}^{\mathrm{st}}(\vec{y})}{P^{\mathrm{st}}(\vec{y})} \cdot [\vec{\nabla}_{\vec{y}} P(\vec{x} t'|\vec{y} t'') ]P^{\mathrm{st}}(\vec{y})
\end{eqnarray}
where in the last term we expressed the average over the noise $\xi(t'')$ conditioned at future time using the Doob transform, see Section \ref{sec:doob}. The quantity $\psi$ is defined as
\begin{equation}\label{eq:phi}
\psi(\vec{x})=\vec{c} (\vec{x})\cdot \vec{J}^{\mathrm{st}}(\vec{x}) + \vec{\nabla}_{\vec{x}} \cdot [\hat{D}(\vec{x}) \cdot \vec{c}(\vec{x}) P^{\mathrm{st}}(\vec{x})].
\end{equation}
We now rewrite the last expression in Eq. \eqref{eq:covar} 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 
\begin{equation}
\sigma^2_{R_\mathrm{bound},R_\mathrm{dev}} = - 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' \int_0^t d t'' \int ~d \vec{x} ~d \vec{y}~ \frac{\psi(\vec{x})}{P^{\mathrm{st}}(\vec{x})} P(\vec{x}; t'|\vec{y}; t'') \vec{\nabla}_{\vec{y}} \cdot \vec{J}^{\mathrm{st}}(\vec{y}) = 0
\end{equation}
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.
\section{Derivation of the Euler-Lagrange equations}
\label{III}
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
\begin{equation}
\langle R\rangle=t\left\langle \frac{\vec{c} \cdot \vec{J}^{\mathrm{st}}}{P^{\mathrm{st}}}\right\rangle
\end{equation}
The second moment can be expressed as
\begin{eqnarray}\label{eq:tominimize}
\langle R^2\rangle &=& 2  t\langle \vec{c} \cdot \hat{D} \cdot \vec{c} \rangle + 2 \int_0^t d t' \int_{0}^{t'} d t''\left[\left\langle \frac{\psi(\vec{x})}{P^{\mathrm{st}}(\vec{x})} \frac{\psi(\vec{y})}{P^{\mathrm{st}}(\vec{y})} \right\rangle + \sqrt{2}\left\langle \frac{\psi(\vec{x})}{P^{\mathrm{st}}(\vec{x})} \vec{c}(\vec{y}) \cdot \hat{\sigma}^{T}(\vec{y})\cdot \vec{\xi}(t'') \right\rangle\right].
\end{eqnarray}
To evaluate the last term on the right hand side of Eq. \eqref{eq:tominimize} we again apply the Doob transform, obtaining
\begin{eqnarray}
\langle R^2\rangle &=& 2 t\langle \vec{c} \cdot \hat{D} \cdot \vec{c} \rangle + \\
&+&
2\int_0^t d t' \int_{0}^{t'}d t''\intd \vec{x}~d \vec{y}~ P(\vec{x},t'|\vec{y},t'')P^\mathrm{st}(\vec{y},t'') \left[\frac{\psi(\vec{x})}{P^{\mathrm{st}}(\vec{x})}  \frac{\psi(\vec{y})}{P^{\mathrm{st}}(\vec{y})} +  2 \frac{\psi(\vec{x})}{P^{\mathrm{st}}(\vec{x})} \vec{c}(\vec{y})  \cdot \hat{D}(\vec{y})\cdot \frac{\nabla_{\vec{y}} P(\vec{x},t'|\vec{y},t'') }{P(\vec{x},t'|\vec{y},t'')}\right]\nonumber
\end{eqnarray}
Integrating by parts the last term over $\vec{y}$ and substituting the definition of $\psi$ we obtain
\begin{eqnarray}\label{eq:R23}
\langle R^2\rangle &=& 2t \langle \vec{c} \cdot \hat{D} \cdot \vec{c} \rangle +\\
&+&2 \int_0^t d t' \int_{0}^{t'} d t'' \left\langle \frac{\vec{c}(\vec{x}) \cdot \vec{J}^{\mathrm{st}}(\vec{x}) + \vec{\nabla}_{\vec{x}} \cdot [\hat{D}(\vec{x}) \cdot \vec{c}(\vec{x}) P^{\mathrm{st}}(\vec{x})]}{P^{\mathrm{st}}(\vec{x})} \frac{\vec{c}(\vec{y}) \cdot \vec{J}^{\mathrm{st}}(\vec{y}) - \vec{\nabla}_{\vec{y}} \cdot [\hat{D}(\vec{y})\cdot \vec{c}(\vec{y}) P^{\mathrm{st}}(\vec{y})]}{P^{\mathrm{st}}(\vec{y})} \right\rangle \nonumber \\
\label{varRdev}
\end{eqnarray}
We introduce a new variable
\begin{equation}
\vec{\Gamma}(\vec{x}) = P^{\mathrm{st}}(\vec{x})\hat{D}(\vec{x}) \cdot \vec{c}(\vec{x}).
\end{equation}
 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
\begin{equation}
\langle R \rangle = t \int d\vec{x} ~\vec{\Gamma}(\vec{x}) \cdot \hat{D}^{-1} \cdot \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})}
\end{equation}
and, using Eq. \eqref{eq:R23}, the second moment can be expressed as
\begin{eqnarray}
\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] 
\label{R2K}
\end{eqnarray}
with
\begin{eqnarray}
\hat{K}^{(1)} &=& \frac{2 \hat{D}^{-1}(\vec{x})}{P^{\mathrm{st}}(x)} t \delta(\vec{x}-\vec{y}) \nonumber\\
&+& 2 \int_0^t d t' \int_0^{t'} d t'' P(x; t'|\vec{y}; t'') P^{\mathrm{st}}(\vec{y}) \frac{ \hat{D}^{-1}(\vec{x}) \cdot \vec{J}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})^2} \frac{\vec{J}^{\mathrm{st}}(\vec{y}) \cdot \hat{D}^{-1}(\vec{y})}{P^{\mathrm{st}}(\vec{y})^2}  \\
K^{(2)} &=& - 2 \int_0^t d t' \int_0^{t'} d t'' P(x; t'|\vec{y}; t'') P^{\mathrm{st}}(\vec{y}) \frac{1}{P^{\mathrm{st}}(\vec{x}) P^{\mathrm{st}}(\vec{y})} \nonumber \\
\vec{K}^{(3)} &=& 2 \int_0^t d t' \int_0^{t'} d t'' P(x; t'|\vec{y}; t'') P^{\mathrm{st}}(\vec{y}) \frac{1}{P^{\mathrm{st}}(\vec{x})} \frac{\vec{J}^{\mathrm{st}}(\vec{y}) \cdot \hat{D}^{-1}(\vec{y}) }{P^{\mathrm{st}}(\vec{y})^2} \nonumber \\
\vec{K}^{(4)} &=& - 2 \int_0^t d t' \int_0^{t'} d t'' P(x; t'|\vec{y}; t'') P^{\mathrm{st}}(\vec{y}) \frac{\vec{J}^{\mathrm{st}}(\vec{x})\cdot \hat{D}^{-1}(\vec{x})}{P^{\mathrm{st}}(\vec{x})^2}  \frac{1}{P^{\mathrm{st}}(\vec{y})} \nonumber
\end{eqnarray}
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})$
\begin{equation}\label{eq:firstvariation}
\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}
\end{equation}
where
\begin{equation}\label{eq:avr}
\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}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})}
\end{equation}
while
\begin{eqnarray}
\frac{\delta \langle R^2\rangle}
{\delta \vec{\Gamma}}\delta \vec{\eta} &=&\int d\vec{x} d\vec{y} ~ \Bigg[  \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})  \nonumber\\
&-&\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}) \nonumber\\
&-&\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})\nonumber\\
&+& \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)\Bigg]
\end{eqnarray}
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
\begin{eqnarray}\label{eq:avr2}
\frac{\delta \langle R^2\rangle}
{\delta \vec{\Gamma}}\cdot\delta \vec{\eta} &=& \int d\vec{x}~ \vec{\eta}(\vec{x}) \cdot \vec{G}(\vec{x},t) \nonumber \\
\end{eqnarray}
with
\begin{equation}
\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}).
\end{equation}
Imposing from Eqs.~\eqref{eq:firstvariation} that the first variation vanishes for any choice of $\eta$ and using  Eqs. \eqref{eq:avr}, and \eqref{eq:avr2} leads to the condition
\begin{equation}
\vec{G}(\vec{x},t) = 2 t \hat{D}^{-1}(\vec{x}) \cdot \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})} \frac{\langle R^2\rangle}{\langle R \rangle}
\label{firstvariationnull}
\end{equation}

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.

\section{Long-time limit}
\label{IV}

We now derive the long-time limit of Eq.~(13) in the Main Text, that we rewrite as
\begin{eqnarray}\label{eq:ltlSI}
P^{\mathrm{st}}(\vec{x}) \hat{D}(\vec{x})\cdot\vec{c_\mathrm{h}}(\vec{x}) &+& 
 \frac{P^{\mathrm{st}}(\vec{x})}{t}\hat{D}(\vec{x})\cdot \vec{\nabla}_{\vec{x}} \left\{ \int d\vec{y} \int_0^t d t' \int_0^{t'} d t'' ~\frac{P(\vec{x}; t'|\vec{y}; t'')}{P^{\mathrm{st}}(\vec{x})} \vec{\nabla}_{\vec{y}} \cdot \left[ P^{\mathrm{st}}(\vec{y}) \hat{D}(\vec{y}) \cdot \vec{c}_\mathrm{h}(\vec{y}) \right] \right\} +\\
&+& \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{t~P^{\mathrm{st}}(\vec{x})} \int d\vec{y} \int_0^t d t' \int_0^{t'} d t'' P(\vec{x}; t'|\vec{y};t'') \vec{J}^{\mathrm{st}}(\vec{y}) \vec{c}_\mathrm{h}(\vec{y}) 
=\vec{J}^{\mathrm{st}}(\vec{x}) \frac{\langle R_\mathrm{h}^2 \rangle}{2 \langle R_\mathrm{h} \rangle}\nonumber
\end{eqnarray}
In the limit $t\rightarrow\infty$, the first two term on the left-hand side of Eq.~\eqref{eq:ltlSI} 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
\begin{equation}
\frac{\vec{J}^{\mathrm{st}}(\vec{x}) \langle R_\mathrm{h} \rangle}{2} = \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{t} \int_0^t d t' \int_0^{t'} d t'' \left\langle \frac{\vec{c}_\mathrm{h} \cdot \vec{J}^{\mathrm{st}} }{P^\mathrm{st}}\right\rangle = \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{t~P^{\mathrm{st}}(\vec{x})} \int_0^t d t' \int_0^{t'} d t'' \int d \vec{y}~ P^{\mathrm{st}}(\vec{x}) \vec{c}_\mathrm{h}(\vec{y}) \cdot \vec{J}^{\mathrm{st}}(\vec{y})
 \end{equation}
Equation~\eqref{eq:ltlSI} then becomes
\begin{eqnarray}
P^{\mathrm{st}}(\vec{x}) \hat{D}(\vec{x})\cdot\vec{c_\mathrm{h}}(\vec{x})  &+& 
\frac{P^{\mathrm{st}}(\vec{x})}{t}\hat{D}(\vec{x})\cdot \vec{\nabla}_{\vec{x}} \left\{ \int d\vec{y} \int_0^t d t' \int_0^{t'} d t'' ~\frac{P(\vec{x}; t'|\vec{y}; t'')}{P^{\mathrm{st}}(\vec{x})} \vec{\nabla}_{\vec{y}} \cdot \left[ P^{\mathrm{st}}(\vec{y}) \hat{D}(\vec{y}) \cdot \vec{c}_\mathrm{h}(\vec{y}) \right] \right\} +\\
&+& \frac{\vec{J}^{\mathrm{st}}(\vec{x})}{t P^{\mathrm{st}}(\vec{x})} \int d\vec{y} \int_0^t d t' \int_0^{t'} d t'' \left[ P(\vec{x}; t'|\vec{y}; t'') - P^{\mathrm{st}}(\vec{x}) \right] \vec{J}^{\mathrm{st}}(\vec{y}) \cdot \vec{c}_\mathrm{h}(\vec{y})  = \vec{J}^{\mathrm{st}}(\vec{x}) \frac{\sigma^2_{R_\mathrm{h}}}{2 \langle R_\mathrm{h} \rangle}.\nonumber
\end{eqnarray}
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'-t''$ we obtain
\begin{eqnarray}
P^{\mathrm{st}}(\vec{x}) \hat{D}(\vec{x})\cdot\vec{c_\mathrm{h}}(\vec{x}) &+& P^{\mathrm{st}}(\vec{x}) \left\{ \int d\vec{y} ~ \hat{D}(\vec{x})\cdot\vec{\nabla}_{\vec{x}}\left[\frac{\phi(\vec{x},\vec{y})}{P^{\mathrm{st}}(\vec{x})} \right]\vec{\nabla}_{\vec{y}} \cdot \left[ P^{\mathrm{st}}(\vec{y}) \hat{D}(\vec{y}) \cdot \vec{c}_\mathrm{h}(\vec{y}) \right] \right\}+\\
&+&\frac{\vec{J}^{\mathrm{st}}(\vec{x})}{P^{\mathrm{st}}(\vec{x})} \int d\vec{y} ~\phi(\vec{x},\vec{y}) \vec{J}^{\mathrm{st}}(\vec{y}) \cdot \vec{c}_\mathrm{h}(\vec{y})  = \vec{J}^{\mathrm{st}}(\vec{x}) \frac{\sigma^2_{R_\mathrm{h}}}{2 \langle R_\mathrm{h} \rangle}
\label{eq:ltlFIN}
\end{eqnarray}
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.~\eqref{eq:ltlFIN} by parts and rewriting the whole expression as an integral equation leads to Eqs.~(14) and (15) of the Main Text.

\section{Discretization of the integral kernel in one dimension}
\label{V}

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
\begin{eqnarray}
K(x,y) = \frac{(J^{\mathrm{st}})^2}{P^{\mathrm{st}}(x)} \phi(x,y) + P^{\mathrm{st}}(x) D(x) \delta(x-y) - D(x) P^{\mathrm{st}}(x) \nabla_{xy} \left[ \frac{\phi(x,y)}{P^{\mathrm{st}}(x)} \right] P^{\mathrm{st}}(y) D(y)
\label{LTL}
\end{eqnarray}
where we used that the the stationary flux is constant in one dimension. 

Writing explicitly the spatial derivatives results in
\begin{eqnarray}
K(x,y) = \frac{(J^{\mathrm{st}})^2}{P^{\mathrm{st}}(x)} \phi(x,y) + P^{\mathrm{st}}(x) D(x) \delta(x-y) + \left[ - \nabla_{xy} \phi(x,y) + \frac{\nabla_x P^{\mathrm{st}}(x)}{P^{\mathrm{st}}(x)} \nabla_y \phi(x,y) \right] D(x) P^{\mathrm{st}}(y) D(y)
\label{LTL2}
\end{eqnarray}

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
\begin{equation}
\sum_j M_{ij}(c_\mathrm{h})_j=J^\mathrm{st}
\end{equation}
where we call $M_{ij}$ the matrix obtained by discretizing the kernel of the integral equation. Such matrix reads
\begin{eqnarray}
M_{ij}  = \frac{(J^{\mathrm{st}})^2}{P^{\mathrm{st}}_i} (\phi)_{ij} + P^{\mathrm{st}}_i D_i \delta_{ij} + \left( - (\nabla_{xy} \phi)_{ij} + \frac{(\nabla_x P^{\mathrm{st}})_i}{P^{\mathrm{st}}_i} (\nabla_y \phi)_{ij} \right) D_i P^{\mathrm{st}}_j D_j
\label{LTL2b}
\end{eqnarray}
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
\begin{eqnarray}
(\nabla_x P^{\mathrm{st}})_{i} &=& \frac{-P^{\mathrm{st}}_{i+2} + 8 P^{\mathrm{st}}_{i+1} - 8 P^{\mathrm{st}}_{i-1} + P^{\mathrm{st}}_{i-2}}{12 \Delta} + \mathcal{O}(\Delta^4) \nonumber \\
(\nabla_y \phi)_{ij} &=& \frac{-\phi_{i,j+2} + 8 \phi_{i,j+1} - 8 \phi_{i,j-1} + \phi_{i,j-2}}{12 \Delta} + \mathcal{O}(\Delta^4) \nonumber \\
(\nabla_{xy} \phi)_{ij} &=& \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}}{12 \Delta} + \mathcal{O}(\Delta^4)
\end{eqnarray}
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 
\begin{equation}
\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).
\end{equation}
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})$.

\section{Hyperaccurate current for the two dimensional model}
\label{VI}

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
\begin{gather}
\label{tosolve}
c_{\mathrm{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_{\mathrm{h},1}(y_2) = f \cos(2 \pi x_2) \\
\label{2component}
c_{\mathrm{h},2}(x_2) = - \int d y_2 ~\partial_{x_2,y_2} \phi(x_2, y_2) ~c_{\mathrm{h},2}(y_2)
\end{gather}
where $\phi(x_2,y_2)$ can be computed explicitly from the diffusion equation 
\begin{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))
\end{equation}
Eq.~\eqref{2component} 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.~\eqref{tosolve}.
We now expand the solution in a Fourier basis
\begin{equation}
c_{\mathrm{h},1}(y) = \kappa_{\mathrm{h},0} + \chi_{\mathrm{h},0} + \sum_{n=1}^{+\infty} \kappa_{\mathrm{h},n} \cos(2 \pi n y) + \sum_{n=1}^{+\infty} \chi_{\mathrm{h},n} \sin(2 \pi n y)
\label{solution}
\end{equation}
Inserting Eq.~\eqref{solution} in Eq.~\eqref{tosolve}, we get
\begin{gather}
\kappa_{\mathrm{h},0} + \chi_{\mathrm{h},0} + \sum_{n=1}^{+\infty} \kappa_{\mathrm{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_{\mathrm{h},n+1}+\kappa_{\mathrm{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_{\mathrm{h},2}}{4} + \frac{\kappa_{\mathrm{h},0}}{2} \right) + \nonumber \\
+ \sum_{n=1}^{+\infty} \chi_{\mathrm{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_{\mathrm{h},n+1}+\chi_{\mathrm{h},n-1}}{4} \sin(2 \pi n y_2) + \frac{f^2}{2 \pi^2} \frac{\chi_{\mathrm{h},2}}{4} \cos(2 \pi y_2) \sin(2 \pi y_2) = f \cos(2 \pi y_2)
\end{gather}
To make progress, we use the properties  of trigonometric functions
\begin{eqnarray}
\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} \\
\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} \bigg(\cos(2 \pi (n+1) y_2) + \cos(2 \pi (n-1) y_2) \bigg) \\
\sin(2 \pi n y_2)\cos(2 \pi y_2) &=& \frac{1}{2} \bigg(\sin(2 \pi (n+1) y_2) + \sin(2 \pi (n-1) y_2) \bigg)
\end{eqnarray}
Expressing the products of trigonometric functions using these relations we obtain
\begin{gather}
\kappa_{\mathrm{h},0} + \chi_{\mathrm{h},0} + \sum_{n=1}^{+\infty} \kappa_{\mathrm{h},n} \cos(2 \pi n y_2) + f^2 \sum_{n=2}^{+\infty} \frac{\kappa_{\mathrm{h},n+1}+\kappa_{\mathrm{h},n-1}}{16 \pi^2 n^2} \bigg(\cos(2 \pi (n+1) y_2) + \cos(2 \pi (n-1) y_2) \bigg) + \nonumber \\ + \frac{f^2}{4 \pi^2} \bigg(\cos(4 \pi y_2) + 1 \bigg) \left( \frac{\kappa_{\mathrm{h},2}}{4} + \frac{\kappa_{\mathrm{h},0}}{2} \right) + \nonumber \\
+ \sum_{n=1}^{+\infty} \chi_{\mathrm{h},n} \sin(2 \pi n y_2) + f^2 \sum_{n=2}^{+\infty} \frac{\chi_{\mathrm{h},n+1}+\chi_{\mathrm{h},n-1}}{16 \pi^2 n^2} \bigg(\sin(2 \pi (n+1) y_2) + \sin(2 \pi (n-1) y_2) \bigg) + \nonumber \\
+ \frac{f^2}{4 \pi^2} \frac{\chi_{\mathrm{h},2}}{4} \bigg(\sin(4 \pi y_2) + 1 \bigg) = f \cos(2 \pi y_2)
\end{gather}
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
\begin{eqnarray}
\kappa_{\mathrm{h},0} &+& \frac{f^2}{4 \pi^2} \left( \frac{\kappa_{\mathrm{h},2}}{4} + \frac{\kappa_{\mathrm{h},0}}{2} \right) = 0 \nonumber \\
\kappa_{\mathrm{h},1} &+& \frac{f^2}{8 \pi^2} \frac{\kappa_{\mathrm{h},3} + \kappa_{\mathrm{h},1}}{8} = f \nonumber \\
\kappa_{\mathrm{h},2} &+& \frac{f^2}{18 \pi^2} \frac{\kappa_{\mathrm{h},4}+\kappa_{\mathrm{h},2}}{8} + \frac{f^2}{4 \pi^2} \left( \frac{\kappa_{\mathrm{h},2}}{4} + \frac{\kappa_{\mathrm{h},0}}{2} \right) = 0 \nonumber \\
\kappa_{\mathrm{h},n} &+& f^2 \frac{\kappa_{\mathrm{h},n} + \kappa_{\mathrm{h},n-2}}{16 \pi^2 (n-1)^2} + f^2 \frac{\kappa_{\mathrm{h},n+2} + \kappa_{\mathrm{h},n}}{16 \pi^2 (n+1)^2} = 0 \;\;\;\;\ \mbox{for } n > 2
\label{setofeq}
\end{eqnarray}
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]$.

\end{widetext}

\end{document}
