Computer simulation of the stochastic dynamics of super-paramagnetic particles in ferrofluids

Computer simulation of the stochastic dynamics of super-paramagnetic particles in ferrofluids

Claudio Scherer

Institute of Physics, Federal University of Rio Grande do Sul, Porto Alegre, Brazil

ABSTRACT

Several papers have been written on the complex problem of the stochastic dynamics of the magnetic moments of super-paramagnetic particles, simultaneously with the stochastic rotation of these colloidal particles in a ferrofluid [1-3]. None of these works, however, is sufficiently general and conveniently simple and clear to be used in sumulational works to appropriately describe the experimental super-paramagntetic resonance lines. We have a new proposal for the equations of rotational motion, which is appropriate for simulations. Those equations are stochastic differential equations with multiplicative noise. Therefore, they have to be interpreted as Stratonovich-Langevin equations and the roles of stochastic calculus have to be used in the simulations. For this reason we will briefly present the essence of the numerical algorithms used in the solutions of Stratonovich equations. Finally, the simulational results for the magnetic response functions, and the corresponding dynamic susceptibilities, will be shown and their consequences will be analyzed.

Keywords: Ferrofluid, Magnetic fluid; Computer simulation; Stochastic dynamic; Super-paramagnetic

I. INTRODUCTION

The dynamics of the magnetic moments of super-paramagnetic particles constitute a nice example of problem in the Physics of condensed matter where to apply the tools of stochastic calculus to perform numerical simulation. If those particles are in a ferrofluid, then the problem becomes still more interesting, because one has a coupled rotational motion of the moments with respect to the particles and of the particles with respect to the fluid. Several physical variables may be calculated from the numerical solutions of the equations of motion for an ensemble of super-paramagnetic particles, which can then be compared to experimental results. We give particular attention to the linear response functions, from which the dynamical susceptibility is obtained by Fourier transform, and then the magnetic resonance absorption lines.

In section II we introduce some hypothesis, which are very realistic for real ferrofluids, defining in this way the essence of the model. The equations of motion are also introduced in this section. Since the simulation is done for an equilibrium ensemble of particles, the probability distributions of the relevant variables is discussed in section III. General properties of stochastic differential equations and numerical procedures to get solutions for their realizations are treated in section IV. Finally, in section V we solve the set of equations for the super-paramagnetic particles, for different values of the physical parameters, calculate the response functions and the corresponding dynamical susceptibilities and present the results in figures. Some consequences of our results are discussed.

II. THE MODEL AND THE EQUATIONS OF MOTION

Some hypothesis used in this work are:

1) The magnetic particle has a symmetry axis of easy magnetization, which will be characterized by the unit vector c;

2) The particle's magnetic moment µ has constant modulus µ and can rotate inside the particle under the influence of an interaction potential modeled by V = - K(µ·c)², where K is a constant, called änisotropy constant"; we will use units such that µ = 1;

3) The particle's rotational inertia is negligible in the equations of motion, in comparison with the Brownian and dissipative terms;

These three hypothesis are very realistic for the normal conditions of super-paramagnetic particles in ferrofluids above the blocking temperature.

The equations of motion for µ are obtained from the classical equation for the rotation of a system with angular momentum J in presence of a torque N,

The existence of a magnetic moment µ implies the existence of an angular momentum, such that µ = gJ. In the general case the "gyromagnetic factor" g is a tensor, but frequently it is simply a scalar. We will assume this to be the case and chose units such that g = 1. Therefore Eq.(1) may be written as

The torques on µ are:

a) The torque due to the interaction potential V of hypothesis 2 above is obtained by deriving V with respect to the angle between µ and c, which results in 2K (µ·c) (m×c);

b) The Brownian torque due to the atomic vibrations of thermal origin (phonons) of the particle: a_mm×x_m;

c) The dissipative torque of resistance to the Brownian rotation (fluctuation dissipation theorem): -g_mm× (throughout this paper means time derivative of µ);

d) The torque of interaction of µ with a magnetic field H, applied or due to the other particles in the ferrofluid, which happens to be present at the particles position: µ×H.

Therefore the equation of motion for µ, Eq.(2), may be written as

Substituting at the RHS of Eq.(3) by the whole of the RHS, and using the properties of the triple vector product we get

If we call H_eff = H + 2K (µ·c)c +a_mx_m, Eq.(4) becomes the equation of Landau-Lifshitz.

The only motion of the particle which will be considered is the rotation of the symmetry axis, c, since the rotation around c is irrelevant for the motion of µ. The last three terms of Eq.(3) are due to the interaction of µ with the particle and, therefore, they have to be present, with opposite sign, in the equation for c. Besides those terms there is the Brownian torque on the particle, due to the thermal motion of the liquid's molecules and the corresponding dissipative torque, according to the fluctuation-dissipation theorem. Since the rotational inertia of the particle will be neglected (hypothesis 3 above), the equation of motion for c reduces to a balance between the torques,

which may be written more conveniently as

Making the vector product of Eq.(6), at the left, by c and using the properties of the triple vector product, follows

where, in the last term, we substituted c ×(c ×x_l) by c ×x_l , because both expressions are white noise perpendicular to c and, therefore, have the same statistical properties and the last expression is simpler.

III. EQUILIBRIUM DISTRIBUTIONS

Before performing the simulations of the dynamics of µ and c we need to get the equilibrium distributions of these vectors, which are the initial distributions from which we will simulate the dynamics. The independent stochastic variables are the polar coordinates q and f of µ and the polar coordinates J and j of c.

The equilibrium distribution for these four independent variables, in presence of a magnetic field H in the z direction, is given by Boltzmann's distribution,

where µ·c = sin q sin J cos(f - j) + cos(q)cos(J), µ·H = H cos(q) and C is a normalization constant.

To obtain, the equilibrium distribution for one of the variables, for example Pº(q), for q, we perform the numerical integration of the other independent variables, f, J and j (see Fig.1, lines a and b).

For the simulation of the dynamics we will need to generate an equilibrium ensemble of particles with coordinates distributed according to Eq.(8). There are several ways of doing this [4, 5]. Fig.1, line c, shows Pº(q) obtained form the histogram for the q's of 10000 particles of an equilibrium ensemble, for different parameters, as described in the figure caption. When the histogram is done for 1000000 particles, as in the dynamic simulation described below, the theoretical and simulated curves become indistinguishable.

It is interesting to note that the equilibrium distribution of the orientation of µ is independent of the anisotropy K. The orientation of µ depends on the ratio H/T while the orientation of the particle, specified by c, depends on H, K and T, as can be seen in Fig.2.

We will see bellow that for the dynamic response of the magnetization this independence of the orientation of µ with respect to the anisotropy does not happen, i.e., the dynamic behavior of µ depends on the anisotropy.

IV. NUMERICAL SIMULATION OF STOCHASTIC DIFFERENTIAL EQUATIONS WITH MULTIPLICATE NOISE

Equations 4 and 7, taken together, may be written in the form of a Stratonovich stochastic differential equation [5, 6],

where X(t) is an n-component vector (in our case n=6, X = [µ, c] = [µ_x, µ_y, µ_z, c_x, c_y, c_z]), A(X) is an n-component vector (RHS of equations 4 and 7, except for the terms with the white-noise), dW(t) is an m-component "Wiener differential" (see reference [6]), which corresponds to the integral of the white-noise over the time interval (t, t+dt), and (X) is an n × m matrix. Each component dW_j of dW is a Gaussian stochastic process with zero mean and variance dt, i.e.,

In some cases is a constant matrix, i.e., does not depend on X, and then Eq.(9) is said to be of ädditive noise". Otherwise, namely when depends on X, the equation is said to be of "multiplicative noise". In the former case Eq.(9) may be solved by ördinary stochastic calculus", while in the last case one have to use "Stratonovich" or "Ito" calculus. This is the case of equations (4) and (7), since contains vector products with µ and c, which are components of X.

We write Eq.(9) in terms of "finite differences", according to Stratonovich calculus. Calling Dt, DX and DW the finite differences of t, X and W between two consecutive instants of the discretized time, we have

The term DX/2 in the argument essential for Stratonovich calculus, whereas it is optional, but convenient, in the argument of A.

To generate, in computer simulation, the independent components DW_j of the "Wiener increments" DW , we proceed in the following way: generate a Gaussian random number R_G, with zero mean and unit variance and write

since the width of the distribution for DW_j is .

It is in general, and in our case in particular, difficult, or even impossible, to isolate DX in Eq.(11). For this reason we will use an algorithm which we call "second order Stratonovich", as follows:

1) Calculate DX⁰ equal to the RHS of Eq.(11) without the terms DX/2 in the arguments of A and ;

2) Use DX⁰ instead of DX in the arguments of A and to obtain DX.

This procedure is analogous to Runge-Kutta 2^nd order in the case of ordinary differential equations. It gives the same level of precision as we would have by transforming Eq.(9) into an Ito differential equation and integrating it by Ito calculus, but it is simpler to implement.

V. SOLUTION OF EQUATIONS (4) AND (7) BY STRATONOVICH CALCULUS

Let us write equations (4) and (7) explicitly in the form of Stratonovich stochastic differential equations:

Now we write equations (13) and (14) in form of finite differences according to Stratonovich integration rule. For the time step (t, t+Dt) we use the following notation: µ = µ(t); c = c(t); and Dm and Dc the corresponding increments in the interval.

We obtain µ(t+Dt) and c(t+Dt) proceeding as follows, according to the 2^nd order Stratonovich algorithm, as described in section IV:

1) Calculate Dm⁰ from Eq.(15), by neglecting everywhere at the RHS the terms Dm/2 and Dc/2;

2) Calculate Dc⁰ from Eq.(16), by neglecting everywhere at the RHS the terms Dm/2 and Dc/2 and using Dm⁰, just calculated, instead of Dm in the term g_mDm;

3) Calculate Dm from Eq.(15), by using Dm⁰/2 and Dc⁰/2 everywhere at the RHS instead of Dm/2 and Dc/2;

4) Calculate Dc from Eq.(16), by using Dm⁰/2 and Dc⁰/2 everywhere at the RHS instead of Dm/2 and Dc/2 and using Dm, just calculated, in the term g_mDm;

5) Obtain µ(t+Dt) = µ+Dm and c(t+Dt) = c +Dc and correct the norms of those quantities from the small errors (2^nd order in Dt), which may be produced by the above calculations, by dividing both by their corresponding norms.

VI. LINEAR RESPONSE FUNCTIONS AND DYNAMIC SUSCEPTIBILITIES

Assume that a weak field F(t) is applied on a system and an observable variable B(t) changes its expectation value from the equilibrium value B⁰ to á B(t) ñ = B⁰ + dB(t), where dB(t) is a linear functional of F(t). Assume also that, as is usually the case, the interaction energy of F with the system may be written in the form (t) = - A(t) F(t), where A(t) is a dynamical variable of the system. The linear response function F_BA is implicitly defined by

Assuming, for simplicity, that the equilibrium values A⁰ and B⁰ are zero, the "Kubo formula" (see equation 4.2.20 of reference [7]), which relates the correlation function of A and B to F_BA, in the classical case, may be written as

where b = 1/k_BT and á ¼ ñ means average over an equilibrium distribution for a canonical ensemble of systems.

In the case of a magnetic field H_i(t), applied in direction i, the response function F_ji of the magnetization in direction j is

The numerical simulation starts from an equilibrium ensemble of N particles, like introduced in section III, follows the evolution of µ_j(t) as explained in section V, makes the product _j(t)µ_i(0) for each particle, sums over all particles and divides by N.

In what follows we apply the above "recipe" to calculate the magnetic response functions of the super-paramagnetic particles in ferrofluids, for several sets of material parameters. From the response functions, by Fourier transform, we obtain the corresponding dynamic susceptibilities,

The absorption power in magnetic resonance is proportional to the imaginary part of the susceptibility,

which makes the susceptibility a very important function.

We chose a set of values for the material parameters, in adimensional units, as our standard. Then, by varying one parameter in each calculation, we isolate its effect on the dynamics. For our standard we chose:

Figure 3 shows the results of the calculations for the response functions F_xx(t) and F_yx(t), for the case of the standard values of the parameters. We used an ensemble of 1000000 particles.

Figure 4 shows the imaginary part of the susceptibility, (w), for the standard case and several values of the anisotropy, keeping the standard values for the other parameters. One sees that in absence of anisotropy there is a very narrow, well pronounced, resonance line, while for increasing anisotropy the resonance becomes broader and less well pronounced. For very high K, which corresponds to "blocked" magnetic moment, the resonance disappears.

Figure 5 shows the effect of both viscosities on the susceptibility. The doted and dashed lines correspond to increase g_m and g_l, respectively, by a factor of 10. One sees that in both cases the line becomes broader, less intense and centered on a slightly higher resonance frequency.

Figure 6 shows the effect of changing temperature. At higher temperature the particle's random rotation becomes faster and, in consequence, the broadening effect of the anisotropy field decreases.

VII. CONCLUSIONS

In conclusion one can say that:

1) Simulation of the dynamics of super-paramagnetic moments in ferrofluids, when compared to resonance experiments, can give important information about their properties:

a) Bigger K: broader resonance line, same resonance frequency;

b) Very big K (blocked): no resonance;

c) Bigger g_m or g_l : broader resonance line and slightly higher resonance frequency;

d) Higher temperature: narrower line, same resonance frequency;

2) Second order Stratonovich is a very good algorithm to simulate realizations of the dynamics of super-paramagnetic moments in ferrofluids.

Received on 7 September, 2005

• [1] M.I. Shliomis and V.I. Stepanov, Adv. Chem. Phys. 87, 1 (1994).
• [2] C. Scherer and H. G. Matuttis, Phys. Rev. E 63, 0115504 (2001).
• [3] C. Scherer and T. F. Ricci Braz. J. Phys. 31, 380 (2001).
• [4] W.H. Press et. al., Numerical Recipies, Cambridge University Press, New York, 1992
• [5] C. Scherer, Métodos Computacionais da Física, ed. Livraria da Física, Universidade de Săo Paulo (2005);
• [6] C.W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and Natural Sciences, Springer-Verlag, Berlin, 1999;
• [7] R. Kubo, M. Toda and N. Hashitsume, Statistical Physics II: Nonequilibrium Statistical Mechanics, Springer-Verlag, Berlin, Heidelberg, 1985;

Publication Dates

• Publication in this collection
  16 Oct 2006
• Date of issue
  Sept 2006

History

• Received
  07 Sept 2005