Stochastic mechanics of nonequilibrium systems

Stochastic Mechanics of Nonequilibrium Systems^*

Tânia Tomé and Mário J. de Oliveira

Instituto de Física, Universidade de São Paulo

Caixa Postal 66318, 05315-970, São Paulo, SP, Brazil

Received April 10, 1997

A nonequilibrium system of classical particles interacting through nonconservative forces is studied by a stochastic Markovian process defined over the space of the particle positions. The velocities of the particles are treated as independent stochocastic variables with a given probability distribution. We deduce expressions for the rate in which energy is exchanged with the surrounding as well as the rate in which energy is dissipated. These results are applied to the special case in which the irrotational and solenoidal parts of the forces acting on the particles are orthogonal. We also solve exactly a model in which these two forces are linear functions of positions.

I. Introduction

We consider here the stochastic Markovian approach in the study of nonequilibrium system of classical particles interacting through forces that may be nonconservative. According to this approach, the deterministic trajectory of the system in phase space is replaced by a random trajectory governed by a stochastic or Langevin equation [1] [2] [3]. Our starting point is a Langevin equation of motion defined over the space of the particle positions only. The probability distribution defined on this space will evolve in time according to a Smoluchowski equation which is a Fokker-Planck equation in position space [1] [2] [3]. Such a stochastic approach can be understood as a coarse grained description in which the velocities become stochastic variables independent of positions, with a certain given probability distribution. The positions are also stochastic variables but depend on the velocities.

The coarse grained description used here is obtained by assuming that within a characteristic time interval t the velocities of the particles become independent of the positions and acquire a distribution whose average speed of the particles is u. The description is then characterized by the velocity relaxation time t and by the average speed u of the particles after relaxation. We assume also that the forces among particles vary so slowly over a distance ut that may be considered constant within the time interval t.

Within the above assumptions we obtain the Langevin equation and the corresponding Fokker-Planck equation in position space. We consider systems described by nonconservative forces and deduce from these equations expressions for the rate in which energy is exchanged with the surrounding as well as the rate in which energy is dissipated. These results are applied to the special case in which the irrotational and solenoidal parts of the forces acting on the particles are orthogonal. We solve exactly a model in which these two orthogonal forces are linear functions of positions.

There are other ways to justify the position space description used here. For instance, one could start from a Kramers equation, which is a Fokker-Planck equation in position and velocity space, and perform the elimination of velocities treated as being fast variables [1] [2] [3]. But in this case we would have to introduce first the stochastic equation that describes the trajectory of the system in the position and velocity space and then proceed to the elimination of the fast variables. Moreover this procedure involves the use of a friction term that we avoided in our approach.

The point of view we use here in the study of nonequilibrium system described by continuous variables is the same one adopted by Tomé [4] in the study of nonequilibrium stochastic lattice gas models. In this case the analysis of lattice gas models by a stochastic dynamics [5] [6] [7] has been successful in describing the phase transitions in system far from equilibrium [7-15]. Other approaches, macroscopic and microscopic, in which fluctuations play an important role, have been used in explaining self-organization and dissipative structures. These include the pioneering work of Nicolis and Prigogine[16] and the one based on the nonequilibrium statistical operator method introduced by Zubarev [17]. This last formalism has been explored extensively and successfully by Luzzi and Vasconcellos [18] in the study of dissipative structures in condensed matter.

II. Stochastic equation of motion

Consider a classical system composed of a certain number of identical particles of mass m each one interacting with each other and possibly with the external environment. Let x_i(t) and v_i(t) be the position and the velocity of the i-th particle at time t and let f_i(x) be the total force on particle i due to the other particles which is supposed to depend only on the positions. The forces acting on particles may be conservative or nonconservative. The notation x stands for the collection {x_i} of positions of the particles and v will denote the set {v_i} of their velocities . Assuming that the forces f_i(x) can be considered constant during a time interval t, the position of the i-th particle at time t+t will be given by

where f_i(x(t))/m is the acceleration and v_i(t) is the velocity of the particle at the beginning of the time interval, and m is the mass of the particle.

According to our assumptions, the velocity v(t+t) of the particles at time t+t will not depend on the positions x(t) of the particles at time t but could depend on the velocity v(t) at time t. However, we assume that it will not depend on the velocity v(t) of the particles on the previous time so that < v_i(t+t)v_j(t) > = 0 for any pair i, j. The velocity probability distribution of distinct particles are identical and are characterized by having zero mean and dispersion u², that is,

< v_i(t) > = 0 and < v_i(t)v_j(t) > = u²d_ij . (2)

Let us define a = 2m/t and z_i(t) = av_i(t) so that equation (1) is written in the form of a stochastic difference equation

where the noise z_i(t) has the properties

where Q = mu² is proportional to the average kinetic energy of the particles

Notice that there are only three independent parameters: the relaxation time t, the average speed u and the mass m of the particles. The others, Q and a, are defined in terms of these three. In what follows however we will adopt as independent parameters t, Q and a so that m = ta/2 and . Formally we may take the limit t® 0, with Q and a constants to obtain the following stochastic differential equation or Langevin equation,

where z_i(t) is a white noise with the properties

< z_i(t) > = 0 , < z_i(t)z_j(t^¢) > = 2aQd_ijd(t-t^¢) . (6)

In this case we are left with only two parameters which are Q and a. In any case, we will always consider equations (5) and (6) as meaning equation (3) and (4) with small t. Eventually we may take the limit t® 0 in case it makes sense.

The Fokker-Planck equation corresponding to the Langevin equation (5) with properties (6), that is, the equation that governs the time evolution of the probability P(x,t) of the system being at position x at time t, is given by

which we call Smoluchowski equation.

By solving this equation we may obtain the average < F > of any state function F(x) by

< F > = ó
õ F(x)P(x,t)dx. (8)

Quantities that are not functions of positions cannot, in principle, be calculated from this formula. We shall see however that there are some quantities of this type that can be written as an average of a certain function of positions. One purpose of this paper is to show how this can be done for dissipation of energy in nonequilibrium systems.

In the case of a free particle (Brownian motion) where the forces vanish, the equation

with the properties (6) can be exactly solved and we get the following result

The Smoluchowski equation becomes the pure diffusion equation

The diffusion constant can also be written as D = tu²/2 = l²/2t where l = ut.

III. Equilibrium state

According to classical statistical mechanics the probability distribution P(x) of a system interacting through a potential energy E(x) and in equilibrium with a heat reservoir at an absolute temperature T is given by

where b = 1/k_BT and Q is a normalization factor. The potential energy E(x) is defined over the position space x = {x_i} composed of the positions of the particles. Here we show that if the force f_i(x) is conservative, that is, if f_i(x) is the gradient of a potential energy E(x), then (12) is the stationary solution of the Smoluchowski equation (7).

Consider the Langevin equation (5) with the noise having the properties (6). The corresponding Fokker-Planck equation (7) can be written in the form

where J_i(x,t) is the current defined by

The stationary solution P(x) satisfies the equation

where the stationary J_i(x) current is given by

If the force f_i(x) is the gradient of a potential energy E(x), that is, if

then one can write

from which one easily sees that the stationary solution is

where A is a normalization constant. This solution not only satisfies equation (15) but also satisfies J_i(x) = 0. That is, the current vanishes at each point of position space.

We conclude that the probability distribution (19) becomes identical to the probability distribution given by equation (12) of classical equilibrium system at temperature T such that

Q = k_BT . (20)

Hence, may say that the stationary state of a system interacting through conservative forces is an equilibrium state. To known a priori whether the stationary state of a given system is an equilibrium state we just have to check whether the condition

is satisfied for any pair (i,j). When this condition is satisfied the system has microscopic reversibility.

The results of this section serve as the basis of the numerical simulation method used in equilibrium statistical mechanics known as the Langevin method. Suppose we wish to estimate numerically averages of functions that depend only on the positions of a system described by the potencial energy E(x). If we were able to generate a set of M points x⁽¹⁾, x⁽²⁾, x⁽³⁾, ...x^(M) in position space, each one sampled according to the probability distribution given by equation (12) then an estimative of the average

< F > = ó
õ F(x)P(x)dx (22)

will be simply given by

The numerical problem is then reduced to devise methods of sampling points in the position space according to a probability distribution (12) known a priori. These methods are englobed in basically two types that are known as the Monte Carlo methods and the method of the Langevin equation. Both approaches use a stochastic Markov process defined over the position space whose stationary probability distribution is P(x). The stochastic process can be viewed as a trajectory in position space such that, for long times, the points of the trajectory will be emerging with the desired probability P(x). If we simulate equation (5), with f_i(x) = -¶E(x)/¶x_i, we will generate, for long times, a set of points that are distributed according to equation (12).

IV. Power of a conservative force

The average rate in which any given force F_i(x) does work, that is, its average power, is determined by

We want to find a state function W(x) such that

W = < W(x) > . (25)

To this end we consider first the case in which the force is conservative, that is, F_i(x) = -¶V(x)/¶x_i. In this case

so that the average power of a conservative force is the rate of the average potential energy < V(x) > . It is clear that in the stationary state < V(x) > does not depend on time and the power W vanishes.

Now, expression (26) gives

Using the equation (13) we obtain that

After an integration by parts and using the definition (16) of current and equation (17) we get

Performing once more an integration by parts in the last integral we further obtain that

and then

as desired.

V. Power of a generic force

Let us now consider the power developed by a generic force F_i(x) that might not be conservative. In this case we cannot use the procedure of the preceding section since the force cannot in general be written as the gradient of a certain function. We proceed then as follows. We start by using the discretize form

of the velocity dx_i/dt, where 0 £ q £ 1. The expression for the power of the force F_i is then

Using the discretize form (3) of the Langevin equation this expression can be written as

Now, in the limit t® 0,

< {qf_i(x(t))+(1-q)f_i(x(t-t))}F_i(x(t)) > ® < f_i(x)F_i(x) > (35)

Also

< z_i(t)F_i(x(t)) > = 0 , (36)

since z_i(t) is independent of x_i(t). Next we perform the following expansion

F_i(x(t+t)) = F_i(x(t))+ å
j F_ij(x(t))[x_j(t+t)-x_j(t)] , (37)

where F_ij(x) = ¶F_i/¶x_j, or, using again the discretize form (3) of the Langevin equation, we find that

From this expression we obtain that

< z_i(t)F_i(x(t+t)) > = 2Q < F_ii(x(t)) > (39)

when t® 0. Therefore in the limit t® 0 it follows that

Since this expression must be valid also for conservative forces we should choose q = 1/2 so that it becomes identical with expression (30).

For any force F_i(x), conservative or nonconservative, we have then

W = < W(x) > , (41)

where

The quantity vanishes in the case of a conservative force F_i no matter whether the force f_i(x) that describes the system is conservative or nonconservative.

VI. Exchanged and dissipated power

Let us consider a system described by a force that might not be conservative, that is, one such that equation (21) might not be satisfied. When the force f_i(x) is nonconservative the stationary state will not be an equilibrium state and we expect that energy will be dissipated continuously. From the result of the previous section, the power of the force f_i(x) will be

In the stationary state this expression will give the dissipated power since the exchanged power must vanish. However, if the system is in a transient state this expression gives the total power, that is, the sum of the exchanged and dissipated power. Let us try next to calculate each of these two contributions.

We start by splitting the force f_i(x) into two parts

f_i(x) = f_i^C(x)+f_i^D(x) , (44)

such that f_i^C(x) is conservative, that is, f_i^C(x) is the gradient of a potential E(x), or

and f_i^D(x) is a nonconservative force with a vanishing divergence (solenoidal force), that is,

That such a splitting can always be done can be seen as follows. This last equation implies that

or yet

Therefore, if we are given a generic force f_i(x), this last equation can be solved to obtain E(x). From E(x) we get f_i^C by means of equation (45) and f_i^D by using equation (44).

The total power is then splitted into two parts

W = W_C+W_D , (49)

where

is the power of force f_i^C(x) and

is the power of force f_i^D(x). These expression were obtained from the expression by setting f_i = f_i^C+f_i^D and F_i equal to f_i^C and f_i^D respectively. We have also used relation (46). Now W_C vanishes in the stationary state since it is the power of a conservative force and it is identified with the energy exchanged per unit time. The quantity W_D vanishes identically when the system is conservative (f_i^D = 0) and it is identified with the energy dissipated per unit time. It also must be positive. For the special case in which the two components of the force f_i^C(x) and f_i^D(x) are orthogonal this is indeed the case. A system described by a force whose components are of this type is considered in the next section.

Equation (49) may be interpreted by assuming, for instance, that the system is an open system in contact with two reservoir. One of them is the source of energy and the other is the recipient of the dissipated energy. In a given time interval Dt the first reservoir releases an amount of energy DE to the system. A part of it DE_C, is retained by the system whereas the remaining part DE_D is dissipated and released to the second reservoir. In the stationary state DE_C = 0 (W_C = 0) and all the energy released to the system is dissipated. If, moreover, the stationary state is an equilibrium state, then also DE_D = 0 (W_D = 0) and there is no exchange of energy between the source reservoir and the system.

VII. Orthogonal forces

Consider the special case where the irrotational forces f_i^C and the solenoidal forces f_i^D are orthogonal that is

å
i f_i^C(x)f_i^D(x) = 0 . (52)

In this case

so that the dissipated power is positive and

Although the system is described by force f_i = f_i^C+f_i^D which is not conservative since f_i^D is not conservative we show that in this case it is possible to find the nonequilibrium stationary state. To this end, consider the Fokker-Planck equation (7) in the stationary regime, that is,

If the conditions (52) is fulfilled and assuming that f_i^D(x) is orthogonal to the gradient of P(x), that is

then the second term in equation (55) vanishes and what is left is a Fokker-Planck equation for a system described by a conservative force f_i^C(x) = -¶E/¶x_i. The solution is therefore

We can check from this solution that the assumption (56) is indeed correct. Using (57) we can check explicitly that W_C vanishes.

VIII. Simple model

Let us consider a simple model of just two variables. The force is given by

where K and A are two parameters, and r and f are the two-dimensional vectors r = x₁+x₂ and f = f₁+f₂. As long as A ¹ 0 the force is nonconservative since Ñ×f = -2A¹ 0. In this case

and

The associate Fokker-Planck equation is

that is,

which can be written as

Although this equation does not have radial symmetry, it does possess a solution with radial symmetry. Indeed if we assume that P is a function of |r| = r only then the vector product term in (63) vanishes and the equation becomes a Fokker-Planck equation for a conservative force given by f = -Kr. The time dependent solution of such an equation is

where

For simplicity we used initial conditions such that < r > = 0.

To obtain W_C and W_D we insert the result < r² > = 2B(t), obtained from (64), into equations (59) and (60) to get

and

Therefore, W_C and W_D decay exponentially to their stationary values. Notice that W_D is always nonnegative since B(t) is nonnegative. The quantity W_C on the other hand can be positive or negative depending on the sign of B₀.

The stationary probability P(r) is given by

where the normalization factor is P₀ = K/2pQ, is obtained by taking the limit t ® ¥ in equation (64). The stationary values W_C and W_D are

W_C = 0 , (69)

and

so that the rate of dissipation of energy W_D is proportional to A² and, of course, it vanishes when A goes to zero, that is, in the equilibrium case.

We may also calculate the stationary current. From its definition

J(x) = f(x)P(x)-QÑP(x) , (71)

we obtain the result

Hence, the stationary current is proportional to A and also vanishes in the equilibrium state, as expected.

IX. Conclusion

We have studied a nonequilibrium classical system of interacting particles by a stochastic approach whose starting point is a Langevin equation of motion in position space. From this equation we have obtained expressions for the rate in which energy is exchanged as well as the rate in which energy is dissipated. The exchanged and dissipated power are related to the irrotational and solenoidal parts of the forces, respectively. We have applied the results to the case where these two parts are orthogonal and solved exactly a model in which the forces are linear functions of positions.

References

[1]
• N. G. van Kampen, Stochastic Processes in Physics and Chemistry, North-Holland, Amsterdam, 1981.
[2]
• C. W. Gardiner, Handbook of Stochastic Methods, Springer-Verlag, Berlin, 1983.
[3]
• H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications, Springer-Verlag, Berlin, 1984.
[4]
• T. Tomé, Tese de Livre-Docência, Universidade de São Paulo, São Paulo, 1996.
[5]
• T. M. Liggett, Interacting Particle Systems, Springer-Verlag, New York, 1985.
[6]
• N. Konno, Phase Transitions of Interacting Particle Systems, World Scientific, Singapore, 1994.
[7]
• J. Marro and R. Dickman, Nonequilibrium Phase Transitions, Cambridge University Press, Boston, 1996.
[8]
• T. Tomé and M. J. de Oliveira, Phys. Rev. A 40, 6643 (1989).
[9]
• T. Tomé, M. J. de Oliveira and M. A. Santos, J. Phys. A 24, 3677 (1991).
[10]
• M. J. Oliveira, J. Stat. Phys. 66, 273 (1992).
[11]
• M. J. de Oliveira, J. F. Ferreira and M. A. Santos, J. Phys. A 26, 2317 (1993).
[12]
• T. Tomé and R. Dickman, Phys. Rev. E 47, 948 (1993).
[13]
• M. J. de Oliveira, Physica A 203, 13 (1994).
[14]
• J. Satulovsky and T. Tomé, Phys. Rev. E 49, 5073 (1994).
[15]
• T. Tomé and J. R. Drugowich de Felício, Phys. Rev. E 53, 3976 (1996).
[16]
• G. Nicolis and I. Prigogine, Self-Organization in Nonequilibrium Systems, Wiley-Interscience, New York, 1977.
[17]
• D. N. Zubarev, Nonequilibrium Statistical Thermodynamics, Consultants Bureau, New York, 1974.
[18]
• R. Luzzi and A. R. Vasconcellos, Ciência e Cultura 43, 423 (1991).

*

This paper is dedicated to Prof. Roberto Luzzi in his 60th anniversary.

Publication Dates

• Publication in this collection
  14 June 1999
• Date of issue
  Dec 1997

History

• Received
  10 Apr 1997