Brownian motion in short range random potentials

Department of Physics, University of California, San Diego, La Jolla, California 92093-0354 Department of Chemistry and Biochemistry, University of California, San Diego, La Jolla, California 92093-0340 Institute for Nonlinear Science, Department 0407, University of California, San Diego, 9500 Gilman Drive, La Jolla, California 92093-0407 and Departament d’Estructura i Constituents de la Mate ́ria, Avenida Diagonal 647, 08028 Barcelona, Spain ~Received 9 March 1998 !


I. INTRODUCTION
Since the pioneering work of Einstein ͓1͔, it has been well known that the motion of a free Brownian particle is random and isotropic.It is of diffusive nature, with a bare diffusion coefficient , determined by the mean square displacement where d is the dimension of the space.The diffusion coefficient is proportional to the absolute temperature of the system, indicating which is the origin of this type of random motion.The brackets indicate statistical ensemble average of many independent particles.The distribution of particles spreads homogeneously in all directions, as in a normal diffusion process.A very simple way to describe this motion is to use a Langevin equation in the Schmoluchowsky limit.
In this paper we consider the motion of the same Brownian particles in a random potential formed of wells and hills whose location and magnitude are random quantities.In this case, particles do not move randomly in all directions as in the former example, but they follow, mostly, the easier paths connecting wells, remaining long time intervals in the deepest wells.The mechanism of the motion is barrier crossing.So, when we consider an ensemble of these noninteracting particles, instead of a homogeneous distribution spreading through all space, what we find is a set of localized regions ͑wells͒ where particles spend most of the time.
For very large temperature values (), with respect to the relative magnitude of the potential wells, one can define an effective diffusion coefficient such as in Refs.͓2,3͔.Nevertheless, a consistent analysis of the numerical trajectories seems to indicate that the regime is subdiffusive with a better description in terms of effective exponents, ͗⌬r 2 ͑ t ͒͘ϳt z eff .

͑2͒
Diffusion in random static potentials has been a subject of a rich variety of studies, most of them using master equations ͓4-6͔.Very powerful theoretical tools have been applied for this problem where Langevin equations have been used as the basic model ͓6-9͔, but only a few numerical simulations have been reported up to now ͓6͔.Thus there is very little information on how the theoretical results can be compared with real simulations.
In this work, we present an extended numerical study of the random diffusion of independent particles embedded in a random correlated potential.A key point, which is the correct statistical generation of the potential, is also discussed.The roles of the potential parameters, intensity and characteristic length are studied.Numerical simulations for one and two dimensions are performed.
In Sec.II, the model and noise characteristics are presented.In Sec.III, we explain the numerical algorithms used and the simulation set up. Results are presented and discussed in Sec.IV.

II. MODEL
Our starting point is a Langevin equation for Brownian particles in the high friction limit ͑the Smoluchouski limit͒ x ˙ϭϪVЈ͑x͒ϩ͑t͒, ͑3͒ where the Cartesian components of the Gaussian white noise () have zero mean value and correlation

͑4͒
The potential we will consider here V(x) is a Gaussian short ranged correlated variable, with zero mean value and correlation ͗V͑x͒V͑xЈ͒͘ϭg͑xϪxЈ͒ d .

͑5͒
This is a nonidealized potential such as white correlated noise, and it can represent more realistic situations.The potential is described by only two parameters: its effective intensity and its correlation length , In this way, for →0, it is possible to recover the white noise limit for the same intensity ⑀.
In what follows we will assume g to be of the form This is a short range potential with characteristic length and strength controlled by the parameter ⑀.
In Refs.͓6,8,9͔ a different prescription was taken.Instead of the potential V(x), they took the force F(x).The equivalence between both descriptions is easily seen, Nevertheless, no correlation length appears explicitly in their approach.We will comment on this point below.

III. NUMERICAL ALGORITHM AND POTENTIAL GENERATION
Now we will review the main steps in the numerical implementation of the problem we are studying.The equation of motion ͑3͒ is numerically solved by means of a second order predictor-corrector algorithm ͓10͔.The motion will take place in a d-dimensional space discretized in N d cells of size ⌬x d .The linear size of the system is LϭN⌬x.
Initially, 250 particles are uniformly randomly located in each one of four domains of linear size L/4, covering the whole system.The simulation is stopped when a first particle runs a distance of the order of L/4.
Here we present details for the one-dimensional case ͑see Fig. 2͒.The two-dimensional case is easily generalized ͑see Fig. 1͒.The position of each particle for each time step of integration ⌬ is obtained as where F(x)ϭϪdV(x)/dx is the force, and 's are Gaussian random numbers of zero mean and variance (2⌬) 1/2 , constructed with standard Gaussian random number generators.
The generation of a random potential with prescribed characteristics is not as standard as the generation of white noises such .We assume that the potential is defined in the corners of each cell and we have N d values of it, with periodic boundary conditions ͓11͔.The potential is constructed in the corresponding Fourier space V(k), V͑k͒ϭ"g͑k͒ d … 1/2 ͑k͒, ͑12͒ where g(k) d is the Fourier transform of the correlation function ͑5͒, and (k) are Gaussian random numbers of zero mean and correlation, ͗͑k͒͑kЈ͒͘ϭ␦ kϩk Ј ,0 .

͑13͒
In order to obtain this type of anticorrelation, variables are constructed in the following way: where are Gaussian random numbers of zero mean and variance equal to 1, and the subindexes i and r indicate the real and imaginary parts of .The Fourier inverse transform of V(k) will give the potential values at the corner cells.This procedure guarantees the statistical properties ͓12͔ of definition ͑5͒.
The force values in the corner cells are calculated from the discrete symmetric derivative of V(x), Comparison of the correlation function, g(r), from the predicted theoretical form and our simulation for the one-and twodimensional cases 1D ͑a͒ and 2D ͑b͒.However, during the numerical integration of Eq. ͑3͒, we need to evaluate the force values for any point inside the simulation region.These values are generated by linear interpolation from the last obtained discrete forces as This numerical scheme can be improved, but as long as ⌬x Ӷ no important errors are made.

IV. NUMERICAL RESULTS AND THEORETICAL PREDICTIONS
For the stochastic dynamics of an overdamped Brownian particle in a random potential, in our numerical simulations, we have observed several time regimes.In the early stages of the evolution, one can identify two different regimes as function of the noise: a deterministic regime for very small noise intensities, and a bare diffusion regime for large noise.In later stages, one can observe two other regimes: either a subdiffusive regime for noise intensities that are not too small, or a frozen regime if noise is very small.These different dynamical regimes can be clearly seen with an appropriate selection of the parameters.Here we will review the theoretical and numerical aspects of the first three regimes.The last one is commented upon later in this section.

A. Deterministic short time regime
For a very small noise intensity ͓a small bare diffusion compared with "g(0) d … 1/2 ͔ particles behave deterministically, relaxing from their initial position (x 0 ) to the bottom of the nearest well of the stochastic potential.This behavior can be understood if every particle is assumed to relax to a potential well approximated by a parabola.The parameters defining the parabolic geometry are taken from the potential average height ͓␣ϭg(0) d 1/2 ͔ and potential average width (2): The solution of the deterministic equation of motion, when (t)ϭ0 in Eq. ͑3͒, shows that the mean square displacement, over the initial density (0), is ͗⌬x͑t͒ 2 ͘ϭ͗x 0 2 ͑͘e Ϫ␣t/ 2 Ϫ1͒ 2 .

͑20͒
In Fig. 3, our numerical simulations show these two regimes.In particular, the early regime described by Eq. ͑19͒ is seen in the trajectories label by b and bЈ of this figure.Both cases have the same value for ⑀/ 3 , and therefore they exhibit the same short time behavior.Also, the t 2 power law is also clear in these trajectories.
The long time regime of Eq. ͑20͒ is also present in the trajectories labeled by b and bЈ in Fig. 3.It is clear that the final frozen mean square displacement values are controlled only by .In particular, with the data used, the final value of ͗⌬r 2 (t)͘ for case b (ϭ8) is four times the one correspond- ing to case bЈ (ϭ4).

B. Bare diffusion regime
Here we consider the case of the short time regime for a noise intensity that is not very small.Within the same approximation used in the deterministic regime, a parabolic potential, particles relax to the bottom of the well, but the influence of the noise is manifested and the particles perform a Brownian motion.The equation of motion with the assumed potential ͑17͒ is now x ˙ϭϪ ␣x The mean square displacement, over the noise realizations and the initial density distribution, is then

͑22͒
and two more limiting behaviors, can be appreciated.͑a͒ For times larger than ␣t/ 2 , we find that ͗⌬x͑t͒ 2 ͘ϳ2t.

͑23͒
͑b͒ For very large times, we recover the same saturation values as in Eq. ͑20͒.
So, independently of the potential parameters (⑀,), for neither small nor large times there exists a diffusive regime with the bare diffusion .
The regime described by Eq. ͑23͒ is seen in cases labeled by a and aЈ in Fig. 3.We see how the early stages are clearly diffusive with the bare diffusion parameter .The theoretical result ͑22͒ saturates as in Eq. ͑20͒, but this regime is not seen now because a new physical mechanism starts, which is barrier crossing.The subdiffusive character of the new regime is manifested for large times of the trajectories in Fig. 3.

C. Subdiffusive regime
This is the most important regime in our study.When the noise intensity is not too small ͓g(0) d / 2 Ͻ40.0͔, the system clearly exhibits, in our numerical simulations, a subdiffusive regime.It starts when the bare diffusion regime saturates as seen in Fig. 3.In other words, it starts when particles are able to overpass the barriers limiting the well in which they were initially located.In this regime, depending on the noise intensity, particles can visit several wells by a barrier crossing process.We expect that this regime is not of a diffusive nature with effective exponents depending on the system parameters: the bare diffusion coefficient (), the potential strength (⑀), and its characteristic length ().With these three parameters, the only possible adimensional combination is ͑24͒ This quantity has to be considered as the most relevant parameter for this problem.Numerical simulations have been performed, in one and two dimensions, for different values of the system parameters.The effective exponents have been extracted from the trajectories of the mean square displacement by numerical fitting of a power law AϩBt z eff .The set of exponents versus the adimensional quantity ͓Eq.͑24͔͒ is presented in Fig. 4.
According to renormalization group calculations ͓6,8,9͔, a linear dependence of the type is expected for dϭ2.Numerical fitting in the dϭ2 case gives the following values for the free parameters: aϭ1.01 and bϭ0.023.The value obtained for b is almost two times smaller than the theoretical calculated value bϭ1/8 ͓6,8͔.
Equation ͑25͒ looks very similar to the theoretical prediction expression, with (0) instead of g(0) ͓see Eq. ͑9͔͒; however the theoretical result does not make any reference to the characteristic length ͓6,8͔.This could explain the difference between the numerical value and the one obtained from the theory.The characteristic length should be explicit in any short range force or potential.
The data dispersion in Fig. 4 gives an idea of the error bars.In general, the errors are lower than 5% in dϭ2, and all our data correspond to a unique potential.So our statistical average is over an ensemble of ''particles'' ͑1000͒ and not over an ensemble of different ''potentials.''If the system size is large enough, with many wells and hills, a second average over different potentials is not a crucial point.Nevertheless, for very small values of the intensity of the noise, particles ''see'' very few wells and the statistic is poorer.
In tion are shown with high and low noise intensities ͓compared with g(0) 2 ͔.In Fig. 5͑a͒, particles spread out, covering a considerable space distance from the initial location as in the homogeneous diffusion case.In Fig. 5͑b͒, this quasihomogeneity is lost, and particles mimic the underlining potential pattern.In this last case, many more statistics would be needed to obtain reliable values for the effective exponents due to the small number of wells explored.In any case, the general trend of the nondiffusive behavior is clear in our simulations.
Simulation results for dϭ1, plotted for reference in Fig. 4, deserve an explicit discussion.A linear fitting, like in Eq. ͑25͒, for parameters a and b gives the values 1.01 and 0.038, respectively.Moreover, theoretical studies ͓6͔ predict that, for this dimensionality, no expression like Eq. ͑25͒ is to be expected.Instead of the power law with an effective exponent, a logarithmic power law ͗⌬x 2 ͘ϳ(lnt) 4 has been proved.This last behavior is not clear in our simulations, due to the special restrictive characteristic of the motion.Indeed, in dϭ1, a very deep well would act as a trap or reflecting wall, bounding the motion.On the other hand, for dϭ2, the particle would have more possible ways to escape from this trapping case.This contrasts with the dϭ1 case, where larger simulations would be needed to discriminate the regime.
Summarizing the main results of this work, we have done numerical simulations of Brownian particles under a short ranged potential.The earlier stage of the motion has been satisfactorily explained as a function of the system parameters.The later subdiffusive character of the motion is clearly seen in two dimensions.The effective exponents obtained agree with existing theoretical predictions, provided the different system parameters ͑noise intensity, potential strength, and characteristic length͒ are properly taken into account.
With respect to the one-dimensional case, numerical simulations do not match theoretical predictions very well.To see the effects predicted by the theory, very long time simulations in larger systems and ensemble potential averages would be required.

FIG. 2 .
FIG. 2. Potential in one dimension ͑dashed line͒ and the density of particles after some time.(0) is the uniform initial particle density.

Fig. 5 ,
FIG. 4. Effective exponents vs the adimensional parameter.Solid symbols refer to the 2D case, and open symbols to the 1D case.