Work fluctuations for diffusion dynamics submitted to stochastic return

Returning a system to a desired state under a force field involves a thermodynamic cost, i.e. work. This cost fluctuates for a small-scale system from one experimental realization to another. We introduce a general framework to determine the work distribution for returning a system facilitated by a confining potential with its minimum at the restart location. The general strategy, based on average over resetting pathways, constitutes a robust method to gain access to the statistical information of observables from resetting systems. We exploit paradigmatic setups, where explicit computations are attainable, to illustrate the theory. Numerical simulations validate our theoretical predictions. For some of these examples, a non-trivial behavior of the work fluctuations opens a door to optimization problems. Specifically, work fluctuations can be minimized by an appropriate tuning of the return rate.


Introduction
Stochastic resetting is a relatively recent, fast-growing field within the realm of nonequilibrium statistical mechanics.Despite the simplicity of the concept, outstanding properties emerge when, on top of its natural dynamics, a stochastic system is submitted to a random reset to a certain state.Due to its rich phenomenology, stochastic resetting still represents an appealing test bench for researchers interested in non-equilibrium even a decade after its original modern formulation [1].
Stochastic resetting exhibits remarkable theoretical properties that have made the former an excellent tool to analyze non-equilibrium stationary states [2,3,4,5], thermodynamic speed limit [6], thermodynamic uncertainty relation [7], inducing and optimizing the Markovian Mpemba effect [8], antiviral therapies [9], modelling cell division [10], tax dynamics [11], just to cite a few hot topics in non-equilibrium statistical mechanics assisted by stochastic resetting.However, it is not only in the theoretical perspective that stochastic resetting becomes popular, there are certainly a myriad of applications, from search and foraging processes [1,12,13,14] to ecological disasters [15] passing though enzymatic reactions [16], see [17] for a thorough review.
In its original formulation [1], stochastic resetting was considered instantaneous, i.e., when the reset occurs the state of the system is instantaneously changed to a resetting state and the natural dynamics starts afresh.From a theoretical perspective, the simple renovation of the dynamics allows to resort to a renewal formulation of the dynamics [18,19].Nevertheless, instantaneous reset is an ideal limit, which may be inconvenient to characterize real systems.For instance, a clear drawback of instantaneous resets is the thermodynamic cost.The cost required to drastically change the state of a physical state in a vanishing interval of time is arbitrarily large.Indeed, some proposals have been introduced to circumvent this unsuitable feature: refractory periods [20,21] and intermittent potentials [22,23].The former considers an instantaneous reset that is followed by a residence time in the resetting state, whereas the latter considers a confining potential that randomly switches on/off in order to avoid the system to get away from the resetting state.Notwithstanding, both models fails someway.In a reset with refractory periods, the reset event itself is still instantaneous, while intermittent potential does not guarantee the reset.
Return dynamics solves the problem [13,24,25,26,27,28].In stochastic returns, the dynamics comprises two evolution phases: the natural dynamics and the return dynamics.The first one is randomly interrupted by switching it to the latter, which ends only when the resetting state is reached.Return phases have been considered to be either dynamically deterministic or stochastic.This strategy guarantees the effectiveness of realistic reset within a finite time.For the moment, experimentation has been scarce, although recent developments have allowed to start implementing stochastic resetting concepts in real experiments [29,30].
Stochastic thermodynamics [31] addresses the study of thermodynamic quantities (work, heat, entropy, . . . ) from the point of view of statistical mechanics, usually within a mesoscopic scale where fluctuations may be relevant.Surprisingly, to the best of our knowledge, there is not much literature on the stochastic resetting viewed under the lens of stochastic thermodynamics, apart a few exceptions [32,33,34,35,36].Specifically, the cost in terms of mechanical work employed to implement the resetting dynamics have never been looked into deeply, despite its evident interest from both theoretical and experimental perspectives.To fill in this gap, in this paper we aim to present a general method to compute the work fluctuations associated with stochastic return process facilitated by a confining potential.
Optimality have been addressed in the context of stochastic resetting [12,37].Nevertheless, optimization problems have been posed traditionally for minimizing first passage time, usually motivated by applications to search processes.This is reasonable since no notion of cost was available before.Nevertheless, once the cost is introduced into this work, it is interesting to pose questions regarding optimization of such a cost or its fluctuations with respect to the resetting parameter.
The rest of the articles is organized as follows.Section 2 is devoted to the presentation of the model system that is a one-dimensional Brownian particle submitted to stochastic return.The central result of the article, which is the computation of the distribution of the mechanical work carried out the system, is derived in Sec. 3. The results are explicitly particularized in Sec. 4. Therein, paradigmatic setups, Poissonian return either with V-shaped or harmonic potential, are used to illustrate the excellent agreement between theory and simulations.Finally, we deliver conclusions and some perspectives in Sec. 5. Some detailed calculations are relegated to the Appendix.

The model
We consider a Brownian particle freely diffusing in a one-dimensional space with diffusion constant D. Thus, the time evolution of the particle position, x(t), follows the Langevin equation: where the dot indicates a time-derivative, and η(t) is a Gaussian white delta-correlated thermal noise with zero mean and unit variance, which fulfills η(t)η(t ) = δ(t − t ).At a random interval of time, drawn from a specified distribution, f (t), a confining potential, U r (x), with its minimum at the origin, is switched on.Once the potential is on, the dynamics reads For the sake of simplicity, the same diffusion constant D as in Eq. ( 1) has been considered in Eq. (2).Nonetheless, results shown below can also be carried out for two different coefficients in the spirit of Ref. [26] along the same line.The external potential is switched off when the system does the first passage to the origin, which is the minimum of the potential U r (x), and then, the system starts afresh the dynamics generated through Eq. (1) during a random time distributed according to f (t).Then, dynamics in Eq. ( 2) starts up to the return to the origin and the same game is played iteratively.In summary, the system explores the space according to Eq. (1) (exploration phase), whereas it returns to the origin (the resetting location) following Eq.(2) (return phase).Notice that the above dynamics can be cast using only one Langevin equation: for U (x, λ(t)) ≡ λ(t)U r (x), where λ is a dichotomous controlled variable that switches between zero and unity for exploration and return phases, respectively.The stationary probability density function for the position emerging from this implementation of stochastic resetting has been studied in Ref. [26], and the relaxation to the steady state is discussed in [38].In this article instead, we focus on the thermodynamic properties of the system.Specifically, the distribution of the work, W tot , required to reset the system under the non-instantaneous resetting protocol discussed above up to a certain observation time t is looked into.
On a general basis, within the framework of stochastic thermodynamics [31], the work in a dynamical process where the potential U (x; λ(t)) is varied through a control parameter λ(t) can be simply defined, For the model under consideration: (i) the aim of the external potential is to bring the system to the resetting location, which coincides with the minimum of the potential, (ii) contributions to the work are instantaneous since λ is different from zero just in a null set of times, which are those where the exploration phase is switched to the return phase or vice versa.For convenience, we consider min{U (x; λ)} = 0, that has no physical implications since it just defines the energy origin.
As introduced above, the function λ(t) is piece-wise constant.Let us define t i and τ i as the times when, respectively, the i-th exploration and return phase end.Consistently, these times are stochastic.Specifically, the duration of the i-th exploration phase, t i − τ i−1 (τ 0 ≡ 0) has to be distributed according f ; whereas the duration of the i-th return phase τ i − t i , follows a first passage distribution, which will be discussed later in detail.Hence, λ(t) jumps from 0 to 1 at times t i and does the opposite at times τ i .Therefore, it is possible to explicitly write down: where Θ(•) is the Heaviside theta function.The sum takes into account formally infinite events but if one is interested in a finite time window, it suffices to sum those within it.It is important to recall that the system starts in exploration phase, thus Differentiating Eq. ( 5) with respect to time, and substituting λ into Eq.( 4) yields where n r (t) counts the number of return phases started up to time t, and we have taken into account that U r (x(τ j )) = 0 since x(τ j ) is the location of the minimum of the potential.Notice that x(t i ), in Eq. ( 6), is the position of particle where the resetting phase starts.Equation ( 6) provides the work performed in a single trajectory.Of course, this is a stochastic quantity because of the underlying stochasticity of the dynamics.The rest of this article is devoted to obtain the distribution of the work.

General theoretical framework
Since the work in a single trajectory is written in terms of the position of the particle when return phases start, it will be required to analyze the statistics of what we call the resetting pathway.We define the resetting pathway as the specific resetting history followed by the Brownian particle, comprising the information concerning times where .Sample trajectory of a Brownian particle submitted to stochastic return to location x 0 .The evolution alternates between the free exploration phase (blue solid line) where the system follows its natural dynamics, and the return phase (red solid line), in which a biased drift towards x 0 is introduced.The resetting pathway is defined by the set of times where a specific phase ends, either exploration, t i or return τ i , along with the positions where the return phases start, the phase is switched and position of the particle therein.Figure 1 schematically illustrates the resetting pathway of a sample evolution.

Statistical weight of resetting pathway
This article aims at deriving the work distribution of the resetting system under study.The strategy to do so is to sum the contributions stemming from all possible resetting pathways weighted consistently with their probability.For the sake of clarity, let us first introduce some useful notation: (i) f (t) is the probability density function of the resetting time intervals.Defining is the probability of not having started the return phase after an exploration phase of duration t.
(iii) p(x, t|x 0 ) is the probability density function of the position of a particle evolving an amount of time t exclusively through the diffusion defined by the exploration phase given that it started at x 0 .This is the solution of the Fokker-Planck equation associated to Eq. ( 1) with an initial condition p(x, 0|x 0 ) = δ(x − x 0 ).
(iv) π(x, τ |x r ; x 0 ) is the probability density function of the position of a particle evolving an amount of time τ exclusively through the return phase given that it started at x r and there is an absorbing boundary at x 0 .This is the solution of the Fokker-Planck equation associated to Eq. ( 2) in the returning phase with initial condition π(x, 0|x r ; x 0 ) = δ(x − x r ) with absorbing boundary condition at x 0 .To avoid cluttered formulae, in what follows we get rid of x 0 explicitly in the notation.For instance, this implies that for π and p, that is, we will write simply π(x, τ |x r ) and p(x, t) respectively.Note that we have chosen x 0 as the origin of coordinates indeed.
(v) φ(τ |x r ) is the first passage time distribution to x 0 in the return phase given that the system starts from x r .If a return phase starts at time t i from position x r , the variable τ i − t i is distributed according to φ.
) is the probability of not having ended the return phase started a time τ ago, that is, the survival probability, Φ( Note that π is zero for sign(x − x 0 ) = sign(x r − x 0 ), since the absorbing boundary forbids the transfer of the particle to the other side of the boundary.
The above statistical objects are the building blocks to characterize the statistical weight of a given resetting pathway.Let ψ p n (t) be the probability of having n finished trial and being currently in the phase of kind p = {diff, ret}.(A finished trial means the completion of a return phase followed by an exploration phase.)Combining carefully the stay probabilities of each phase, and defining χ(x, t) ≡ f (t)p(x, t) for the sake of convenience, the general formulas for these probabilities read: and The integrand on the right-hand side of the above expressions are the probability density function for a specific resetting pathway.Of course, normalization holds since the sum, plus all integrals inside the explicit form of ψ, cover all possible resetting pathways.For a detailed discussion regarding the construction of the above probabilities as well as for a demonstration of its normalization, see Appendix A. Functions ψ can be written in a compact form taking into account that they have a convolution structure.Using ' * ' to denote the convolution, i.e., [A * B](t) ≡ t 0 dt A(t )B(t − t ), between functions and as exponent, e.g., A * 2 (t) = [A * A](t), for convolution power, it is possible to get and

Work distribution
The convolution structure found out in the previous subsection 3.1 makes especially convenient to study the distribution of the work through the Laplace transform of its moment generating function.Let us start by defining the moment generating function: where the notation • stands for the weighted average over all possible resetting pathway.
Since the total work is just the sum of contributions of U r evaluated at the particle's positions where the return phases start, we have that with χ W (k, x, t) ≡ χ(x, t)e kUr(x) .
Introducing the Laplace transform, we can write explicitly the Laplace transform of the moment generating functions in terms of the Laplace transforms of the functions characterizing the probabilistic ingredients of the model.Above, we have used the convolution theorem for Laplace transform, carried out the sum of the geometric series and recalled the relation 15) is the central result of this study.It is an exact general result for the distribution work of a system submitted to return dynamics with arbitrary return potential U r and arbitrary return rate f .

Results for diffusion under Poissonian return with paradigmatic potentials
From now on, we consider the case of Poissonian return, f (t) = re −rt , with r being a constant resetting rate.This choice allows to give some explicit simplifications in Eq. (15).Specifically, we get where is the Laplace transform of the free propagator of the diffusion phase, p(x, t) = exp [−x 2 /(4Dt)] / √ 4πDt.In the following, we consider two specific paradigmatic potentials: V-shaped and harmonic.

V-shaped potential
We consider first a V-shaped potential, U r (x) = γ V |x| in the return phase.This is one of the quite few cases, where the first passage density, φ(t|x), for the particle to reach the origin starting from location x can be obtained exactly, see section 3.2.2.2 from Ref. [39], and its Laplace transform reads Substituting Eqs. ( 17) and ( 19) into Eq.( 16) and performing the integrals over x, where, for the sake of simplicity, we have introduced the parameter α V ≡ γ 2 V /4D, which characterizes the inverse of a time-scale, defined by the relative intensity of the confining potential with respect to the diffusion.
Moments can be obtained from the moment generating function (15).Specifically the n-th moment of the work in the Laplace space is We have carried out the explicit calculation for the first and second moments.On the one hand, the expression obtained exactly are not especially illuminating and then we relegate them to Appendix B. On the other hand, the limit of large times is quite informative and fairly simple.Thus, we develop it here.To look into large times, it is needed to study the Laplace transform for small s.Specifically, the Laurent series for the two first moments read: with µ Application of the generalized final value theorem allows to study the asymptotic behavior of the first moments.The mean work diverges linearly with time, whereas the second moment does it quadratically Nevertheless, the variance of the work, /2. Provided the asymptotic behavior for the mean and the variance, one trivially gets that, for the square of the coefficient of variation, Off course, the coefficient of variation without the t-rescaling, goes to zero for large times.Thus, the relative fluctuations being negligible in such regime.The theoretical predictions derived for the V-shaped potential are compared to results obtained from averaging over simulated trajectories in Fig. 2. Therein, we have obtained the theoretical evolution of the moments by numerical inversion of the Laplace transform of the moments, whose analytical expression can be found inAppendix B. The Rescaled mean (a), variance (b), and the square of the coefficient of variation (c) for a particle submitted to V-shaped potential under return dynamics.Theoretical predictions for the time evolution stems from numerical inversion of the Laplace transform of the moments (see Appendix B for their analytical forms) whereas the asymptotic values comes from Eqs. ( 28), ( 30) and (31).The parameters considered here are r = 0.3, D = 2, and α V = 1.Simulation results stand for average over 10 5 trajectories using a time step dt = 10 −3 .
prediction of the asymptotic value for large times is obtained directly from Eqs. ( 28), (30) and (31).The agreement is excellent, not only for the asymptotic regime but also in the time evolution predicted by the numerical inversion.Let us discuss in detail the dependence of the asymptotic values obtained above on the return rate r.During the whole discussion, we make reference to the t−rescaled moments appearing in Eqs.(28), (30) and (31).The rescaled mean work is an increasing function of r that goes from W /t = 0 for r = 0, and possesses a horizontal asymptote at W /t = 4Dα V for r → ∞.Both, the variance and the coefficient of variation are non-monotonous functions of r.Both of them start decreasing from r = 0 up to find a minimum value and then increase up to reach a horizontal asymptote (see Fig. 3) The values of the return rate r min at which the minimum is attained are for σ 2 W /t and tσ 2 W / W 2 , respectively.This is a remarkable feature from the optimization point of view.The dependence of the work fluctuations on the return rate are non-trivial, exhibiting a minimum value for the aforementioned values either in absolute or relative terms.

Harmonic potential
We consider now a harmonic potential, U r (x) = γ h x 2 /2 in the return phase.This is another simple case, where the first passage density, φ(t|x), for the particle to reach the origin starting from location x can be obtained (see for instance Appendix A.2 in  30) and ( 31) respectively.The vertical red dashed lines stand for the optimal values minimizing the fluctuations in Eqs. ( 32) and (33).As in Fig. 2, D = 2 and α V = 1.Regardless the parameters, the optimal resetting rate minimizing the coefficient of variation is roughly four times the one minimizing the variance.
Ref. [26]).The result in the Laplace domain is where Γ and U stand for the gamma function and the confluent hypergeometric function, also known as Tricomi's function, respectively.Obtaining the Laplace transform of the moment generating function for the work reduces to substitute Eqs. ( 17) and (34) into Eq.( 16) with U r (x) = γ h x 2 /2 and performing the integrals over x.
Unlike the case of V-shaped potential discussed in the previous subsection 4.1, the analytical computation of the first and second moment is not straightforward.Nevertheless, it is possible to compute these quantities in the long-time limit numerically.Let us briefly sketch the procedure for the computation of the scaled moments.For a given Laplace transform variable s, we compute the first and second moment (in the Laplace space), respectively, by taking first and second order derivative of the numerical expression G W (k, s) with respect to k and set those derivatives equal to zero [see Eq. ( 21)].By numerically inverting the Laplace transform, we can compute the first and second scaled cumulant in the long-time limit.
In Fig. 4(a-d), we show the comparison of these scaled cumulants obtained using Langevin simulations for different sets of parameters.Furthermore, we show the comparison of the long-time analytical results with long-time numerical simulations results.Clearly, we see that the scaled cumulants becomes independent of time similar to what we have observed in the previous subsection 4.1 for the case of V-shaped potential.Further, we notice from the right panel that the scaled cumulants of the work increase (decreases) with the stiffness parameter, λ (resetting rate, r) in the long-time limit.

Conclusions
As far as we are aware, this is the first time that the distribution of the mechanical work performed by a resetting mechanism is worked out.Our framework, relying on the average over resetting pathways, leads to write the distribution-the Laplace transform of the moment generating function to be more accurate-in terms of statistical information of the processes that constitute the return dynamics, see Eq. (15).Note that this approach remains completely valid for both arbitrary resetting time distributions f (t) and arbitrary dynamics for the diffusion and return phases.Furthermore, not only is the framework useful for computing the work distribution but any physical observable depending on the resetting pathway.
The power to obtain exact results through this method has been shown through explicit computation in paradigmatic cases.Specifically, a V-shaped potential, which confine particles by submitting to constant force to both sides of the confining point, and harmonic potential have been considered under Poissonian resetting times.
For the first one, explicit calculations are manageable, being possible to reach a closed expression for the Laplace transform of the moment generating function.Hence, Laplace transform of the moments is obtained exactly and validated with simulation results pleasurably after inverting numerically the Laplace transform.The asymptotic behavior for large times is fully resolved analytically for the first moments.An optimal resetting rate emerges when analyzing the fluctuations of the work.This is a remarkable feature characterizing the non-trivial behavior of the work.
Regarding the harmonic case, the explicit calculation is not so straightforward as for the V-shaped one.Nevertheless, numerical evaluation of the predictions and comparison with simulation results have validated our approach in this more involved case.
The methods introduced here are expected to open new insight to stochastic thermodynamics in resetting systems, allowing to look into the energetics of the resetting mechanism and its cost.Optimal features, as the one shown for the V-shaped potential, are especially appealing since they are the hallmark for engineering efficient resetting mechanisms.

Figure 1
Figure1.Sample trajectory of a Brownian particle submitted to stochastic return to location x 0 .The evolution alternates between the free exploration phase (blue solid line) where the system follows its natural dynamics, and the return phase (red solid line), in which a biased drift towards x 0 is introduced.The resetting pathway is defined by the set of times where a specific phase ends, either exploration, t i or return τ i , along with the positions where the return phases start, x

Figure 3 .
Figure 3. Asymptotic rescaled fluctuations as a function of the resetting rate.Variance (a) and the square of the coefficient of variation (b) are non-monotonous functions as given by Eqs.(30) and (31) respectively.The vertical red dashed lines stand for the optimal values minimizing the fluctuations in Eqs.(32) and(33).As in Fig.2, D = 2 and α V = 1.Regardless the parameters, the optimal resetting rate minimizing the coefficient of variation is roughly four times the one minimizing the variance.

Figure 4 .
Figure 4. First and second scaled cumulant of work for stochastic returns facilitated by a harmonic trap.Left panel (a-d): Figures are obtained using Langevin simulations.Right panel (a-d): Figures show the plateau points of each figure in the left panel, and show the comparison of analytical results with Langevin simulation in the long-time time.(a-b): Resetting rate r = 0.5.(c-d): Trap stiffness γ h = 0.5.In all plots, the diffusion constant D = 0.75 is kept fixed.The numerical simulation is performed for a time increment dt = 10 −3 , and the averaging is performed over 10 4 realizations.