text:
# Anomalous diffusion in nonlinear oscillators with multiplicative noise

## I Introduction

Randomness in the external conditions entails the parameters of a dynamical system to fluctuate. The extent of these fluctuations is independent of any thermodynamic characteristic of the system in contrast to intrinsic fluctuations the amplitude of which is proportional to the equilibrium temperature, in accordance with the fluctuation-dissipation theorem vankampen ; risken . Usually, external randomness appears as a multiplicative noise in the dynamical equations. The interplay of noise and nonlinearity in a system far from equilibrium results in some unusual phenomena lefever . In fact, the presence of noise dramatically alters the properties of a nonlinear dynamical system both qualitatively and quantitatively (for a recent review see landaMc ). For example, it was shown recently that in a spatially-extended system, a multiplicative noise, white in space and time, generates an ordered symmetry-breaking state through a non-equilibrium phase transition whereas no such transition exists in the absence of noise vandenb1 ; vandenb2 . Noise can also induce spatial patterns vandenb3 ; toral  or improve the performance of a nonlinear device through stochastic resonance barbay . Furthermore, even if some important qualitative features of a deterministic system survive to external noise, their quantitative characteristics may change: a stable fixed point may become unstable bourret , a bifurcation may be delayed (noise-induced stabilization) lucke ; roder  and scale-invariant properties which manifest themselves as power-laws may be altered with the appearance of non-classical scaling exponents genovese .

The discovery of Brownian motors which are able to rectify random fluctuations into a directed motion (noise-induced transport) has triggered a renewed interest in the study of simple one-dimensional mechanical models of particles in a potential with random parameters reimann . It is well known that a linear oscillator submitted to parametric noise can be unstable even if damping is taken into account bourret ; linden . This noise-induced energetic instability has been observed in diverse experimental contexts such as electronic oscillators strato ; kabashima , nematic liquid crystals kawakubo  and surface waves (Faraday instability) fauve . In engineering fields this instability plays a crucial role in the study of the dynamic response of flexible structures to random environmental loading such as the wave-induced motion of off-shore structures or the vibration of tall buildings in a turbulent wind roberts . The presence of nonlinear friction tends to limit the oscillation amplitude: the pendulum with randomly vibrating suspension axis and undergoing nonlinear friction, known as the van der Pol oscillator, has been studied in the small noise limit using perturbative expansions strato ; landa .

In the present work, we consider the motion of an undamped nonlinear oscillator trapped in a general confining potential and submitted to parametric random fluctuations. Because there is no dissipation, the energy of the system increases with time and we shall show that the position, the momentum and the energy grow as power-laws of time with nontrivial exponents that depend on the behavior of the confining potential at infinity bouchaud . A key feature of our method is to use the integrability properties of the associated deterministic nonlinear oscillator in order to derive exact stochastic equations in action-angle variables. We then use the averaging technique of classical mechanics landau , together with a reduction procedure vandenbroeck1 ; drolet , to calculate exactly the anomalous scaling exponents, irrespective of the amplitude of the noise. Some of our results were derived before in the particular case of a cubic nonlinearity using an energy-enveloppe equation lind2 . Our method enables us to derive the numerical prefactors appearing in the scaling laws (generalized diffusion constants), and our analytical predictions compare very satisfactorily with the numerical results. In the case of noise correlated in time, the anomalous diffusion exponents are modified: they can be obtained by dimensional analysis arguments and the values thus found also agree with numerical results. Throughout this work, crossover phenomena between different scaling regimes are emphasized.

This article is organized as follows. In section II, we recall that the energy of a linear oscillator with multiplicative noise grows exponentially with time and that this growth may be characterized by a Lyapunov exponent. In section III, we analyze the classical Duffing oscillator in presence of parametric noise. Our technique allows us to study precisely the long-time behavior of the system. In section IV, we consider a particle in an arbitrary confining potential that grows as a polynomial at large distances. In section V, we discuss the case of colored noise where the presence of a new timescale (the correlation time) leads to a nontrivial crossover from the white noise regime to another scaling regime. Our conclusions are presented in section VI. In Appendix A, the nonlinear oscillator in presence of both additive and multiplicative noise is briefly studied: we show that at long times the effect of additive noise is irrelevant. Appendix B is devoted to numerical methods and in Appendix C some useful mathematical relations are recalled.

## II The linear oscillator with parametric noise

In this section we recall known results for an undamped linear oscillator submitted to parametric noise, a generic and widely studied model, in order to understand the role of external multiplicative noise lefever ; bourret ; strato . The dynamical equation for such a system is

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+(\omega ^2+\xi (t))x(t)=0,$$
(1)
where $`x(t)`$ represents the position of the oscillator at time $`t`$ and $`\omega `$ its frequency. The random noise $`\xi (t)`$ is a Gaussian white noise of zero mean-value and of amplitude $`𝒟`$:

$`\xi (t)`$ $`=`$ $`0,`$
$`\xi (t)\xi (t^{})`$ $`=`$ $`𝒟\delta (tt^{}).`$ (2)
The physical interpretation of Eq. (1) is that the frequency of the oscillator is not constant in time but fluctuates around its mean value $`\omega `$ because of randomness in the external conditions (external noise). When these fluctuations are deterministic and periodic in time, Eq. (1) is a Mathieu equation which has been extensively studied landau . Here, we are interested in the case where these fluctuations are random with no deterministic part. The origin, $`x=0`$ and $`\text{d}x/\text{d}t=0`$, is an unstable stationary solution of Eq. (1). As shown in Ref. bourret , this instability can be studied from the dynamical evolution of the Probability Distribution Function $`P(x,v,t)`$ of $`x`$ and $`v`$ (with $`v=\dot{x}=\text{d}x/\text{d}t)`$. This P.D.F. obeys the Fokker-Planck equation vankampen ; risken  associated with Eq. (1):

$$\frac{P}{t}=v\frac{P}{x}+\omega ^2x\frac{P}{v}+\frac{𝒟}{2}\frac{^2}{v^2}\left(x^2P\right),$$
(3)
where Eq. (1) is understood according to Stratonovich rules.

This Fokker-Planck equation leads to a closed system of ordinary differential equations that couple the $`n+1`$ moments of order $`n`$, i.e., moments of the type $`x^{nk}v^k`$, where $`n`$ and $`k`$ are positive integers and $`0kn`$:

$$\frac{\text{d}}{\text{d}t}x^{nk}v^k=(nk)x^{nk1}v^{k+1}\omega ^2kx^{nk+1}v^{k1}+\frac{𝒟}{2}k(k1)x^{nk+2}v^{k2}.$$
(4)
The divergence of the moments with time results from the existence of at least one positive eigenvalue of the linear system (4). In particular, the mean value of the mechanical energy $`E`$ of the system (i.e. the sum of its kinetic and potential energies) grows exponentially with time:

$$E=\frac{1}{2}v^2+\frac{1}{2}\omega ^2x^2e^{\mu t},$$
(5)
where the growth rate $`\mu `$ is the positive real root of the equation:

$$\mu ^3+4\omega ^2\mu =2𝒟.$$
(6)
It has also been proved that the quenched average of the energy $`E`$ grows linearly with time, hence the Lyapunov exponent $`\mathrm{\Lambda }`$, defined as

$$\mathrm{\Lambda }=\underset{t\mathrm{}}{lim}\frac{1}{t}\mathrm{log}E,$$
(7)
is finite and strictly positive hansel ; tessieri . The positivity of the Lyapunov exponent implies the instability of all moments at long times. Note that the growth rate $`\mu `$, defined in Eq. (5), is larger than the Lyapunov exponent because of the convexity inequality, $`\mathrm{log}E\mathrm{log}E`$.

In Figure 1(a), we present the numerical solution of Eq. (1) averaged over a large number of realizations of the noise, where the pulsation is $`\omega =1`$. The algorithm used to solve this stochastic differential equation with multiplicative noise is inspired from mannella  and explained in Appendix B. A numerical estimate of the Lyapunov exponent, given in Fig. 1(a), agrees very well with the analytic expression of Refs. hansel ; tessieri . The usual statistical equipartition of the total energy between kinetic and potential contributions is satisfied: $`E=\omega ^2x^2=v^2`$.

In Figure 1(b), we show the same quantities for the degenerate linear oscillator obtained by taking $`\omega `$ equal to 0. This degenerate case exhibits the same behavior as the generic case and the Lyapunov exponent can be calculated by taking the $`\omega 0`$ limit in the formulas of Refs. hansel ; tessieri . We conclude that the instability triggered by the noise is the dominant effect and that the presence of the linear restoring force $`\omega ^2x`$ is irrelevant.

Hence, in order to avoid an exponential increase of the energy and the amplitude of the oscillator, it is necessary to go beyond the linear approximation and consider the effect of nonlinear restoring forces degli .

## III Duffing oscillator with multiplicative noise

We now analyze the effect of a nonlinear restoring force in an oscillator submitted to an external multiplicative noise. In order to preserve the $`xx`$ symmetry, the nonlinear term has to be odd in the amplitude $`x`$. In this section, we study the particular case of a cubic nonlinearity

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+(\omega ^2+\xi (t))x(t)+\lambda x(t)^3=0.$$
(8)
The coefficient of the nonlinear term is set equal to unity by rescaling the variable $`x(t)`$ to $`x(t)\sqrt{\lambda }`$. The random noise $`\xi `$ is a Gaussian white noise of amplitude $`𝒟`$ as defined in Eq. (2). The deterministic nonlinear mechanical system correponding to Eq. (8) is known as the Duffing oscillator. We shall prove that the cubic term is relevant and prevents the average amplitude from growing exponentially. Instead of an exponential behavior, the average energy of the oscillator as well as the variances of its position and velocity exhibit a power-law behavior with time. We shall calculate exactly the associated scaling exponents.

### III.1 The degenerate cubic oscillator

The linear part of the restoring force, $`\omega ^2x`$, is negligible in comparison to the cubic term when the amplitude of the oscillator is large. In order to study the long-time behavior of the oscillator, we therefore simplify Eq. (8) to that of a degenerate cubic oscillator:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\xi (t)x(t)+x(t)^3=0.$$
(9)
We first study the deterministic part of Eq. (9) and shall add the noise term afterwards roberts . In one dimension, the deterministic cubic oscillator is integrable because the energy $`E`$ defined as

$$E=\frac{1}{2}\dot{x}^2+\frac{1}{4}x^4,$$
(10)
is a conserved quantity. The exact solution of the mechanical system, $`\ddot{x}+x^3=0`$, for a fixed value of $`E`$, is given by

$`x`$ $`=`$ $`E^{1/4}{\displaystyle \frac{\mathrm{sn}((4E)^{\frac{1}{4}}t;\frac{1}{\sqrt{2}})}{\mathrm{dn}((4E)^{\frac{1}{4}}t;\frac{1}{\sqrt{2}})}},`$ (11)
$`\dot{x}`$ $`=`$ $`(2E)^{1/2}{\displaystyle \frac{\mathrm{cn}((4E)^{\frac{1}{4}}t;\frac{1}{\sqrt{2}})}{\mathrm{dn}^2((4E)^{\frac{1}{4}}t;\frac{1}{\sqrt{2}})}},`$ (12)
where sn, cn and dn are the classical Jacobi elliptic functions of modulus $`1/\sqrt{2}`$ abram ; byrd . The properties of elliptic functions that will be useful for our purpose are recalled in Appendix C. From Eqs. (11) and (12), we define the action-angle variables of the cubic oscillator landau ; lichtenberg . The action variable $`I`$ corresponds to the area under a constant energy curve in phase space:

$$I=4_0^{(4E)^{\frac{1}{4}}}\sqrt{2E\frac{x^4}{2}}\text{d}x=8\left(_0^1\sqrt{1u^4}\text{d}u\right)E^{3/4}E^{3/4}.$$
(13)
The angle variable $`\varphi `$, canonically conjugate to the action $`I`$, is equal to $`\omega (E)t`$ (but for an unimportant additive constant) where $`\omega (E)`$ is the frequency corresponding to the energy $`E`$. From Eqs. (11) and (12), the phase is identified as

$$\varphi =(4E)^{\frac{1}{4}}t.$$
(14)
This phase variable is defined modulo the common real period $`4K`$ of the elliptic functions $`\mathrm{sn},\mathrm{cn},\mathrm{dn}`$ that appear in Eqs. (11) and (12). The quarter of the period is given by

$$K=K(\frac{1}{\sqrt{2}})=\sqrt{2}_0^1\frac{\text{d}u}{\sqrt{1u^4}}1.854.$$
(15)
In terms of the energy-angle coordinates $`(E,\varphi )`$, the original variables $`(x,\dot{x})`$ read as follows:

$$x=E^{1/4}\mathrm{sd}(\varphi ;1/\sqrt{2})\text{ and }\dot{x}=(2E)^{1/2}\frac{\mathrm{cn}(\varphi ;1/\sqrt{2})}{\mathrm{dn}^2(\varphi ;1/\sqrt{2})},$$
(16)
where we have introduced the function sd = sn/dn. Hence, for the free system (without noise) the second-order dynamical equation, $`\ddot{x}+x^3=0`$, is equivalent to the following two first-order equations, the first one representing energy conservation:

$`\dot{E}`$ $`=`$ $`0,`$
$`\dot{\varphi }`$ $`=`$ $`(4E)^{\frac{1}{4}}.`$ (17)
The presence of external noise spoils the integrability of the dynamical system (9) and causes $`E`$ to grow with time by continuously injecting energy into the system. We now take into account the external noise and rewrite the system (17) in $`(E,\varphi )`$ coordinates. The stochastic evolution equation for the energy is given by

$$\dot{E}=x\dot{x}\xi (t)=\sqrt{2}E^{3/4}\frac{\mathrm{sn}(\varphi ;1/\sqrt{2})\mathrm{cn}(\varphi ;1/\sqrt{2})}{\mathrm{dn}^3(\varphi ;1/\sqrt{2})}\xi (t).$$
(18)
Inverting the relation (16) between $`\varphi `$ and $`x`$,

$$\varphi =\mathrm{sd}^1(\frac{x}{E^{1/4}},\frac{1}{\sqrt{2}}),$$
(19)
and using Eqs. (16) and (18), we obtain the stochastic evolution of the phase variable:

$$\dot{\varphi }=(4E)^{\frac{1}{4}}\frac{\mathrm{sd}^2(\varphi ;1/\sqrt{2})}{2\sqrt{2}E^{1/4}}\xi (t).$$
(20)
With the help of the auxiliary variable $`\mathrm{\Omega }`$ defined as

$$\mathrm{\Omega }=2\sqrt{2}E^{1/4},$$
(21)
we derive a compact and symmetric form for the two stochastic evolution equations (18) and (20):

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`{\displaystyle \frac{\mathrm{sn}(\varphi )\mathrm{cn}(\varphi )}{\mathrm{dn}^3(\varphi )}}\xi (t),`$ (22)
$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{2}}{\displaystyle \frac{\mathrm{sd}^2(\varphi )}{\mathrm{\Omega }}}\xi (t).`$ (23)
where the elliptic modulus $`1/\sqrt{2}`$ common to all the elliptic functions has been omitted. We emphasize that the coupled equations (22) and (23) are mathematically equivalent to the initial system and have been derived without any approximation. Moreover, the nature of the parametric perturbation has played no role in the derivation: the function $`\xi (t)`$ can be a deterministic or a stochastic function with arbitrary statistical properties.

We now perform a precise analysis of the long-time behavior of the cubic oscillator driven by a multiplicative Gaussian white noise. From a heuristic point of view, we observe from Eq. (22) that $`\mathrm{\Omega }`$ undergoes a diffusion process and should scale typically as $`t^{1/2}`$. We also notice from Eq. (23) that, as $`\mathrm{\Omega }`$ grows, the phase $`\varphi `$ varies more and more rapidly with time. Hence the phase $`\varphi `$ is a fast variable and it is natural to average the dynamics over its rapid variations. This averaging process leads to some remarkable and general identities between different physical quantitites. Thus, we obtain the average of $`\dot{x}^2`$ from Eq. (16),

$$\dot{x}^2=2\overline{\frac{\mathrm{cn}^2(\varphi )}{\mathrm{dn}^4(\varphi )}}E,$$
(24)
where the overline indicates the mean-value over a period $`4K`$. The coefficient in Eq. (24) can be calculated explicitly:

$$2\overline{\frac{\mathrm{cn}^2(\varphi )}{\mathrm{dn}^4(\varphi )}}=\frac{4}{3};$$
(25)
the proof of this identity (and of some forthcoming similar identities) is presented in Appendix C. Equation (24) is nothing but a statistical equipartition  relation between the mean values of the total energy and kinetic energy:

$$E=\frac{3}{4}\dot{x}^2.$$
(26)
Such an identity provides a generalization of the usual equipartition theorem valid for quadratic degrees of freedom. Another statistical equality follows from Eq. (26) and the definition (10) of the energy $`E:`$

$$\dot{x}^2=x^4.$$
(27)
We emphasize that these equalities are ‘universal’ in the sense that they are independent of the form of the noise we consider. The derivation of these equipartition relations relies only on the hypothesis that the distribution of the angle $`\varphi `$ is uniform over the interval $`[0,4K]`$ when $`t\mathrm{}`$. In particular, identities (26) and (27) are valid for multiplicative as well as for additive noise (see Appendix A). Our numerical results (Fig. 2) confirm these relations and that the averaging procedure is therefore valid.

The same averaging procedure allows us to derive a closed equation for the stochastic evolution of the slow variable $`\mathrm{\Omega }.`$ We start by writing the Fokker-Planck equation governing the evolution of the P.D.F. $`P_t(\mathrm{\Omega },\varphi )`$ associated with the system (22,23). Since the noise term appears as a multiplicative factor, one must be cautious about the convention used to define stochastic calculus. Here as well as in the following, we shall use Stratonovich rules because they are obtained naturally when white noise is considered as a limit of colored noise with vanishing correlation time vankampen ; risken . The Fokker-Planck equation corresponding to Eqs. (22) and (23) reads:

$`_tP`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{2}}_\varphi P`$ (28)
$`+`$ $`{\displaystyle \frac{𝒟}{2}}\left(_\varphi \left({\displaystyle \frac{g(\varphi )}{\mathrm{\Omega }}}_\varphi {\displaystyle \frac{g(\varphi )}{\mathrm{\Omega }}}P\right)_\varphi \left({\displaystyle \frac{g(\varphi )}{\mathrm{\Omega }}}_\mathrm{\Omega }Pf(\varphi )\right)_\mathrm{\Omega }\left(f(\varphi )_\varphi {\displaystyle \frac{g(\varphi )}{\mathrm{\Omega }}}P\right)+_\mathrm{\Omega }\left(f(\varphi )_\mathrm{\Omega }f(\varphi )P\right)\right),`$
where we have defined:

$$f(\varphi )=\frac{\mathrm{sn}(\varphi )\mathrm{cn}(\varphi )}{\mathrm{dn}^3(\varphi )}\text{ and }g(\varphi )=\mathrm{sd}^2(\varphi ).$$
(29)
The Fokker-Planck equation (28) written in the variables $`(\mathrm{\Omega },\varphi )`$ is exact because we study the case of a Gaussian white noise. In order to pursue our calculations, we assume that $`P_t(\mathrm{\Omega },\varphi )`$ becomes independent of $`\varphi `$ when $`t\mathrm{}`$, i.e., that the probability measure for $`\varphi `$ is uniform over the interval $`[0,4K]`$. We now average the Fokker-Planck equation (28) over the angular variable. We shall use the fact that the average of the derivative of any function is zero:

$$\overline{_\varphi (\mathrm{})}=0.$$
(30)
This implies that:

$$\overline{_\varphi \left(g(\varphi )_\varphi g(\varphi )\right)}=0\text{ and }\overline{f(\varphi )}=0\text{ because }f(\varphi )=\frac{1}{2}_\varphi g(\varphi ).$$
(31)
Using these properties, and in particular the last identity in (31), we derive the phase-averaged Fokker-Planck equation:

$$_t\stackrel{~}{P}=\frac{𝒟}{2}\left(_\mathrm{\Omega }^2\stackrel{~}{P}2_\mathrm{\Omega }\frac{\stackrel{~}{P}}{\mathrm{\Omega }}\right),$$
(32)
where $`\stackrel{~}{P}_t(\mathrm{\Omega })`$ is now a function of $`\mathrm{\Omega }`$ and $`t`$ only, and where $``$ is given by

$$=\overline{f^2(\varphi )}=\frac{1}{K}_0^K\frac{\mathrm{sn}^2(\varphi )\mathrm{cn}^2(\varphi )}{\mathrm{dn}^6(\varphi )}\text{d}\varphi =2\frac{_0^1u^2\sqrt{1u^4}\text{d}u}{_0^1\frac{\text{d}u}{\sqrt{1u^4}}}=\frac{96\pi ^2}{15\left(\mathrm{\Gamma }(\frac{1}{4})\right)^4}0.3655\mathrm{}$$
(33)
The third equality is obtained by setting $`u=\mathrm{sd}(\varphi )/\sqrt{2}`$; the expression in terms of the Euler $`\mathrm{\Gamma }`$ function can be found in abram . The averaged Fokker-Planck equation (32) corresponds to the following effective Langevin dynamics for the variable $`\mathrm{\Omega }`$,

$$\dot{\mathrm{\Omega }}=\frac{𝒟}{\mathrm{\Omega }}+\mathrm{\Xi }(t),$$
(34)
where the effective noise $`\mathrm{\Xi }(t)`$ satisfies the relation

$$\mathrm{\Xi }(t)\mathrm{\Xi }(t^{})=𝒟\delta (tt^{}).$$
(35)
We conclude from Eq. (34) that $`\mathrm{\Omega }`$ exhibits a diffusive behavior. The dynamics of all moments of $`\mathrm{\Omega }`$ can be deduced from Eq. (32), as for example:

$`_t\mathrm{\Omega }^2`$ $`=`$ $`3𝒟,`$ (36)
$`_t\mathrm{\Omega }^4`$ $`=`$ $`10𝒟\mathrm{\Omega }^2.`$ (37)
From Eq. (36), we calculate the effective diffusion constant $`D_\mathrm{\Omega }`$ of $`\mathrm{\Omega }`$:

$$D_\mathrm{\Omega }=3𝒟1.096𝒟,$$
(38)
where $`\mathrm{\Omega }^2D_\mathrm{\Omega }t`$ when $`t\mathrm{}`$. We conclude from Eq. (37) that $`\mathrm{\Omega }`$ has a non-trivial flatness:

$$\frac{\mathrm{\Omega }^4}{\mathrm{\Omega }^2^2}=\frac{5}{3}.$$
(39)
This last equality shows that although the variable $`\mathrm{\Omega }`$ diffuses as $`t^{1/2}`$, it is not Gaussian (and hence it is not a Brownian variable).

From this statistical information about $`\mathrm{\Omega }`$, we now derive the scaling with time of the average of $`E`$, $`x^2`$, $`v^2`$ and calculate the numerical prefactors (generalized diffusion constants). From Eqs. (21), (38) and (39), we deduce that:

$$E=\frac{\mathrm{\Omega }^4}{64}=\frac{5}{3}\frac{\mathrm{\Omega }^2^2}{64}=\frac{15}{64}^2𝒟^2t^2=\frac{48\pi ^4}{5\left(\mathrm{\Gamma }(\frac{1}{4})\right)^8}𝒟^2t^20.0313𝒟^2t^2.$$
(40)
Using the equipartition relations (26) and (27), we find that:

$`\dot{x}^2`$ $`=`$ $`{\displaystyle \frac{4}{3}}E={\displaystyle \frac{64\pi ^4}{5\left(\mathrm{\Gamma }(\frac{1}{4})\right)^8}}𝒟^2t^20.0417𝒟^2t^2,`$ (41)
$`x^2`$ $`=`$ $`{\displaystyle \frac{\overline{\mathrm{sd}^2(\varphi )}}{8}}\mathrm{\Omega }^2={\displaystyle \frac{15}{16}}^2𝒟t={\displaystyle \frac{192\pi ^4}{5\left(\mathrm{\Gamma }(\frac{1}{4})\right)^8}}𝒟t0.125𝒟t.`$ (42)
The prefactor in the mean value of the energy, Eq. (40), agrees with that of lind2 . We have verified numerically the power-laws predicted by Eqs. (40), (41) and (42) and determined the prefactors. The results are displayed in Fig. 3. The numerical values for the exponents and the prefactors are in perfect agreement with the analytical calculations: $`D_E=0.0313𝒟^2`$, $`D_x=0.125𝒟`$ and $`D_v=0.0417𝒟^2`$.

The averaged distribution fonction $`\stackrel{~}{P}_t(\mathrm{\Omega })`$ can be calculated because Eq. (32) is exactly solvable due to its invariance under rescalings $`t\lambda ^2t`$ , $`\mathrm{\Omega }\lambda \mathrm{\Omega },`$ $`\lambda `$ being an arbitrary real number (this invariance is the same as that of the heat equation). Equation (32) is solved by using the self-similar Ansatz:

$$\stackrel{~}{P}_t(\mathrm{\Omega })=\frac{1}{\sqrt{t}}\mathrm{\Pi }(\frac{\mathrm{\Omega }}{\sqrt{t}}),$$

and the P.D.F of $`\mathrm{\Omega }`$ is found to be

$$\stackrel{~}{P}_t(\mathrm{\Omega })=\sqrt{\frac{2}{\pi }}\mathrm{\Omega }^2\frac{\mathrm{e}^{\mathrm{\Omega }^2/(2𝒟t)}}{(𝒟t)^{\frac{3}{2}}}.$$
(43)
This distribution function allows us to calculate the P.D.F. in the $`(x,\dot{x})`$ variables in the time-asymptotic regime. In particular, the skewness and the flatness factors of the position and of the velocity can be calculated analytically. Since both variables $`x`$ and $`\dot{x}`$ are parity symmetric, their skewness vanishes. The flatness is given by

$`{\displaystyle \frac{x^4}{x^2^2}}`$ $`=`$ $`{\displaystyle \frac{4}{3\left(\overline{\mathrm{sd}^2(\varphi )}\right)^2}}{\displaystyle \frac{\mathrm{\Omega }^4}{\mathrm{\Omega }^2^2}}={\displaystyle \frac{16}{45^2}}2.66,`$ (44)
$`{\displaystyle \frac{v^4}{v^2^2}}`$ $`=`$ $`{\displaystyle \frac{9}{4}}\overline{{\displaystyle \frac{\mathrm{cn}^4(\varphi )}{\mathrm{dn}^8(\varphi )}}}{\displaystyle \frac{\mathrm{\Omega }^8}{\mathrm{\Omega }^4^2}}={\displaystyle \frac{27}{5}}=5.4.`$ (45)
These values are in excellent agreement with the numerical computations shown in Fig. 4. We notice that the variables $`x`$ and $`\dot{x}`$ are non-Gaussian because their flatness differs from the Gaussian value 3.

To summarize, the following scalings have been derived:

$`E`$ $``$ $`t^2,`$
$`x`$ $``$ $`t^{\frac{1}{2}},`$
$`\dot{x}`$ $``$ $`t.`$ (46)
In particular, we conclude that the amplitude $`x`$ of a degenerate cubic oscillator does not grow exponentially with time but has a normal diffusive behavior.

### III.2 The Duffing oscillator

We now study the general case of a non-zero pulsation $`\omega `$:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+(\omega ^2+\xi (t))x(t)+x(t)^3=0.$$
(47)
Here the coefficient of the nonlinear term is taken to be unity and the random noise is Gaussian and white as defined in Eq. (2).

The results of sections II and III.1 show two regimes: starting from a small initial condition, the amplitude of the oscillator grows exponentially with time until $`x\omega `$ where the linear and nonlinear terms are of the same order and then the amplitude grows as the square-root of time according to Eq. (46). Because the deterministic system corresponding to Eq. (47) is integrable with the help of elliptic functions, this crossover from exponential to algebraic can be derived in a quantitative manner.

When $`\omega `$ is non zero, Eqs. (11) and (12) become:

$`x`$ $`=`$ $`\left({\displaystyle \frac{4E^2}{4E+\omega ^4}}\right)^{1/4}{\displaystyle \frac{\mathrm{sn}((4E+\omega ^4)^{\frac{1}{4}}t;k)}{\mathrm{dn}((4E+\omega ^4)^{\frac{1}{4}}t;k)}},`$ (48)
$`\dot{x}`$ $`=`$ $`(2E)^{1/2}{\displaystyle \frac{\mathrm{cn}((4E+\omega ^4)^{\frac{1}{4}}t;k)}{\mathrm{dn}^2((4E+\omega ^4)^{\frac{1}{4}}t;k)}},`$ (49)
where the elliptic modulus $`k`$ varies with the energy and is given by

$$k^2=\frac{\sqrt{4E+\omega ^4}\omega ^2}{2\sqrt{4E+\omega ^4}}.$$
(50)
We notice that $`k`$ tends to the limiting value $`1/\sqrt{2}`$ when the energy goes to infinity. Defining the angle variable as:

$$\varphi =\mathrm{sd}^1(\frac{(4E+\omega ^4)^{\frac{1}{4}}}{\sqrt{2E}}x;k),$$
(51)
we rewrite the dynamical equation in energy-angle coordinates. However, while deriving the dynamical equation for $`\varphi `$ we must remember that the elliptic modulus $`k`$ depends on the energy $`E`$ and is, therefore, a function of time. After reintroducing the multiplicative noise term, the stochastic Duffing oscillator in the energy-angle variables becomes:

$`\dot{E}`$ $`=`$ $`{\displaystyle \frac{2E}{(4E+\omega ^4)^{\frac{1}{4}}}}{\displaystyle \frac{\mathrm{sn}(\varphi ;k)\mathrm{cn}(\varphi ;k)}{\mathrm{dn}^3(\varphi ;k)}}\xi (t);`$ (52)
$`\dot{\varphi }`$ $`=`$ $`(4E+\omega ^4)^{\frac{1}{4}}\xi (t){\displaystyle \frac{(2E+\omega ^4)\mathrm{sd}^2(\varphi ;k)}{(4E+\omega ^4)^{\frac{5}{4}}}}\omega ^2\xi (t)\left({\displaystyle \frac{E\mathrm{sd}^4(\varphi ;k)}{(4E+\omega ^4)^{\frac{7}{4}}}}{\displaystyle \frac{E\mathrm{sn}(\varphi ;k)\mathrm{cn}(\varphi ;k)_0^\varphi \mathrm{sd}^2(\theta ;k)𝑑\theta }{(4E+\omega ^4)^{\frac{7}{4}}\mathrm{dn}^3(\varphi ;k)}}\right).`$ (53)
As compared to Eq. (20), two supplementary terms appear in Eq. (53). These terms are related to $`dk/dt`$ and are proportional to $`\omega ^2`$.

Although the dynamical equations are more complicated than those of the purely cubic case, the analysis can be performed as above. We shall however simplify our discussion here by taking $`k`$ equal to its asymptotic value $`1/\sqrt{2}`$. This approximation is justified as soon as the energy is large. We also replace the noise $`\xi (t)`$ by the effective noise $`\mathrm{\Xi }(t)`$ defined in Eq. (35). This second approximation is only qualitatively correct since it amounts to neglecting a deterministic force in the effective Langevin dynamics for $`E`$. We thus obtain:

$$\dot{E}\frac{2E}{(4E+\omega ^4)^{\frac{1}{4}}}\mathrm{\Xi }(t).$$
(54)
We deduce from Eq. (54) that as long as $`E\omega ^4`$, the energy behaves as the exponential of a Brownian motion and, therefore, increases exponentially with time. However, when $`E>\omega ^4`$, the nonlinear term becomes important, Eq. (52) reduces to Eq. (18), and the energy grows as the square of time.

We expect that the crossover from exponential to algebraic growth will appear when $`E\omega ^4`$, or $`x\omega `$. Using unscaled variables, the balance between linear and non-linear terms in Eq. (8) is obtained when $`x=x_c\omega /\sqrt{\lambda }`$. Fig. 5 demonstrates that the two regimes are observed numerically when the nonlinear coefficient $`\lambda `$ is very small compared to $`\omega ^2`$: we use the numerical values $`\omega =1`$ and $`\lambda =10^6`$. We notice that in the short-time linear regime, the usual equipartition relation for a quadratic potential is verified ($`E\dot{x}^2`$), while the exponential growth of the energy is characterized by the Lyapunov exponent $`\mathrm{\Lambda }`$ defined in Eq. (7). In the long-time regime, the equipartition ratio reaches its nonlinear value $`3/4`$, while the energy growth becomes algebraic, with a generalized diffusion constant $`D_E(\lambda )`$ in good agreement with Eq. (40), up to the expected scaling factor: $`D_E(\lambda )=D_E/\lambda `$.

## IV The general nonlinear oscillator with parametric noise

We now consider the case of a particle trapped in a confining potential $`𝒰(x)`$ and subject to an external noise. This mechanical system generalizes the Duffing oscillator studied in the previous section and its dynamics is given by

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\xi (t)x(t)=\frac{𝒰(x)}{x},$$
(55)
where $`\xi (t)`$ is the Gaussian white noise of Eq. (2). For the potential to be confining, we must have $`𝒰+\mathrm{}`$ when $`|x|\mathrm{}`$. We restrict our analysis to the case where $`𝒰`$ is a polynomial; in order to respect the $`xx`$ symmetry $`𝒰`$ is even in $`x`$. Hence, when $`|x|\mathrm{}`$,

$$𝒰\frac{x^{2n}}{2n}\text{ with }n2,$$
(56)
the coefficient of $`x^{2n}`$ is chosen to be $`1/2n`$ by a suitable rescaling. As before, we expect the amplitude of the oscillator to grow without bounds at large times. Keeping only the relevant terms in Eq. (55) leads to

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\xi (t)x(t)+x(t)^{2n1}=0.$$
(57)
The deterministic version of Eq. (57), $`\ddot{x}+x^{2n1}=0`$, is integrable because of energy conservation. For a given value of $`E`$ the deterministic solution is

$`x`$ $`=`$ $`E^{1/2n}𝒮_n\left((2nE)^{\frac{n1}{2n}}t\right),`$ (58)
$`\dot{x}`$ $`=`$ $`(2n)^{\frac{n1}{2n}}E^{1/2}𝒮_n^{}\left((2nE)^{\frac{n1}{2n}}t\right),`$ (59)
where the energy reads

$$E=\frac{1}{2}\dot{x}^2+\frac{1}{2n}x^{2n}.$$
(60)
The function $`𝒮_n`$ is defined as the inverse function of an hyperelliptic integral as follows

$$𝒮_n(X)=YX=\sqrt{n}_0^{\frac{Y}{(2n)^{1/2n}}}\frac{\text{d}u}{\sqrt{1u^{2n}}}=\frac{\sqrt{n}}{(2n)^{1/2n}}_0^Y\frac{\text{d}u}{\sqrt{1\frac{u^{2n}}{2n}}}.$$
(61)
From this definition, we find a relation between $`𝒮_n`$ and its derivative $`𝒮_n^{}`$

$$𝒮_n^{}(X)=\frac{(2n)^{\frac{1}{2n}}}{\sqrt{n}}\left(1\frac{(𝒮_n(X))^{2n}}{2n}\right)^{\frac{1}{2}}.$$
(62)
The action-angle variables for the deterministic system are given by

$`I`$ $`=`$ $`4{\displaystyle _0^{(2nE)^{\frac{1}{2n}}}}\sqrt{2E{\displaystyle \frac{x^{2n}}{n}}}\text{d}x=4(2^{n+1}n)^{\frac{1}{2n}}E^{\frac{n+1}{2n}}{\displaystyle _0^1}\sqrt{1u^{2n}}\text{d}uE^{\frac{n+1}{2n}},`$ (63)
$`\varphi `$ $`=`$ $`𝒮_n^1\left({\displaystyle \frac{x}{E^{\frac{1}{2n}}}}\right)=\sqrt{n}{\displaystyle _0^{\frac{x}{(2nE)^{1/2n}}}}{\displaystyle \frac{\text{d}u}{\sqrt{1u^{2n}}}}.`$ (64)
The angle variable $`\varphi `$ is well defined modulo the period $`4K_n`$ of the function $`𝒮_n`$ where

$$K_n=\sqrt{n}_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}.$$
(65)
When the noise term is taken into account, the energy is not conserved and the stochastic evolution of the variables $`E`$ and $`\varphi `$ is given by

$`\dot{E}`$ $`=`$ $`x\dot{x}\xi (t)=(2n)^{\frac{n1}{2n}}E^{\frac{n+1}{2n}}𝒮_n(\varphi )𝒮_n^{}(\varphi )\xi (t),`$ (66)
$`\dot{\varphi }`$ $`=`$ $`(2nE)^{\frac{n1}{2n}}{\displaystyle \frac{1}{(2n)^{\frac{1}{n}}}}{\displaystyle \frac{𝒮_n(\varphi )^2}{(2nE)^{\frac{n1}{2n}}}}\xi (t).`$ (67)
Introducing the auxiliary variable $`\mathrm{\Omega }`$,

$$\mathrm{\Omega }=(2n)^{\frac{n+1}{2n}}E^{\frac{n1}{2n}},$$
(68)
the Eqs. (66) and (67) can be written in the simpler form:

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`(n1)𝒮_n(\varphi )𝒮_n^{}(\varphi )\xi (t),`$ (69)
$`\dot{\varphi }`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }}{(2n)^{\frac{1}{n}}}}{\displaystyle \frac{𝒮_n(\varphi )^2}{\mathrm{\Omega }}}\xi (t).`$ (70)
As in the case of the Duffing oscillator, we observe that $`\varphi `$ is a fast variable and therefore we shall average out its rapid variations. Starting from (59), we write

$$\dot{x}^2=(2n)^{\frac{n1}{n}}\overline{𝒮_n^{}(\varphi )^2}E=2\frac{_0^1\text{d}u\sqrt{1u^{2n}}}{_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}}E,$$
(71)
the last equality is derived by writing $`u=𝒮_n(\varphi )`$, and using Eqs. (62) and (65). Moreover, the following identity is true (as can be shown by integrating $`_0^11.\sqrt{1u^{2n}}\text{d}u`$ by parts):

$$_0^1\text{d}u\sqrt{1u^{2n}}=n_0^1\text{d}u\frac{u^{2n}}{\sqrt{1u^{2n}}}=n_0^1\text{d}u\sqrt{1u^{2n}}+n_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}.$$
(72)
Substituting this identity in Eq. (71) leads to

$$E=\frac{n+1}{2n}\dot{x}^2.$$
(73)
From the definition (60) of the energy, we derive another statistical equality:

$$\dot{x}^2=x^{2n}.$$
(74)
We emphasize again that these generalized equipartition relations are valid for any type of noise, colored or white, additive or multiplicative. The only hypothesis is that the probability distribution function $`P_t(\mathrm{\Omega },\varphi )`$ becomes uniform in $`\varphi `$ over the interval $`[0,4K_n]`$ when $`t\mathrm{}`$. We observe from the numerical simulations presented in Fig. 6.c that this condition is very well satisfied.

The averaging process with respect to the fast variable generates an effective Langevin dynamics for the slow variable $`\mathrm{\Omega }(t)`$. Starting from the Fokker-Planck equation for the total P.D.F., $`P_t(\mathrm{\Omega },\varphi ),`$ and averaging over $`\varphi `$ leads to the following evolution equation for the averaged probability distribution $`\stackrel{~}{P}_t(\mathrm{\Omega })`$

$$_t\stackrel{~}{P}=\frac{_n𝒟}{2}\left(_\mathrm{\Omega }^2\stackrel{~}{P}\frac{2}{n1}_\mathrm{\Omega }\frac{\stackrel{~}{P}}{\mathrm{\Omega }}\right),$$
(75)
where $`_n`$ is given by

$$_n=\frac{(n1)^2}{n}(2n)^{\frac{2}{n}}\frac{_0^1\text{d}uu^2\sqrt{1u^{2n}}}{_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}}=\frac{(n1)^2}{n+1}(2n)^{\frac{2}{n}}\frac{\mathrm{\Gamma }(\frac{3}{2n})}{\mathrm{\Gamma }(\frac{1}{2n})}\frac{\mathrm{\Gamma }(\frac{3n+1}{2n})}{\mathrm{\Gamma }(\frac{3n+3}{2n})},$$
(76)
$`\mathrm{\Gamma }(.)`$ being the Euler Gamma function. The effective Langevin dynamics for the variable $`\mathrm{\Omega }`$ is thus

$$\dot{\mathrm{\Omega }}=\frac{_n𝒟}{n1}\frac{1}{\mathrm{\Omega }}+\mathrm{\Xi }_n(t),$$
(77)
where the effective Gaussian white noise $`\mathrm{\Xi }_n(t)`$ satisfies the relation

$$\mathrm{\Xi }_n(t)\mathrm{\Xi }_n(t^{})=_n𝒟\delta (tt^{}).$$
(78)
The averaged distribution fonction $`\stackrel{~}{P}_t(\mathrm{\Omega })`$ can be calculated because Eq. (75) is exactly solvable due to its invariance under rescalings $`t\lambda ^2t`$ , $`\mathrm{\Omega }\lambda \mathrm{\Omega },`$ $`\lambda `$ being an arbitrary real number. The P.D.F. of $`\mathrm{\Omega }`$ is found to be

$$\stackrel{~}{P}_t(\mathrm{\Omega })=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2(n1)}\right)}\frac{\mathrm{\Omega }^{\frac{2}{n1}}}{\left(2_n𝒟t\right)^{\frac{n+1}{2(n1)}}}\mathrm{exp}\left\{\frac{\mathrm{\Omega }^2}{2_n𝒟t}\right\},$$
(79)
from which we obtain the P.D.F. of the energy:

$$\stackrel{~}{P}_t(E)=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2(n1)}\right)}\frac{n1}{nE}\left(\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{2_n𝒟t}\right)^{\frac{n+1}{2(n1)}}\mathrm{exp}\left\{\frac{(2n)^{\frac{n+1}{n}}E^{\frac{n1}{n}}}{2_n𝒟t}\right\}.$$
(80)
The long-time behavior of the amplitude, velocity and energy of the general nonlinear oscillator can now be derived. Using the equipartition identity (73) and Eq. (80), we obtain

$`E`$ $`=`$ $`{\displaystyle \frac{1}{(2n)^{\frac{n+1}{n1}}}}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3n+1}{2n2}\right)}{\mathrm{\Gamma }\left(\frac{n+1}{2n2}\right)}}\left(2𝒟_nt\right)^{\frac{n}{n1}}t^{\frac{n}{n1}},`$ (81)
$`\dot{x}^2`$ $`=`$ $`{\displaystyle \frac{2n}{n+1}}Et^{\frac{n}{n1}}.`$ (82)
Using Eq (76), we find that

$$\overline{𝒮_n^2(\varphi )}=\frac{1}{K_n}_0^{K_n}𝒮_n^2(\varphi )\text{d}\varphi =(2n)^{\frac{1}{n}}\frac{_0^1\frac{u^2\text{d}u}{\sqrt{1u^{2n}}}}{_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}}=\frac{n+3}{(2n)^{\frac{1}{n}}(n1)^2}_n,$$
(83)
where we made the change of variables $`u=𝒮_n(\varphi ),`$ and used the following identity (obtained by integrating by parts):

$$_0^1u^2\sqrt{1u^{2n}}\text{d}u=\frac{n}{n+3}_0^1\frac{u^2\text{d}u}{\sqrt{1u^{2n}}}.$$
(84)
Finally, we deduce from Eqs. (58), (80) and (83)

$$x^2=\overline{𝒮_n^2(\varphi )}E^{1/n}=\frac{2}{n1}\frac{_n}{(2n)^{\frac{1}{n}}}\frac{\mathrm{\Gamma }\left(\frac{3n+1}{2n2}\right)}{\mathrm{\Gamma }\left(\frac{n+1}{2n2}\right)}\left(\frac{2_n𝒟t}{(2n)^{\frac{n+1}{n}}}\right)^{\frac{1}{n1}}t^{\frac{1}{n1}}.$$
(85)
The analytical results for the nonlinear oscillators with quintic $`x^5`$ ($`n=3`$) and heptic nonlinearity $`x^7`$ ($`n=4`$) are as follows:

$`\text{For }n=3,`$ $`E=0.097t^{3/2},\dot{x}^2=0.145t^{3/2},x^2=0.290t^{1/2}.`$ (86)
$`\text{For }n=4,`$ $`E=0.130t^{4/3},\dot{x}^2=0.208t^{4/3},x^2=0.347t^{1/3}.`$ (87)
The formulas (81), (82) and (85) were verified numerically. The scaling exponents and the prefactors given in Eqs. (86) and (87) are in excellent agreement with the numerical values, as shown in Fig. 6.

In conclusion, we have derived the following scaling relations:

$`E`$ $``$ $`t^{\frac{n}{n1}},`$
$`x`$ $``$ $`t^{\frac{1}{2(n1)}},`$
$`\dot{x}`$ $``$ $`t^{\frac{n}{2(n1)}}.`$ (88)
In particular, it should be remarked that $`x`$ undergoes an anomalous diffusion with time with exponent $`1/(2n2)`$. If we make $`n1`$ formally, this exponent diverges to infinity: this is consistent with the exponential growth of the linear oscillator (see Sec. II).

We end this section by considering the case of a general confining potential energy $`𝒰`$ neither necessarily polynomial in $`x`$ nor even in $`x`$. The only requirement is that $`𝒰+\mathrm{}`$ when $`|x|\mathrm{}`$. We discuss the qualitative behavior of $`E`$, $`\dot{x}`$ and $`x`$ at large times from elementary scaling considerations. Suppose first that $`𝒰|x|^r`$ for large values of $`|x|`$, $`r`$ being an arbitrary real number.

* If $`r>2`$, then balance between kinetic and potential energies leads to $`E\dot{x}^2x^r`$; thus the time evolution of the energy is given by $`\dot{E}x\dot{x}\xi E^{\frac{1}{r}\frac{1}{2}}\xi `$. From the scaling relations $`\dot{E}E/t`$ and $`\xi t^{1/2}`$, we conclude that

$$Et^{\frac{r}{r2}},xt^{\frac{1}{r2}},\dot{x}t^{\frac{r}{2(r2)}}.$$

This qualitative argument can be made rigorous by generalizing the results obtained above to non-integer values of $`r`$.
* If $`r2`$, then the potential $`𝒰`$ is negligible with respect to the multiplicative noise term and we are back to the case of the degenerate linear oscillator. Therefore $`E`$, $`\dot{x}`$ and $`x`$ grow exponentially with time.

If the potential grows exponentially, i.e. $`𝒰\mathrm{e}^{x^\beta }`$, $`\beta `$ being a positive real number, then similar considerations lead to $`Et`$ and $`\dot{x}t^{1/2}`$ (leaving aside logarithmic corrections). We then conjecture that the amplitude $`x`$ diffuses in a logarithmically slow manner: $`x(\mathrm{ln}t)^{1/\beta }`$.

## V Colored Gaussian noise

We now consider the case where the Gaussian noise has a non-zero correlation time, and discuss how the previously found scalings are modified. The system we want to study satisfies the dynamical equation:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)\eta (t)+x(t)^{2n1}=0,$$
(89)
where $`\eta `$ is a colored Gaussian noise of zero mean-value. The statistical properties of $`\eta `$ are determined by

$`\eta (t)`$ $`=`$ $`0,`$
$`\eta (t)\eta (t^{})`$ $`=`$ $`{\displaystyle \frac{𝒟}{2\tau }}\mathrm{e}^{|tt^{}|/\tau },`$ (90)
where $`\tau `$ is the correlation time of the noise. The noise $`\eta `$ can be obtained from white noise by solving the Ornstein-Uhlenbeck equation:

$$\frac{\text{d}\eta (t)}{\text{d}t}=\frac{1}{\tau }\eta (t)+\frac{1}{\tau }\xi (t),$$
(91)
where $`\xi (t)`$ is the Gaussian white noise defined in Eq. (2), and $`t,t^{}\tau `$.

Introducing action-angle variables as in Sec. IV, we obtain the same set of coupled Langevin equations, Eqs. (69,70), where $`\xi (t)`$ is replaced by the colored noise $`\eta (t)`$. As emphasized previously, the generalized equipartition relations are independent of the nature of the noise: Eqs. (73) and (74) remain valid when the noise is correlated in time. This is confirmed by numerical simulations (see Fig. 7.b for $`n=2`$).

The scalings found in sections III and IV are deduced by averaging the Fokker-Planck equation. Here, we must write the evolution equation for the joint P.D.F. of $`x`$, $`v=\dot{x}`$ and $`\eta `$, $`P_t(x,v,\eta )`$:

$$\frac{P}{t}=v\frac{P}{x}+(x^{2n1}x\eta )\frac{P}{v}+\frac{1}{\tau }\frac{\eta P}{\eta }+\frac{𝒟}{2\tau ^2}\frac{^2P}{\eta ^2}.$$
(92)
We perform a scaling analysis of this equation in the spirit of yp . Balancing the diffusion term with the time derivative leads to $`\eta t^{1/2}`$. Then we compare the terms of probability current $`\frac{vP}{x}`$ and $`\frac{(x^{2n1}x\eta )P}{v}.`$ A consistent balance between these terms is possible only if $`v^2x^{2n}`$ and $`x^{2n2}\eta `$. We thus find the following scaling relations

$`E`$ $``$ $`t^{\frac{n}{2(n1)}},`$
$`x`$ $``$ $`t^{\frac{1}{4(n1)}},`$
$`\dot{x}`$ $``$ $`t^{\frac{n}{4(n1)}}.`$ (93)
Thus, we predict that the scaling exponents for colored noise are half the exponents calculated for white noise (88). Numerical simulations indeed confirm that the average energy of a cubic oscillator ($`n=2`$) with colored multiplicative noise grows linearly with time (see Fig. 7.a).

The period $`T`$ of a deterministic oscillator (without noise) decreases as the energy increases: $`TE^{\frac{n1}{2n}}`$ from Eq. (67). When the equations are written in the energy-angle coordinates, two time scales $`T`$ and $`\tau `$ appear. In the regime where $`\tau T`$, the correlation time of the noise is much smaller than the typical variation time of the angle. Hence the noise can be considered to be white, the averaging procedure can be applied as in Section IV and the scalings found in (88) are correct. When $`T\tau `$ the noise becomes correlated over a period of the free system and can not be treated as white anymore. Now, $`T\tau `$ corresponds to $`E\tau ^{\frac{2n}{n1}}`$ which leads to crossover time $`t_c`$ of the order $`t_c^{\frac{n}{n1}}\tau ^{\frac{2n}{n1}}`$ i.e. $`t_c\tau ^2`$. For times larger than $`t_c`$, the scalings (93) are observed.

## VI Conclusion

A particle trapped in a confining potential with white multiplicative noise undergoes anomalous diffusion: if the confining potential grows as $`x^{2n}`$ at infinity, the particle diffuses as $`D_xt^{\beta _n}`$. We have calculated the anomalous diffusion exponent $`\beta _n=1/(2n2)`$, and the coefficient $`D_x`$. We have found similar laws for the diffusion of velocity and energy. Thanks to generalized equipartition identities, we have derived universal relations between the exponents and between the prefactors. Our calculations are based on the assumption that in the long-time limit the probability distribution function becomes uniform in the phase variable. By averaging out the phase variations, an effective, projected, dynamics for the action (or energy) can be defined. This technique enabled us to derive the asymptotic distribution law of the energy in the $`t\mathrm{}`$ limit, and to calculate its non-Gaussian features (skewness and flatness). Our analytical results agree with the numerical simulations within the numerical error bars. Thus, the averaging procedure produces very accurate results; it would be an interesting and challenging problem to characterize deviations from our results and to calculate subleading corrections.

In the case of colored multiplicative noise, we have deduced the anomalous diffusion exponents from an elementary scaling argument. Our result, supported by numerics, shows that the exponents are halved in presence of time correlations. The efficiency of parametric amplification decreases if the noise is coherent over a period of the system and therefore the particle diffuses at a much slower rate. In this case, however, the averaging technique is harder to apply because the noise itself is averaged out to the leading order. A precise calculation in the case of colored Gaussian noise still remains to be done.

We have considered only Hamiltonian systems, i.e. systems where no friction is present. Nevertheless, if the damping is small, the results we have derived for the undamped oscillator remain valid until the crossover time (identical to the typical decay time of the energy) is reached. The general case of a non-linear oscillator with (linear) friction leads to interesting results and is currently under study phkir .

###### Acknowledgements.

It is a pleasure to thank Yves Pomeau for encouraging us to work in this field of nonlinear stochastic equations and for his advice. K.M. is grateful to Michel Bauer for many useful discussions.

## Appendix A Nonlinear oscillator with external and internal noise

In this Appendix we discuss the behavior of an oscillator submitted to both additive and multiplicative noises. Because we are considering non-dissipative systems, there is no stationary probability distribution; the position, velocity and energy of the system satisfy scaling laws.

We first consider the case where the noise is only additive. The dynamics of a linear oscillator with additive noise is exactly solvable and the solution of the equation:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\omega ^2x(t)=\xi (t),$$
(94)
is given by:

$$x(t)=\frac{1}{\omega }_0^t\mathrm{sin}\left(\omega (ts)\right)\xi (s)\text{d}s.$$
(95)
We conclude from this expression that $`x(t)`$ grows diffusively with time:

$$x^2(t)=\frac{𝒟}{2\omega ^2}\left(t\frac{\mathrm{sin}2\omega t}{2\omega }\right).$$
(96)
Such a diffusive behavior is entirely different from the exponential growth of a linear oscillator subject to parametric noise.

When the nonlinear restoring force is cubic, the linear term can be neglected in the long-time limit and the dynamics is given by

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)^3=\xi (t),$$
(97)
where $`\xi (t)`$ is the Gaussian white noise defined in Eq. (2). This equation can be analyzed as in Sec. III.1. After introducing the energy and angle variables, defined in Eqs. (10) and (19), we obtain

$`\dot{E}`$ $`=`$ $`\dot{x}\xi (t)=(2E)^{\frac{1}{2}}{\displaystyle \frac{\mathrm{cn}(\varphi ;1/\sqrt{2})}{\mathrm{dn}^2(\varphi ;1/\sqrt{2})}}\xi (t);`$ (98)
$`\dot{\varphi }`$ $`=`$ $`(4E)^{\frac{1}{4}}{\displaystyle \frac{\mathrm{sd}(\varphi ;1/\sqrt{2})}{2\sqrt{2}E^{1/2}}}\xi (t).`$ (99)
We can show that the auxiliary variable $`\mathrm{\Theta }`$ defined as

$$\mathrm{\Theta }=E^{\frac{1}{2}},$$
(100)
undergoes a normal diffusion process that can be precisely studied after averaging out the fast variable $`\varphi `$. The P.D.F. of the variable $`\mathrm{\Theta }`$, the scaling laws for $`E`$, $`\dot{x}`$ and $`x`$, as well as the generalized diffusion constants can be calculated exactly, as before. We indicate here only the asymptotic power-laws:

$`E`$ $``$ $`t,`$
$`x`$ $``$ $`t^{\frac{1}{4}},`$
$`\dot{x}`$ $``$ $`t^{\frac{1}{2}}.`$ (101)
The scaling exponents for additive noise are different from those for multiplicative noise. For exemple, we observe that a particle submitted to an additive noise in a quartic potential is subdiffusive with an anomalous exponent equal to $`1/2`$, whereas in the presence of multiplicative noise it behaves diffusively.

In the general case, with a nonlinearity of the type $`x^{2n1}`$, with $`n1`$, one can prove that the variable $`\mathrm{\Theta }`$, defined in Eq. (100), diffuses as $`t^{1/2}`$ in the long time limit. Thus the following scalings are satisfied when $`t\mathrm{}`$:

$`E`$ $``$ $`t,`$
$`x`$ $``$ $`t^{\frac{1}{2n}},`$
$`\dot{x}`$ $``$ $`t^{\frac{1}{2}}.`$ (102)
All these exponents are smaller than the exponents found for a multiplicative white noise in Eq. (88).

Finally, when both additive and multiplicative noises are present, the oscillator is governed by the equation:

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)\xi _{\text{mult}}(t)+x(t)^{2n1}=\xi _{\text{add}}(t),$$
(103)
where $`\xi _{\text{mult}}`$ and $`\xi _{\text{add}}`$ are independent white noises of amplitude $`𝒟_{\text{mult}}`$ and $`𝒟_{\text{add}},`$ respectively. If we study the energy variation due to noise, $`\dot{E}x\dot{x}\xi _{\text{mult}}+\dot{x}\xi _{\text{add}}E^{\frac{n+1}{2n}}\xi _{\text{mult}}+E^{1/2}\xi _{\text{add}}`$, we observe that the first term has a dominant effect. From this simple argument, we conclude that the multiplicative noise is expected to dominate over the additive noise and therefore, asymptotically, the scaling laws will be those derived for the multiplicative noise alone. However, a crossover between the two scalings (88) and (102) should be observed by choosing $`𝒟_{mult}𝒟_{add}`$. Comparing Eqs. (88) and (102), we find that the effect of the multiplicative noise starts to dominate after a crossover time of the order of $`(𝒟_{mult}t_c)^{\frac{1}{2n2}}(𝒟_{add}t_c)^{\frac{1}{2n}}`$ i.e., $`t_c𝒟_{add}^{n1}/𝒟_{mult}^n`$.

## Appendix B Numerical Algorithm

The algorithm used to integrate numerically the stochastic ordinary differential equations studied in this article is the one-step collocation method advocated in mannella . In this Appendix, we recall the general principles underlying this method, and give the algorithms we used to integrate Eqs. (57) and (89)-(91), for white and correlated noise respectively. All stochastic equations are understood according to Stratonovich rules.

Let $`\{x_i\}_{i=1,\mathrm{},N}`$ be $`N`$ real variables of time $`t`$, and $`\xi (t)`$ a stochastic process assumed to be Gaussian and white. We wish to solve systems of $`N`$ coupled Langevin equations of the form

$$\dot{x}_i=f_i(\{x_j(t)\})+g_i(\{x_j(t)\})\xi (t),$$
(104)
where $`f_i`$ and $`g_i`$ are $`N`$ (smooth) functions of the $`x_i`$’s. Let $`\delta t`$ be the integration timestep. Upon integrating formally Eq. (104) between $`0`$ and $`\delta t`$, we obtain the following set of coupled equations implicit in $`\{x_i(t)\}`$

$$x_i(\delta t)x_i(0)=_0^{\delta t}f_i(\{x_j(s)\})ds+_0^{\delta t}g_i(\{x_j(s)\})\xi (s)ds.$$
(105)
For $`\delta t`$ small enough, the functions $`f_i`$ and $`g_i`$ may be Taylor-expanded in the vicinity of $`t=0`$, *e.g.*,

$$f_i(\{x_j(s)\})=f_i(\{x_j(0)\})+_kf_i(\{x_j(0)\})\delta x_k(s)+\frac{1}{2}\left(_k_lf_i(\{x_j(0)\})\right)\delta x_k(s)\delta x_l(s)+\mathrm{},$$
(106)
where $`\delta x_k(s)=x_k(s)x_k(0)`$. Replacing the functions $`f_i`$ and $`g_i`$ in Eq. (105) by the expansion (106), we obtain a new set of coupled integral equations implicit in $`\{\delta x_i\}`$. Upon solving these equations up to an arbitrary order in $`\delta t`$, we obtain $`\delta x_k(\delta t)`$. In practice, we implemented an algorithm exact to $`O(\delta t^2)`$.

### B.1 White noise

We integrate the following set of first-order differential equations

$`\dot{x}`$ $`=`$ $`v,`$ (107)
$`\dot{v}`$ $`=`$ $`\omega ^2x\lambda x^{2n1}+x\xi ,`$ (108)
where the Gaussian white noise $`\xi (t)`$ verifies Eq. (2). The exact evolution equations for $`\delta x(\delta t)=x(\delta t)x(0)`$ and $`\delta v(\delta t)=v(\delta t)v(0)`$ read

$`\delta x(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}ds\left(v(0)+\delta v(s)\right),`$ (109)
$`\delta v(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}ds\left\{\omega ^2\left(x(0)+\delta x(s)\right)\lambda \left(x(0)+\delta x(s)\right)^{2n1}+\left(x(0)+\delta x(s)\right)\xi (s)\right\}.`$ (110)
The auxiliary variables $`Z_1`$ and $`Z_2`$, defined by

$`Z_1(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}\xi (s)ds,`$ (111)
$`Z_2(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}Z_1(s)ds,`$ (112)
are Gaussian random variables, with zero average, and the following correlations: $`Z_1(\delta t)^2=𝒟\delta t`$, $`Z_2(\delta t)^2=𝒟\delta t^3/3`$, $`Z_1(\delta t)Z_2(\delta t)=𝒟\delta t^2/2`$. Up to order $`(\delta t)^2`$, we find

$`x(\delta t)`$ $`=`$ $`x(0)+v(0)\delta t+Z_2(\delta t)+{\displaystyle \frac{1}{2}}\delta t^2\left(\omega ^2x(0)\lambda x(0)^{2n1}\right),`$ (113)
$`v(\delta t)`$ $`=`$ $`v(0)+Z_1(\delta t)+\delta t\left(\omega ^2x(0)\lambda x(0)^{2n1}\right)+{\displaystyle \frac{1}{2}}\delta t^2v(0)\left(\omega ^2x(0)(2n1)\lambda x(0)^{2n2}\right).`$ (114)
In practice, we use two independent Gaussian random noises $`Z_1`$ and $`Y_1`$ with zero mean and correlations $`Z_1(\delta t)^2=𝒟\delta t`$, $`Y_1(\delta t)^2=𝒟\delta t`$. As shown in Ref. mannella , the variable $`Z_2`$ may be approximated by the expression:

$$Z_2(\delta t)=\delta t\left(\frac{1}{2}Z_1(\delta t)+\frac{1}{2\sqrt{3}}Y_1(\delta t)\right),$$
(115)
when the algorithm is exact up to order $`\delta t^2`$.

### B.2 Colored noise

When the noise $`\eta (t)`$ is correlated, we must solve a set of three coupled equations

$`\dot{x}`$ $`=`$ $`v,`$ (116)
$`\dot{v}`$ $`=`$ $`\omega ^2x\lambda x^{2n1}+x\eta ,`$ (117)
$`\dot{\eta }`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}\eta +{\displaystyle \frac{1}{\tau }}\xi .`$ (118)
The set of exact integral equations becomes

$`\delta x(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}ds\left(v(0)+\delta v(s)\right),`$ (119)
$`\delta v(\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}ds\left\{\omega ^2\left(x(0)+\delta x(s)\right)\lambda \left(x(0)+\delta x(s)\right)^{2n1}+\left(x(0)+\delta x(s)\right)\left(\eta (0)+\delta \eta (s)\right)\right\},`$ (120)
$`\delta \eta (\delta t)`$ $`=`$ $`{\displaystyle _0^{\delta t}}ds\left\{{\displaystyle \frac{1}{\tau }}\left(\eta (0)+\delta \eta (s)\right)+{\displaystyle \frac{1}{\tau }}\xi (t)\right\}.`$ (121)
The algorithm to order $`\delta t^2`$ reads

$`x(\delta t)`$ $`=`$ $`x(0)+v(0)\delta t+{\displaystyle \frac{1}{2}}v(0)\delta t^2,`$ (122)
$`v(\delta t)`$ $`=`$ $`v(0)+\delta t\left(\omega ^2x(0)\lambda x(0)^{2n1}+x(0)\eta (0)\right)+{\displaystyle \frac{1}{\tau }}x(0)Z_2(\delta t)`$ (123)
$`+{\displaystyle \frac{1}{2}}\delta t^2\left(v(0)\left(\omega ^2x(0)(2n1)\lambda x(0)^{2n2}+\eta (0)\right){\displaystyle \frac{1}{\tau }}x(0)\eta (0)\right),`$
$`\eta (\delta t)`$ $`=`$ $`\eta (0)+{\displaystyle \frac{1}{\tau }}Z_1(\delta t)+{\displaystyle \frac{1}{\tau }}\eta (0)\delta t{\displaystyle \frac{1}{\tau ^2}}Z_2(\delta t)+{\displaystyle \frac{1}{2}}\delta t^2{\displaystyle \frac{\eta (0)}{\tau ^2}}.`$ (124)
## Appendix C Some useful Properties of Elliptic Functions

The Jacobian elliptic functions $`\mathrm{sn}(u;k),\mathrm{cn}(u;k),`$ and $`\mathrm{dn}(u;k)`$ are meromorphic functions of the complex variable $`u`$. They depend on a parameter $`k`$ called the elliptic modulus abram ; byrd . These functions, like the trigonometric functions, have a real period and, like the hyperbolic functions, have an imaginary period. They are thus doubly periodic in the complex plane and their periods are $`(4K,2iK^{})`$, $`(4K,2K+2iK^{})`$ and $`(2K,4iK^{})`$, respectively, where

$$K(k)=_0^1\frac{\mathrm{d}t}{\sqrt{(1t^2)(1k^2t^2)}}=_0^{\frac{1}{k^{}}}\frac{\mathrm{d}y}{\sqrt{(1k_{}^{}{}_{}{}^{2}y^2)(1+k^2y^2)}}\text{ with }t^2=\frac{y^2}{1+k^2y^2}\text{ and }k_{}^{}{}_{}{}^{2}=1k^2.$$
(125)
The imaginary period is given by $`K^{}(k)=K(k^{})=K\left(\sqrt{1k^2}\right).`$ The Jacobian elliptic functions satisfy the fundamental relations:

$`\mathrm{sn}^2(\varphi ;k)+\mathrm{cn}^2(\varphi ;k)`$ $`=`$ $`1,`$
$`k^2\mathrm{sn}^2(\varphi ;k)+\mathrm{dn}^2(\varphi ;k)`$ $`=`$ $`1,`$ (126)
and their derivatives are

$`{\displaystyle \frac{\text{d}}{\text{d}\varphi }}\mathrm{sn}(\varphi ;k)`$ $`=`$ $`\mathrm{cn}(\varphi ;k)\mathrm{dn}(\varphi ;k),`$
$`{\displaystyle \frac{\text{d}}{\text{d}\varphi }}\mathrm{cn}(\varphi ;k)`$ $`=`$ $`\mathrm{sn}(\varphi ;k)\mathrm{dn}(\varphi ;k),`$
$`{\displaystyle \frac{\text{d}}{\text{d}\varphi }}\mathrm{dn}(\varphi ;k)`$ $`=`$ $`k^2\mathrm{sn}(\varphi ;k)\mathrm{cn}(\varphi ;k).`$ (127)
If we define the function $`\mathrm{sd}(\varphi ;k)=\mathrm{sn}(\varphi ;k)/\mathrm{dn}(\varphi ;k),`$ the relations (126) and (127) lead to two new formulas

$$\frac{\text{d}}{\text{d}\varphi }\mathrm{sd}(\varphi ;k)=\frac{\mathrm{cn}(\varphi ;k)}{\mathrm{dn}^2(\varphi ;k)}\text{ and }\frac{\mathrm{cn}^2(\varphi ;k)}{\mathrm{dn}^4(\varphi ;k)}=(1k_{}^{}{}_{}{}^{2}\mathrm{sd}^2(\varphi ;k))(1+k^2\mathrm{sd}^2(\varphi ;k)).$$
(128)
Moreover, from Eqs. (127) and (128), we obtain

$`\mathrm{sn}^1(y;k)`$ $`=`$ $`{\displaystyle _0^y}{\displaystyle \frac{\text{d}t}{\sqrt{(1t^2)(1k^2t^2)}}},`$ (129)
$`\mathrm{sd}^1(y;k)`$ $`=`$ $`{\displaystyle _0^y}{\displaystyle \frac{\text{d}u}{\sqrt{(1k_{}^{}{}_{}{}^{2}u^2)(1+k^2u^2)}}}.`$ (130)
Similar relations can be derived for $`\mathrm{cn}^1`$ and $`\mathrm{dn}^1.`$ The expressions of inverse elliptic functions in terms of quartic integrals provide a classic way of defining elliptic functions. The change of variables $`t=u/\sqrt{1+k^2u^2}`$ in Eq. (129) leads to the identity

$$\mathrm{sd}^1(y;k)=\mathrm{sn}^1(\frac{y}{\sqrt{1+k^2y^2}};k)$$
(131)
We deduce from Eq. (129)

$$\frac{\text{d}}{\text{d}k}\mathrm{sn}^1(y;k)=k_0^y\frac{t^2\text{d}t}{\sqrt{(1t^2)(1k^2t^2)}}=k_0^{y/\sqrt{1k^2y^2}}\frac{u^2\text{d}u}{\sqrt{(1k_{}^{}{}_{}{}^{2}u^2)(1+k^2u^2)}}=k_0^{\varphi _1}\mathrm{sd}^2(\varphi ;k)\text{d}\varphi ,$$
(132)
where the last equality was obtained by setting $`u=\mathrm{sd}(\varphi ;k)`$ with $`\varphi _1=\mathrm{sn}^1(y;k).`$

We verify from Eqs. (126) and (127) that $`x`$ and $`\dot{x}`$ defined in Eqs. (11) and (12) are solutions of the cubic oscillator $`\ddot{x}+x^3=0.`$ Similarly, we verify that $`x`$ and $`\dot{x}`$ defined in Eqs. (48) and (49) satisfy the differential equation for the Duffing oscillator, $`\ddot{x}+\omega ^2x+x^3=0.`$ We now prove the identities used to derive the equipartition relations and the diffusion constants in section III.1. In the case of a purely cubic oscillator, $`k=k^{}=1/\sqrt{2},`$ and we obtain from Eq. (130)

$$\frac{\text{d}}{\text{d}y}\mathrm{sd}^1(y;\frac{1}{\sqrt{2}})=\frac{1}{\sqrt{1\frac{y^2}{4}}}.$$
(133)
From Eq. (128) we deduce

$$\overline{\frac{\mathrm{cn}^2(\varphi )}{\mathrm{dn}^4(\varphi )}}=\frac{1}{K}_0^Kd\varphi \frac{\mathrm{cn}^2(\varphi )}{\mathrm{dn}^4(\varphi )}=\frac{_0^1\sqrt{1u^4}\text{d}u}{_0^1\frac{\text{d}u}{\sqrt{1u^4}}}=\frac{2}{3},$$
(134)
where we have set $`u=\mathrm{sd}(\varphi )/\sqrt{2}`$ to obtain the second equality; the final equality results from Eq. (72). Similarly, we obtain:

$$\overline{\frac{\mathrm{cn}^4(\varphi )}{\mathrm{dn}^8(\varphi )}}=\frac{1}{K}_0^Kd\varphi \frac{\mathrm{cn}^4(\varphi )}{\mathrm{dn}^8(\varphi )}=\frac{_0^1(1u^4)^{\frac{3}{2}}\text{d}u}{_0^1\frac{\text{d}u}{\sqrt{1u^4}}}=\frac{4}{7}.$$
(135)
The last equality is obtained by the following integration by parts:

$$_0^11.(1u^4)^{\frac{3}{2}}\text{d}u=_0^16u^4(1u^4)^{\frac{1}{2}}\text{d}u=_0^16(1u^4)^{\frac{1}{2}}\text{d}u_0^16(1u^4)^{\frac{3}{2}}\text{d}u;$$

thus

$$_0^1(1u^4)^{\frac{3}{2}}\text{d}u=\frac{6}{7}_0^1(1u^4)^{\frac{1}{2}}\text{d}u.$$

Now, using Eq. (134), we obtain the numerical value 4/7 in Eq. (135). We also have:

$$\overline{\mathrm{sd}^2(\varphi )}=\frac{1}{K}_0^Kd\varphi \mathrm{sd}^2(\varphi )=2\frac{_0^1\frac{u^2\text{d}u}{\sqrt{1u^4}}}{_0^1\frac{\text{d}u}{\sqrt{1u^4}}}=\frac{5}{2},$$
(136)
where $``$ was defined in Eq. (33). The last equality is obtained by using the following identities

$$_0^1\frac{u^6\text{d}u}{\sqrt{1u^4}}=_0^1\frac{u^3}{2}\frac{2u^3\text{d}u}{\sqrt{1u^4}}=\frac{3}{2}_0^1u^2\sqrt{1u^4}\text{d}u,$$

$$_0^1\frac{u^2\text{d}u}{\sqrt{1u^4}}\text{d}u=_0^1\left(u^2\sqrt{1u^4}+\frac{u^6}{\sqrt{1u^4}}\right)\text{d}u=\frac{5}{2}_0^1u^2\sqrt{1u^4}\text{d}u.$$

We end this Appendix by explaining how the equations for the phase dynamics are obtained. Eq. (20) is derived from Eqs. (133) and (18):

$$\frac{\text{d}}{\text{d}t}\varphi =\frac{\text{d}}{\text{d}t}\mathrm{sd}^1(\frac{x}{E^{1/4}},\frac{1}{\sqrt{2}})=\frac{\frac{\dot{x}}{E^{1/4}}\frac{x\dot{E}}{4E^{5/4}}}{\sqrt{1\frac{x^4}{4E}}}=\left(\frac{\dot{x}}{E^{1/4}}\frac{x^2\dot{x}\xi }{4E^{5/4}}\right)\frac{\sqrt{2E}}{\dot{x}}=(4E)^{\frac{1}{4}}\frac{\mathrm{sd}^2(\varphi ;1/\sqrt{2})}{2\sqrt{2}E^{1/4}}\xi (t).$$
(137)
In order to derive the dynamics of the phase of a Duffing oscillator (53), we must differentiate Eq. (51) with respect to time. We must bear in mind that the elliptic modulus $`k`$ defined in Eq. (50) is now a function of time. Therefore, when differentiating the function $`\mathrm{sd}^1`$, we obtain two contributions: $`(i)`$ the derivative with respect to the argument produces a result similar to the one for the cubic case, $`(ii)`$ the derivative with respect to the elliptic modulus, calculated using Eqs. (131) and (132), generates the two last terms in Eq. (53).


format: text

text:
# Renormalization of the Minimal Supersymmetric Standard Model

## Abstract

The renormalization of the Minimal Supersymmetric Standard Model (MSSM) is presented. We describe symmetry identities that constitute a framework in which the MSSM is completely characterized and renormalizability can be proven. Furthermore, we discuss applications of this framework for the determination of symmetry-restoring counterterms, the gauge dependence of $`\mathrm{tan}\beta `$ and the derivation of non-renormalization theorems.

In this talk the renormalization of the Minimal Supersymmetric Standard Model (MSSM) is presented . A framework is set up where all counterterms are uniquely determined. This comprises a set of symmetry identities providing a complete characterization of the MSSM and a set of on-shell renormalization conditions that forbid on-shell mixing between different physical fields. In this framework it has been shown that the MSSM is multiplicatively renormalizable, infrared finite, and that all on-shell conditions can be satisfied simultaneously.

This study is motivated by the fact that no satisfactory regularization for supersymmetric gauge theories is known. In particular, dimensional regularization breaks supersymmetry; hence, supersymmetry-restoring counterterms have to be calculated and added. On the other hand, dimensional reduction is mathematically inconsistent, and therefore its area of validity is unclear. In practice, renormalization of the MSSM means first to check whether the chosen regularization preserves the symmetries and to add — if necessary — symmetry-restoring counterterms, and second to add the usual symmetric counterterms (corresponding to field and parameter renormalization) in order to cancel divergences and satisfy renormalization conditions.

This raises the deeper question of how to formulate the symmetries of the MSSM at all on the quantum level and on the level of Green functions. For the Standard Model  and general supersymmetric Yang–Mills theories  the answer is known: the Slavnov–Taylor and Ward identities provide a complete characterization of the symmetries. Similar identities should also be formulated in the MSSM.

The outline of the talk is as follows. First the symmetry identities of the MSSM are presented. Then we draw important conclusions of these identities, in particular on the proof of renormalizability and the practical determination of symmetry-restoring counterterms. Finally, two applications of the symmetry identities are discussed, concerning the gauge dependence of the parameter $`\mathrm{tan}\beta `$ and a new approach to the non-renormalization theorems in supersymmetric gauge theories.

The basic symmetries of the (electroweak part of the) MSSM are spontaneously broken SU(2)$`\times `$U(1) gauge invariance and softly broken supersymmetry. Clearly, the basic structure of the symmetry identities in the MSSM can be obtained by combining the results for the Standard Model  and for general supersymmetric Yang–Mills theories . The main symmetry content is described by a Slavnov–Taylor identity

$$S(\mathrm{\Gamma })=0,$$
(1)
where $`\mathrm{\Gamma }`$ denotes the generating functional of the one-particle irreducible Green functions. It combines all information on gauge invariance and supersymmetry including the quantum corrections to the transformations and the commutation relations of the generators. However, the Slavnov–Taylor identity does not fix the values of the hypercharges $`Y_i`$, which however are crucial in order to fix the electric charges $`Q_i=T_i^3+Y_i/2`$ correctly. The $`Y_i`$ are fixed by a local Ward identity for the U(1)-symmetry:

$`^\mu {\displaystyle \frac{\delta \mathrm{\Gamma }}{\delta V^{}^\mu }}=ig^{}{\displaystyle \underset{\mathrm{Fields}\phi _i}{}}{\displaystyle \frac{Y_i}{2}}\phi _i{\displaystyle \frac{\delta \mathrm{\Gamma }}{\phi _i}}+\begin{array}{c}\text{gauge-fixing}\\ \text{terms}\end{array}`$ (4)
Along with this local Ward identity, global Ward identities describing global SU(2)$`\times `$U(1) and $`R`$ invariance are formulated. The symmetry breakings are introduced by using external fields with a constant shift. For soft supersymmetry breaking a chiral supermultiplet $`(A,a_\alpha ,F_A+v_A)`$ with a shift in its highest component is used, and for spontaneous breaking of gauge invariance an SU(2)$`\times `$U(1)-multiplet $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ is used. Using these fields, the symmetry identities take the same form as in the cases with unbroken symmetries, but when the external fields are set to their constant values, symmetry breaking is described.

While the basic structure of the symmetry identities seems obvious, the difficulty lies in the detailed implementation. In fact, it turns out that the $`R_\xi `$-gauge requires that in the MSSM the detailed structure of the external fields $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ appearing in the symmetry identities must differ from the one in the literature.

The problem is that the MSSM contains an extended Higgs sector, and even if CP-conservation is assumed the physical CP-odd Higgs boson $`A^0`$ and the charged Higgs bosons $`H^\pm `$ can mix with unphysical degrees of freedom:

$`A^0`$ $``$ $`(G^0,A_{\mathrm{Long}}^\mu ,Z_{\mathrm{Long}}^\mu ),`$
$`H^\pm `$ $``$ $`(G^\pm ,W^\pm {}_{\mathrm{Long}}{}^{\mu }).`$ (5)
However, in the $`R_\xi `$-gauge fixing terms only the unphysical fields appear in the gauge-fixing functions (we restrict ourselves to the case of the neutral fields for simplicity):

$$F^A=^\mu A_\mu ,F^Z=^\mu Z_\mu +\xi M_ZG^0.$$
(6)
In this form, the gauge fixing would even break global SU(2)$`\times `$U(1)-invariance and the U(1)-Ward identity. As proposed in , the U(1)-Ward identity can be restored by using the external field $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ for writing the gauge-fixing functions. However, if $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ is chosen as an SU(2)-doublet like in the Standard Model-case, it turns out that necessarily $`A^0`$ and/or $`H^\pm `$ appear in the gauge fixing. Thus we are lead to the question which multiplet assignment to choose for $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$. The answer has been found in . The multiplet structure of $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ is chosen as the product of the adjoint and doublet representation of SU(2)$`\times `$U(1), and there are two of these 8-component multiplets, one for each Higgs doublet $`H_{1,2}`$. Then the gauge-fixing functions can be written as

$$F^a=^\mu V_\mu ^a2\mathrm{I}\mathrm{m}((\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})_i^a{}_{}{}^{}H_{i}^{}),$$
(7)
so that $`F^a`$ transforms in the adjoint representation and is compatible with the U(1)-Ward identity. At the same time, $`\widehat{\mathrm{v}}`$ has enough components that can be adjusted such that the $`R_\xi `$-gauge conditions are reproduced for $`\widehat{\mathrm{\Phi }}=0`$, i.e. $`F^a|_{\widehat{\mathrm{\Phi }}0}`$ does not contain the fields $`A^0`$, $`H^\pm `$, and the $`F^{A,Z}`$-components coincide with (6).

A more complicated multiplet structure of this external field, however, also complicates the proof of infrared finiteness of the MSSM. Consider the counterterm $`_{\mathrm{ct}}=\widehat{\mathrm{\Phi }}_jA^\mu A_\mu `$ that might be necessary to restore symmetries as an example. In the calculation of global Ward identities there appear terms like $`W\mathrm{\Gamma }=\mathrm{}+d^4x\widehat{\mathrm{v}}_i\frac{\delta \mathrm{\Gamma }}{\delta \widehat{\mathrm{\Phi }}_j}`$, leading to diagrams of the type shown in Fig. 1. Since the $`\widehat{\mathrm{\Phi }}_j`$-field carries no momentum, all diagrams of this type contain the infrared divergent integral $`\frac{d^4k}{k^4}`$, no matter what other external lines and momenta are present. It has to be shown that such situations cannot arise and that $`\mathrm{\Gamma }`$ itself as well as all Ward and Slavnov–Taylor identities are infrared finite.

The proof of infrared finiteness can be carried out by identifying the dangerous components of $`\widehat{\mathrm{\Phi }}`$ and also of $`(A,a_\alpha ,F_A+v_A)`$ using infrared powercounting, finding an optimal assignment of the infrared quantum numbers and checking that every dangerous term is in fact excluded for a symmetry reason.

Combining the mentioned elements yields the detailed structure of the symmetry identities in the MSSM. The main identities are the Slavnov–Taylor and Ward identities in presence of the external $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$ and $`(A,a_\alpha ,F_A+v_A)`$ fields, whose multiplet structure and infrared powercounting has been determined. There are further identities corresponding to gauge-fixing conditions, but their discussion is beyond the scope of the present talk (see sec. 3.2 of ). The set of symmetry identities constitutes a full definition of the MSSM at the quantum level.

Using the symmetry identities we can draw the following conclusions of basic importance: the symmetry-restoring counterterms are uniquely determined, the remaining freedom consists of adding symmetric counterterms — which can be shown to correspond to multiplicative renormalization, and infrared finiteness can be completely proven.

Furthermore, the Slavnov–Taylor identity can be used to show that all on-shell renormalization conditions can be satisfied simultaneously . Since the MSSM is multiplicatively renormalizable, it is quite obvious that most on-shell conditions can be satisfied be choosing the field renormalization constants appropriately. An exception arises from the mixings (5). The question is whether it is possible to establish on-shell renormalization conditions that characterize the fields $`A^0`$, $`H^\pm `$ as mass eigenstates. Using field renormalization it is always possible to satisfy the on-shell conditions (neglecting the finite widths of the particles)

$`\mathrm{\Gamma }_{A^0A^0}(M_A^2)=\mathrm{\Gamma }_{A^0G^0}(M_A^2)=0,`$
$`\mathrm{\Gamma }_{H^+H^{}}(M_{H^\pm }^2)=\mathrm{\Gamma }_{H^+G^{}}(M_{H^\pm }^2)=0.`$ (8)
The Slavnov–Taylor identity yields

$`0`$ $`=`$ $`{\displaystyle \underset{\phi =A^0,G^0}{}}\mathrm{\Gamma }_{c^aY_\phi }\mathrm{\Gamma }_{A^0\phi }\mathrm{\Gamma }_{c^a\overline{c}^b}{\displaystyle \frac{1}{\xi }}{\displaystyle \frac{\delta F^b}{\delta A^0}}`$ (9)
$`+{\displaystyle \underset{V=A,Z}{}}\mathrm{\Gamma }_{c^aY_{V_\mu }}\mathrm{\Gamma }_{A^0V_\mu },`$
and a similar identity for $`A^0H^\pm `$, where $`c^a,\overline{c}^a`$ are the Faddeev–Popov (anti-)ghosts and $`Y_\phi `$ denotes the source of the BRS transformation of $`\phi `$. Because of (8) and because the $`R_\xi `$-gauge can be realized, the first two terms of (9) vanish for on-shell momenta, and we obtain

$`\mathrm{\Gamma }_{A^0A^\mu }(M_A^2)=\mathrm{\Gamma }_{A^0Z^\mu }(M_A^2)=0,`$
$`\mathrm{\Gamma }_{H^+W^{}^\mu }(M_{H^\pm }^2)=0.`$ (10)
Hence, indeed no on-shell mixing between the physical and unphysical degrees of freedom occurs and all on-shell conditions can be satisfied.

This completes the general discussion of the renormalization and provides the basis for practical applications.

As the first application we will briefly discuss how supersymmetry-restoring counterterms can be determined in practice. Generally, the counterterms have to be chosen such that the symmetry identities hold after renormalization. For supersymmetry, to kinds of identities are important. Using111There are fields $`\phi `$ with no corresponding source $`Y_\phi `$. For them, $`\frac{\delta \mathrm{\Gamma }}{\delta ϵ\delta Y_{\phi _i}}`$ has to be replaced by $`_ϵs\phi _i`$, where $`s\phi _i`$ is the BRS variation of $`\phi _i`$.

$$0=\frac{\delta S(\mathrm{\Gamma })}{\delta ϵ}=\underset{\mathrm{Fields}\phi _i}{}\frac{\delta \mathrm{\Gamma }}{\delta ϵ\delta Y_{\phi _i}}\frac{\delta \mathrm{\Gamma }}{\delta \phi _i},$$
(11)
where $`ϵ`$ is the ghost corresponding to supersymmetry, we obtain identities corresponding to supersymmetry relations, where the prefactors are quantum corrected and renormalized supersymmetry transformations. Using

$`0`$ $`=`$ $`{\displaystyle \frac{\delta S(\mathrm{\Gamma })}{\delta ϵ\delta \overline{ϵ}\delta Y_{\phi _j}}}`$ (12)
$`=`$ $`{\displaystyle \underset{\mathrm{Fields}\phi _i}{}}({\displaystyle \frac{\delta \mathrm{\Gamma }}{\delta \overline{ϵ}\delta Y_{\phi _i}}}{\displaystyle \frac{\delta \mathrm{\Gamma }}{\delta ϵ\delta Y_{\phi _j}\phi _i}}+(ϵ\overline{ϵ}))`$
$`\pm 2i\sigma ^\mu _\mu \phi _j+\mathrm{},`$
identities corresponding to the supersymmetry algebra are obtained. Here $`\pm `$ holds for bosonic/fermionic $`\phi _j`$, respectively, and the dots denote calculable terms corresponding to gauge transformations and equations of motion in the supersymmetry algebra. By the requirement that these identities should be satisfied after renormalization, it is possible to determine first the counterterms for the renormalized supersymmetry transformations, and then the counterterms to the vertex functions without BRS insertions. The results are unique up to symmetric counterterms, which correspond to multiplicative renormalization. Identities of the first kind have already been considered in , but they alone do not lead to unique results for symmetry-restoring counterterms and cannot serve as tests of the symmetry of a given regularization scheme. Both kinds of identities have been considered in  for supersymmetric QED and QCD at the one-loop level. In all cases it has been found that dimensional reduction preserves the identities at the regularized level. Thus this scheme is indeed supersymmetric in the cases considered up to now.

In the remainder of this talk we want to discuss two applications of the Slavnov–Taylor identities of the MSSM or other supersymmetric models: calculating the gauge dependence of the parameter $`\mathrm{tan}\beta `$  and deriving non-renormalization theorems . Both use certain extensions of the Slavnov–Taylor identities as tools. The gauge dependence is calculated using a Slavnov–Taylor identity containing an additional BRS transformation of the gauge parameter, and the non-renormalization theorems are derived by introducing BRS transformations of the coupling constants.

The quantity $`\mathrm{tan}\beta `$ is one of the main input parameters of the MSSM. At tree level,

$$\mathrm{tan}\beta =\frac{v_2}{v_1}.$$
(13)
In , different renormalization schemes for $`\mathrm{tan}\beta `$ were analyzed with the aim to find a scheme that defines $`\mathrm{tan}\beta `$ at the same time in a gauge-independent and process-independent way. Using the extended Slavnov–Taylor identity, the gauge dependence of $`\mathrm{tan}\beta `$ can be explicitly calculated for any given scheme. It was found that in several well-known schemes, the $`\overline{DR}`$-scheme and the schemes introduced in , $`\mathrm{tan}\beta `$ is defined in a gauge-dependent way, i.e. the relation between $`\mathrm{tan}\beta `$ and observable quantities that can be used to extract $`\mathrm{tan}\beta `$ from experiment is gauge dependent. The only exception is the $`\overline{DR}`$-scheme, which is gauge independent if its application is restricted to $`R_\xi `$-gauges and to the one-loop level.

Hence, a large class of new schemes was considered where $`\mathrm{tan}\beta `$ is defined in the Higgs sector, as at the tree level. All gauge-independent schemes in this class were identified. Unfortunately it turned out that each of these schemes leads to severe numerical instabilities in the perturbative expansion, so these schemes are not useful in practice. Given these results, the $`\overline{DR}`$-scheme seems to be the best choice of all process-independent definitions for $`\mathrm{tan}\beta `$.

Non-renormalization theorems are among the deepest and most exciting properties of supersymmetric models. They state the absence of certain divergences, e.g. of quadratic divergences, and thus provide in particular a solution to the naturalness problem. In , a new approach towards these theorems has been developed, which has been applied to supersymmetric QED and QCD in . The origin of the non-renormalization theorems is identified as follows. Every term in a supersymmetric Lagrangian is the highest component of a supermultiplet and thus related to lower-dimensional field polynomials. Similarly, supersymmetry relates diagrams to diagrams with a lower degree of divergence.

This fact can be implemented into an extended Slavnov–Taylor identity by replacing the coupling constants by full supermultiplets. In this way, every supersymmetric term in the Lagrangian is replaced by a sum of the form (considering the case of chiral multiplets for simplicity)

$$g_{\mathrm{susy}}g_{\mathrm{susy}}\chi \mathrm{\Xi }+fA,$$
(14)
where $`(g,\chi ,f)`$ is the multiplet of the coupling and $`(A,\mathrm{\Xi },_{\mathrm{susy}})`$ the multiplet of the respective term in the Lagrangian. The higher components of the coupling thus couple to the lower components of the Lagrangian. The supermultiplet structure of the couplings implies BRS transformations of the couplings. If these are included into the Slavnov–Taylor identity, an extended identity is obtained that can be used to derive non-renormalization theorems.

The advantages of this approach are that the non-renormalization theorems can be derived without assuming a supersymmetric regularization, in the context of the Wess–Zumino gauge, and that it makes the underlying algebraic origin apparent. As an illustration, we list the results for supersymmetric QED. In this model, there are two independent divergences, conveniently expressed in terms of divergent renormalization constants for the charge and electron mass:

$`\delta Z_e^{(1)}:\text{ only one-loop},\delta Z_m^{(l)}.`$ (15)
The renormalization constants for all soft-breaking parameters can be expressed in terms of $`\delta Z_e`$ and $`\delta Z_m`$:

$`\delta Z_{M_\lambda }^{(1)}`$ $`=`$ $`2\delta Z_e^{(1)}:\text{ only one-loop},`$ (16)
$`\delta Z_b^{(l)}`$ $`=`$ $`(2lM_\lambda m/b+1)\delta Z_m^{(l)},`$ (17)
$`\delta Z_M^{(l)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}l(l+1)(M_\lambda ^2/b)\delta Z_m^{(l)}.`$ (18)
In addition to relating all renormalization constants, these equations imply that the charge and the photino mass counterterms are finite from the two-loop level on and that the two scalar mass counterterms are only logarithmically divergent.

Similar results can also be derived for non-abelian supersymmetric gauge theories. In this case a deep connection between the form of the non-renormalization theorems and two anomalies — the Adler–Bardeen anomaly and a supersymmetry anomaly in presence of the supercoupling  — is exhibited. As a byproduct, also the non-renormalization of the Adler–Bardeen anomaly coefficient can be proven in a simple way.

To summarize, we started with a discussion of the symmetry identities of the MSSM. Particular attention was payed to the mixing of physical and unphysical fields and the implications on the external fields $`(\widehat{\mathrm{\Phi }}+\widehat{\mathrm{v}})`$, $`(A,a_\alpha ,F_A+v_A)`$ and their structure. Once the detailed form of the fields and the symmetry identities was established, it was possible to prove the renormalizability of the MSSM. On the practical side, the symmetry identities constitute an important tool. We have shown that they can be used for the unambiguous determination of possible symmetry-restoring counterterms, for calculating the gauge dependence of $`\mathrm{tan}\beta `$, and in a new approach to the non-renormalization theorems.