Explicit solution of the generalised Langevin equation

Generating an initial condition for a Langevin equation with memory is not trivial. We introduce a generalisation of the Laplace transform as a useful tool for solving this problem, in which a limit procedure may send the extension of memory effects to arbitrary times in the past. This machinery allows us to compute average position, work, its variance and the entropy production rate for a particle dragged in a complex fluid by a harmonic potential, which could represent the effect of moving optical tweezers. For initial conditions in equilibrium we generalise the results by van Zon and Cohen, finding the variance of the work for generic protocols of the trap. In addition, we study a particle dragged since a long time by a trap with constant velocity, hence in an established steady state. Our formulas open the door to thermodynamic uncertainty relations for systems with memory.


Introduction
The driven diffusion process of a colloidal bead immersed in a fluid is one of the paradigms of nonequilibrium statistical physics [1,2,3,4,5,6,7]. Fluctuations play a prominent role for this mesoscopic system due to the multitude of random hits on the bead by the molecules of the surrounding fluid. If these molecules are tinier and faster than the colloid, a net separation between the timescales of fast and slow degrees of freedom occurs and the colloidal particle undergoes Markovian dynamics. In this case, the motion of the bead can be equivalently described by using the Langevin equation, path integrals and the Fokker-Plank equation [8]. Historically, the Langevin approach came first and arguably remains the most intuitive. In fact, for a one dimensional system, incorporating the effects of the fluid in Newton's second law, one may write a Langevin equation of motion for the position x(t) of a bead of mass m as a second order stochastic differential equation, The random force is generated by a Gaussian white noise ξ (t), with average ξ (t) = 0 and correlation ξ (t )ξ (t ) = 2γ 0 k B T δ (t −t ). The prefactor of the delta function ensures thermodynamic consistency according to the (second) fluctuation-dissipation theorem [9], linking the drag coefficient γ 0 of the dissipative term −γ 0ẋ to the strength of the noisy term. As a deterministic force not due to the fluid we focus on the case F (x,t) = −∂ x U(x,t) with a time-dependent potential energy U(x,t).
If the bead is immersed in a solution containing for example long and complex polymers [10,11], the above-mentioned separation of time scales is no longer possible and memory effects occur. One may then consider a generalised Langevin equation (GLE) with constant diffusion coefficient, whose formal derivation can be found in [12,13,14]. For t ≥ 0 this equation reads where Γ (t) is the memory kernel, t m ≤ 0 is the time to which the memory effects extend and η(t) is a coloured Gaussian noise obeying η(t) = 0. The above equation could also describe the motion of a particle under the effect of hydrodynamic backflow [15]. The fluctuation-dissipation relation [9] is still valid in the more general form η(t )η(t ) = k B T Γ (|t − t |): thermodynamic equilibrium is present in the medium if its two effects (dissipation and noise) are proportional at all times. Note that a Markovian memory kernel Γ Markov (t) = 2k B T δ (t) would lead to the usual fluctuation-dissipation theorem. The aim of this paper is to solve the GLE with a parabolic confinement potential U(x,t) = κ 2 (x − λ (t)) 2 , obtained for instance by using optical tweezers centred on a moving coordinate λ (t) The non-dynamical case was already discussed for example in [16]. Moreover, we will restrict ourselves to the case of a non-divergent time dependent friction coefficient γ = lim t→∞ t 0 dt Γ (t ) < ∞, which is a sensible physical requirement [17,18].
One of the first analytical solutions for the GLE with κ = 0 and no external force can be found in [19]. It is obtained through the use of Laplace transforms and is expressed in terms of the velocity susceptibility χ v (t), a key quantity discussed in the next sections. In this paper we obtain a more general solution in terms of the susceptibility and its integrals. This enables us to calculate averages and variances of relevant quantities such as position, thermodynamic work and entropy production, with a dynamics starting from different initial conditions. Some of these results are already known in the literature, especially for equilibrium initial conditions, see for example [20]. However, imposing a nonequilibrium steady state as initial condition is not trivial for the GLE, due to its memory. A scheme for achieving an initial condition with memory extending far in the past seems required. To this end, we introduce a modified version of Laplace transforms with arbitrary initial time t m , which is then shifted back to minus infinity by a limit procedure. The explicit dependence of the solution on t m along with the well defined limits of susceptibilities will make this procedure very straightforward.
The following section introduces the technical details of the modified Laplace transform. In Sec. 3 we discuss the solution of the GLE and in Sec. 4 we show how to use the solution for computing interesting thermodynamic quantities. We show that the entropy production rate can be expressed in terms of a retarded velocity, coinciding with the usual velocity of the particle in the Markovian case, see (65). Finally in Sec. 5 we apply the obtained results to transient dynamics from equilibrium and to the stationary state generated by a linear dragging protocol λ (t) = vt, also discussed in [21]. In the latter scenario we note that quantities such as average position, velocity, work and entropy production rate have the same structure as for Markov dynamics. Moreover, we discuss the limit t m → −∞ for an arbitrary dragging protocol λ (t) which can be seen as a generalised stationary state, in the sense that the memory of initial conditions is lost. In this case we manage to show that the variance of the thermodynamic work is equal to the one of a system prepared in equilibrium initial conditions (see equation (88)), thus generalising the results by van Zon and Cohen [2].

Modified Laplace transform
A standard way of dealing with the linear GLE uses Laplace transforms. This technique is particularly useful when dealing with an initial condition at finite times, for instance when the system starts from equilibrium at time t = 0. If the initial time is rather taken infinitely back in the past, traditional Laplace transforms are no longer suitable to find a solution for the GLE. However, it is well known that, for Markovian dynamics, non-equilibrium steady states can be obtained from this limit. Hence, we would find it useful to have a framework in which Laplace transforms are available and steady states may be considered.
Our way to tackle this problem is to introduce a modified Laplace transform with an arbitrary initial time t m ≤ 0 that acts on a given function g(t) as followŝ The standard Laplace transform of course is recovered for t m 0. The aim is to solve the GLE finding the explicit dependence of the solution on t m and then, if interested in steady states, eventually take the limit t m → −∞. For our purposes, we just need to know the effect of such modified transform on first and second derivatives of a function. They can be readily expressed as Note thatĝ t m (k) stands for the modified Laplace transform of the function g(t) while g t m ≡ g(t m ) is the function calculated at time t m . Furthermore, it is not hard to show that the action of the modified Laplace transform on integrals is equal to the action of the standard transform, namely We also need to know the effect of such transform on the convolution of a causal function (in our case the memory kernel Γ (t) = 0 for t < 0) with an arbitrary g(t), First, to compute an explicit version of this equation, we note that i.e. these integrals define the same region of integration so that (7) becomes which is a generalisation of the convolution theorem. It states that the modified Laplace transform of the convolution of a causal function with an arbitrary function g(t) is equal to the product of the standard Laplace transform of the causal function, i.e.Γ (k), and the modified Laplace transform of g(t), that isĝ t m (k). We conclude this section by remarking that, of course, the modified Laplace transform of a causal function is equal to the standard Laplace transform of that function.

GLE solution
By applying the modified Laplace transform (4) to the GLE (3) and by using the results obtained above we get Furthermore, with a bit of algebra we can isolate the position x from the other quantities obtaininĝ where we introduced the "position susceptibility" χ x (t), a key quantity of this paper, defined via its Laplace transform In the following we will also use its integral χ(t) and its derivative χ v (t) ("velocity susceptibility"), In Appendix A we discuss the limits of these susceptibilities for t → 0 and t → ∞. Two examples are shown in figure 1.
We stress that all the susceptibilities are of course causal functions. Noting that L t m k e −kt m = δ (t − t m ) and L t m k e −ktm k = θ (t − t m ), where θ (t) is the Heaviside step function, and transforming back equation (11) to real time, we obtain that is the solution to the Langevin equation. The average of the position can now be obtained, given that η(t) = 0: (17) From this the average velocity can be calculated using that the system is linear and Here we exploit that in the underdamped limit χ(0) = 0. For overdamped dynamics, instead the velocity of the particle is not defined.

Variance of the position and correlations
Another important quantity we are interested in is the variance of the position at time t. Given that the system started at time t m with position x t m and velocity v t m , we have that Using the previously obtained expression for the position (16) and defining we find that (19) becomes Focusing on the the first term on the right hand side, we further define and calculate the following quantity (also for future convenience): In the last line we used the second fluctuation-dissipation theorem η(t )η (t ) = k B T Γ (|t − t |) that relates the correlation of the noise to the memory kernel. Taking the double Laplace transform of both sides of equation (22) we get where β = 1/k B T as usual. Moreover, we used (8) between the second and the third line and then we made the change of variable u = t − s. Focusing on the remaining integrals, we have that where in the last line we recognised the Laplace transform of Γ (t) and used that As for the second term in the last line of equation (24), using integration by parts we get where we noted that the first term in the second line is equal to zero. Going back to equation (24) and remembering that we started from (23) we finally obtain Recalling the definition of the position susceptibility via its Laplace transform and its relation with the memory kernelχ x (k) = [mk 2 + kΓ (k) + κ] −1 and doing some algebra it is possible to show that The inverse transformation back to real time yields Using that along with the generalised convolution theorem, we are able to show that (29) becomes To finally obtain an expression for the variance of the position, we evaluate this quantity at equal times (i.e. t = t = t ) and then plug it into equation (21). Since Finally, by using (21), we obtain an expression for the variance of the position from arbitrary initial conditions, Note that, being the GLE linear, if the initial probability distribution function (PDF) P(x t m , v t m ,t m ) is a (bivariate) Gaussian, so will be the P t m (x t , v t ,t) at time t > t m . In particular, for overdamped dynamics, in order to fully characterise the PDF P t m (x t ,t) (velocity is not defined in this limit) one only needs the average and variance of the position, which can be readily obtained from (17) and (33) by taking the massless limit m → 0. Instead, for underdamped dynamics, one also needs the variance of the velocity and the covariance between position and velocity They can be calculated similarly to the variance of the position (33), where, in order to calculate (36), we used that χ v (0) = 1/m and the convention for the Heaviside step function for which θ (0) = 1/2. Equations (33), (36) and (37) are the explicit expressions of the components of the covariance matrix used to define the PDF of the process at time t with Gaussian initial conditions at time t m with x t = (x t , v t ) and x t m ,t = ( x t m ,t , v t m ,t ) given by (17) and (18).

Thermodynamic quantities
This section is devoted to the analysis of relevant thermodynamic quantities such as work, entropy production and entropy production rate.

Work
We consider the definition according to stochastic energetics [5,6] of work done on a particle by a time dependent external potential during the time interval [0,t], where we used λ (0) = 0. We can calculate the work as a function of the external protocol and the susceptibilities (13) and (14) by just plugging the explicit solution for the position of the particle (16) into (40), which reads Its average can be obtained, again by noting that η(t) = 0, as It is well known that, for such linear systems, the PDF P(W ) of the work is Gaussian.
Hence, in addition to the average W t m ,t , again we need its variance to completely characterise the PDF. It can be calculated similarly to the variance of the position (19), starting from the definition of work (40), where C (t ,t ) was defined in (31). By calculating the first term in the second line we get that is the expression for the variance of the work for an arbitrary initial distribution of position and velocities. Although it might look rather complicated, in the next section we will see that the above equation simplifies significantly for some usual initial distributions.

Entropy production and entropy production rate
Entropy production does not need any introduction, it is a crucial quantity in stochastic thermodynamics that encodes the information about the irreversibility of a given process. In particular, for a colloidal bead in contact with a heat bath, the entropy production during a time interval [0,t] can be split into two parts where Q is the heat injected into the heat reservoir, β is the inverse temperature (hence Σ med (x,t) is the entropy change in the reservoir) and Σ sys (x,t) is the difference between the Shannon entropy of the final and initial states of the system. In particular, for Gaussian PDFs, it holds that where |S t m ,t | is the determinant of the covariance matrix (38) at time t.
As for the heat absorbed from the bath, it can be defined through the Stratonovich integral [5,6] where F bath (x,t) is the force exerted from the particle on the bath, i.e., using the GLE (2), Equation (48) thus becomes Taking its average, since r 2 = ∆ 2 r + r 2 for an arbitrary stochastic variable r, we get At this stage one cannot simplify further this expression for the entropy production.
On the other hand, we can obtain a much more compact form for the entropy production rate, defined as For the system entropy production rate we immediately see from (47) that From equation (51) instead we get that where again we used that v t m ,t = ẋ t m ,t = ∂ t x t m ,t . Consider now the term between square brackets on the right hand side of equation (54) and name it Taking its modified Laplace transform we obtain where we used the formula for the modified Laplace transform of a second derivative (5). Moreover, looking back to the expression for the average of the position (17) we note that it can be effectively written as where I (t,t m ) = x t m (1 − κ χ(t − t m )) + m v t m χ x (t − t m ) contains the information relative to initial conditions, in particular I (t m ,t m ) = x t m . Going back to equation (56), recalling the definition of the position susceptibility via its Laplace transform (χ x (k) = [mk 2 + kΓ (k) + κ] −1 ) and using that for the first term on the right hand side of equation (57) we have that along with the generalised convolution theorem for the second one, we get Its inverse can be calculated using again the convolution theorem is the time dependent friction coefficient and γ = lim t→∞ γ(t) is its asymptotic limit for long times. As for v ret t m ,t , it can be interpreted as a retarded velocity that converges to the real velocity for t → ∞. In fact The same decoupling between the kernel and the average velocity can be obtained for t m → −∞ if one is able to show that v t m ,t = v t−t m . It will be for example the case of a trapped particle dragged at a constant velocity, i.e. λ (t) = vt. Under these hypothesis, with a calculation analogous to that of equation (62), one can show that the above relation can be more generally written as Moreover, note that for Markovian dynamics defined by a memory kernel Γ Markov (t) = 2γ 0 δ (t) it holds that γ = γ(t) = γ 0 and v ret t m ,t = v t m ,t for every t.
Finally, putting together equation (54) and (60), we get while for the total entropy production rate (assuming that P t m (x t , v t ,t) is Gaussian) we have that

Applications
In this paragraph we apply the general formulas derived in the previous sections to specific initial conditions. In particular, we will discuss a transient dynamics from equilibrium and a steady state dynamics.

Transient from equilibrium
For a colloidal particle trapped in a parabolic potential with stiffness κ, the equilibrium PDF at time t m = 0 has a Gaussian shape, with parameters given by Using equations (67) and (68), we can evaluate the evolution of all the quantities discussed in the previous section, starting from the probability distribution defined above and for an arbitrary λ (t). Starting from the average of the position (17) and velocity (18) we find that while for the covariance matrix, using equations (33), (36) and (37) we get that i.e. if we start from equilibrium and the trap stiffness κ does not change, then the covariance matrix remains constant in time.
Going forward to the estimate of thermodynamic work, from (42) and (44) and again using that λ (0) = 0 along with χ(t) = t 0 dt χ x (t), we get that Being the PDF of the work P(W t ) Gaussian, an integral fluctuation theorem for the thermodynamic work W (x t , v t ,t) holds (see [20] for details) and a Jarzynski equality would follow [22]. Finally, being the covariance matrix and its determinant both constants, a very simple expression can be found for the rate of entropy production

Stationary State
The first thing we do in this section is to clarify under which circumstances the system reaches a stationary state when trapping is active and either the initial time t m → −∞ or the observation time t → +∞. One usually defines the stationary distribution as the solution of the Fokker Planck equation when the PDF does not depend explicitly on time. Nevertheless, this definition becomes problematic when the drift term or the diffusion coefficient of the associated Langevin equation depend explicitly on time, as in the cases we are considering in this paper. In order to tackle this problem, first of all we note that if a sufficiently large time has passed from the beginning of the dynamics, i.e. t − t m → ∞, the PDF P t m (x t , v t ,t) at time t will be a bivariate Gaussian with the usual form depending on time via the averages of position and velocity and the covariance matrix. For the latter we again use the expressions for the variance of position and velocity (33) and (36) alongside with their covariance (37) and the limits of susceptibilities that can be found in Appendix A, finding that As in the previous example starting form equilibrium, also for this sort of steady state we have that the covariance matrix does not depend on time. However, the averages of position of velocity will be in general dependent on time, as it can be seen from their expression (17) and (18) and proceeding similarly as above We outflank this problem by moving the centre of the harmonic trap at constant speed, i.e. λ (t) = vt, so that we get the following GLE Performing the change of variable y(t) = x(t) − vt, we see that the system can be mapped through a Galilean transformation to the centre of the trap reference frame. Note that this does not change the covariance matrix and that the new PDF P t m (y t ,ẏ t ,t) will be defined from the same matrix and by y t m ,t and ẏ t m ,t . The transformed GLE becomes and its solution can be found similarly to that for the original GLE. In particular we find that (80) Taking the limit t − t m → ∞, which corresponds to either t m → −∞ or t → +∞, and using the limits derived in Appendix A, we see that which are both constant. We conclude that for a harmonic potential with constant strength and with centre travelling at constant speed (λ (t) = vt) it is possible, through a Galilean transformation, to map the system to another one for which an equilibrium distribution exists. In fact, the PDF P t m (y t ,ẏ t ,t) inherits the Gaussian character from the PDF of the original variable x(t). This in turn implies that, being the covariance matrix constant as well as the averages of the dynamical variables (81), the PDF for y(t) becomes time independent. In this sense we mean that P t m (x t , v t ,t) becomes stationary as t − t m → +∞.
In addition, we note that i.e. they do not depend on the specific form of the memory kernel but only on the limit of its time integral. Moreover, note that the expression above has the same structure as in the usual Markov case where instead of γ there appears the conventional Stokes friction coefficient γ 0 . Consider now the thermodynamic work, in particular equations (40) and (44) for the specific case of λ (t) = vt. For the average work we find that again has the same form as the well known Markov case. For the variance of the work, instead, we use the limits of susceptibilities discussed in Appendix A i.e. from (72) and for λ (t) = vt we get that ∆ 2 W ss t = ∆ 2 W eq t . As for the entropy production rate we immediately see that it has the same form as for Markov dynamics with the usual substitution γ 0 → γ σ tot ss t = γv 2 (85) We conclude this section by highlighting that if the dragging protocol λ (t) = vt a "stationary state" in a generalised sense can be achieved as t m → −∞, meaning that memory of initial conditions is lost and that the covariance matrix has become constant, see (76). This can be easily seen again by considering the limits of the susceptibilities discussed in the appendix. For position and velocity, using equation (77) what we get is while for the average work it holds that As for its variance, from (44) and using that χ(∞) = 1/κ along with χ x (∞) = 0, we can easily see that  meaning that the equivalence for this quantity starting from equilibrium initial conditions and the stationary state in the strict sense, also holds for this generalised stationary state and for arbitrary driving protocol λ (t). Finally, for the entropy production rate we use equation (65) along with the fact that the covariance matrix is constant in order to obtain 5.3 Example: exponentially decaying memory kernel As a standard example for non-Markovian dynamics, we consider a GLE with exponentially decaying memory kernel, and, for causality, Γ exp (t) = 0 for t < 0. The characteristic time τ emerges from the relaxation of the molecules in the reservoir. In the limit τ → 0 the exponential memory kernel tends to a Dirac delta  We set m = 1, κ = 1, γ = 1, A = 1 and ω = 1. In both scenarios we compare the Markov case (τ → 0) with an exponential memory case (τ = 2). The amplitude of the oscillations is wider in the latter case. and the Markovian limit is recovered.
For finite τ the susceptibilities display oscillations, as shown in figure 1. Hence, quantities considered in the previous section may also exhibit oscillations. This is the case for all quantities displayed in figure 2(a), for a system starting from an equilibrium condition, even if the dragging protocol λ (t) = vt is linear. In the stationary state, memory effects are not visible anymore in the averages of position, work and entropy production rate (they grow linearly, see figure 2(b)) but oscillations are still present in the variance of work, which we have shown to follow the same formula for transient dynamics and for the stationary state. The non-monotonicity with time of the work variance is clearly due to the memory stored by the complex fluid. The variance of position and velocity are not shown in the figure as they are constant in both cases.
Finally, if we we consider an intrinsically oscillating driving protocol of the form λ (t) = A sin(ωt), the effects of memory may determine an increase of the amplitude of the already present oscillations, both from equilibrium ( figure 3(a)), and in the steady state ( figure 3(b)). Panels on the left in figure 3 represent the Markovian limit while panels on the right show an example for an exponential memory kernel with τ = 2. In the latter case, the average position fluctuates more, and the entropy production rate can become negative (while having a positive average over one cycle in the steady oscillatory regime).

Conclusions
The Gaussian process with memory is a classic in statistical mechanics. Yet, we have shown that further results can be derived for this process realised by a generalised Langevin equation for a particle driven by a harmonic strap with constant strength in a complex fluid. An explicit solution of the GLE is based on computing susceptibilities. In terms of these important dynamical quantities, several other expressions are derived.
For generic protocols and initial Gaussian conditions, the quantities we computed for every time t ≥ 0 are the average particle position (17), its autocorrelation function (31) and hence its variance (33), the average work done on the system (42), its variance (44), and the entropy production rate (65). These formulas can be simplified in some standard scenarios, e.g. starting from equilibrium or in steady states. Moreover, the variance of the work starting from equilibrium is equal to that for a steady state in a generalised sense and is proportional to the average of work starting from the same initial conditions. Since we can deal with various dragging protocols, this means that the two cumulants for the work (73) generalise formulas by van Zon and Cohen [2].
Especially aiming at dealing with steady states, everything starts by introducing a new Laplace transform with arbitrary initial time t m . The explicit dependence of the solution on t m along with the well defined behaviour of the susceptibilities for the limits t → 0 and t → ∞ allow us to recognise a steady state for a linear dragging protocol λ (t) = vt as t m → −∞. More in general, for an arbitrary protocol, this limit leads to a loss of the information about the initial state. We can interpret it as a generalised steady state.
Going into some more details about the quantities calculated throughout the paper, for a steady state generated by a linear dragging protocol we recognise the same structure of the average of position and velocity as for usual Markov dynamics and same thing for the covariance matrix. Finally, we are able to write the entropy production rate in terms of a quantity that we termed the retarded velocity, matching the usual velocity if no memory effects are included in the kernel.
In conclusion, we note that this framework yields average quantities but also their variances. Hence it will be used to derive one of the first examples of thermodynamic uncertainty relation [23,24,25,26,27,28,29,30] for systems with memory [31,32].

A Appendix: Limits of susceptibilities
In this section we discuss the limits of the position susceptibility defined in Laplace space aŝ χ x (k) = [mk 2 + kΓ (k) + κ] −1 (92) along with the limit of its integral and and of its derivative, To this end we use that for a given function g(t) it holds where the last line immediately follows from the first line. Note that all this limits do not depend on m and hence they are valid also for overdamped dynamics. Things become different in the limit of t → 0. where we used that lim k→∞ mk 2 kΓ (k) 1. In factΓ (k) k→∞ = k 2 would correspond to ballistic motion which we do not consider, see [17] for more details. As for its integral and derivative of course we have that lim t→0 We see that this result does not depend on the kernel form, in fact inertial effects dominate the particle behaviour in the small time limit.