text:
# Universally diverging Grüneisen parameter and the magnetocaloric effect close to quantum critical points

\[

## Abstract

At a generic quantum critical point, the thermal expansion $`\alpha `$ is more singular than the specific heat $`c_p`$. Consequently, the “Grüneisen ratio”, $`\mathrm{\Gamma }=\alpha /c_p`$, diverges. When scaling applies, $`\mathrm{\Gamma }T^{1/(\nu z)}`$ at the critical pressure $`p=p_c`$, providing a means to measure the scaling dimension of the most relevant operator that pressure couples to; in the alternative limit $`T0`$ and $`pp_c`$, $`\mathrm{\Gamma }\frac{1}{pp_c}`$ with a prefactor that is, up to the molar volume, a simple universal combination of critical exponents. For a magnetic-field driven transition, similar relations hold for the magnetocaloric effect $`(1/T)T/H|_S`$. Finally, we determine the corrections to scaling in a class of metallic quantum critical points.

\]

Introduction: The anomalous behavior observed in an increasing number of systems, ranging from insulating magnets and heavy fermion compounds to cuprate superconductors, has been attributed to the presence of quantum critical points (QCPs). These occur in systems where a continuous quantum phase transition (QPT) at $`T=0`$ is induced by tuning some control parameter like pressure $`p`$, doping or magnetic field $`H`$. Such zero-temperature critical points can determine the properties of materials in a wide range of temperatures. In general, quantum critical points are more difficult to characterize compared to their classical counterparts. At a classical critical point, thermodynamic quantities typically diverge; the associated critical exponents historically played a central role in our eventual understanding of scaling and universality. Some of these divergences, however, have to disappear at a QCP: there are constraints placed by the third law of thermodynamics due to the very fact that the transition takes place at zero temperature. Here we show that the Grüneisen ratio  diverges at any QCP, in a way that provides a novel thermodynamic means of probing quantum phase transitions.

We define the Grüneisen ratio $`\mathrm{\Gamma }`$  in terms of the molar specific heat $`c_p=\frac{T}{N}\frac{S}{T}|_p`$ and the thermal expansion $`\alpha =\frac{1}{V}\frac{V}{T}|_{p,N}=\frac{1}{V}\frac{S}{p}|_{T,N}`$

$`\mathrm{\Gamma }={\displaystyle \frac{\alpha }{c_p}}={\displaystyle \frac{1}{V_mT}}{\displaystyle \frac{S/p}{S/T}}`$ (1)
where $`S`$ is the entropy and $`V_m=V/N`$ the molar volume. In ordinary situations, pressure dependences are regular and a finite Grüneisen ratio is expected as is indeed observed in all previous measurements of this quantity in the literature. Such a regular dependence is typically described by assuming that the system is dominated by a single energy scale $`E^{}`$ (e.g. the Fermi energy in a metal or the Debye energy if acoustic phonons dominate) so the molar entropy takes the form $`S/N=f(T/E^{})`$. The Grüneisen ratio is then temperature independent and given by $`\mathrm{\Gamma }=\frac{1}{V_mE^{}}\frac{E^{}}{p}`$. However, this formula already suggests that a diverging $`\mathrm{\Gamma }`$ can be expected when some energy scale $`E^{}`$ vanishes as it happens at a QCP.

Divergence of the Grüneisen ratio at QCPs: A quantum critical point is reached in a singular fashion by tuning some external parameter and, in general, this external parameter is thermodynamically coupled to pressure. In the low temperature limit, the singular terms of $`S`$ and $`T`$ in Eq. (1) cancel out leaving $`\mathrm{\Gamma }`$ to depend only on singularities associated with the pressure $`p`$. As the pressure controls the QPT, such a singularity always exists and the Grüneisen ratio diverges at any QCP. This divergence is entirely determined by the scaling dimension of the control parameter, which is the most relevant operator to which the pressure couples. As shown below, this leads to a $`T`$ dependence, $`\mathrm{\Gamma }1/T^{1/\nu z}`$ \[see Eq. (10)\]. In other words, the temperature exponent of the Grüneisen ratio provides a direct means to measure $`\nu z`$ and, as a result, characterize a QCP. Put in a slightly different way, the thermal expansion contains valuable information complementary to that obtained from the specific heat: while $`c_p`$ measures the response to $`T`$ (y-axis in Fig. 1), $`\alpha `$ describes the response to the tuning parameter of the QPT, the second relevant variable at a QCP (x-axis in Fig. 1). This has to be contrasted with a classical phase transition. There, generically, only one relevant operator exists to which both $`T`$ and $`p`$ couple. Accordingly, $`\mathrm{\Gamma }`$ will be constant close to a classical transition.

To observe the singular behavior of $`\mathrm{\Gamma }`$ or the thermal expansion the pressure has to couple sufficiently strongly to the critical dynamics. This is for example the case in heavy fermion compounds where the intricate competition between magnetic interactions and the Kondo effect can be tuned by pressure, doping or magnetic field to yield a QPT, typically from a metallic antiferromagnet to a metallic paramagnet. The high sensitivity to pressure arises from the exponential dependence of the Kondo temperature on system parameters. Whether the transitions in these systems conform to the Gaussian picture associated with $`T=0`$ spin-density wave (SDW) transitions or are non-Gaussian as in a locally quantum critical point is a question of great current interest.

If the control parameter of the QPT is not pressure but an external magnetic field $`H`$, the role of the Grüneisen ratio is played by the ratio of the $`T`$-derivative of the magnetization $`M`$ (per mole) to the molar specific heat, for either fixed pressure or fixed volume,

$`\mathrm{\Gamma }_H={\displaystyle \frac{(M/T)_H}{c_H}}={\displaystyle \frac{1}{T}}{\displaystyle \frac{(S/H)_T}{(S/T)_H}}={\displaystyle \frac{1}{T}}{\displaystyle \frac{T}{H}}|_S.`$ (2)
It can be determined directly from the magnetocaloric effect – the change of temperature in response to an adiabatic ($`S=const.`$) change of $`H`$.

In the following, we go beyond these general considerations by carrying through a) a more detailed analysis based on the assumption of scaling and b) a model study on the spin-density-wave (SDW) QCPs in metallic systems. The latter is a model system that is above or equal to the upper critical dimension, so corrections to scaling are important; it is of direct interest in comparing with experiments in heavy fermion compounds.

Scaling analysis: Close to any QCP, the correlation length $`\xi `$ diverges as a function of a control parameter $`r`$, $`\xi |r|^\nu `$, where e.g. $`r=(pp_c)/p_c`$ or $`r=(HH_c)/H_c`$. Correspondingly, a typical correlation (imaginary) time, $`\xi _\tau \xi ^z`$, diverges as the QCP is approached. The “dynamical critical exponent” $`z`$ depends on the dynamics of the order parameter and relates time and length scales.

If one assumes that the critical behavior is governed by $`\xi `$ and $`\xi _\tau `$ (a more careful discussion of this assumption is given below), the critical contribution to the free energy per mole, $`F_{cr}=FF_{\text{reg}}`$, can be cast into the following standard scaling Ansatz (using hyperscaling)

$`{\displaystyle \frac{F_{cr}}{N}}`$ $`=`$ $`\rho _0r^{\nu (d+z)}\stackrel{~}{f}\left({\displaystyle \frac{T}{T_0r^{\nu z}}}\right)`$ (3)
$`=`$ $`\rho _0\left({\displaystyle \frac{T}{T_0}}\right)^{(d+z)/z}f\left({\displaystyle \frac{r}{(T/T_0)^{1/(\nu z)}}}\right),`$ (4)
where $`\rho _0`$ and $`T_0`$ are non-universal constants, while $`f(x)`$ and $`\stackrel{~}{f}(x)`$ are universal scaling functions. Obviously, $`f(x0)f(0)+xf^{}(0)+\mathrm{}`$ is regular as there is no phase transition at $`r=0,T>0`$ (see Fig. 1). The limit $`\stackrel{~}{f}(x0)=\stackrel{~}{f}(0)+cx^{y_0+1}`$ describes the low temperature behavior of the phases to the left or right side of the QCP (in general different for $`r>0`$ and $`r<0`$). Note that the exponent $`y_0>0`$ has to be positive due to the third law of thermodynamics. It characterizes the power-law behavior of the specific heat $`c_pT^{y_0}`$, e.g. $`y_0=1`$ for a Fermi liquid, $`y_0=2`$ for a d-wave superconductor in $`d=2`$, or $`y_0=d`$ and $`d/2`$ for an insulating antiferromagnet and ferromagnet, respectively.

Thermodynamical quantities are easily obtained from (4). The critical contribution $`c_{cr}`$ to the specific heat at $`r=0`$ is given by

$$c_{cr}(T,r=0)=\frac{(d+z)d}{z^2}\frac{\rho _0}{T_0}f(0)\left(\frac{T}{T_0}\right)^{d/z}$$
(5)
and for $`T0`$, $`r0`$

$$c_{cr}(T0,r)=\frac{\rho _0cy_0(y_0+1)}{T_0}\left(\frac{T}{T_0}\right)^{y_0}r^{\nu (dy_0z)}.$$
(6)
Similarly, in the case of a pressure tuned QCP with $`r=(pp_c)/p_c`$ the critical contribution $`\alpha _{cr}`$ to the thermal expansion reads

$`\alpha _{cr}(T,r=0)`$ $`=`$ $`{\displaystyle \frac{d+z\frac{1}{\nu }}{z}}{\displaystyle \frac{\rho _0f^{}(0)}{T_0p_cV_m}}\left({\displaystyle \frac{T}{T_0}}\right)^{(d\frac{1}{\nu })/z}`$ (7)
and for $`r0`$

$`\alpha _{cr}(T0)`$ $`=`$ $`{\displaystyle \frac{\rho _0(y_0+1)c\nu (dy_0z)}{T_0V_m}}{\displaystyle \frac{r^{\nu (dy_0z)}}{p_cr}}\left({\displaystyle \frac{T}{T_0}}\right)^{y_0}.`$ (8)
The thermal expansion is more singular than the specific heat leading to a Grüneisen ratio

$`\mathrm{\Gamma }_{cr}(T,r=0)={\displaystyle \frac{\alpha _{cr}}{c_{cr}}}=G_TT^{1/(\nu z)},`$ (10)
where the prefactor $`G_T=\frac{(d+z1/\nu )zf^{}(0)}{(d+z)df(0)}\frac{T_0^{1/(\nu z)}}{p_cV_m}`$ contains some non-universal parameters ($`p_c`$ and $`T_0`$). We reach the important conclusion that the temperature exponent of the Grüneisen ratio is equal to $`1/\nu z`$.

In the other limit $`T0`$, $`r0`$ we obtain the universal result

$`\mathrm{\Gamma }_{cr}(T0,r)=G_r{\displaystyle \frac{1}{V_m(pp_c)}}.`$ (11)
Remarkably, even the (generally unknown) scaling functions cancel out in the amplitude $`G_r`$, leaving only a combination of critical exponents and the dimensionality:

$`G_r={\displaystyle \frac{\nu (dy_0z)}{y_0}}.`$ (12)
Note that the universality of this prefactor is connected to the third law of thermodynamics – a finite residual entropy per volume ($`y_0=0`$) would spoil this result.

It is rather difficult to measure thermal expansion inside a pressure cell. However, in many systems doping acts like “chemical pressure”. If doping $`x`$ and pressure $`p`$ can be quantitatively related, $`pp_c=c(xx_c)`$, a measurement of $`\mathrm{\Gamma }`$ for different samples at ambient pressure can be used to check the prediction (11) quantitatively. For generic tuning parameters, we need to substitute $`(r/p)`$ for $`1/p_c`$ in Eqs. (7,LABEL:alpha-r) and modify Eqs. (10,11) accordingly.

Similarly, for a QCP tuned by magnetic field \[$`r=(HH_c)/H_c`$\] one obtains for the magnetocaloric effect

$`\mathrm{\Gamma }_{H,cr}(T0,r)={\displaystyle \frac{(M/T)_H}{c_{cr}}}=G_r{\displaystyle \frac{1}{HH_c}}.`$ (13)
Again, in the $`T0`$ limit the prefactor (12) is universal. The $`T`$-dependence of $`\mathrm{\Gamma }_{H,cr}`$ at $`r=0`$ is also given by (10).

It is interesting to compare the above with the case of a quantum critical line, where the critical behavior is not restricted to a single point but to a finite (pressure) interval. Here, since only marginal and irrelevant operators exist for $`T=0`$, $`\mathrm{\Gamma }_{cr}`$ can diverge at most logarithmically

$`\mathrm{\Gamma }_{cr}\pm \mathrm{log}T.`$ (14)
Conversely, if $`\mathrm{\Gamma }`$ diverges algebraically for $`T0`$ a critical line scenario can be excluded.

Applicability of scaling: The applicability of the scaling results (413) depends on a number of assumptions. Most importantly, in an actual experiment not $`\mathrm{\Gamma }_{cr}=\alpha _{cr}/c_{cr}`$ but $`\mathrm{\Gamma }=\alpha /c`$ is measured and sometimes non-critical contributions can dominate (for an example see below). To verify Eqs. (10) and (11) in a situation where $`c_{cr}`$ is subleading, one would have to subtract carefully the non-critical contributions to the specific heat.

Generally, the scaling Ansatz (4) holds only below the upper critical dimension ($`d+z<4`$ within $`\mathrm{\Phi }^4`$ theories). At the upper critical dimension, logarithmic corrections to scaling arise. Above the upper critical dimension, the scaling argument can be spoiled by the presence of so called “dangerously irrelevant operators”: the free energy is a singular function of irrelevant variables. Explicit calculations (see below) for the case of an SDW transition  show that on the paramagnetic side the irrelevant operator at most leads to logarithmic corrections.

A more subtle question is whether one of the basic assumptions underlying the scaling approach (4) holds: Is there a single diverging time-scale close to the QCP? For example, in a nearly magnetic metal the answer to this question is not obvious as there are at least two types of low-energy degrees of freedom: magnetic fluctuations and fermionic quasiparticles. This can indeed lead to a breakdown of simple scaling relations as shown e.g. by Belitz et al. . In the case of a local critical point induced by a (non-local) magnetic transition, as has been suggested by one of the authors in , two scaling dimensions need to be considered: one associated with the tuning of the long-wavelength fluctuations and the other with the tuning of the local fluctuations.

SDW transitions: We now turn to more specific calculations at SDW quantum critical points, for two reasons. First, they allow us to address a number of questions concerning the scaling results: How do corrections to scaling arise at the upper critical dimension? Are the scaling results valid above this dimension? What happens, if the prefactor $`dy_0z`$ in Eqs. (11,13) vanishes? Second, our calculations are important for the purpose of assessing the relevance of SDW QCPs to the magnetic quantum phase transitions in heavy fermion compounds.

Our starting point is the Ginzburg–Landau–Wilson functional of Hertz :

$`S[\varphi ]`$ $`=`$ $`{\displaystyle \underset{𝐪,i\omega _n}{}}\left(\delta +q^2+{\displaystyle \frac{|\omega _n|}{\mathrm{\Gamma }_q}}\right)|\varphi _{𝐪,i\omega _n}|^2+S^{(4)},`$ (15)
$`S^{(4)}`$ $`=`$ $`u{\displaystyle _0^\beta }𝑑\tau {\displaystyle d^d𝐫[\varphi (𝐫,\tau )]^4},`$ (16)
with $`\mathrm{\Gamma }_q=\mathrm{\Gamma }_0q^{z2}`$, where $`z=2`$ for an antiferromagnetic SDW transition in a metal. The $`z=3`$ theory may be used to describe the critical endpoint of a metamagnetic first order transition  in $`d=2,3`$. In the case of a ferromagnetic QCP in $`d=3`$, the model (15) with $`z=3`$ is only valid up to logarithmic corrections and breaks down in $`d=2`$ . For commensurate 2D magnetism coupled to 2D fermions there are additional singularities in the fermion-collective-mode coupling; these singularities are absent when the fermions are taken to be 3D. Following the renormalization group scheme adopted by Millis  we have calculated the thermal expansion and the Grüneisen ratio for $`d=2,3`$ and $`z=2,3`$ on the non-magnetic side of the phase diagram, $`\delta \delta _c`$. Details of the calculation will be reported elsewhere.

The results are summarized in Tables I and II. Up to logarithmic corrections the results obey the scaling forms (413) with $`\nu =1/2,y_0=1`$. Note that for $`d=z`$ the prefactor in (8) and (11,13) vanishes. The $`1/r`$ dependence of $`\alpha _{cr}`$ for $`d=z`$ arises from a $`T^2\mathrm{log}1/r`$ correction to $`F_{cr}`$ not captured by scaling. For the quantum critical regime in $`d=1/\nu =2`$ the thermal expansion is logarithmic. The argument of the logarithm is a power of $`T`$ for $`d+z>4`$ and is itself logarithmically dependent on $`T`$ for $`d+z=4`$; these features reflect the dangerously irrelevant or marginal nature of the quartic coupling $`u`$.

In addition to the critical contributions, the measured quantities also contain non-critical background components. We list here the full results for the purpose of comparisons with experiments in heavy fermion compounds undergoing an antiferromagnetic transition ($`z=2`$). Consider first $`d=3`$. At the QCP ($`r=0`$)

$`\alpha `$ $`=`$ $`a_1T^{1/2}+a_2T,`$ (17)
where the $`a_2`$ term comes from the (fermionic) background contribution. However, approaching the QCP in the Fermi-liquid regime

$`\alpha `$ $`=`$ $`(a_1/r^{1/2}+a_2)T.`$ (18)
For $`d=2`$ and $`z=2`$, we have at the QCP ($`r=0`$)

$`\alpha `$ $`=`$ $`a_1\mathrm{log}[b\mathrm{log}{\displaystyle \frac{T_0}{T}}]+a_2T,`$ (19)
and in the Fermi-liquid regime approaching the QCP:

$`\alpha `$ $`=`$ $`(a_1/r+a_2)T.`$ (20)
In two dimensions, the thermal expansion at $`r=0`$ diverges in the zero temperature limit in sharp contrast to the textbook statement that $`\alpha (T0)=0`$. Still, it is straightforward to show that our results satisfy the third law of thermodynamics. As $`\alpha =(1/V)S/p`$, we can write, for generic pressure,

$`S(p,T)=S(p_c,T){\displaystyle _{p_c}^p_{}}\alpha V𝑑p{\displaystyle _p_{}^p}\alpha V𝑑p,`$ (21)
where $`p_{}`$ characterizes the crossover between the QC and FL regimes. $`S(p,T0)0`$ due to a vanishing integration region \[$`(p^{}p_c)T`$\] over which $`\alpha `$ is divergent.

We now briefly discuss the experimental implications of our results. Many heavy fermion compounds have long been known to show a Grüneisen ratio that increases to a very large value as temperature is lowered. Experiments are also becoming available in the heavy fermion metals tuned to an antiferromagnetic quantum critical point, making possible a systematic comparison with our theory. It is hoped that the present paper will stimulate similar measurements in other kinds of (real and putative) quantum critical materials.

In conclusion, we argue that the Grüneisen ratio and the magnetocaloric effect are divergent at any QCP. In addition, they can be used to measure the scaling dimensions and to check the very existence of a quantum critical point.

We would like to acknowledge helpful discussions with P. Gegenwart, H.v. Löhneysen, J. Mydosh, C. Pfleiderer, F. Steglich, and J. D. Thompson. This work has been supported in part by NSF Grant No. DMR-0090071, TCSAM, the Welch foundation (L.Z. and Q.S.) and the Emmy-Noether program of the Deutsche Forschungsgemeinschaft (M.G. and A.R.).


format: text

text:
# Scaling behavior of a nonlinear oscillator with additive noise, white and colored

## 1 Introduction

The concept of noise in Physics stems from the theory of Brownian motion in which the force exerted on a macroscopic impurity by the surrounding molecules is represented by a random function $`(t)`$ wax . This model provides a microscopic explanation of the law of diffusion and yields the celebrated fluctuation-dissipation relation vankampen . Assuming that the Brownian particle of position $`x(t)`$ and velocity $`v(t)`$ is also submitted to a harmonic potential $`𝒰(x)=\frac{1}{2}\omega ^2x^2`$, the dynamic equation associated with this system is the linear Langevin equation:

$$\frac{\text{d}}{\text{d}t}v(t)+\gamma v(t)+\omega ^2x(t)=(t),$$
(1)
where $`\gamma `$ is the effective damping coefficient.

The interplay between noise and nonlinearity generates many original phenomena that make the study of random dynamical systems a rich and fascinating field. Noise in a nonlinear system may induce non-equilibrium phase transitions lefever ; landaMc  or improve the performance of a device via stochastic resonance marchesoni . Random fluctuations may be rectified into a directed motion when a particle in a fluctuating environment is submitted to a “ratchet-like” potential reimann .

A simple nonlinear noisy dynamical system is obtained by considering Eq. (1) with an anharmonic potential $`𝒰(x)`$. If the random force is approximated by a Gaussian white noise, the stationary probability distribution function (PDF) of the energy, valid when $`t\gamma ^1`$, is given by the canonical distribution. To the best of our knowledge, the closed form of the distribution function of the nonlinear oscillator’s energy is not known when the correlation time $`\tau `$ of the random force is non-zero (see colored  and references therein for approximate results valid in the small $`\tau `$ limit).

In this work, we study the motion of an *undamped* nonlinear oscillator submitted to an additive noise. For Gaussian white noise (Section 2), we apply the method presented in a recent article philkir , where we studied the dynamical behavior of an undamped nonlinear oscillator with a fluctuating frequency represented as a parametric noise. We calculate the probability distribution function of the energy in the long-time limit and match it to the distribution obtained in presence of damping. The average energy, root-mean-square amplitude and velocity of the oscillator grow algebraically with time, with anomalous diffusion exponents different from those obtained for parametric noise. We show in Section 3 that diffusion is reduced at large times when the additive noise has a non-zero correlation time. Anomalous diffusion exponents are calculated in a self-consistent manner for an Ornstein-Uhlenbeck random force. In a somewhat counter-intuitive fashion, the white-noise behavior is observed for sufficiently *short* times. The crossover from white to colored noise can be embodied in a single scaling function as shown by dimensional analysis. All the results obtained for a nonlinear oscillator submitted to additive or multiplicative, white or colored Gaussian noise are summarized in Table 1.

## 2 The nonlinear oscillator with additive white noise

We study a particle trapped in a confining potential $`𝒰(x)`$ and subject to additive noise. We neglect dissipative effects unless stated otherwise. The dynamics of the nonlinear oscillator is given by

$$\frac{\text{d}^2}{\text{d}t^2}x(t)=\frac{𝒰(x)}{x}+\xi (t),$$
(2)
where $`x(t)`$ represents the position of the oscillator at time $`t`$. In this Section, 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^{}).`$ (3)
For analytical as well as numerical calculations, we shall interpret stochastic differential equations such as (2) according to Stratonovich’s convention. Indeed, Stratonovich calculus appears naturally in the limit of a vanishingly small correlation time vankampen . It must be emphasized that this convention (as opposed to Ito calculus) introduces correlations between the noise and the dynamical variables at the same time $`t`$ itostrato . Moreover, the mechanical oscillator defined by Equations (2,3) is *not* conservative: noise feeds energy into the system (see Eqs. (15) and (35)). In the absence of damping, the oscillator’s amplitude increases with time.

For the potential $`𝒰(x)`$ to be confining, we must have $`𝒰+\mathrm{}`$ when $`|x|\mathrm{}`$. We shall restrict our analysis to the case where $`𝒰(x)`$ is a polynomial, and in order to respect the $`xx`$ symmetry, $`𝒰(x)`$ is taken to be even in $`x`$. Hence, when $`|x|\mathrm{}`$,

$$𝒰\frac{x^{2n}}{2n}\text{ with }n2,$$
(4)
where the coefficient of $`x^{2n}`$ has been set to $`1/(2n)`$ after a suitable rescaling. As the amplitude $`x`$ of the oscillator grows at large times, only the asymptotic behavior of $`𝒰(x)`$ for $`|x|\mathrm{}`$ is relevant and thus Eq. (2) reduces to

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)^{2n1}=\xi (t).$$
(5)
The method we shall use to study Eqs. (3,5) is akin to that presented in our previous work on nonlinear oscillators with multiplicative noise philkir . In Section 2.1, we use the integrability properties of the deterministic nonlinear oscillator to write exact stochastic differential equations in energy-angle variables. We then derive equipartition relations in Section 2.2. Averaging out the fast angular variable strato ; lindenberg , we calculate the energy’s probability distribution function in the long-time limit (Section 2.3). In Section 2.4, we obtain the associated anomalous diffusion exponents and constants and study their dependence on the stiffness of the confining potential at infinity. Our numerical simulations are based on a time discretization presented in mannella , and described in detail in philkir .

### 2.1 Energy-angle coordinates

We rewrite the second-order stochastic differential equation (5) as a first order system in energy (or action) and angle variables. Without noise, the deterministic version of Eq. (5):

$$\ddot{x}+x^{2n1}=0,$$
(6)
is integrable because of energy conservation, where the energy is defined by

$$E=\frac{1}{2}\dot{x}^2+\frac{1}{2n}x^{2n}.$$
(7)
The action-angle variables $`(I,\varphi )`$ associated with Eq. (6) are

$`I`$ $`=`$ $`4{\displaystyle _0^{(2nE)^{\frac{1}{2n}}}}\sqrt{2E{\displaystyle \frac{x^{2n}}{n}}}\text{d}x`$ (8)
$`=`$ $`4(2^{n+1}n)^{\frac{1}{2n}}E^{\frac{n+1}{2n}}{\displaystyle _0^1}\sqrt{1u^{2n}}\text{d}u`$
$``$ $`E^{\frac{n+1}{2n}},`$
$`\varphi `$ $`=`$ $`\sqrt{n}{\displaystyle _0^{x/(2nE)^{1/2n}}}{\displaystyle \frac{\text{d}u}{\sqrt{1u^{2n}}}},`$ (9)
where the angle variable $`\varphi `$ is defined modulo the oscillation period $`4K_n`$, with

$$K_n=\sqrt{n}_0^1\frac{\text{d}u}{\sqrt{1u^{2n}}}.$$
(10)
The solution of Eq. (6) can be parametrized as a function of energy $`E`$ and angle $`\varphi `$

$`x`$ $`=`$ $`E^{1/2n}𝒮_n(\varphi ),`$ (11)
$`\dot{x}`$ $`=`$ $`(2n)^{\frac{n1}{2n}}E^{1/2}𝒮_n^{}(\varphi ),`$ (12)
where the hyper-elliptic function $`𝒮_n`$ is defined as

$`𝒮_n(\varphi )=Y\varphi `$ $`=`$ $`\sqrt{n}{\displaystyle _0^{\frac{Y}{(2n)^{1/2n}}}}{\displaystyle \frac{\text{d}u}{\sqrt{1u^{2n}}}}`$ (13)
$`=`$ $`{\displaystyle \frac{\sqrt{n}}{(2n)^{1/2n}}}{\displaystyle _0^Y}{\displaystyle \frac{\text{d}u}{\sqrt{1\frac{u^{2n}}{2n}}}}.`$
From this definition, we obtain the following relation between $`𝒮_n`$ and its derivative $`𝒮_n^{}`$

$$𝒮_n^{}(\varphi )=\frac{(2n)^{\frac{1}{2n}}}{\sqrt{n}}\left(1\frac{(𝒮_n(\varphi ))^{2n}}{2n}\right)^{\frac{1}{2}}.$$
(14)
Writing Eq. (5) in terms of energy and angle variables, we obtain the following stochastic dynamical system

$`\dot{E}`$ $`=`$ $`\dot{x}\xi (t)=(2n)^{\frac{n1}{2n}}E^{\frac{1}{2}}𝒮_n^{}(\varphi )\xi (t),`$ (15)
$`\dot{\varphi }`$ $`=`$ $`(2nE)^{\frac{n1}{2n}}{\displaystyle \frac{1}{(2n)^{\frac{1}{2n}}}}{\displaystyle \frac{𝒮_n(\varphi )}{(2nE)^{\frac{1}{2}}}}\xi (t).`$ (16)
With the help of the auxiliary variable $`\mathrm{\Omega }`$

$$\mathrm{\Omega }=(2n)^{\frac{n+1}{2n}}E^{\frac{1}{2}},$$
(17)
Eqs. (15) and (16) are written in the simpler form:

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`n𝒮_n^{}(\varphi )\xi (t),`$ (18)
$`\dot{\varphi }`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{(2n)^{\frac{1}{2n}}}}\right)^{\frac{n1}{n}}{\displaystyle \frac{𝒮_n(\varphi )}{\mathrm{\Omega }}}\xi (t).`$ (19)
These equations have been derived without any hypothesis on the function $`\xi `$ and are rigorously equivalent to the original equation (5).

Although the method is identical, this set of equations differs from that obtained for parametric noise (See Eqs. (27)-(28) in philkir ). As a consequence, additive and multiplicative random noises lead to different long-time behaviors of the observables of the system (See table 1).

### 2.2 Equipartition relations

Noting that $`𝒮_n`$ and $`𝒮_n^{}`$ are bounded functions of $`\varphi `$, we deduce from Eqs. (18) and (19) that $`\mathrm{\Omega }`$ is a diffusive variable that scales as $`t^{1/2}`$, whereas $`\varphi `$ is a fast variable that scales like $`t^{3/21/2n}`$ (these assertions will be justified rigorously in the next section). In the long time limit, it is therefore justified to average the dynamics over the variations of the fast variable $`\varphi `$. Using Eqs. (11)-(12), and assuming that $`\varphi `$ is uniformly distributed over the interval $`[0,4K_n]`$ of a period, we find the following relations:

$`E`$ $`=`$ $`{\displaystyle \frac{n+1}{2n}}\dot{x}^2,`$ (20)
$`\dot{x}^2`$ $`=`$ $`x^{2n}.`$ (21)
Figure 1 shows that Eqs. (20)-(21) are indeed verified with excellent accuracy.

Because the transformation of variables (11)-(12) is independent of the nature of the noise $`\xi (t)`$, the statistical identities (20)-(21) are the same as those obtained for parametric noise philkir . Such generalized equipartition relations are valid whenever the asymptotic probability distribution function is independent of the angle variable and is a function of the energy only. This is the case, in particular, of the canonical Boltzmann-Gibbs distribution: Eqs. (20)-(21) are valid at thermodynamic equilibrium, as can be shown using standard techniques of statistical mechanics in the canonical ensemble.

### 2.3 Calculation of the effective PDF

The Fokker-Planck equation corresponding to the system (18, 19) reads

$`_t`$ $`P`$ $`=\left({\displaystyle \frac{\mathrm{\Omega }}{(2n)^{\frac{1}{2n}}}}\right)^{(n1)/n}_\varphi P`$
$`+`$ $`{\displaystyle \frac{𝒟}{2}}(_\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 }f(\varphi )P\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))`$
where we have defined

$$f(\varphi )=n𝒮_n^{}(\varphi )\text{ and }g(\varphi )=𝒮_n(\varphi ).$$
(23)
This equation governs the dynamics of the probability distribution function $`P_t(\mathrm{\Omega },\varphi )`$. As we already noticed, the angle variable $`\varphi `$ varies rapidly as $`\mathrm{\Omega }`$ grows. We carry out the averaging of Eq. (2.3) under the hypothesis that the probability density, $`P_t(\mathrm{\Omega },\varphi )`$, becomes uniform in $`\varphi `$ when $`t\mathrm{}`$. The averaged evolution equation of the reduced probability density $`\stackrel{~}{P}_t(\mathrm{\Omega })`$ is then

$$_t\stackrel{~}{P}=\overline{f^2(\varphi )}\frac{𝒟}{2}\left(_\mathrm{\Omega }^2\stackrel{~}{P}\frac{1}{n}_\mathrm{\Omega }\frac{\stackrel{~}{P}}{\mathrm{\Omega }}\right),$$
(24)
where we have used $`f=(_\varphi g)/n`$. The notation $`\overline{f^2(\varphi )}`$ represents the mean value of the function $`f^2`$ over one period $`4K_n`$ of $`\varphi `$. This constant has been calculated explicitly in philkir  and is given by

$$\overline{f^2(\varphi )}=n^2\overline{𝒮_n^{}(\varphi )^2}=\frac{n^2(2n)^{\frac{1}{n}}}{n+1}.$$
(25)
The Fokker-Planck equation (24) is equivalent to an effective first-order stochastic differential equation for the slow variable $`\mathrm{\Omega }`$. Defining an effective Gaussian white noise $`\stackrel{~}{\xi }(t)`$ such that

$$\stackrel{~}{\xi }(t)\stackrel{~}{\xi }(t^{})=\stackrel{~}{𝒟}\delta (tt^{})\text{ with }\stackrel{~}{𝒟}=\frac{n^2(2n)^{\frac{1}{n}}}{n+1}𝒟,$$
(26)
we deduce from Eq. (24) the effective Langevin equation for $`\mathrm{\Omega }`$

$$\dot{\mathrm{\Omega }}=\frac{\stackrel{~}{𝒟}}{2n}\frac{1}{\mathrm{\Omega }}+\stackrel{~}{\xi }(t).$$
(27)
Thus $`\mathrm{\Omega }`$ can be reinterpreted as a Brownian variable in a logarithmic potential.

The scale-invariant effective Fokker-Planck equation (24) can be solved explicitly. We find

$$\stackrel{~}{P}_t(\mathrm{\Omega })=\frac{2}{\mathrm{\Gamma }\left(\frac{n+1}{2n}\right)}\frac{\mathrm{\Omega }^{\frac{1}{n}}\mathrm{e}^{\mathrm{\Omega }^2/(2\stackrel{~}{𝒟}t)}}{(2\stackrel{~}{𝒟}t)^{\frac{n+1}{2n}}},$$
(28)
where $`\mathrm{\Gamma }`$ is the Euler Gamma function. Using Eq. (17), we deduce from Eq. (28) the asymptotic probability distribution function of the energy variable:

$$\stackrel{~}{P}_t(E)=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2n}\right)}\frac{1}{E}\left(\frac{(2n)^{\frac{n+1}{n}}E}{2\stackrel{~}{𝒟}t}\right)^{\frac{n+1}{2n}}\mathrm{exp}\left\{\frac{(2n)^{\frac{n+1}{n}}E}{2\stackrel{~}{𝒟}t}\right\}$$
(29)
Moreover, the P.D.F. (29) can be matched to the canonical distribution. Suppose that some damping is present in the system. The dynamical equation (5) now reads

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+\gamma \frac{\text{d}}{\text{d}t}x(t)+x(t)^{2n1}=\xi (t),$$
(30)
where $`\gamma `$ is the friction coefficient. The stationary distribution associated with the stochastic equation (30) is:

$$\stackrel{~}{P}_{\text{can}}(E)=\frac{1}{\mathrm{\Gamma }\left(\frac{n+1}{2n}\right)}\frac{1}{E}\left(\frac{2\gamma E}{𝒟}\right)^{\frac{n+1}{2n}}\mathrm{exp}\left\{\frac{2\gamma E}{𝒟}\right\}.$$
(31)
This PDF becomes the Boltzmann-Gibbs canonical measure when supplemented with the fluctuation-dissipation relation vankampen . The crossover time $`t_c`$ from the asymptotic distribution function for the undamped oscillator to the Gibbs measure is found by matching Eq. (29) with Eq. (31):

$$t_c=\frac{n+1}{2n\gamma }.$$
(32)
Thus, there are three distinct time scales involved: a fast time scale over which the angle variable $`\varphi `$ gets uniformly distributed, a ‘long’ time scale over which the effective Fokker-Planck equation for $`\mathrm{\Omega }`$ (24) is relevant, and finally times longer than the crossover time $`t_c`$ beyond which the inevitable damping in the system takes over.

### 2.4 Time-asymptotic behavior of observables

With the probability distribution function (29), we can calculate the long time behavior, in the absence of damping, of the expectation value of any observable that depends on the position and velocity of the particule. In particular, we obtain

$`x^2`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3}{2n}\right)}{\mathrm{\Gamma }\left(\frac{1}{2n}\right)}}\left({\displaystyle \frac{2n^2}{n+1}}𝒟t\right)^{\frac{1}{n}},`$ (33)
$`\dot{x}^2`$ $`=`$ $`{\displaystyle \frac{n𝒟}{n+1}}t,`$ (34)
$`E`$ $`=`$ $`{\displaystyle \frac{𝒟}{2}}t.`$ (35)
From Eq. (33), we notice that the particle diffuses more slowly in a stiffer potential well, as one would intuitively expect. The equipartition relation between energy and velocity obtained in Section 2.2 is confirmed by Eqs. (34,35). For *all* potentials growing algebraically at infinity, we find that the mean value of the energy grows linearly with time, with a diffusion constant equal to $`𝒟/2`$ (Eq. (35)). Such a remarkable and simple behavior was not obvious a priori.

These predictions are confirmed by numerical simulations. The universal behavior of $`E`$ is represented in Fi–gure 2.a. Using a numerical value of the noise amplitude $`𝒟=1`$ in Eq. (33), we obtain

$`\text{For }n=2,`$ $`x^2=0.552t^{1/2},`$ (36)
$`\text{For }n=3,`$ $`x^2=0.526t^{1/3},`$ (37)
$`\text{For }n=4,`$ $`x^2=0.500t^{1/4}.`$ (38)
These formulae are in excellent agreement with the numerical results displayed in Figures 2.b and 2.c.

The validity of Eqs. (33,34,35) can be checked in the linear case, $`n=1`$. The exact solution of the equation

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

$`x(t)`$ $`=`$ $`{\displaystyle \frac{1}{\omega }}{\displaystyle _0^t}\mathrm{sin}\left(\omega (tu)\right)\xi (u)\mathrm{du},`$ (40)
$`\dot{x}(t)`$ $`=`$ $`{\displaystyle _0^t}\mathrm{cos}\left(\omega (tu)\right)\xi (u)\mathrm{du}.`$ (41)
We obtain the following exact formulae vankampen :

$`x^2`$ $`=`$ $`{\displaystyle \frac{𝒟}{2\omega ^2}}t\left(1{\displaystyle \frac{\mathrm{sin}2\omega t}{2\omega t}}\right),`$ (42)
$`\dot{x}^2`$ $`=`$ $`{\displaystyle \frac{𝒟}{2}}t\left(1+{\displaystyle \frac{\mathrm{sin}2\omega t}{2\omega t}}\right),`$ (43)
$`E`$ $`=`$ $`{\displaystyle \frac{𝒟}{2}}t.`$ (44)
In the long time limit $`t\mathrm{}`$, these expressions agree with Eqs. (33,35) when one substitutes $`n=1`$ (and $`\omega =1`$).

## 3 The nonlinear oscillator with additive colored noise

In order to model an additive Gaussian noise with a non vanishing correlation time $`\tau `$, we use an Ornstein-Uhlenbeck process, denoted by $`\eta (t)`$. We thus replace Eq. (5) by the following equation

$$\frac{\text{d}^2}{\text{d}t^2}x(t)+x(t)^{2n1}=\eta (t).$$
(45)
where the random noise $`\eta (t)`$ is obtained from

$$\frac{\text{d}\eta (t)}{\text{d}t}=\frac{1}{\tau }\eta (t)+\frac{1}{\tau }\xi (t),$$
(46)
and $`\xi (t)`$ is the Gaussian white noise defined in Eq. (3). The stationary statistical properties of $`\eta `$ are given by

$`\eta (t)`$ $`=`$ $`0,`$
$`\eta (t)\eta (t^{})`$ $`=`$ $`{\displaystyle \frac{𝒟}{2\tau }}\mathrm{e}^{|tt^{}|/\tau }.`$ (47)
We rewrite Eq. (45) in $`(\mathrm{\Omega },\varphi )`$ coordinates and obtain as above \[see Eqs. (18) and (19)\]

$`\dot{\mathrm{\Omega }}`$ $`=`$ $`n𝒮_n^{}(\varphi )\eta (t),`$ (48)
$`\dot{\varphi }`$ $`=`$ $`\left({\displaystyle \frac{\mathrm{\Omega }}{(2n)^{\frac{1}{2n}}}}\right)^{\frac{n1}{n}}{\displaystyle \frac{𝒮_n(\varphi )}{\mathrm{\Omega }}}\eta (t).`$ (49)
Again, the angle $`\varphi `$ is a fast variable: we expect that the equipartition relations (20) and (21) are also valid for colored additive noise. This is indeed confirmed by numerical simulations (See Figure 3).

The system (46,48,49) is equivalent to a Fokker-Planck equation for the P.D.F. $`P_t(\mathrm{\Omega },\varphi ,\eta )`$. However, averaging this equation over $`\varphi `$ does not lead to any conclusive result. The qualitative reason for this failure is as follows: when the period of the angular variable becomes smaller than the coherence time $`\tau `$ of the noise, the noise itself is averaged out.

The goal of this section is to introduce alternative methods to determine the scaling behavior of nonlinear oscillators with colored additive noise, first qualitatively in Section 3.1 thanks to dimensional analysis, then in a more rigorous manner in Section 3.2 using a self-consistent argument. Finally, we show in Section 3.3 that the crossover between the behaviors typical of white noise and colored noise can be embodied in a single scaling function.

### 3.1 Scaling analysis

The equations (45) and (46) can be interpreted as a system of three coupled first-order stochastic differential equations in presence of white noise. The corresponding Fokker-Planck equation for the P.D.F. $`P_t(x,v,\eta )`$ is given by

$$\frac{P}{t}=v\frac{P}{x}+(x^{2n1}\eta )\frac{P}{v}+\frac{1}{\tau }\frac{\eta P}{\eta }+\frac{𝒟}{2\tau ^2}\frac{^2P}{\eta ^2}.$$
(50)
If we perform a scaling analysis in the spirit of yp  by comparing terms in this equation two by two, we find that a consistent balance is found if $`v^2x^{2n}`$, $`x^{2n1}\eta `$ and $`\eta ^2𝒟t/2\tau ^2`$. Eliminating $`\eta `$, we obtain the following scaling laws

$`E`$ $``$ $`\left({\displaystyle \frac{𝒟t}{2\tau ^2}}\right)^{\frac{n}{(2n1)}},`$
$`\dot{x}`$ $``$ $`\left({\displaystyle \frac{𝒟t}{2\tau ^2}}\right)^{\frac{n}{2(2n1)}},`$
$`x`$ $``$ $`\left({\displaystyle \frac{𝒟t}{2\tau ^2}}\right)^{\frac{1}{2(2n1)}}.`$ (51)
A similar argument was used in philkir  in the case of colored *multiplicative* noise, and yielded different exponents (See Table 1).

Numerical simulations confirm the scalings predicted by Eqs. (51) for general $`n`$ (See Figure 4). Besides, these scalings are consistent with exact results found in the linear case $`n=1`$. In the next section we present an analytic derivation of the exponents appearing in Eq. (51).

### 3.2 A self-consistent calculation

We make the a priori Ansatz that, in the long-time limit, $`\mathrm{\Omega }`$ grows algebraically with time as

$$\mathrm{\Omega }t^\alpha ,$$
(52)
and we will determine the scaling exponent $`\alpha `$ from a self-consistent calculation. We deduce from Eq. (49) that $`\varphi t^\nu `$ where

$$\nu =\frac{n1}{n}\alpha +1.$$
(53)
Substituting this scaling of $`\varphi `$ in Eq. (48), we obtain (hereafter we shall leave aside all proportionality constants)

$$\mathrm{\Omega }_0^t\mathrm{dz}𝒮_\mathrm{n}^{}(\mathrm{z}^\nu )\eta (\mathrm{z}).$$
(54)
We must now determine the asymptotic statistical behavior of this expression in the $`t\mathrm{}`$ limit. In order to simplify the calculation, we replace the Ornstein-Uhlenbeck noise by a discrete dichotomous noise, where time is discretized in steps of duration $`\tau `$ (see chris  for a general discussion). This approximation amounts to replacing exponentially decaying correlations, by finite time correlations and therefore leaves the diffusion exponents unchanged. Equation (54) then reduces to

$$\mathrm{\Omega }\underset{k=0}{\overset{t/\tau }{}}ϵ_k_{k\tau }^{(k+1)\tau }\mathrm{dz}𝒮_\mathrm{n}^{}(\mathrm{z}^\nu ),$$
(55)
where the variable $`ϵ_k`$ takes the value $`\pm \sqrt{𝒟/(2\tau )}`$ randomly during the time interval $`[k\tau ,(k+1)\tau ]`$; in other terms,

$$ϵ_kϵ_l=\frac{𝒟}{2\tau }\delta _{kl}.$$
(56)
From Eqs. (55) and (56), we deduce that

$$\mathrm{\Omega }^2\underset{k=0}{\overset{t/\tau }{}}\left(_{k\tau }^{(k+1)\tau }\mathrm{dz}𝒮_\mathrm{n}^{}(\mathrm{z}^\nu )\right)^2.$$
(57)
Integrating by parts we obtain

$`{\displaystyle _{k\tau }^{(k+1)\tau }}`$ $`\mathrm{dz}𝒮_\mathrm{n}^{}(\mathrm{z}^\nu )={\displaystyle \frac{𝒮_\mathrm{n}\left((\mathrm{k}+1)^\nu \tau ^\nu \right)}{\nu \left((\mathrm{k}+1)\tau \right)^{\nu 1}}}`$
$`{\displaystyle \frac{𝒮_n\left(k^\nu \tau ^\nu \right)}{\nu \left(k\tau \right)^{\nu 1}}}+{\displaystyle \frac{\nu 1}{\nu }}{\displaystyle _{k\tau }^{(k+1)\tau }}\mathrm{dz}{\displaystyle \frac{𝒮_\mathrm{n}(\mathrm{z}^\nu )}{\mathrm{z}^\nu }}`$
The integral term on the r.h.s. of (3.2) is of order $`𝒪(k^\nu )`$. Since $`𝒮_n`$ is a bounded function, the first two terms are of order $`𝒪(k^{\nu +1})`$, and will dominate when $`k1`$ since $`\nu >1`$. We thus obtain

$$\mathrm{\Omega }^2\underset{k=1}{\overset{t/\tau }{}}\frac{1}{\left(k\tau \right)^{2\nu 2}}t^{32\nu }.$$
(59)
Assuming that the variable $`\mathrm{\Omega }`$ is not multifractal, we know from the scaling hypothesis, Eq. (52), that $`\mathrm{\Omega }^2t^{2\alpha }`$ Comparing with Eq. (59), we deduce that

$$2\alpha =32\nu .$$
(60)
Equations (53) and (60) provide the required self-consistent condition and allow to calculate the exponents:

$$\alpha =\frac{n}{2(2n1)}\text{ and }\nu =\frac{5n3}{2(2n1)}.$$
(61)
We can now deduce from Eq. (17) and from the equipartition relations, Eqs. (20) and (21), the scaling behavior of the dynamical variables

$`E`$ $``$ $`t^{\frac{n}{2n1}},`$
$`v`$ $``$ $`t^{\frac{n}{2(2n1)}},`$
$`x`$ $``$ $`t^{\frac{1}{2(2n1)}}.`$ (62)
This result is in agreement with the scaling analysis of the previous section \[see Equation (51)\].

### 3.3 Crossover from white to colored noise

We shall now match the white and the colored noise behaviors obtained in the long-time limit through a scaling function and will describe the crossover from the first regime to the second one. In the system (45) with Ornstein-Uhlenbeck noise, one can define two dimensionless variables: $`t/\tau `$ and $`𝒟\tau ^{\frac{3n1}{n1}}`$. Any dimensionless quantity can be written as a function of these two variables (or of any two independent combinations of these two variables). For example, the mean value of the energy can be written in full generality, as a function of time, as

$$E=\frac{𝒟t}{2}\mathrm{\Phi }(𝒟\tau ^{\frac{2n}{n1}}t,𝒟\tau ^{\frac{3n1}{n1}}),$$
(63)
where the prefactor $`𝒟t/2`$ ensures that the equation is dimensionally correct and has been chosen to match the white noise limit. The scaling function $`\mathrm{\Phi }`$ allows to interpolate between white and colored noise and exhibits the following asymptotic behavior:

* The white noise limit corresponds to $`\tau =0`$, therefore $`\mathrm{\Phi }(0,0)=1`$.
* When $`\tau `$ is finite (colored noise case) and $`t\mathrm{}`$, we deduce from Eq. (51) that, for $`u\mathrm{}`$ and for $`v`$ finite, $`\mathrm{\Phi }(u,v)u^{\frac{1n}{2n1}}\varphi (v)`$. The prefactor $`\varphi (v)`$ is independent of time and cannot be determined by dimensional analysis.

The crossover between white noise and colored noise is obtained when $`u=𝒟\tau ^{\frac{2n}{n1}}t1,`$ i.e., for a typical time of the order of

$$t_{\mathrm{cross}}\frac{\tau ^{\frac{2n}{n1}}}{𝒟}.$$
(64)
At this crossover time the energy is of the order of $`E𝒟t_{\mathrm{cross}}\tau ^{\frac{2n}{n1}}`$. The period of the deterministic oscillator for such an energy is given by $`TE^{\frac{n1}{2n}}\tau `$ (see Eq. (16)). Thus, the crossover occurs when the period of the deterministic oscillator is of the order of the correlation time of the noise. This fact has the following intuitive explanation. When $`tt_{\mathrm{cross}}`$, the period of the deterministic oscillator is very large compared to the correlation time of the noise: the noise is totally uncorrelated over one period and can be viewed as a sequence of uncorrelated kicks and thus acts as a white noise. However, when $`tt_{\mathrm{cross}}`$, the period of the deterministic oscillator is small compared to the correlation time of the noise: the noise is highly correlated over a period and can be viewed as an almost constant quantity over this time scale. Many periods must be added before the effect of the noise becomes perceptible: hence the diffusion slows down.

We remark that in Eq. (63) we have not used $`t/\tau `$ as the dimensionless time variable which would have been the most natural choice. In fact, the white noise limit $`\tau 0`$, and the long time limit $`t\mathrm{}`$ with finite $`\tau `$ are indistinguishable if one uses the variable $`t/\tau `$. For this reason, it may be incorrectly stated that ‘in the long time limit, colored noise appears as white’. We have seen in our problem that exactly the opposite is true: the effect of a non-zero correlation time becomes relevant at long times. The limit $`t/\tau \mathrm{}`$ is singular and its value depends on whether the second dimensionless variable $`v=𝒟\tau ^{\frac{3n1}{n1}}`$ is equal to zero or is strictly positive. By choosing the variable $`u=𝒟\tau ^{\frac{2n}{n1}}t`$ as the dimensionless time variable in Eq. (63) rather than $`t/\tau `$ we avoid this difficulty because the white noise and the colored noise cases now appear as different limits.

## 4 Conclusion

A Hamiltonian nonlinear oscillator submitted to random internal noise (additive noise) exhibits anomalous diffusion phenomena. We have calculated analytically the associated scaling exponents in the case of white noise and have obtained, in the long time limit, an explicit expression for the probability distribution function in phase space. The energy of such an oscillator grows linearly with time irrespective of the form of the confining potential. It is remarkable that even the rate of energy growth does not depend on the stiffness of the potential and is simply proportional to the amplitude of the noise. Our results also describe the behavior of an oscillator with small damping rate $`\gamma `$ for times less than the dissipation time $`\gamma ^1`$. We have shown that the probability distribution function of the energy of the undamped oscillator can be matched with the canonical Boltzmann-Gibbs distribution when $`t\gamma ^1`$.

In the case of colored noise, the scaling exponents are modified. Their values have been determined by a self-consistent calculation and the predictions of dimensional analysis have been confirmed. The diffusion exponents are reduced because the coherence of the noise over a period of the system makes the energy transfer less efficient. The growth of energy is slower than linear and in the limit of a very stiff potential, the energy grows only as the square-root of time.

In Table 1, we summarize and compare the results that we have derived for additive noise and multiplicative noise, in this work and in philkir , respectively. In the case of white noise, precise results are available: thanks to the averaging technique, the asymptotic probability distribution functions are known. For colored noise, only the anomalous diffusion exponents have so far been calculated. We notice that in all cases the exponent for the velocity is half the exponent for the energy and $`n`$ times the exponent for the amplitude i.e., the ratios between corresponding exponents are independent of the problem considered. This is the consequence of the generalized equipartition relations, which are independent of the nature of the noise (in fact, equipartition also provides universal relations between the generalized diffusion constants appearing as prefactors in the scaling laws). Besides, it was argued in philkir  that, in the long time limit, the multiplicative noise always dominates over the additive noise and indeed we observe that the diffusion exponents for the multiplicative noise are always larger than those for the additive noise.

Although we have been able here to derive the exponents for additive colored noise, the generalized diffusion constants, as well as the long time asymptotics of the probability distribution function were not found by our approach and their calculation remains an open problem. Besides, the precise calculation of the scaling function that describes the crossover between white and colored noise is certainly a challenging problem.

###### Acknowledgements.

We would like to thank Jean-Marc Luck for his constant interest in our work, as well as Michel Bauer, Peter Holdsworth and Kunihiko Kaneko for useful comments.


text:
# Rapidity and Transverse Momentum Distributions of DIS 1-jet Inclusive Cross Section — A Next-to-Leading Order Monte-Carlo Prediction

## 1 Introduction

In contrast to “fixed order” calculations, Monte-Carlo (MC) event generators attempt to predict the bulk of the exclusive components of perturbative QCD subprocesses and to give the details of the final states. While this is at the expense of model dependence in the non-perturbative hadronization, it is also the case that MC event generators incorporate a useful approximation to the resummation of large logarithms in their parton showering. This provides a successful phenomenology of infrared sensitive cross sections.

Now, it is quite hard to treat non-leading corrections in MC event generators. However, this is also quite essential. For example, in deep inelastic scattering (DIS) at small $`x_{\mathrm{bj}}`$ or in diffractive DIS, the gluon density is substantially larger than the quark densities, and the next-to-leading order (NLO) gluon-induced process is comparable to the leading order (LO) quark-induced process. Therefore, to get a sensible phenomenology, the gluon-induced NLO corrections must be included.

Precisely because MC event generators treat exclusive components of the cross section, standard “fixed order” calculations of hard scattering coefficients are not well adapted for use in event generators. Some symptoms of this are that in LO, the coefficients contain $`\delta `$-functions of parton kinematics, and that in NLO the coefficients are singular distributions (generalized functions). They only give results that can be compared with hadronic level cross sections after an integration over a wide range of parton kinematics, which provides a kind of infrared safety. Standard fixed-order calculations approximate the partons as on-shell with zero transverse momentum. For cross sections that do not include a wide integration over parton kinematics, only a correct treatment of parton kinematics allows a correct conversion of the generalized functions to ordinary functions of kinematic variables. In fact, the standard approximations used in hard-scattering coefficients cause complications and difficulties in the derivation of the algorithms used in event generators.

In this paper, we will examine the rapidity and transverse momentum distribution of the 1-jet inclusive cross section in DIS. Even though this is a jet cross section, it is not infra-red safe because its LO hard scattering coefficient is a $`\delta `$-function, and in the approximation that parton virtuality and transverse momentum are ignored, the jet momentum can be calculated from the measured 4-momentum of the scattered electron. Thus the $`\delta `$-function appears in the hadron-level jet cross section, in this simple parton-model approximation. Similarly, the conventional NLO corrections are singular at the parton-model values. Because we only integrate over a narrow range of parton kinematics around these singularities, it is not clear how to apply these conventional NLO hard scattering coefficients to this cross section. However, it is crucial for a good NLO MC method to consistently treat this kind of cross section with NLO accuracy.

As we argued in , a proper treatment of higher order (NLO and beyond) corrections in event generators requires (a) that the parton kinematics should be treated more exactly and (b) that the hard scattering coefficients should be ordinary functions, not generalized functions. This is in contrast to other proposals  for matrix element corrections to event generators. In  our approach is proposed in general and is designed to be generalizable to all orders in perturbation theory.

More concretely, in  this method was applied to DIS in a form specifically designed for use in event generators using the Bengtsson–Sjöstrand  algorithm for showering. Such event generators include RAPGAP  and PYTHIA . We have now implemented111Only gluon-induced NLO corrections are implemented, as in . The quark-induced NLO terms have the complication of a soft divergence, which we plan to deal with in a later paper. Gluon-induced NLO corrections are important in DIS, because they are not suppressed compared to the LO terms whenever the gluon distribution is large.  this method in RAPGAP, a widely used event generator for both inclusive and diffractive DIS. The implementation is based on earlier work by Schilling .

In this paper we show the results of the 1-jet inclusive cross section calculated by our subtraction method. We also compare our method with a method previously implemented in RAPGAP, which uses a $`k_t`$-cut applied to the simple unsubtracted NLO cross section. The $`k_t`$-cut method is intended mainly for the di-jet rate calculation of high $`p_t`$ jet events.

We also compare our results to those of a LO calculation from RAPGAP, where the $`\delta `$-function in the 1-jet inclusive cross section is smeared out only by the parton showering and the hadronization. At small $`x_{\mathrm{bj}}`$ the gluon-induced NLO contributions are enhanced, and we expect the cross section with NLO accuracy to be a very broad distribution, and the three methods should give very different results. At large $`x_{\mathrm{bj}}`$ the NLO corrections are suppressed by $`O(\alpha _s)`$, we expect to see similar results by these methods.

The organization of this paper is the following. In Sec. 2 we briefly describe the subtraction method and the $`k_t`$-cut method in RAPGAP. In Sec. 3 we show and discuss the results of the 1-jet inclusive cross sections for DIS at small $`x_{\mathrm{bj}}`$ and large $`x_{\mathrm{bj}}`$ predicted by the three different methods. In Sec. 4, we discuss the possible applications and future work of the NLO RAPGAP with subtraction method.

## 2 NLO hard scattering

In , we devised a subtraction method for NLO corrections, and obtained the coefficient function for the theoretically simple yet phenomenologically significant gluon-induced correction to DIS. It was assumed that the Bengtsson-Sjöstrand (BS)  algorithm is used for the initial state parton shower, as in RAPGAP. A point-by-point subtraction in the phase space ensures that the hard-scattering coefficient is non-singular . The pdfs needed are specific to the showering algorithm , and those needed to match the BS algorithm were calculated in . The quark-induced NLO corrections have the complication of a soft divergence, and we plan to treat this case, matched to the BS algorithm, in a later paper; the quark-induced term is a genuine correction, order $`\alpha _s`$ compared with the LO process, whereas the gluon-induced term is not suppressed at all when the gluon density is large.

In RAPGAP, the previously implemented method for NLO was the $`k_t`$-cut method, designed mainly to give the correct di-jet rate at large $`Q^2`$. Here the unsubtracted parton-level matrix elements of the NLO processes are used with a cut on the transverse momentum $`k_t`$ of the out going partons in the center of mass frame of $`\gamma ^{}+\mathrm{parton}`$ scattering. $`k_{t\mathrm{cut}}`$ is required to be large enough so that the NLO contributions are smaller than the total cross section. Usually $`k_{t\mathrm{cut}}^2=4\mathrm{GeV}^2`$. This gives us a large logarithm $`\mathrm{log}(k_{t\mathrm{cut}}^2/Q^2)`$ in the integrated NLO cross section when $`Q^2`$ is large.

## 3 Rapidity and Transverse Momentum Distribution of DIS 1-jet Inclusive Cross Section

### 3.1 Theoretical and experimental considerations

To probe the differences between different treatments of NLO corrections in an event generator, we need an infrared sensitive cross section that potentially gives us the maximum difference between the different methods. We choose the rapidity and transverse momentum distribution of the 1-jet inclusive cross section, $`d^2\sigma /(d\widehat{\eta }dx_t)`$ in the laboratory frame, with the rapidity and transverse momentum measured relative to the simplest parton-model values. At the LO parton level, the scattered quark has a transverse momentum of $`E_t^{\mathrm{PM}}`$ and pseudorapidity $`\eta ^{\mathrm{PM}}`$, both of which can be calculated from the momentum of the scattered electron, and so we define

$`x_t`$ $``$ $`{\displaystyle \frac{E_t^{\mathrm{jet}}}{E_t^{\mathrm{PM}}}},`$ (1)
$`\widehat{\eta }`$ $``$ $`\eta ^{\mathrm{jet}}\eta ^{\mathrm{PM}}.`$ (2)
If we ignore parton virtuality and transverse momentum, the LO cross section is a simple $`\delta `$-function

$$\frac{d^2\sigma }{d\widehat{\eta }dx_t}\delta (x_t1)\delta (\widehat{\eta }).$$
(3)
In real QCD, with effects from higher order QCD subprocesses and the hadronization, the $`\delta `$-function will be smeared out. The shape of the smeared $`\delta `$-function predicted by a MC event generator will provide us information about its treatment of the NLO corrections.

We calculate the cross section from the event generator RAPGAP by three methods:

1. LO hard scattering, where jet structure is generated exclusively by parton showering and hadronization.
2. With NLO corrections implemented by the $`k_{t\mathrm{cut}}`$ method.
3. With NLO corrections implemented by our subtraction method.

Due to the enhanced gluon density at small $`x_{\mathrm{bj}}`$, the cross section in the small $`x_{\mathrm{bj}}`$ region is likely to give us the maximum difference between the different methods. At large $`x_{\mathrm{bj}}`$ the NLO corrections are suppressed by $`\alpha _s`$, and therefore in this region the cross section is not sensitive to the difference between the NLO methods. So in this paper we calculate this cross section in two different regions, $`10^4<x_{\mathrm{bj}}<10^2`$, where we expect to see very different results from the three calculations, and $`0.02<x_{\mathrm{bj}}<1`$, where similar results are expected.

Our calculations are performed under the following conditions, with hadron kinematics appropriate to the collider experiments at HERA:

* The positron is moving in the $`z`$ direction with $`P_z^e=27.5`$ GeV.
* The proton is moving in the $`+z`$ direction with $`P_z^p=820`$ GeV.
* $`20<Q^2<50`$ GeV2.
* In LO and $`k_t`$-cut calculations, the LO pdfs Cteq4L are used.
* In RAPGAP with our subtraction method, we use the BS-algorithm specific pdf calculated by a scheme change  from the $`\overline{\mathrm{MS}}`$ pdf, Cteq4M.
* The jet-finder is PXCONE, with $`E_t`$-mode and cone radius$`=1`$.
* The running scale of the pdfs and $`\alpha _s`$ is $`Q`$.
* Initial state showering is turned on, with $`Q^2`$ ordering.
* Final state showering is on.
* Hadronization and proton remnant are required.
* The LO subprocess and (for the NLO calculations) the gluon-induced NLO subprocesses are chosen.
* Only light quarks are considered.

Figs. 1, 2 and 3 give the results for $`10^4<x_{\mathrm{bj}}<10^2`$, while Figs. 4 and 5 give the results for $`0.02<x_{\mathrm{bj}}<1`$.

### 3.2 Cross Section at small $`x_{\mathrm{bj}}`$

In Fig. 1 we plot $`d^2\sigma /(d\widehat{\eta }dx_t)`$ as a function of $`x_t`$ and $`\widehat{\eta }`$ for the small $`x_{\mathrm{bj}}`$ region, $`10^4<x_{\mathrm{bj}}<10^2`$. It clearly shows that the LO (graph c) prediction is strongly peaked at the PM value $`x_t=1`$ and $`\widehat{\eta }=0`$, as expected. The cross section predicted by the $`k_t`$-cut method (graph b) is somewhat broader but still peaked at the PM values. In contrast, the prediction from our subtraction method (graph a) is broadly distributed; although it has a maximum close to the PM position, the peak is not at all $`\delta `$-function like. This indicates that when the gluonic NLO term is comparable to the LO term, a correct treatment of off-shell parton kinematics has a substantial effect on the shape of the cross section.

For a more quantitative comparison between these methods, we now concentrate on the one parameter cross sections $`d\sigma /d\widehat{\eta }`$ and $`d\sigma /dx_t`$, shown in Figs. 2 and 3. In each figure, we show in graphs (a), (b) and (c) the components of the 1-jet inclusive cross section separated according to the number of jets in the final state: 1 jet, 2 jets and $`3`$ jets. Then in graph (d), we give the total 1-jet inclusive cross section. The solid curves, dashed curves and the dotted curves are the predictions of our subtraction method, the $`k_t`$-cut method, and the LO calculation, respectively.

The differences between the predictions are evident. The subtraction method gives the broadest distribution, while the $`k_t`$-cut method and the LO results show progressively more strongly peaked distributions. A comparison with data should fairly easily test our claim that the subtraction method is more correct. We also see that the number of 1-jet events and the number of 2-jets events predicted by our subtraction method are of the same order, while the $`k_t`$-cut method and LO predict more dominance of 1-jet events222The number of events is roughly proportional to the area under the curve divided by the number of jets in an event.. Also note that, in Fig. 2, the peaks of the curves are not at the PM position $`\widehat{\eta }=0`$. Even though the peak corresponds to a smeared out $`\delta (\widehat{\eta })`$, the peak has shifted to $`\widehat{\eta }`$ around $`1/2`$. The reason is that showering and hadronization of the jet changes a massless on-shell quark into a massive quark, and the kinematics of 2-body collisions shows that this moves the rapidity of the quark towards the proton’s rapidity. Since the proton is chosen to move in the $`+z`$ direction, this means that the rapidity of the jet is biased to a positive value relative to the simple parton model value.

### 3.3 Cross Section at large $`x_{\mathrm{bj}}`$

In Figs. 4 and 5 we show the results for large $`x_{\mathrm{bj}}`$. The parameters are the same as before, except that we choose to analyze events with $`0.02<x_{\mathrm{bj}}<1`$. In contrast to the small $`x_{\mathrm{bj}}`$ results, the three different methods give quite similar results, all peaked near $`\widehat{\eta }=0`$ and $`x_t=1`$. In addition, most of the events are 1-jet events. This is because gluon distribution is now much smaller, and the LO process dominates.

## 4 Applications and Future Work

We have implemented our subtraction method  in RAPGAP for DIS processes. The subtraction method is designed to give perturbative QCD predictions for the hard scattering with a correct treatment of parton kinematics. We chose to calculate an infra-red sensitive observable $`d^2\sigma /dx_td\widehat{\eta }`$ using three different methods, and to perform the calculation in a region where the NLO correction is large because of a large gluon density. We see large differences between the subtraction method and the $`k_t`$-cut method, thereby demonstrating the importance of the treatment of parton kinematics.

There is no experimental analysis on $`d^2\sigma /dx_td\widehat{\eta }`$ available yet, current analyses  being in the Breit frame. Since our subtraction method, combined with the parton shower, is intended to give more accurate perturbative QCD predictions for infrared sensitive cross sections, it would be interesting to do the experimental analysis of the jet cross section in the laboratory frame relative to the position of parton-model jet, and to compare it with the MC calculations using our subtraction method.

In the future we would like to work on two directions. One is to implement the subtraction method for diffractive DIS where there is also a large gluon distribution. Another is to work on the quark-induced NLO subprocess in order to get a complete NLO MC event generator.

###### Acknowledgments.

We would like to thank Professor Hannes Jung for his generous help, various discussion and helpful suggestions. This work was supported in part by the U.S. Department of Energy under grant number DE-FG02-90ER-40577.