Asymptotics of work distributions in a stochastically driven system

We determine the asymptotic forms of work distributions at arbitrary times $T$, in a class of driven stochastic systems using a theory developed by Engel and Nickelsen (EN theory) (arXiv:1102.4505v1 [cond-mat.stat-mech]), which is based on the contraction principle of large deviation theory. In this paper, we extend the theory, previously applied in the context of deterministically driven systems, to a model in which the driving is stochastic. The models we study are described by overdamped Langevin equations and the work distributions in the path integral form, are characterised by having quadratic actions. We first illustrate EN theory, for a deterministically driven system - the breathing parabola model, and show that within its framework, the Crooks flucutation theorem manifests itself as a reflection symmetry property of a certain characteristic polynomial function. We then extend our analysis to a stochastically driven system, studied in ( arXiv:1212.0704v2 [cond-mat.stat-mech], arXiv:1402.5777v1 [cond-mat.stat-mech]) using a moment-generating-function method, for both equilibrium and non - equilibrium steady state initial distributions. In both cases we obtain new analytic solutions for the asymptotic forms of (dissipated) work distributions at arbitrary $T$. For dissipated work in the steady state, we compare the large $T$ asymptotic behaviour of our solution to that already obtained in ( arXiv:1402.5777v1 [cond-mat.stat-mech]). In all cases, special emphasis is placed on the computation of the pre-exponential factor and the results show excellent agreement with the numerical simulations. Our solutions are exact in the low noise limit.


Introduction
Stochastic thermodynamics extends the definition of thermodynamic quantities such as entropy, heat and work, to the level of stochastic trajectories, and has become an important area of research in non-equilibrium statistical mechanics of small systems [5]. Non-equilibrium fluctuations are particularly relevant for systems with a small number of degrees of freedom, where there are large fluctuations around the average and considerable contributions from rare events. Hence a major advancement in the field was due to the discovery of fluctuation theorems (FTs) which extend equilibrium fluctuation-dissipation relations to far-from-equilibrium regimes. Depending on the initial conditions of the system under consideration, FTs put constraints on various thermodynamic distributions and demonstrate that the second law holds statistically in systems with stochastic dynamics; a positive entropy production is exponentially more likely than its negative counterpart. Among the various versions of the FT, the notable ones are the Crooks Fluctuation Theorem (CFT) [6], and the Jarzynski Equality (JE) [7,8] which relate non-equilibrium quantities such as the thermodynamic (Jarzynski) work to equilibrium quantities such as the free energy difference. Advancements in technology have allowed experimental verifications of these fluctuation theorems in a variety of systems (see references in [5]).
Relatedly, there has been significant interest in calculating work distributions for different models. A well studied example is that of a colloidal particle in a harmonic optical trap, where either the mean position of the trap [9,10,11,12] or the stiffness [13,1,14,15,16] is externally modulated. In the literature these models are referred to as the sliding parabola or the breathing parabola respectively. For the sliding parabola, if the driving is deterministic, the work distribution is known to be a Gaussian [9]. For the breathing parabola, the solution given in [13] is formally exact for an arbitrary driving protocol. In [17], an exact calculation of the work distribution has been carried out for the Brownian particle in a logarithmic-harmonic potential. In [18], exact work statistics have been obtained for a Brownian particle in the presence of non-conservative forces such as torques. Table 1: Asymptotic behaviour of work distributions. List of asymptotic forms of P (W ) known so far adapted from [22]. A, B, C are constants that depend upon the explicit form of the driving protocol and the duration of the protocol T .
It is however hard to find an exact expression for the full work distribution P (W ) except in the few cases mentioned above. This problem is resolved to a certain extent by using techniques based on large deviation theory [19]. One of the recent developments in this regard is a theory developed in [1] by Engel and Nickelsen (EN Theory), which is used to compute the tail forms of work distributions including the pre-exponential factor at arbitrary times T . To derive the asymptotic probability for a certain rare work value, the probability of the most likely trajectory that gives rise to this specific work value is considered. The constraint to rare work values is effectively equivalent to a low temperature limit. The advantage of this method over other similar large deviation techniques [20,3] is that the results for the tails are also valid for very small time durations of the process which are typical for experimental situations. EN theory also has the additional benefit that it reduces the problem of computing the asymptotic form of the work distributions to solving a system of ordinary differential equations. One can therefore use available BVP solvers [21] to obtain the full asymptotic form [1,22]. In addition to the problems studied in [1], a few more systems have been studies using EN theory [23,24]. In [22], Holubec et al put forward a functional-form conjecture to classify the exact asymptotic form of P (W ) based on the form of the driving used. In addition, they also obtain the analytic solution to the asymptotic form of P (W ) in an absolute value potential (V-potential). Table 1 is adapted from [22], where the so far known results for the asymptotic forms of work distributions for various driving protocols are listed. In all the above cases the driving protocol λ(t) considered, is a deterministic function of time.
In this paper, we will first revisit a model with a deterministic driving protocol, the breathing parabola [1], to familiarize the reader with the methods of EN theory. Then, in this simple model, by considering a particular forward and reverse protocol, we will illustrate the fluctuation theorem for the Jarzynski work. We will show that, within the framework of the EN theory, the fluctuation theorem manifests itself as a symmetry property of a certain characteristic polynomial, which also determines the exact moment generating function at finite times. We then extend this analysis to a stochastically driven system (i.e. λ(t) is stochastic). Stochastic driving protocols have been looked at previously in [20,3,4] and recently in [25]. In [20], λ(t) is taken to be a Gaussian random noise, and the large-time asymptotic form of the distribution of work done by the stochastic force λ(t) (the product of the stochastic force λ times the displacement dx) has been computed, and fluctuation theorems analysed. In [3], λ(t) is considered instead to be the Ornstein Uhlenbeck process. This is also the stochastic protocol that we study in this paper. The model describes the dynamics of a colloidal particle in a harmonic potential whose mean position is stochastically modulated; we call this model the stochastic sliding parabola (SSP) model. The model studied in [25] in the context of the recently discussed finite time thermodynamics uncertainty relation, is equivalent to a discrete version of the SSP. In [3], it was claimed that the work fluctuation theorem was valid only in certain regions of the parameter space. Subsequently it was shown in [4], that when the dissipated work, W d ( defined as the difference between the Jarzynski work, that differs from the definition in [20,3] by a boundary term, and the equilibrium free energy difference ∆F ) is considered, the Crooks fluctuation theorem is satisfied in all regions of the parameter space. In all these cases, the authors used a moment generating function method, and the P (W ) obtained is valid only in the large T limit. In this paper, we study the SSP for both equilibrium and non-equilibrium initial conditions using EN theory. We show that when the initial points are sampled from an equilibrium distribution, the Jarzynski work (W ) satisfies a transient fluctuation theorem. To our knowledge, this has not been noted before in the context of this model. We compute the exact forms of the tails of the PDF for different time durations, and compare them with numerical simulations. For the special case when λ(t = 0) is unconstrained, we are able to obtain the closed asymptotic form of P (W ) as a function of T . For non-equilibrium steady state initial conditions, we obtain the asymptotic form of P (W d ) in a similar manner and compare its large-time behaviour with the asymptotic form obtained using the results in [4]. Our comparison shows that, for very large times, the leadingorder asymptotic form of both methods match but they disagree in the sub-leading pre-exponential behaviour. We validate all our results using numerical simulations. In all cases, special emphasis is put on the computation of the pre-exponential factor, which involves calculating a fluctuation determinant [26,27] that is a functional determinant ( of a matrix differential operator [28,29] in the case of the SSP). We do this by using a technique developed in [28], which is based on the spectral -ζ functions of Sturm-Liouville type operators.
The paper is organized as follows. In Section 2, we introduce various techniques and methods used in this article and illustrate EN theory in the context of the breathing parabola problem. Particularly, in Section 2.2, we discuss the FT for a specific choice of forward and reverse protocols and obtain the asymptotic forms of P (W ) in each case. In Section 3 we extend this analysis to the SSP model. Exact asymptotic forms of PDFs of the Jarzynski work (W ) and the dissipated work (W d ) are obtained for equilibrium and non-equilibrium steady state initial conditions in Sections 3.1 and 3.2 respectively. In Section 3.2.2, we compare the asymptotic forms of P (W d ) with the results from [4]. In Section 4 we present our conclusions. Various technical details of the paper including the computation of the pre-exponential factor, are presented in Appendices A -E.

Basic methods
For systems modelled using overdamped Langevin equations, EN theory may be used to compute the asymptotic form of the work probability distribution analytically [1,30]. The theory was initially developed for equilibrium initial conditions, and later generalized to non-equilibrium initial conditions as well [24]. The method can be summarized as follows: 1. Let W [x(·)] be any functional of the stochastic trajectory x(·). Write down P (W [x(·)] = W ) as a path integral, with the corresponding action S[x(·)], by constraining the trajectories to have a work value W .
2. The first order approximation of P (W ) for rare W is obtained as whereS is the action S evaluated along the optimal trajectory that minimizes the action, found by solving the corresponding Euler Lagrange equations together with natural boundary conditions.
3. An improved estimate is then obtained by including the pre-exponential factor, which also takes into account the contributions from the trajectories lying close to the optimal trajectory. This is done by expanding S to second order in variations around the optimal trajectory and by performing Gaussian integrations over the variations.
As is usual for Gaussian integrals, the computation of the pre-exponential factor requires the evaluation of the ratio of determinants of certain differential operators (functional determinants). In this article we use the generalization of a method developed in [28], based on the spectral -ζ function of Sturm-Liouville operators to determine this ratio. In section 2.1, we illustrate EN theory using the breathing parabola model. This model has been studied by Engel and Nickelsen [1] and the full asymptotic form of the work distribution including the pre-exponential factor has been computed. This model serves as a good starting point to the discussions that follow in the paper since the techniques we will use later are the generalizations of the techniques presented here. We will stick to the notations used in [1] unless required otherwise.

Illustrative example: The breathing parabola
The breathing parabola potential is given by where the stiffness of the trap λ(t) is the driving protocol which is varied deterministically from a value λ(t = 0) = λ 0 to λ(t = T ) = λ T during each realization of the process. The motion of a colloidal particle in this potential can be described by the overdamped Langevin equation, where V (x) = dV /dx and β is the inverse temperature defined as 1/(k B T ). η(t) is a Gaussian white noise, with η(t) = 0, and η(t) η(s) = δ(t − s). The Jarzynski work done along a trajectory x(·) of this system for a time interval [0, T ] is defined [31,32] as, W [x(·)] being a functional of a stochastic trajectory, is a random variable by itself. Using the Onsager Machlup formalism [33,34], the probability of W [x(·)] taking a value W can be written down as a path integral [1], where the augmented action S is given by, N is the normalization constant corresponding to mid point discretization in the functional integral and Z 0 is the initial equilibrium partition function. EN theory uses the contraction principle of large deviation theory [19] to calculate the asymptotic behaviour of P (W ) for large W . This is implemented by first evaluating the integrals in Eq. (6) using the saddle point approximation, and finding the optimal trajectory (x(·),q) which minimizes S. To find the optimal trajectory, S is studied near the vicinity of a trajectoryx(t) and a valueq of q by writing x(t) =x(t) + y(t) and q =q + r and by expanding S to second order 1 in y(·) and r as, where A is a second order Sturm-Liouville type differential operator,

Leading-order behaviour
The leading-order asymptotic behaviour of P (W ) is obtained by approximating, whereS is the action evaluated along the optimal trajectoryx(·) andq, obtained by demanding that S lin vanishes for an arbitrary variation (y(·), r). This results in the Euler-Lagrange equations, with boundary conditions,ẋ The above equations constitute a Sturm-Liouville eigenvalue problem with the parameter iq. Therefore there are infinitely many values iq (n) for which a non trivial solution exists, all of which are saddle points of the action S. It can be verified that each iq (n) value corresponds to two solutions, ±x n (t). It is possible to show that all saddle points except the one corresponding to the smallest iq(≡ iq * ) value are unstable because they fail the Hessian (second derivative) test 2 . Therefore they do not contribute to the asymptote of P (W ). The term proportional to r in Eq. (8b) gives us the constraint equation, Using Eq. (11), (12) and Eq. (13), it can be shown that along the optimal trajectory, Therefore using Eq. (10), the leading-order asymptotic behaviour of P (W ) is determined to be,

The pre-exponential factor
The next step is to improve this estimate by also taking into account the fluctuations around the optimal trajectory. This is done by retaining S quad (Eq. (8c)) in the exponent of Eq. (5) by computing the Gaussian integrals in As already mentioned, since the action is quadratic, the fluctuation governing operator A that appear in S quad is the same as the operator for the optimal trajectory. The operator A therefore has a zero-mode which is nothing but the optimal trajectory itself. Writing integration variable y as a series expansion in terms of the normalized eigenfunctions φ n (t) of A as, y(·) = Σ n c n φ n (t) (17) and then performing the Gaussian integrals over the expansion parameters c n and r, it can be shown that the zero-mode of A gets omitted naturally and does not cause any problems to the integral in Eq. (16). This gives us a compact expression for the pre-exponential factor, The notation det A iq=iq * stands for a determinant omitting the zero-mode. The factor d 0 is defined as, wherex(t) is the zero-mode and ||x|| stands for its norm. J is a factor stemming from the Jacobian of the transformation of integration variables. The value of J for the breathing parabola problem and also a class of similar potentials has been determined in [1], and can be shown to be equal to, A derivation of Eq. (18) and (19) can be found in [1]. Now using Eqns. (20), (18), (15) and (5) we finally obtain, The factor 2 corresponds to the two equipotent saddle points (iq * , ±x(t)). Computing the pre-exponential factor now reduces to evaluating d 0 and the ratio of functional determinants appearing in the expression above. As we will see, computing d 0 is rather straightforward. The evaluation of the determinant ratio is however more involved and is carried out using a technique developed in [28]. Notice also that in Eq. (21), the protocol λ(t) is not specified. In the next Section, we apply the method to a specific forward and the corresponding reverse protocol and exactly compute the asymptotic form of P (W ) including the pre-exponential factor, for both cases. We also analyse the FT within the framework of EN theory.

Forward and reverse protocols: Fluctuation relation
In this Section, we will compute the exact asymptotic form of the work distribution in the breathing parabola problem, for a specific choice of forward (F) and reverse (R) protocols. The work distributions are then known to satisfy the Crooks fluctuation theorem: We consider the following forward protocol, and the corresponding reverse protocol, This particular reverse protocol was studied as an illustrative example in [1]. In order to find the leading-order behaviour of P (W ) in each case, one has to obtain the smallest iq values for which the ELE Eqns. (11), (12) have a non trivial solution. Here we use a formalism used in [28]. As we will see, this method will turn out to be useful for the computation of the pre-exponential factor as well. First we write the boundary conditions as two matrix equations (for notational simplicity, we will use x(t) instead ofx(t).), In this case, M and N are matrices, It can then be shown that the iq values for which the equations (11), (12) have a non trivial solution, can be obtained as the roots of a function (the characteristic polynomial), where H k (t) is the matrix of fundamental solutions of ELE Eqns. (11), (12) defined as [28,36], We will also make a particular choice of H k (t), namely that H k (0) = I 2 . For convenience, we have defined a variable, The iq * value may then be found by looking at the smallest roots of the characteristic polynomial F (k). In Fig.  1 we plot F (k) vs. k for both the forward and reverse protocols for a time period T = 1. Notice the interesting symmetry of the characteristic polynomials namely that, F (k) for the forward protocol is the mirror image of F (k) of the corresponding reverse protocol. As we will show in Appendix D, this symmetry is a consequence of the FT [6,16,37]. As a consequence of this form of F (k), we have For T = 1, we get k * F = 3.66 = −k * R and therefore P F/R (W ) has the leading-order asymptotic form (β = 1), Now we move on to the computation of the pre-exponential factor. Obtaining the factor d 2 0 is straightforward, using Eq. (13) in Eq. (19) we see that, where ||x|| is the norm of the zero-mode, and is defined as the usual inner product, The * in the equation above stands for complex conjugation. The computation of the second factor, which is the ratio of the functional determinants, is rather involved, and may be evaluated using a technique that is developed based on the spectral ζ function of Sturm-Liouville type operators [28]. The power of the method is that it enables one to compute the determinant ratio in terms of the zero-mode itself. In terms of suitably normalized zero-mode solutions, the determinant ratio become The subscript N denotes that a particular normalization is chosen. Explicit details of the calculation including the choice of normalization will be discussed in Appendix B and C. Using Eq. (21), Eq. (32) and Eq. (34), the exact asymptotic form of P (W ) with the pre-exponential factor can now be determined for both forward and reverse protocols. Doing explicit computations for T = 1, (see Appendix C) we get, The asymptotic forms obtained above are consistent with the Crooks fluctuation relation (Eq. (22)), with β∆F = 1 2 ln k f ki = 1 2 ln 1 2 for the reverse process [13,16]. We have also calculated P (W ) numerically, and the results are in good agreement with the theoretical predictions.

Remarks
The main points in the calculation of the asymptotic form of P (W ) for the breathing parabola hold in general for PDFs with a class 3 of quadratic augmented action S W (Eq. (6)), and may be summarised as follows.
• P (W ) has, to leading-order, the functional form, where iq * is the smallest root of some characteristic polynomial F (k ≡ 1 − iq). The Crook's fluctuation theorem is seen to manifest itself as a reflection-symmetry property of this characteristic polynomial F around k = 0. This symmetry is no surprise when we realize that the exact moment generating function of work in this case is given by, The iq * value which determines the asymptotics correspond to the singularity of this moment generating function lying close to zero. We will derive Eq. (37) in Appendix D.
• Including the pre-exponential factor, P (W ) takes the form, where all the factors may be computed in terms of the fundamental solutions of the Euler Lagrange equations. This form for the tail is exact in the low noise limit β → ∞, since there are no higher order expansion terms in (7) that we are neglecting.
In the following Sections, we will apply the techniques developed here to a stochastically driven system.

The stochastic sliding parabola
Consider the dynamics of a colloidal particle in a harmonic trap where the mean position of the trap is externally modulated [9,10,30,3,4,20]. Such potentials go by the name sliding parabola, and have the general form : where x(t) is the position variable and λ(t) is the externally modulated mean position. We have set the stiffness of the trap to 1. The dynamics of the colloidal particle in this potential can be described by the Langevin equation,ẋ where η(t) is a thermal noise and D is the diffusion coefficient. η is assumed to be Gaussian with η(t) = 0, and η(t) η(s) = δ(t − s). One of the natural time scales in the system is given by the relaxation time in the harmonic trap τ r = γ/κ where γ is the friction coefficient and κ is the stiffness of the trap. In the case of a deterministic driving protocol, the exact statistics of the Jarzynski work done on the colloidal particle is known [9]. The work distribution is a Gaussian and satisfies the transient fluctuation theorem. The sliding parabola with a deterministic driving was also looked at in [1,30], as a test example for the EN Theory, and in this case the method gives the full probability distribution (not just the large-W form). Eq. (40) has also been studied both experimentally [2] and analytically [3,20,4] when λ(t) is a stochastic driving protocol. One of the cases studied is when λ(t) is the Ornstein-Uhlenbeck process given by, ξ(t) is again assumed to be a Gaussian noise with ξ = 0 and ξ(t)ξ(s) = δ(t − s). The noise ξ is usually athermal in origin with a diffusion coefficient A as given in Eq. (41). τ 0 gives the second natural time scale in the system in terms of the relaxation time of λ correlations. The two noises are assumed to not have cross correlations, i.e. η(t) ξ(s) = 0. Hereafter, we will refer to the coupled equations (40) and (41) as the Stochastic Sliding Parabola (SSP).
In the remaining Sections of this paper, we study the SSP model for both equilibrium and non-equilibrium steady state initial conditions using EN theory. For equilibrium initial conditions, the dissipation function that satisfies a fluctuation theorem of the SSP can be identified with the Jarzynski work 4 . The form of the dissipated work in the steady state has been obtained in [4]. For both situations, the exact form of work distributions at arbitrary times is not known. Here we show that the discussions in Section 2.3 can be applied to this system, and we can hence compute the exact asymptotic form of both the transient and steady state work distributions including the pre-exponential factor, at arbitrary times T . For steady state initial conditions, we compare our results with [4], in the appropriate limits. Without loss of generality, for the calculations that follow, we set A = D = k B T and τ 0 = τ r = 1.

Equilibrium initial condition. : Transient fluctuations
In the first case that we look at, we compute the asymptotic form of the distribution of the Jarzynski work done on the colloidal particle (W ) by the stochastic force Eq. (41) starting from an initial equilibrium distribution given by, The particle is assumed be in thermal equilibrium initially for a fixed value of λ 0 and the partition function Z 0 is computed accordingly. The Jarzynski work done on the colloidal particle along each trajectory is defined in the same way as before : In terms of the joint probability density functional of trajectories {x (·) , λ (·)} T 0 , the probability density function of work can be written down as, with the action The normalization constant for this case is [39], In order to find the large W asymptotic behaviour of P (W ), we will adopt the methods discussed in Section 2.1. Here that means, we need to identify the optimal choice of bothx(t) andλ(t) that minimizes the action S for a given value of W . We follow the same procedure as before and put x(t) =x(t) + y(t), λ(t) =λ(t) + z(t) and q =q + r and expand S to second order in y(·), z(·) and r.
Notice again that since S is quadratic, the expansion will not contain terms of order greater than quadratic in y(·), z(·) and r.

Leading-order form of P (W )
As in the case of the Breathing Parabola problem, The leading-order form of P (W ) can be computed in terms of the optimal trajectory (x(t),λ(t)) as, where (x(t),λ(t)) solve the following Euler Lagrange equations, together with Robin-type boundary conditions,x and the constraint equation, Notice that (as we have seen in section 2.1.1,) a non trivial solution to the ELEs exists only for some specific values of iq ≡ iq, and only the smallest iq value (≡ iq * ) is relevant. Based on the discussion in Section 2.3, we can conclude that, to leading order, The iq * value may be found by looking at the roots of the function, corresponding to this problem (see Appendix C). In Fig. 3 we plot F (k) vs k for T = 1. It may be verified that This form of the Crooks fluctuation theorem is a consequence of the fact that ∆F = 0 for this particular choice of initial conditions. The smallest roots of F (k) are found to be k * = ± 2.300. Solving Eq. (49) with Eq. (50) for k = ± 2.3, and then using the constraint equation (51), it can be verified that k = + 2.3 corresponds to the positive tail of P (W ) and k = − 2.3 corresponds to the negative tail of P (W ). Hence using Eq.(52), the leading-order asymptotic forms of P (W ) may be written down as (For β = 1), Notice that the asymptotic forms are again consistent with the fluctuation theorem (Eq. (54)).

P (W ) including the pre-exponential factor
As in Section 2.1.2, the asymptotic estimate for P (W ) obtained above can be improved by also taking into account the contributions coming from S quad in Eq. (47). Since the action in Eq. (45) is quadratic, it can be shown that the fluctuations around the optimal trajectory are again governed by the same operator A (Eq. (49)) which determines the optimal trajectory. P (W ) including the pre-exponential factor therefore takes the form: The form of the Jacobian that is required to derive Eq. (56) is obtained in Appendix A. As in the previous case, det A in Eq. (56) is the determinant of the operator A omitting the zero-mode. d 0 is given by, where || x λ || is the norm of the zero-mode. The factor 2 in Eq. (56) again accounts for the two equipotent saddle points, (±x(t), ±λ(t), iq * ± ). Notice that one significant difference from the case of the breathing parabola is that the functional operators appearing in the determinant ratio in Eq. (56) are 2D functional operators. We hence generalise the method of functional determinants [28] that we used for the 1D Sturm-Liouville operator in the previous case, for this 2D case, to compute the ratio of the determinants appearing in Eq. (56). The explicit details are given in Appendices B and C. For the case T = 1 one can hence compute, and therefore using Eq. (56), the improved estimate to the positive and negative tails are, In Table 2 we give the exact asymptotic forms of P (W ) for different values of T .  Table 2: Asymptotic Forms of P (W, T ). The exact asymptotic form of P (W ) can be determined for any value of T . Table 2 gives the exactly computed form for a few T values, which may be compared with experiments / numerical simulations.
We have also numerically integrated the Langevin equations to get an estimate for P (W ). For a step size of ∆t = 10 −3 , averaging over 10 6 realizations gives a reasonably good agreement with our theoretical predictions. The results are plotted in Figure 4.  Table 2 using EN theory.
The results agree for large values of W , except for the very far tails where there are not sufficiently many sample points. In Figure 5 we present the numerical result, verifying the fluctuation theorem. The case when we leave λ 0 unconstrained, (equivalent to sampling λ 0 from a Gaussian distribution with very large variance), is a special case in which one can obtain a closed asymptotic form of P (W ) as a function of T -the time period of duration of the protocol. Following similar calculations as for the previous case (Appendix C), but with modified boundary conditions (50), we find for the positive tail, and for the negative tail, For T → ∞ we find that the time-dependent part of the the pre-exponential factor converges to a value 7 2 . Notice that the positive tail of the probability distribution decays as a power-law in W , and therefore the mean and higher moments do not exist. This behaviour of the tails can be attributed to the large fluctuations in the system which contribute to positive work values. However, for any finite variance of the initial Gaussian distribution of λ 0 , the work distributions can be shown to have tails of the form (59), and well defined moments.

The steady state fluctuations.
In this section we study the steady state work fluctuations in the SSP. The functional form of the dissipated work (W d ) in the steady state of the SSP was identified in [4]. It was then shown to satisfy the FT by computing the corresponding moment generating function in the large T limit. Here we look at the exact asymptotic form of P (W d ) using EN theory. We will then make a comparison with results from [4] in the appropriate limits.

Asymptotic form of
The dissipated work in the steady state in the time interval [0, T ] for the SSP is given by [4], The form of the dissipation function may also be identified from the formalism presented in [38]. Sampling the initial points from the stationary probability distribution (from [3]), one can write down the probability distribution for dissipated work in the steady state as, with the augmented action Notice that the action in Eq. (65) is again quadratic. Therefore one can infer that in the asymptotic regime, P (W d ) must again be of the form (Section 2.3): It can be verified that the functional operator that determines the optimal trajectory, as well as the fluctuations around the optimal trajectory, is again A, as in the previous Sections 2.1.2 and 3.1.2. The difference comes in the boundary conditions; matrices M and N get modified accordingly. As before, the leading-order behaviour of P (W d ) can be found from the smallest roots of the corresponding characteristic polynomial F (k) ( Eq. (53) ). In Fig. 6 is a plot of the characteristic polynomial for T = 1. Notice that F (k) is again symmetric about k = 0, and this indicates that the dissipated work given by Eq. (62) indeed satisfies the fluctuation relation The smallest roots of F (k), are found to be k * = ±3. 16 and therefore the leading-order asymptotic form is given by (for β = 1): where the ± denotes the positive and negative tails respectively. As done earlier, one can also improve this estimate by including the pre-exponential factor. The steps involved are identical to the previous cases considered. The quadratic nature of the augmented action in Eq. (65) leads to a sub-leading power-law behaviour ∼ |W d | − 1 2 , and a numerical factor which is completely determined in terms of the zero-mode of the operator A, and depends only on the time duration T of the protocol. Explicit calculations may be carried out in the same way as in the previous case (see Appendix C). For T = 1 we find: In  Table 3: Asymptotic Forms for P (W d , T ). The exact asymptotic forms of P (W d ) for different values of T , obtained using Eq. (66).

Comparison with results in [4].
Verley et al [4], obtain the exact form of the generating function of the probability distribution of the dissipation function for large values of the time duration T of the protocol, and show that the Crooks fluctuation theorem is satisfied in this limit [4]. In order to compare the two methods, we have first inverted the generating function from [4] (details in Appendix E). The exact form of P (W d , T ) for T 1 is obtained as : In order to compare this result with our calculations, we do an asymptotic expansion of Eq. (70) for large W d . To leading-order we find that, We compare this with the leading-order form in Eq. (66) by checking how k * (T ) behaves as T becomes large. Knowing the exact form F (k) for any value of T , this behaviour may be readily found. As we show in Fig. 7, we find, and this leads to the asymptotic forms given in Eq. (71). Therefore in the large-T limit, the leading-order behaviour predicted by both methods agree. Now we look at the sub-leading pre-exponential behaviour predicted by both methods. The asymptotic expansion of the pre-exponential factor of P (W d , T ) in Eq. (70) gives a sub-leading pre-exponential behaviour ∼ |W d | −5/2 , which is a much faster decay than predicted by Eq. (66) which suggests a sub-leading behaviour ∼ |W d | −1/2 for any value of T . In Figure 8, we compare the results from both methods with our simulations of the Langevin dynamics. The asymptotic forms obtained in Table 3.2.1 are good fits to the tails for all values of T . As T gets larger, Eq. (70) becomes a good fit to the numerical data, and the two methods agree to a good extent at the tails (as expected from the same leading-order behaviour). In order to see the difference coming from the disagreement in the pre-exponential factors of both methods, extensive Langevin simulations will be needed.

Conclusion.
In this paper we have determined the exact asymptotic form of the work distribution, including the preexponential factor, in a class of stochastically driven systems, using a theory developed by Engel and Nickelsen (EN theory) [1]. In cases where the exact finite time work distribution is not known, EN theory can be applied to obtain good analytical approximations for the tails of the distributions, which are generically hard to observe in experiments or numerical simulations. This can then be combined with data from experimentally / numerically viable regimes to construct the full probability distribution. The extension of EN theory to stochastically driven systems has not been done previously and the analytic solutions that we obtain here are new. EN theory involves writing an augmented action which carries all the information about initial conditions as well as the functional -W [x(·)] whose probability distribution is to be computed. The asymptotic form of work distributions are then computed by using the saddle point approximation, which formally corresponds to the small noise limit. For a class of quadratic (augmented) actions, we have shown that the smallest roots of a certain characteristic polynomial function F (k) determines both the leading-order asymptotic form as well as the pre-exponential factor. The Crooks fluctuation theorem is then shown to manifest itself as the reflection symmetry property of this function. These features can be explained using the relation (proved in Appendix D), for the exact moment generating function (MGF). For a colloidal particle in a harmonic trap, where the mean position of the trap is modulated according to the Ornstein-Uhlenbeck process, we have shown that the (dissipated) work distributions have the asymptotic behaviour, in both the transient and the steady state. Here the constants, C 1 and C 2 are fixed by the time duration of the driving and the noise coefficients, and are explicitly determined for all cases. The asymptotic form given by Eq. (74) can be shown to be universal for quadratic augmented actions. A rigorous discussion of this point and also the relation of the asymptotic form to the singularities of the MGF may be found in [37,16]. In [37], the authors have shown that for the same class of systems for which the above discussions apply, determination of various PDFs and MGFs reduce to finding solutions of certain nonlinear differential equations (NLDEs), which in many cases need to be solved numerically. The method of functional determinants simplifies this problem to instead determining the solutions of ELEs, which are linear ordinary differential equations. In [17], for the special case of a Brownian particle in a logarithmic harmonic potential, the authors obtained the asymptotic form of work distribution in terms of the solution of a Riccati differential equation. We are not aware of any other analytic methods for performing the exact finite-time computation of the asymptotic form including the pre-factor, in Langevin systems. Although we have restricted ourselves to the computation of asymptotics of work distributions in this work, Eq. (73) contains more information, such as the complete moment hierarchy. These aspects will be discussed in a future publication. We believe that the methods that we discuss here has potential applications in the context of finite time thermodynamics of stochastic systems. It should be interesting to look into more applications of this theory, particularly in other stochastic potentials typical for experimental situations. For example, the theory can be developed to include stochastic driving governed by discrete time stochastic processes such as the one studied in [25], by using an appropriate path integral representation [40,41,42]. However, identifying universality classes for asymptotics of work distributions still remains an open problem.

A The Jacobian
In this Appendix, we obtain the exact form of the Jacobian of the transformations that is required in deriving Eq. (56) in Section 3.1. We generalise the derivation of Engel and Nickelsen in [1], to the SSP studied in Section 3 described by the Langevin equations,ẋ The propagator of the corresponding Fokker-Planck equation may be written down as, From the normalization condition where,Ṽ we have, After several partial integrations in the RHS of Eq. (79) we obtain, The above Gaussian integral is straightforward to compute and gives, where J is the Jacobian to be determined. Note that A = A iq=0 as given by Eqns. (49) and (50). Therefore we have, This is the result used in deriving Eq. (56). In a similar manner one can show that all the other Jacobians appearing in the main text, have the above form.

B Functional determinants
The necessity of computing functional determinants of certain differential operators arises in many different situations. As we have seen in the main text, computing the leading-order contribution to path integrals is one of them. In many cases it is not an absolute functional determinant that is required, but a ratio where one or both of the operators can in principle have a zero-mode. Profound mathematical techniques have been developed to compute functional determinants (or the ratio of determinants) even in situations when the operators have zero-modes. Here we apply the contour integration method suggested in [28] to the SSP problem discussed in Section 3. As we have already seen, the operators that appear in the ratio of determinants in Eq. (56) are 2 × 2 matrix differential operators. In [28] Kirsten et al, have discussed the possible generalization of their techniques to such matrix differential operators as well. Recently in [29], Falco et al have also looked at a similar generalization. In contrast to these previous studies, the operator A that we study here has differential operator entries in the off diagonal terms as well. However, the methods discussed in [28] can also be generalized to this situation. In the case of the SSP, the matrix differential operator that we need to find the determinant of, is defined by the following problem: together with Robin-type boundary conditions: The form of the matrices M and N can be deduced from Eq. (50). Using the results from [28], one can then write down the determinant ratio as, Here H k is the matrix of fundamental solutions of the homogeneous equation, defined as, The function B appear due to the presence of the zero-mode and need to be determined using the self adjointness property of the differential operator A in each case. The normalized zero-mode, u N (t) is defined as, where the constants are determined by, The inner product is the usual one, given by In case of the SSP discussed in Section 3.1, we find, In the next Section, we show how this theory can be used to compute the functional determinants that appear in the main text.

C Explicit Computations
Here we provide the explicit calculations using EN theory for two of the cases considered in the main text, the breathing parabola in Section 2.1 and the SSP in Section 3.1. For the breathing parabola problem, in [1], for a specific choice of (reverse) protocol, EN theory was used to compute the exact asymptotic form of P (W ). Here we present the calculations for a particular choice of forward protocol, and give only the final solution for the corresponding reverse protocol. The solutions we obtain for the SSP problem in Section C.2 are new, and are generalizations of the calculations in Section C.1.

C.1 Breathing Parabola: P F/R (W )
For the breathing parabola, we have considered the specific forward protocol, In terms of the shifted variable k = 1 − iq, the corresponding ELE reads (for simplicity we will use x(t) instead ofx(t) everywhere.),ẍ Two independent solutions of this 2nd order differential equations are given by, Together with the boundary conditions, where M and N are matrices, The above system constitutes a second order Sturm-Liouville eigenvalue problem in k. A non trivial solution exists only for some specific values of k, which are given by the roots of the corresponding characteristic polynomial, As we have seen in the main text, the asymptotic behaviour of P (W ) is determined by the smallest value of k(≡ k * ) for which F (k) = 0. The value of k * may be obtained numerically. For the case T = 1 we find k * = 3.67 (see Figure 1 in the main text). The leading-order asymptotic form of P (W ) is therefore, In order to improve this estimate one has to compute the pre-exponential factor as well. The computation of the factor d 2 0 is rather straightforward. From (19) we see that, where x(t) is the zero-mode. For k = −3.67, we find that the zero-mode is Here C 1 is the undetermined constant in the solution, which is to be fixed using the constraint equation (13).
Using the above form of x(t) and the constraint equation (13), we find, The other factor which appears in the calculation of the pre-exponential factor is the determinant ratio of the two functional differential operators, Let us first compute the factor appearing in the numerator of the RHS. Using Eq. (94), the matrix of normalized fundamental solutions (H(0) = I 2 ) when iq = 0 may be found as, Notice that, this is also the limiting value of F (k) as k → 1 in Figure 1. Let us now look at the term in the denominator of the RHS of Eq. (102). The appropriately normalized solutions x N (t) can be identified using (94) and Eq. (89) adapted to this problem. We obtain, x N (t) = (−0.0011 + 0.0035 i)(−2 + t) (0.5−1.55 i) + (19.92 + 61.97 i)(−2 + t) (0.5+1.55 i) .
Using the methods discussed in [28] we find that for this problem, B = − 1 x N (T ) . Together with this, one finds Putting the factors together in Eq. (21), we finally get, Similarly for the reverse protocol, it can be shown that, We compare these results with numerical simulations, and as we show in the main text, Figure 2, the prediction for the tail region is in excellent agreement with the theoretical predictions.

C.2 The Stochastic Sliding Parabola
In this Appendix, we do explicit computations to obtain the asymptotic form of the probability distribution Eq. (59), found in Section 3.1 for equilibrium initial conditions and T = 1 (explicit calculations for the other cases discussed in section 3.1.3 and section 3.2 can be carried out in a similar manner). For simplicity we will use the notation x(t) and λ(t) instead ofx(t) andλ(t). First we note that when |k| = 1, following some algebra, the system of Euler-Lagrange equations for (x, λ) given by Eq. (49) in the main text can be reduced to a fourth order ordinary differential equation for one of the variables (for example, λ) as, In terms of the solution λ(t), x(t) is then given by, Eq. (110) has four independent solutions given by, A general solution for a specific optimal trajectory (k = k * ) can always be written as a linear combination of these four independent solutions, where the coefficients are fixed by the boundary conditions and the constraint equation. As we discussed in Section 3.1.1, in order to compute the leading-order behaviour of P (W ), we only require the k * values and not the explicit solution. In order to find k * , we look at the smallest roots of the function: M and N corresponding to the boundary conditions in (50) can be written down as, H k (t) has the form as given in Eq. (87). From Figure 3 in the main text, we see that the relevant k * values are given by k * = ±2.3. Solving the ELEs (49) along with the boundary conditions (50), for k = 2.3 yields, λ(t) = C 1 (0.0043 sin(0.76t) − 1.30 sin(1.30t) − 0.0056 cos(0.76t) + cos(1.30t)), x(t) = C 1 (−0.0035 sin(0.76t) + 0.62 sin(1.30t) − 0.0083 cos(0.76t) + 1.82 cos(1.30t)).
The work done along this trajectory become, This value is negative, and therefore k * = −2.3 corresponds to the negative tail of P (W ). With this we find that to leading-order, the positive and negative tails of P(W) have the functional form, respectively. In order to improve this estimate, we next include the pre-exponential factor. The first factor which goes into the pre-exponential is d 2 0 defined in Eq. (57). This can be computed relatively easily as in the case of the breathing parabola. For both k = ± 2.3 we find using Eq. (115) and (118), (The superscript ± is used to denote the solutions for positive or negative tails). The next factor to be computed is the square root of the ratio of two functional determinants, for which we will use the formula, where, is the appropriately normalized solution to the ELEs, which have to be found using Eq. (89). For this particular case, We find that, Using the self-adjointness property of A, one can again compute, .
(127) Also using Figure 3 to compute F (1), we find, Therefore the full pre-exponential factor is, and the asymptotic form of P (W ), including the pre-exponential factor becomes, In a similar manner, the exact asymptotic forms discussed in Section 3.1.3 and Section 3.2 can be computed for any value of T . In particular, for the case that we discussed in Section 3.1.3, the exact asymptotic form can be obtained as a function of T .

D Origin of the symmetry of F (k)
In this section we will show the relation between F (k) and the exact moment generating function of dissipated work, which explains the reflection symmetry of F (k). First, using the path integral representation, an exact relation for the moment generating function can be written down as, In all the cases we have considered, the augmented action S[ x, λ, q ] is quadratic, therefore by doing several partial integrations, it can be shown that it reduces to S[ x, λ, q ] = 1 4 x λ A k x λ + Boundary terms in (x, λ, k), k = 1 − iq.
where the kernel A k is defined by the same operator that determines the optimal trajectory (Eq. (11), (49)) with the same boundary terms. Therefore the integral in Eq. (131) is a standard Gaussian integral which can be computed as, This determinant ratio can then be computed using the techniques developed in [28]. In terms of the function F (k), we find, Due to Crooks fluctuation theorem [6], the moment generating function of dissipated work (G) must satisfy the relation, By writing iq as 1 − k and using Eq. (134), it can be immediately verified that for the above relation to hold, the function F must satisfy, F (k) = F (−k).
Hence the symmetry property of F (k) is a consequence of the fluctuation theorem.
E Comparison with the results in [4] In this Appendix, we will invert the generating function of P (W d ) obtained in [4] using the methods discussed in [3] and [20]. The solution that we obtain here will be used for comparison with the asymptotic form of P (W d ) calculated using the EN theory in Section 3.2.
Verley et al, in [4], have shown that for very large T , the probability generating function of the dissipated work has the form, The notation W T stands for the dissipated work W d over a time duration T of the driving. For the SSP considered in Section3.2, the functions φ and g are given by [4], φ(µ) = 1 − ν(µ), where ν(µ) = 1 − µ(1 + µ), g(µ) = 4ν(µ) The probability density function can be obtained from the moment generating function Z(µ, t) by taking the inverse Fourier (two-sided Laplace) transform: where the integration is done along the imaginary axis in the complex-µ plane [20]. Using the large T form of Z(µ, T ) given by Eq. (137) and (138) we write, where The large-T form of P (W T ) can be obtained from Eq. (140) by using the method of steepest descent (for completeness, we reproduce the method and discussion from [20] here ). The saddle point µ * is obtained from the solution of the condition f w (µ * ) = 0 as µ * (w) = 1 2 From the above expression, one finds that µ * (w → + − ∞) → µ ± , where Therefore µ * ∈ (µ − , µ + ). It is useful to notice that in terms of µ ± , On the real axis, outside the interval [µ − , µ + ], ν(µ) is therefore imaginary. However, in order for the the integral in the definition of Z (Eq. (137)) to converge, Z(µ, T ) must be real for real values of µ. For this reason, it is only within the range µ − < µ < µ + ( for which ν(µ) is real and analytic), that the analytic continuation of Z(µ, T ) to real µ is allowed. We hence expect the saddle point to also lie between these values. As we have already seen in Eq. (143), this is indeed the case. Since for µ − < µ < µ + , g(µ) is analytic (the denominator is positive for all µ in the specified range), it can be neglected in the saddle point calculation as a sub-leading contribution.
The saddle point calculation relates φ(µ) to the large deviation function h s (w) by the Legendre transform, We also see that, This means that f w (µ) has a minimum at µ * along real µ. Now since g(µ) is analytic, the usual saddle point approximation method [3] gives, P (W T = w T ) ∼ g(µ * )e T hs(w) 2πT f w (µ * ) .