- Email: [email protected]
Stochastic continuity equation and related processes
ARTICLE IN PRESS
Physica A 347 (2005) 17–37 www.elsevier.com/locate/physa
Stochastic continuity equation and related processes Gabriele Bassia, Armando Bazzanib, Helmut Maisc, Giorgio Turchettib, a Department of Mathematics and Statistics, University of New Mexico, USA Department of Physics, University of Bologna, INFN Sezione di Bologna, Via Irnerio 46, Bologna 40127, Italy c DESY, Hamburg, Germany
b
Received 13 December 2003; received in revised form 7 July 2004 Available online 27 October 2004
Abstract We consider the processes defined by a Langevin equation and the associated continuity equation. The average of the density function, solution of the continuity equation, satisfies the Fokker–Planck equation. For a volume preserving vector field the same equation is satisfied by the average of the integer powers of the density, which are the moments of the related probability density. For a generic vector field the Fokker–Planck equation for the moments is slightly modified. We first illustrate the problem in the simple case of a free particle subject to a white noise, since the averages can be computed by an elementary procedure using the factorization property of the correlation functions of the noise. The probabilistic meaning of the moments is discussed and the comparison between the analytical results and the numerical simulation is shown. The case of a generic Langevin equation is treated by computing the averages via a Dyson expansion after observing that, for a volume preserving vector field, any power of the density function satisfies the same continuity equation with the appropriate initial conditions. As a consequence the results obtained for the free particle are easily extended. An alternative approach is based on the characteristics of the continuity equation; the probability density of this process in an extended phase
Corresponding author. Tel.: +39 0512091095.
E-mail address: [email protected] (G. Turchetti). 0378-4371/$ - see front matter r 2004 Published by Elsevier B.V. doi:10.1016/j.physa.2004.08.083
ARTICLE IN PRESS 18
G. Bassi et al. / Physica A 347 (2005) 17–37
space still satisfies a Fokker–Planck equation and its moments coincide with the previous definition. r 2004 Published by Elsevier B.V. Keywords: Stochastic processes; Continuity equation; Moments
1. Introduction Stochastic processes are commonly used to build phenomenological models of physical and biological systems, in which a small subsystem is isolated taking into account its interaction with the fluctuating environment. The Langevin equation, continuous version of a random walk, is an ordinary first-order differential equation with a deterministic vector field plus a white noise component. The probability density of the process satisfies the Fokker–Planck equation; its fundamental solution is the weighted sum over all paths connecting the initial and final points in the extended phase space [7]. We consider here the continuity equation associated to the Langevin equation. This is a first-order linear partial differential equation, which defines a stochastic process whose properties are analyzed by considering the average of its integer powers. These are the moments of the probability density associated to the stochastic process, and determine it completely. What we propose is an extension of the stochastic Liouville equation previously considered [1–3], and discussed in the context of beam dynamics [4]. It is important to remark that there are physical processes, such as the Brownian motion, in which each particle experiences a different realization of the noise (describing random collisions with the medium). In this case the physical process is described by the Langevin equation and, even though the stochastic continuity equation for the density r is mathematically well defined, only the average hri; satisfying the Fokker–Planck equation has a physical meaning as the probability density of the process. The evaluation of r at a given location is useful from the numerical viewpoint since it allows to compute the average hri there with a high accuracy. However, there are physical processes in which all the particles experience the same noise: consider for instance a bunch of particles crossing a magnetic element in an accelerator, where the current has small random fluctuations. In this case the density r; satisfying the continuity equation, describes the physical process and hri has to be interpreted as the average over many distinct runs, each one providing a different realization of the stochastic process, for the same initial bunch. In a very high-energy low-intensity beam the electromagnetic interaction between the particles is negligible and the single particle Langevin equation, describing the motion of a charge in external fluctuating field, is fully justified. The variance of r and more generally all its cumulants provide the required information on the physical process. The investigation of the stochastic continuity equation (SCE), we propose, extends the previous works [1,2,4] on the stochastic Liouville equation and shows that the moments hrn i satisfy a simple generalization of the Fokker–Planck equation satisfied by hri:
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
19
We first analyze the simple case of a free particle with white noise and Gaussian initial distribution. The average of the density and any of its integer powers hrn i satisfies the same Fokker–Planck equation. The direct proof is based on a series expansion and the average of even powers of a Wiener process, which have a simple form. The extension from the free particle case to any vector field with an additive or multiplicative white noise is based on a Dyson expansion in the interaction representation, a method first suggested in Ref. [1] and used in Ref. [5] to obtain a generalized Fokker–Planck equation in the case of correlated noise. When the vector field is volume preserving, a nontrivial case in dimension dX2; the average of any integer power of the fluctuating density hrn i satisfies the Fokker–Planck equation, with the appropriate initial condition rnin : Finally, we consider an extended phase space ðx; rÞ and we investigate the evolution of a probability density Pðx; r; tÞ associated to the fluctuating variables x; r: The probability density Pðr; x; tÞ has the property that its moments of order n with respect to r coincide with the averages of the nth power of the fluctuating density rðx; t; xÞ: The function P satisfies the Fokker–Planck equation associated to the generalized Langevin equation in the extended ðx; rÞ phase space. Concerning the mathematical properties of the fluctuating density we notice that at any point x of the phase space the average hri is the transition probability from the initial point to the point x of the process defined by Langevin equation. However, the mean square deviation s2 ¼ hr2 i hri2 and higher cumulants of the fluctuating density are not zero. This is easy to verify in the free particle case where s2 is a bimodal distribution vanishing for t ! 1: A fluctuating density r; which evolves according to the SCE, differs from the standard probability density associated to Langevin equation. The difference is evident since in the first case for any realization of the stochastic process there is a rigid transport of all the points in the support of the initial distribution, whereas the Langevin equation propagates any initial point with an independent realization of the stochastic process. In some problems related to accelerators and plasma physics one can be interested in evaluating the probability that some particles reach a specific region in phase space. A Brownian molecular dynamics (BMD) approach to the Langevin equation usually requires many realizations of the noise to keep the statistical error small, since the whole distribution function is computed. A direct integration of the Fokker–Planck equation using a fixed grid in phase space (Liouville approach) can introduce interpolating errors that destroy some structures of the distribution due for instance to the Hamiltonian character of the deterministic dynamics. Once the initial distribution is known analytically the BMD approach allows to compute the average of the fluctuating density, solution of the SCE, in the region of interest, without any grid interpolation. The density at the given point is equal to the initial density at the back propagated point for the volume preserving flows and the error due to the averaging process is inversely proportional to the number of realizations. This paper is organized as follows: in Section 1 we consider the SCE for the density function rðx; t; xÞ and its integer powers, and the equation satisfied by their averages; in Section 2 we analyze in detail the case of a free particle subject to a white noise; in Section 3 we compare the BMD solution of the Fokker–Planck equation
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
20
with the average hrðx; t; xÞi of the SCE, analyzing the related errors; in Section 4 we discuss the space correlations and the equations they satisfy; in Section 5 we introduce the probability density Pðr; x; tÞ; the equations it satisfies and its moments; in Appendix B we recall the basic properties of the Wiener and Uhlenbeck processes; in Appendix C we compute the space correlations for the free particle and damped particle with white noise; and in Appendix D we compute the average hr2 i for a damped particle with a white or correlated noise.
2. Stochastic continuity equation The Langevin equation describes the evolution of a point subject to a vector field F ¼ a þ bxðtÞ in a phase space Rd where a is the deterministic component and xðtÞ is a white noise, formal time derivative of a Wiener process wðtÞ dx ¼ Fðx; t; xÞ aðx; tÞ þ bðx; tÞxðtÞ : dt
(1)
In (1) x; F; a; b denote vectors in Rd ; whereas xðtÞ is the usual scalar white noise. A regularization of the stochastic process, achieved by replacing the white noise xðtÞ with a Uhlenbeck process, see Appendix A, would allow to make our formal derivation rigorous by letting the correlation time vanish as a final step. To the stochastic differential equation we associate the stochastic continuity equation, qr qðFrÞ þ ¼0; qt qx
(2)
which is well defined for any realization of a regular stochastic process. Given an initial distribution rin ðxÞ and letting x ¼ S t;t0 ðx0 ; xÞ be the evolution semi-group, the solution can be written as rðx; t; xÞ ¼
rin ðx0 Þ ; mðx; t; xÞ
(3)
where m is the Jacobian of the flow from x0 at time t0 to x at time t: In the volume preserving case qF=qx ¼ 0 the solution simplifies, since r is constant along any trajectory, and reads rðx; t; xÞ ¼ rin ðx0 Þ :
(4)
This is the case of the free particle F ¼ xðtÞ; we shall analyze in detail, where rðx; t; xÞ ¼ rin ðx wðtÞ þ wðt0 ÞÞ :
(5)
In the general case we write the continuity equation (2) in the following operator form: qr ¼ ðAðtÞ þ xðtÞBðtÞÞr; qt
A¼
q aðx; tÞ; qx
B¼
q bðx; tÞ qx
(6)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
and its solution as a Dyson series according to ! Z t rðx; t; xÞ ¼ T exp ððAðsÞ þ xðsÞBðsÞÞ ds rin ðxÞ ;
21
(7)
t0
where T is the time ordering operator defined by TðOðt1 ÞOðt2 ÞÞ ¼ Oðt2 ÞOðt1 Þ if t2 Xt1 and TðOðt1 ÞOðt2 ÞÞ ¼ Oðt1 ÞOðt2 Þ if t1 4t2 : In Ref. [5] it was proved that density averaged with respect to the white noise x satisfies the following equation: q hri ¼ qt
B2 Aþ hri 2
(8)
which is just the Fokker–Planck equation corresponding to the Stratonovich interpretation of the Langevin equation (1). In order to characterize the stochastic process rðx; t; xÞ the average is clearly not sufficient and to this end we consider the integer powers of the density and their averages hrn ðx; t; xÞi: 2.1. The moments In the case of a volume preserving vector field qF=qx ¼ 0 the integer powers of the density rn ðx; t; xÞ for n41 satisfy the same equation as for n ¼ 1; qrn qrn þF ¼0 qt qx
(9)
the initial conditions being rn ðx; t0 ; xÞ ¼ rnin ðxÞ
(10)
as a consequence the solution can be written as a Dyson series [6] ! Z t n r ðx; t; xÞ ¼ T exp ððAðsÞ þ xðsÞBðsÞÞ ds rnin ðxÞ
(11)
t0
and their averages still satisfy the equation q n B2 hr i ¼ A þ hrn i qt 2
(12)
with initial conditions hrn ðx; t0 ; xÞi ¼ rnin ðxÞ
(13)
q q and B ¼ b qx for a volume preserving flow. In the case Note that in (12) A ¼ a qx of a generic flow equation (9) is replaced by
qrn qrn qF ¼0 þF þ nrn qx qt qx
(14)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
22
and the equation for the averages is still given by (11) where the operators A and B are given by An ¼ ðn 1Þ
qa q a; qx qx
Bn ¼ ðn 1Þ
qb q b: qx qx
(15)
3. The free particle example We consider the simplest Langevin equation defined by dx ¼ xðtÞ dt
xðtÞ ¼ x0 þ wðtÞ wðt0 Þ
and a Gaussian initial distribution 1 x2 rin ðxÞ ¼ pffiffiffiffiffiffiffiffiffi exp 2t0 2pt0
(16)
(17)
In this case we know that the density rðx; t; xÞ is given by (5) namely 1 ðx wðtÞ þ wðt0 ÞÞ2 : rðx; t; xÞ ¼ pffiffiffiffiffiffiffiffiffi exp 2t0 2pt0
(18)
Letting the averages of the nth power be labeled according to rðnÞ ðx; tÞ ¼ hrn ðx; t; xÞi
(19)
by a straightforward, even though tedious calculation reported in Appendix A, we show that 1 x2 ð1Þ p ffiffiffiffiffiffiffi r ðx; tÞ ¼ exp : (20) 2t 2pt By the same procedure we prove that rðnÞ ðx; tÞ satisfy the Fokker–Planck equation qrðnÞ 1 q2 rðnÞ ¼ 2 qx2 qt
(21)
with initial conditions ðnÞ
r ðx; t0 Þ ¼
rnin ðxÞ
¼
1 ð2pt0 Þn=2
nx2 exp 2t0
(22)
The computation of hrð2Þ i follows a similar procedure and the result is given by 1 x2 rð2Þ ðx; tÞ ¼ pffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp (23) 2t t0 2p t0 2t t0 It is not hard to verify that it satisfies the Fokker–Planck equation (21) with initial condition (22), where n ¼ 2: In a similar way we can prove that rðnÞ ðx; tÞ is
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
23
given by pffiffiffiffi t0 x2 rðnÞ ðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ; 2ðt tn Þ ð2pt0 Þ nðt tn Þ where tn ¼ t0
(24)
n1 n :
4. Numerical results and error analysis It is immediate to check that the mean square deviation of the fluctuating density rðx; t; xÞ is not zero S2 ðx; tÞ ¼ rð2Þ ðx; tÞ ðrð1Þ ðx; tÞÞ2 a0 :
(25)
A simple computation shows that the first two even moments with respect to x (the odd ones disappear) of S2 ðx; tÞ are given by Z þ1 1 1 1 2 dx S ðx; tÞ ¼ pffiffiffi pffiffiffiffi pffiffi h1iS2
(26) t0 2 p t 1 and hx2 iS2
Z
þ1 1
1 2t t0 pffiffi pffiffiffiffi t : dx x2 S2 ðx; tÞ ¼ pffiffiffi t0 4 p
(27)
Both disappear for t ¼ t0 ; whereas for t ! 1 the former h1iS2 reaches a constant limit, the latter hx2 iS2 diverges linearly. In Fig. 1 we show the behavior of rð1Þ ðx; tÞ and the SCE variance Sðx; tÞ for a Gaussian initial distribution (17) with t0 ¼ 1 using the analytical results (20), (23), (25). It is interesting to observe that the SCE variance is a bimodal distribution with a minimum at the origin. In Fig. 2 we plot the initial density corresponding to the Gaussian initial distribution rðx; t0 Þ given by (17) for t0 ¼ 12; BMD simulated with N ¼ 2 104 points, and its evolution hrðx; t0 Þi at time t ¼ 2; obtained by propagating each initial
Fig. 1. Left panel: plot of the mean rð1Þ ðx; tÞ of the density rðx; t; xÞ solution of the stochastic Liouville 2 equation for a free particle with white noise for an initial Gaussian density rin ðxÞ ¼ ex =2 =ð2pÞ1=2 : The three curves correspond to t ¼ 1:01 (green), 1.1 (red), and 2 (blue), respectively. Right panel: plot of the variance sðx; tÞ for the same values of t and the same initial condition as the right panel.
ARTICLE IN PRESS 24
G. Bassi et al. / Physica A 347 (2005) 17–37
Fig. 2. Left: one-dimensional diffusion for different realizations of the noise (Brownian motion): comparison between the exact and numerical solution at time t ¼ 0:5; 2 with N ¼ 2 104 : Right: onedimensional diffusion for same realizations of the noise (stochastic Liouville equation): initial distribution at t ¼ 0:5 (red), mean distribution at t ¼ 100 (blue), distributions obtained for one realization of the noise (black).
Fig. 3. Left: one-dimensional diffusion: comparison between the exact and numerical mean density at time t ¼ 2: In red the initial distribution at t ¼ 0:5 with N ¼ 2 104 : Right: the same for the variance of the density.
point with a distinct Wiener trajectory. With such a value of N the discrepancy, with respect to the exact result due to statistical error, is visible at the initial time and at the next one. In the next frame of Fig. 2 five realizations of the process rðx; t; xÞ are shown; they correspond to the analytical initial distribution translated according to five realizations of the Wiener process. The relative error is proportional to N 1=2 where N is the number of realizations, as one can show using the central limit theorem. In Fig. 3 the SCE average density is compared with the exact result. The numerical procedure consists in computing the fluctuation density rðx; t; xÞ according to Eq. (5) on N ¼ 2 104 distinct Wiener trajectories and taking the average. The error cannot be appreciated within the graphical resolution; this is also true for the mean square deviation S2 ðx; tÞ which is compared with the exact result in the right frame.
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
25
4.1. Error on hri from SCE For any volume preserving vector field the fluctuating density rðx; t; xÞ is given by the initial density rin ; we suppose to be analytically known, at the initial point x0 ¼ S t0 ;t ðx; xÞ obtained by back propagating from t to t0 the point x with a suitable stochastic integrator. The relative error on the average density rð1Þ hrðx; tÞi; obtained by averaging over N stochastic trajectories, is proportional to N 1=2 : If the field is not volume preserving the knowledge of the Jacobian function mðx; tÞ is also required. The use of the SCE is convenient if we are interested in the average density at a given point x; since we can choose N large enough to reach the desired accuracy. This advantage is lost if we need to evaluate the average density on a domain whose measure is close to 1, since the computational load is comparable with BMD, unless the noise is very weak. Indeed, letting be the noise amplitude we have hri ¼ r0 þ 2 r2 and using the SCE we introduce an error only in average r2 of the fluctuating component, supposing we can neglect the error of our integrator on the deterministic trajectory. 4.2. Error on hri from BMD The solution of the Langevin equation is obtained by approximating the initial density with N particles and letting each of them evolve in time with an independent realization of the noise. The error on the density is due to the statistical fluctuation pffiffiffi Dn of the number n particles in a given cell and it is well known that Dn ¼ n: This error depends both on the initial distribution and on the stochastic trajectories; the error on the initial distribution, which does not affect the SCE if rin is known analytically, is unavoidable. Let us consider for example the one-dimensional case and suppose that the support of the particles density is an interval of length L in the time interval ½t0 ; t: We subdivide the space interval into equal subintervals of length Dx: The accuracy in x is determined by the space resolution Zx ¼ Dx=L: The numerical density at x is proportional to the number n of points falling in the interval which contains x (we may suppose it to be the midpoint) and the error Zr ¼ Dhri=hri is given by the fluctuation of n Zr ’ Zx ’ Z;
Zr ¼
Dn 1 ¼ pffiffiffi : n n
To have an explicit estimate we suppose that hri is slowly varying; in this case we may replace r with its space average 1=L and write n ’ N Dx=L ¼ NZx : Since the above equation yields n ’ 1=Z2r we finally obtain 1 ’ NZx Z2r ’ NZ3 ;
Zr ’ Zx ’ Z
if we impose that the space and density resolutions agree with a unique given accuracy Z in order to have a good balance. To resume with the BMD we need N ¼ 1=Z3 initial points, to which correspond N independent stochastic trajectories, in order to reach an accuracy Z in a grid with 1=Z
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
26
intervals. With the SCE we need N ¼ 1=Z2 points to reach the same accuracy for the average density at a single point. However, if we wish to evaluate the density at all the 1=Z grid points, the number of trajectories to be evaluated is 1=Z3 just as with the BMD method. For the free particle case, using (18) for the SCE solution, we have analyzed the ð1Þ average density error jrð1Þ num ðx; tÞ=r ðx; tÞ 1j taking a further average on the time or space interval considered. Choosing t0 ¼ 1; t ¼ 2 and x 2 ½3; 3 we found an error cN 1=2 with c of order 1 ð0:5oco0:8Þ: For the weakly noisy pendulum we have c51 as long as the average of the fluctuating component is a small fraction of the deterministic part.
5. The noisy pendulum The Hamiltonian systems with weak noise are physically significant, and mathematically challenging since the diffusion process occurs jointly with the filamentation process and the distributions are extremely hard to describe using the standard BMD methods, because the statistical error usually makes the fine structure of the probability distribution difficult to study. On the contrary the SLE allows to obtain very accurate results and hence can be used for instance to check the numerical integration schemes used to solve the Fokker–Planck equation [8]. The Langevin equations of our noisy pendulum, which describe the longitudinal motion in an RF cavity of HERAp, are dx ¼ p; dt
pffiffiffiffiffiffiffi dp ¼ A sin x B sinð4xÞ þ 2DxðtÞ ; dt
(28)
where we choose A ¼ 1:2875 and B ¼ 4:1201: In Fig. 4 we compare the phase space density for an initial Gaussian distribution in the noiseless case and the case of a
2
2
p x Fig. 4. Pendulum with an initial Gaussian distribution rðx; p; 0Þ ¼ expð 2s 2 2s2 Þ=ð2psx sp Þ ðsx ¼ x p 0:2; sp ¼ 0:4Þ at time t ¼ 30 without noise and with noise of small amplitude ¼ ð2DÞ1=2 with D ¼ 4 10 : Left: space distribution rðx; 0; tÞ: Right: momentum distribution rð0; p; tÞ: The results are obtained with an accurate integration of the Fokker–Planck equation for Da0:
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
27
Fig. 5. Phase space density distribution for the double pendulum with an initial Gaussian distribution 2
2
p x rðx; p; 0Þ ¼ expð 2s 2 2s2 Þ=ð2psx sp Þ ðsx ¼ 0:2; sp ¼ 0:4Þ at time t ¼ 30 with noise of small amplitude x
p
¼ ð2DÞ1=2 where D ¼ 104 : Left: comparison of the ‘‘exact’’ space distribution rðx; 0; tÞ with the solution obtained by solving the Langevin equation with the BMD method with N ¼ 106 particles. Right: comparison of the exact solution with the SCE obtained with 104 realizations of the noise.
weak noise. As it can be seen at time t ¼ 30 the structures due to filamentation are not yet destroyed by the noise. These structures are almost impossible to reproduce with the BMD solution of Eq. (28) unless the number of particles is large ð106 –107 Þ; whereas they are reproduced with a very small error by the SCE with few realizations of the noise. In Fig. 5 the comparison of the ‘‘exact’’ space distribution rðx; 0; tÞ with the solution obtained by solving the Langevin equation with the BMD method with N ¼ 106 particles (right frame) and with the SCE obtained with 104 realizations of the noise (left frame) is shown. While the difference in the distributions plotted in the right frame cannot be appreciated within the graphical resolution, the difference in the left frame is quite remarkable. In Fig. 5 we have done an error analysis for the space distributions rðx; 0; tÞ:
6. Space correlations The fluctuating densities are space correlated. The two and more generally the n points space correlation function for the stochastic continuity equations (SCE) satisfy Fokker–Planck-like equations which are obtained in a way very similar to the one followed for the average of the density and its integer powers. Indeed, starting from Eq. (6) for the SCE we write q rðx; tÞ ¼ ðAx ðtÞ þ xðtÞBx ðtÞÞrðx; tÞ; qt q bðx; tÞ Bx ðtÞ ¼ qx
Ax ðtÞ ¼
q aðx; tÞ ; qx ð29Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
28
it is not hard to show that the product of the density at two distinct points x; y satisfies the following equation: q rðx; tÞrðy; tÞ ¼ ðAðtÞ þ xðtÞBðtÞÞrðx; tÞrðy; tÞ; qt
A ¼ Ax þ Ay ;
B ¼ Bx þ By : (30)
q 2 qx r ðxÞ;
We should notice that ArðxÞrðyÞ in the limit y ! x becomes taking the limit separately on the operator and the function gives a wrong result. The solution can be written in the same way as (7) where the operators A; B are defined according to (30) and the initial condition is replaced by rin ðxÞrin ðyÞ: Taking the average and denoting for brevity r2 ðx; y; tÞ ¼ hrðx; tÞrðy; tÞi the two-point space correlation function, the equation it satisfies is q B2 r ðx; y; tÞ ¼ A þ (31) r2 ðx; y; tÞ qt 2 2 with initial condition r2 ðx; y; 0Þ ¼ rin ðxÞrin ðyÞ: More generally denoting rn ðx1 ; . . . ; xn ; tÞ ¼ hr1 ðx; tÞ rn ðxn ; tÞi the n points correlation function, the equation it satisfies is similar to (31) and reads q B2 rn ðx1 ; . . . ; xn ; tÞ ¼ A þ rn ðx1 ; . . . ; xn ; tÞ ; qt 2 A¼
n X q aðxi ; tÞ; qxi i¼1
B¼
n X q bðxi ; tÞ qxi i¼1
(32)
with initial conditions rn ðx1 ; . . . ; xn ; 0Þ ¼ rin ðx1 Þ rin ðxn Þ:
7. The probability distribution The process associated to the fluctuating density rðx; t; xÞ can be determined from a probabilistic viewpoint by associating a probability density Pðx; r; tÞ: At each time t indeed our field is specified by a space coordinate x and the value r it takes there. We may also introduce a conditional probability Pðx; r; t j x0 ; r0 ; t0 Þ equal to the weighted sum of all the paths joining in an extended phase space the initial point ðx0 ; r0 ; t0 Þ with the final point ðx; r; tÞ: The phase space ðx; rÞ 2 Rdþ1 has one extra dimension with respect to the ordinary phase space x 2 Rd where the Langevin equation is defined. The initial probability distribution is P0 ðx; rÞ ¼ Pðx; r; t0 Þ; not to be confused with the initial density rin ðxÞ: If the initial density is rin ðxÞ with probability 1 then P0 ðx; rÞ ¼ dðr rin ðxÞÞ : The normalization condition reads for any tXt0 Z þ1 Z þ1 dr dx rPðx; r; tÞ ¼ 1 : 0
1
(33)
(34)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
The relation between the initial probability and its value at time t is Z þ1 Z þ1 Pðx; r; tÞ ¼ dr0 dx0 Pðx; r; t j x0 ; r0 ; t0 ÞP0 ðx0 ; r0 Þ ; 0
29
(35)
1
where the transition probability satisfies the standard normalization condition Z þ1 Z þ1 dr dx Pðx; r; t j x0 ; r0 ; t0 Þ ¼ 1 0
1
and the group property Z þ1 Z þ1 dr1 dx1 Pðx; r; t j x1 ; r1 ; t1 ÞPðx1 ; r1 ; t1 j x0 ; r0 ; t0 Þ 0
1
¼ Pðx; r; t j x0 ; r0 ; t0 Þ :
ð36Þ
To determine the transition probability we have to consider the characteristics of the SCE in the extended phase space ðx; rÞ: We assume for simplicity that qb=qx ¼ 0 so that the divergence of the vector field f is qa=qx ¼ 0: In this case the evolution equations read dx dr qa ¼ aðx; tÞ þ xðtÞbðx; tÞ; ¼ r : (37) dt dt qx We consider the continuity equation in the extended phase space for the fluctuating density Pðx; r; t; xÞ which satisfies the equation qP ¼ ðA þ Ar þ xðtÞBÞP ; qt where q qa q q aðx; tÞ; Ar ¼ r; B ¼ bðx; tÞ : qx qx qr qx Denoting by P¯ the average of P with respect to x; ¯ qP 1 2 ¯ ¼ A þ Ar þ Bn P : qt 2 A¼
(38)
(39)
(40)
¯ In the volume preservingR case Ar ¼ 0 and consequently the equation satisfied by P n¯ n and its moments hr i ¼ r P dr with respect to r is the Fokker–Planck equation (12). As a consequence we may identify P¯ with the probability density P in extended phase space. In the volume non preserving case it is not difficult to see that the equation satisfied by the moments hrn i is q n qa n 1 2 hr i þ hr i ¼ An þ B hrn i ; (41) qt qx 2 where An ¼ A ðn 1Þqa=qx is the operator defined by (15) and BRn ¼ B having ¯ dr ¼ assumed that qb=qx ¼ 0: This result easily follows by observing that rn Br P n n n qa=qx nhr i after an integration by parts. The relation with hr i is simple only if qa=qx is constant. In the one-dimensional case letting a be this scalar constant it
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
30
turns out that hrn i ¼ eaðtt0 Þ hrn i so that the probability density in extended phase ¯ space is given by P ¼ eaðtt0 Þ P: 7.1. The free particle with damping In order to clarify the above discussion we analyze in detail the case of the free particle with damping in which Eq. (37) reduces to x_ ¼ ax þ xðtÞ;
r_ ¼ ar
(42)
¯ is and the equation for P 2¯ q ¯ q ¯ a q ðrPÞ ¯ þ1 q P : ðxPÞ P¼a qt qx qr 2 qx2
(43)
The equation is solved by separation of variables and we write its fundamental solutions as ¯ r; t j x0 ; r0 ; t0 Þ ¼ Gðx; t j x0 ; t0 Þdðr r0 eaðtt0 Þ Þ ; Pðx;
(44)
where Gðx; t j x0 ; t0 Þ is the fundamental solution of the Fokker–Planck equation for the Uhlenbeck process 1 ðx x0 eaðtt0 Þ Þ2 Gðx; t j x0 ; t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ; 2s2 ðt t0 Þ 2ps2 ðt t0 Þ s2 ðtÞ ¼
1 e2at : 2a
ð45Þ
¯ the probability density in the extended As a consequence recalling that P ¼ eaðtt0 Þ P phase space at time t is given by Z þ1 Z þ1 aðtt0 Þ dr0 dx0 Gðx; t j x0 ; t0 Þ Pðx; r; tÞ ¼ e 0
1
dðr r0 eaðtt0 Þ ÞP0 ðx0 ; r0 Þ :
ð46Þ
For an initial probability distribution with support on a unique distribution rin ðxÞ we have P0 ðx; rÞ ¼ dðr rin ðxÞÞ: In Fig. 6 we plot Pðx; r; tÞ as a function of r having chosen x ¼ 1 and a ¼ 1 and rin ðxÞ equal to the Gaussian distribution (17). The moments of this probability distribution are easily computed Z þ1 Z þ1 Z þ1 rn dr dr0 dx0 hrn i ¼ eaðtt0 Þ 0
0
1 aðtt0 Þ
Gðx; t j x0 ; t0 Þdðr r0 e ÞP0 ðx0 ; r0 Þ Z þ1 Z þ1 dr0 enaðtt0 Þ rn0 dx0 Gðx; t j x0 ; t0 Þdðr0 rin ðx0 ÞÞ ¼ eaðtt0 Þ 0 1 Z þ1 dx0 Gðx; t j x0 ; t0 Þrnin ðx0 Þ : ¼ eðn1Þaðtt0 Þ 1
ð47Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
31
Fig. 6. Plot of the probability density Pðx; r; tÞ for the free particle with damping ða ¼ 1Þ and white noise as a function of r at x ¼ 1 and t ¼ 0:501; 0.51, 0.6, and 0.8 for an initial probability density at t ¼ 0:5 concentrated on a Gaussian. The circles correspond to the numerical result while the lines correspond to the analytical result.
In Appendix C we carry out the explicit calculation of the second moment hr2 i using (47) also in the case where the white noise xðtÞ is replaced by a correlated noise.
8. Conclusions The SCE associated to a Langevin equation is a more general process in which a given initial density function is propagated by any realization of the process. The probability density associated to the Langevin equation is the average of the density function r computed over all the realizations. Naturally, the average of the density function does not exhaust all the information, which is captured instead by the average of all the integer powers of the density. This density evolves following the Fokker–Planck equation associated to the Langevin equation for a volume preserving vector field. The averages are equivalent to the moments of a probability density in the phase space whose dimension is increased by one, by including the r coordinate. This probability also satisfies a Fokker–Planck equation associated to a new Langevin equation in the extended phase space, describing the characteristics of the continuity equation. The meaning of the process corresponding to the SCE has been illustrated in the case of the free particle where all the computations are easily
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
32
carried out. From the numerical viewpoint the calculation of the probability density, given initially in analytic form, is fast and accurate at a given space point, since one can average the value there of the fluctuating density. The standard numerical procedure requires in any case the propagation of all the points of the initial probability density, which is usually done with a BMD procedure computing for each point a stochastic trajectory, which is an autonomous realization. We did not analyze the case of a correlated noise, since the extension is rather straightforward, at least conceptually. One can treat the correlated noise as an additional variable, by including the white noise-driven differential equation into the previous set of differential equations. It is possible to avoid the increase of the phase space dimensionality; however, the probability densities no longer satisfy an equation in closed form and the approximation schemes required are applicable only to noise of small amplitude and with rapidly decaying correlations.
Acknowledgements We wish to thank Prof. Ellison for useful discussions. A. Bazzani and G. Turchetti wish to thank DESY for hospitality during the preparation of this work, which was completed within the framework of the COFIN project 2003 Order and chaos in nonlinear extended systems: coherent structures, weak stochasticity and anomalous transport. G. Bassi wishes to thank DESY for a doctoral fellowship and the US Department of Energy for a partial support by Contract DE-FG03-99ER41104.
Appendix A. Direct computation of rð1Þ ðx; tÞ This is achieved by expanding the exponential (18) and performing the averages by using (B.6) 1 1 X ð1Þn hðx wðtÞ wðt0 ÞÞ2n i rð1Þ ¼ pffiffiffiffiffiffiffiffiffi n! 2pt0 n¼0 ð2t0 Þn 1 n 2n 1 X ð1Þn X x2ðnkÞ k ðt t ¼ pffiffiffiffiffiffiffiffiffi Þ ð2k 1Þ!! 0 2k 2pt0 n¼0 ð2t0 Þn k¼0 n! 1 1 X 1 X ð1Þkþm k ð2k 1Þ!! 2m þ 2k ¼ pffiffiffiffiffiffiffiffiffi x2m ðt t Þ 0 kþm ðm þ kÞ! 2k 2pt0 m¼0 k¼0 ð2t0 Þ 1 1 X 1 X x2m t0 t k 1 ð2k 1Þ!! ð2m þ 2kÞ! m! ¼ pffiffiffi ð1Þm t0 ð2kÞ!ð2mÞ! p m¼0 2k ðm þ kÞ! ð2t0 Þmþ1=2 m! k¼0 1 1 X 1 X x2m t0 t k 1 ð2m þ 2k 1Þ!! 2m m! ¼ pffiffiffi ð1Þm t0 k! ð2mÞ! p m¼0 2k ð2t0 Þmþ1=2 m! k¼0
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
1 1 X 1 X x2m t0 t k 1 ð2m þ 2k 1Þ!! ¼ pffiffiffi ð1Þm t0 p m¼0 2k ð2m 1Þ!!k! ð2t0 Þmþ1=2 m! k¼0 1 1 X x2m t0 t m1=2 m ¼ pffiffiffi 1 ð1Þ t0 p m¼0 ð2t0 Þmþ1=2 m! 2 1 x ¼ pffiffiffiffiffiffiffi exp : 2t 2pt
33
ðA:1Þ
In summing the series we have used the following identity: 1 X
an
n¼0
n X
bnk ck ¼
1 X
bm
an cnm ;
(A.2)
n¼m
m¼0
k¼0
1 X
where we have set m ¼ n k and interchanged the summation order. From k ¼ n mX0 it follows mpno1: Changing again the summation index in the second sum from n to k ¼ n m it follows now that 0pko1 and consequently we have 1 X n¼0
an
n X
bnk ck ¼
1 X m¼0
k¼0
bm
1 X
akþm ck
(A.3)
k¼0
We have used the following formula for the Taylor expansion of an algebraic function: ð1 xÞm1=2 ¼
1 X k¼0
xk
1 1 ð2m þ 2k 1Þ!! : 2k k! ð2m 1Þ!!
(A.4)
The series converges within jyjo1; namely our previous series converges for
t0 t
(A.5)
t o1; 0oto2t0 : 0 Obviously for t42t0 the result still holds with an analytic continuation performed after summing the series.
Appendix B. Wiener and Uhlenbeck processes The time evolution of a stochastic process is determined by a deterministic field and its fluctuations (noise). The white noise xðtÞ; formal time derivative of a Wiener process wðtÞ; describes a fluctuating field, whose correlations have an instantaneous decay hxðtÞi ¼ 0;
hxðtÞxðt0 Þi ¼ dðt t0 Þ :
(B.1)
The presence of memory corresponds to a non instantaneous decay of correlations and the decay rate characterizes the fluctuating field. The fluctuating part wðtÞ of the
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
34
process xðtÞ defined by dx ¼ axðtÞ þ xðtÞ dt
(B.2)
is asymptotic to the Ornstein–Ulhenbeck process for large t and reduces to the Wiener process wðtÞ for a ¼ 0: As a consequence letting xðt0 Þ ¼ x0 be the initial condition we have Z t aðtt0 Þ x0 þ wðtÞ; wðtÞ ¼ eaðtsÞ xðsÞ ds (B.3) xðtÞ ¼ e t0
and for a ¼ 0 we have Z
t
wðtÞ ¼ wðtÞ wðt0 Þ
xðsÞ ds
(B.4)
t0
By using the following property of the white noise correlations: X hxðs1 Þxðs2 Þ xðs2n1 Þxðs2n Þi ¼ dðsi1 si2 Þ dðsi2n1 si2n Þ ;
(B.5)
i1 ...i2n
where the sum runs over all the ð2n 1Þ!! pairings of the indices, it is not difficult to show that the following relations hold for tXt0 : n 1 e2aðtt0 Þ hw2n ðtÞi ¼ ð2n 1Þ!! 2a hðwðtÞ wðt0 ÞÞ2n i ¼ ðt t0 Þn ð2n 1Þ!! :
ðB:6Þ
The most general Langevin equation can be written as dx ¼ Fðx; t; xÞ aðx; tÞ þ bðx; tÞxðtÞ ; dt
(B.7)
where the stochastic field may be regularized by replacing xðtÞ with awðtÞ; a regular process whose limit is the white noise when a ! 1:
Appendix C. Space correlations for a free and damped particle In order to illustrate the computation of space correlations, we consider the free particle case in which Eq. (31) explicitly reads q 1 q q 2 r2 ðx; y; tÞ ¼ þ r2 ðx; y; tÞ : (C.1) qt 2 qx qy In order to find the solution corresponding to rin ðxÞ ¼ dðxÞ we perform a Fourier decomposition Z 2 1 r2 ðx; ty; tÞ ¼ ð2pÞ2 e2ðkx þky Þ tikx xiky y dkx dky : (C.2) R2
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
35
It is evident that pffiffiffi this is a solution pofffiffiffi Eq. (C.1). Rotating the coordinates with u ¼ ðkx þ ky Þ= 2 and v ¼ ðkx ky Þ= 2 we obtain xþy 2 1 Z Z xþy xy eð 2 Þ 2t u2 tiu pffi i pffi v 2 2 2 r2 ðx; y; tÞ ¼ ð2pÞ e du e dv ¼ pffiffiffiffiffiffiffi dðx yÞ : (C.3) 2pt R R 2
In the initial distribution is a Gaussian r0 ðxÞ ¼ ð2pt0 Þ1=2 expð 2tx Þ we may calculate 0 the space correlation as rð2Þ ¼ hrin ðx wðtÞ þ wðt0 ÞÞrin ðy wðtÞ þ wðt0 ÞÞi rather than solving (C.1) Z 1 ð2Þ rin ðx wÞrin ðy wÞpðw; t t0 Þ dw ; (C.4) r ðx; y; tÞ ¼ 1
where 2 1 w pðw; t t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2ðtt0 Þ : 2pðt t0 Þ
(C.5)
It follows that Z hrðx; tÞrðy; tÞi ¼
1
2 ðxwÞ2 ðywÞ2 1 1 1 w dw pffiffiffiffiffiffiffiffiffi e 2t0 pffiffiffiffiffiffiffiffiffi e 2t0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2ðtt0 Þ 2pt0 2pt0 2pðt t0 Þ 1
tðx2 þy2 Þ2ðtt0 Þxy 1 ð Þ 2t0 ð2tt0 Þ pffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2p t0 2t t0 xþy 2 1 ðxyÞ2 1 1 ð Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2 2tt0 pffiffiffiffiffiffiffiffiffi e 4t0 : 4pt0 pð2t t0 Þ
¼
ðC:6Þ
Note that in the limit t0 ! 0 the new result (C.6) agrees with the previous one (C.3). Moreover, for y ¼ x we recover rð2Þ ðx; tÞ hr2 ðx; tÞi given by (23). C.1. Damped particle with white noise 2
x If the initial distribution is a Gaussian rin ðxÞ ¼ ð2ps20 Þ1=2 expð 2s 2 Þ we calculate 0
the space correlation for x_ ¼ ax þ xðtÞ which is given by (C.6) where 2 1 w pðw; t t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2sðtt0 Þ ; 2psðt t0 Þ
sðtÞ ¼
1 ð1 e2at Þ : 2a
(C.7)
It follows that hrðx; tÞrðy; tÞi ¼
e2aðtt0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ps20 2ps2 ðt t0 Þ
Z
1
dw e
12 ðe2aðtt0 Þ ½ðxwÞ2 þðywÞ2 Þ 2s 0
e
w2 2s2 ðtt0 Þ
1
aðtt0 Þ
¼
e qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 ðx yÞ2 s2 ðt t0 Þe2aðtt0 Þ þ ðx2 þ y2 Þs20 exp : 2s20 ð2s2 ðt t0 Þ þ e2aðtt0 Þ s20 Þ
ðC:8Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
36
For x ¼ y we obtain hr2 ðx; tÞi ¼
eaðtt0 Þ qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 x2 exp 2 2s ðt t0 Þ þ e2aðtt0 Þ s20
ðC:9Þ
in agreement with the result (D.4) obtained in Appendix D.
Appendix D. The second moment hr2 i for a damped particle In order to give a simple explicit example we consider a point subjected to a damping force and random uncorrelated collisions. The corresponding process for the point velocity, denoted by x; is defined by Eq. (B.2). The average density hri and the average of its square hr2 i satisfy the equations defined by (12) where A; B are replaced by An ; Bn defined by (15) for n ¼ 1; 2: These equations explicitly read q q 1 q2 hri ¼ a ðxhriÞ þ hri qt qx 2 qx2
(D.1)
and q 2 q 1 q2 2 hr i ¼ a ðxhr2 iÞ þ ahr2 i þ hr i ; (D.2) qt qx 2 qx2 which differs from the previous one by the divergence of the deterministic field. In this case there is no Langevin equation associated to the process defined by the fluctuating field r2 and in order to compute the mean hr2 i we evaluate the second moment with respect to r of the probability density Pðx; r; tÞ according to (47). Choosing an initial Gaussian density with variance s20 1 x2 rin ðxÞ ¼ qffiffiffiffiffiffiffiffiffiffi exp 2 (D.3) 2s0 2ps20 and inserting into (46) we obtain Z þ1 eaðtt0 Þ ðx x0 eaðtt0 Þ Þ2 x2 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p dx0 exp hr i ¼ exp 2 2s2 ðt t0 Þ 2s0 2ps20 2ps2 ðt t0 Þ 1 ¼
eaðtt0 Þ qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 x2 exp 2 : 2s ðt t0 Þ þ e2aðtt0 Þ s20
ðD:4Þ
If we choose s20 ¼ s2 ðt0 Þ then the denominators in the exponent and the argument of the second square root becomes s2 ðt t0 Þ þ s2 ðtÞ since the following identity holds
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
37
s2 ðt t0 Þ þ e2aðtt0 Þ s2 ðt0 Þ ¼ s2 ðtÞ: In the limit a ! 0 the result (30) for the second moment of the free particle is recovered. For the moment hrn i the result is obtained in a similar way. If the white noise is replaced by a Uhlenbeck noise according to dx dy ¼ ax þ yðtÞ; ¼ byðtÞ þ bxðtÞ (D.5) dt dt choosing yð0Þ ¼ 0 so that hyðtÞi ¼ 0 we have a new process, which tends to the previous one in the limit b ! 1: The process xðtÞ is Gaussian and writing the solution of (D.5) as Z t Z t at aðttÞ xðtÞ ¼ x0 e þ e yðtÞ dt yðtÞ ¼ b dtebðttÞ xðtÞ (D.6) 0
0 at
Rt the variance is easily computed after writing x x0 e ¼ 0 Kðt; sÞxðsÞ ds according to Z t s2 ðtÞ ¼ K 2 ðt; sÞ ds 0 2
b ð1 e2at Þ 2b3 ð1 eðaþbÞt Þ b3 ð1 e2bt Þ þ 2 ¼ : ðD:7Þ 2aðb aÞ2 ðb a2 Þðb2 abÞ 2ðb2 abÞ2 The evolution for r is still the same dr=dt ¼ ar and the probability density P is still ¯ is given by (44) and in (45) s2 ðtÞ is given by expressed by P ¼ eaðtt0 Þ P¯ where P (D.7). As a consequence the second moment hr2 i is still given by (D.4) where the variance is defined by (D.7). References [1] R. Kubo, in: D. Ter Haar (Ed.), Fluctuation, Relaxation and Resonance in Magnetic Systems, Oliver and Boyd, Edinburgh, 1962. [2] R. Kubo, Stochastic Liouville equations, J. Math. Phys. 4 (1963) 174–183. [3] R. Graham, Covariant formulation of non equilibrium statistical thermodynamics, Z. Phys. 26 (1977) 397–405. [4] J. Ellison, Accelerators and probability: the special effect of noise in beam dynamics, nonlinear and stochastic beam dynamics in accelerators—a challenge to theoretical and computational physics, DESY Proceedings, 1998. [5] A. Bazzani, G. Bassi, G. Turchetti, Diffusion and memory effects for stochastic processes and fractional Langevin equations, Physica A 324 (2003) 530–550. [6] F. Dyson, The radiation theories of Tomonaga, Schwinger and Feynmann, Phys. Rev. 75 (1949) 486–502. [7] R. Kubo, M. Toda, N. Hashitsume, Statistical Physics II, Non equilibrium Statistical Mechanics, Springer Verlag Series in Solid State Science 31 (1978) 25. [8] G. Bassi, A. Bazzani, H. Mais, G. Turchetti, Hamiltonian dynamics with weak noise and the echo effect, Physica A, submitted for publication.
Physica A 347 (2005) 17–37 www.elsevier.com/locate/physa
Stochastic continuity equation and related processes Gabriele Bassia, Armando Bazzanib, Helmut Maisc, Giorgio Turchettib, a Department of Mathematics and Statistics, University of New Mexico, USA Department of Physics, University of Bologna, INFN Sezione di Bologna, Via Irnerio 46, Bologna 40127, Italy c DESY, Hamburg, Germany
b
Received 13 December 2003; received in revised form 7 July 2004 Available online 27 October 2004
Abstract We consider the processes defined by a Langevin equation and the associated continuity equation. The average of the density function, solution of the continuity equation, satisfies the Fokker–Planck equation. For a volume preserving vector field the same equation is satisfied by the average of the integer powers of the density, which are the moments of the related probability density. For a generic vector field the Fokker–Planck equation for the moments is slightly modified. We first illustrate the problem in the simple case of a free particle subject to a white noise, since the averages can be computed by an elementary procedure using the factorization property of the correlation functions of the noise. The probabilistic meaning of the moments is discussed and the comparison between the analytical results and the numerical simulation is shown. The case of a generic Langevin equation is treated by computing the averages via a Dyson expansion after observing that, for a volume preserving vector field, any power of the density function satisfies the same continuity equation with the appropriate initial conditions. As a consequence the results obtained for the free particle are easily extended. An alternative approach is based on the characteristics of the continuity equation; the probability density of this process in an extended phase
Corresponding author. Tel.: +39 0512091095.
E-mail address: [email protected] (G. Turchetti). 0378-4371/$ - see front matter r 2004 Published by Elsevier B.V. doi:10.1016/j.physa.2004.08.083
ARTICLE IN PRESS 18
G. Bassi et al. / Physica A 347 (2005) 17–37
space still satisfies a Fokker–Planck equation and its moments coincide with the previous definition. r 2004 Published by Elsevier B.V. Keywords: Stochastic processes; Continuity equation; Moments
1. Introduction Stochastic processes are commonly used to build phenomenological models of physical and biological systems, in which a small subsystem is isolated taking into account its interaction with the fluctuating environment. The Langevin equation, continuous version of a random walk, is an ordinary first-order differential equation with a deterministic vector field plus a white noise component. The probability density of the process satisfies the Fokker–Planck equation; its fundamental solution is the weighted sum over all paths connecting the initial and final points in the extended phase space [7]. We consider here the continuity equation associated to the Langevin equation. This is a first-order linear partial differential equation, which defines a stochastic process whose properties are analyzed by considering the average of its integer powers. These are the moments of the probability density associated to the stochastic process, and determine it completely. What we propose is an extension of the stochastic Liouville equation previously considered [1–3], and discussed in the context of beam dynamics [4]. It is important to remark that there are physical processes, such as the Brownian motion, in which each particle experiences a different realization of the noise (describing random collisions with the medium). In this case the physical process is described by the Langevin equation and, even though the stochastic continuity equation for the density r is mathematically well defined, only the average hri; satisfying the Fokker–Planck equation has a physical meaning as the probability density of the process. The evaluation of r at a given location is useful from the numerical viewpoint since it allows to compute the average hri there with a high accuracy. However, there are physical processes in which all the particles experience the same noise: consider for instance a bunch of particles crossing a magnetic element in an accelerator, where the current has small random fluctuations. In this case the density r; satisfying the continuity equation, describes the physical process and hri has to be interpreted as the average over many distinct runs, each one providing a different realization of the stochastic process, for the same initial bunch. In a very high-energy low-intensity beam the electromagnetic interaction between the particles is negligible and the single particle Langevin equation, describing the motion of a charge in external fluctuating field, is fully justified. The variance of r and more generally all its cumulants provide the required information on the physical process. The investigation of the stochastic continuity equation (SCE), we propose, extends the previous works [1,2,4] on the stochastic Liouville equation and shows that the moments hrn i satisfy a simple generalization of the Fokker–Planck equation satisfied by hri:
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
19
We first analyze the simple case of a free particle with white noise and Gaussian initial distribution. The average of the density and any of its integer powers hrn i satisfies the same Fokker–Planck equation. The direct proof is based on a series expansion and the average of even powers of a Wiener process, which have a simple form. The extension from the free particle case to any vector field with an additive or multiplicative white noise is based on a Dyson expansion in the interaction representation, a method first suggested in Ref. [1] and used in Ref. [5] to obtain a generalized Fokker–Planck equation in the case of correlated noise. When the vector field is volume preserving, a nontrivial case in dimension dX2; the average of any integer power of the fluctuating density hrn i satisfies the Fokker–Planck equation, with the appropriate initial condition rnin : Finally, we consider an extended phase space ðx; rÞ and we investigate the evolution of a probability density Pðx; r; tÞ associated to the fluctuating variables x; r: The probability density Pðr; x; tÞ has the property that its moments of order n with respect to r coincide with the averages of the nth power of the fluctuating density rðx; t; xÞ: The function P satisfies the Fokker–Planck equation associated to the generalized Langevin equation in the extended ðx; rÞ phase space. Concerning the mathematical properties of the fluctuating density we notice that at any point x of the phase space the average hri is the transition probability from the initial point to the point x of the process defined by Langevin equation. However, the mean square deviation s2 ¼ hr2 i hri2 and higher cumulants of the fluctuating density are not zero. This is easy to verify in the free particle case where s2 is a bimodal distribution vanishing for t ! 1: A fluctuating density r; which evolves according to the SCE, differs from the standard probability density associated to Langevin equation. The difference is evident since in the first case for any realization of the stochastic process there is a rigid transport of all the points in the support of the initial distribution, whereas the Langevin equation propagates any initial point with an independent realization of the stochastic process. In some problems related to accelerators and plasma physics one can be interested in evaluating the probability that some particles reach a specific region in phase space. A Brownian molecular dynamics (BMD) approach to the Langevin equation usually requires many realizations of the noise to keep the statistical error small, since the whole distribution function is computed. A direct integration of the Fokker–Planck equation using a fixed grid in phase space (Liouville approach) can introduce interpolating errors that destroy some structures of the distribution due for instance to the Hamiltonian character of the deterministic dynamics. Once the initial distribution is known analytically the BMD approach allows to compute the average of the fluctuating density, solution of the SCE, in the region of interest, without any grid interpolation. The density at the given point is equal to the initial density at the back propagated point for the volume preserving flows and the error due to the averaging process is inversely proportional to the number of realizations. This paper is organized as follows: in Section 1 we consider the SCE for the density function rðx; t; xÞ and its integer powers, and the equation satisfied by their averages; in Section 2 we analyze in detail the case of a free particle subject to a white noise; in Section 3 we compare the BMD solution of the Fokker–Planck equation
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
20
with the average hrðx; t; xÞi of the SCE, analyzing the related errors; in Section 4 we discuss the space correlations and the equations they satisfy; in Section 5 we introduce the probability density Pðr; x; tÞ; the equations it satisfies and its moments; in Appendix B we recall the basic properties of the Wiener and Uhlenbeck processes; in Appendix C we compute the space correlations for the free particle and damped particle with white noise; and in Appendix D we compute the average hr2 i for a damped particle with a white or correlated noise.
2. Stochastic continuity equation The Langevin equation describes the evolution of a point subject to a vector field F ¼ a þ bxðtÞ in a phase space Rd where a is the deterministic component and xðtÞ is a white noise, formal time derivative of a Wiener process wðtÞ dx ¼ Fðx; t; xÞ aðx; tÞ þ bðx; tÞxðtÞ : dt
(1)
In (1) x; F; a; b denote vectors in Rd ; whereas xðtÞ is the usual scalar white noise. A regularization of the stochastic process, achieved by replacing the white noise xðtÞ with a Uhlenbeck process, see Appendix A, would allow to make our formal derivation rigorous by letting the correlation time vanish as a final step. To the stochastic differential equation we associate the stochastic continuity equation, qr qðFrÞ þ ¼0; qt qx
(2)
which is well defined for any realization of a regular stochastic process. Given an initial distribution rin ðxÞ and letting x ¼ S t;t0 ðx0 ; xÞ be the evolution semi-group, the solution can be written as rðx; t; xÞ ¼
rin ðx0 Þ ; mðx; t; xÞ
(3)
where m is the Jacobian of the flow from x0 at time t0 to x at time t: In the volume preserving case qF=qx ¼ 0 the solution simplifies, since r is constant along any trajectory, and reads rðx; t; xÞ ¼ rin ðx0 Þ :
(4)
This is the case of the free particle F ¼ xðtÞ; we shall analyze in detail, where rðx; t; xÞ ¼ rin ðx wðtÞ þ wðt0 ÞÞ :
(5)
In the general case we write the continuity equation (2) in the following operator form: qr ¼ ðAðtÞ þ xðtÞBðtÞÞr; qt
A¼
q aðx; tÞ; qx
B¼
q bðx; tÞ qx
(6)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
and its solution as a Dyson series according to ! Z t rðx; t; xÞ ¼ T exp ððAðsÞ þ xðsÞBðsÞÞ ds rin ðxÞ ;
21
(7)
t0
where T is the time ordering operator defined by TðOðt1 ÞOðt2 ÞÞ ¼ Oðt2 ÞOðt1 Þ if t2 Xt1 and TðOðt1 ÞOðt2 ÞÞ ¼ Oðt1 ÞOðt2 Þ if t1 4t2 : In Ref. [5] it was proved that density averaged with respect to the white noise x satisfies the following equation: q hri ¼ qt
B2 Aþ hri 2
(8)
which is just the Fokker–Planck equation corresponding to the Stratonovich interpretation of the Langevin equation (1). In order to characterize the stochastic process rðx; t; xÞ the average is clearly not sufficient and to this end we consider the integer powers of the density and their averages hrn ðx; t; xÞi: 2.1. The moments In the case of a volume preserving vector field qF=qx ¼ 0 the integer powers of the density rn ðx; t; xÞ for n41 satisfy the same equation as for n ¼ 1; qrn qrn þF ¼0 qt qx
(9)
the initial conditions being rn ðx; t0 ; xÞ ¼ rnin ðxÞ
(10)
as a consequence the solution can be written as a Dyson series [6] ! Z t n r ðx; t; xÞ ¼ T exp ððAðsÞ þ xðsÞBðsÞÞ ds rnin ðxÞ
(11)
t0
and their averages still satisfy the equation q n B2 hr i ¼ A þ hrn i qt 2
(12)
with initial conditions hrn ðx; t0 ; xÞi ¼ rnin ðxÞ
(13)
q q and B ¼ b qx for a volume preserving flow. In the case Note that in (12) A ¼ a qx of a generic flow equation (9) is replaced by
qrn qrn qF ¼0 þF þ nrn qx qt qx
(14)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
22
and the equation for the averages is still given by (11) where the operators A and B are given by An ¼ ðn 1Þ
qa q a; qx qx
Bn ¼ ðn 1Þ
qb q b: qx qx
(15)
3. The free particle example We consider the simplest Langevin equation defined by dx ¼ xðtÞ dt
xðtÞ ¼ x0 þ wðtÞ wðt0 Þ
and a Gaussian initial distribution 1 x2 rin ðxÞ ¼ pffiffiffiffiffiffiffiffiffi exp 2t0 2pt0
(16)
(17)
In this case we know that the density rðx; t; xÞ is given by (5) namely 1 ðx wðtÞ þ wðt0 ÞÞ2 : rðx; t; xÞ ¼ pffiffiffiffiffiffiffiffiffi exp 2t0 2pt0
(18)
Letting the averages of the nth power be labeled according to rðnÞ ðx; tÞ ¼ hrn ðx; t; xÞi
(19)
by a straightforward, even though tedious calculation reported in Appendix A, we show that 1 x2 ð1Þ p ffiffiffiffiffiffiffi r ðx; tÞ ¼ exp : (20) 2t 2pt By the same procedure we prove that rðnÞ ðx; tÞ satisfy the Fokker–Planck equation qrðnÞ 1 q2 rðnÞ ¼ 2 qx2 qt
(21)
with initial conditions ðnÞ
r ðx; t0 Þ ¼
rnin ðxÞ
¼
1 ð2pt0 Þn=2
nx2 exp 2t0
(22)
The computation of hrð2Þ i follows a similar procedure and the result is given by 1 x2 rð2Þ ðx; tÞ ¼ pffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp (23) 2t t0 2p t0 2t t0 It is not hard to verify that it satisfies the Fokker–Planck equation (21) with initial condition (22), where n ¼ 2: In a similar way we can prove that rðnÞ ðx; tÞ is
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
23
given by pffiffiffiffi t0 x2 rðnÞ ðx; tÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffin pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ; 2ðt tn Þ ð2pt0 Þ nðt tn Þ where tn ¼ t0
(24)
n1 n :
4. Numerical results and error analysis It is immediate to check that the mean square deviation of the fluctuating density rðx; t; xÞ is not zero S2 ðx; tÞ ¼ rð2Þ ðx; tÞ ðrð1Þ ðx; tÞÞ2 a0 :
(25)
A simple computation shows that the first two even moments with respect to x (the odd ones disappear) of S2 ðx; tÞ are given by Z þ1 1 1 1 2 dx S ðx; tÞ ¼ pffiffiffi pffiffiffiffi pffiffi h1iS2
(26) t0 2 p t 1 and hx2 iS2
Z
þ1 1
1 2t t0 pffiffi pffiffiffiffi t : dx x2 S2 ðx; tÞ ¼ pffiffiffi t0 4 p
(27)
Both disappear for t ¼ t0 ; whereas for t ! 1 the former h1iS2 reaches a constant limit, the latter hx2 iS2 diverges linearly. In Fig. 1 we show the behavior of rð1Þ ðx; tÞ and the SCE variance Sðx; tÞ for a Gaussian initial distribution (17) with t0 ¼ 1 using the analytical results (20), (23), (25). It is interesting to observe that the SCE variance is a bimodal distribution with a minimum at the origin. In Fig. 2 we plot the initial density corresponding to the Gaussian initial distribution rðx; t0 Þ given by (17) for t0 ¼ 12; BMD simulated with N ¼ 2 104 points, and its evolution hrðx; t0 Þi at time t ¼ 2; obtained by propagating each initial
Fig. 1. Left panel: plot of the mean rð1Þ ðx; tÞ of the density rðx; t; xÞ solution of the stochastic Liouville 2 equation for a free particle with white noise for an initial Gaussian density rin ðxÞ ¼ ex =2 =ð2pÞ1=2 : The three curves correspond to t ¼ 1:01 (green), 1.1 (red), and 2 (blue), respectively. Right panel: plot of the variance sðx; tÞ for the same values of t and the same initial condition as the right panel.
ARTICLE IN PRESS 24
G. Bassi et al. / Physica A 347 (2005) 17–37
Fig. 2. Left: one-dimensional diffusion for different realizations of the noise (Brownian motion): comparison between the exact and numerical solution at time t ¼ 0:5; 2 with N ¼ 2 104 : Right: onedimensional diffusion for same realizations of the noise (stochastic Liouville equation): initial distribution at t ¼ 0:5 (red), mean distribution at t ¼ 100 (blue), distributions obtained for one realization of the noise (black).
Fig. 3. Left: one-dimensional diffusion: comparison between the exact and numerical mean density at time t ¼ 2: In red the initial distribution at t ¼ 0:5 with N ¼ 2 104 : Right: the same for the variance of the density.
point with a distinct Wiener trajectory. With such a value of N the discrepancy, with respect to the exact result due to statistical error, is visible at the initial time and at the next one. In the next frame of Fig. 2 five realizations of the process rðx; t; xÞ are shown; they correspond to the analytical initial distribution translated according to five realizations of the Wiener process. The relative error is proportional to N 1=2 where N is the number of realizations, as one can show using the central limit theorem. In Fig. 3 the SCE average density is compared with the exact result. The numerical procedure consists in computing the fluctuation density rðx; t; xÞ according to Eq. (5) on N ¼ 2 104 distinct Wiener trajectories and taking the average. The error cannot be appreciated within the graphical resolution; this is also true for the mean square deviation S2 ðx; tÞ which is compared with the exact result in the right frame.
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
25
4.1. Error on hri from SCE For any volume preserving vector field the fluctuating density rðx; t; xÞ is given by the initial density rin ; we suppose to be analytically known, at the initial point x0 ¼ S t0 ;t ðx; xÞ obtained by back propagating from t to t0 the point x with a suitable stochastic integrator. The relative error on the average density rð1Þ hrðx; tÞi; obtained by averaging over N stochastic trajectories, is proportional to N 1=2 : If the field is not volume preserving the knowledge of the Jacobian function mðx; tÞ is also required. The use of the SCE is convenient if we are interested in the average density at a given point x; since we can choose N large enough to reach the desired accuracy. This advantage is lost if we need to evaluate the average density on a domain whose measure is close to 1, since the computational load is comparable with BMD, unless the noise is very weak. Indeed, letting be the noise amplitude we have hri ¼ r0 þ 2 r2 and using the SCE we introduce an error only in average r2 of the fluctuating component, supposing we can neglect the error of our integrator on the deterministic trajectory. 4.2. Error on hri from BMD The solution of the Langevin equation is obtained by approximating the initial density with N particles and letting each of them evolve in time with an independent realization of the noise. The error on the density is due to the statistical fluctuation pffiffiffi Dn of the number n particles in a given cell and it is well known that Dn ¼ n: This error depends both on the initial distribution and on the stochastic trajectories; the error on the initial distribution, which does not affect the SCE if rin is known analytically, is unavoidable. Let us consider for example the one-dimensional case and suppose that the support of the particles density is an interval of length L in the time interval ½t0 ; t: We subdivide the space interval into equal subintervals of length Dx: The accuracy in x is determined by the space resolution Zx ¼ Dx=L: The numerical density at x is proportional to the number n of points falling in the interval which contains x (we may suppose it to be the midpoint) and the error Zr ¼ Dhri=hri is given by the fluctuation of n Zr ’ Zx ’ Z;
Zr ¼
Dn 1 ¼ pffiffiffi : n n
To have an explicit estimate we suppose that hri is slowly varying; in this case we may replace r with its space average 1=L and write n ’ N Dx=L ¼ NZx : Since the above equation yields n ’ 1=Z2r we finally obtain 1 ’ NZx Z2r ’ NZ3 ;
Zr ’ Zx ’ Z
if we impose that the space and density resolutions agree with a unique given accuracy Z in order to have a good balance. To resume with the BMD we need N ¼ 1=Z3 initial points, to which correspond N independent stochastic trajectories, in order to reach an accuracy Z in a grid with 1=Z
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
26
intervals. With the SCE we need N ¼ 1=Z2 points to reach the same accuracy for the average density at a single point. However, if we wish to evaluate the density at all the 1=Z grid points, the number of trajectories to be evaluated is 1=Z3 just as with the BMD method. For the free particle case, using (18) for the SCE solution, we have analyzed the ð1Þ average density error jrð1Þ num ðx; tÞ=r ðx; tÞ 1j taking a further average on the time or space interval considered. Choosing t0 ¼ 1; t ¼ 2 and x 2 ½3; 3 we found an error cN 1=2 with c of order 1 ð0:5oco0:8Þ: For the weakly noisy pendulum we have c51 as long as the average of the fluctuating component is a small fraction of the deterministic part.
5. The noisy pendulum The Hamiltonian systems with weak noise are physically significant, and mathematically challenging since the diffusion process occurs jointly with the filamentation process and the distributions are extremely hard to describe using the standard BMD methods, because the statistical error usually makes the fine structure of the probability distribution difficult to study. On the contrary the SLE allows to obtain very accurate results and hence can be used for instance to check the numerical integration schemes used to solve the Fokker–Planck equation [8]. The Langevin equations of our noisy pendulum, which describe the longitudinal motion in an RF cavity of HERAp, are dx ¼ p; dt
pffiffiffiffiffiffiffi dp ¼ A sin x B sinð4xÞ þ 2DxðtÞ ; dt
(28)
where we choose A ¼ 1:2875 and B ¼ 4:1201: In Fig. 4 we compare the phase space density for an initial Gaussian distribution in the noiseless case and the case of a
2
2
p x Fig. 4. Pendulum with an initial Gaussian distribution rðx; p; 0Þ ¼ expð 2s 2 2s2 Þ=ð2psx sp Þ ðsx ¼ x p 0:2; sp ¼ 0:4Þ at time t ¼ 30 without noise and with noise of small amplitude ¼ ð2DÞ1=2 with D ¼ 4 10 : Left: space distribution rðx; 0; tÞ: Right: momentum distribution rð0; p; tÞ: The results are obtained with an accurate integration of the Fokker–Planck equation for Da0:
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
27
Fig. 5. Phase space density distribution for the double pendulum with an initial Gaussian distribution 2
2
p x rðx; p; 0Þ ¼ expð 2s 2 2s2 Þ=ð2psx sp Þ ðsx ¼ 0:2; sp ¼ 0:4Þ at time t ¼ 30 with noise of small amplitude x
p
¼ ð2DÞ1=2 where D ¼ 104 : Left: comparison of the ‘‘exact’’ space distribution rðx; 0; tÞ with the solution obtained by solving the Langevin equation with the BMD method with N ¼ 106 particles. Right: comparison of the exact solution with the SCE obtained with 104 realizations of the noise.
weak noise. As it can be seen at time t ¼ 30 the structures due to filamentation are not yet destroyed by the noise. These structures are almost impossible to reproduce with the BMD solution of Eq. (28) unless the number of particles is large ð106 –107 Þ; whereas they are reproduced with a very small error by the SCE with few realizations of the noise. In Fig. 5 the comparison of the ‘‘exact’’ space distribution rðx; 0; tÞ with the solution obtained by solving the Langevin equation with the BMD method with N ¼ 106 particles (right frame) and with the SCE obtained with 104 realizations of the noise (left frame) is shown. While the difference in the distributions plotted in the right frame cannot be appreciated within the graphical resolution, the difference in the left frame is quite remarkable. In Fig. 5 we have done an error analysis for the space distributions rðx; 0; tÞ:
6. Space correlations The fluctuating densities are space correlated. The two and more generally the n points space correlation function for the stochastic continuity equations (SCE) satisfy Fokker–Planck-like equations which are obtained in a way very similar to the one followed for the average of the density and its integer powers. Indeed, starting from Eq. (6) for the SCE we write q rðx; tÞ ¼ ðAx ðtÞ þ xðtÞBx ðtÞÞrðx; tÞ; qt q bðx; tÞ Bx ðtÞ ¼ qx
Ax ðtÞ ¼
q aðx; tÞ ; qx ð29Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
28
it is not hard to show that the product of the density at two distinct points x; y satisfies the following equation: q rðx; tÞrðy; tÞ ¼ ðAðtÞ þ xðtÞBðtÞÞrðx; tÞrðy; tÞ; qt
A ¼ Ax þ Ay ;
B ¼ Bx þ By : (30)
q 2 qx r ðxÞ;
We should notice that ArðxÞrðyÞ in the limit y ! x becomes taking the limit separately on the operator and the function gives a wrong result. The solution can be written in the same way as (7) where the operators A; B are defined according to (30) and the initial condition is replaced by rin ðxÞrin ðyÞ: Taking the average and denoting for brevity r2 ðx; y; tÞ ¼ hrðx; tÞrðy; tÞi the two-point space correlation function, the equation it satisfies is q B2 r ðx; y; tÞ ¼ A þ (31) r2 ðx; y; tÞ qt 2 2 with initial condition r2 ðx; y; 0Þ ¼ rin ðxÞrin ðyÞ: More generally denoting rn ðx1 ; . . . ; xn ; tÞ ¼ hr1 ðx; tÞ rn ðxn ; tÞi the n points correlation function, the equation it satisfies is similar to (31) and reads q B2 rn ðx1 ; . . . ; xn ; tÞ ¼ A þ rn ðx1 ; . . . ; xn ; tÞ ; qt 2 A¼
n X q aðxi ; tÞ; qxi i¼1
B¼
n X q bðxi ; tÞ qxi i¼1
(32)
with initial conditions rn ðx1 ; . . . ; xn ; 0Þ ¼ rin ðx1 Þ rin ðxn Þ:
7. The probability distribution The process associated to the fluctuating density rðx; t; xÞ can be determined from a probabilistic viewpoint by associating a probability density Pðx; r; tÞ: At each time t indeed our field is specified by a space coordinate x and the value r it takes there. We may also introduce a conditional probability Pðx; r; t j x0 ; r0 ; t0 Þ equal to the weighted sum of all the paths joining in an extended phase space the initial point ðx0 ; r0 ; t0 Þ with the final point ðx; r; tÞ: The phase space ðx; rÞ 2 Rdþ1 has one extra dimension with respect to the ordinary phase space x 2 Rd where the Langevin equation is defined. The initial probability distribution is P0 ðx; rÞ ¼ Pðx; r; t0 Þ; not to be confused with the initial density rin ðxÞ: If the initial density is rin ðxÞ with probability 1 then P0 ðx; rÞ ¼ dðr rin ðxÞÞ : The normalization condition reads for any tXt0 Z þ1 Z þ1 dr dx rPðx; r; tÞ ¼ 1 : 0
1
(33)
(34)
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
The relation between the initial probability and its value at time t is Z þ1 Z þ1 Pðx; r; tÞ ¼ dr0 dx0 Pðx; r; t j x0 ; r0 ; t0 ÞP0 ðx0 ; r0 Þ ; 0
29
(35)
1
where the transition probability satisfies the standard normalization condition Z þ1 Z þ1 dr dx Pðx; r; t j x0 ; r0 ; t0 Þ ¼ 1 0
1
and the group property Z þ1 Z þ1 dr1 dx1 Pðx; r; t j x1 ; r1 ; t1 ÞPðx1 ; r1 ; t1 j x0 ; r0 ; t0 Þ 0
1
¼ Pðx; r; t j x0 ; r0 ; t0 Þ :
ð36Þ
To determine the transition probability we have to consider the characteristics of the SCE in the extended phase space ðx; rÞ: We assume for simplicity that qb=qx ¼ 0 so that the divergence of the vector field f is qa=qx ¼ 0: In this case the evolution equations read dx dr qa ¼ aðx; tÞ þ xðtÞbðx; tÞ; ¼ r : (37) dt dt qx We consider the continuity equation in the extended phase space for the fluctuating density Pðx; r; t; xÞ which satisfies the equation qP ¼ ðA þ Ar þ xðtÞBÞP ; qt where q qa q q aðx; tÞ; Ar ¼ r; B ¼ bðx; tÞ : qx qx qr qx Denoting by P¯ the average of P with respect to x; ¯ qP 1 2 ¯ ¼ A þ Ar þ Bn P : qt 2 A¼
(38)
(39)
(40)
¯ In the volume preservingR case Ar ¼ 0 and consequently the equation satisfied by P n¯ n and its moments hr i ¼ r P dr with respect to r is the Fokker–Planck equation (12). As a consequence we may identify P¯ with the probability density P in extended phase space. In the volume non preserving case it is not difficult to see that the equation satisfied by the moments hrn i is q n qa n 1 2 hr i þ hr i ¼ An þ B hrn i ; (41) qt qx 2 where An ¼ A ðn 1Þqa=qx is the operator defined by (15) and BRn ¼ B having ¯ dr ¼ assumed that qb=qx ¼ 0: This result easily follows by observing that rn Br P n n n qa=qx nhr i after an integration by parts. The relation with hr i is simple only if qa=qx is constant. In the one-dimensional case letting a be this scalar constant it
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
30
turns out that hrn i ¼ eaðtt0 Þ hrn i so that the probability density in extended phase ¯ space is given by P ¼ eaðtt0 Þ P: 7.1. The free particle with damping In order to clarify the above discussion we analyze in detail the case of the free particle with damping in which Eq. (37) reduces to x_ ¼ ax þ xðtÞ;
r_ ¼ ar
(42)
¯ is and the equation for P 2¯ q ¯ q ¯ a q ðrPÞ ¯ þ1 q P : ðxPÞ P¼a qt qx qr 2 qx2
(43)
The equation is solved by separation of variables and we write its fundamental solutions as ¯ r; t j x0 ; r0 ; t0 Þ ¼ Gðx; t j x0 ; t0 Þdðr r0 eaðtt0 Þ Þ ; Pðx;
(44)
where Gðx; t j x0 ; t0 Þ is the fundamental solution of the Fokker–Planck equation for the Uhlenbeck process 1 ðx x0 eaðtt0 Þ Þ2 Gðx; t j x0 ; t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi exp ; 2s2 ðt t0 Þ 2ps2 ðt t0 Þ s2 ðtÞ ¼
1 e2at : 2a
ð45Þ
¯ the probability density in the extended As a consequence recalling that P ¼ eaðtt0 Þ P phase space at time t is given by Z þ1 Z þ1 aðtt0 Þ dr0 dx0 Gðx; t j x0 ; t0 Þ Pðx; r; tÞ ¼ e 0
1
dðr r0 eaðtt0 Þ ÞP0 ðx0 ; r0 Þ :
ð46Þ
For an initial probability distribution with support on a unique distribution rin ðxÞ we have P0 ðx; rÞ ¼ dðr rin ðxÞÞ: In Fig. 6 we plot Pðx; r; tÞ as a function of r having chosen x ¼ 1 and a ¼ 1 and rin ðxÞ equal to the Gaussian distribution (17). The moments of this probability distribution are easily computed Z þ1 Z þ1 Z þ1 rn dr dr0 dx0 hrn i ¼ eaðtt0 Þ 0
0
1 aðtt0 Þ
Gðx; t j x0 ; t0 Þdðr r0 e ÞP0 ðx0 ; r0 Þ Z þ1 Z þ1 dr0 enaðtt0 Þ rn0 dx0 Gðx; t j x0 ; t0 Þdðr0 rin ðx0 ÞÞ ¼ eaðtt0 Þ 0 1 Z þ1 dx0 Gðx; t j x0 ; t0 Þrnin ðx0 Þ : ¼ eðn1Þaðtt0 Þ 1
ð47Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
31
Fig. 6. Plot of the probability density Pðx; r; tÞ for the free particle with damping ða ¼ 1Þ and white noise as a function of r at x ¼ 1 and t ¼ 0:501; 0.51, 0.6, and 0.8 for an initial probability density at t ¼ 0:5 concentrated on a Gaussian. The circles correspond to the numerical result while the lines correspond to the analytical result.
In Appendix C we carry out the explicit calculation of the second moment hr2 i using (47) also in the case where the white noise xðtÞ is replaced by a correlated noise.
8. Conclusions The SCE associated to a Langevin equation is a more general process in which a given initial density function is propagated by any realization of the process. The probability density associated to the Langevin equation is the average of the density function r computed over all the realizations. Naturally, the average of the density function does not exhaust all the information, which is captured instead by the average of all the integer powers of the density. This density evolves following the Fokker–Planck equation associated to the Langevin equation for a volume preserving vector field. The averages are equivalent to the moments of a probability density in the phase space whose dimension is increased by one, by including the r coordinate. This probability also satisfies a Fokker–Planck equation associated to a new Langevin equation in the extended phase space, describing the characteristics of the continuity equation. The meaning of the process corresponding to the SCE has been illustrated in the case of the free particle where all the computations are easily
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
32
carried out. From the numerical viewpoint the calculation of the probability density, given initially in analytic form, is fast and accurate at a given space point, since one can average the value there of the fluctuating density. The standard numerical procedure requires in any case the propagation of all the points of the initial probability density, which is usually done with a BMD procedure computing for each point a stochastic trajectory, which is an autonomous realization. We did not analyze the case of a correlated noise, since the extension is rather straightforward, at least conceptually. One can treat the correlated noise as an additional variable, by including the white noise-driven differential equation into the previous set of differential equations. It is possible to avoid the increase of the phase space dimensionality; however, the probability densities no longer satisfy an equation in closed form and the approximation schemes required are applicable only to noise of small amplitude and with rapidly decaying correlations.
Acknowledgements We wish to thank Prof. Ellison for useful discussions. A. Bazzani and G. Turchetti wish to thank DESY for hospitality during the preparation of this work, which was completed within the framework of the COFIN project 2003 Order and chaos in nonlinear extended systems: coherent structures, weak stochasticity and anomalous transport. G. Bassi wishes to thank DESY for a doctoral fellowship and the US Department of Energy for a partial support by Contract DE-FG03-99ER41104.
Appendix A. Direct computation of rð1Þ ðx; tÞ This is achieved by expanding the exponential (18) and performing the averages by using (B.6) 1 1 X ð1Þn hðx wðtÞ wðt0 ÞÞ2n i rð1Þ ¼ pffiffiffiffiffiffiffiffiffi n! 2pt0 n¼0 ð2t0 Þn 1 n 2n 1 X ð1Þn X x2ðnkÞ k ðt t ¼ pffiffiffiffiffiffiffiffiffi Þ ð2k 1Þ!! 0 2k 2pt0 n¼0 ð2t0 Þn k¼0 n! 1 1 X 1 X ð1Þkþm k ð2k 1Þ!! 2m þ 2k ¼ pffiffiffiffiffiffiffiffiffi x2m ðt t Þ 0 kþm ðm þ kÞ! 2k 2pt0 m¼0 k¼0 ð2t0 Þ 1 1 X 1 X x2m t0 t k 1 ð2k 1Þ!! ð2m þ 2kÞ! m! ¼ pffiffiffi ð1Þm t0 ð2kÞ!ð2mÞ! p m¼0 2k ðm þ kÞ! ð2t0 Þmþ1=2 m! k¼0 1 1 X 1 X x2m t0 t k 1 ð2m þ 2k 1Þ!! 2m m! ¼ pffiffiffi ð1Þm t0 k! ð2mÞ! p m¼0 2k ð2t0 Þmþ1=2 m! k¼0
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
1 1 X 1 X x2m t0 t k 1 ð2m þ 2k 1Þ!! ¼ pffiffiffi ð1Þm t0 p m¼0 2k ð2m 1Þ!!k! ð2t0 Þmþ1=2 m! k¼0 1 1 X x2m t0 t m1=2 m ¼ pffiffiffi 1 ð1Þ t0 p m¼0 ð2t0 Þmþ1=2 m! 2 1 x ¼ pffiffiffiffiffiffiffi exp : 2t 2pt
33
ðA:1Þ
In summing the series we have used the following identity: 1 X
an
n¼0
n X
bnk ck ¼
1 X
bm
an cnm ;
(A.2)
n¼m
m¼0
k¼0
1 X
where we have set m ¼ n k and interchanged the summation order. From k ¼ n mX0 it follows mpno1: Changing again the summation index in the second sum from n to k ¼ n m it follows now that 0pko1 and consequently we have 1 X n¼0
an
n X
bnk ck ¼
1 X m¼0
k¼0
bm
1 X
akþm ck
(A.3)
k¼0
We have used the following formula for the Taylor expansion of an algebraic function: ð1 xÞm1=2 ¼
1 X k¼0
xk
1 1 ð2m þ 2k 1Þ!! : 2k k! ð2m 1Þ!!
(A.4)
The series converges within jyjo1; namely our previous series converges for
t0 t
(A.5)
t o1; 0oto2t0 : 0 Obviously for t42t0 the result still holds with an analytic continuation performed after summing the series.
Appendix B. Wiener and Uhlenbeck processes The time evolution of a stochastic process is determined by a deterministic field and its fluctuations (noise). The white noise xðtÞ; formal time derivative of a Wiener process wðtÞ; describes a fluctuating field, whose correlations have an instantaneous decay hxðtÞi ¼ 0;
hxðtÞxðt0 Þi ¼ dðt t0 Þ :
(B.1)
The presence of memory corresponds to a non instantaneous decay of correlations and the decay rate characterizes the fluctuating field. The fluctuating part wðtÞ of the
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
34
process xðtÞ defined by dx ¼ axðtÞ þ xðtÞ dt
(B.2)
is asymptotic to the Ornstein–Ulhenbeck process for large t and reduces to the Wiener process wðtÞ for a ¼ 0: As a consequence letting xðt0 Þ ¼ x0 be the initial condition we have Z t aðtt0 Þ x0 þ wðtÞ; wðtÞ ¼ eaðtsÞ xðsÞ ds (B.3) xðtÞ ¼ e t0
and for a ¼ 0 we have Z
t
wðtÞ ¼ wðtÞ wðt0 Þ
xðsÞ ds
(B.4)
t0
By using the following property of the white noise correlations: X hxðs1 Þxðs2 Þ xðs2n1 Þxðs2n Þi ¼ dðsi1 si2 Þ dðsi2n1 si2n Þ ;
(B.5)
i1 ...i2n
where the sum runs over all the ð2n 1Þ!! pairings of the indices, it is not difficult to show that the following relations hold for tXt0 : n 1 e2aðtt0 Þ hw2n ðtÞi ¼ ð2n 1Þ!! 2a hðwðtÞ wðt0 ÞÞ2n i ¼ ðt t0 Þn ð2n 1Þ!! :
ðB:6Þ
The most general Langevin equation can be written as dx ¼ Fðx; t; xÞ aðx; tÞ þ bðx; tÞxðtÞ ; dt
(B.7)
where the stochastic field may be regularized by replacing xðtÞ with awðtÞ; a regular process whose limit is the white noise when a ! 1:
Appendix C. Space correlations for a free and damped particle In order to illustrate the computation of space correlations, we consider the free particle case in which Eq. (31) explicitly reads q 1 q q 2 r2 ðx; y; tÞ ¼ þ r2 ðx; y; tÞ : (C.1) qt 2 qx qy In order to find the solution corresponding to rin ðxÞ ¼ dðxÞ we perform a Fourier decomposition Z 2 1 r2 ðx; ty; tÞ ¼ ð2pÞ2 e2ðkx þky Þ tikx xiky y dkx dky : (C.2) R2
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
35
It is evident that pffiffiffi this is a solution pofffiffiffi Eq. (C.1). Rotating the coordinates with u ¼ ðkx þ ky Þ= 2 and v ¼ ðkx ky Þ= 2 we obtain xþy 2 1 Z Z xþy xy eð 2 Þ 2t u2 tiu pffi i pffi v 2 2 2 r2 ðx; y; tÞ ¼ ð2pÞ e du e dv ¼ pffiffiffiffiffiffiffi dðx yÞ : (C.3) 2pt R R 2
In the initial distribution is a Gaussian r0 ðxÞ ¼ ð2pt0 Þ1=2 expð 2tx Þ we may calculate 0 the space correlation as rð2Þ ¼ hrin ðx wðtÞ þ wðt0 ÞÞrin ðy wðtÞ þ wðt0 ÞÞi rather than solving (C.1) Z 1 ð2Þ rin ðx wÞrin ðy wÞpðw; t t0 Þ dw ; (C.4) r ðx; y; tÞ ¼ 1
where 2 1 w pðw; t t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2ðtt0 Þ : 2pðt t0 Þ
(C.5)
It follows that Z hrðx; tÞrðy; tÞi ¼
1
2 ðxwÞ2 ðywÞ2 1 1 1 w dw pffiffiffiffiffiffiffiffiffi e 2t0 pffiffiffiffiffiffiffiffiffi e 2t0 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2ðtt0 Þ 2pt0 2pt0 2pðt t0 Þ 1
tðx2 þy2 Þ2ðtt0 Þxy 1 ð Þ 2t0 ð2tt0 Þ pffiffiffiffipffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2p t0 2t t0 xþy 2 1 ðxyÞ2 1 1 ð Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2 2tt0 pffiffiffiffiffiffiffiffiffi e 4t0 : 4pt0 pð2t t0 Þ
¼
ðC:6Þ
Note that in the limit t0 ! 0 the new result (C.6) agrees with the previous one (C.3). Moreover, for y ¼ x we recover rð2Þ ðx; tÞ hr2 ðx; tÞi given by (23). C.1. Damped particle with white noise 2
x If the initial distribution is a Gaussian rin ðxÞ ¼ ð2ps20 Þ1=2 expð 2s 2 Þ we calculate 0
the space correlation for x_ ¼ ax þ xðtÞ which is given by (C.6) where 2 1 w pðw; t t0 Þ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi e 2sðtt0 Þ ; 2psðt t0 Þ
sðtÞ ¼
1 ð1 e2at Þ : 2a
(C.7)
It follows that hrðx; tÞrðy; tÞi ¼
e2aðtt0 Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2ps20 2ps2 ðt t0 Þ
Z
1
dw e
12 ðe2aðtt0 Þ ½ðxwÞ2 þðywÞ2 Þ 2s 0
e
w2 2s2 ðtt0 Þ
1
aðtt0 Þ
¼
e qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 ðx yÞ2 s2 ðt t0 Þe2aðtt0 Þ þ ðx2 þ y2 Þs20 exp : 2s20 ð2s2 ðt t0 Þ þ e2aðtt0 Þ s20 Þ
ðC:8Þ
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
36
For x ¼ y we obtain hr2 ðx; tÞi ¼
eaðtt0 Þ qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 x2 exp 2 2s ðt t0 Þ þ e2aðtt0 Þ s20
ðC:9Þ
in agreement with the result (D.4) obtained in Appendix D.
Appendix D. The second moment hr2 i for a damped particle In order to give a simple explicit example we consider a point subjected to a damping force and random uncorrelated collisions. The corresponding process for the point velocity, denoted by x; is defined by Eq. (B.2). The average density hri and the average of its square hr2 i satisfy the equations defined by (12) where A; B are replaced by An ; Bn defined by (15) for n ¼ 1; 2: These equations explicitly read q q 1 q2 hri ¼ a ðxhriÞ þ hri qt qx 2 qx2
(D.1)
and q 2 q 1 q2 2 hr i ¼ a ðxhr2 iÞ þ ahr2 i þ hr i ; (D.2) qt qx 2 qx2 which differs from the previous one by the divergence of the deterministic field. In this case there is no Langevin equation associated to the process defined by the fluctuating field r2 and in order to compute the mean hr2 i we evaluate the second moment with respect to r of the probability density Pðx; r; tÞ according to (47). Choosing an initial Gaussian density with variance s20 1 x2 rin ðxÞ ¼ qffiffiffiffiffiffiffiffiffiffi exp 2 (D.3) 2s0 2ps20 and inserting into (46) we obtain Z þ1 eaðtt0 Þ ðx x0 eaðtt0 Þ Þ2 x2 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p dx0 exp hr i ¼ exp 2 2s2 ðt t0 Þ 2s0 2ps20 2ps2 ðt t0 Þ 1 ¼
eaðtt0 Þ qffiffiffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2p s20 2s2 ðt t0 Þ þ e2aðtt0 Þ s20 x2 exp 2 : 2s ðt t0 Þ þ e2aðtt0 Þ s20
ðD:4Þ
If we choose s20 ¼ s2 ðt0 Þ then the denominators in the exponent and the argument of the second square root becomes s2 ðt t0 Þ þ s2 ðtÞ since the following identity holds
ARTICLE IN PRESS G. Bassi et al. / Physica A 347 (2005) 17–37
37
s2 ðt t0 Þ þ e2aðtt0 Þ s2 ðt0 Þ ¼ s2 ðtÞ: In the limit a ! 0 the result (30) for the second moment of the free particle is recovered. For the moment hrn i the result is obtained in a similar way. If the white noise is replaced by a Uhlenbeck noise according to dx dy ¼ ax þ yðtÞ; ¼ byðtÞ þ bxðtÞ (D.5) dt dt choosing yð0Þ ¼ 0 so that hyðtÞi ¼ 0 we have a new process, which tends to the previous one in the limit b ! 1: The process xðtÞ is Gaussian and writing the solution of (D.5) as Z t Z t at aðttÞ xðtÞ ¼ x0 e þ e yðtÞ dt yðtÞ ¼ b dtebðttÞ xðtÞ (D.6) 0
0 at
Rt the variance is easily computed after writing x x0 e ¼ 0 Kðt; sÞxðsÞ ds according to Z t s2 ðtÞ ¼ K 2 ðt; sÞ ds 0 2
b ð1 e2at Þ 2b3 ð1 eðaþbÞt Þ b3 ð1 e2bt Þ þ 2 ¼ : ðD:7Þ 2aðb aÞ2 ðb a2 Þðb2 abÞ 2ðb2 abÞ2 The evolution for r is still the same dr=dt ¼ ar and the probability density P is still ¯ is given by (44) and in (45) s2 ðtÞ is given by expressed by P ¼ eaðtt0 Þ P¯ where P (D.7). As a consequence the second moment hr2 i is still given by (D.4) where the variance is defined by (D.7). References [1] R. Kubo, in: D. Ter Haar (Ed.), Fluctuation, Relaxation and Resonance in Magnetic Systems, Oliver and Boyd, Edinburgh, 1962. [2] R. Kubo, Stochastic Liouville equations, J. Math. Phys. 4 (1963) 174–183. [3] R. Graham, Covariant formulation of non equilibrium statistical thermodynamics, Z. Phys. 26 (1977) 397–405. [4] J. Ellison, Accelerators and probability: the special effect of noise in beam dynamics, nonlinear and stochastic beam dynamics in accelerators—a challenge to theoretical and computational physics, DESY Proceedings, 1998. [5] A. Bazzani, G. Bassi, G. Turchetti, Diffusion and memory effects for stochastic processes and fractional Langevin equations, Physica A 324 (2003) 530–550. [6] F. Dyson, The radiation theories of Tomonaga, Schwinger and Feynmann, Phys. Rev. 75 (1949) 486–502. [7] R. Kubo, M. Toda, N. Hashitsume, Statistical Physics II, Non equilibrium Statistical Mechanics, Springer Verlag Series in Solid State Science 31 (1978) 25. [8] G. Bassi, A. Bazzani, H. Mais, G. Turchetti, Hamiltonian dynamics with weak noise and the echo effect, Physica A, submitted for publication.