Stationary energy probability density of oscillators driven by a random external force

We derive rigorous analytical results for the stationary energy probability density function of linear and nonlinear oscillators driven by additive Gaussian noise. Our study focuses on two cases: (i) a harmonic oscillator subjected to Gaussian colored noise with an arbitrary correlation function and (ii) nonlinear oscillators with a general potential driven by Gaussian white noise. We also derive analytical expressions for the stationary moments of the energy and investigate the partition of the mean energy between kinetic and potential energy. To illustrate our general results, we consider specifically the case of exponentially correlated noise for (i) and power-law and bistable potentials for (ii). Our theoretical results are substantiated by Langevin simulations.


I. INTRODUCTION
The study of the effects of noise on oscillators has a long and distinguished history [1].Early studies focused on the motion of a harmonically bound Brownian particle in a viscous fluid.The Brownian particle experiences thermal noise, modeled as a Gaussian white noise process.The evolution equation for the position is then a Langevin equation, a stochastic differential equation.The corresponding Fokker-Planck equation for the probability density of the position and velocity can be readily derived and solved [2][3][4][5].Kramers considered a nonlinear oscillator, a system with a double-well potential, to study the noise-activated escape from a potential well or the crossing of an energy barrier [6,7].The Kramers oscillator is essentially the Duffing oscillator [8,9], which has played the role of a paradigm in nonlinear dynamics [10].The Duffing oscillator has found many applications to mechanical systems (see, e.g., [9]), but it also has been employed to describe nonlinear electrical circuits [11,12] and to model certain aspects of Rayleigh-Bénard convection [13,14].More recently, the Duffing oscillator has served as a model for nanoand micro-mechanical oscillators [15][16][17][18].
In the 1980s, studies of linear and nonlinear oscillators were extended to the case of colored noise driving forces.The effect of correlations in the random driving force on the stationary probability density and its moments were investigated [19][20][21][22], as were related problems, such as the escape from a potential well and the mean first passage time [7], stochastic resonance [23][24][25], the stability of metastable states [26], and the dynamics of excitable systems [27].Much less attention has been paid to how the color of the noise and the nonlinearities of the potential affect the statistical properties of the particle's energy.This is a important topic, since oscillators have recently been considered for harvesting energy from ambient fluctuations to power small self-contained sensors and actuators [28][29][30][31][32][33][34][35].The aim of this paper is to begin addressing this gap.We derive analytical expressions for the stationary probability density function (PDF) of the particle's energy.We determine the mean energy and its variance and compare our theoretical predictions with Monte Carlo simulations.We further investigate the partition of the energy between the kinetic energy and the potential energy of the oscillator.Our studies focus on two cases.First, we consider a harmonically bound particle driven by general additive Gaussian colored noise.Second, we consider a particle in a general potential driven by Gaussian white noise.In both cases we develop a general theory and study specific illustrative examples.
The paper is organized as follows.In Sec.II we formulate the Langevin equations for a particle with inertia in a general potential driven by an additive random force.We also derive the general relation between the stationary energy PDF and the stationary PDF for the particle's position and velocity.Section III is devoted to the motion of a particle in a harmonic potential driven by colored Gaussian noise with an arbitrary correlation function.The stationary energy PDF, as well as the first and second moments of the energy, is derived analytically.We apply these results to the case of exponentially correlated Gaussian noise and to Gaussian white noise.For these specific cases, we also determine the partition of the energy between kinetic and potential energy.We analyze the effect of the potential shape on the stationary energy PDF for oscillators driven by an additive Gaussian white noise in Sec.IV.We apply our analytical results to two specific cases: a general power-law potential and a bistable potential.The first type of oscillator, without friction and with additive noise [36] and multiplicative noise [37], has recently been studied in the context of anomalous diffusion.The second type, the Duffing oscillator, has been investigated as an alternative to the linear oscillator for harvesting energy from external random sources [28,29,35,38,39].

II. STATIONARY ENERGY PROBABILITY DENSITY
While linear and nonlinear oscillators are encountered in a wide variety of applications as discussed in Sec.I, we adopt a mechanical picture for the sake of concreteness in the following.To be specific, we consider the motion of a particle in a one-dimensional fluctuating environment.The state of the particle is fully described by its position x and its velocity v.Both variables are random and evolve according to the system of stochastic differential equations (Langevin equations), Without loss of generality, we assume that the particle has unit mass.It moves in the potential V (x) whose gradient gives rise to a force acting on the particle.The potential is determined only up to an additive constant.We choose that constant such that the deepest potential well corresponds to an energy level of zero.The particle also experiences friction.
We assume the friction force to be a linear function of the velocity, and b is the friction coefficient.The fluctuations of the environment are represented by the external random force ξ (t) with mean zero, ξ (t) = 0. We consider only nonequilibrium situations.For such systems, the fluctuating force ξ (t) and the friction coefficient b are not related by a fluctuation-dissipation theorem.Our aim is to describe the evolution of the particle's energy.Since position and velocity are random quantities, the energy of the particle, is also a random quantity.(We use T to denote the kinetic energy of the particle.)The temporal evolution of its mean value is given by 3) can be interpreted as the balance equation for the power of the particle.Two processes contribute to the rate of change of the particle's energy.The first term is negative and represents the rate of energy dissipation due to friction.The second term corresponds to the rate of exchange of energy with the external noise.We expect that the particle eventually reaches a steady state where the dissipation of the energy is balanced by the injected energy; i.e., E becomes constant.In the steady state where the subscript s denotes "steady state."Our goal is to find the PDF of the particle's energy in the steady state.This PDF can be obtained from the stationary density for the position and velocity of the particle, P s (x,v), by applying the following relation between the PDFs of related random variables [40]: where The composition δ[g(x)] for continuously differentiable functions g is defined by where x i are the roots of g(x), which are assumed to be simple, and the sum runs over all roots [41].We apply this property to where v 0 (x) is the solution of the equation and taking into account (2.9), we obtain from (2.5) . (2.10)

III. LINEAR OSCILLATOR
In this section we analyze the effect of Gaussian noise on the stationary energy PDF of a particle in a harmonic potential, V (x) = ω 2 x 2 /2.We assume that the correlation function of the Gaussian noise ξ (t) is arbitrary: (3.1) Since the system (2.1) is linear in this case and since the external driving force is Gaussian, the stationary PDF for the position and velocity of the particle will also be Gaussian and reads, with u = ( x v ), where m = x s v s (3.3) denotes the mean, and A represents the covariance matrix, Taking the average of (2.1) in the steady state, we find that The elements of the covariance matrix A are given by [42][43][44]] where Here 2 = ω 2 − b 2 /4, and we assume that the oscillator operates in the underdamped regime, i.e., 4ω 2 > b 2 .Since we are interested in the stationary state, this assumption does not imply any lack of generality.The final steady state is the same, regardless of whether the transient behavior is overdamped or underdamped.Given the form of G(t) and H (t) and using Euler's formula to represent the sine and cosine function, we find that the elements of the covariance matrix (3.6) correspond to Laplace transforms with a complex argument, s = b/2 ± i , of the correlation function C(t).Defining Ĉ(s) to be the Laplace transform of C(t) with argument s, after some calculations we obtain the following expressions: (3.8a) where Re[•] and Im[•] represent the real and imaginary parts.
Using the elements of the covariance matrix (3.8) in the stationary PDF (3.2) for the position and velocity of the particle, we find Substituting (3.9) into (2.10),we obtain, after some calculations, where and I 0 (•) is the modified Bessel function.In (3.10) we display explicitly the dependence on E, i.e., the parameters λ ± depend on b, ω, and the parameters involved in the correlation function C(t), but they do not depend on the energy.For large z the function I 0 (z) behaves asymptotically as I 0 (z) ∼ e z/2 / √ z, and the tail of the stationary energy PDF decays as The average of the energy in the stationary state is given by where we have used (3.8).Similarly, (2.4) implies that the power in the stationary state is given by The mean squared energy can be calculated with the help of (3.10) and is (3.15) We apply our results for general external Gaussian forces to two specific cases.

A. Exponentially correlated Gaussian noise
As a first example, we consider the case of Gaussian colored noise with an exponential correlation function: The Laplace transform of the correlation function, (3.17b) (3.18b) From (3.10) we obtain the stationary energy PDF, for E ∈ [0,∞).In Fig. 1 we plot the stationary PDF (3.19), using solid curves for different values of τ , and compare the results with the corresponding simulations.As can be appreciated the stationary energy PDF shows a monotonic behavior with respect to the correlation time τ .It decays faster as τ increases.The numerical simulations are performed by solving numerically the Langevin equations (2.1).We have discretized the equations by using the so-called Heun algorithm [45].At each time step we generate Gaussian white random numbers for ξ (t) by means of the Box-Mueller-Wiener algorithm [45].For the case of Fig. 1 we have added to the system (2.1) a third equation for the variable ξ (t), namely, dξ/dt = −ξ/τ + aη(t)/ √ τ , where η(t) is Gaussian white noise with correlation η(t)η(t ) = 2a 2 τ δ(t − t ).In this case the system of stochastic equations have been discretized again by using the Heun method and the random numbers for η(t) have been generated by the Box-Mueller-Wiener algorithm.
To obtain mean values we have averaged over 10 5 realizations of the noise.
The stationary mean energy can be calculated directly from (3.19), from (2.2) using the second moments (3.18), or from (3.13): (3.20) Equation (3.14) implies that the power dissipated by the particle and the power supplied by the Gaussian colored noise in the steady state are given by It is proportional to the noise intensity D = a 2 τ and inversely proportional to the damping coefficient b and the square of the oscillation frequency.We obtain the expression for the mean squared energy from (3.15), and the variance of the energy fluctuations is given by The coefficient of variation or relative standard deviation, a measure for the dispersion of a PDF, is defined as the ratio of the standard deviation to the mean.Keeping the amplitude of the noise a fixed and varying the correlation time τ from 0 to ∞ we find The coefficient of variation increases as the correlation time τ increases from 1 at τ = 0 to √ 2 as τ → ∞.For any τ , the fluctuations of the energy are of the same order as the mean energy due to the Gaussian nature of the noise.
Note that (3.18) implies that the stationary mean potential energy and the stationary mean kinetic energy of the oscillator are given by Consequently, we obtain for the ratio of the mean potential energy to the mean kinetic energy This ratio does not depend on ω −1 , the characteristic time of the undamped oscillator, but it depends on τ , the characteristic time of the external random driving force.More precisely, it depends on the ratio of the correlation time to the characteristic time of the dissipative force, b −1 .Note that correlations in the external Gaussian noise prevent equipartition of the energy.
For nonvanishing τ , the mean potential energy always exceeds the mean kinetic energy.

B. Gaussian white noise
As a second example, we consider the case of Gaussian white noise.The white noise limit corresponds to a → ∞ and τ → 0 keeping D = a 2 τ < ∞ fixed.In this limit, (3.19) yields In other words, the energy is exponentially distributed in the steady state.In Fig. 2 we compare the theoretical result of (3.27) with numerical Monte Carlo simulations.
Taking the Gaussian-white-noise limit of (3.20)-(3.23),we obtain ) It is easily verified that the mean energy E s and the power vξ s of the particle are larger for Gaussian white noise than for exponentially correlated Gaussian colored noise.Note that the stationary energy PDF does not depend on the frequency ω for Gaussian white noise, in contrast to the case of colored noise.As noted above, the ratio of mean potential energy to the mean kinetic energy does not depend on ω, but solely on τ .For Gaussian white noise, the system achieves equipartition of the energy, (3.33)

IV. NONLINEAR OSCILLATOR
In this section we analyze the effect of the shape of the potential on the stationary energy PDF for systems driven by Gaussian white noise, The stationary PDF P s (x,v) for the system (2.1) with Gaussian white noise is given by (see, for example, [46][47][48]) where Z is the normalization constant, Note that the stationary PDF (4.2) factorizes as In other words, the position and the velocity of a particle driven by an additive external Gaussian white noise force are independent random variables in the steady state for the motion in an arbitrary potential.Equations (2.10) and (4.2) imply that If we define μ = b/D, then the stationary PDF (4.2) formally looks like a Boltzmann equilibrium distribution, We proceed in a similar way to calculate For the mixed expectations values, we exploit the stochastic independence of x and v and find Equation (4.9) implies that the stationary mean kinetic energy is given by for any potential.

A. Power-law potential
In this section we consider the case of a general powerlaw potential, V (x) = kx 2n , n 1 [37].The normalization constant (4.3) is given by and (4.5) provides the stationary energy PDF,

.14)
Note that for n = 1 we recover the case of the linear oscillator (3.27).For E → 0 + , the PDF (4.14) diverges if n > 1, while for the case n = 1 (harmonic oscillator) it approaches b/D.In Fig. 3 we compare the analytical result (4.This result allows us to compare the linear and nonlinear oscillators in terms of their mean energy and their energy fluctuations, i.e., their variance.For the mean energies we find In other words, the mean energy of the linear oscillator exceeds the mean energy of nonlinear oscillators.The variance of the energy in the steady state is given by Note that the variance of the energy for a nonlinear oscillator is smaller than that for the linear harmonic oscillator.This suggests that a nonlinear potential "absorbs" the energy fluctuations by reducing their amplitudes.Combining (4.12), i.e., the stationary mean kinetic energy is independent of the potential, T s = D/(2b), and (4.16), we determine the mean potential energy as The mean kinetic energy of a particle driven by Gaussian white noise in a power-law potential V (x) = kx 2n that is softer than the harmonic potential, i.e., n > 1, always exceeds the potential energy for such a system, n V (x) s = T s (4.19) and Q = 1/n.Equation (4.19) shows that the power-law oscillator with additive Gaussian white noise obeys the virial theorem [49].As the potential becomes softer, i.e., as n becomes larger, more of the energy of the particle will on average be kinetic energy.

B. Bistable potential
In this section we consider case of the bistable potential System (2.1) with this potential is the well-known Duffing oscillator.We are interested in analyzing the effect of the barrier height on the stationary energy PDF, in particular on how the mean energy of the oscillator depends on .Simple analysis shows that the potential (4.20) has two minima at x = ± √ a/δ and a maximum at x = 0.The height of the potential barrier is V (0) − V (± √ a/δ) = a 2 /4δ.As stated in Sec.II, we choose the additive constant of the potential such that the deepest well corresponds to an energy level of zero.Therefore we set = a 2 /4δ, which results in V (± √ a/δ) = 0 and V (0) = .
To evaluate the integral in (4.5), we have to determine the domain of integration V (x) E, which is defined as the interior of the boundary V (x) = E.This boundary will be different for E > and E < .The integral needs to be split into two parts, one for 0 < E < and one for E > .Since E = corresponds to a separatrix of the deterministic problem, we expect the stationary energy PDF to diverge at E = .
Within the region 0 < E < , the boundary V (x) = E has four real roots for x.Since the potential is an even function, we integrate over x ∈ [x ,x + ], where and multiply by 2. Equation (4.5) yields where K(•) is the complete elliptic integral of the first kind.Within the region E > , the boundary V (x) = E has only two real roots for x.Due to the symmetry of the potential we integrate over x from 0 to x + and multiply by 2. Equation (4.5) yields P (2)  s We need to determine the normalization constant by integrating P s (E) over E for E ∈ [0,∞), 0 P (1)  s (E)dE The final result reads Here In Fig. 4 we compare the analytical result (4.25) with numerical simulations.As we expected, the stationary energy PDF diverges at E = .
By virtue of (4.12), the stationary mean kinetic energy is independent of the potential, that is, T s = D/2b for the Duffing oscillator.To evaluate the stationary mean potential energy, we use (4.8) and find This implies that and the stationary mean total energy is given by We evaluate the second moment x 2 s directly from (4.2) and (4.3), obtaining Introducing the new variable u = δbx 4 /4D, we obtain To evaluate the integrals, we employ the following series expansion: and obtain, using (4.27), where we have defined Here I ν (•) is the modified Bessel function of order ν.We substitute this result into (4.30) to find As in the previous sections, we investigate the partition of the energy in the stationary state between the potential energy and the kinetic energy.For the Duffing oscillator, the ratio of the mean potential energy to the mean kinetic energy in the stationary state is given by Substituting the expression (4.34) into (4.37) and using (4.27), we obtain Note that Q, in contrast to E s , depends only on the dimensionless potential barrier β and not on the individual parameters of the random Duffing oscillator.On the other hand, (4.36) and (4.38) show that the behavior of the stationary mean energy and of the ratio of the potential energy to the kinetic energy is qualitatively the same as the dimensionless barrier β is varied.In Fig. 5 we compare the stationary mean energy (4.36) with numerical simulations.In this case the fluctuations are even greater than for the previous oscillators studied here, and it is difficult to obtain a mean value with sufficient accuracy.The behavior of Q as a function of the dimensionless potential barrier β is shown in the inset of Fig. 5.
Using the limiting form of the modified Bessel function I ν for small arguments [50], On the other hand, using the asymptotic expansion of the modified Bessel function I ν for large arguments [50], This maximum value has been predicted numerically [28,29,32] but has never been established analytically.At this point, the mean potential energy exceeds the mean kinetic energy by only about 16%, and the mean total energy is only about 8% larger than the mean energy for the case of a harmonic potential.A further increase of the potential barrier leads to a monotonic decrease of Q and E s toward the asymptotic values of 1 and D/b, respectively.At β = 20.0,V (x) s is already only about 2% larger than T s .
As in previous sections we study now the energy dispersion.To compute the energy variance we have to compute first E 2 s .Taking into account (2.2) and (4.12) we find On the other hand, we can make use of (4.20) and (4.28) to get The moments x 6 s and x 8 s can be computed from (4.2) straightforwardly as we did to find (4.32).The expression for the relative energy dispersion can be finally obtained as a function of the dimensionless potential barrier β only:

.50)
Since β ranges from 0 (where the potential resembles the quartic potential) to ∞ (where the potential is harmonic) it can be easily shown from (4.50) that 0.9

.51)
The maximum value for the relative energy dispersion, 2/ √ 3, is reached when β → 0 + , i.e, when the potential is quartic.Note that this result equals that obtained from (4.16) and (4.17) with n 2. On the other hand, the minimum is not reached for β → ∞ where σ EE / E s = 1 but at β = 1.587where σ EE / E s reaches a local minimum.

V. CONCLUSIONS
The problem of harvesting energy from ambient fluctuations to power small devices in a self-contained manner has attracted much attention.Linear and nonlinear oscillators have been considered for this purpose [29,30,32,35,39].It is therefore important to understand the interplay of the potential shape of an oscillator with the nature of the external random driving force on the energy of oscillators in the steady state.
We have addressed this problem by analyzing the stationary energy PDF of oscillators driven by additive Gaussian noise and focused on two cases, namely, linear oscillators driven by Gaussian noise with an arbitrary correlation function and oscillators with a general nonlinear potential driven by Gaussian white noise.We have applied our analytical results to specific examples, such as exponentially correlated external noise and the Duffing oscillator, and substantiated our theoretical results by Langevin simulations.
Besides obtaining rigorous expressions for the stationary energy PDF and its moments, we have also derived analytical results for the mean kinetic and potential energy.We find that exponential correlations in the Gaussian driving force destroy the equipartition of kinetic and potential energy in the harmonic oscillator.The larger the correlation time compared to the dissipative time scale, the larger the mean potential energy compared to the mean kinetic energy.We find that for the power-law oscillator, V (x) = kx 2n , the virial theorem holds for additive Gaussian white noise.
The mean total energy of a Gaussian-white-noise-driven Duffing oscillator, as well as the ratio of mean potential energy to mean kinetic energy, depends in a nonmonotonic way on the barrier height.We have provided an analytical derivation that these quantities pass through a minimum and then a maximum as the barrier height between the two potential wells is increased.The oscillator reaches a state of equipartition between kinetic and potential energy as the dimensionless potential barrier β goes to infinity.For β < 1.524 664, the mean kinetic energy exceeds the mean potential energy and achieves its largest value relative to the latter at β min = 0.0654 649 6.This is also the point of minimum mean total energy.For β > 1.524 664, the mean potential energy exceeds the mean kinetic energy and achieves its largest value relative to the latter at β max = 3.325 336, which is also the point of maximum total energy.Our results demonstrate that the mean total energy of the Gaussianwhite-noise-driven Duffing oscillator is minimal when the mean kinetic energy is maximal relative to the mean potential energy.Conversely, the mean total energy is maximal when the mean kinetic energy is minimal relative to the mean potential energy.
There are three main directions that we will explore in future work.The first is the effect of correlations in the random driving force on the stationary energy PDF and the mean kinetic and potential energy of the power-law and Duffing oscillators.The second is the effect of non-Gaussian noise.While the assumption of a normal distribution for the ambient fluctuations is often justified due to the central limit theorem, the external force on the oscillator may be of an impulsive type, which is better modeled by Poisson noise [32].Finally, we will explore the effect of the asymmetry of the potential shape and the presence of periodic forcing on the stationary energy PDF.

FIG. 1 .
FIG. 1. Stationary energy PDF of a particle in a harmonic potential and driven by external Gaussian noise with an exponential correlation function (3.16).The values of the amplitude a and the correlation time τ are listed on the figure.The other parameters are ω = b = 1.The symbols correspond to numerical simulations and the curves represent results obtained from (3.19).

. 30 )FIG. 2 .
FIG.2.Stationary energy PDF of a particle in a harmonic potential and driven by external Gaussian white noise.The parameter values are ω = 0.63, b = 0.2, and D = 0.5.The jagged curve corresponds to the results of numerical simulations and the solid curve corresponds to(3.27).Inset: The same results plotted on a semilog scale.
14) with numerical simulations for different values of n.The rth-order moment of the stationary energy PDF (4

FIG. 3 .
FIG. 3. Stationary energy PDF of a particle in a potential V (x) = kx 2n and driven by external Gaussian white noise.The parameter values are k = 0.2, b = 0.2, and D = 0.5.Each panel corresponds to a different value of the exponent n.The symbols correspond to numerical simulations and the curves are obtained from (4.14).

. 26 )FIG. 4 .
FIG. 4. Stationary energy PDF of a particle in the bistable potential (4.20) and driven by external Gaussian white noise.The parameter values are a = b = D = 1.Each plot corresponds to a different value of the potential barrier .The symbols correspond to numerical simulations, while the solid curves are obtained from (4.25).

41 )FIG. 5 .
FIG. 5. Mean energy of a particle in the potential (4.20) and driven by external Gaussian white noise.The parameter values are a = b = D = 1.The symbols correspond to numerical simulations, while the solid curve represents (4.36).Inset: Ratio of the mean potential energy to the mean kinetic energy Q of a particle in the bistable potential (4.20) and driven by external Gaussian white noise vs the dimensionless potential barrier β.The solid curve represents Eq. (4.38) and the symbols correspond to numerical results.andaccording to(4.19)