Magnetization and susceptibility of ferrofluids
A second-order Taylor series expansion of the free energy functional provides analytical expressions for the magnetic field dependence of the free energy and of the magnetization of ferrofluids, here modelled by dipolar Yukawa interaction potentials. The corresponding hard core dipolar Yukawa reference fluid is studied within the framework of the mean spherical approximation. Our findings for the magnetic and phase equilibrium properties are in quantitative agreement with previously published and new Monte Carlo simulation data.
Ferrofluids are colloidal suspensions of ferromagnetic grains dispersed in a solvent. To keep such dispersions stable, one must avoid aggregation and counterbalance the attractive van der Waals interaction by repulsive forces, either electrostatic ones, if the particles are charged, or steric ones, if their surface is coated with polymers or surfactants. Each particle of a ferrofluid possesses a permanent magnetic dipole moment so that the ferrocolloids have much in common with dipolar liquids [1, 2]. The effective interaction of such magnetic particles can be modeled by the Stockmayer potential, which is the sum of an isotropic Lennard-Jones (LJ) and the anisotropic dipole-dipole interaction [3, 4]. Within this approach the LJ interaction parameters are considered to be effective ones, incorporating the effects of the solvent. (In a more refined approach one would model the suspension as a binary mixture of a nonpolar solvent and a polar solute component [5, 6].)
Based on previous experience [7, 8], here we substitute the Lennard-Jones potential by a hard core Yukawa potential. For a homogeneous isotropic dipolar Yukawa (DY) fluid this substitution provides the valuable benefit that one can solve analytically the corresponding Ornstein-Zernike equation within the mean spherical approximation (MSA). It turns out that the thermodynamic properties of a DY fluid can be calculated easily from the corresponding results for a hard core Yukawa and a dipolar hard sphere (DHS) fluid. The thermodynamic properties of a DY fluid differ from those of the DHS fluid but the magnetic susceptibility of the DY fluid is the same as that of the DHS fluid, which was predicted first by Wertheim . The disadvantage of the MSA is that it can predict the magnetization only for weak magnetic fields because it is a linear response theory. In order to eliminate this shortcoming, first Lebedev  suggested a semi-empirical method for using MSA at finite magnetic fields. Later Morozov and Lebedev  extended the applicability of DHS MSA to the case of arbitrary strengths of the magnetic fields on the basis of the Lovett-Mou-Buff equation , relating the single particle distribution function and the direct correlation function. Within the framework of density functional theory, here we also propose an equation for the magnetization of ferrofluids, which turns out to be simpler than the corresponding equation of Morozov and Lebedev . In order to facilitate the calculation of the phase diagrams, we present also the magnetic field dependence of the chemical potential and of the pressure.
2 The model
We study so-called hard core dipolar Yukawa fluids which consist of spherically symmetric particles interacting via a Yukawa (Y) potential characterized by parameters , , and :
In addition there is a dipolar interaction due to point dipoles embedded at the centres of the particles:
with the rotationally invariant function
where particle 1 (2) is located at () and carries a dipole moment of strength with an orientation given by the unit vector () with polar angles (); is the difference vector between the centres of particle 1 and 2, , and is a unit vector with orientation . The hard core dipolar Yukawa potential is defined by the aforementioned interaction potentials as
In the following we consider the DY fluid in a homogeneous external magnetic field (the direction of which is taken to coincide with the direction of the axis). The contribution of the latter to the interaction potential is
where the angle measures the orientation of the dipole of a single particle relative to the field direction.
3 Homogeneous magnetization in an external field
In general, an inhomogeneous and anisotropic dipolar fluid is described by the one particle distribution function . For bulk fluid phases the number density is spatially constant and is the orientational distribution function at the point . In a homogeneously magnetized bulk phase the orientational distribution function depends only on () and it is normalized to 1, i.e., . In these terms our subsequent bulk analysis is based on the minimization of the following grand canonical variational functional:
where is the Helmholtz free energy functional of the DY fluid, is the temperature and is the chemical potential. The Helmholtz free energy functional consists of an ideal gas and an excess contribution:
For the system under study the ideal contribution has the form
where is the de Broglie wavelength, is the inverse temperature and is the volume of the fluid. The central quantity in this theory is the excess free energy functional , which originates from interparticle interactions. The excess DY free energy functional for an anisotropic system is not known. However, it can be approximated by a functional Taylor series, expanded around a homogeneous isotropic reference system with bulk density . Neglecting all terms beyond second order one has
where is the difference between the anisotropic and the isotropic orientational distribution function. The first- () and second-order () direct correlation functions appearing in Eq. (9) are the first and second functional derivatives of () evaluated at the homogeneous and isotropic density . Since only the second- and higher-order direct correlation functions provide nonzero contributions to the functional . In the following we approximate the second-order direct correlation function of a DY fluid by the corresponding MSA solution due to Henderson et al :
(see Eq. (3)) with the rotationally invariant function
In Eq. (10), where is the hard core Yukawa direct correlation function , and are spherical harmonic components of the corresponding direct correlation function of a DHS fluid in MSA . The latter two functions can be expressed in terms of the Percus-Yevic hard sphere direct correlation function. If the volume of the fluid has the shape of a rotational ellipsoid (elongated around the magnetic field direction) and the magnetization is homogeneous, depends only on the angle relative to the long axis of the ellipsoid, and thus it can be expanded in terms of Legendre polynomials:
Due to one has . In order to avoid depolarization effects, we consider sample shapes of infinitely prolate ellipsoids, i.e., needle-shaped volumes . Since in an expansion in terms of the orthogonal functions and (see Eq. (10)) have contributions only of the order , is the only term in which contributes to Eq. (9). Thus after some elementary calculations one finds for the free energy functional of Eq. (9)
where is the Helmholtz free energy density of the isotropic DY fluid and
is the reduced inverse compressibility function of hard spheres within the Percus-Yevic approximation . The parameter stems from the DHS MSA and is given by the implicit equation
where is the Langevin susceptibility. Minimization of the grand canonical functional with respect to the orientational distribution function yields
with the normalization constant
The coefficient is determined by the implicit equation
where is the well known Langevin function. Further details of this minimization scheme can be found in Ref. . Minimization of with respect to yields
where is the pressure of the DY fluid. Each particle of the magnetic fluid carries a dipole moment which may be preferentially aligned in a particular direction. This gives rise to a magnetization (the direction of which coincides with the direction of the external magnetic field):
Accordingly, the zero field magnetic susceptibility is
4 Monte Carlo simulations
We carried out Monte Carlo simulations for DY fluids using NVT and NpT ensembles in order to verify the predictions of the present DFT. The initial configurations in a cubic box were created by random insertions of the particles, excluding particle overlaps. Boltzmann sampling and periodic boundary conditions with the minimum-image convention were applied. A spherical cutoff of the pair interaction potentials at half of the cell length was applied, and long-ranged corrections (LRC) were taken into account. For the dispersion part of the potential, the usual Yukawa-tail LRC was used while for the dipole-dipole interaction a reaction field LRC and a conducting boundary condition were applied. For the simulation of the magnetization after 40000 equilibration cycles 0.6-0.8 million production cycles were used involving 512 particles. In the simulations the equilibrium magnetization can be obtained from the expression
where the brackets denotes the ensemble average.
We obtained the external field dependence of the numerical data for liquid-vapor coexistence of a DY fluid by using the extended NpT plus test particle (NpT+TP) method, which is described in detail in Ref.  so that here we present only an outline. Starting from a point in the one-phase region of the parameter plane, the chemical potential can be expanded into a Taylor series about the point up to third order:
On the basis of basic thermodynamic relations the coefficients of this series can be expressed in terms of the derivatives of the enthalpy and the volume of the system with respect to and and they can be calculated from fluctuation formulas by performing an NpT+TP Monte Carlo simulation at the point . These derivatives and fluctuation formulas are given in Ref. . Carrying out this scheme for a point in the liquid and in the vapor phase, respectively, and rewriting the third-order Taylor series of for these points accordingly, the vapor pressure curve as well as other equilibrium data can be obtained from the intersection of these curves in the appropriate temperature range within an accuracy depending on the number of terms taken into account in the Taylor expansion. The NpT ensemble MC simulations involved 512 particles and about 0.8 million cycles were performed. The chemical potential was calculated by using Widom’s test particle method.
5 Results and discussion
In the following we shall use reduced quantities: as reduced temperature, as reduced density, as reduced dipole moment, as reduced magnetic field strength, and as the reduced magnetization. For the reference system, all results are obtained for (see Eq. (1)), and different from Ref.  the so-called first-order MSA  is used to describe the Yukawa fluid. The magnetization curves obtained from the numerical solution of Eqs. (23) and (15)
for a needle-shaped sample are displayed in Fig. 1(a) together with the result of the Langevin theory for . Figure 1(a) shows that the interparticle interaction enhances the magnetization relative to the corresponding values of the Langevin theory. The theoretical data of the present theory are consistent with our new simulation data. For the reduced density we find both DFT and MC results to be in excellent, quantitative agreement for all field strengths. For the level of quantitative agreement reduces for larger field strengths.
The density dependence of the magnetic susceptibility obtained from the numerical solution of Eqs. (23) and (15) in comparison with earlier simulation data of Szalai et al  is displayed in Fig. 1(b). For the reduced dipole moment the agreement between the theoretical and simulation results is quite good, especially at lower densities. For the level of quantitative agreement reduces for larger reduced densities.
The liquid-vapor phase diagrams at a given temperature follow from requiring the equality of the chemical potential and of the pressure for the densities of coexisting phases. The phase diagrams in the density-temperature plane are presented in Fig. 2(a) for various fixed values of the magnetic field and for the reduced dipole moment . Figure 2(a) shows that for and there is no spontaneous magnetization within the framework of this theory, i.e., there is no ferromagnetic liquid phase. For there is a first-order phase transition between a low density gas and a high density liquid phase and both phases exhibit a nonzero magnetization. The agreement between the NpT+TP simulation and the DFT data is reasonable. For Fig. 2(b) represents the magnetic field dependence of the critical temperature corresponding to the liquid-gas phase transitions, obtained from the present DFT and liquid-gas equilibrium NpT+TP simulation data. For small both the theoretical and the simulation shifts increase quadratically as a function of the field strength. Similar results have been found for Stockmayer fluids in Refs. [17, 18].
The following main results have been obtained.
(1) Based on a second-order Taylor series expansion of the free energy functional of the dipolar Yukawa fluids, within the mean spherical approximation we have derived an analytical expression for the magnetization of ferrofluids. There is good agreement between the results of the density functional theory and the Monte Carlo simulation data for reduced dipole moments . The proposed implicit equation for the magnetization (Eq. (23)) extends the applicability of MSA for arbitrary strengths of the magnetic fields.
-  Buyevich Y A and Ivanov A O 1992 Physica A 190 276
-  Huke B and Lücke 2004 Rep. Prog. Phys. 67 1731
-  Russier V and Douzi M 1994 J. Colloid Interface Sci. 162 356
-  Groh B and Dietrich S 1994 Phys. Rev. E 50 3814
-  Szalai I and Dietrich S 2005 Mol. Phys. 103 2873
-  Range G M and Klapp S H L 2004 Phys. Rev. E 70 031201
-  Henderson D, Boda D, Chan K-Y and Szalai I 1999 J. Chem. Phys. 110 7348
-  Szalai I, Henderson D, Boda D and Chan K-Y 1999 J. Chem. Phys. 111 337
-  Wertheim M S 1971 J. Chem. Phys. 55 4291
-  Lebedev A V 1989 Magnetohydrodynamics 25 520
-  Morozov K I and Lebedev A V 1990 J. Magn. Magn. Mater. 85 51
-  Lovett R, Mou C Y and Buff F P 1976 J. Chem. Phys. 65 570
-  Hausleitner C Hafner J 1988 J. Phys. F: Met. Phys. 18 1013
-  Percus J K and Yevick G J 1956 Phys. Rev. 110 1
-  Tang Y 2003 J. Chem. Phys. 118 4141
-  Boda D, Liszi J and Szalai I 1995 Chem. Phys. Lett. 235 140
-  Groh B and Dietrich S 1996 Phys. Rev. E 53 2509
-  Boda D, Winkelmann J, Liszi J and Szalai I 1996 Mol. Phys. 87 601