Ferrofluids are colloidal suspensions of ferromagnetic single-domain particles in a liquid carrier. As the single particles contain only one magnetic domain, they can be seen as small permanent magnets. To prevent agglomeration of the particles, due to van-der-Waals or magnetic attraction, they are usually sterically or electrostatically stabilized (see figure 1). The former is achieved by adsorption of long chain molecules onto the particle surface, the latter by adsorption of charged coating particles. The size of the ferromagnetic particles are in the region of 10 nm. With the surfactant layer added they can reach a size of a few hundred nanometers. Have in mind that if we refer to the particle diameter $\sigma$ we mean the diameter of the magnetic core plus two times the thickness of the surfactant layer.

Some of the liquid properties, like the viscosity, the phase behavior or the optical birefringence can be altered via an external magnetic field or simply the fluid can be guided by such an field. Thus ferrofluids possess a wide range of biomedical applications like magnetic drug targeting or magnetic thermoablation and technical applications like fine positioning systems or adaptive bearings and dampers. In figure 2 the picture of a ferrofluid exposed to the magnetic field of a permanent magnet is shown. The famous energy minimizing thorn-like surface is clearly visible.

For simplicity in this tutorial we simulate spherical particles in a monodisperse ferrofluid system which means all particles have the same diameter $\sigma$ and dipole moment $\mu$. The point dipole moment is placed at the center of the particles and is constant both in magnitude and direction (in the coordinate system of the particle). This can be justified as the Néel relaxation times are usually negligible for the usual sizes of ferrofluid particles. Thus the magnetic interaction potential between two single particles is the dipole-dipole interaction potential which reads

\begin{equation*} u_{\text{DD}}(\vec{r}_{ij}, \vec{\mu}_i, \vec{\mu}_j) = \gamma \left(\frac{\vec{\mu}_i \cdot \vec{\mu}_j}{r_{ij}^3} - 3\frac{(\vec{\mu}_i \cdot \vec{r}_{ij}) \cdot (\vec{\mu}_j \cdot \vec{r}_{ij})}{r_{ij}^5}\right) \end{equation*}with $\gamma = \frac{\mu_0}{4 \pi}$ and $\mu_0$ the vacuum permeability.

For the steric interaction in this tutorial we use the purely repulsive Weeks-Chandler-Andersen (WCA) potential which is a Lennard-Jones potential with cut-off radius $r_{\text{cut}}$ at the minimum of the potential $r_{\text{cut}} = r_{\text{min}} = 2^{\frac{1}{6}}\cdot \sigma$ and shifted by $\varepsilon_{ij}$ such that the potential is continuous at the cut-off radius. Thus the potential has the shape

\begin{equation*} u_{\text{sr}}^{\text{WCA}}(r_{ij}) = \left\{ \begin{array}{ll} 4\varepsilon_{ij}\left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}} \right)^6 \right] + \varepsilon_{ij} & r_{ij} < r_{\text{cut}} \\ 0 & r_{ij} \geq r_{\text{cut}} \\ \end{array} \right. \end{equation*}where $r_{ij}$ are the distances between two particles. The purely repulsive character of the potential can be justified by the fact that the particles in real ferrofluids are sterically or electrostatically stabilized against agglomeration.

The whole interaction potential reads

\begin{equation*} u(\vec{r}_{ij}, \vec{\mu}_i, \vec{\mu}_j) = u_{\text{sr}}(\vec{r}_{ij}) + u_{\text{DD}}(\vec{r}_{ij}, \vec{\mu}_i, \vec{\mu}_j) \end{equation*}The liquid carrier of the system is simulated through a Langevin thermostat.

For ferrofluid systems there are three important parameters. The first is the volume fraction in three dimensions or the area fraction in two dimensions or quasi two dimensions. The second is the dipolar interaction parameter $\lambda$

\begin{equation} \lambda = \frac{\tilde{u}_{\text{DD}}}{u_T} = \gamma \frac{\mu^2}{k_{\text{B}}T\sigma^3} \end{equation}where $u_\mathrm{T} = k_{\text{B}}T$ is the thermal energy and $\tilde{u}_{DD}$ is the absolute value of the dipole-dipole interaction energy at close contact (cc) and head-to-tail configuration (htt) (see figure 3) per particle, i.e. in formulas it reads

\begin{equation} \tilde{u}_{\text{DD}} = \frac{ \left| u_{\text{DD}}^{\text{htt, cc}} \right| }{2} \end{equation}The third parameter takes a possible external magnetic field into account and is called Langevin parameter $\alpha$. It is the ratio between the energy of a dipole moment in the external magnetic field $B$ and the thermal energy

\begin{equation} \alpha = \frac{\mu_0 \mu}{k_{\text{B}} T}B \end{equation}The aim of this tutorial is to introduce the basic features of **ESPResSo** for ferrofluids or dipolar fluids in general. In **part I** and **part II** we will do this for a monolayer-ferrofluid, in **part III** for a three dimensional system. In **part I** we will examine the clusters which are present in all interesting ferrofluid systems. In **part II** we will examine the influence of the dipole-dipole-interaction on the magnetization curve of a ferrofluid. In **part III** we calculate estimators for the initial susceptibility using fluctuation formulas and sample the magnetization curve.

We assume the reader is familiar with the basic concepts of Python and MD simulations.

**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.

For this tutorial the following features of **ESPResSo** are needed

```
#define DIPOLES
#define DP3M
#define LENNARD_JONES
```

Please uncomment them in the `myconfig.hpp` and compile **ESPResSo** using this `myconfig.hpp`.

We start with checking for the presence of ESPResSo features and importing all necessary packages.

In [1]:

```
import espressomd
import espressomd.magnetostatics
import espressomd.cluster_analysis
import espressomd.pair_criteria
import espressomd.zn
espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])
import numpy as np
import tqdm
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
plt.rcParams.update({'font.size': 18})
```

Now we set up all simulation parameters.

In [2]:

```
# Lennard-Jones parameters
LJ_SIGMA = 1
LJ_EPSILON = 1
LJ_CUT = 2**(1. / 6.) * LJ_SIGMA
# Particles
N_PART = 1200
MASS=1.0
# Area fraction of the mono-layer
PHI = 0.1
# Dipolar interaction parameter lambda = mu_0 m^2 /(4 pi sigma^3 KT)
DIP_LAMBDA = 4
# Temperature
KT = 1.0
# Friction coefficient
GAMMA = 1.0
# Time step
TIME_STEP = 0.01
```

Note that we declared a `lj_cut`. This will be used as the cut-off radius of the Lennard-Jones potential to obtain a purely repulsive WCA potential.

Now we set up the system. The length of the simulation box is calculated using the desired area fraction and the area all particles occupy. Then we create the **ESPResSo** system and pass the simulation step. For the Verlet list skin parameter we use the built-in tuning algorithm of **ESPResSo**.

**Exercise:**

How large does `BOX_SIZE`

have to be for a system of `N_PART`

particles with a volume (area) fraction `PHI`

?
Define `BOX_SIZE`

.

$$
L_{\text{box}} = \sqrt{\frac{N A_{\text{sphere}}}{\varphi}}
$$

In [3]:

```
# SOLUTION CELL
BOX_SIZE = (N_PART * np.pi * (LJ_SIGMA / 2.)**2 / PHI)**0.5
```

In [4]:

```
# System setup
# BOX_SIZE = ...
print("Box size", BOX_SIZE)
# Note that the dipolar P3M and dipolar layer correction need a cubic
# simulation box for technical reasons.
system = espressomd.System(box_l=(BOX_SIZE, BOX_SIZE, BOX_SIZE))
system.time_step = TIME_STEP
```

Box size 97.08129562778495

In [5]:

```
# Lennard-Jones interaction
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=LJ_EPSILON, sigma=LJ_SIGMA, cutoff=LJ_CUT, shift="auto")
```

Now we generate random positions and orientations of the particles and their dipole moments.

Hint: It should be noted that we seed the random number generator of numpy. Thus the initial configuration of our system is the same every time this script will be executed. You can change it to another one to simulate with a different initial configuration.

**Exercise:**

How does one set up randomly oriented dipole moments?

Hint: Think of the way that different methods could introduce a bias in the distribution of the orientations.

Create a variable `dip`

as a `N_PART x 3`

numpy array, which contains the randomly distributed dipole moments.

In [6]:

```
# SOLUTION CELL
# Random dipole moments
np.random.seed(seed=1)
dip_phi = 2. * np.pi * np.random.random((N_PART, 1))
dip_cos_theta = 2 * np.random.random((N_PART, 1)) - 1
dip_sin_theta = np.sin(np.arccos(dip_cos_theta))
dip = np.hstack((
dip_sin_theta * np.sin(dip_phi),
dip_sin_theta * np.cos(dip_phi),
dip_cos_theta))
```

In [7]:

```
# Random dipole moments
# ...
# dip = ...
# Random positions in the monolayer
pos = BOX_SIZE * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))
```

`fix` argument.

In [8]:

```
# Add particles
system.part.add(pos=pos, rotation=N_PART * [(True, True, True)],
dip=dip, fix=N_PART * [(False, False, True)])
```

Out[8]:

<espressomd.particle_data.ParticleSlice at 0x7ae8e401a080>

Be aware that we do not set the magnitude of the magnetic dipole moments to the particles. As in our case all particles have the same dipole moment it is possible to rewrite the dipole-dipole interaction potential to

\begin{equation} u_{\text{DD}}(\vec{r}_{ij}, \vec{\mu}_i, \vec{\mu}_j) = \gamma \mu^2 \left(\frac{\vec{\hat{\mu}}_i \cdot \vec{\hat{\mu}}_j}{r_{ij}^3} - 3\frac{(\vec{\hat{\mu}}_i \cdot \vec{r}_{ij}) \cdot (\vec{\hat{\mu}}_j \cdot \vec{r}_{ij})}{r_{ij}^5}\right) \end{equation}where $\vec{\hat{\mu}}_i$ is the unit vector of the dipole moment $i$ and $\mu$ is the magnitude of the dipole moments. Thus we can only prescribe the initial orientation of the dipole moment to the particles and take the magnitude of the moments into account when calculating the dipole-dipole interaction with Dipolar P3M, by modifying the original Dipolar P3M prefactor $\gamma$ such that

\begin{equation} \tilde{\gamma} = \gamma \mu^2 = \frac{\mu_0}{4\pi}\mu^2 = \lambda \sigma^3 k_{\text{B}}T \end{equation}Of course it would also be possible to prescribe the whole dipole moment vectors to every particle and leave the prefactor of Dipolar P3M unchanged ($\gamma$). In fact we have to do this if we want to simulate polydisperse systems.

Now we choose the steepest descent integrator to remove possible overlaps of the particles.

In [9]:

```
# Choice of Fmax: Force acting on a particle moving 0.01*Sigma in a single time step
FMAX = 0.01 * LJ_SIGMA * MASS / system.time_step**2
system.integrator.set_steepest_descent(
f_max=FMAX,
gamma=0.1,
max_displacement=0.05)
system.integrator.run(5000)
assert np.all(np.abs(system.part.all().f) < FMAX), "Overlap removal did not converge!"
```

For the simulation of our system we choose the velocity Verlet integrator. After that we set up the thermostat which is, in our case, a Langevin thermostat to simulate in an NVT ensemble.

Hint: It should be noted that we seed the Langevin thermostat, thus the time evolution of the system is partly predefined. Partly because of the numeric accuracy and the automatic tuning algorithms of Dipolar P3M and DLC where the resulting parameters are slightly different every time. You can change the seed to get a guaranteed different time evolution.

In [10]:

```
# Switch to velocity Verlet integrator
system.integrator.set_vv()
system.thermostat.set_langevin(kT=KT, gamma=GAMMA, seed=1)
```

`DipolarP3M` to our system, a tuning function is started automatically
which tries to find the optimal parameters for dipolar P3M and prints them to the screen.
The last line of the output is the value of the tuned skin.

In [11]:

```
CI_DP3M_PARAMS = {'cao':3,'r_cut':8.34,'mesh':[8,8,8],'alpha':0.2115,'tune':False}
# Setup dipolar P3M and dipolar layer correction
dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT, **CI_DP3M_PARAMS)
mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=BOX_SIZE - LJ_SIGMA)
system.magnetostatics.solver = mdlc
# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)
# print skin value
print(f'tuned skin = {system.cell_system.skin:.2f}')
```

tuned skin = 0.55

Now we equilibrate the dipole-dipole interaction for some time

In [12]:

```
# equilibrate
EQUIL_ROUNDS = 10
EQUIL_STEPS = 100
for i in tqdm.trange(EQUIL_ROUNDS):
system.integrator.run(EQUIL_STEPS)
```

100%|██████████| 10/10 [00:02<00:00, 4.69it/s]

In [13]:

```
LOOPS = 100
```

**Exercise:**

Setup a cluster analysis object (`ClusterStructure`

class) and assign its instance to the variable `cluster_structure`

.
As criterion for the cluster analysis use a distance criterion where particles are assumed to be
part of a cluster if the nearest neighbors are closer than $1.3\sigma_{\text{LJ}}$.

In [14]:

```
# SOLUTION CELL
# Setup cluster analysis
cluster_structure = espressomd.cluster_analysis.ClusterStructure(
system=system,
pair_criterion=espressomd.pair_criteria.DistanceCriterion(cut_off=1.3 * LJ_SIGMA))
```

In [15]:

```
n_clusters = []
cluster_sizes = []
```

In [16]:

```
color = {0: "#7fc454"}
radii = {0: LJ_SIGMA/2}
vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)
```

**Exercise:**

Write an integration loop which runs a cluster analysis on the system, saving the number of clusters `n_clusters`

and the size distribution `cluster_sizes`

.
Take the following as a starting point:

```
for i in tqdm.trange(LOOPS):
# Run cluster analysis
cluster_structure.run_for_all_pairs()
# Gather statistics:
n_clusters.append(# < excercise >)
for c in cluster_structure.clusters:
cluster_sizes.append(# < excercise >)
system.integrator.run(100)
vis.update()
```

In [17]:

```
# SOLUTION CELL
for i in tqdm.trange(LOOPS):
# Run cluster analysis
cluster_structure.run_for_all_pairs()
# Gather statistics:
n_clusters.append(len(cluster_structure.clusters))
for c in cluster_structure.clusters:
cluster_sizes.append(c[1].size())
system.integrator.run(100)
vis.update()
```

100%|██████████| 100/100 [00:24<00:00, 4.16it/s]

`matplotlib`.

In [18]:

```
plt.figure(figsize=(10, 10))
plt.xlim(0, BOX_SIZE)
plt.ylim(0, BOX_SIZE)
plt.xlabel('x-position', fontsize=20)
plt.ylabel('y-position', fontsize=20)
plt.plot(system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1], 'o')
plt.show()
```

**Exercise:**

Use `numpy`

to calculate a histogram of the cluster sizes and assign it to the variable `size_dist`

.
Take only clusters up to a size of 19 particles into account.

Hint: In order not to count clusters with size 20 or more, one may include an additional bin containing these.
The reason for that is that `numpy`

defines the histogram bins as half-open intervals with the open border at the upper bin edge.
Consequently clusters of larger sizes are attributed to the last bin.
By not using the last bin in the plot below, these clusters can effectively be neglected.

In [19]:

```
# SOLUTION CELL
size_dist = np.histogram(cluster_sizes, range=(2, 21), bins=19)
```

In [20]:

```
import scipy.optimize
def kernel(x, a, b):
return a * np.exp(-np.abs(b) * x)
xdata = size_dist[1][:-2]
ydata = size_dist[0][:-1] / float(LOOPS)
popt, _ = scipy.optimize.curve_fit(kernel, xdata[2:], ydata[2:])
xdata_fit = np.linspace(np.min(xdata), np.max(xdata), 100)
ydata_fit = kernel(xdata_fit, *popt)
plt.figure(figsize=(10, 6))
plt.plot(xdata_fit, ydata_fit, label='Exponential fit')
plt.plot(xdata, ydata, 'o', label='Simulation results')
plt.xlabel('Cluster size', fontsize=20)
plt.ylabel('Probability density function', fontsize=20)
plt.legend()
plt.gca().get_xaxis().set_major_locator(ticker.MaxNLocator(integer=True))
plt.show()
```

[1] Juan J. Cerdà, V. Ballenegger, O. Lenz, and Ch. Holm. P3M algorithm for dipolar interactions. *Journal of Chemical Physics*, 129:234104, 2008.

[2] A. Bródka. Ewald summation method with electrostatic layer correction for
interactions of point dipoles in slab geometry. *Chemical Physics Letters* 400(1–3): 62–67, 2004. DOI:10.1016/j.cplett.2004.10.086

Image sources:

[3] Ayouril, Electro-Steric Stabilization, CC BY-SA 3.0

[4] Gregory F. Maxwell <gmaxwell@gmail.com> PGP:0xB0413BFA, Ferrofluid Magnet under glass edit, CC BY-SA 3.0