Work fluctuations for a confined Brownian particle: the role of initial conditions

We study the large fluctuations of the work injected by the random force into a Brownian particle under the action of a confining harmonic potential. In particular, we compute analytically the rate function for generic uncorrelated initial conditions, showing that, depending on the initial spread, it can exhibit no, one, or two singularities associated to the onset of linear tails. A dependence on the potential strength is observed for large initial spreads (entailing two singularities), which is lost for stationary initial conditions (giving one singularity) and concentrated initial values (no singularity). We discuss the mechanism responsible for the singularities of the rate function, identifying it as a big jump in the initial values. Analytical results are corroborated by numerical simulations.


Introduction
Rare events represent an active research field in physics, mathematics and natural sciences [1,2].One way to describe them is suggested by large deviation theory [3,4], which provides * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. a quantitative evaluation of probabilities of events beyond the regime of normal fluctuations.This theory also offers the possibility to extend the usual equilibrium free-energy approach of statistical physics to non-equilibrium and dynamical contexts [5].In fact, given an extensive physical observable W τ computed by cumulating a large number τ of microscopic events, if a large deviation principle holds, then the asymptotics of the probability distribution P(W τ /τ = w) can be characterised by the rate function I(w) = − lim τ ↑∞ 1 τ log P(W τ /τ = w), playing indeed a role similar to a free-energy.
A singularity in the rate function can separate the regime of typical fluctuations around the mean from the regime of far rare events.In dynamical contexts, singularities of this kind are associated with dynamical phase transitions [6][7][8][9][10][11][12][13][14][15][16][17][18].In particular, if W τ is an observable measured along the trajectories of a given system, a singularity in the graph of the rate function I(w) can mark a phase separation in trajectory space.A relevant question is which are the trajectories contributing to the different regimes of I(w).A possible answer could come from the so-called single big-jump principle [19][20][21][22], which explains rare events not in terms of an accumulation of many small microscopic events but solely as an effect of the biggest one.It would be interesting to analyse this principle in the context of dynamical phase transitions.
The dynamics of Brownian particles offers the possibility to deepen our understanding of the above subjects.Farago [23] studied the work done by the thermal bath force on a single underdamped Brownian particle.In particular, the rate function of the work for a free particle was computed and, in the stationary regime, a singularity and an associated linear tail were found at negative values of the work.The consequences on the non validity of the Fluctuation Relation were discussed.Then, the effects of a confining potential were considered.A confining potential at the level of a single particle, in addition to having an intrinsic interest as experimentally realisable by means of an optical trap, can mimic the trapping created by other particles at finite densities, and has been studied with also this purpose [24][25][26].Arguments were given in [23] to conclude that the presence of a harmonic potential would not modify the rate function for the thermal bath work with respect to the free case.Singular rate functions for trajectory dependent quantities have also been found for Brownian particles in a moving potential [27,28] or in contact with several baths [29], for Brownian particles under the action of an additional Gaussian force [30,31], for single harmonically confined active particles [32] and for active Brownian particles at finite densities [33][34][35][36][37].
In this paper, we consider again a harmonically confined Brownian particle and compute analytically the rate function of the work done by the random force of the thermal bath.We extend the analysis to the case of generic uncorrelated initial conditions.Our results prove the claims made in [23], as they show that the rate function does not depend on the strength of the potential in the stationary regime and for concentrated initial values.For generic initial conditions, we demonstrate that the rate function always possesses a singularity at a value of the work smaller than the mean work, except for the limiting case of concentrated initial values where such singularity moves to the boundary.This singularity corresponds to the singularity observed in [23] for the stationary, free particle.In addition, we demonstrate that a second singularity emerges above the mean work when the spread of the initial conditions is large enough.Each singularity is the beginning point of an associated linear tail in the graph of the rate function.We provide an interpretation of the singularities in terms of the single big-jump principle, showing that trajectories in the linear tail regimes are characterised by a big jump in the initial values.
The paper is organised as follows.In section 2 we present the model and introduce the work injected by the random force.In this section, we also report a brief summary of the analytical approach developed in [38], which we use to compute the rate function.Section 3 is devoted to the computation of the scaled cumulant generating function (SCGF) and the rate function, the latter being the Legendre transform of the former.In section 4 we investigate on the trajectories phenomenology in order to explain the mechanism that causes the singularities of the rate function.Finally, section 5 summarises our findings and future perspectives.

Definition of the model
We consider a single unit-mass Brownian particle in one dimension under the action of a harmonic potential described by the following Langevin equation where x(t) is the particle position, γ is the viscous friction coefficient, k (the potential strength) is the elastic constant, η(t) (the random force) is a Gaussian white noise with ⟨η(t)⟩ = 0 and ⟨η(t)η(s)⟩ = δ(t − s) and D is the diffusion coefficient.Following [23], we interpret −γ ẋ(t) and √ 2D η(t) as generic energy dissipation and injection channels, respectively, so that throughout calculations we keep D unspecified.The usual case of an equilibrium thermal bath at temperature T is recovered by fixing D = γk B T with k B the Boltzmann constant, as prescribed by the Einstein relation.Introducing the velocity v(t), equation ( 1) is recast into the following system of first-order differential equations with initial conditions x(0) and v(0), to which we will refer as the equations of motion.We consider uncorrelated Gaussian initial conditions with zero mean and standard deviations σ x for the position and σ v for the velocity.The case of a stationary process is obtained with σ 2 x = D/kγ and σ 2 v = D/γ [39].We will study the probability distribution of the work injected by the random force into the Brownian particle up to time τ , defined as where in the second row we have used the η(t) expression extracted from equation (2).Our goal is in fact to compute analytically the rate function I(w) of the work injected per unit time: ) , where P(W τ /τ = w) is the probability distribution expressed by the path integral with path probability which weighs each trajectory realisation combining the distribution of the initial values with the Onsager-Machlup weight [40].To compute the rate function I(w), we follow the approach developed in [38], which is briefly summarised in the next subsection.Application of this approach to our problem is made in section 3. We mention that an alternative approach to deal with the work injected by a Gaussian external force into an underdamped Brownian particle is suggested in [28,30,31].Although both methods are ultimately traced back to Gaussian integrals for computational purposes, the theory developed in [38] provides, with mathematical rigour, a large deviation principle for a general class of quadratic functionals and a comprehensive recipe for evaluating their rate functions.

Overview of the analytical approach
In order to evaluate the rate function I(w), we resort to the approach recently proposed in [38] for quadratic functionals of stable Gauss-Markov chains, whose applicability to our problem is feasible upon a convenient time discretization.Results for the original problem are later obtained by performing a continuum limit.Let {X n } n⩾0 be a Markov chain taking values in R d (in our case d = 2, corresponding to position and velocity), and assume that there exist a drift matrix S with spectral radius ρ(S) < 1 and an invertible, diagonal diffusion matrix C such that where {G n } n⩾0 is a sequence of i.i.d.standard Gaussian random vectors valued in R d .Suppose that the initial state X 0 is a Gaussian random vector independent of {G n } n⩾0 with zero mean and positive-definite covariance matrix Σ 0 .The chain ( 4) is stationary if and only if ) m , → p denoting transposition, and the request on the spectral radius ensures that Σ s actually exists.Let also be a quadratic functional of the chain, where ⟨•, •⟩ is the Euclidean scalar product and L, U, R, and V are d × d real matrices with L, U, and R symmetric.With full probability, the typical value of W N is given by the following law of large numbers [38] lim The theory developed in [38] describes the large deviations of W N .To expose it, we need to formulate some mathematical details.For each real numbers µ and θ we introduce the Hermitian matrix where ı is the imaginary unit, and we term primary domain the set O of µ for which The matrices with I the identity matrix, can be shown to be invertible [38] and can be proved to be Hermitian [38].The interval E ≡ (µ − , µ + ) ⊆ O in which L µ and R µ are simultaneously positive definite is termed effective domain.We are now in the position to state the following large deviation result [38].The SCGF of W N /N in the large N limit, in symbols lim N↑∞ 1 N ln⟨e µWN ⟩, turns out to be the function (8) with domain Moreover, the quadratic functional W N /N satisfies a large deviation principle with rate function J(w) given by the Legendre transform of the function ( 8) in E, i.e.
Within the effective domain E, the SCGF is differentiable, and the limits lim µ↓µ− φ ′ (µ) ≡ w − and lim µ↑µ+ φ ′ (µ) ≡ w + exist by convexity.When μ− < µ − and/or µ + < μ+ , the SCGF is non-steep at the boundaries, i.e. w ± are finite, and the general expression for the corresponding rate function is where j(w) is a convex function connecting continuously at w ± along with its first derivative.The rate function is characterised by second-order singularities at w ± and linear tails outside the interval (w − , w + ).Thus, whenever the effective domain is smaller than the primary one, at least one linear stretch in the rate function occurs.We stress that the matrices Σ 0 , L, and R neither affect the typical value of W N nor the SCGF function expression and the primary domain.Instead, they play a crucial role in determining the effective domain extension, and hence in determining the occurrence of singularities and linear tails in the graph of the rate function.

SCGF and rate function
In this section we apply the above method to our problem.Section 3.1 introduces a convenient time discretization.Section 3.2 gives the SCGF, whereas the effective domain is discussed in section 3.3.Finally, the rate function is investigated in section 3.4.The reader interested only in the latter can directly go over to the last subsection.

Discretization
By time discretization, we recast the equations of motion (2) and the work (3) as a Gauss-Markov chain and a quadratic functional, respectively.To make a contact with the framework of [38], time discretization must provide an invertible, diagonal diffusion matrix C. To this aim we introduce a fictitious noise in the upper equation of equation ( 2) with vanishing diffusion coefficient in the continuum limit.Specifically, we divide the time interval of duration τ in equation ( 3) into the sum of N time steps of size ϵ, so that τ = Nϵ, and approximate the position and velocity derivatives as dx (t)  dt ≃ xn+1−xn ϵ and dv(t) dt ≃ vn+1−vn ϵ , with discrete-time position x n ≡ x(nϵ) and discrete-time velocity v n ≡ v(nϵ).In this way, we turn the equations of motion (2) into the following Markov chain where {ξ n } n⩾0 and {η n } n⩾0 are two independent sequences of i.i.d.standard Gaussian random variables and D F is a fictitious diffusion coefficient that goes to zero when ϵ ↓ 0. The precise dependence of D F on ϵ is irrelevant, but we suppose that D F ∝ ϵ for simplicity.Results for the original continuum system will be obtained by performing the continuum limit ϵ ↓ 0. The drift and diffusion matrix extracted from equation ( 13) are The eigenvalues of S are 1 − ϵ 2 (γ ± √ γ 2 − 4k), so that ρ(S) < 1 for all sufficiently small ϵ.Thus, the stationary covariance matrix exists for all sufficiently small ϵ and reads The leading order corresponds to the stationary covariance matrix for the continuum system [39], as expected.
Concerning the work injected by the random force into the Brownian particle, the integral contribution in equation ( 3) is discretized according to the trapezoidal rule for convenience as Once W N is recast as in equation ( 5), we find From equation ( 6), with full probability the typical value of the discrete-time work turns out to be lim ) .
This formula and τ = Nϵ give the typical value of the work for the continuum system: ⟨w⟩ ≡ lim τ ↑∞ W τ /τ = D. To conclude, we observe that if J(w) ≡ − lim N↑∞ 1 N log P(W N /N = w) denotes the rate function of W N /N in equation (12), then the rate function I(w) of W τ /τ is

SCGF
Here we aim to evaluate the SCGF φ(µ) in equation (17).The first step is the knowledge of the Hermitian matrix F µ (θ) defined by equation (7).From equations ( 14) and ( 16), we find where the explicit expressions of the coefficients are We note that only the coefficient p 12 depends on the variable µ.The determinant of F µ (θ) is given by where According to Sylvester's criterion, the Hermitian matrix F µ (θ) is positive definite if and only if its upper left entry and its determinant are positive.The former is positive for all θ, µ, and ϵ > 0. For ϵ > 0 small enough, the latter is positive for all θ provided that µ ∈ (μ − , μ+ ) with Therefore, recalling the discussion in section 2.2, the primary domain of the SCGF φ(µ) is O = (μ − , μ+ ).Bearing in mind that D F ∝ ϵ, we note that We can now compute the function φ(µ).For µ ∈ O, we have the identity with Thus, combining equation ( 8) with equation ( 19), we get the explicit expression ) ) .
We define to which we refer as the SCGF for the continuum system.This function is defined on the primary domain (−∞, γ/4D), which, manifestly, is the continuum limit of equation ( 18).

Effective domain
Here we establish the effective domain E in equation (17).To this aim, the explicit expressions of the matrices L µ and R µ in equation (11) are needed.This requires at first to compute the matrices Φ µ (0) and Φ µ (1) in equation ( 9) for µ ∈ O.All the involved integrals in Φ µ (0) and Φ µ (1) fall within the following general case: where g 1 and g 2 are as in equation ( 20), whereas ζ 0 , ζ 1 , and ζ 2 are generic numbers.Thus, for µ ∈ O, we find ) .
Combining these expressions with equations ( 10) and ( 11), expanding around ϵ = 0 with D F ∝ ϵ, a lengthy calculation that we omit gives the following simple result for the leading order of the matrices L µ and R µ : The interval E ≡ (µ − , µ + ) accounts for those µ ∈ O that make L µ and R µ positive-definite at the same time.It is straightforward to check that, at small ϵ, the matrix R µ is positive definite for all µ in the primary domain O. Instead, at small ϵ, the matrix L µ is positive definite only if both the diagonal entries of the leading order are positive, that is only if In this way, the effective domain of the SCGF in the continuum limit, which is ϕ(µ) given by equation ( 21), is the set of numbers µ < γ/4D satisfying equation (22).Solving equation ( 22), for the continuum system we finally find E = (µ − , µ + ) with We note that µ − > μ− = −∞ for every choice of the parameters, so that ϕ(µ) is non-steep at µ − in any situation, except for the limiting case of concentrated initial values defined by the limits σ x , σ v ↓ 0. With ϕ(µ) given by equation ( 21) and µ ± given by equation ( 23), from equation ( 17) we get

Rate function
Here we investigate I(w) given by equation (24), which constitutes the rate function for the continuum system and is one of the main objectives of the paper.Using equation ( 21) for ϕ(µ) and equation ( 23) for µ ± , when M ≡ max{kσ 2 x , σ 2 v } ⩽ 4D/γ we find If instead M ≡ max{kσ 2 x , σ 2 v } > 4D/γ, then we get , where w − is as before and We realise that the rate function always has a singularity at w − below the mean D with an associated left linear tail.Moreover, when M > 4D/γ the model exhibits a second singularity at w + above D and a corresponding right linear tail.We stress that the parameter M is the only place where the elastic constant k appears.
In the special case of concentrated initial values defined by the limits σ x , σ v ↓ 0, i.e. when the Brownian particle starts fluctuating from x(0) = v(0) = 0, we have M = 0 ⩽ 4D/γ, so that the rate function I(w) reads In this limiting case the rate function does not display linear tails since the assumption x(0) = v(0) = 0 rules out the possibility of negative fluctuations for the work, as it becomes a sum of positive contributions according to equation (3).We note that the rate function is independent of the elastic constant k in this limiting case.The rate function is independent of k also in the case of stationary system, corresponding to σ 2 x = D/kγ and σ 2 v = D/γ.In fact, for a stationary system we find M = D/γ ⩽ 4D/γ, giving The left linear tail, and that alone, is now present.Equations ( 25) and ( 26) demonstrate the claim made in [23] that the injected power is completely insensitive to the presence of a harmonic confinement under concentrated initial values or stationary initial conditions.Consistently, they reproduce the analytic results of [23] for the free system corresponding to k = 0.As observed in [23], the rate function satisfies no Fluctuation Relation [41][42][43][44].Figure 1 reports the SCGF ϕ(µ) (top) and the rate function I(w) (bottom) for the concentrated initial values x(0) = v(0) = 0, for stationary initial conditions and for generic initial conditions entailing both left and right linear tails.The figure summarises the ability of the model to exhibit no linear tail, one linear tail, or two linear tails.In order to confirm analytical results, figure 1(bottom) also compares the analytical rate function with rate functions estimated by means of numerical simulations of the original model.In fact, we simulated the trajectories of the discrete-time model (13) with time step ϵ = 10 −2 and fictitious diffusion coefficient D F set equal to 0. Then, we sampled the injected work W τ for different times τ = Nϵ as prescribed by equation (15).In the regions of w around D that we are able to explore numerically, we find a good agreement between analytic results and simulations.Recall that D is the typical value of W τ , with D = 1.0 (in arbitrary units) in figure 1.
Figure 2 provides some phase diagrams of the system as deduced by the ratios r − ≡ w − /D and r + ≡ D/w + between the rate function singularities w ± , where linear tails begin, and the mean work D. We have set r + ≡ 0 when the right linear tail does not occur, that is when M ≡ max{kσ 2 x , σ 2 v } ⩽ 4D/γ.In this way, we have r − ∈ (0, 1) and r + ∈ [0, 1), with r − ↑ 1 when w − approaches D from below, r + ↑ 1 when w + approaches D from above, and r + ↓ 0 when the region M ⩽ 4D/γ is approached from outside.Panels (a) and (b) of figure 2 report the phase diagrams for r − and r + , respectively, in the T − k plane at fixed γ, σ x and σ v with D = γk B T. We see that w − and w + approach D at low temperature T and/or large elastic constant k.At higher T and finite k, w − decreases to zero, whereas w + increases to infinity up to a critical value of T where the right linear tail ceases to exist.This critical value is identified by the condition 4k B T = max{kσ 2 x , σ 2 v }.Panels (c) and (d) of figure 2 report the phase diagrams for r − and r + in the σ x − σ v plane at fixed system parameters.We note that both w − and w + approach D at large σ x and/or large σ v .On the contrary, when σ x and σ v become small, the point w − decreases to zero, whereas w + ceases to exist as soon as max{kσ 2 x , σ 2 v } = 4D/γ.

Phenomenology of the trajectories
In this section we shed light on the physical mechanisms that originate singular trajectories, i.e. trajectories giving values w of the work in a linear tail of the RF I(w).We pursue the goal by analysing the trajectories of the model with system parameters as in figure 1(f), so that the rate function has both left and right linear tails.The simulation time τ we consider is much larger than the inertial time t I ≡ γ −1 .We fixed γ = 1.0,D = 1.0, k = 1.0, σ x = 10, σ v = 10, and τ = 50 ≫ t I = 1.0 (in arbitrary units).Inspection of singular trajectories reveals the presence of some big jump phenomena.Figure 3 reports three trajectories of the system, in terms of the position x(t) in figure 3(a) and An analogous behaviour characterises also the trajectories corresponding to w ≪ w − in the far left linear tail.Such big jumps are confirmed and characterised by figure 3(c), where the initial position and velocity are reported for a large pool of trajectories conditional on W τ = wτ with w ≪ w − (red), w ≫ w + (blue) and w − ⩽ w ⩽ w + (green).This panel clearly shows two different patterns for singular and non-singular trajectories, the former roughly localising on an annulus, so that at least one among x(0) and v(0) makes a big jump, the latter concentrating around the origin.At larger times τ , green points tend to accumulate more and more densely around the origin, whereas red and blue points remain essentially unaltered.This implies that big jumps in the initial conditions are of order √ τ .Therefore, we conclude that big jumps of order √ τ in the initial values of position and velocity are an essential ingredient for singular trajectories to occur.Coherently with what we find here, we note that in [23] the physical origin of the singularity in the rate function of the stationary, free particle was ascribed to a very energetic initial condition.Similarly, a big jump mechanism causing a particle to roll uphill in the potential was suggested in [32] to explain a singular rate function.
It is interesting to observe that the patterns for w ≪ w − and those for w ≫ w + in figure 3(c) overlap, meaning that same initial big jumps can produce trajectories associated to values of the work in both the left and the right linear tail of the rate function.In order to distinguish such trajectories in terms of the injected work, we discuss the action of the random force η(t).In fact, once a singular initial configuration is given, the linear tail where the value of the work falls is determined by the sign of the random force realisation, as demonstrated by figure 3(d).Specifically, consider a trajectory with w ≫ w + .This trajectory involves big jumps in the initial conditions and is characterised by a certain random force path.Figure 3(d) shows that if we now construct a trajectory with same initial conditions and opposite random force, then the corresponding value of the work falls in the left linear tail.Similarly, if we consider a trajectory with w ≪ w − and change the sign of the random force, then we end up with a value of the work in the right linear tail of the rate function.In conclusion, changing the sign of the random force systematically turns a trajectory with w ≫ w + into a trajectory with w ≪ w − , and vice versa.

Conclusions
In this paper we have studied the large deviations of the power injected by the random force of the thermal bath on a Brownian particle under the action of an external harmonic potential with generic uncorrelated Gaussian initial conditions.
In particular, we have computed analytically the SCGF and the rate function using the approach recently proposed in [38] for quadratic functionals of stable Gauss-Markov chains.The SCGF is found to be non-steep at the left boundary, except from the limiting case of concentrated initial values, i.e. x(0) = v(0) = 0.It is also non-steep at the right boundary when the variance of the initial position or initial velocity is large enough.The corresponding rate function therefore can exhibit no singularity and no linear tail, or one second-order singularity and one linear tail on the left, or two second-order singularities and two linear tails.Interestingly, the dependence of the rate function on the elastic constant k is found for large enough initial variances and lost for concentrated initial values and stationary initial condition, thus confirming the claim made in [23] for a confined Brownian particle.
In order to understand the physical mechanism originating the singularities, we have analysed the singular trajectories, i.e. the trajectories that give values of the work either in the left or in the right tail of the rate function.By resorting to numerical simulations, we have shown that singular trajectories are characterised by big jumps of order ∼ √ τ in the initial position and/or the initial velocity.Once the big jumps occur, the random force path determines if the value of the work falls within the left linear tail or the right linear tail.
A natural extension of this research will be to consider correlations in the initial condition or coloured thermal noise, such as the Ornstein-Uhlenbeck process [31,45] or the fractional Brownian motion [46].Leaving the field of Gaussian processes, another possible extension will be to consider an anharmonic potential, which could disrupt the singularities of the rate function, and continuous-time random walks, for which large deviation principles are known in great generality [47,48].On a more practical note, in analogy with experiments on Fluctuation Relations [49,50], we argue that our findings could be tested by means of optical trap experiments with µm diameter beads.

Figure 1 .
Figure 1.SCGF ϕ(µ) (top row) and rate function I(w) (bottom row) with system parameters γ = 1.0,D = 1.0, and k = 1.0 (in arbitrary units) for several initial conditions.(a) and (d): concentrated initial values x(0) = v(0) = 0 given by the limits σx, σv ↓ 0. (b) and (e): stationary system corresponding to σ 2 x = D/kγ and σ 2 v = D/γ.(c) and (f ): generic initial conditions defined by σx = 10 and σv = 10 (in arbitrary units).In the top panels, the blue line shows the SCGF in the effective domain, whereas the purple dashed curve in (b) and (c) depicts its trend in the primary domain (see sections 3.2 and 3.3 for details).In the bottom panels, the red solid line is the analytical rate function, whereas the coloured dotted lines are the numerical rate functions at different values of τ .In (e) and (f ), the red dashed lines mark the beginning of the left linear tail at w− and of the right linear tail at w+.

Figure 2 .
Figure 2. Phase diagrams as deduced by the ratios r− ≡ w−/D < 1 and r+ ≡ D/w+ < 1 between the rate function singularities w± and the mean value of the work D. Regarding r+, white regions delimited by black dashed lines are the regions of parameters where the right linear tail does not occur, that is where M ≡ max{kσ 2 x , σ 2 v } ⩽ 4D/γ.(a) and (b): r− and r+ for D = γk B T in the − k plane with γ = 1.0, k B = 1.0, σx = 10 and σv = 10 (in arbitrary units).(c) and (d): r− and r+ in the σx − σv plane with γ = 1.0,D = 1.0 and k = 1.0 (in arbitrary units).