Ferrofluid - Part III

Table of Contents

  1. Susceptibility with fluctuation formulas
    1. Derivation of the fluctuation formula
    2. Simulation
  2. Magnetization curve of a 3D system

Remark: The equilibration and sampling times used in this tutorial would be not sufficient for scientific purposes, but they are long enough to get at least a qualitative insight of the behaviour of ferrofluids. They have been shortened so we achieve reasonable computation times for the purpose of a tutorial.

Susceptibility with fluctuation formulas

In this part we want to calculate estimators for the initial susceptibility, i.e. the susceptibility at zero external magnetic field. One could carry out several simulations with different external magnetic field strengths and get the initial susceptibility by fitting a line to the results. We want to go a more elegant way by using fluctuation formulas known from statistical mechanics. In three dimensions the initial susceptibility $\chi_{init}$ can be calculated with zero field simulations through

\begin{equation} \chi_{init} = \frac{V \cdot \mu_0}{3 \cdot k_B T} \left( \langle \boldsymbol{M}^2 \rangle - \langle \boldsymbol{M} \rangle^2 \right) = \frac{\mu_0}{3 \cdot k_B T \cdot V} \left( \langle \boldsymbol{\mu}^2 \rangle - \langle \boldsymbol{\mu} \rangle^2 \right) \end{equation}

where $\boldsymbol{M}$ is the magnetization vector and $\boldsymbol{\mu}$ is the total magnetic dipole moment of the system. In direction $i$ it reads

\begin{equation} M_i = \frac{1}{V} \Bigg\langle \sum_{j=1}^N \tilde{\mu}_j^i \Bigg\rangle \end{equation}

where $\tilde{\mu}_j^i$ is the $j$ th dipole moment in direction $i$.

Derivation of the fluctuation formula

We want to derive the fluctuation formula. We start with the definition of the magnetic susceptibility. In general this reads

\begin{equation} \chi \equiv \frac{\partial}{\partial H} \langle M_{\boldsymbol{H}} \rangle \end{equation}

with $\langle M_{\boldsymbol{H}} \rangle$ the ensemble averaged magnetization in direction of a homogeneous external magnetic field $\boldsymbol{H}$.

In thermal equilibrium the ensemble average of the magnetization reads

\begin{equation} \langle M_{\boldsymbol{H}} \rangle = \frac{1}{V Z_c} \left \lbrack \sum_{\alpha} \mu_{\boldsymbol{H},\alpha} e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }\right \rbrack \end{equation}

with $Z_c$ the canonical partition function, $E_{\alpha}(H=0)$ the energy without an external magnetic field $\boldsymbol{H}$, $\beta$ the inverse thermal energy $\frac{1}{k_BT}$, $\mu_{\boldsymbol{H},\alpha}$ the total magnetic dipole moment of the system in direction of the external magnetic field $\boldsymbol{H}$ in microstate $\alpha$ and $V$ the system volume.

Now we insert the magnetization $\langle M_{\boldsymbol{H}} \rangle$ in the definition of the magnetic susceptibility $\chi$ and let the derivative operate on the ensemble average. We get the fluctuation formula

\begin{equation} \chi = \frac{\beta\mu_0}{V} \left \lbrack \frac{1}{Z_c}\sum_{\alpha} \mu_{\alpha}^2~ e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H } - \frac{1}{Z_c}\sum_{\alpha} \mu_{\alpha}~ e^{ -\beta E_{\alpha}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }~~ \frac{1}{Z_c}\sum_{\alpha'}\mu_{\alpha'}~ e^{ -\beta E_{\alpha'}(H=0) + \beta\mu_0\mu_{\boldsymbol{H},\alpha}H }\right \rbrack = \frac{\beta\mu_0}{V} \left \lbrack \langle \mu_{\boldsymbol{H}}^2 \rangle - \langle \mu_{\boldsymbol{H}} \rangle^2 \right \rbrack = \frac{\beta\mu_0}{V} \left(\Delta \mu_{\boldsymbol{H}}\right)^2 \end{equation}

At zero external magnetic field ($H = 0$) there is no distinct direction for the system, so we can take the fluctuations $\Delta \mu$ in all directions and divide it by the dimension. Thus we can use more data points of our simulation for the average and get a more precise estimator for the susceptibility. Thus finally the fluctuation formula for the initial susceptibility in three dimensions reads

\begin{equation} \chi_{init} = \frac{\beta\mu_0}{3V} \left \lbrack \langle \boldsymbol{\mu}^2 \rangle - \langle \boldsymbol{\mu} \rangle^2 \right \rbrack = \frac{V\beta\mu_0}{3} \left \lbrack \langle \boldsymbol{M}^2 \rangle - \langle \boldsymbol{M} \rangle^2 \right \rbrack \end{equation}

where $\boldsymbol{\mu}$ and $\boldsymbol{M}$ are defined above.


In this part we want to consider a three dimensional ferrofluid system and compare our result for the initial susceptibility $\chi_{init}$ with them of Ref. [1].

First we import all necessary packages and check for the required ESPResSo features

Now we set up all necessary simulation parameters

and the system, where we, as we did in part I, only commit the orientation of the dipole moment to the particles and take the magnitude into account in the prefactor of Dipolar P3M (for more details see part I).

Hint: It should be noted that we seed both the Langevin thermostat and the random number generator of numpy. Latter means that the initial configuration of our system is the same every time this script is executed. As the time evolution of the system depends not solely on the Langevin thermostat but also on the numeric accuracy and DP3M (the tuned parameters are slightly different every time) it is only partly predefined. You can change the seeds to simulate with a different initial configuration and a guaranteed different time evolution.

Now we equilibrate for a while

As we need the magnetization of our system, we import MagneticDipoleMoment from observables which returns us the total dipole moment of the system which is the magnetization times the volume of the system.

Now we set the desired number of loops for the sampling

and sample the first and second moment of the magnetization or total dipole moment, by averaging over all total dipole moments occurring during the simulation

For the estimator of the initial susceptibility $\chi_{init}$ we need the magnitude of one single dipole moment

Now we can calculate $\chi_{init}$ from our simulation data

and print the result

Compared with the value $\chi = 0.822 \pm 0.017$ of Ref. [1] (see table 1) it should be very similar.

Now we want to compare the result with the theoretical expectations. At first with the simple Langevin susceptibility

and at second with the more advanced one (see Ref. [1] eq. (6)) which has a cubic accuracy in $\chi_L$ and reads

\begin{equation} \chi = \chi_L \left( 1 + \frac{\chi_L}{3} + \frac{\chi_L^2}{144} \right) \end{equation}

Both of them should be smaller than our result, but the second one should be closer to our one. The deviation of the theoretical results to our simulation result can be explained by the fact that in the Langevin model there are no interactions between the particles incorporated at all and the more advanced (mean-field-type) one of Ref. [1] do not take occurring cluster formations into account but assumes a homogeneous distribution of the particles. For higher values of the volume fraction $\phi$ and the dipolar interaction parameter $\lambda$ the deviations will increase as the cluster formation will become more pronounced.

Magnetization curve of a 3D system

At the end of this tutorial we now want to sample the magnetization curve of a three dimensional system and compare the results with analytical solutions. Again we will compare with the Langevin function but also with the approximation of Ref. [2] (see also Ref. [1] for the right coefficients) which takes the dipole-dipole interaction into account. For this approximation, which is a modified mean-field theory based on the pair correlation function, the Langevin parameter $\alpha$ is replaced by

\begin{equation} \alpha' = \alpha + \chi_L~L(\alpha) + \frac{\chi_L^{2}}{16} L(\alpha) \frac{d L(\alpha)}{d\alpha} \end{equation}

where $\chi_L$ is the Langevin susceptibility

\begin{equation} \chi_L = \frac{N}{V}\frac{\mu_0 \mu^2}{3k_BT} = 8 \cdot \lambda \cdot \phi \end{equation}

Analogous to part II we start at zero external magnetic field and increase the external field successively. At every value of the external field we sample the total dipole moment which is proportional to the magnetization as we have a fixed volume.

First we create a list of values of the Langevin parameter $\alpha$. As we already sampled the magnetization at zero external field in the last section, we take this value and continue with the sampling of an external field unequal zero

Now for each value in this list we sample the total dipole moment / magnetization of the system for a while. Have in mind that we have only committed the orientation of the dipole moments to the particles. Thus we have to commit $H\cdot \mu$ as the external magnetic field to ESPResSo, where $\mu$ is the magnitude of a single magnetic dipole moment.

As in part II we use the same system for every value of the Langevin parameter $\alpha$. Thus we use that the system is already pre-equilibrated from the previous run so we save some equilibration time. For scientific purposes one would use a new system for every value for the Langevin parameter to ensure that the systems are independent and no correlation effects are measured. Also one would perform more than just one simulation for each value of $\alpha$ to increase the precision of the results.

Now we define the Langevin function and the modified mean-field-approximation of the Langevin parameter of Ref. [2]

We also want to plot the linear approximation at $\alpha = 0$ to see for which values of $\alpha$ this approximation holds. We use the initial susceptibility calculated in the first chapter of this part as the gradient. As we want the gradient of $M^*$ with respect to $\alpha$ which fulfills the relation

\begin{equation} \frac{\partial M^*}{\partial \alpha} = \frac{1}{M_{sat}}\frac{\partial M}{\partial \left( \frac{\mu_0\mu}{k_BT} H\right)} = \frac{k_BT~V}{\mu_0\mu^2N}\frac{\partial M}{\partial H} = \frac{k_BT~V}{\mu_0\mu^2N}~\chi \end{equation}

we have to scale our calculated initial susceptibility $\chi_{init}$ by a factor to get it in our dimensionless units.

Now we plot the resulting curves together with our simulation results and the linear approximation

We can see that the magnetization curve where we used the Langevin parameter of the modified mean-field-theory is closer to our simulation results. Also we can clearly see that the linear approximation holds only for very small values of $\alpha$.

At this point is should be mentioned, that the modified mean-field-model assumes a spatially homogeneous system which is not the case at higher volume fractions $\phi$ and dipolar interaction parameters $\lambda$ as the particles form clusters. We can already see this with our simulation results as they visibly deviate from the modified mean-field-model.

At sufficiently high volume fractions $\phi$ and dipolar interaction parameters $\lambda$ these clusters can be thus rigid, that simulation with normal methods are impossible as the relaxation times exceeds normal simulation times by far resulting in strongly correlated configurations and thus measurements.

[1] Zuowei Wang, Christian Holm, and Hanns Walter Müller. “Molecular dynamics study on the equilibrium magnetization properties and structure of ferrofluids”. In: Phys. Rev. E 66 (2 Aug. 2002), p. 021405. doi: 10.1103/PhysRevE.66.021405. url: https://link.aps.org/doi/10.1103/PhysRevE.66.021405.

[2] Alexey O. Ivanov and Olga B. Kuznetsova. “Magnetic properties of dense ferrofluids: An influence of interparticle correlations”. In: Phys. Rev. E 64 (4 Sept. 2001), p. 041405. doi: 10.1103/PhysRevE.64.041405. url: https://link.aps.org/doi/10.1103/PhysRevE.64.041405.