{'fraction_non_alphanumeric': 0.05689892883656803, 'fraction_numerical': 0.015749284123448934, 'mean_word_length': 4.0927626248987306, 'pattern_counts': {'":': 0, '<': 21, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Work fluctuations in a time-dependent harmonic potential: rigorous results and beyond the overdamped limit

We investigate the stochastic motion of a Brownian particle in the harmonic potential with a time-dependent force constant. It may describe the motion of a colloidal particle in an optical trap where the potential well is formed by a time-dependent field. We use the path integral formalism to solve the Langevin equation and the associated Fokker-Planck (Kramers) equation. Rigorous relations are derived to generate the probability density function for the time-dependent nonequilibrium work production beyond the overdamped limit. We find that the work distribution exhibits an exponential tail with a power-law prefactor, accompanied by an interesting oscillatory feature (multiple pseudo locking-unlocking transitions) due to the inertial effect. Some exactly solvable cases are also discussed in the overdamped limit.


I. INTRODUCTION
Nonequilibrium (NEQ) fluctuation has been an important issue in the field of statistical mechanics for the last two decades since the first discovery of the fluctuation theorem for the entropy production [1][2][3]. Fluctuation theorems (FT's) [4][5][6][7][8][9][10][11] and the Jarzynski equality (JE) [12,13] are the central theoretical relations governing NEQ fluctuation phenomena, widely valid for many NEQ systems, deterministic or stochastic, thermostatted with heat baths. It basically deals with a large fluctuation around the average measurement in theory and experiment with considerable contribution from rare events. Such phenomena become dominant for a system with small degrees of freedom. There have been extensive studies of small experimental systems such as a microscopic bead dragging in a viscous fluid [14], a single molecule of the RNA under mechanical stretch [15,16], an oscillating bead under the translating center of the optical trap [17], the circuit of an electric dipole in electric potential bias [18], an ultra light metallic wire under torsion [19], etc.
External bias was considered as a typical underlying mechanism for NEQ systems such as the nano circuit device with potential bias [18,20], the harmonic oscillator under constant torque applied [19,21,22], one dimensional lattice gas in contact at boundaries with different heat or particle baths [23,24]. A nonconservative force was also recognized as a source for the entropy production [6] such as in a nano heat engine in contact with multiple reservoirs for a circulating current in high dimensional systems [25][26][27]. Non-Markovian nature caused by memory effect or colored noise is another source for NEQ [28,29]. In these examples, the system reaches a NEQ steady state (NESS) after a transient period, where a persistent nonzero current, directed or circulated, generates the incessant work production. The probability distribution function (PDF) of the work production exhibits an exponential decay with a power-law prefactor in the rare-event region [27], along with interesting unusual features such as initial condition dependency of the large deviation function [30][31][32][33] and multiple dynamic transitions in reaching the NESS [27,34].
On the other hand, a time-dependent perturbation on external parameters such as electric field, magnetic field, volume, force constant, etc. generates a genuine time-dependent NEQ state where the system never maintains the NESS. For example, the stochastic motion of a Brownian particle was studied in the harmonic potential moving with a constant velocity (sliding parabola potential) [35][36][37][38][39][40], and also in the harmonic potential with a time-dependent force constant (breathing parabola potential) [41][42][43][44]. In these cases, the work PDF also shows an exponential decay with a power-law prefactor in the rare-event region, along with a time-dependent characteristic value for the work production determining the exponential decay shape. Most of previous studies considered the overdamped limit, partly because the experimental situation for a colloidal particle in a harmonic trap can be well approximated in the overdamped limit and also partly because the analytical treatment is much easier [35][36][37][38][39][40][41][42][43][44].
In this paper, we generalize these results beyond the overdamped limit (underdamped case) for a Brownian particle in a breathing parabola potential with the momentum variable kept intact. We focus on the inertial effect on the time-dependent characteristic value for the work production. Our model may also serve as a soft-wall version of the box expansion or compression with a single Brownian particle inside, in contact with a thermal reservoir [45]. The experimental setup is also feasible: In a molecular tweezer or an optical trap experiment, the potential well can be approximated by the harmonic potential. The shape of the potential well, so the force constant of the harmonic potential is set to vary with a time-dependent external field. The stochastic motion is described by the Langevin equation and the corresponding Fokker-Planck (Kramers) equation. We use the path integral formalism to derive rigorous relations from which the time-dependent work PDF can be easily calculated with machine accuracy and also its cumulants at all temperatures with any time-dependent force constant. We find an interesting oscillatory feature of the work PDF shape, solely due to the inertial effect (absent in the overdamped limit), which resembles multiple dynamic transitions found in the linear diffusion system [27,34], but shows smooth crossovers rather than sharp transitions. Thus we call this crossover as a pseudo locking-unlocking transition. The existence of multiple pseudo dynamic transitions may be related to the existence of the phase-space circulating current in the underdamped case, but full intuitive understanding calls for further investigation in future.
In section II, we introduce the breathing harmonic potential function and discuss the FT's. In section III, we derive the equations for the PDF and the cumulants for the work production using the path integral formalism. Our formalism is tested for the systems in the sudden change limit. In section IV, we present the analysis for the work PDF and find the exponential tail with a powerlaw prefactor. In section V, we study the overdamped limit for exactly solvable cases. In section VI, we discuss the main results of our work and the perspective to future works.

II. TIME-DEPENDENT HARMONIC POTENTIAL
We consider the Brownian motion of a particle in one dimension under the breathing harmonic potential with a time-dependent force constant k = k(t) and in contact with a heat bath. The equations of motion are given bẏ where γ is a damping coefficient and ξ is white noise with zero mean satisfying ξ(t)ξ(t ′ ) = 2d δ(t − t ′ ). The diffusion coefficient d is chosen to satisfy the Einstein relation, d = β −1 γ, which guarantees the equilibrium (EQ) Boltzmann distribution at inverse temperature β in the steady state, if k is constant in time.
The equations of motion can be rewritten aṡ where q ≡ (x, p) T and ζ ≡ (0, ξ) T . Here, the superscript T denotes the transpose of a vector or a matrix. The force matrix F is given by The energy of a particle is given by E(q; k) = p 2 2m + kx 2 2 , which is written as E = 1 2 q T · H · q with a Hamiltonian matrix Let P (q, t) be the probability density function for finding a particle at state q and at time t. Then it satisfies the Fokker-Planck equation, specially called the Kramers equation, where ∇ = (∂ x , ∂ p ) T and the diffusion matrix is given by where a small positive parameter ǫ is introduced for convenience, making possible the inversion of the diffusion matrix D during a formal manipulation in the path integral formulation. In the end, we take the ǫ → 0 limit to recover the delta function constraint δ(ẋ − p/m) for position and momentum. With ǫ, position and momentum can be treated on the same footing, which gives us the formal advantage over the usual path integral with the δ-function constraint. This approach works well. For instance, one can reproduce the expected results for the EQ PDF when k is a time-independent constant [46]. In this case, the EQ Boltzmann distribution becomes the stationary solution of the Kramers equation in the limit ǫ → 0. The partition function is given by Z(k) = dq e −βE(q;k) = (4π 2 m/(β 2 k)) 1/2 , so that the free energy is given by F (k) = − 1 2β ln(4π 2 m/(β 2 k)). When the force constant k varies in time, the system is driven into a NEQ state. It belongs to the Jarzynski's criterion for NEQ, where the rate of work production is given byẆ =k(∂E/∂k). Then the NEQ work W done on the particle moving along a path q(τ ) for 0 < τ < t is given by This is rewritten in a matrix form as where The system is assumed to be initially in EQ at β with k(0) = k i and will reach a final state with k(t) = k f , which is certainly far from EQ. In this situation, the JE states where · · · denotes the average over all possible paths q(τ ) and ∆F is the free energy difference F (k f ) − F (k i ) at β. The JE can be trivially derived from the Crooks relation [5] where P F (W ) = δ(W − W[q]) F is the PDF for the work production W during the forward process with the change from k i to k f and vice versa P R (W) for the reverse process. The JE and the Crooks relation can be proved for a general form of energy and external perturbation for the Langevin dynamics, if the system is initially in EQ at β [11]. However, the explicit expression for P (W ) is not generally known. There are not many stochastic models which can be solved analytically for the PDF of fluctuating quantities. The breathing harmonic potential with a time-dependent force constant is not only analytically tractable, but can also serve an appropriate model for the potential well in an optical tweezer or trap.

III. PATH INTEGRAL FORMALISM WITH TIME-DEPENDENT FORCE
The Fokker-Planck equation for a multivariate system with a linear drift force, known as the high dimensional Ornstein-Uhlenbeck process, is solvable, i.e., the time-dependent PDF P (q, t) can be obtained analytically [47,48]. NEQ properties for this process were investigated in detail from the view of the circulating NESS current [25] and the violation of the fluctuationdissipation relation [49]. Recently we have revisited this system in the light of the fluctuation theorem [27] in the case that the (non-conservative) drift force does not vary with time. The path integral formalism developed in that study can be extended to the present problem with timedependent drift represented by the force matrix To describe the NEQ fluctuations, it is convenient to introduce a path integral during time period t as The initial PDF for q 0 is chosen to follow the EQ Boltzmann distribution P eq (q 0 ; k(0)) as in Eq. (7), and D[q](· · ·) denotes the integration over all possible paths connecting q(0) = q 0 and q(t) = q 1 for 0 < τ < t. The Lagrangian L is chosen to read The source term ( dτ l T · q) is introduced for a later use. Note that the exponent of the integrand is at most quadratic in q. Hence the path integration can be computed exactly by Gaussian integrations. The quantity I is useful in calculating physical quantities of interest. For example, the PDF P (q, t) is given by [50] P (q, t) = I(q, λ = 0; l(τ ) = 0) .
The PDF for the NEQ work production can be also calculated from I. First, we define a dimensionless quantity for the work as w = βW for simplicity and introduce its generating function which can be obtained as Note that the JE, G(1) = exp[−β∆F ], can be proven explicitly in this path integral formalism as well as the generalized Crooks relation as process. The PDF for the dimensionless work w is then obtained by the inverse Fourier transformation as For an arbitrary functional A[q(τ )], one can also calculate its ensemble-averaged value from I. Defining the cumulant generating function as one finds that We will use this relation to calculate the cumulants of the work. The path integral, Eq. (13), can be evaluated by using the methods developed in our recent study [27]. Here, we will present the results without showing detailed calculation steps.

A. Probability distribution function
The PDF P (q, t) is given by where the kernel A(t) is a symmetric matrix, satisfying the differential equation as The formal solution is given by with the initial condition A(0) = βH(0). Here the evolution operator U is given by where the subscript denotes the time-ordered product, and satisfies the differential equation with U(t ′ ; t ′ ) = I (the identity matrix).
In the absence of noises (D = 0), U(t; t ′ ) describes the deterministic evolution by q(t) = U(t; t ′ ) q(t ′ ). When the force matrix is constant in time so that U(t, t ′ ) = e −(t−t ′ )F , one can do the integral in Eq. (23) and find the explicit solution for A(t) (see Eq. (21) in [27]). For a general time-dependent F(t), it is difficult to treat U analytically. However, Eqs. (23) and (25) can be solved very precisely by numerical integrations.

B. Work distribution function
The generating function for the work distribution in Eq. (17) involves the integration of the quantity I with nonzero λ. The work W[q] coupled to λ is also quadratic in q (see Eq. (9)), hence the integration can be performed in the same way as was done for the probability distribution. After some algebra, one can derive where Λ = βḢ in Eq. (10) andÃ(τ ; λ) is the modified kernel due to −λβW in Eq. (13). It is found to satisfy the nonlinear differential equation where the initial condition is given byÃ(0; λ) = βH(0). This nonlinear differential equation can be solved easily for λ = 0 and 1. The solution isÃ(τ ; 0) = A(τ ) in Eq. (23), whileÃ(τ ; 1) = βH(τ ). Interestingly,Ã(τ ; 1) corresponds to the kernel for P (q, τ ) in the quasi-static process. Inserting this into Eq. (26), we find which verifies the JE.
It is not efficient to integrate Eq. (27) numerically, be-causeÃ sometimes becomes singular, as τ increases in some range of λ. Then,Ã −1 cannot be defined any more. Therefore it is more convenient to rewrite Eq. (27) in terms ofÃ(τ ; λ) as along with an equivalent and more effective expression for the generating function replacing Eq. (26) as (30) A similar result has been found in the time-independent case [27]. Equations (29) and (30) are ingredients for numerical study of the work production distribution P (w). Especially, the exponentially decaying tail behavior of P (w) is manifested by the divergence of G(λ), which turned out to be fully captured by the singularity in the logarithmic boundary term in Eq. (30). Thus we will focus on the behavior of detÃ(t; λ) in the next section.
One can observe that Eq. (29) becomes independent of β ifÃ is scaled by β. This proves that G(λ) as well as P (w) is independent of β. Therefore, P (W ) is simply equal to βP (w) with w = βW . In the weak noise (large β) limit [41,42], the tail behavior of P (w) for large |w| determines exactly and fully the work distribution P (W ) except for a narrow central region, |W | < β −1 . We will come back to this issue later.

C. Cumulants of work production
The cumulant generating function in Eq. (19) is found as This form is expected because the Lagrangian is quadratic in q and the source field l(τ ) is linearly coupled to q. The kernel Γ(τ, τ ′ ) is given as and Γ(τ, τ ′ ) = Γ T (τ ′ , τ ). Using Eqs. (20) and (31), one can express any functional of path q(τ ) in terms of Γ(τ, τ ′ ). For example, the first and the second cumulant of the work are given by where W 2 c = W 2 − W 2 . Note thatḢ 11 =k anḋ H ab = 0, otherwise. Then, the expressions become simpler: One can also find higher-order cumulants in terms of Γ(τ, τ ′ ), which are non-zero in all orders. This implies that the PDF P (w) should have a non-Gaussian form. The PDF shape will be discussed further in the next session.

D. Sudden change limit
A sudden change is a rare case where one can calculate the work PDF exactly, even in the underdamped case. Suppose that the particle is in EQ under the harmonic potential with the force constant k i , and the force constant is changed abruptly to k f at time t = 0 [41,42]. If the particle is in a state q = (x, p) T just before the change (t = 0 − ), its state still remains the same right after the change (t = 0 + ) as well as the PDF P (q, t). Only change occurs in the potential energy, which results in the Then the work production W(q) = ∆E.
As the initial distribution is given by the EQ Boltzmann distribution, the PDF P (x) = dp P eq (q; k i ) = βk i /(2π)e −βkix 2 /2 . Then, the PDF P (w) of the dimensionless work w = βW can be easily derived, using P (w)dw = 2P (x)|dx|, which yields that with a = k i /(k f − k i ) and θ(w) is the Heaviside step function.
The generating function G can be easily calculated, using Eq. (16), as The JE is also seen from G(1) = (k f /k i ) −1/2 . Our analytic formalism in the previous subsections can also reproduce G(λ). The sudden change in the potential function can be studied by considering Integrating Eq. (29) from τ = 0 − to τ = 0 + , one get whereÃ(0 − ; λ) = βH(0 − ) is used. Then, from Eq. (30), one can easily reproduce the result in Eq. (36). The cumulants of the work can be also easily calculated as

IV. ANALYSIS OF WORK DISTRIBUTION
The analytic formalism developed in this paper is very useful to investigate the work distribution P (w) numerically, in particular, its tail behavior, in contrast to direct numerical integration of the equations of motion where we always face with a statistics problem, becoming serious in rare-event regions. In this session, we first present numerical data from the latter method to check fluctuation relations and get some insights on the nature of the work production distribution. Then, the tail behavior of P (w) is carefully examined by the former method.
It is convenient to work with dimensionless variables in numerical calculations. Rescaling γ , one finds that the dimensionless variablesx andp satisfy the same equations of motion as in Eq. (1) with respect to the dimensionless timet with dimensionless parameters m =γ =β = 1 andk = ( m γ 2 )k. Hence, without loss of generality, we will set m = γ = β = 1. The only independent parameter is the force constant k. We will consider a special case wherek is a time-independent constant, i.e., k(t) = k i (1 + αt) for convenience.
First, we check the JE and the Crooks relation from direct numerical integration of the time-discretized equations of motion. We adopt a notation X n = X(t = t n ) for a time-dependent quantity X(t) where t n = n∆t (n = 0, 1, 2, · · ·) are discretized times in unit of ∆t. Then, the equations of motion are solved from the difference equations where η n are independent Gaussian-distributed random variables with zero mean and unit variance. An initial configuration q 0 = (x 0 , p 0 ) is drawn from the EQ distribution of Eq. (7). The dimensionless NEQ work production w n = βW n up to time t n is evaluated from the recursion relation with w 0 = 0. Repeating the simulations N S times, one can measure the PDF P (w) and the generating function G(λ) numerically. Note that the work production w is always positive for α > 0 and negative for α < 0, independent of noise realizations.
In simulations, we take ∆t = 10 −3 and N S = 10 7 . The force constant k(t) is taken to vary linearly from k i = 1 to k f = 4 and k i = 4 to k f = 1, which will be referred to as a forward (F) and a reverse (R) process, respectively. Figure 1(a) shows the P F (w) for the F process till t = 3 with α = 1 and P R (w) for the R process with α = −1/4. We compare P F (w) and e w−β∆F P R (−w) with β∆F = 1 2 ln k f ki = ln 2. They seem to overlap each other well (except for the region with very small P (w)), which supports the validity of the Crooks relation in Eq. (12).
In order to examine the PDF in detail, we compute the generating function G(λ) = e −λw . These are plotted in Fig. 1(b). The JE states that G(λ = 1) = e −β∆F where β∆F = ln 2 for the F process and − ln 2 for the R process. Indeed, the numerical curves pass through the JE points. We also compare G F (λ) and e −β∆F G R (1 − λ), to check the generalized Crooks relation. For moderate values of λ, both data align along a single curve. However, there is a slight discrepancy for large |λ| where rare fluctuations with large values of |w| are important. This reflects a statistical uncertainty due to limited samplings. Even with N S = 10 7 samples, statistics is poor for those rare fluctuations. Now, we utilize the analytic results in Eqs. (29) and (30) which are free from statistical errors, in order to determine the tail part of P (w) precisely. In discretized times in unit of ∆t = 10 −3 , the nonlinear differential equation (29) forÃ(τ ; λ) is solved with the initial condi-tionÃ(0; λ) = βH(0) and the integration in Eq. (30) is performed numerically. We present the numerical results for the F and R processes (dashed lines) in Fig. 1(b). As expected, the previous simulation results deviate significantly from our new improved numerical results in rare- In fact, our numerical data in Fig. 1(b) shows that G(λ) is divergent at a threshold λ 0 (F) and 1 − λ 0 (R) with λ 0 ≃ −0.84713 < 0. The divergence occurs when detÃ(t; λ) = 0 as seen in Eq. (30). Figure 2 shows the time evolution of detÃ(t; λ) at several values of λ in the case with k i = 1 and α = 1. To a given value of t, detÃ becomes smaller as λ decreases and vanishes at a threshold λ 0 . One can solve the equation detÃ(t; λ) = 0 numerically to obtain the t-dependent threshold λ 0 . Figure 3 shows the numerical results for the system with k i = 1 and α = 0.5, 1, and 2. The threshold depends on k i and α, and increases monotonically and converges to a finite limiting value λ ∞ 0 ≃ −1.39162, −0.81311, and −0.48110 in the t → ∞ limit.
The singular behavior of G(λ) reveals the asymptotic behavior of the tail shape of P (w) for large |w|. Due to the generalized Crooks relation, it suffices to consider the F process with positive α (compression). Figure 2 suggests that detÃ(t; λ) is regular near λ = λ 0 (t) so that one can write as detÃ ≃ c(λ − λ 0 (t)) with a positive constant c. Then, from Eq. (30), G(λ) diverges as and its square-root singularity indicates that [27] P (w) ∼ w −1/2 e −|λ0(t)|w (42) for large positive w with the characteristic work w 0 ≡ 1/|λ 0 |. Our analytic formalism, combined with the numerical analysis, predicts that there is the exponential tail in P (w) with the power-law prefactor with the exponent −1/2. Note that the abrupt change of k (sudden change limit) also yields the same tail, as was shown in the preceding section. We test the tail shape by direct numerical integration of the equations of motion, using    Fig. 4, the PDF P (w) multiplied with e |λ0|w follows a power-law scaling for large w, which confirms the tail shape. Huge fluctuations for large w are due to statistical errors in sampling rare events.
The tail shape of P (w) in Eq. (42) is consistent with previous findings in the overdamped limit [41][42][43][44]. It is also not surprising to find that |λ 0 (t)| decreases with t because one can easily expect that the work PDF should be more distributed (flatter) as t increases. However, the monotonic behavior of |λ 0 (t)| is not trivially smooth, but has an interesting repeating structure.
In Fig. 3, we observe a stepwise change of |λ 0 (t)| in time, composed of a rather fast linear change followed by quite slow plateau-type change, which repeats itself but with decreasing size both in magnitude and time period, and finally converges to the limiting value of |λ ∞ 0 |. This implies that the exponential tail of P (w) relaxes into the limiting distribution via multiple (possibly infinitely many) fast-slow-type relaxation dynamics. These repeated fast-slow-type dynamics resemble multiple locking-unlocking dynamic transitions found in twodimensional linear diffusion systems in the overdamped limit [34]. However, our case shows rather smooth crossovers between fast and slow dynamics, in contrast to sharp transitions with completely flat plateaus in |λ 0 (t)| [34]. We call the stepwise changes in our case as pseudo locking-unlocking transitions. In the mathematical language, we can not find any detÃ(t; λ) curve tangential to the t axis (detÃ = 0) in Fig. 2, which prohibits a completely flat plateau, so no sharp transition is realized.
It is easy to recognize that the oscillatory feature of detÃ in Fig. 2 evokes the stepwise change of |λ 0 (t)|. First, note that all detÃ curves show oscillatory wiggles almost simultaneously in time and the oscillation frequency grows as t increases. So we can define a set of characteristic times (t ± 1 , t ± 2 , · · ·) where all curves show a local minimum (+) or maximum (−) simultaneously, at least, approximately. The oscillatory behavior is related to the increasing frequency of the harmonic oscillator caused by the increasing force constant k(t) = k i (1 + αt).
Due to this oscillatory feature of the detÃ curves, one can easily figure out that the curves cross the t axis sparsely right after t = t + 1 until t = t − 1 , and densely for t − 1 < t < t + 2 , and so on. Therefore, λ 0 increases very fast during 0 < t < t + 1 and very slow during t + 1 < t < t − 1 , and this fast-slow relaxation dynamics repeats itself with increasing frequency.
The underlying mechanism of these pseudo lockingunlocking transitions should be similar to one for the sharp transitions found in the two-dimensional linear diffusion systems [34]. Only differences are the nature of the rotational current, which exists here only in the phase space of (x, p) and the time-dependent external force, which acts a role of the rotational driving force as well as the (time-dependent) anisotropic harmonic potential in the phase space. However, we could not find a sharp dynamic transition in our model with an arbitrary choice of parameters (k i , α). Recalling what we learned in [34], we guess that the anisotropy may be always small in our model, compared to the driving force magnitude, in order to avoid a sharp transition. For full understanding, however, a further investigation is necessary.
For the overdamped one-dimensional case, we cannot have any rotational current, so the oscillatory behavior is expected to be absent at all, which is confirmed rigorously in the next session. Therefore, it can be concluded that the pseudo locking-unlocking transitions found in the underdamped case originate from the existence of the rotational current in the phase space.

V. OVERDAMPED LIMIT
In the overdamped limit, the usual Fokker-Planck equation, replacing the Kramers equation, reads Then one can use the analytic formalism developed so far for the Brownian dynamics by replacing D with (γβ) −1 , F with γ −1 k, H with k, andḢ withk. Then, the work generating function G(λ) is given by where the scalar quantityÃ(t; λ) satisfies a nonlinear differential equation with the initial conditionÃ(0; λ) = βk(0) = βk i . We will set γ = β = 1 without loss of generality. The overdampled limit is investigated in detail with different choices of k(τ ).
A. k(τ ) = ki(1 + ατ ) All the relevant informations are obtained from Eqs. (44) and (45). Unfortunately, the closed-form solution forÃ(t; λ) or G(λ) is not available. However, the highly accurate numerical solution is possible, which is shown in Fig. 5(a) forÃ(t; λ) with k i = 1 and α = 1. As in the Brownian dynamics, it becomes zero at a tdependent threshold λ 0 . The threshold is plotted in Fig. 5(b). ¿From the similar analysis for Eq. (41), the PDF P (w) can be found to have the same tail shape as in Eq. (42). Note that the particle dynamics does not display any oscillatory motion in the overdamped limit, henceÃ does not either as seen in Fig 5(a). Thus, the threshold λ 0 in Fig. 5(b) varies in time smoothly, showing no stepwise change at all.
Engel and Nickelsen studied the same harmonic potential problem in the overdamped limit [41,42]. In their studies, they evaluated the path integral using the saddle point method in the low noise (large β) limit. However, as pointed out in section III, P (w) is independent of β, so P (W ) with w = βW can be exactly determined from the tail behavior of P (w) except for |W | < β −1 .
In order to calculate the cumulants of the work production, we need to evaluate the PDF kernel A(t), first, using Eqs. (23) and (24). As we do not need to deal with the time ordered product, we simply gets U (t, t ′ ) = exp[− t t ′ dτ k(τ )]. Then, one can explicitly calculate the cumulants for the work production. Using Eqs. (23) and (32), we find and where we use and where The first and second cumulants are found as where h(τ, τ ′ ) = e −2ki(τ +τ ′ )−kiα(τ 2 +τ ′2 ) .
There are two extreme cases; quasi-static and sudden processes. For the quasi-static process, one can take the limit: α → 0, t → ∞ with a finite value of αt = (k f − k i )/k i . Changing variable as u = k i ατ /(k f − k i ), one can get an approximate result for the time integral. The important ingredient for the integration is b a du e cu 2 → 1 2c which agree with the results by Speck [43]. It indicates that the work distribution function is perfectly Gaussian centered around w = β∆F up to O(α), and the non-Gaussianity starts to appear in O(α 2 ) [51]. In the quasistatic limit (α → ∞), the work distribution function becomes a delta function, as expected for EQ processes. For a sudden process, one can take the opposite limit; α → ∞, t → 0 with a finite value of αt = (k f − k i )/k i . Also using the same variable u, the integrand of b a du e cu 2 can be expanded in orders of c in the c → 0 limit. As a result, one can get Note that w and w 2 c are finite even for an instantaneous change (α = ∞), which agrees with the sudden change limit for the underdamped case in Eq. (39). It is different from the case for the rigid wall moving with speed v in the v → ∞ limit, where we expect w → 0 [45]. The difference is due to the distinctive situations. For the former, the collision occurs everywhere with the harmonic potential, while for the latter the collision occurs only at the descending wall. The similarity lies in the non-trivial fluctuation around the average value.
With this specific form, one can find the closed-form solution for the work generating function G(λ). For α > 0, the harmonic potential becomes flatter with time τ ≥ 0 and the work w done on the particle is always negative. On the other hand, for α < 0, the harmonic potential becomes stiffer with time τ (0 ≤ τ < 1/|α|) and w is always positive.
It is convenient to change variables as The time-like variable u is monotonically increasing with τ , starting from 0 to ∞ for any nonzero α. ¿From Eq. (45), one obtains a differential equation for f λ (u): Note that κ may be either positive real or pure imaginary depending on the range ofλ. In either case, the solution is given by with cos(ix) = cosh x and sin(ix) = i sinh x for any x.
With the solution for f λ (u) or equivalently forÃ(τ ; λ), one can obtain the work generating function using Eq. (44). It is useful to note that f λ (u) = c + 1 2 d du ln[cos(2κu) + (k i − c)sin(2κu)/κ]. After a straightforward algebra, we find that with u t = 1 α ln(1 + αt). The work PDF P (w) can be obtained by the inverse Fourier transformation of G(λ) in Eq. (18). Note that the generating function has an inverse square-root singularity at a particular value of λ = λ 0 (u) at which the denominator in Eq. (60) vanishes. In fact, the singularity occurs whenÃ(τ ; λ) = 0 or equivalently f λ (u) = 0, as seen in Eq. (44). The inverse square-root singularity in G(λ) at λ = λ 0 implies an exponential tail with a power-law prefactor [27] in P (w) in the ω → −∞ limit for α > 0 (λ 0 > 0) or in the ω → ∞ limit for α < 0 (λ 0 < 0). The singular point satisfies the relation withλ 0 = αk i λ 0 /2. It should be understood that tan −1 (ix) = i tanh −1 (x) and that 0 ≤ tan −1 x < π for a real x. Figure 6 shows the plots for the solution of Eq. (62), where the divergence of u t is observed as λ 0 approaches the limiting value from above. Interestingly, the time dependence and the limiting value are very different depending on whether c < 0 (α ≤ −2k i ), Especially, the limiting value λ ∞ 0 = lim u→∞ λ 0 is given by We present the plot of the λ ∞ 0 as a function of α/k i in Fig. 7.
There is an interesting symmetry of λ ∞ 0 (−α) = 1 − λ ∞ 0 (α). This comes from the Crooks relation. The reverse protocol with respect to the forward protocol k(τ ) = k i /(1 + ατ ) should be given as k r (τ ) = k(t − τ ) = k f /(1 + α r τ ) with α r = −αk f /k i and k f = k i /(1 + αt). If the system starts with the EQ distribution with k r (0) = k f , all results derived here can be also applied for the reverse process by replacing k i by k f and α by −αk f /k i . Then, Eq. (63) gives us λ ∞ 0,r (α r ) = λ ∞ 0 (−α). The Crooks relation of Eq. (12) yields λ ∞ 0 = 1 − λ ∞ 0,r in the large w limit, which leads to our symmetry of λ ∞ 0 (−α) = 1 − λ ∞ 0 (α). We add a few remarks on the interesting α dependence of λ ∞ 0 : (i) For α ≥ 2k i , the tail shape of P (w) does not change with α with fixed λ ∞ 0 = 1. When α is large enough, the harmonic potential flattens very fast. Then, the particle dynamics starting from the EQ distribution with k i would be rather localized and not fully relaxed into a flattened harmonic potential. So the fluctuation in w may be dominated by an initial transient behavior even in the long-time limit (t → ∞), independent of the detailed shape of k(t). The sudden change limit discussed in subsection III D corresponds to the α = ∞ limit with (35) is consistent with the result for α ≥ 2k i . Nevertheless, it is still quite remarkable to find λ ∞ 0 = 1 for large but finite α. Similar features of initial-distribution dominance in the large deviation function in the long-time limit have been found in various different situations [30][31][32][33].
(ii) For |α| < 2k i , |λ ∞ | decreases monotonically as |α| increases. This behavior is compatible with the common wisdom that the fluctuation gets stronger (longer tail in P (ω)) as the rate of the change in driving increases.
(iii) When α ≤ −2k i (or c ≤ 0), we obtain that λ ∞ 0 = 0. This implies that P (ω) has a pure power-law tail in the positive-w region in the u → ∞ (t → 1/|α|) limit. In this case, the driving is strong enough to generate huge fluctuations.
The generating function also produces the cumulants of the work production by w n c = d n ln G/d(−λ) n | λ=0 . We focus on the mean value of the work, which is given by The quasi-static process corresponds to the limiting case where α → 0, t → ∞ with fixed αt = (k i − k f )/k f . In this limit, we find which agrees with Eq. (54). For a sudden process, we take the opposite limit where α → ∞, t → 0 with the same fixed value of αt in the above. Eq. (65) approaches w = (k f − k i )/2k i , which agrees with Eq. (55).

VI. DISCUSSION
The Brownian dynamics with mixture of position and momentum variables, having the even and odd parity respectively in time reversal, has not been fully scrutinized. In many literatures, the overdamped limit was taken for simplicity or as regards a light mass in experiment; otherwise the Kramers equation should be investigated, which is a nontrivial task. Putting the position and momentum on the same footing and introducing the singular diffusion matrix as in Eq. (6), we converted it into the usual Fokker-Planck equation, where the strict constraint on the position and momentum as δ(ẋ − p/m) is relaxed, so it is easier to handle. It is a well-known fact from the classical text books, but there are not many examples exploiting this method. In our study, we have proven this approach to be successful in finding the results analytically and numerically beyond the overdamped limit.
For the motion under the harmonic potential with a time-dependent force constant k(t), we have succeeded in examining the PDF of the work production rigorously, the most important quantity in NEQ fluctuations. As a result, we have found the exponential tail with a power-law prefactor in the PDF P (w) and |λ 0 (t)|, the characteristic constant in the exponential tail, to decrease with time t showing a fine structure of infinite but not sharply-edged staircase. By comparing the multiple locking-unlocking transitions (sharply-edged stair-case) found in the two-dimensional linear diffusion system [34], we call these rather smooth staircase as a manifestation of multiple pseudo locking-unlocking transitions. These pseudo transitions completely go away in the overdamped limit where no rotational current exists even in the phase space. We also consider some exactly solvable models in the overdamped limit and found an interesting power-law (not exponential) tail in P (w) for the case of rather fast compression (α ≤ −2k i ) with the protocol k(t) = k i /(1 + αt), which implies huge NEQ fluctuations.
The potential well in optical tweezers or an optical trap experiment is controlled by an external field, so can have the shape change in time due to a time-varying external field. Therefore, our study can serve a theoretical basis for such experiments, since the potential well may be approximated to be harmonic. The perturbation theory might be exploited to investigate an anharmonic effect. Our recent study of the multi-dimensional diffusion dynamics for a linear drift force [27,34] in the overdamped limit can also be realized in such experiments. We suggest many interesting experiments to be carried out in this direction.