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 4) 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`.

For interesting ferrofluid systems, where the fraction of ferromagnetic particles in the liquid carrier and their dipole moment are not vanishingly small, the ferromagnetic particles form clusters of different shapes and sizes. If the fraction and/or dipole moments are big enough the clusters can interconnect with each other and form a whole space occupying network. In this part we want to investigate the number of clusters as well as their shape and size in our simulated monolayer ferrofluid system. It should be noted that a monolayer is a quasi three dimensional system (q2D), i.e. two dimensional for the positions and three dimensional for the orientation of the dipole moments.

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
espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])
import numpy as np
import tqdm
%matplotlib inline
import matplotlib.pyplot as plt
```

In [2]:

```
import matplotlib.ticker as ticker
```

In [ ]:

```
plt.rcParams.update({'font.size': 18})
```

Now we set up all simulation parameters.

In [3]:

```
# Lennard-Jones parameters
LJ_SIGMA = 1
LJ_EPSILON = 1
LJ_CUT = 2**(1. / 6.) * LJ_SIGMA
# Particles
N_PART = 1200
# 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**.

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 [4]:

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

In [5]:

```
# 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
```

Now we set up the interaction between the particles as a non-bonded interaction and use the Lennard-Jones potential as the interaction potential. Here we use the above mentioned cut-off radius to get a purely repulsive interaction.

In [6]:

```
# 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.

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 [7]:

```
# 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 [8]:

```
# 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))))
```

Now we add the particles with their positions and orientations to our system. Thereby we activate all degrees of freedom for the orientation of the dipole moments. As we want a two dimensional system we only allow the particles to translate in $x$- and $y$-direction and not in $z$-direction by using the `fix` argument.

In [9]:

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

Out[9]:

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 [10]:

```
# Set integrator to steepest descent method
system.integrator.set_steepest_descent(
f_max=0, gamma=0.1, max_displacement=0.05)
```

Perform a steepest descent energy minimization. Track the relative energy change $E_{\text{rel}}$ per minimization loop (where the integrator is run for 10 steps) and terminate once $E_{\text{rel}} \le 0.05$, i.e. when there is less than a 5% difference in the relative energy change in between iterations.

In [11]:

```
import sys
energy = system.analysis.energy()['total']
relative_energy_change = 1.0
while relative_energy_change > 0.05:
system.integrator.run(10)
energy_new = system.analysis.energy()['total']
# Prevent division by zero errors:
if energy < sys.float_info.epsilon:
break
relative_energy_change = (energy - energy_new) / energy
print(f'Minimization, relative change in energy: {relative_energy_change:.4f}')
energy = energy_new
```

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 [12]:

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

To calculate the dipole-dipole interaction we use the Dipolar P3M method (see Ref. [1]) which is based on the Ewald summation. By default the boundary conditions of the system are set to conducting which means the dielectric constant is set to infinity for the surrounding medium. As we want to simulate a two dimensional system we additionally use the dipolar layer correction (DLC) (see Ref. [2]). As we add `DipolarP3M` to our system as an actor, 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 [13]:

```
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.actors.add(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}')
```

Now we equilibrate the dipole-dipole interaction for some time

In [14]:

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

The system will be sampled over 100 loops.

In [15]:

```
LOOPS = 100
```

As the system is two dimensional, we can simply do a scatter plot to get a visual representation of a system state. To get a better insight of how a ferrofluid system develops during time we will create a video of the development of our system during the sampling. If you only want to sample the system simply go to Sampling without animation

To get an animation of the system development we have to create a function which will save the video and embed it in an html string.

In [16]:

```
import matplotlib.animation as animation
```

In [17]:

```
import tempfile
import base64
VIDEO_TAG = """<video controls>
<source src="data:video/x-m4v;base64,{0}" type="video/mp4">
Your browser does not support the video tag.
</video>"""
def anim_to_html(anim):
if not hasattr(anim, '_encoded_video'):
with tempfile.NamedTemporaryFile(suffix='.mp4') as f:
anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264'])
with open(f.name, "rb") as g:
video = g.read()
anim._encoded_video = base64.b64encode(video).decode('ascii')
plt.close(anim._fig)
return VIDEO_TAG.format(anim._encoded_video)
animation.Animation._repr_html_ = anim_to_html
def init():
# Set x and y range
ax.set_ylim(0, BOX_SIZE)
ax.set_xlim(0, BOX_SIZE)
x_data, y_data = [], []
part.set_data(x_data, y_data)
return part,
```

In the following an animation loop is defined, however it is incomplete.
Extend the code such that in every loop the system is integrated for 100 steps.
Afterwards `x_data`

and `y_data`

have to be populated by the folded $x$- and $y$- positions of the particles.

(You may copy and paste the incomplete code template to the empty cell below.)

```
def run(i):
# < excercise >
# Save current system state as a plot
x_data, y_data = # < excercise >
ax.figure.canvas.draw()
part.set_data(x_data, y_data)
print(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', end='\r')
return part,
```

In [18]:

```
def run(i):
system.integrator.run(100)
# Save current system state as a plot
x_data, y_data = system.part.all().pos_folded[:, 0], system.part.all().pos_folded[:, 1]
ax.figure.canvas.draw()
part.set_data(x_data, y_data)
print(f'progress: {(i + 1) * 100. / LOOPS:3.0f}%', end='\r')
return part,
```

Now we use the `animation` class of `matplotlib` to save snapshots of the system as frames of a video which is then displayed after the sampling is finished. Between two frames are 100 integration steps.

In the video chain-like and ring-like clusters should be visible, as well as some isolated monomers.

In [19]:

```
fig, ax = plt.subplots(figsize=(10, 10))
part, = ax.plot([], [], 'o')
animation.FuncAnimation(fig, run, frames=LOOPS, blit=True, interval=0, repeat=False, init_func=init)
```

Out[19]: