Non-monotonic skewness of currents in non-equilibrium steady states

Measurements of any property of a microscopic system are bound to show significant deviations from the average, due to thermal fluctuations. For time-integrated currents such as heat, work or entropy production in a steady state, it is in fact known that there will be long stretches of fluctuations both above as well as below the average, occurring equally likely at large times. In this paper we show that for any finite-time measurement in a non-equilibrium steady state - rather counter-intuitively - fluctuations below the average are more probable. This discrepancy is higher when the system is further away from equilibrium. For overdamped diffusive processes, there is even an optimal time when time-integrated current fluctuations mostly lie below the average. We demonstrate that these effects result from the non-monotonic skewness of current fluctuations and provide evidence that they are easily observable in experiments. We also discuss their extensions to discrete space Markov jump processes and implications to biological and synthetic microscopic engines.

Terminating the expansion of φ to the second order in λ, leads to the thermodynamic uncertainty relations [12][13][14] which are trade-off relations connecting the precision of arbitrary current measurements to the entropy production rate. In the t → ∞ limit, the LHS of Eq. (1) converges to a time independent function [9] referred to as the large deviation rate function [15][16][17][18][19], knowing which helps fully characterize the fluctuations in the long-time limit. In general, such long-time results can be obtained within the mathematical framework of large deviation theory and many such results have been obtained for the statistics of the fluctuations of entropy production [20][21][22][23][24][25], efficiency distributions [6][7][8], first passage problems [26][27][28] and current fluctuations in general [29,30]. An interesting addition to this class of results was obtained in Ref. [31], where it was shown that the fraction of time that a current spends above its average value follows the arcsine law in the long time limit [32]. As a consequence, stochastic currents with long streaks above or below their average value are much more and equally likely than those that spend similar fractions of time above and below their average.
The other extreme of very short-time fluctuations of currents, is also surprisingly non-trivial. It has been shown that Eq. (1) saturates for J = ∆S tot in the limit t → 0 for overdamped diffusive processes [33][34][35][36]. As a consequence, σ can be exactly inferred for such systems by studying the mean and variance of current fluctuations at short times [33][34][35][36], even for non-stationary systems [37]. The saturation of the bound also implies that the fluctuations of ∆S tot are Gaussian in these systems in the short-time limit even when arbitrarily far from equilibrium [33]. In fact, as was recently pointed out, the Gaussianity in the t → 0 limit holds for any arbitrary current in overdamped diffusive processes [37]. As we show below, this short-time behaviour combined with the large-deviation results mentioned above hold important clues for interesting finite-time fluctuation properties.
Finite-time fluctuations are clearly of interest since this is most often what is observed in experiments. However, when neither the t → ∞ nor the t → 0 limit can be taken, generic features of such fluctuations are harder to identify because of the prevalence of transient effects and time-correlations. If the t → ∞ limit can be thought in terms of applying a thermodynamic limit [38], then finitetime fluctuations include effects which vanish in the thermodynamic limit, making them harder to access. Some important general results in this category include the integral fluctuation theorem for stochastic entropy production [39], statistical properties of entropy production derivable from the fluctuation theorem [40], the finitetime versions of the thermodynamic uncertainty relations [11,[41][42][43], universal results known for the statistics of infima, stopping times, and first-passage probabilities of entropy production [44], a generic equation describing the time evolution of the stochastic entropy production [45], statistics of the time of the maximum of a one dimensional stationary process [46], and bounds on first passage times of current fluctuations [47].
In this paper, we unravel a previously unnoticed property of current fluctuations at finite and short-times. We demonstrate that the skewness of current fluctuations is positive, and non-monotonic in time and argue that this behaviour is generic for non-equilibrium steady states generated by any overdamped diffusive process. As a consequence, we find that at all finite times, current fluctuations below the long-time average are more probable than those above. For a single realization of the process, interestingly, this implies that a below average outcomes will be typical. Moreover, due to the non-monotonicity, there is an optimal time when this discrepancy is the highest. We show that these features of current fluctuations are easily visible in non-trivial models studied numerically and experimentally. We also discuss their extensions to discrete space Markov jump processes and implications to biological and synthetic microscopic engines. In all cases, in the limit of large t, we recover results consistent with Ref. [31].
The central results we present in this manuscript apply to non-equilibrium systems in a stationary state. We first consider generic overdamped diffusive processes of the form,ẋ where A(x) is the drift vector, and B(x, t) is a d × d matrix, and η(t) represents a Gaussian white noise sat- . Consider a current J in the stationary state of this system defined as J = x(t) is any arbitrary function of x, and • corresponds to the Stratanovich product. First we look at the implications of Eq. (1) for any such current. Without loss of generality, we consider currents for which J ≥ 0. Let [J(t)] k be the k-th cumulant of J [48]. Expanding the inequality in Eq. (1) in powers of (−λ) i , it can be shown that [J(t)] k ≥ 0 for k ≥ 2 for any t. For 1 < k ≤ 3, the cumulants coincide with the k−th central moment of J. Here · corresponds to an ensemble average over steady state trajectories of length t. An important standardized moment which can then be constructed is the skewness, . Skew-ness quantifies the asymmetries of the fluctuations about the average value. Here we focus on the properties of the skewness as a function of t. Using the results in Ref. [37] which proved the Gaussianity of general current fluctuations in the t → 0 limit (see Eq. S13 -S15 in the Supplementary Note 1 of [37], where it is shown that J(t) , and J(t) 2 are ∝ t for small t, but J(t) 3 ∝ t 2 for small t.
For J = ∆S tot , this behaviour was shown already in [33]), we obtain that S ∝ t 1/2 for small t. Further, the existence of the large deviation function Φ(λ) = lim t→∞ φ(λ, t) ensures that all the cumulants scale linearly in time as t → ∞. As a result, S ∝ t −1/2 for large t. Combining these two limiting behaviours with the positivity of the cumulants, we obtain that S is a positive, non-monotonic function of time that vanishes both in the t → ∞ limit as well as the t → 0 limit. As a consequence, there will also be a special time, where the skewness attains a maximum value. This is the first central observation we make in this paper. From the generality of the above arguments, we expect this behaviour to be generic for any non-equilibrium steady state current if the short-time and long-time behaviour are as detailed above. But to be more concrete, we now take the example of the non-equilibrium steady state of a colloidal system [24,36,[49][50][51][52]. The model consists of a single colloidal particle in a harmonic trap with stiffness κ, whose mean position is modulated according to the Ornstein-Uhlenbeck process. The dynamics of the system with position variable x(t), and the trap center x 0 (t) ≡ λ(t) can be described using a system of overdamped Langevin equations as: Here D is the diffusion constant at room temperature (T ), τ = γ/κ is the relaxation time of the harmonic trap, κ is the trap stiffness and γ is the drag coefficient related to D by the Stokes-Einstein relation as Dγ = k B T . Similarly, τ 0 is the relaxation time of the OU process and A corresponds to its strength. The noise terms ζ(t) and ξ(t) are Gaussian white noises obeying ξ(t) = 0, and ξ(t)ζ(t ) = 0. In Fig. 1, we plot the skewness S of several arbitrary currents as a function of time, using numerical data. Currents are constructed using J = Here {c i } ∈ R are constants which can be varied to construct different currents. Without loss of generality, we consider currents with J ≥ 0. We find that the skewness of arbitrary currents is a positive, bounded, and non-monotonic function of t which vanishes both in the short -time limit as t 1/2 and in the large-time limit as t −1/2 . Next, we investigate these properties for a particular choice of the current, which is J = ∆S tot [53]. The corre- sponding entropy production rate is given by σ = δ 2 θ (δ+1)τ0 , where δ = τ0 τ and θ = A D [24,50,52]. In Fig. 2a, we plot the analytically computed skewness of ∆S tot (t), for different values of θ [54]. As expected, we find that S(t) ≥ 0 and is non-monotonic in time, featuring a maximum at an intermediate time t ≡ τ E . We also find that the skewness increases with θ, which also increases the entropy production rate of the system. Now we look at the measurable consequences of this non-monotonicity. We first consider the mean fraction of the time when the measured entropy stays above the average value, denoted by T + = 1 t t 0 Θ(J(s) − J(s) ) ds, where Θ(x) is the Heaviside function. If the fluctuations of ∆S tot (t) were symmetric about the average, then it would have implied, T + = 1 2 . Indeed, this is known to be the case in the t → ∞ limit [31]. In Fig. 2b, we plot T + , obtained from the numerical data as a function of t. We find that, T + < 1 2 for all t and tends to 1 2 in the limits t → 0 and t → ∞. We also find that T + is non-monotonic in time, and attains a minimum close to when S(t) is the highest.
Interestingly, the behaviour of T + implies that it is more likely that a current measured for a finite t duration stays below the average value for most part of the measurement. This discrepancy is the highest when the skewness peaks. We also find that T + is monotonically decreasing as a function of θ for all t. Hence, the further away the system is from equilibrium, the more likely that arbitrary currents or entropy production measured along a single trajectory stay below the average value for most of the time.
where t sup is the time of global maximum of ∆S tot (t)−σt. The results are shown in Fig. 2c. We find that T max has a similar time dependence as T + and stays below 1 2 for all t and tends to 1 2 in the limits t → 0 and t → ∞. T max (t) is also found to be non-monotonic in time featuring a minimum at t ∼ τ E . This means, for currents measured for a finite time, it is more likely that it's maximum deviation above the mean will be found before the half-length of the measurement. Again, the discrepancy will be the highest when t ∼ τ E .
As seen in Fig. 1, the time at which the skewness has the highest value varies from current to current. In particular, for J = ∆S tot , from the analytical solutions, we find that τ E monotonically decreases with increasing θ (See the supplementary information where the dependence of τ E on θ and τ 0 is analytically obtained). As a result, in order to capture the non monotonic nature of S(t), the further a system is away from equilibrium, the finer the resolution needs to be. We demonstrate below that this is, however, still very detectable in experimental data.
Experiments: The model in Eq. (3) was first realized in an Optical Tweezers setup in Ref. [49] and was also studied recently in Ref. [36]. To realize this system, we trap a 3 µm polystyrene particle in an aqueous solution in a harmonic potential well given by is the time dependent mean position of the trap which is modulated using an acousto-optic modulator according to an Ornstein-Uhlenbeck process (see Eq 3) with τ 0 = 2.5 ms, and A = [0.1, 0.2, 0.3] × (0.6 × 10 −6 ) 2 m 2 /s. We sample the one-directional trajectory of the probe at a spatiotemporal resolution of ∼ 1 nm -10 kHz for 100 second. We then use the autocorrelation of the time series of a trapped particle and the noise to calibrate the fluctuation of the probe from volts to nm [55]. We plot the experimental results in Fig. 2d, e and f. We find that the numerical results in Fig. 2 (a)-(c) are reproduced by our experiments and that the non-monotonicities in S(t), T + and T max are clearly visible [56].
So far, we have considered a linear Langevin model. However, our results are generic and are expected to hold for any overdamped diffusive process. To substantiate this, we numerically explore the non-monotonic nature of fluctuations of entropy currents of a system with non-linearities. For this, we consider an anharmonic Brownian gyrator with a double-well confining potential (studied recently in [57]; see the supplemental material for details and for the additional example of a gyrator with a quartic confining potential). We vary the nonequilibrium conditions by changing the ratio of the temperatures (α = T2 T1 ) along the two orthogonal directions. In Fig. 3, we show that all the features that we showed in Fig. 2 are also present in this non-linear model.
Apart from the generic bound in Eq. (1) which holds for any continuous-time Markov process, a crucial ingre-   dient in our results is the emergence of Gaussian fluctuations in the short-time limit of overdamped diffusive systems, recently proved in [37]. For discrete-space systems, which evolve according to a continuous-time Markov process, current fluctuations are not necessarily Gaussian at t → 0; in fact it can be shown that all the moments scale ∝ t at short times [58]. Thus, the skewness will scale ∝ t −1/2 , diverging as t → 0, and will not necessarily be non-monotonic in time. As a consequence, shorter current measurements in such systems are more likely to lie below the average. In all the cases these effects increase when the system is further away from equilibrium. In Fig. 4, we demonstrate this for a three-level system, coupled to two thermal reservoirs at unequal temperatures [59][60][61]. Variants of this system have recently appeared in a number of different contexts [62][63][64], an important example being classical/ quantum clocks [65]. It is known that clocks need to run for a long time to give reliable estimates [66] which is also the limit in which the skewness vanishes.
In summary, we have unravelled and quantitatively characterized the universal properties of skewness of current fluctuations, and their measurable consequences in non-equilibrium steady states. The results imply that, rather counter intuitively, it is more probable that currents such as entropy production will mostly lie below average in a single realization of finite-time duration. Our results can be potentially verified in molecular motors such as kinesin [67] by looking at the statistics of it's steps [68] or energy dissipation [69]. Our results also show that the nature of the fluctuations of currents, and thus the most probable outcome, crucially depend on the time duration of the measurement. It would be of substantial interest to investigate whether biological motors have optimized their timescales considering such constraints (through molecular evolution). For artificial microscopic engines such as the ones studied in Refs. [70][71][72][73][74][75], this implies that it might be possible to choose cycling times to optimize fluctuations, and to get a reliable performance in a limited number of runs. We plan to attempt answering these questions in our future research.

SUPPLEMENTAL MATERIAL
In this section, we reproduce the exact calculation of the moment generating function of ∆S tot (t) for the Stochastic Sliding Parabola model, previously obtained in [52].

I. EXACT CALCULATION OF THE MGF OF ∆Stot(t) FOR THE STOCHASTIC SLIDING PARABOLA (SSP) MODEL
The stationary probability distribution for x and λ is given by [50], the MGF of total entropy production ∆S tot (t) can be written down in the following manner. First, the joint probability density functional of trajectories starting at t = 0 at (x 0 , λ 0 ) and ending at t = τ at (x t , λ t ) may be written as, ds L(ẋ(s), x(s),λ(s), λ(s), s) (5) with the Lagrangian, The normalization constant N for this case is [78], The entropy production in the steady-state in the time interval [0, τ ] for the SSP is then, This form of the entropy production can easily be obtained by equating it to the ratio of the probabilities of forward and time-reversed trajectories using Eq. (4) and Eq. (5) and the form of the Lagrangian Eq. (6). Hence, upto a normalization factor C (determined by Eq. (4) and (7)), we have the following expression for the MGF of ∆S tot (t), with the augmented action S[ x(·), λ(·), u ] = (δ + 1) δ 2 θ(x 0 − λ 0 ) 2 + δ θx 2 0 + λ 2 0 + λ 2 0 2Dτ 0 θ (δ 2 (θ + 1) + 2δ + 1) after several partial integrations, it can be shown that the above quadratic action reduces to where the kernel is defined by the operator: Carrying out the Gaussian integral, and requiring the boundary terms to vanish, the generating function at arbitrary times τ can be written down as a ratio of functional determinants, This ratio can be computed using a technique described in [79] and used in [51], which is based on the spectral -ζ functions of Sturm-Liouville type operators. Applying this method, it can be shown that this ratio can be obtained in terms of a characteristic polynomial function F as, where, H is the matrix of suitably normalized fundamental solutions of the homogeneous equation, A u x = 0, and is defined as, M and N have information about the boundary conditions from Eq. (11) and we require, A derivation of Eq. (14), applicable to a class of driven Langevin systems with quadratic actions is given in [51]. We would also like to stress that, the expression given in Eq. (14) is valid only for u ∈ [u − (τ ), u + (τ )] for which the operator A u doesn't have negative eigenvalues. The MGF is not analytic outside this interval. For the SSP in the steady state, we find, the four independent solutions of A u x = 0 to be where, Matrices M and N are given by, Using these, the MGF can be computed exactly using Eq. (14), and various moments of the probability distribution can also be exactly obtained for any t.

A. Dependence of τE on system parameters
For the stochastic sliding parabola model, it is possible to analytically compute how τ E , the time at which Skewness peaks, depends on the system parameters, using the exact expressions. In Fig. S1, we demonstrate how τ E depends on θ and τ 0 . We find that τ E monotonically decreases with increasing θ, which is the limit at which the system goes further away from equilibrium. This means, further away the system is from equilibrium, a higher resolution in time will be required to see the nonmonotonicity in current fluctuations. The dependence of τ E on τ 0 is found to be more non-trivial. We find that τ E is a non-monotonic function of τ 0 and features a minimum at a particular τ 0 value. FIG. S1. Dependence of τE on the relative magnitude of the Ornstein Uhlenbeck driving θ. We find that τE monotonically decreases with increasing θ, which is the limit at which the system goes further away from equilibrium. b) Dependence of tE on τ0. We find that τE is a non-monotonic function of τ0 and features a minimum at a particular τ0 value.

II. ENTROPY CURRENTS OF ANHARMONIC BROWNIAN GYRATORS
The brownian gyrator is one of the minimal prototypes of a microscopic heat engine [72,73,75]. It consists of a micron-sized particle trapped in a generic potential well, coupled to two heat reservoirs with different temperatures along two orthogonal directions. When the two degrees of freedom of the trapped particle are coupled, it can be shown that the system reaches a non-equilibrium stationary state, where the particle starts gyrating around the minima of the potential. The dynamics of the system, in the overdamped limit, can be expressed in terms of coupled langevin equations: Here, U (x, y) denotes the confining potential in the x − y plane. The x axis is coupled to a thermal reservoir at temperature T 1 and the y axis is coupled to another thermal reservoir at temperature T 2 . The corresponding thermal noises η i (t) are of Gaussian nature and white in time, such that η i (t) = 0 and η i (t)η j (t ) = δ ij δ(t − t ). The viscous drag of the medium is denoted by γ, which is related to the temperatures of the reservoirs through the Einstein relation, D i γ = k B T i , k B is the Boltzmann constant (We set k B = 1 for simplicity). In this work, we choose the anharmonic Brownian gyrator, where the confining potential is anharmonic [57,74], as an example of a non-linear diffusive model.
We first consider a Brownian gyrator with a double-well confining potential of the form, Where axes of the potential x and y are rotated by an angle θ with respect to the temperature axes (x, y) [ axes of the co-ordinate frame] as, The parameter 'b' can be used to tune the bi-stable nature of the potential along the x direction as the barrier height (= b 2 ) and the position of the minima ((= ± √ b) ) of the potential are dependent on it. The stiffness constant 'k' characterises the harmonic part of the potential along the y direction. We choose b = 1, k = 2, and θ = 45 • for the analysis performed in this work.
We also consider another anharmonic Brownian gyrator with a quartic confining potential given by, The parameters 'k 1 ' and 'k 2 ' can be considered as the stiffness constants along the respective directions. We set k 1 = 1, k 2 = 2 along with θ = 45 • to construct an anisotropic quartic potential for the analysis performed in this work.
The entropy currents for both gyrators are constructed numerically with third-order polynomial basis functions using the TUR-based short-time inference scheme [33][34][35], as described in Ref. [57] [80]. We have looked into the properties of entropy currents corresponding to the different non-equilibrium conditions controlled by the ratios of the temperatures (α = T 2 /T 1 ) of the thermal reservoirs of the systems, where α = 0.1 corresponds to the most non-equilibrium configuration we considered, and α = 0.3 is the least non-equilibrium configuration of both systems. The non-monotonic features of the entropy current fluctuations for the gyrator system with double-well confining potential are shown in Fig. 3 of the main text and in Fig. S2 we show the similar results for the quartic well confining potential.

III. SKEWNESS OF CURRENTS IN CONTINUOUS TIME MARKOV JUMP PROCESSES
Here we demonstrate that Skewness ∝ t −1/2 in the t → 0 limit for discrete space system evolving according to a continuous time Markov jump process. We consider a set of M number of states {i}, i = 1, 2, ...M , and transition rates Γ ij ≥. Let π(i) correspond to the steady state probability of finding the system in the state i. A stochastic realization of the system is denoted by x(s) ∈ {i} with s ∈ [0, t). The fluctuating current between any two pairs of states i and j can be computed as, where x(s + k ) (x(s − k )) corresponds to the state of the system immediately after (before) the transition at times s = s k . A generalized time-integrated current in this system is given by, , where d ij = −d ji are weighting factors which are constants. The steady state average of this current is given by, A particular choice, d ij = F ij = log Γij π(i) Γjiπ(j) corresponds to the current J = ∆S tot (t).
The fluctuations of any such current J can be calculated using it's moment generating function, where L(λ) is the tilted transition matrix with elements, We are particularly interested in the small time properties of the moments of the current J(t). This can be obtained by Taylor expanding G as near t = 0. Keeping to first order in t, we obtain, For a generic choice of Γ ij and d ij , it is possible to verify that all the moments of J, obtained by Taylor expanding G as a function of λ, will be proportional to t for small t. Thus the skewness as defined in the maintext, will be proportional to t − 1 2 for t → 0. In Fig. (4)a of the maintext, we show this for a three-level system x ∈ {0, 1, 2}, with energy levels E i = {0, E 1 , E 2 }. We have assumed that the transitions between the levels 0 and 1, and the levels 0 and 2 are mediated by a hot reservoir at inverse temperature β 1 = 1 k B T1 , where k B is the Boltzmann constant and T 1 is the temperature of the hot reservoir. The corresponding transition rates obey the local detailed balance condition: The transitions between the levels 1 and 2 are assumed to be mediated by a cold reservoir at inverse temperature β 2 = 1 k B T2 > β 1 . The corresponding transition rates obey The skewness of arbitrary currents, and in particular entropy production can be straightforwardly computed for this model using the expressions given above. There are also standard techniques available to numerically simulate this model as a