Work statistics at first-passage times

We investigate the work fluctuations in an overdamped non-equilibrium process that is stopped at a stochastic time. The latter is characterised by a first passage event that marks the completion of the non-equilibrium process. In particular, we consider a particle diffusing in one dimension in the presence of a time-dependent potential U(x,t)=k|x−vt|n/n , where k > 0 is the stiffness and n > 0 is the order of the potential. Moreover, the particle is confined between two absorbing walls, located at L±(t) , that move with a constant velocity v and are initially located at L±(0)=±L . As soon as the particle reaches any of the boundaries, the process is said to be completed and here, we compute the work done W by the particle in the modulated trap upto this random time. Employing the Feynman–Kac path integral approach, we find that the typical values of the work scale with L with a crucial dependence on the order n. While for n > 1, we show that ⟨W⟩∼L1−n expkLn/n−vL/D for large L, we get an algebraic scaling of the form ⟨W⟩∼Ln for the n < 1 case. The marginal case of n = 1 is exactly solvable and our analysis unravels three distinct scaling behaviours: (i) ⟨W⟩∼L for v > k, (ii) ⟨W⟩∼L2 for v = k and (iii) ⟨W⟩∼exp−(v−k)L for v < k. For all cases, we also obtain the probability distribution associated with the typical values of W. Finally, we observe an interesting set of relations between the relative fluctuations of the work done and the first-passage time for different n—which we argue physically. Our results are well supported by the numerical simulations.

Recently, there has been a growing interest in the study of the thermodynamic quantities until a particular event of interest has occurred for the first time [23,[37][38][39][40][41][42][43][44][45][46][47][48].Examples include the escape of a colloidal particle from a metastable state or stretching of a polymer till a certain length is attained [41].Since the underlying dynamics is stochastic, the time at which these events take place also varies from realisation to realisation.This drastically changes the properties of the thermodynamic quantities compared to situations where the observation time is fixed for all realisations.For instance, the bound on the average work picks up a non-trivial correction term due to the fact that the system is generally out of equilibrium at the end of the first-passage time [41].Similar observation has also been made for the efficiency of heat engines [40].In fact, based on the martingale theory, many new results on the integral fluctuation theorem and stopping times for entropy production have been analytically obtained [49].Despite these general results, the number of thermodynamic first-passage problems that have been solved exactly, seems to be very limited.In this paper, we aim to partially fill this gap by studying an analytically tractable class of models, where one can get exact results for the moment generating function associated with the mechanical work.
In order to study this, we reformulate the work done as a first-passage functional of the stochastic process [50].Statistical properties of such functionals have been quite extensively studied in the literature by deploying the celebrated Feynman-Kac formalism suitably adapted for first-passage problems [51].Based on this formalism, these functionals have been shown to have many applications in fields ranging from queue theory [52], sandpile and percolation models [53,54] to disordered systems [55], among others [56][57][58].Moreover, they have been studied for diverse stochastic processes such as diffusion and drift-diffusion [59][60][61][62], random acceleration [63], Lévy process [64], Ornstein-Uhlenbeck particle [65,66] and resetting processes [67,68].In this paper, we are interested in using these tools and techniques to calculate the moments and the distribution of the work done by a diffusing particle subjected to a time-dependent potential U (x, t) = k n |x − vt| n with order n > 0. In particular, we show that the work can be reformulated as a first-passage functional and one can then employ the Feynman-Kac formalism to investigate its statistical properties.
The rest of the paper is structured as follows: In Section 2, we introduce our model and present a general formalism to study the work done.We first focus on the analytically tractable cases of the linear potential (n = 1) in Section 3 and the harmonic potential (n = 2) in Section 4. The insights gained from these two cases become instrumental in dealing with the general n in Section 5. We then discuss an interesting relation between the work done and the first-passage time in Section 6.Finally, in Section 7, we present the conclusion and make some future remarks.

Model and general formulation of the problem
We consider a one-dimensional diffusing particle whose position, in presence of a timedependent potential U (x, t), evolves according to an overdamped Langevin equation where η(t) is the Gaussian white noise satisfying the properties ⟨η(t)⟩ = 0 and ⟨η(t)η(t ′ )⟩ = δ(t − t ′ ).Also D and γ denote the diffusion and damping coefficients, respectively, and are related by the Einstein relation D = k B T /γ.Throughout this paper, we will fix γ = 1.The initial position of the particle will be denoted as x 0 .
Moreover, we choose the potential U (x, t) to be a confining potential with stiffness k (> 0) and we let it move with a constant velocity v (> 0): Finally, we will assume the existence of two absorbing walls located at positions L ± (t) that themselves move with the same constant velocity v, i.e.L ± (t) = ±L + vt, with −L < x 0 < L. Plots of the potential U (x, t) and the associated force F (x, t) are provided in Figure 1.We observe that for n < 1, the magnitude of the force decays with increasing x, while for n > 1, it increases with x.Intriguingly, we observe different striking behaviour depending on the order of the potential strength as will be illustrated in this paper.
As mentioned in the introduction, we are interested in the work distribution of first-passage time problems.Here, we allow the particle to move until it touches one of the absorbing boundaries at L ± (t).Let us denote this time by t f and the corresponding trajectory as {x(τ )} with 0 ≤ τ ≤ t f .We then calculate the work done up to time t f as [4,[69][70][71] with Since the motion is random, we get different values of t f for different realisations.This indicates that in contrast to the usual case of fixed observation time, the definition of work in Eq. ( 3) possesses two sources of stochasticity -one arises due to the random trajectory {x(τ )} and the other stems from random first-passage time t f .Notice that we have considered the absorbing boundaries to be dynamical in our model.Such scenarios might arise in experiments in which a Brownian particle is dragged using a moving potential [69,70].Typically, they have been carried out for a fixed observation time.However, one can also consider scenarios where one waits till the particle comes out of the field of view of the camera or escapes from the trap.Now if the camera/trap is moving, the field of view also changes with time and these, in turn, mimic a moving boundary scenario.Another example might be investigating the statistics of work for stretching a polymer to a certain threshold length that itself depends on time [41].This can also be formulated as a moving boundary problem.
In order to derive the statistics of W (x 0 ), we reformulate the work done as a firstpassage functional which has been extensively studied in the literature [50].To see this, let us first perform a change of variable as and recast the Langevin equation and the work done in terms of this new co-moving variable: Meanwhile, the first-passage time t f , in the co-moving frame, simply becomes the time at which the particle moves out of the fixed interval [−L, L] for the first time.Therefore, by suitably defining the variable ξ(t), we have been able to recast the entire problem with fixed absorbing walls.Following [50,51], one can now derive a backward differential equation for the moment generating function with ξ 0 = ξ(t = 0) = x 0 and the averaging ⟨...⟩ involves averaging both with respect to the trajectories as well as the first-passage time t f .To do this, we now look at a trajectory of ξ(τ ) from τ = 0 to τ = t f and break it into two parts: (i) a left interval [0, ∆t] and (ii) a right interval [∆t, t f ] with ∆t → 0. At the end of the left interval, the variable ξ(t) takes the value ξ(∆t In the remaining time interval (t f − ∆t), the particle starting from ξ(∆t) moves out of the interval [−L, L].Therefore, decomposing the integral in Eq. ( 8) as ´tf 0 = ´∆t 0 + ´tf ∆t , we obtain Inserting the expression of ξ(∆t) and performing the expansion in ∆t, we obtain the following differential equation for Q(p, ξ 0 ) as ∆t → 0: where f n (ξ 0 ) = sgn(ξ 0 ) ξ 0 n−1 .This equation has to be solved in the domain −L ≤ ξ 0 ≤ L and it a differential equation in the initial position ξ 0 .This approach is referred to as a"backward" approach as it involves varying the initial position instead of the final one.Moreover it is a second order differential equation.Thus, we need two boundary conditions to solve it.These conditions can be obtained by analysing the behaviour of Q(p, ξ 0 ) at ξ 0 = ±L.For these two cases, the particle is already at the edge of the interval and gets instantly absorbed.Thus, for both cases, we have t f = 0 which implies W (±L) = 0 from Eq. (7).Plugging this in Eq. ( 8), we obtain The goal now is to solve the backward Eq. ( 11) with these boundary conditions and then utilize it to obtain the moments of W (x 0 ) as One can also obtain the full probability distribution by performing the inverse transformation with respect to the p-variable in Eq. (8).In what follows, we first illustrate this rigorously for the analytically tractable cases n = 1 and n = 2.Then, we will combine the intuition gained in these solvable cases with Eq. (11) to derive some general results for arbitrary n.

Linear potential (n = 1)
Let us first look at the n = 1 case for which we can solve the backward Eq. ( 11) analytically.For this case, the Langevin equation ( 6) takes the form In essence, this is a drift-diffusion process in which the particle experiences two distinct drifts.The first one involves a drift with a magnitude k, that tries to pull the particle back to the origin.The second one involves a drift of magnitude v that tends to take the particle away from the origin along the negative x direction.The behaviour of the particle varies depending on which of the two terms dominate.To see this in the context of work done, let us proceed to solve the backward Eq. (11).Depending on the sign of the initial position ξ 0 , we can split the backward equation as where µ = (v + k)/D, γ = (v − k)/D.Solving this set of equations yields To compute the K(p)-functions, we will need four conditions.Two of these are the boundary conditions Q(p, ξ 0 = ±L) which are given in Eq. (12).The other two can be obtained by looking at the behaviour of Q(p, ξ 0 ) and its derivative in the vicinity of ξ 0 = 0. Integrating the backward Eq. ( 11) from −ϵ to +ϵ and then taking ϵ → 0 + , we obtain Using these conditions, it is possible to compute all K(p)-functions as , and inserting them in Eq. ( 16), we get for ξ 0 = 0 This gives us the exact form of the moment-generating function for the case of linear potential.By taking derivative of Eq. ( 19) suitably with respect to p, we can now obtain all moments of the work W .For instance, we find the mean as For all panels, the red curve corresponds to the exact result in Eq. ( 20) and the green curve is the asymptotic expression in Eq. ( 21).We have chosen D = 1.
⟨W ⟩ = kv γDµ Here, we have introduced the simplified notation ⟨W m ⟩ instead of ⟨W m (0)⟩.To gain some physical insights, let us look at the asymptotic behaviour of this mean as L becomes large.This behaviour turns out to crucially depend on the signature of the parameter γ and we find Physically, for γ < 0 (or equivalently v < k), the drift v in Eq. ( 6) is not strong enough for the particle to overcome the force k.Thus, the particle typically takes exponentially large time to reach the absorbing walls which, in turn, gives rise to exponentially large value of the work.On the other hand, for γ > 0 (or equivalently v > k), the particle can easily overcome the constant force and reach the absorbing walls with relatively high probability.Therefore, we get smaller values of the work which scales linearly with L.
For the marginal case γ = 0, we anticipate the scaling to lie somewhere in between γ > 0 and γ < 0 cases.In Figure 2, we have compared these asymptotic results and their exact expressions in Eq. ( 20) with the numerical simulation of the Langevin equation (1).
After analysing the mean, we now proceed to compute the full probability distribution of the work.For this, we need to perform the inverse Laplace transformation of Q(p, 0) in Eq. (19).For large L, we saw that the work typically takes a large positive value for all values of γ.In terms of p-variable, this translates to taking small p limit of the moment generating function.Taking p → 0 limit in Eq. ( 19), we find that Q(p, 0) also takes different forms depending on the signature of γ.In Appendix A, we have  23) and (25).We have chosen D = 1.
rigorously derived these forms for large L and found where α = v/ √ D and β = α e µL − 1 .For γ ̸ = 0, we can perform the inverse Laplace transformation explicitly and obtain the distribution of work to be On the other hand, for γ = 0, we perform the expansion and then carry out the inversion with respect to p to obtain In Figure 3, we have compared our analytic results with the same obtained from numerical simulations.We observe excellent match between the derived results and numerics for γ ̸ = 0. Contrarily for γ = 0, we see a departure at smaller values of W [see middle panel in Figure 3].This stems from the fact that we have truncated the series in Eq. ( 24) at p 2 and ignored the higher order terms.While this is valid for large values of W (which in the Laplace domain corresponds to small p), it becomes less accurate for smaller values of W (which corresponds to large p).Thus, for a more accurate match at small W , one needs to consider higher order terms in Eq. (24).Furthermore, the distribution of W for these different cases turn out to be completely different.For instance, as seen in panel (a) in Figure 3 for γ > 0, the distribution close to the mean can be effectively approximated by a Gaussian form.This can also be seen in Eq. ( 23) where we can replace W ≃ ⟨W ⟩ in the vicinity of the mean and the distribution then simply becomes a Gaussian distribution.However, as we move towards the tail, this approximation ceases to remain valid and one needs to consider the full non-Gaussian form of the distribution.Contrarily, for γ ≤ 0, the distribution is strictly (always) non-Gaussian with exponentially decaying tails.In fact, as illustrated in Figure 3, the distributions for γ ≤ 0 become highly skewed compared to the γ > 0 case.
One can also use the simplified expressions of Q(p, 0) in Eq. ( 24) to calculate the higher moments of work.For example, using Eq. ( 13), we obtain the second moment of W to be Once again, we see the emergence of different large-L behaviours for different signatures of γ, stemming essentially from the same underlying physical reasoning as discussed in the context of the mean.

Harmonic potential (n = 2)
We now consider the other solvable case of harmonic potential.For this case, the backward equation ( 11) takes the form Solving this, we obtain where H m (x) stands for the generalised Hermite polynomial with degree m (where m can be a real number) and 1 F 1 (m; 1/2; x) is the Kummer confluent hypergeometric function.
The other two functions K 1 (p) and K 2 (p) can be evaluated using the absorbing boundary In both panels, green curves are obtained by numerically differentiating the moment-generating function in Eq. ( 28), while red curves are the asymptotic results in Eqs. ( 31) and (32).We have conditions in Eq. ( 12) and we find Now by taking the derivative of Q(p, ξ 0 ), we can find all moments of work done.In Figure 4, we have illustrated this for the mean and the second moment of W by numerically carrying out the derivative of Eq. (28).As seen from this plot, our results are consistent with the numerical simulations.Later, we will derive these moments exactly by using a slightly different (but related) method and show that for large L while the second moment of W behaves as We have verified these scaling behaviours in Figure 4 along with their exact counterparts.
Having looked at the first two moments, we next proceed to calculate the distribution of W . Obtaining the distribution using Q(p, ξ 0 ) in Eq. ( 28) analytically seems daunting.Nevertheless, in Appendix B we have provided a heuristic calculation that correctly gives the distribution for typical values of W at large L. In particular, we find as L becomes large.We have compared this expression with the numerical simulations in Figure 5 and we observe an excellent match between them.For large W , the distribution decays exponentially as ∼ exp (−W/ζ W ) with the decay constant . This is in stark contrast with the case of the fixed observation time, where distribution of W for n = 2 case just turns out to be a Gaussian function [69,72].In fact, the exponential decay turns out to be the hallmark property for the work distributions of all potentials with n > 1 under a first-passage time as we discuss later.

General n
Having examined the analytically tractable cases, we now shift our focus to deriving general results applicable to arbitrary n (> 0).While solving the backward equation (11) directly for a general n is a challenging task, we can utilize it to derive a differential equation governing the moments of W , which we can subsequently solve.These moments can then be employed to deduce the behaviour of the distribution of W for large L as done for the linear and harmonic cases previously.Writing the m-th moment of W as ⟨W m (ξ 0 )⟩ = W m (ξ 0 ), we take the derivative with respect to p in Eq. ( 11) to obtain where W 0 (ξ 0 ) = 1 and f n (ξ 0 ) = sgn(ξ 0 ) ξ 0 n−1 .For m = 1, the solution of this equation is given by and the constants A and B can be computed using the condition W 1 (ξ 0 = ±L) = 0.For ξ 0 = 0, the solution in Eq. ( 35) reduces to This is an exact expression of the mean valid for all n.Although analytically integrating in this expression turns out to be challenging in general, one can simplify the expression and obtain the leading order behaviour of ⟨W ⟩ in the large-L limit.As shown in Appendix C, this turns out to depend on the exponent n and we get qualitatively different results depending on whether n is greater or smaller than 1.In particular, we have shown in Appendix C that For the marginal case n = 1, these results are shown in Eq. ( 21).This can also be obtained from Eq. (37) where n → 1 − gives the γ > 0 case in Eq. ( 21) while n → 1 + leads to the γ < 0 case.Next let us try to understand the difference between n > 1 and n < 1 cases physically.When the potential is sufficiently confining (n > 1 case), the motion of the particle is heavily restricted due to the confining potential [see Figure 1].Thus, typically the particle will be located away from the absorbing walls at ±L and it will hit them very rarely.Therefore, for n > 1, the first-passage time to reach the walls is typically high which in turn gives large value for the work.On the other hand, for n < 1, while the potential is still confining, the corresponding force decays with increasing distance.Hence, in Eq. ( 6), the drift v term is always dominating beyond a certain distance.As a result, we get relatively smaller values of the first-passage time and hence the work.
In Figure 6, we have compared the exact result in Eq. ( 36) with the same obtained from the numerical simulations for n < 1 (left panel) and n > 1 (right panel).We observe an excellent agreement between the theory and the numerics for both cases.In addition to this, we have also presented a comparison of the large-L expressions in Eq. ( 37) for the two cases.We find that Eq. ( 37) converges to the exact result as we take larger values of L.
So far in this section, we have solved Eq. ( 34) only for the mean (m = 1).But one can also utilize it to obtain the higher moments of W .A similar analysis for m = 2 gives with the asymptotic behaviour derived in Appendix D to be Note that these expressions only provide the leading order behaviour in L and there are still sub-leading corrections that depend on model parameters v and k.Interestingly, it turns out that one can use the first two moments to heuristically calculate the probability distribution describing the typical fluctuations of the work [as done before for n = 2].We refer to Appendix B for this calculation where we show We emphasize that this expression only describes the typical fluctuations of W around the mean for n > 1 and will not capture the regimes far away from the mean.In Figure 7 (left panel), we have compared this with the numerics for n = 1.5.We notice an excellent agreement between them.Next, we address the case where n < 1.As observed in the instance of n = 1 with γ > 0, one needs to solve the complete backward equation (11) in order to obtain Q(p, 0).While achieving this for n = 1 is feasible, doing the same analytically for arbitrary n < 1 turns out to be difficult.However, based on extensive numerical simulation, we find that for n < 1 also, the distribution has the same form (around the mean) as for n = 1 with γ > 0 in Eq. ( 23).This leads us to assume for large L.Here A 1 and A 2 and A 3 are constants that depend on the parameters of the model.To compute them, we use the normalisation condition and the fact that the first two moments are given.We then obtain Plugging these constants in Eq. (41) gives The right panel in Figure 7 shows a comparison of Eq. ( 43) with the same obtained from the numerical simulations for n = 0.75.An excellent match between them validates our results.In summary, we have derived the exact expressions for the first two moments and their large-L forms in this section, applicable for all values of n.We then combined these expressions with a heuristic analysis to obtain the probability distribution of W for large values of L across all n values.
Before ending this section, it is worth noting that although we have focused on the large-L behaviour of the work done, one can also use the exact expressions given in Eqs.(36) and (38) for computing the first two moments for small values of L.
To calculate this, we first observe that for small L, one can perform the expansion G ± (y) ≃ 1 + 1 D ky n n ± vy inside integrals in Eqs.(36) and (38).Now the integrations can be carried out analytically and we obtain Compared to the large-L expressions, we find that both ⟨W ⟩ and ⟨W 2 ⟩ exhibit algebraic scaling with L for smaller values of L.

Relation with the first-passage time
So far, we have studied the statistics of the work done till the particle reaches one of the walls at L ± (t) for the first time.Deploying a backward approach in terms of the co-moving variable ξ(t), we have been able to compute the first two moments and the probability distribution of W for large L. In the remaining part of our paper, we show that the fluctuations of the work done W bear an interesting relation with the fluctuations of the first-passage time t f defined in Eq. ( 4).In order to see this, we first write a backward equation for the moments T m (ξ 0 ) = t m f (ξ 0 ) of the first-passage time as done before in Eq. ( 34) for the work done [50,51] where T 0 (ξ 0 ) = 1 and f n (ξ 0 ) = sgn(ξ 0 ) ξ 0 n−1 .Meanwhile the moments satisfy the boundary conditions T m (ξ 0 = ±L) = 0, since the particle gets instantly absorbed if its initial position coincides with one of the absorbing walls.Solving Eq. ( 45) in the same way as done in Section 5, we obtain for m = 1 with functions G ± (y) given in Eq. (35).Constants B 1 and B 2 can be evaluated using the boundary conditions T 1 (ξ 0 = ±L) = 0.For ξ 0 = 0, the mean ⟨t f ⟩ = T 1 (0) is given by For large L, the integrations in this expression can be approximately carried out as done for the work done.Proceeding similarly in Appendix C, we find while for the marginal case of n = 1, we obtain   54) for three cases of n > 1, n = 1 and n < 1.For every n, we have considered scenarios with v < k (pink), v = k (blue) and v > k (green).
Both Eqs. ( 48) and (49) give the leading order behaviours in L. Interestingly, for n > 1, we see that the large-L scaling of ⟨t f ⟩ is similar to ⟨W ⟩ in Eq. (37).We believe this occurs because for large t f , the rate of work done W/t f approaches a constant value (by the law of large numbers), and hence W ∼ t f .Therefore, the stochasticity of the total work is dominated by the stochasticity of the first-passage time.
Next, one can solve Eq. ( 45) for for m = 2 to show that the second moment For large L, one can again simplify this expression [see Appendix D] and obtain and for n = 1, one gets Once again in Eq. ( 51), we only obtain the leading order asymptotic expression in L.
However, there will be sub-leading corrections that depend on parameters v and k and we do not delve into these corrections in Eq. ( 51).
Next, we set out to examine an interesting connection between the first-passage time and the work done.To this end, we define their respective coefficents of variation or the variabilities as and then study the following ratio As we show below, this "variability ratio" has intricate behaviors as a function of system size L and the potential order n.
For n = 1, one can plug the analytic expressions of first and second moments of W and t f derived in the previous sections and show that the ratio R(L) converges to the value 1 as the length L becomes large for all values of v and k.Similarly, for n > 1, one can use the asymptotic results ⟨W 2 ⟩ ≃ 2 ⟨W ⟩ 2 and t 2 f ≃ 2 ⟨t f ⟩ 2 to show that both CV W and CV t f individually go to 1.This indicates that the ratio R(L) also converges to 1 for larger L for n > 1. However for n < 1, an analytical calculation of R(L) remains unlikely.As indicated in Eqs. ( 39) and ( 51), for the case of n < 1, one needs to have the knowledge of the sub-leading corrections in L within the expressions of ⟨W 2 ⟩ and ⟨t 2 f ⟩ in order to compute the CV -s.Moreover, extracting these corrections analytically from the exact expression turns out to be difficult and our analysis in Eqs. ( 39) and ( 51) cannot capture these details.
In Figure 8, we have plotted the ratio R(L) as a function of the initial wall separation L for different values of the exponent n.For n ≥ 1, we find that R(L) indeed saturates to the universal value equal to one at large L independent of the model parameters k, v and D. On the other hand, for n < 1, the large-L behaviour of R(L) depends on these parameters.For instance, in Figure 8 [right panel], we see that the saturation values for v = 3, k = 3 (shown in blue) and v = 3, k = 2 (shown in green) are completely different for n = 0.8.On the other hand, for v = 2, k = 3 (shown in pink), R(L) increases monotonically at large L and does not saturate to a constant value.This is possibly due to the fact that the saturation of R(L) only takes place at large L, i.e. for L ≫ L th , where L th is some threshold length.For a given n, the threshold L th also depends on the parameters v and k.In Figure 8 [right panel], we find that L th ∼ 50 for the green curve (k/v < 1) whereas L th ∼ 500 for the blue one (k/v = 1).Thus, we see a consistent rise in the value of L th as the ratio k/v starts to increase.Following this trend, we anticipate L th to be very large for the pink curve (k/v > 1).Probing such large L values in simulations is computationally expensive.Thus, we observe that R(L) increases for the pink curve since we are still in the regime L ≪ L th .However, if we go to very large L values (beyond what is considered in the figure), we anticipate R(L) to converge to a constant value even for the pink curve (v = 2, k = 3).

Conclusion
In summary, we have analyzed the work statistics in a driven overdamped system in one dimension.In contrast to the conventional cases where the observables are usually measured upto a fixed time, here we measure the work upto a random first passage time that is conditioned on a certain criterion.We employ the Feynman-Kac path integral approach to compute the work functional in various set-ups consisting of potentials of different configurations.We provide a comprehensive analysis of the work fluctuations which shows a rich behaviour as a function of the potential strength, the external drive and the interval over which the particle moves.Furthermore, we showed an interesting symmetry relation between the signal-to-noise relation or the variability of the work and that of the first-passage time.To understand this phenomena better, we also delved deeper into the attributing physical conditions.Notably, our results illustrate a marked difference in the work statistics between systems driven up to a fixed time and a firstpassage time.Although we considered a fixed initial condition x 0 = 0 to calculate these results, our method can, in principle, be extended to a general distribution P (x 0 ).In particular, we anticipate the large-L expressions in Eqs.(37) and (39) for work functionals and in Eqs.(48) and (51) for first-passage time to be valid even for a general P (x 0 ) as long as the variance of x 0 does not scale with the initial separation L.
Moving forward, it would be interesting to study the "variability ratio" in other systems and to identify any similar pattern as unveiled here.There has been a myriad of studies in recent times in the field of stochastic thermodynamics involving various other thermodynamic observables such as injected power, dissipated heat or entropy production.It is only natural to study the same also within our set-up and will be pursued elsewhere.Furthermore, our model could also be useful to study existing thermodynamic bounds for first-passage-time problems [24,37,41,73,74].Another interesting direction is to investigate how our method can be extended to incorporate other potentials.For instance, we assumed that both the potential and the boundaries move at the same velocity v.It would be interesting to explore what happens if we relax this condition and allow them to move with different velocities.Concluding, we believe that our results can be tested using controlled optical trap experiments [72,75,76] which, hopefully, will also unfold new research directions to this problem.
from the vicinity of ω = 1.Performing an expansion of the integrand in Eq. (C.4) around y = 1 − u and keeping the leading order in u, we get the large-L behaviour of J ± (L) as We next turn to I ± (L).Rewriting their forms from Eq. (C.2) Let us first analyse the integration over z for large values of L. By changing the variable u = kz n /nD, this integral can be rewritten as ˆωL The sum converges for n > 1.This implies that the integration over z in Eq. (C.8) becomes independent of ωL as L becomes very large (for a given ω) and I ± (L) accordingly takes the form Let us now derive the large-L expression of the mean work ⟨W ⟩ for n < 1.Looking at Eq. (C.1), this again translates to calculating the functions J ± (L) and I ± (L).We first notice that even for n < 1, G + (y) has a faster than exponential growth for large y.Therefore, one can follow a similar approach as described in Section Appendix C.1 for J + (L) and I + (L) and obtain On the other hand, G − (y) for n < 1 is a decaying function for large y and requires a different approach.To this end, we recast J − (L) as and carry out the integration over y for L → ∞ to yield For n < 1, this is a convergent sum.Hence J − (L) becomes independent of L for this case.
Next we turn to I − (L) whose explicit expression follows from Eq. (C.2) as By changing the variables y = ωL and z = uωL, this expression becomes For large L and fixed ω and u, one can simplify this expression as with G ± (y) defined in Eq. (35).In order to find the large-L expression of W 2 , we have to calculate the corresponding form of these functions.Meanwhile, we computed J ± (L) in Appendix C and showed that they take different forms depending on whether the exponent n is greater than or smaller than one.In the following, we carry out the same analysis for Y ± (L).Looking at this, we first notice that for y integration, the integrand has a faster than exponential growth.This implies that at large L, this integration will be dominated predominantly by the larger values of y.On the other hand, the integration over z is damped due to the presence of the decaying exponential term.Therefore, even though y takes larger values, the integration over z will still be dominated by small values of z (more precisely z ≪ L).This allows us to replace W 1 (±z) ≃ W 1 (0) in Eq. (D.In the remaining part of this section, we provide a derivation of Eq. (D.11) which was instrumental in deriving W 2 for n < 1 case.Following Eq. ( 35), we can write the mean work as

Figure 1 .
Figure 1.Left and right panels: The depiction of the potential U (x, t) (via Eq. (2)) and the corresponding force F (x, t) for different values of n.Parameters: k = v = t = 1.

Figure 4 .
Figure 4. First two moments of the work done by the particle in presence of the harmonic potential (n = 2).In both panels, green curves are obtained by numerically differentiating the moment-generating function in Eq. (28), while red curves are the asymptotic results in Eqs.(31) and(32).We have chosen v = k = D = 1.

Figure 5 .
Figure 5. Probability distribution of the work done W for the time-dependent harmonic potential (n = 2).Analytic expression in Eq. (33) is plotted in black and the simulation data are shown in red.Inset shows the comparison for small values of W . Parameters chosen are v = k = D = 1, L = 4.

Figure 6 .
Figure 6.Illustration of the average work ⟨W ⟩ performed by the particle in potentials with n values 0.75 and 1.5.The theoretical plot is based on the expression in Eq. (36), presenting the exact calculation.In addition, the asymptotic large-L results are given in Eq.(37).Parameters chosen are k = 0.5, v = D = 1.

Figure 7 .
Figure 7.Comparison of the probability distribution of W outlined in Eqs.(40) and (43) is presented in the left panel for n > 1 and L = 10 and the right panel for n < 1 and L = 100.We have set the parameters as k = 0.5, v = D = 1 for this comparison using numerical simulations.

Figure 8 .
Figure 8. Simulation results for the ratio R(L) defined in Eq. (54) for three cases of n > 1, n = 1 and n < 1.For every n, we have considered scenarios with v < k (pink), v = k (blue) and v > k (green).

C. 13 )
where the second line follows from Eq. (C.6).Combining this result with the expressions of J ± (L) in Eq. (C.6) and plugging all of them in Eq. (C.1), we find⟨W ⟩ ∼ L 1−n exp 1 D kL n n − vL ,for n > 1 and large L. (C.14) Appendix C.2. n < 1 case