Hyperaccurate bounds in discrete-state Markovian systems

Generalized empirical currents represent a vast class of thermodynamic observables of mesoscopic systems. Their fluctuations satisfy the thermodynamic uncertainty relations (TURs), as they can be bounded by the average entropy production. Here, we derive a general closed expression for the hyperaccurate current in discrete-state Markovian systems, i.e., the one with the least fluctuations, for both discrete- and continuous-time evolution. We show that its associated hyperaccurate bound is generally much tighter than the one given by the TURs, and might be crucial to providing a reliable estimation of the average entropy production. We also show that one-loop systems (rings) exhibit a hyperaccurate current only for finite times, highlighting the importance of short-time observations. Additionally, we derive two novel bounds for the efficiency of work-to-work converters, solely as a function of either the input or the output power. Finally, our theoretical results are employed to analyze a 6-state model network for kinesin, and a chemical system in a thermal gradient exhibiting a dissipation-driven selection of states.

Generalized empirical currents represent a vast class of thermodynamic observables of mesoscopic systems. Their fluctuations satisfy the thermodynamic uncertainty relations (TURs), as they can be bounded by the average entropy production. Here, we derive a general closed expression for the hyperaccurate current in discrete-state Markovian systems, i.e., the one with the least fluctuations, for both discrete-and continuous-time evolution. We show that its associated hyperaccurate bound is generally much tighter than the one given by the TURs, and might be crucial to providing a reliable estimation of the average entropy production. We also show that one-loop systems (rings) exhibit a hyperaccurate current only for finite times, highlighting the importance of short-time observations. Additionally, we derive two novel bounds for the efficiency of work-to-work converters, solely as a function of either the input or the output power. Finally, our theoretical results are employed to analyze a 6-state model network for kinesin, and a chemical system in a thermal gradient exhibiting a dissipation-driven selection of states.
The breaking of detailed balance, a positive total entropy production rate, the presence of steady probability currents, and a limited efficiency (in the case of thermal engines) are only a few possible fingerprints of a nonequilibrium picture. Some of these features are also intimately connected through the celebrated fluctuation theorems [2,21] and satisfy universal bounds known as thermodynamic uncertainty relations (TURs), generally dictating that the dissipation constraints current fluctuations out of equilibrium. TURs have attracted increasing attention in recent years, hinting at the fascinating perspective of estimating the entropy production by measuring stochastic currents [22][23][24].
In its original formulation, the TUR relates fluctuations of any stochastic currents in steady state arbitrarily far from equilibrium to the average total entropy production rate, Σ [25]: where σ 2 J and J are the variance and mean of the * daniel.busiello@epfl.ch current J over the ensemble of stochastic trajectories, respectively, being such left side known as coefficient of variation squared (CV 2 ) of J. The TUR expressed in Eq. (1) has been proven to hold both for Markovian discrete-state systems [26] and for continuous-state systems following a Langevin dynamics [27][28][29][30]. Subsequently, TURs have been extended to several cases, such as periodically-driven systems [31] and discretetime processes [32,33], in turn generating a wealth of novel bounds in stochastic thermodynamics [34,35] and highlighting their connection with fluctuation theorems [36,37]. Additionally, several TURs have been recently unified under a geometric interpretation [38]. As stated before, besides the richness of its physical content, i.e. the minimum amount of dissipation required to have a current of a desired precision, TURs also play a leading role in estimating the average entropy production, Σ , by inverting Eq. (1). Some works in this direction exploited the saturation of the bound in short-time experiments [24], even if the bottleneck of this inference problem relies on the ability to identify a current approaching the bound, so to provide a reliable estimate of Σ . In [39], a closed expression for the hyperaccurate current, i.e. the one minimizing the CV 2 , is derived for a set of overdamped Langevin equations. This is clearly the best observable to bound the average entropy production rate using Eq. (1).
Here, we generalize the concept of hyperaccurate current to Markovian discrete-state systems. We derive a general closed expression for the hyperaccurate current in the case of both discrete (Markov chains) and continuous (Master Equation) time evolution. For systems with only one loop in the transition network (rings), we show that all currents have the same CV 2 in the longtime limit, while finite-time hyperaccurate currents can be defined. Conversely, in the presence of more than one loop, we derive the hyperaccurate current and its associated bound, both for finite times and in the long-time regime. The knowledge of the hyperaccurate current can also provide two novel bounds for the efficiency of general work-to-work converters, respectively as a function solely of the output or input work. We then illustrate our theory for two paradigmatic master equation systems, a six-state model for kinesin moving along a microtubule [40,41], and a chemical system in a thermal gradient exhibiting a dissipation-driven selection of states [10].

II. GENERALIZED EMPIRICAL CURRENTS
Consider a stochastic trajectory performed by a discrete-state system in the time interval t ∈ [0, t f ]. This is characterized by the set of visited states, {x i } i=0,...,N . The generalized empirical currents associated with this trajectory are defined as [26]: where n ml is the number of jumps from the state l to m up to time t f , and d ml the element (ml) of an antisymmetric matrixd. The specific form ofd determines the current. Clearly, J is a trajectory-dependent quantity, since n ml depends on the set {x i } i=0,...,N as follows: where δ i,j attempts to the Kronecker delta.
To evaluate the CV 2 for discrete-state systems, we compute average and variance of a generalized empirical current over all stochastic trajectories with the same duration t f . From Eq. (2) and by employing the antisymmetric property ofd together with the fact that it does not depend on the trajectory, the average current is given by where the sum now runs over all indices m < l, and j ml is the average current from the state l to m: Analogously, the variance of J reads as follows: where C mlm l = n ml n m l − n ml n m l . Using again the anti-symmetry of d ml , we can restrict the summation in Eq. (6) over all indices m < l and m < l , obtaining: with M mlm l = C mlm l + C lml m − C mll m − C lmm l . Later on we will determine the explicit form of j ml and σ 2 J for Markov chains and master equation systems.

III. HYPERACCURATE CURRENTS AND BOUND
The hyperaccurate current is determined by the matrix d that minimizes the CV 2 , namelyd (h) . Hence, for each element d ij we have to solve the following equation: where J (h) and σ 2 J (h) correspond to the mean and variance of the hyperaccurate current, respectively, and {m, l} is a short notation to denote that the sum is constrained to m < l and m < l . Analogously, σ 2 J (h) is the variance of the hyperaccurate current. We can exploit the fact that CV 2 does not change when multiplyingd (h) by a constant. The solution of Eq. (8) is therefore defined up to an arbitrary factor. From now on, we shall fix this constant by setting σ 2 J (h) / J (h) = 1. This procedure is analogous to the one employed in [39]. Moreover, since both M mlij and j ij diverge linearly with t f in the longtime limit, it is convenient to introduce the scaled quan-titiesM mlm l = M mlm l /t f andj ml = j ml /t f that stay finite when t f → ∞. With these choices, from Eq. (8), the hyperaccurate coefficients d Moreover, we arrive at the general expression for the hyperaccurate bound, that is the minimum possible value of the CV 2 of any generalized empirical current: It is worth mentioning that Eqs. (9) and (10) hold for any discrete-state system, whether it is described by a Markov chain or evolves according to a master equation, both at and out of the steady-state. In the next sections, we shall derive the statistics of the currents, namely mean and variance, for Markov chains and master equation systems, restricting ourselves to the relatively simple, yet quite general, case of stationary processes for simplicity.

A. Statistics of currents for Markov chains
Markov chains are characterized by a discrete-time evolution. Indeed, a transition between discrete states can only happen at a definite time interval, ∆t. As a consequence, a trajectory of length t f will be necessarily constituted by N = t f /∆t transitions. Let p i;t the probability to be in the state i at time t, the dynamics of a Markov chain is given by being indeed fully specified by the transition matrix A ij = p i;t+∆t|pj ;t . Hence, given a stochastic trajectory {x i } i=0,...,N := x, taking place in the time interval t ∈ [0, t f ], its path probability reads: This probability is properly normalized, i.e.
For stationary processes, the initial condition p x0;0 is equal to the steady-state probability distribution p st x0 . Therefore, the average number of jumps from the state l to the state m is given by: where we used the following properties of the transition matrix A xixi−1 : A xixi−1 = 1, and Following a similar procedure, it is possible to compute the second moment, n ml n m l , which is given by: where summations over k and k in the equation above includes three kinds of terms: a first one including only trajectories in which k > k, a second one taking contributions from trajectories in which k < k, and a third one accounting for the cases k = k. By evaluating each term separately, we arrive at the following expression: Hyperaccurate coefficients, and thus the hyperaccurate bound, can be readily obtained substituting Eqs. (13) and (14) into Eq. (9) for stationary processes described by a Markov chain.

B. Statistic of currents for master equation systems
To adapt the formalism developed in the previous subsection to master equation systems, it is sufficient to modify the transition matrix as follows: where W ij is now the transition rate from the state j to the state i and corresponds to the (ij)-th element of the transition rate matrixŴ . As in the previous section, we consider time-independent transition rates to ensure that the system will eventually reach a unique stationary state. With this form of A ij , the moments of the generalized currents can be obtained by performing the continuous-time limit, that is ∆t → 0. Hence, we obtain: It is worth pointing out that the number of jumps occurring in a single trajectory is not fixed a-priori by its duration for master equation systems. As outlined above, Eqs. (16) and (17) determine the hyperaccurate current and its associated bound. However, it is instructive to derive the equation for hyperaccurate coefficients explicitly in the long-time limit, i.e. t f → +∞. Noting that: and following the same procedure outlined in [39], we arrive at the following expression forC mlm l : As expected, there are no divergences in the long-time limit and the propagators only depend on time differences since we are considering stationary processes. Since the expression forC mlm l enters into the definition ofM mlij , we can write the equation for the hyperaccurate coefficients as follows: where J st ij = W ij p st j − W ji p st i are the steady-state probability currents obtained from the master equation.
A useful way to handle Eq. (19) is to expand all probabilities in terms of the eigenvalues and eigenvectors of transition matrixŴ , i.e.
where M is the total number of accessible states in the system, v (i) l is the l-th component of the i-th eigenvector, and λ i its associated eigenvalue. Note that the eigenvalues are enumerated in descending order so that λ 1 = 0 > λ 2 > · · · > λ M . Here, the initial conditions, i.e. the fact the system is in the state m at time 0, are encoded in the coefficients a (m) i satisfying the following equations: Hence,C mlij takes the following form: As a starting point, consider a discrete-state system constituted by one single loop (a ring) in the transition network. At stationarity, all edges not belonging to the loop satisfy the detailed balance, and hence such a system can only support a unique nonzero steady current, J st . In the long-time limit, all possible currents of the system have to be proportional to J st , so that their CV 2 is equal to the one of J st , CV 2 st : where α is the proportionality factor. However, due to the transient regime of the propagators, it is possible to define a finite-time hyperaccurate current for any oneloop system, even in the presence of stationary processes. As illustrated in Fig. 1 for a simple 4-state ring, the hyperaccurate bound changes over time (dot-dashed vertical lines with increasing opacity). As time increases, the probability distribution function of the CV 2 becomes narrower and narrower, eventually becoming a Dirac δ centered at B h (red star) when t → ∞, since all currents have the same CV 2 in the long-time limit. In the presented example, we have 2/ Σ < B h (t) ≤ B h , where the first inequality comes from the definition of the hyperaccurate bound. Hence, by inverting the finite-time hyperaccurate bound, one can get an improved estimation of the average entropy production at stationarity, namely Σ hyp (t). Indeed, we have: In the inset of Fig. 1, we illustrate that the short-time estimation of the entropy production gets closer to the actual value of Σ than in the long-time limit. This finding has been found to be robust across all one-loop systems we numerically studied, suggesting that finite-time observations might be beneficial to estimate the dissipation of discrete-state Markovian systems, in line with a recent result obtained for overdamped systems [24]. A proof of the generality of this property, along with the subsequent design of a feasible experimental procedure to take advantage of it, might be a fascinating topic for future works.

IV. HYPERACCURATE EFFICIENCY BOUNDS
Hyperaccurate currents set the tightest possible bound to the entropy production rate of stochastic systems. As a consequence, they also provide us with general bounds on the efficiency of work-to-work converters that do not require a specific knowledge of system features, being instead directly linked to the hyperaccurate bound B h .
Work-to-work converters are a broad class of molecular engines operating at the nanoscale at a constant temperature (isothermal). They convert a given form of input work, e.g., chemical, into a different form of output work, e.g. mechanical, with a limited efficiency [42]. They substantially differ from heat engines, as working substances do not undergo cyclic transformation between two different temperatures. Cellular transporters [40,41,43], catalytic enzymes [18,44], Hsp70 chaperones [45] are a few prominent examples of these machines in biology.
The most relevant quantities in this scenario are the input, W in , and the output power, W out . Here, we derive two general upper bounds for the efficiency, each one depending on either W in or W out . We start with the expression the steady-state average entropy production: which can be rewritten in the usual bilinear form Σ = e J e F e , where F e and J e are thermodynamic forces and fluxes, respectively, and the sum runs over all fundamental cycles [46]. Clearly, J e is in general a linear combination of some microscopic stationary fluxes, J st mn . Analogously, F e will correspond to a combination of some microscopic forces, F mn = log(W mn /W nm ).
To define an operating work-to-work converter, we consider the presence of a load and a drive force, respectively F l and F d , with their corresponding fluxes, J d and J l . Hence, we have: where the right-hand side of this equation corresponds to the first law of thermodynamics with no variations of internal energy. To operate as an engine, one necessarily requires that the dissipation provided by the driving force, W in ≥ 0, generates a work that counteracts the external load, i.e. W out ≤ 0. Hence, the efficiency can be defined as follows: By combining Eqs. (25) and (26) together with the fact that B h ≥ 2/ Σ by construction, we obtain: Analogously, a second bound involving W in is derived: Notice that both these bounds do not require specific knowledge of system features, and they have to hold simultaneously at steady state. We stress the fact that η out b requires, in principle, the measurement of the output power, whereas η in b is solely based on the a-priori knowledge of the input work, e.g., the available chemical energy from ATP [43,45]. Moreover, it is possible to show that η out b ≤ η in b , meaning that the knowledge of the output power provides a tighter bound to the efficiency. This finding also agrees with the naive expectation that W out is more informative than W in to predict the efficiency of a work-to-work converter.

A. Hyperaccurate current and efficiency in a model network for kinesin
Kinesin is a molecular motor playing a fundamental role in biological processes, including mitosis, meiosis, and the transport of cellular cargo [43,47]. It consists of two amino acid chains forming a coiled coil with two motor heads on one end that are able to bind to microtubules. The other end of the dimer binds to cellular organelles. Kinesin performs processive walks on microtubules by subsequent binding and unbinding of the two heads. The hydrolysis of one ATP (adenosine triphosphate) into an ADP (adenosine diphosphate) and an inorganic phosphate (P) in the catalytic site placed in the motor head drives conformational changes that make the walk possible [40,41]. It constitutes a remarkable example of a work-to-work converter since it transduces chemical energy into mechanical work.
Here, we calculate both the hyperaccurate current and the efficiency bounds for kinesin, by applying the developed framework to the six-state transition network introduced in [40,41]. Let us start with a brief introduction to the model. Each state is determined by the chemical composition of the two motor heads, e.g., ATP, ADP, or empty. Since we aim at describing processive motion, we ignore states in which both heads have the same composition [40]. The network of all possible transitions is sketched in Fig. 2. The system moves from the state 1 to 2 and from 4 to 5 via ATP binding; ADP binding drives the transition from the state 6 to 5 and from 3 to 2; the transition from the state 6 to 1 and from 3 to 4 are associated with ATP hydrolysis. Moreover, the dashed arrows identify the forward (from 2 to 5) and backward (from 5 to 2) mechanical steps.
There are six possible cycles in this network. F + = 1 → 2 → 5 → 6 → 1, indicated in Fig. 2, encompasses the ATP hydrolysis and the subsequent forward step. Conversely, B + = 2 → 3 → 4 → 5 → 2 (see Fig. 2) converts the energy from ATP hydrolysis into a backward step. Additionally, the system can also hydrolyze two ATP molecules, while performing no steps, following the purely dissipative cycle D + = 1 → 2 → · · · → 6 → 1, not reported in figure. Clearly, also the opposite cycles involving ADP synthesis can be performed, namely F − , B − , and D − . The net processive walk is given by a competition between forward and backward cycles.
The dynamics of this system is controlled by two in-  [40,41]. F + and B + denote the forward and backward cycles, respectively. ATP states are indicated in red, while ADP states are in blue. The dashed arrows indicates the mechanical step, taken forward from 2 to 5, backward viceversa. dependent parameters, the dimensionless load force and f = LF /k B T (F and L being the load force and the step size, respectively), and the chemical energy available from ATP hydrolysis, in dilute conditions. The dynamics of kinesin can be by a master equation in which the transition rate W ij from the state j to i is given by: where κ ij is a constant, I ij ([X]) = [X] only if the reaction involves binding of X, otherwise it is equal to 1, and Φ(f ) takes the following form: with θ and χ ij additional constant factors. This choice of the transition rates agrees with the experimental observations, and also satisfies the energetic balance for each cycles [40]. This condition states that detailed balance holds if no energy is available from ATP, otherwise chemical energy is converted into mechanical motion. For example, by inspecting the cycle F + , we have the following energetic balance: We notice that K eq in Eq. (29) can be written as: K eq = κ 52 κ 21 κ 65 κ 16 κ 25 κ 12 κ 56 κ 61 = κ 25 κ 54 κ 32 κ 43 κ 52 κ 45 κ 23 κ 34 .
From the transition matrix, hyperaccurate coefficients and bound can be readily obtained by employing the framework outlined above. In Fig. 3, we report B h (solid red) together with the bound provided by the TUR (dotdashed black) and an ensemble of CV 2 of random currents (black dots), as a function of the dimensionless load force f . By construction, B h provides the tighest possible bound to the CV 2 and is markedly tighter than 2/ Σ . Indeed, we also report in the inset the ratio B h Σ /2 for different values of ∆µ, which quantifies the difference between B h and the TUR bound. As the system approaches equilibrium, this ratio decreases, as expected [39,48]. Moreover, its behavior as a function of f is non-monotonous, exhibiting also the presence of a peak corresponding to the value k B T f = ∆µ, where mechanical cycles become futile and the motor only dissipates energy (see Eq. (32)).
To quantify the performance of the hyperaccurate efficiency bounds, we explicitly write down the steady-state entropy production. From Eq. (34), Σ reads where J st B and J st F are the steady-state probability fluxes associated with the cycles B + and F + , respectively. Notice that J st F + J st B = J st 16 + J st 43 , which is the total thermodynamic flux associated with ATP hydrolysis, ∆µ. Anal- For a given force f , large values of ∆µ allow the kinesin to work as a motor, converting chemical energy into mechanical motion. However, when ∆µ is small, the available energy is not sufficient to displace the kinesin, hence it effectively uses mechanical energy to produce ATP or ADP (depending on the sign of f and ∆µ) [43]. For simplicity, we perform the numerical analysis for f > 0 and ∆µ > 0, although it can be straightforwardly extended to all other cases. From Eq. (26), when the kinesin converts chemical into mechanical energy, the efficiency reads where the numerator is the output work W out = −(J st F − J st B )k B T f , since kinesin operates against the external load force. Fig. 4 shows the efficiency and its associated bounds for two representative values of ∆µ. Results for other values of ∆µ (not shown) exhibit similar features. For ∆µ = 5.52 (Fig. 4a), η out b provides a very tight bound for small values of f , while both η out b and η in b converges to the actual efficiency when η approaches its maximum. We also report a change of behavior of the kinesin before and after the maximum efficiency, highlighting the change of regime from molecular transporter to an ATP synthesizer, respectively. In these conditions, high efficiencies are also associated with small fluxes since for small ∆µ the system is close to equilibrium. Conversely, when ∆µ = 19.34 (Fig. 4b), the system is far from equilibrium and exhibits large probability fluxes. Although hyperaccurate efficiency bounds are less tight than for ∆µ = 5.52, η out b still provides a tighter bound for the efficiency than the one derived using TUR [49]. We also notice that the system stops operating as a work-to-work converter when the flux in the backward cycle is equal to the one in the forward cycle, i.e., J st B = J st F .

B. Hyperaccurate currents and dissipation-driven selection of states
As a second application, we study a three-state chemical system diffusing in a temperature gradient. This model has been introduced in Ref. [10] as a paradigmatic example of a selection of chemical states driven by internal dissipation processes. Later on, a possible solution to the furanose conundrum has been proposed starting from an analogous modelization [50]. These studies have been stimulated by, and in turn fueled, the idea that life might have been an inevitable consequence of nonequilibrium thermodynamics [10].
The system consists of three chemical states, A, B and C, living in two different boxes at two different temperatures, T 1 and T 2 with T 1 > T 2 . Moreover, each chemical species can diffusively move between the boxes, leading to a 6-state model, as depicted in Fig. 5. All possible internal transitions among states are: where X i indicates the species X in the box i.
To determine the transition matrix governing the system dynamics, we write the chemical rates in the standard Kramers' form [10]: where k XiYi is the chemical rates associated with the reaction from X i to Y i , and ∆E = E A − E B = E A − E C for simplicity. We introduce a kinetic asymmetry by setting two different energetic barriers in going from A i to B i , ∆ B , and from A i to C i , ∆ C , so that: with ∆ = ∆ B − ∆ C > 0, which means that the state C is kinetically favorable with respect to B [10].
Three-state chemical system in a temperature gradient, T1 > T2. Each state can diffuse between the boxes with the same diffusion rate, d, and A can convert into B or C in both boxes. The vertical position of a state is proportional to its energy, i.e., EA > EB = EC . Red arrows indicate the net steady-state probability flux flowing through C states only. This is the fastest net flux in the system. Blue arrows are associated to the slow net flux that only visit B states.
We also assume that all species move between boxes with the same (symmetric) diffusive rate, , which can be interpreted as the stabilization energy of C with respect to B. When T 1 = T 2 = T , the system is at thermodynamic equilibrium and eventually reaches a Boltzmann distribution in which the states B and C are equally populated since they have the same energy. Analogously, when d = 0, each box will relax to its own Boltzmann distribution with temperature T i , and the total population of B will be identical to the one of C. However, in the presence of diffusion and a temperature gradient, the system dissipates thermal energy performing diffusive cycles between boxes. In particular, two stationary fluxes emerge: one only flows through B states (blue arrows in Fig. 5) and exhibits slow dissipation, while the other only flows through C states and dissipates faster (red arrows in Fig. 5). This symmetry breaking is associated with the kinetic asymmetry in the energetic barriers and will result in a steady-state population [C] st higher than [B] st , i.e., R CB > 1.
It is possible to show that, in the limit of fast diffusion d → +∞, and small gradient T 1 T 2 , we have: where P eq m (A) is the equilibrium probability distribution of A at the average temperature T m = (T 1 + T 2 )/2 and Σ is the entropy production. Moreover, a positive correlation between R CB and Σ has been shown to hold even beyond this limit [10].
Providing a lower bound for the average entropy production, Σ , we obtain a lower bound for R CB , namely CB , provides an accurate prediction for the selection parameter close to the equilibrium, while deviations arise as the temperature gradient increases.

VI. CONCLUSIONS
Thermodynamic uncertainty relations (TURs) set universal bounds for the precision of a stochastic current, quantified as the ratio between its variance and squared mean (CV 2 ), in terms of the dissipation of the system. Thus, by inverting this inequality, it is possible to provide a lower bound to the entropy production and estimate the distance from thermodynamic equilibrium. The main advantage of this indirect approach is that it does not require the large sample sizes and observational times that are commonly required to provide a direct estimation of the entropy production. In [39], it has been pointed out how the knowledge of the hyperaccurate current, i.e., the one with the minimum CV 2 , might greatly improve our predicting power. Here, we introduced and derived the hyperaccurate current for generic discretestate Markovian systems, both evolving in discrete and continuous time. Our analytical closed formula has been tested against two models of chemical systems. As a side result, we also provided two hyperaccurate bounds for the efficiency of work-to-work converters, as a function of either the input or the output power. Possible future extensions might include the study of nonintegrated currents and nonstationary dynamics, and a connection among TURs, hyperaccuracy, and information theory, whose role is becoming dominant in understanding stochastic systems [51,52].
Additionally, we employed our framework to compute the finite-time hyperaccurate bounds for one-loop discrete-state systems (rings). Our results suggest that short-time experiments might be much more informative than the long-time limit to estimate the average dissipation. A formal proof and generalization of this statement, along with a feasible experimental approach for discretestate Markovian dynamics, is left as an intriguing perspective for the future.