# Generation of short and long range temporal correlated noises

## I Introduction

In stochastic process simulations, the generation of noise with some specific statistical properties has been of great importance. In most cases, the required noise is Gaussian and white (delta correlated); in many other cases, as in noises of real systems, a specific correlation is needed and an appropriate noise generation has to be implemented which can be carried out in standard computer facilities. Despite the extensive work in developing algorithms to generate noises with particular correlations, there is still a lag concerning those noises with any given temporal, spatial or spatiotemporal correlation.

Some algorithms have been proposed in the last few years for noises which have been proved to obey a linear Langevin equation with a linear, Gaussian, white, noise term. A formal integration of this linear equation is enough to generate a random process with a very particular time correlation. Examples of this approach are the so called Ornstein-Uhlenbeck (O-U) and Wiener (W) processes. The physical meaning of the O-U process is the velocity of a Brownian particle under the influence of friction and immersed in a heath bath. W-process is the archetype for the dynamical behavior of change in position of a free Brownian particle in the high friction limit. The correlations of these two processes are well known and can be obtained analytically  and numerically  by standard methods of stochastic processes.

Nevertheless, quite often in numerical stochastic simulations, we are faced with processes characterized by a Gaussian noise with a specific correlation and an unknown Langevin-like equation dynamics. Therefore, a more general algorithm that only depends on the knowledge of the temporal correlation is necessary .

The purpose of this work is to present the implementation of an algorithm which allows to simulate Gaussian noises with almost any given temporal correlation function. The only requirement for the algorithm to work is that the Fourier transform of the temporal correlation function should be known. In cases where the Fourier transform of the required temporal correlation function can not be properly defined, the algorithm can be used with a suitable cutoff.

In Section II, algorithms based in Langevin equations are revised and the new algorithm of noise generation is presented in practical terms. The algorithm has the property of simulating different correlation regimes, from short to long ranged. Section III considers different applications for which the dynamics of the system are very sensitive to the noise correlation. A summary of results and some conclusions are presented in the last section.

## II Algorithms to Generate Gaussian Noises

### A The Linear Langevin Equation Method

For the sake of comparison with the method we want to present here, let us revise briefly the main points of the generation of a Gaussian noise which follows a linear Langevin equation in terms of a Gaussian white noise. The stochastic discretized trajectories are generated by formally integrating the corresponding Langevin equation and using a set of Gaussian random numbers for a given realization. The statistical properties are calculated from many realizations of these trajectories (“an ensemble”). In the case of Gaussian noise, only two moments are necessary to characterize the statistical properties of the random process. Due to the linear character of the Langevin equation the Gaussian property of the white noise is transmitted to the generated noise.

The simplest case is the W-process, which follows the Langevin equation

$$\dot{\eta }=\xi (t),$$
(1)   
where $`\xi (t)`$ is a Gaussian white noise of zero mean and delta correlation:

$$\xi (t)\xi (t^{})=2ϵ\delta (tt^{}).$$
(2)   
Here, $`ϵ`$ is the noise intensity.

The statistical properties are easily evaluated  and a trajectory for this type of noise can be obtained by formally integrating Eq.(1) over time:

$$\eta (t+\mathrm{\Delta }t)=\eta (t)+_t^{t+\mathrm{\Delta }t}\xi (t^{})𝑑t^{}=\eta (t)+X(t),$$
(3)   
where the noisy term $`X(t)`$, is constructed from

$$X(t)=\sqrt{2ϵ\mathrm{\Delta }t}\alpha (t).$$
(4)   
Here, $`\alpha (t)`$ is a set of Gaussian independent random numbers of zero mean and variance equal to unity, obtained from any reliable Gaussian random generator .

As a second example, we consider the O-U process. This process contains a temporal memory and obeys the Langevin equation

$$\dot{\eta }=\frac{\eta }{\tau }+\frac{\xi (t)}{\tau },$$
(5)   
where $`\tau `$ is the characteristic time memory. As in the former case, a formal integration gives:

$$\eta (t+\mathrm{\Delta }t)=\eta (t)e^{\frac{\mathrm{\Delta }t}{\tau }}+\frac{ϵ}{\tau }_t^{t+\mathrm{\Delta }t}e^{\frac{tt^{}}{\tau }}\xi (t^{})𝑑t^{}.$$
(6)   
Studying the statistical properties of the noisy term, the algorithm reads,

$$\eta (t+\mathrm{\Delta }t)=\eta (t)e^{\frac{\mathrm{\Delta }t}{\tau }}+\sqrt{\frac{ϵ}{\tau }(1e^{\frac{2\mathrm{\Delta }t}{\tau }})}\alpha (t).$$
(7)   
In this way, a time stochastic trajectory is generated step by step as in the former example. Moreover, in the algorithm we introduce below, the whole trajectory is constructed in a single calculation step.

### B Spectral method

As we noted in the introduction, this method starts from the knowledge of the time correlation function. Probably the basis of this algorithm has been rediscovered and implemented several times, but now, due to the actual high speed and large storage space in computers, the algorithm can be implemented rather easily. Here we follow the main ideas of Ref.  and first introduce some steps that simplify the calculations.

We want to generate a Gaussian, random correlated, noise, $`\eta (t)`$, whose correlation function $`\gamma (t)`$, is defined by:

$$\eta (t)\eta (t^{})=\gamma (|tt^{}|),$$
(8)   
and its Fourier transform,

$$\gamma (\omega )=e^{i\omega t}\gamma (t)𝑑t.$$
(9)   
is known (to some extent). In the $`\omega `$-Fourier space, this correlation reads:

$$\eta (\omega )\eta (\omega ^{})=2\pi \delta (\omega +\omega ^{})\gamma (\omega ),$$
(10)   
where $`\eta (\omega )`$ is the Fourier transform of $`\eta (t)`$.

With this initial information in mind, the algorithm could be summarized as follows. First we discretize the time in $`N=2^n`$ intervals of mesh size $`\mathrm{\Delta }t`$, and note that this time interval has be much smaller than any other characteristic time of the system for our method to work. Every one of these intervals will be denoted by a Roman-index in real space (time) and by a Greek-index in Fourier space (frequency).

In the discrete Fourier space, the noise has a correlation given by:

$$\eta (\omega _\mu )\eta (\omega _\mu ^{}^{})=N\mathrm{\Delta }t\delta _{\mu +\mu ^{},0}\gamma (\omega _\mu ).$$
(11)   
Now, $`\eta (\omega _\mu )`$ can be constructed from

$`\eta (\omega _\mu )=\sqrt{N\mathrm{\Delta }t\gamma (\omega _\mu )}\alpha (\omega _\mu ),`$ (12)   
$`\mu =0\mathrm{}N,\eta (\omega _0)=\eta (\omega _N),\omega _\mu ={\displaystyle \frac{2\pi \mu }{N\mathrm{\Delta }t}},`$ (13)   
where $`\alpha (\omega _\mu )\alpha _\mu `$ are Gaussian random numbers with zero mean and correlation,

$$\alpha _\mu \alpha _\nu =\delta _{\mu +\nu ,0}.$$
(14)   
This type of delta-anticorrelated noise can be generated rather easily if the symmetry properties of real periodic series in the Fourier space ($`\alpha _i`$) are used . Thus, we avoid the extra work involved in Fourier transforming real random numbers. Consider for example a system of size $`N`$, the Fourier components of a periodic series are then related by:

$$\alpha _\mu =\alpha _{\mu +pN}\alpha _\mu =\alpha _\mu ^{},$$
(15)   
where $`p`$ is an integer number. Since $`\alpha _{\mu =0}`$is real and $`Im(\alpha _{\mu >0})=Im(\alpha _{\mu <0})`$, as can be seen in figure 1, the requested anticorrelated complex random numbers (14) can be constructed as

$$\alpha _\mu =a_\mu +ib_\mu ,b_0=0,$$
(16)   
where $`a_\mu `$ and $`b_\mu `$ represent Gaussian random numbers with zero mean and a variance of one half:

$$a_\mu ^2=b_\mu ^2=\frac{1}{2},\mu 0,a_0^2=1.$$
(17)   
The discrete inverse transform of $`\eta (\omega _\mu )`$ is then numerically calculated by a Fast Fourier Transform algorithm . The result is a string of N numbers, $`\eta (t_i)`$, which by construction, have the proposed time correlation (8). However, due to the symmetries of the Fourier Transform, only $`N/2`$ of these values are actually independent and the remaining numbers are periodically correlated with them.

In order to check the suitability of the procedure, the time correlation of Equation (8) is numerically evaluated by an independent (non Fourier based) method, namely,

$$\gamma (t_i)=\frac{_{j=0}^{N_{max}}\eta (t_j+i\mathrm{\Delta }t)\eta (t_j)}{N_{max}+1}.$$
(18)   
Here, $`N_{max}`$ is a number smaller than $`N/2`$ (in the examples that follow, it is taken equal to $`N/4`$).

We now present several examples where this approach is implemented. Since the number of applications making use of temporal random noises is quite large we have chosen applications with rather different correlation properties in order to illustrate the power of the method. The same procedure has already been used for spatial correlated noises .

### C Short range correlated noises

a.1. A Gaussian Correlation.

Let us consider a noise with a Gaussian correlation function defined as

$$\eta (t)\eta (t^{})=\gamma (|tt^{}|)=\frac{2ϵ}{\tau \sqrt{2\pi }}e^{\frac{|tt^{}|^2}{2\tau ^2}},$$
(19)   
where $`ϵ`$ and $`\tau `$ are the noise intensity and correlation time, respectively. The correlation is normalized in such a way that

$$ϵ=_0^{\mathrm{}}\gamma (t)𝑑t.$$
(20)   
Setting $`\tau 0`$, the white noise limit is recovered. The Fourier transform of this Gaussian correlation(19) is given by:

$$\gamma (\omega )=2ϵe^{\frac{\tau ^2\omega ^2}{4}}.$$
(21)   
According to the prescription (13), we now generate the discrete field $`\eta (\omega _\mu )`$ as:

$$\eta (\omega _\mu )=\left(N\mathrm{\Delta }t2ϵe^{\frac{\tau ^2}{\mathrm{\Delta }t^2}(cos(2\pi \mu /N)1)}\right)^{1/2}\alpha _\mu .$$
(22)   
In Fig 2, an explicit comparison between the numerical results and the expected theoretical prediction is presented.

a.2. The Ornstein-Uhlenbeck process.

The Ornstein-Uhlenbeck process simulates the behavior of the velocity of a Brownian particle under friction and immersed in a thermal bath. Quite often, it is used to represent a real noise with memory, whose intensity is $`ϵ`$ and whose correlation time (or memory intensity) is $`\tau `$. This is a well known Gaussian, Markovian, and stationary noise, which obeys a linear Langevin equation (5). Its well known correlation is given by:

$$\eta (t)\eta (t^{})=\gamma (|tt^{}|)=\frac{ϵ}{\tau }e^{\frac{|tt^{}|}{\tau }},$$
(23)   
with the same normalization used for the previous example. Since this noise obeys a linear Langevin equation, which was not the case for the first example, a realization of this noise could be simulated using the formal solution of the stochastic differential equation given by Eq.(7). Instead, we present in this subsection the results obtained following our spectral method. The Fourier transform of this particular correlation function is:

$$\gamma (\omega )=\frac{2ϵ}{1+(\tau \omega )^2}.$$
(24)   
According to the prescription (13), we generate a discrete field $`\eta (\omega _\mu )`$ as:

$$\eta (\omega _\mu )=\left(\frac{N\mathrm{\Delta }t2ϵ}{1+(\frac{2\tau }{\mathrm{\Delta }t}sin(\pi \mu /N))^2}\right)^{1/2}\alpha _\mu .$$
(25)   
In Fig. 2, we compare the numerical and theoretical results. The discretization of the $`\omega `$–variable using either the function $`\mathrm{cos}(\frac{2\pi \mu }{N})`$ as in Eq.(22) or the function $`\mathrm{sin}(\frac{\pi \mu }{N})`$ as in Eq. (25) does not make any difference.

We have seen that in the short range correlated noises considered here, the agreement between the statistical properties of the noise generated from the spectral method and the statistical properties required from the noise are quite good.

### D Long range correlated noises

A more dificult noise, in the sense that it can not be obtained from any known linear Langevin equation, is the one characterized by a power-law decaying correlation function,

$$\gamma (|tt^{}|)\frac{ϵ}{|tt^{}|^\beta },0<\beta <1.$$
(26)   
This correlation is not well defined in Fourier space (it has a singularity at $`|tt^{}|=0`$). Therefore, in order to implement the spectral method, we start with a guess for $`\gamma (\omega _\mu )`$, and then we look at the dynamics of $`\gamma (t)`$ in real space. The Fourier values of the noise are now discretely generated according to the expression,

$$\eta (\omega _\mu )=\frac{N\mathrm{\Delta }tϵ}{\left[\frac{2\tau }{\mathrm{\Delta }t}\mathrm{sin}(\pi \mu /N)+\omega _0\right]^{(\beta 1)/2}}\alpha _\mu ,$$
(27)   
where $`\omega _0`$ is a predefined cut-off. In Ref. modified Bessel functions are used for the correlation in Fourier space but they have the same long range decay. In our case, we assume the following form for the correlation function in real space,

$$\gamma (t)=\frac{Aϵ}{\pi \beta (t+t_0)^\beta },0<\beta <1.$$
(28)   
where $`t_0=\mathrm{\Delta }t/\pi `$, $`\beta `$ is the parameter describing the power law decay, $`ϵ`$ is the noise intensity, and $`A`$ is a parameter to be fitted from a correlation average.

Fig. 3, presents two examples of the power-law Gaussian noise with $`\beta =1/3`$ and $`\beta =2/3`$. It can be seen that the numerical results are fitted very well by the expected power laws.

The existence of power law noises has been predicted and discussed for quite some time in the literature. Many examples arise from theoretical as well as from experimental studies. A recent reported case is the examination of the temperature fluctuations in climatological data , where a power law functional form is found for the correlation of these fluctuations. In this particular case, the analysis is carried over a temperature autocorrelation function as defined by Eq. ( 18), and it is found that the distribution of temperature fluctuations is well described by a Gaussian distribution with the long-ranged decaying auto-correlation function: $`C(s)s^\gamma `$. The exponent found for their data is around $`2/3`$. Noises with these properties can be generated using our procedure and can be implemented in modelling processes similar to the one reported in this climatological anylisis.

Fig. 4 shows a realization of a particular noise trajectory generated with our method for a power law noise with $`\beta =2/3`$. The variations of the data generated with the algorithm look very similar to figures 1.b, 6.a and 6.b in ref. . The persistence of the series is evident from our simulation where the stochastic trajectory appears as “packets”.

## III Applications

We proceed now to present and discuss two different non-trivial examples in non-equilibrium statistical physics where our algorithm can be applied: (i) the dispersion process of a Brownian particle and (ii) the decay from an unstable state, both cases under the influence of long range noises generated through the spectral method. In the first case the noise is additive whereas in the second case it appears multiplicatively.

### A Superdiffusive motion

In this section we want to focus our attention on the random motion of a Brownian particle which obeys the Langevin equation

$$\frac{dx}{dt}=\eta (t),$$
(29)   
where $`\eta `$ is a Gaussian noise with a given correlation function. We consider three different cases, a short range and two long range noises.

The solution for the relative dispersion is well known:

$$\delta x(t)^2=2_0^t_0^t^{}\gamma (s)𝑑s𝑑t^{}.$$
(30)   
For large times, this expression is either $`t`$ (diffusion) for short range noises, or $`t^{2\beta }`$ (super–diffusion) for long range noises. The explicit form of the time dependence of this quantity appears in Fig. 5, where we also show the variance of the three different examples in which our algorithm is applied. Two super–diffusive cases are clearly seen, with dispersion exponents $`1.75`$ and $`1.25`$, corresponding to noise decay power laws, $`\beta =0.25`$ or $`0.75`$, respectively. A short range noise (exponential) is also included for comparison. We see that in this last case the behavior is ballistic (deterministic) $`t^2`$ for $`t<\tau `$, but diffusive for $`t\tau `$ as predicted by Eq. (30).

### B Decay of an unstable state

The decay of an unstable state is one of many interesting problems which appear in nonequilibrium phenomena and nonlinear relaxation process. The switch-on process of a dye laser has been a prototype of a nonequilibrium situation in which the influence of several sources of noise have been tested. The dynamics of this system can discriminate the effects of both additive or multiplicative noises . It has been established that additive noise is responsible for the short time dynamics but multiplicative noise effects appear in the medium and long term dynamics. In previous studies, white or short range, colored noises have been used, but long range noises have never been considered. Here we will present several numerical simulation results for the decay of an unstable state under the influence of long range multiplicative noises.

In order to simplify the analysis, let us illustrate the situation with the following Langevin equation,

$$\frac{dx}{dt}=axbx^3+x\eta (t),$$
(31)   
where $`\eta `$ is a Gaussian noise whose correlation function is to be specified. For $`a`$ and $`b`$ positive parameters, the initial value $`x(0)=0`$ is an unstable state which does need small perturbation to start relaxing towards its steady state. To trigger this process we need either the additive noise of Eq. (31) or an initial distribution for $`x(0)`$. We have chosen this last assumption and the initial values of $`x(0)`$ are Gaussian distributed with statistical moments $`x(0)=0`$ and $`x(0)^2=\sigma ^2`$. Due to the symmetries of the problem, the mean value is $`x(t)=0`$, and hence we look at the dynamical evolution of the second moment.

To get a precise idea of the short time behavior, it is enough to study Eq. (31) in a linear approximation. Formal integration of this equation gives,

$$x(t)=x(0)exp\left(at+_0^t\eta (t^{})𝑑t^{}\right).$$
(32)   
Using Gaussian properties through the calculation, we obtain for the second moment

$$x^2(t)=\sigma ^2exp\left(2at+4\mathrm{\Omega }(t)\right),$$
(33)   
where $`\mathrm{\Omega }(t)`$ is defined by:

$$\mathrm{\Omega }(t)=_0^t_0^t^{}\gamma (s)𝑑s𝑑t^{}.$$
(34)   
For the deterministic case, $`ϵ=0`$, and the multiplicative white noise case, we have:

$$x^2(t)=\sigma ^2exp\left(2at\right),$$
(35)   
$$x^2(t)=\sigma ^2exp\left(2at+4ϵt\right),$$
(36)   
respectively.

As in the superdiffusion case, we get for long range noises a power time dependence $`t^{2\beta }`$ in the exponential; we therefore expect a larger rate for the relaxation of the initial state.

In Fig. 6 we see how the decay of the unstable initial state is influenced in the different cases discussed. For the values of the parameters we have used, Eqs.(35) and (36) predict that the multiplicative white noise should decay six times faster than that the deterministic case. For long range correlated noises we expect that the smaller the exponent $`\beta `$ is, the faster it will decay. Within statistical errors, these points are actually seen in the simulation data presented in the figure .

It is interesting to note that for intermediate and long times no analytical results can be obtained. However, figure 6 is telling us that the final steady state also depends on the correlation of the noise, taking larger values for smaller $`\beta `$’s. This is not an easy problem and it would need further theoretical analysis.

## IV Summary and Conclusions

We have presented a method to generate Gaussian noises with a prescribed time correlation function which does not depend on any type of dynamics. The algorithm herein incorporates whithin the same framework the generation of long and short ranged noises. In particular, the need for a good algorithm is manifest when long range noises representing physical dynamical process with unknown or very complicated dynamical equations are required. In this sense, this method can numerically simulate the influence of long and short ranged noises in physical systems in a very reliable and controlled way. With the prescription given, many problems with long range realizations can be simplified and the possibility of some kind of analytical approach to specific problems increases (such is the case for the examples we have presented).

The generalization of the algorithm to spaciotemporal noises is straight forward for any number of dimensions. Since by construction, the noise generation is finite, the upper time limit of the simulation has to be selected in advance, but it does not constitute a serious limitation for the algorithm.

## V Acknowledgment

The work was supported in part by the U.S. Department of Energy under Grant No. DE-FG03-86ER13606, and in part by the Comisión Interministerial de Ciencia y Tecnología (Spain) Project No. DGICYT PB96-0241. A.R. thanks K. Lindenberg and A. Sarmiento for support and helpful comments.