Energetics of a driven Brownian harmonic oscillator
Abstract. We provide insights into energetics of a Brownian oscillator in contact with a heat bath and driven by an external unbiased time-periodic force that takes the system out of thermodynamic equilibrium. Solving the corresponding Langevin equation, we compute average kinetic and potential energies in the long-time stationary state. We also derive the energy balance equation and study the energy flow in the system. In particular, we identify the energy delivered by the external force, the energy dissipated by a thermal bath and the energy provided by thermal equilibrium fluctuations. Next, we illustrate Jarzynski work-fluctuation relation and consider the stationary state fluctuation theorem for the total work done on the system by the external force. Finally, by determining time scales in the system, we analyze the strong damping regime and discuss the problem of overdamped dynamics when inertial effects can be neglected.
pacs:05.10.Gg, 05.40.-a, 05.40.Ca, 05.60.-k,
Classical thermodynamics and statistical physics usually deals with systems either isolated or at equilibrium with their environment. For example, if a system is in contact and weakly interacts with a thermostat (heat bath), its stationary state is a Gibbs canonical equilibrium state . However, this “ideal” situation is not always satisfied and real systems often have strong couplings with their environments or are driven far from thermodynamic equilibrium by an external agent. Recently, a lot of theoretical developments have improved our understanding of non-equilibrium thermodynamics such as fluctuation theorems [2, 3, 4, 5, 6], which have been generalized both in classical and quantum statistical physics [2, 7, 8, 9, 10]. Moreover, a number of important results have been obtained in the theory of stochastic thermodynamics [11, 12, 13]. These approaches are particularly relevant for small systems where thermal fluctuations cannot be neglected, such as colloidal particles, polymers, molecular motors, etc. . Far from thermal equilibrium systems exhibit many fascinating and non-intuitive properties like absolute negative mobility , stochastic resonance , non-monotonic temperature dependence of a diffusion coefficient  or transient but long-lasting anomalous diffusion [18, 19]. Interestingly, in this recent framework, the Langevin equation, first introduced in 1908 to describe the Brownian motion , has become a paradigmatic study case.
The intention of this paper is to present in a pedagogical way an example of a non-equilibrium system. It should be useful for students in statistical physics to follow step by step a derivation of some formulas and next to perform detailed analysis of the system properties, in particular its energetics. We consider only one model: the classical harmonic oscillator in contact with a heat bath and driven by an external, unbiased time-periodic force. The system is described by a Langevin equation. It is simple enough so that the main statistical quantities characterizing its dynamics can be calculated exactly in analytical forms. Moreover, it is a relevant model for a lot of experimental applications at small scales. For example, optical tweezers used to trap and manipulate micrometric objects are very well described by a harmonic potential for small displacements [21, 22]. Nearly forty years after their invention , optical tweezers have become a very common tool and are now widely used in biology and physics laboratories [24, 25, 26]. Another classical example of applications is a small torsion pendulum [27, 28, 29]. It has also been shown in a one dimensional system of hard particles that the hard core interaction among particles can lead to a restoring force that can be modelled by a harmonic potential . Moreover, a two-dimensional diffusion of a Brownian particle in a periodic asymmetric channel with soft walls can been modeled by a parabolic potential which mimics a restoring force [31, 32]. Harmonic potential has been also used in modelling confined colloidal systems . More generally, any system that undergoes thermal fluctuations close enough to its stable equilibrium position can be described as a Brownian harmonic oscillator.
In this article, we have chosen to consider the case of an external sinusoidal force acting on the oscillator and we focus on interpretation of the oscillator energetics. Namely, we analyze the long time dynamics and compute the potential and kinetic energies, and we describe the energy balance between the system, the thermal bath and the external force. We also discuss in details the regime of strong damping (frequently named overdamped) dynamics where the inertia term can be neglected. We show that the overdamped regime can be rigorously analyzed by revealing characteristic time scales in the system and then by transforming the Langevin equation to its dimensionless form. The similar procedure can be applied to energetics of the system.
The article is organized as follows: In section 2 we describe the model of a periodically driven Brownian particle in a harmonic potential and we present the Langevin equation that governs its dynamics. We also define three kinds of averages of dynamical variables of the oscillator. In section 3, we present a solution of the Langevin equation for the position and the velocity of the particle. Next, in section 4, we compute the averaged kinetic and potential energies of the system in the long-time regime. In section 5 we derive equations for the energy balance and we discuss the energy exchanges between the system, the heat bath, and the external driving force. In section 6, we analyse period-averaged energy in dependence on the frequency of the external driving. In section 7, we study the Jarzynski-type relation for the work performed by the external time-periodic force. In section 8 , we discuss the regime of strong damping and discuss the problem of overdamped dynamics. Finally, in section 9 we conclude the paper with discussion and summary of the findings. Appendices contain supplementary materials on Green functions (Appendix A), the particle position and velocity mean square deviation (Appendix B) and the first two statistical moments of the work (Appendix C).
2 The forced Brownian harmonic oscillator
2.1 Langevin equation
We consider a classical Brownian particle of mass in contact with a thermal bath of temperature . For simplicity we consider its motion in one dimension in the harmonic potential . It is driven by an unbiased, symmetric time-periodic force of amplitude strength and angular frequency . The dynamics of the particle is described by a Langevin equation in the form :
Here, is the particle coordinate, the dot denotes differentiation with respect to the time and is the friction coefficient. Usually is the Stokes coefficient: , with the radius of the particle and the viscosity of the surrounding fluid. However, in this paper we do not refer to particle’s geometry and will simply take it as point-like one. The restoring potential force is . Thermal fluctuations due to coupling of the particle with the thermal bath are modeled by a -correlated Gaussian white noise of zero mean and unit intensity:
The noise intensity factor (where is the Boltzmann constant) comes from the fluctuation-dissipation theorem [36, 37, 38]. It ensures that the system stationary state is the canonical Gibbs state (thermodynamic equilibrium) under the dynamics governed by equation (1) for vanishing driving, i.e., when .
The following points are noteworthy:
The system can be seen as a Brownian particle moving in a harmonic potential or equivalently as a harmonic oscillator.
The stochastic process determined by equation (1) is non-Markovian. However, the pair constitutes a Markovian process.
Due to the presence of the external time-periodic force the system tends for long time to a unique asymptotic non-equilibrium state characterized by a temporally periodic probability density . In this long-time state, the average values of dynamical variables are periodic functions of time, that have the same period as the acting force .
2.2 Notations for averages
In this problem, we define three kinds of averages:
We denote by the ensemble averaging of a dynamical variable over all realizations of the Gaussian white noise and over all initial conditions .
We denote by the average value in the long-time limit state where all transient effects die out:
In this state, the dynamics is characterized by a temporally periodic probability density with the same period as the acting force: .
We denote by the time averaging of over one period in the long-time regime:
By doing this, we obtain the time-independent asymptotic characteristic average of the dynamical variable .
In the following sections, we will use the three notations: , and .
3 Solution of the Langevin equation
3.1 Position of the Brownian particle
The Langevin equation (1) is a stochastic non-homogeneous linear differential equation and its solution is in the form , where the complementary part satisfies the homogeneous differential equation (i.e., when and ), with the initial condition and . The particular solution satisfies the inhomogeneous equation with the initial conditions and . It is convenient to decompose the particular solution into two additive parts: , where is a particular solution corresponding to the deterministic force and is a solution related to the noise term . We therefore write
The term is a linear combination of two independent solutions of the homogeneous equation (in the case when and ):
where and are determined by initial conditions for the position and velocity of the particle which can be deterministic or random. The remaining parameters are defined as
The second part , related to the deterministic force , can be expressed by the Green function of equation (1) and takes the form 
and stands for the Heaviside step function: for and for . The method of obtaining the Green’s function is explained in details in A. The form of the Green function (9) is convenient when , i.e. when the system is weakly damped. In the opposite strongly damped regime, when , one can use the equivalent form with the hyperbolic sine function ,
Below, in all intermediate calculations, we will use only one form, namely equation (9).
The third term is non-deterministic and is related to thermal equilibrium noise . It can also be expressed by the Green function , namely,
Because is a Gaussian process, as a linear functional of is also Gaussian and in consequence is a Gaussian stochastic process. It means that two statistical moments and are sufficient to determine the probability density of the particle coordinate.
From equation (8) we get
Now, we can calculate the average value of the particle position at time . Because , from Eq.(11) it follows that . Therefore . For any initial conditions, when and in the stationary state, when all decaying terms vanish, we obtain
where is the long time limit of the deterministic solution (3.1), i.e. a purely periodic, non-decaying part of (3.1). Let us note that according to the recipe (4), the period-averaged position in the stationary state is zero, .
For long times, the second moment of the position takes the form
where we exploited equation (2) for statistical moments of thermal noise and neglected all exponentially decaying terms. The integral in equation (15) can be calculated and in long-time limit, we again can neglect the expenentially decaying terms to get the relation in the stationary state,
It can be verified that the mean squared deviation (MSD) of the particle position is . In Appendix B we present details of this calculation. The result turns out to be
Neglecting the exponentially decaying terms, we get the following relation in the stationary state,
which is notably mass independent! In fact the particle’s mass only determines the relaxation to the steady state value . In Fig. 1 we depict the temporal evolution of MSD for several values of the particle mass. As one can see, the larger the particle’s mass, the slower the relaxation towards steady state. The role of the mass will also be discussed later when we consider the overdamped regime.
3.2 Velocity of the Brownian particle
We now come to the particles’s velocity analysis. From equation (5) it follows that the particle velocity can be decomposed as
The first term can be calculated from equation (6) and for any initial conditions its average value if . The second term can be calculated form equation (3.1) and in the long time regime it takes the form
The last term is given by:
Its mean value is zero, , because both and . As a result, the mean velocity in the steady-state is
One can note that is simply the time derivative of given by equation (3.1).
The second moment of the velocity can be obtained by several methods. One of them is presented in Appendix B. In the steady-state it reads
Finally, the mean square deviation of the velocity in the steady-state is
which is remarkably independent of the spring stiffness ! More this value is the one that would be expected in the equilibrium situation (i.e. when ), due to the energy equipartition theorem .
the first part for the position (or for the velocity) is generated by the deterministic driving force ,
and the second part for the position (or for the velocity) is generated by the thermal noise (and is equal to the thermal equilibrium value).
This separation takes place because the system is linear and noise is additive.
4 Average energy in a stationary state
In this section, we compute the kinetic and potential energies of the Brownian harmonic oscillator. The kinetic energy of the oscillator:
is a stochastic quantity. Its potential energy:
is also a stochastic quantity. Naturally, the total energy of the oscillator is the sum . Its average value is determined by the second moments of the coordinate and velocity of the oscillator. Using the relations (24) and (20), the steady-state average kinetic energy can be expressed by the relation
where is defined in equation (13) and
The total averaged energy of the oscillator in the stationary state is
Both kinetic energy and potential energy consist of two terms: thermal equilibrium energy and the driving-generated part proportional to its amplitude squared . In turn, the driving-generated part consists of two terms: the part constant in time and the time-periodic part. For the latter, its mean value over the driving -period according to the rule (4) is zero. In contrast to the position and velocity , which are -periodic function of time, the energy is periodic with respect to time but now the period is .
Time-changes of energy in the stationary state are depicted in Fig. 2. In all panels of this figure, the rescaled kinetic energy is: ; the rescaled potential energy is: ; and the rescaled total energy is: , where the characteristic energy . One can note that the minimal values of the kinetic and potential energies are the equilibrium values . However, their maximal values depend on the ratio of two characteristic frequencies, namely, the frequency of the time-periodic force and the eigenfrequency of the oscillator. If (panel a in Fig. 2) then the amplitude of the potential energy is greater than the amplitude of the kinetic energy. In turn, if (panel b in Fig. 2) then the amplitude of the kinetic energy is greater than the amplitude of the potential energy. In both cases, the total energy oscillates in time. Finally, if (panel c in Fig. 2), the total energy does not depend on time, it is constant! Remember that the system is open, not isolated. There is dissipation of energy due to contact with thermostat of temperature and there is a pumping of energy by the time-periodic force . Under some specific conditions (here ), these two processes can lead to the conservation of total energy in the noisy system.
5 Energy balance
We now derive equations for the energy balance. To do this, we first rewrite equation (1) in the Ito-type form of stochastic differential equations:
where formally and is the Wiener process of Gaussian statistics with zero mean and the second moment . Next, we employ the Ito-differential calculus to the functions and , namely,
where we have only kept the terms that at most are first order in . We then perform the ensemble averaging over all realizations of the Wiener process and over all initial conditions . As a result we obtain the rate of change of the potential and kinetic energies in the form:
where for the part containing the Wiener process in equation (LABEL:dEk) we used the Ito-martingale property, i.e. . By adding these two equations together, we get an equation for the energy balance. It reads:
We conclude that there are three processes responsible for the energy change: (i) the first term in the r.h.s. of equation (38) is always negative and represents the rate of energy loss due to dissipation (friction of the fluid); (ii) the second term describes pumping of energy to the system by the external force and (iii) the third term is always positive and characterizes the energy supplied to the system by thermal equilibrium fluctuations. We recall that in thermal equilibrium (i.e. when ), the theorem on the equipartition of energy states that , from which it follows that . Therefore equation (38) can be rewritten in the form:
Let us note that when , in the stationary state (which is thermal equilibrium), the l.h.s. tends to zero and we retrive for any initial conditions.
Now, we perform the time averaging of equation (39) according to the prescription (4). Let us remember that in the stationary state the average values of any functions of the particle coordinate and its velocity are periodic functions of time and therefore the left hand side equation (39) vanishes after time averaging. Indeed, by definition:
As a consequence, in the stationary regime the mean power delivered to the system by the external force over the period is expressed by the relation
where the characteristic time is the relaxation time of the particle velocity for the system when . This relation has an appealing physical interpretation: the input power from the external force to the system determines the difference between the averaged non-equilibrium and equilibrium values of the kinetic energy of the Brownian particle. Moreover, the input power depends not only on the force itself (i.e. on the amplitude and frequency ) but also on properties and parameters of the system: the friction coefficient and the mass . In contrast, the energy supplied by thermal fluctuations does not depend on the external force but only on and .
6 Period-averaged energy
In the non-equilibrium stationary state all averaged forms of energy are periodic functions of time. Time independent characteristics can be obtained from the formulas (4)-(31) by additional averaging over the period of the time-periodic force . The averaging scheme is according to the rule (4). In this way we obtain the period-averaged kinetic energy in the form
The period-averaged potential energy reads
and the total period-averaged energy is given by the expression
The period-averaged input power can be obtained form equation (41) and the result is
We recall that the time averaged input power is independent of temperature of the system. It is not the case for non-linear systems  for which the input power is temperature dependent.
Now, we analyse properties of these time-independent mean energies. The dependence of energy on the rescaled frequency of the external force is shown in Fig. 3. We observe that:
The kinetic energy is zero at zero driving frequency and approaches zero when . It means that there must be an optimal value of the driving force frequency for which the kinetic energy is maximal. Calculation of derivative shows that the maximal value is for the resonance frequency .
The shape of the potential energy depends on the rescaled friction coefficient : If , the potential energy takes a bell-shape form of the maximal value at . However, if , there is one minimum at and the maximum at .
The shape of the total energy also depends on the rescaled friction coefficient: If , the bell-shaped total energy is maximal at and tends to zero when the frequency increases. In turn, if , it exhibits one minimum at and the maximal value at .
The dependence of the input power on the rescaled friction coefficient and the rescaled mass is shown in Fig. 4. In panel a, the scaled frequency is and the rescaled input power is given by the formula:
where the characteristic power is . The minimum input power is always at , the maximum is always at the resonance and does not depend on . In panel b, the scaled frequency is and the rescaled input power now reads:
where the characteristic power is and the rescaled mass is . The minimum input power is always at . The maximal input power is which is located at .
7 The Jarzynski-type work relation
Many physical observables undergo random fluctuations. The coordinate, velocity, kinetic and potential energies of a Brownian particle are generic examples of such observables. Work and fluctuation relations constitute an additional basic characterization of statistical properties of the system. In this section, we present the counterpart of the Jarzynski work relation for the considered system. We have to stress that what we present is not an exact Jarzynski relation because conditions for the Jarzynski theorem to be hold are not fulfilled and the interpretation is not in terms of a free energy, the notion which is defined for systems in thermal equilibrium but not for non-equilibrium ones.
Although the conception of work seems to be well known, there is no consensus on its precise and general definition, see discussion in Refs. [9, 44, 45]. Here we consider the inclusive work defined as
where is given by equations (5)-(13). It is the work performed by the external force on the system in the time interval . We are interested in the limit of long times when the asymptotic state is reached. Therefore the regime should be considered.
The coordinate , via the relation (11), is a linear functional of the Gaussian process . Therefore for each fixed time the work (48) as a function of time
is a stochastic Gaussian process. However, because of the time-periodic driving , it is a non-stationary random process and therefore the statistical properties of the work can depend on the reference time ! It is a radically different case than for stationary processes.
The probability distribution of the work is a Gaussian function in the form 
are the mean value and variance of the work, respectively.
The characteristic function of the Gaussian stochastic process is well known, namely ,
Now, let us assume that , where . Then it takes the form
This characteristic function of the work for the imaginary argument is just a counterpart of the Jarzynski relation:
where is the free energy difference between the final and initial equilibrium states. The argument of the exponent in r.h.s. of (52) formally corresponds to in (53) but cannot be interpreted as a free energy of the system.
In Appendix C, we calculate both and which allow to calculate . As a result we get
where is defined by equation (45) and the remaining notation is explained in Appendix C.
In the literature [47, 48, 49], the relation (52) is called a stationary state fluctuation relation for the total work performed during the time . However, in Ref. [47, 49], a special case of the stationary state is considered, namely, when both statistical moments and do not depend on the reference time . In our case, both and do depend on , cf. equation (C.1) and (C.7) in Appendix C. Moreover, the well-known equality does not hold, even in the limiting case . It is only one term in and one term in which connotes this equality.
8 Strong damping regime
The regime of strong damping has played an important role in study of deterministic and noise-assisted dynamical systems and is usually relevant for experimental micro-systems. The common way to treat the overdamped dynamics is to put mass zero, , in Eq. 1. However, this approach is not fully rigorous because of several reasons. First, the mass of the particle is never exactly zero: it can be ’small’, but then remains the question ’small’ in comparison to what? Second, the energy balance cannot be fully evaluated since some terms do not exist separately but only in specific combinations. It can give rise to squares of the Dirac -functions which do not exist . In other method the order of magnitude of the dimensionless quantity is considered. The problem with this criterion is that the fraction is not constant in time and it is even worse because it can diverge to infinity. A more correct procedure is to identify regimes where ’something’ is small or ’something’ is large. Moreover, one can use a technique of mathematical sequences and try to obtain an exact limit of a considered sequence. We present the following scheme: We transform equation (1) to its dimensionless form and work in dimensionless variables, because from the physical point of view only relations between scales of time, length and energy play a role and not their absolute values. For the system defined by Eq. 1, one can identify four characteristic times:
The relaxation time of the velocity , which can be extracted from the special case of equation (1), namely, .
The relaxation time of the position deduced from the equation .
The reciprocal of the frequency of small oscillations obtained from the equation .
The periodicity time related to the frequency of the external driving .
To analyze strong damping regime, one has to choose the characteristic time as a scaling time since it does not contain the particle mass and therefore the dimensionless time is defined as . We also introduce the dimensionless position , where is some characteristic distance. In the case of the harmonic potential, the only reasonable length scale is defined by the equipartition theorem . In consequence a characteristic distance is . Then the dimensionless Langevin equation corresponding to Eq. (1) takes the form:
where the rescaled quantities read
The rescaled Gaussian noise has the same statistics as in equation (2). Note that (i): in the left-hand-side of equation (55) there is only one parameter which is the same as in equation (47) and (ii): the intensity of the rescaled thermal noise is 1. The parameter is a ratio of two characteristic times of the system and if then , and we can neglect the inertial term. Observe that the inequality means that and it defines a real physical regime of a strong damping regime. It depends not only on mass and friction coefficient but also on the stiffnes of the potential trap! In literature, this fact is often ignored.
Now, we can rewrite Eqs. (42)-(43) in the dimensionless form . The rescaled kinetic energy is (here we consider only the contribution of the forcing on the energies, so we subtract the contribution due to the thermal agitation):
where the rescaled frequency and . For the potential energy we get
The rescaled input power is given by equation (47). We find that if decreases to zero, the rescaled kinetic energy decreases to zero. In this limit, the potential energy takes the similar form like the Cauchy probability distribution, namely,
Similar behavior exhibits the total energy because kinetic energy is zero. As follows from equation (47), if the input power assumes the form
It is zero for zero frequency (i.e. when the force is constant) and saturates to the value when the frequency . However, it is not a correct physical result. We should stress that from the exact formulas (45)-(47) it follows that for any non-zero but arbitrary small values of parameters, when . In particular, in the strong damping regime, the mass parameter is small but non-zero. So, the reader should take care with the order of limits.
The Brownian harmonic oscillator is an extremely simple and at the same time paradigmatic model with many applications in both classical and quantum physics. In the paper, we considered a less popular problem of energetics of the system, which is permanently moved out of thermal equilibrium, but reaches a stationary state in the long-time regime. Because the model is exactly solved, all interesting characteristics can be analytically evaluated. The first conclusion is that in the stationary state, various forms of averaged energy are additive: one part is related to thermal-equilibrium energy according to the equipartition theorem and the second part is generated by the external unbiased time-periodic force of frequency . Additivity is due to linearity of the system (i.e., the potential force is a linear function of the particle coordinate) and for non-linear systems it does not hold true. In the long time limit, when a non-equilibrium stationary state is reached, energy is a periodic function of time of the same period as the external driving. If (where is the eigenfrequency of the undamped oscillator), the maximal value of the kinetic energy is greater than the maximal value of the potential energy. In turn, if , the maximal value of the kinetic energy is smaller than the maximal value of the potential energy. For the resonance frequency , the maximal value of the potential energy is equal to the maximal value of the kinetic energy and the total average energy is conserved, i.e., it is a constant function of time. It is the second important result.
We analysed time-independent characteristics of energy which are obtained by the method of averaging over the period of the time-periodic force. We observe that in some regimes, determined by the dissipation strength (quantified by the friction coefficient ), the frequency dependence of energy exhibits non-monotonic behaviour. New results are obtained for the work-fluctuation relation of the Jarzynski type. We analysed this problem in the regime of long times, when the non-equlibrium stationary but time-dependent state is reached. Both the mean value and variance of the work performed on the system in the interval grows linearly with time (as it should be) but also depend periodically on the reference time . Finally, we critically discussed the so-called overdamped dynamics. We are of opinion that the correct framework for studying a regime of strong damping has to be based on analysis of time scales of the system. The overdamped dynamics cannot simply be obtained by putting the particle mass (which can bring physical nonsense), but rather requires to identify the correct parameters to compute a dimensionless equation, and to calculate the limits for each quantity of interest. It is what we presented in Sec. VII. This procedure is always correct and we hope that we were able to clarify a physically acceptable method of an approximation of some limiting dynamics on a simple example of the driven Brownian harmonic oscillator.
Appendix A Obtaining the Green’s function
In this appendix we obtain the Green’s function of the full Langevin equation. According to the theory of Green’s function [mintU], we first should find the solutions of the corresponding homogeneous differential equation:
The solution takes the exponential form , where verifies:
Assuming for simplicity that , we have two real solutions and :
which gives us two independent functions and , respectively.
Now, we construct the Green’s function that verifies: