Stochastic Molecular Dynamics of colloidal particles

Stochastic Molecular Dynamics of colloidal particles

Claudio Scherer

Instituto de Física, Universidade Federal do Rio Grande do Sul, CEP 91501-970, Porto Alegre, RS, Brazil

ABSTRACT

Colloidal particles move in the carrier liquid under the action of several forces and torques. When the particles carry a dipole moment, electric or magnetic, as in ferrofluids, the rotational and translational motions are coupled because the field on a particle depends on the spatial and directional distribution of the others and the force and torque on it depends on the field. Moreover, there is Brownian, as well as dissipative forces and torques on each particle. Consequently, the numerical solution of the equations of motion requires, besides the techniques of Classical Molecular Dynamics, those of Stochastic Dynamics. The algorithm is explained in some detail and applied on a typical ferrofluid. For different values of the temperature, the possibility of the formation of structures is examined.

1 Introduction

Colloidal particles have translational and rotational motion, due to several forces and torques. Besides forces and torques due to applied fields and inter-particle interactions, the molecules of the carrier liquid collide incessantly with the particle, causing translational and rotational Brownian motion. Their movement inside the liquid is opposed by viscous, dissipative, forces and torques. For simplicity, we will consider only spherical particles. We also assume that each particle carries a permanent dipole moment, which may be electric or magnetic. By using the appropriate Langevin equations we simulate realizations of the stochastic process which is the coupled motion of a sample of particles during a time which is sufficiently long for the thermodynamic equilibrium to be established. We use in the simulation typical values for the parameters, which correspond to realistic ferrofluids[1]. For recent references on ferrofluids, see the book Ferrofluids: Magnetically Controllable Fluids and Their Applications, edited by Stefan Odenbach[2]. Our procedure and results of the simulations are compared with those of Wang, Holm and Müller[3], recently published.

2 The Equations of Motion

The equations of motion for each particle can be written as one vector equation for its rotation around its center of mass and one vector equation for the translation of the center of mass. There is formally only one difference in the rotational equations of motion of colloidal particles having magnetic or electric dipole moments. This difference is in the existence of an intrinsic angular momentum, S, associated to the magnetic moment m ,

where g is the gyromagnetic factor. There is no intrinsic angular momentum associated with electric dipoles. The rotational equation of motion for a solid particle is given by Classical Mechanics,

where J is the total angular momentum and N is the total torque. The magnetic particles of ferrofluids may be superparamagnetic, when the magnetic moment is free to rotate with respect to the particle, or "blocked", when the magnetic moment is fixed in the particle. The case of superparamagnetic particles in ferrofluids has a more complicated dynamics because the magnetic moment and the particle's Euler angles are independent, but interacting, variables[4, 5, 6]. In this work we consider only blocked magnetic particles. For spherical magnetic particles the angular momentum may be written as

where I is the moment of inertia and w is the angular velocity. For electric dipolar particle the relation is the same, except that S =0. The time derivative of m is related to w , for blocked magnetic or electric dipoles, by

The torque N has several origins. If there is a field F, which may be the magnetic induction B or the electric field E, at the particle's position, there is a torque

N_F = m × F.

The field F is the sum of the applied field with the fields due to the other particles. We assume that the carrier liquid has no magnetic or electric moments. The thermal molecular motion of the liquid causes a stochastic torque, s x(t), on the particle, where s is a constant and x(t) we simulate by normalized white noise, with average zero and delta type correlation function,

where _i,j indicate the Cartesian components.

The rotational motion is opposed by a dissipative torque, –l w , due to the liquid's viscosity h, through the Stokes relation for rotation,

l = phd³

where d is the diameter of the spherical particle. Einstein's relation for the constants s and l and the temperature T reads

where k is Boltzmann's constant.

Summing up, we come to the following equation for the rotational motion,

which, together with Eq.(4), form, if F is known, a complete set of equations for the vectors w and m. The same equations describe also electric dipolar particles, except that the second term at the LHS of Eq.(8) is absent. However, when particle-particle interaction is important, as we want to consider in this work, F is dependent on the particle position r as well as on the other particles positions and moments. Consequently, we have to solve simultaneously the rotational and translational equations of motion for all particles.

The translational motion is described by Newton's equations,

for the position vector r and velocity v. The force f on the particle, like the torque, has several origins. If there is a field gradient ÑF at the particle's position, then the force due to this field is

The particle-particle interactions are partly due to their contribution to the field F, but we also assume a hard core interaction which avoids two particles to come closer than a distance d between their centers. This is not contained in the expression for f, but will be introduced directly in the integration procedure. The stochastic force, due to the collisions of the liquid's molecules with the particle, has the form aG(t), where a is a constant and G(t) is also modeled by normalized white noise, like in Eqs. (5) and (6). There is also a dissipative force –gv opposing the translational motion. The dissipative constant g is related to the viscosity by the following Stokes relation,

Einstein's relation for the constants a and g and the temperature T reads

Summing up, we arrive at the following equation for the translational motion:

which, together with Eq.(9), form a complete set of equations for the variables r and v if F(r, t) and m are known. However the particle-particle interactions make F to depend on the magnetic moments and positions of all particles, so that simultaneous solutions of all equations of motion is necessary.

3 Numerical Solution of Langevin Equations

Before discussing the procedures to solve the equations of motion of section 2, we present a brief introduction to numerical solutions of stochastic differential equations with white noise terms, known as Langevin equations.

Consider an n-dimensional stochastic process X(t), described by the following general Langevin equation,

where A(X,t) is a well behaved n-dimensional vector function, B is an n × m matrix and m is the number of independent components of the normalized white noise x(t). If B does not depend on X we say the noise is additive, otherwise it is called multiplicative. For our present purpose it is not necessary to allow for explicit dependencies of A and B on t, so that A = A(X) and B = B(X)

It is impossible to simulate, in the computer, realizations of the white noise, since its correlation time is zero. Therefore we integrate formally Eq.(15) between t and t + Dt,

Since A is a continuous function of X, which is a continuous function of t, we know, from the mean value theorem, that there is at least one value of t¢ such that the first integral above is

For the second integral at the RHS of Eq.(16) we consider first the case of additive noise. Then B may be taken outside the integral, so that the equation becomes

The integral above is known as Wiener increment,

where

is the Wiener process.

The components of the Wiener increment DW(t) are Gaussian stochastic processes with zero averages and standard deviations equal . These properties of DW are fundamental for simulating realizations of the solution X(t) of Eq.(15).

Eq.(18) has then a form very appropriate for numerical simulation,

In numerical simulation Dt, in Eq.(21), is the length of the time step for integration. If we choose a very small Dt, X() may be substituted simply by X(t) , or we may prefer a more precise algorithm, like Runge-Kutta second order. The components of DW are generated, at each time step, as the product of by random Gaussian numbers with zero average and unit variance.

The case of multiplicative noise is more complex. The rigorous treatments of Ito and of Stratonovich give us the prescriptions to follow, which are equivalent in the limit Dt ® 0, but for finite Dt they have different speed of convergence to the exact result. According to my own experience, by treating several examples, the Stratonovich procedure converges more rapidly than that of Ito, but, in many cases, it is more difficult to implement. We use here the Stratonovich procedure, according to which the equation equivalent to Eq.(21) for multiplicative noise should be written as

When we can isolate DX(t) in Eq.(22), very good, otherwise we may, for example, use simply B(X(t)) to obtain a first approximation for DX(t) and then substitute this value in the RHS of the Equation.

4 Numerical Simulation on Ferrofluids

Initially we consider the orders of magnitude of the terms in Eqs. (8) and (14). We start with the simplest one, which is Eq.(14). Applying the procedure described in section 3 to the system of Eqs. (9) and (14), follows

and

An alternative approximation to those equations consists in neglecting the inertia term in Eq.(14), which leads to

The numerical simulation of the system of Eqs. (23) and (24) and of Eq.(25) was made for the same realization of W(t), with ÑF = 0. As parameters we used those of a magnetite spherical particle of diameter d=10 nm in a liquid whose viscosity is that of water, at T=300 K. The results for one component of r(t) are shown in Fig. 1. We see that the effect of taking into account the mass is, indeed, negligible. This conclusion becomes even more evident in Fig. 2, where we used for the mass a value 100 times bigger, keeping the other parameters unchanged.

The advantage of neglecting the mass is that the convergence to the limit Dt ® 0 is much faster than when the mass term is present, which makes an important difference in CPU time when we simulate a system of many interacting particles.

Now we consider the equation of rotational motion, Eq.(8), using three alternative approximations: 1) neglecting the moment of inertia, I; 2) neglecting the intrinsic angular momentum, S = m/g; 3) neglecting both, I and S.

1) Neglecting I:

Let us define the following symbols:

s = m/m₀, where m₀ = |m| is a constant for blocked magnetic moment; DW_^ = DW – s (s · DW) is the component of DW perpendicular to s. After some vector algebra we get, neglecting terms of order higher than Dt²,

where

being s_x, s_y, s_z the Cartesian components of s, = g l/m₀ and = g s/m₀. Finally,

We improve the quality of the simulation by renormalizing s after each integration step, making |s| = 1.

2) Neglecting the intrinsic angular momentumS = m/g:

Again, after some algebra and neglecting terms of order higher than Dt², we come to the set of equations

and

The procedure is then simple: for given w(t), s(t), F(t) and DW we calculate Dw, and Ds in this order, and then w(t + Dt) = w(t) + Dw and s(t + Dt) = s(t) + Ds.

3) Neglecting both, I and S:

In this case Eq.(8) becomes

It takes again some vector algebra and the use of equation

to transform Eq.(32) to the form

where

The corresponding discretized Langevin equation is

In Fig. 3 we show s_z(t) obtained with the three approaches just described, for the same realization of W(t), using the parameters of a realistic magnetic particle in a ferrofluid: a spherical magnetite particle with diameter of 10 nm in a liquid with the viscosity of water at room temperature and in presence of a weak magnetic field of 10 G.

We see that taking into account the moment of inertia, I, or the intrinsic angular momentum, S, has practically no effect on the result of the simulation. To enlarge the effect of those terms we repeat the simulation with fictitious values of I and S 10 times bigger than those of the standard magnetite particle of 10 nm in diameter, keeping the other parameters unaltered. The result is shown in Fig. 4. In this case S has a considerable effect but I is still absolutely irrelevant, so that, if one of those terms should be taken into account, that should be S and not I. The opposite procedure was made in reference [3], as mentioned in our introduction.

5 The Inter-particle Forces and Torques

The forces between the magnetic particles in a ferrofluid are of two different origins: 1) A short range repulsion, which avoids two particles to overlap in space and 2) magnetic dipole-dipole force. The first we simulate by a hard spherical core with the size of the particle; the dipolar force on a given particle is calculated in the following way: we calculate the field gradient B, where B is the magnetic induction on the particle's position due to the other particles in the sample, then

where m_p is the particle's magnetic moment.

For a particle at position r_p, the magnetic induction due to the other particles is

where n is a unit vector in the direction of r_p – r_j.

The torque on m_p is

The numerical simulations on a random homogeneous distribution of magnetic dipoles have convinced us that the contribution to B(r_p) due to the dipoles which are more than a few inter-particle distances apart from r_p is totally negligible. This is so because when the spacial distribution of the moments is precisely homogeneous, even if there is a preferential direction for their orientation, the field at the center of a cube, due to the other particles in the cube, is zero. Therefore a non zero field comes only from deviations from homogeneity, which is really important only in short distances. Our simulation is done on a cubic box of side L of a ferrofluid containing 1000 magnetic particles. We use periodic boundary conditions. To calculate the induction B at r_p due to the particle at r_j, if the x-coordinate x_p – x_j > L/2 we substitute x_j by x_j + L in Eq.(38), and similarly for the other coordinates; if x_p – x_j < –L/2 we substitute x_j by x_j – L. With this procedure, the field on each particles considers all other particles which are in the cube of side L of which the considered particle is at the center, as is indicated in a two-dimensional projection in Fig. 5. To calculate the field on particle p in Fig. 5, due to all other particles inside the simulation box, in solid line, we use, in the case of particle 1, its actual coordinates but in the cases of particles 2, 3, 4 and 5, the coordinates of their periodic images, which are inside the box in doted line, of which p is the center.

After each integration step we translate each particle which went out of the simulation box like particle 2 in the figure, to inside the box, at the position of its periodic image.

We start the simulation by placing the magnetic particles at the sites of a cubic lattice inside the simulation box, the components of the magnetic moments are chosen at random and the modulus normalized to m₀. The fields and field gradients on the particles are calculated as explained above. For the integration we use the approach of the equations without mass, moment of inertia and intrinsic angular momentum, as explained in section 4. For realistic values of the parameters, as those used for Figs. 1 and 3, room temperature and zero applied field, after sufficiently long simulation time, when the distribution seem not to change significantly anymore, the resulting distribution is shown in Fig. 6 for a slice of the simulation box, parallel to the x-y plane and of width equal to the average inter-particle distance. In Fig. 6 the lines indicate the projection of the magnetic moments on the x-y plane. We see that nothing of the lines of aligned dipoles, as reported in reference [3], can be seen. On the other hand, simulation made on a fictitious ferrofluid at T=3K (not 300K) lead to Fig. 7, where lines similar to those reported in reference [3] are seen. This result makes us suspect that an error in the generation of the noise was made by the authors of reference [3], resulting in noise terms much weaker than what should be the case at 300 K.

Received on 2 September, 2003

• [1] R.E. Rosensweig, Ferrohydrodynamics, Dover Publications, Inc., New York, 1997.
• [2] S. Odenbach (Ed.), Ferrofluids: Magnetically Controllable Fluids and Their Applications, Springer-Verlag, Berlin (2002).
• [3] Z. Wang, C. Holm, and H.W. Müller, Phys. Rev. E66, 021405 (2002).
• [4] M.I. Shliomis and V.I. Stepanov, Adv. Chem. Phys. Series 87, 1 (1994).
• [5] C. Scherer and G. Matuttis, Phys. Rev. E63, 011504 (2001).
• [6] C. Scherer and T. Ricci, Braz. J. Phys. 31, 380 (2001).

Publication Dates

• Publication in this collection
  01 Sept 2004
• Date of issue
  June 2004

History

• Received
  02 Sept 2003
• Accepted
  02 Sept 2003