Ferrofluid - Part 2

Table of Contents

1. Applying an external magnetic field
2. Magnetization curve

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.

Applying an external magnetic field

In this part we want to investigate the influence of a homogeneous external magnetic field exposed to a ferrofluid system.

We import all necessary packages and check for the required ESPResSo features

import espressomd
import espressomd.magnetostatics
import espressomd.zn

espressomd.assert_features(['DIPOLES', 'DP3M', 'LENNARD_JONES'])

%matplotlib inline
import matplotlib.pyplot as plt
plt.rcParams.update({'font.size': 18})
import numpy as np
import tqdm

and set up the simulation parameters where we introduce a new dimensionless parameter

\begin{equation} \alpha = \frac{\mu B}{k_{\text{B}}T} = \frac{\mu \mu_0 H}{k_{\text{B}}T} \end{equation}

which is called Langevin parameter. We intentionally choose a relatively high volume fraction $\phi$ and dipolar interaction parameter $\lambda$ to clearly see the influence of the dipole-dipole interaction

# Lennard-Jones parameters
LJ_SIGMA = 1.
LJ_EPSILON = 1.
LJ_CUT = 2**(1. / 6.) * LJ_SIGMA

# Particles
N_PART = 700

# Area fraction of the mono-layer
PHI = 0.06

# 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

# Langevin parameter ALPHA = MU_0 m H / kT
ALPHA = 10.

# vacuum permeability
MU_0 = 1.

Now we set up the system. As in part I, the orientation of the dipole moments is set directly on the particles, whereas the magnitude of the moments is taken into account when determining the prefactor of the 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. The latter means that the initial configuration of our system is the same every time this script will be executed. As the time evolution of the system depends not solely on the Langevin thermostat but also on the numeric accuracy and DP3M as well as DLC (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.

# System setup
box_size = (N_PART * np.pi * (LJ_SIGMA / 2.)**2. / PHI)**0.5

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

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

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

# Random positions in the monolayer
pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))

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

# Remove overlap between particles by means of the steepest descent method
MASS = 1.0
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!"


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

# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)

# Setup dipolar P3M and dipolar layer correction (DLC)
dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT)
mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)
system.magnetostatics.solver = mdlc

# tune verlet list skin again
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}')

Box size 95.72344839677595
tuned skin = 0.90

We now apply the external magnetic field which is

\begin{equation} B = \mu_0 H = \frac{\alpha~k_{\text{B}}T}{\mu} \end{equation}

As only the current orientation of the dipole moments, i.e. the unit vector of the dipole moments, is saved in the particle list but not their magnitude, we have to use $B\cdot \mu$ as the strength of the external magnetic field. We will apply the field in x-direction using the class constraints of ESPResSo

# magnetic field times dipole moment
H_dipm = ALPHA * KT / MU_0
H_field = [H_dipm, 0, 0]

Exercise:

Define a homogenous magnetic field constraint using H_field and add it to system's contraints.

# SOLUTION CELL
H_constraint = espressomd.constraints.HomogeneousMagneticField(H=H_field)
system.constraints.add(H_constraint)

<espressomd.constraints.HomogeneousMagneticField at 0x770c9bf42440>

Equilibrate the system.

# Equilibrate
print("Equilibration...")
equil_rounds = 10
equil_steps = 100
for i in tqdm.trange(equil_rounds):
    system.integrator.run(equil_steps)

Equilibration...

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

Now we can visualize the current state and see that the particles mostly create chains oriented in the direction of the external magnetic field. Also some monomers should be present.

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(particles.pos_folded[:, 0], particles.pos_folded[:, 1], 'o')
plt.show()

You may want to get an insight of how the system develops in time. We use the ZnDraw visualiser:

color = {0: "#7fc454"}   
radii = {0: LJ_SIGMA/2}
vis = espressomd.zn.Visualizer(system, colors=color, radii=radii)
for _ in range(100):
    system.integrator.run(50)
    vis.update()

In the visualization video we can see that the single chains break and connect to each other during time. Also some monomers are present which break from and connect to chains.

Magnetization curve

An important observable of a ferrofluid system is the magnetization $M$ of the system in direction of an external magnetic field $H$

\begin{equation} M = \frac{\sum_i \mu_i^H}{V} \end{equation}

where the index $H$ means the component of $\mu_i$ in direction of the external magnetic field $H$ and the sum runs over all particles. $V$ is the system's volume.

The magnetization plotted over the external field $H$ is called magnetization curve. For particles with non-interacting dipole moments there is an analytical solution

\begin{equation} M = M_{\text{sat}}\cdot L(\alpha) \end{equation}

with $L(\alpha)$ the Langevin function

\begin{equation} L(\alpha) = \coth(\alpha)-\frac{1}{\alpha} \end{equation}

and $\alpha$ the Langevin parameter

\begin{equation} \alpha=\frac{\mu_0\mu}{k_{\text{B}}T}H = \frac{\mu}{k_{\text{B}}T}B \end{equation}

$M_{sat}$ is the so called saturation magnetization which is the magnetization of a system where all dipole moments are aligned to each other. Thus it is the maximum of the magnetization. In our case all dipole moments are equal, thus

\begin{equation} M_{\text{sat}} = \frac{N_{\text{part}}\cdot\mu}{V} \end{equation}

For better comparability we now introduce a dimensionless magnetization

\begin{equation} M^* = \frac{M}{M_{sat}} = \frac{\sum_i \mu_i^H}{N_{\text{part}}\cdot \mu} \end{equation}

Thus the analytical solution for non-interacting dipole moments $M^*$ is simply the Langevin function.

For interacting dipole moments there are only approximations for the magnetization curve available.

Here we want to use the approximation of Ref. [1] for a quasi two dimensional system, which reads with adjusted coefficients (Ref. [1] used a different dipole-dipole interaction prefactor $\gamma = 1$)

\begin{equation} M_{\parallel}^{\text{q2D}} = M_{\text{sat}} L(\alpha) \left( 1 + \mu_0\frac{1}{8} M_{\text{sat}} \frac{\mathrm{d} L(\alpha)}{\mathrm{d}B} \right) \end{equation}

and

\begin{equation} M_{\perp}^{\text{q2D}} = M_{\text{sat}} L(\alpha) \left( 1 - \mu_0\frac{1}{4} M_{\text{sat}} \frac{\mathrm{d} L(\alpha)}{\mathrm{d}B} \right) \end{equation}

for the magnetization with an external magnetic field parallel and perpendicular to the monolayer plane, respectively. Here the dipole-dipole interaction is approximated as a small perturbation and

\begin{equation} \frac{\mathrm{d} L(\alpha)}{\mathrm{d}B} = \left( \frac{1}{\alpha^2} - \frac{1}{\sinh^2(\alpha)} \right) \cdot \frac{\mu}{k_{\text{B}}T} \end{equation}

By comparing the magnetization curve parallel $M_{\parallel}^{\text{q2D}}$ and perpendicular $M_{\perp}^{\text{q2D}}$ to the monolayer plane we can see that the magnetization is increased in the case of an external field parallel to the monolayer plane and decreased in the case of an external field perpendicular to the monolayer plane. The latter can be explained by the fact that an orientation of all single dipole moments perpendicular to the monolayer plane results in a configuration with a repulsive dipole-dipole interaction as the particles have no freedom of movement in the direction perpendicular to the monolayer plane. This counteracts the magnetization perpendicular to the monolayer plane.

We now want to use ESPResSo to get an estimation of how the magnetization curve is affected by the dipole-dipole interaction parallel and perpendicular to the monolayer plane and compare the results with the Langevin curve and the magnetization curves of Ref. [1].

For the sampling of the magnetization curve we set up a new system, where we decrease the dipolar interaction parameter $\lambda$ drastically. We do this as we want to compare our results with the approximation of Ref. [1] which is only valid for small dipole-dipole interaction between the particles (decreasing the volume fraction $\phi$ would also be an appropriate choice). For smaller dipolar interaction parameters it is possible to increase the time step. We do this to get more uncorrelated measurements.

# Dipolar interaction parameter lambda = MU_0 m^2 /(4 pi sigma^3 kT)
DIP_LAMBDA = 1.

# increase time step
TIME_STEP = 0.02

# dipole moment
dipm = np.sqrt(4 * np.pi * DIP_LAMBDA * LJ_SIGMA**3. * KT / MU_0)

# remove all particles
system.part.clear()
system.magnetostatics.clear()
system.thermostat.turn_off()

# Random dipole moments
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))

# Random positions in the monolayer
pos = box_size * np.hstack((np.random.random((N_PART, 2)), np.zeros((N_PART, 1))))

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

# Remove overlap between particles by means of the steepest descent method
system.integrator.set_steepest_descent(f_max=0, gamma=0.1, max_displacement=0.05)

while system.analysis.energy()["total"] > 5 * KT * N_PART:
    system.integrator.run(20)

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

# tune verlet list skin
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)

# Setup dipolar P3M and dipolar layer correction
dp3m = espressomd.magnetostatics.DipolarP3M(accuracy=5E-4, prefactor=DIP_LAMBDA * LJ_SIGMA**3 * KT)
mdlc = espressomd.magnetostatics.DLC(actor=dp3m, maxPWerror=1E-4, gap_size=box_size - LJ_SIGMA)
system.magnetostatics.solver = mdlc

# tune verlet list skin again
system.cell_system.tune_skin(min_skin=0.4, max_skin=2., tol=0.2, int_steps=100)

DipolarP3M tune parameters: Accuracy goal = 5.00000e-04 prefactor = 4.00000e+00
System: box_l = 9.57234e+01 # magnetic part = 700 Sum[mu_i^2] = 7.00000e+02
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
8    3   7.31480e-02 2.54050e+01 4.966e-04 3.536e-04 3.488e-04 0.93    
8    2   7.97126e-02 2.19599e+01 4.959e-04 3.536e-04 3.478e-04 0.92    
8    4   6.93968e-02 2.76533e+01 4.991e-04 3.536e-04 3.522e-04 1.05    
8    1   9.19039e-02 1.67627e+01 4.947e-04 3.536e-04 3.460e-04 0.99    
10   3   7.97126e-02 2.19599e+01 5.311e-04 3.536e-04 3.963e-04 accuracy not achieved
10   4   7.97126e-02 2.19599e+01 5.057e-04 3.536e-04 3.616e-04 accuracy not achieved
10   5   7.90898e-02 2.22637e+01 4.987e-04 3.536e-04 3.518e-04 1.26    
10   6   7.84671e-02 2.25719e+01 4.981e-04 3.536e-04 3.508e-04 1.52    
10   7   7.84671e-02 2.25719e+01 4.954e-04 3.536e-04 3.470e-04 1.93    
12   3   7.90898e-02 2.22637e+01 5.839e-04 3.536e-04 4.647e-04 accuracy not achieved
12   4   7.90898e-02 2.22637e+01 5.431e-04 3.536e-04 4.123e-04 accuracy not achieved
12   5   7.90898e-02 2.22637e+01 5.210e-04 3.536e-04 3.826e-04 accuracy not achieved
12   6   7.90898e-02 2.22637e+01 5.090e-04 3.536e-04 3.662e-04 accuracy not achieved
12   7   7.90898e-02 2.22637e+01 5.027e-04 3.536e-04 3.573e-04 accuracy not achieved
14   3   7.90898e-02 2.22637e+01 6.003e-04 3.536e-04 4.851e-04 accuracy not achieved
14   4   7.90898e-02 2.22637e+01 5.493e-04 3.536e-04 4.203e-04 accuracy not achieved
14   5   7.90898e-02 2.22637e+01 5.199e-04 3.536e-04 3.812e-04 accuracy not achieved
14   6   7.90898e-02 2.22637e+01 5.027e-04 3.536e-04 3.574e-04 accuracy not achieved
14   7   7.90898e-02 2.22637e+01 4.927e-04 3.536e-04 3.431e-04 2.38    
14   6   7.90898e-02 2.22637e+01 5.027e-04 3.536e-04 3.574e-04 accuracy not achieved
16   3   7.90898e-02 2.22637e+01 5.878e-04 3.536e-04 4.696e-04 accuracy not achieved
16   4   7.90898e-02 2.22637e+01 5.324e-04 3.536e-04 3.980e-04 accuracy not achieved
16   5   7.90898e-02 2.22637e+01 5.003e-04 3.536e-04 3.539e-04 accuracy not achieved
16   6   7.84719e-02 2.25695e+01 4.910e-04 3.536e-04 3.407e-04 2.09    
16   7   7.72362e-02 2.31947e+01 4.990e-04 3.536e-04 3.521e-04 2.55    
18   3   7.84719e-02 2.25695e+01 5.720e-04 3.536e-04 4.496e-04 accuracy not achieved
18   4   7.84719e-02 2.25695e+01 5.142e-04 3.536e-04 3.733e-04 accuracy not achieved
18   5   7.78589e-02 2.28774e+01 4.932e-04 3.536e-04 3.439e-04 2.31    
18   6   7.66328e-02 2.35069e+01 4.944e-04 3.536e-04 3.456e-04 2.55    
18   7   7.60197e-02 2.38289e+01 4.899e-04 3.536e-04 3.391e-04 3.09    
20   3   7.78589e-02 2.28774e+01 5.461e-04 3.536e-04 4.162e-04 accuracy not achieved
20   4   7.78589e-02 2.28774e+01 4.892e-04 3.536e-04 3.382e-04 2.24    
20   5   7.60341e-02 2.38213e+01 4.946e-04 3.536e-04 3.459e-04 2.47    
20   6   7.48175e-02 2.44746e+01 4.957e-04 3.536e-04 3.475e-04 2.70    
20   7   7.42093e-02 2.48090e+01 4.900e-04 3.536e-04 3.392e-04 3.20    
22   3   7.78589e-02 2.28774e+01 5.005e-04 3.536e-04 3.543e-04 accuracy not achieved
22   4   7.60341e-02 2.38213e+01 4.890e-04 3.536e-04 3.378e-04 2.98    
22   5   7.42093e-02 2.48090e+01 4.953e-04 3.536e-04 3.469e-04 3.11    
22   6   7.29927e-02 2.54937e+01 4.968e-04 3.536e-04 3.491e-04 3.40    
22   7   7.23844e-02 2.58444e+01 4.904e-04 3.536e-04 3.398e-04 3.85    

resulting parameters: mesh: (8, 8, 8), cao: 2, r_cut_iL: 7.9713e-02,
                      alpha_L: 2.1960e+01, accuracy: 4.9593e-04, time: 0.92
DipolarP3M tune parameters: Accuracy goal = 5.00000e-04 prefactor = 1.00000e+00
System: box_l = 9.57234e+01 # magnetic part = 700 Sum[mu_i^2] = 7.00000e+02
mesh cao r_cut_iL    alpha_L     err       rs_err    ks_err    time [ms]
8    3   7.31480e-02 2.54050e+01 4.966e-04 3.536e-04 3.488e-04 0.87    
8    2   7.97126e-02 2.19599e+01 4.959e-04 3.536e-04 3.478e-04 0.87    
8    4   6.93968e-02 2.76533e+01 4.991e-04 3.536e-04 3.522e-04 1.05    
8    1   9.19039e-02 1.67627e+01 4.947e-04 3.536e-04 3.460e-04 0.97    
10   3   7.97126e-02 2.19599e+01 5.311e-04 3.536e-04 3.963e-04 accuracy not a

0.7000000000000001

To increase the performance we use the built-in function MagneticDipoleMoment to calculate the dipole moment of the whole system. In our case this is only the orientation as we never set the strength of the dipole moments on our particles.

Exercise:

Import the magnetic dipole moment observable and define an observable object dipm_tot. Use particle slicing to pass all particle ids.

# SOLUTION CELL
import espressomd.observables
dipm_tot = espressomd.observables.MagneticDipoleMoment(ids=particles.id)

We use the dimensionless Langevin parameter $\alpha$ as the parameter for the external magnetic field. As the interesting part of the magnetization curve is the one for small external magnetic field strengths—for large external magnetic fields the magnetization goes into saturation in all cases—we increase the spacing between the Langevin parameters $\alpha$ up to higher values and write them into a list

alphas = [0,1,2,3,4.5,8]

For both the magnetization perpendicular and parallel to the monolayer plane we use the same system for every value of the Langevin parameter $\alpha$. Thus we use that the system is already more or less 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 sample the magnetization for increasing $\alpha$ (increasing magnetic field strength) in direction perpendicular to the monolayer plane.

Exercise:

Complete the loop such that for every alpha a magnetic field of strength of the respective H_dipm is applied:

# sampling with magnetic field perpendicular to monolayer plane (in z-direction)

# remove all constraints
system.constraints.clear()

# list of magnetization in field direction
magnetization_perp = np.full_like(alphas, np.nan)

# number of loops for sampling
loops = 500

for ndx, alpha in enumerate(pbar := tqdm.tqdm(alphas)):
    pbar.set_description(f"Sampling for α={alpha:.2f}")
    # set magnetic field constraint
    H_dipm = alpha * KT
    # < exercise >

    # equilibration
    for i in range(equil_rounds):
        system.integrator.run(equil_steps)

    # sampling
    magn_temp = 0.
    for i in range(loops):
        system.integrator.run(20)
        magn_temp += dipm_tot.calculate()[2]

    # save average magnetization
    magnetization_perp[ndx] = magn_temp / loops

    # remove constraint
    system.constraints.clear()

# SOLUTION CELL
# sampling with magnetic field perpendicular to monolayer plane (in z-direction)

# remove all constraints
system.constraints.clear()

# list of magnetization in field direction
magnetization_perp = np.full_like(alphas, np.nan)

# number of loops for sampling
loops = 100

for ndx, alpha in enumerate(pbar := tqdm.tqdm(alphas)):
    pbar.set_description(f"Sampling for α={alpha:.2f}")
    # set magnetic field constraint
    H_dipm = alpha * KT
    H_field = [0, 0, H_dipm]
    H_constraint = espressomd.constraints.HomogeneousMagneticField(H=H_field)
    system.constraints.add(H_constraint)

    # equilibration
    for i in range(equil_rounds):
        system.integrator.run(equil_steps)

    # sampling
    magn_temp = 0.
    for i in range(loops):
        system.integrator.run(20)
        magn_temp += dipm_tot.calculate()[2]

    # save average magnetization
    magnetization_perp[ndx] = magn_temp / loops

    # remove constraint
    system.constraints.clear()

Sampling for α=8.00: 100%|██████████| 6/6 [00:15<00:00,  2.60s/it]

and now we sample the magnetization for increasing $\alpha$ or increasing magnetic field in direction parallel to the monolayer plane.

Exercise:

Use the code from the previous exercise as a template. Now sample the magnetization curve for a magnetic field parallel to the quasi-2D layer and store the calculated magnetizations in a list named magnetization_para (analogous to magnetization_perp).

Hint: Set up the field in $x$- or $y$-direction and sample the magnetization along the same axis.

# SOLUTION CELL
# sampling with magnetic field parallel to monolayer plane (in x-direction)

# remove all constraints
system.constraints.clear()

# list of magnetization in field direction
magnetization_para = np.full_like(alphas, np.nan)

# number of loops for sampling
loops = 100

for ndx, alpha in enumerate(pbar := tqdm.tqdm(alphas)):
    pbar.set_description(f"Sampling for α={alpha:.2f}")
    # set magnetic field constraint
    H_dipm = alpha * KT
    H_field = [H_dipm, 0, 0]
    H_constraint = espressomd.constraints.HomogeneousMagneticField(H=H_field)
    system.constraints.add(H_constraint)

    # equilibration
    for i in range(equil_rounds):
        system.integrator.run(equil_steps)

    # sampling
    magn_temp = 0
    for i in range(loops):
        system.integrator.run(20)
        magn_temp += dipm_tot.calculate()[0]

    # save average magnetization
    magnetization_para[ndx] = magn_temp / loops

    # remove constraint
    system.constraints.clear()

Sampling for α=8.00: 100%|██████████| 6/6 [00:15<00:00,  2.62s/it]

Now we can compare the resulting magnetization curves with the Langevin curve and the more advanced ones of Ref. [1] by plotting all of them in one figure.

For the approximations of $M_{\parallel}^{\text{q2D}}$ and $M_{\perp}^{\text{q2D}}$ of Ref. [1] we need the dipole moment of a single particle. Thus we calculate it from our dipolar interaction parameter $\lambda$

# dipole moment
dipm = np.sqrt(DIP_LAMBDA * 4 * np.pi * LJ_SIGMA**3. * KT / MU_0)
print(f'dipole moment = {dipm:.4f}')

dipole moment = 3.5449

and the saturation magnetization by using

\begin{equation} M_{\text{sat}} = \rho \mu = \phi \frac{4}{\pi \sigma^2} \mu \end{equation}

thus

M_sat = PHI * 4. / np.pi * 1. / (LJ_SIGMA**2.) * dipm

Further we need the derivation of the Langevin function after the external field $B$ thus we define the function

def dL_dB(alpha):
    return (1. / (alpha**2) - 1. / np.sinh(alpha)**2) * dipm / (KT)

Now we define the approximated magnetization curves parallel and perpendicular to the monolayer plane

# approximated magnetization curve for a field parallel to the monolayer plane
def magnetization_approx_para(alpha):
    return L(alpha) * (1. + MU_0 / 8. * M_sat * dL_dB(alpha))

# approximated magnetization curve for a field perpendicular to the monolayer plane
def magnetization_approx_perp(alpha):
    return L(alpha) * (1. - MU_0 / 4. * M_sat * dL_dB(alpha))

Now we define the Langevin function

# Langevin function
def L(x):
    return (np.cosh(x) / np.sinh(x)) - 1. / x

and plot the three theoretical curves together with our simulation results

Exercise:

The following listing plots the analytical models for the magnetization. Add your sampled magnetization curves (for parallel and perpendicular fields) each normalized by the number of particles $N_{\text{part}}$.

# list of the values for alpha (x-axis)
x = np.arange(0.01,9, 0.1, dtype=float)

plt.figure(figsize=(10,10))
plt.xlabel(r'$\alpha$', fontsize=20)
plt.ylabel(r'$M^*$', fontsize=20)
plt.plot(x, L(x), label='Langevin function')
plt.plot(x, magnetization_approx_perp(x), label=r'q2D approximation $\perp$')
plt.plot(x, magnetization_approx_para(x), label=r'q2D approximation $\parallel$')
# < exercise >
plt.legend(fontsize=20)
plt.show()

# SOLUTION CELL
# list of the values for alpha (x-axis)
x = np.arange(0.01, 9, 0.1, dtype=float)

plt.figure(figsize=(10, 10))
plt.xlabel(r'$\alpha$', fontsize=20)
plt.ylabel(r'$M^*$', fontsize=20)
plt.plot(x, L(x), label='Langevin function')
plt.plot(x, magnetization_approx_perp(x), label=r'q2D approximation $\perp$')
plt.plot(x, magnetization_approx_para(x), label=r'q2D approximation $\parallel$')
plt.plot(alphas, magnetization_perp / N_PART, 'o', label=r'simulation results $\perp$')
plt.plot(alphas, magnetization_para / N_PART, 'o', label=r'simulation results $\parallel$')
plt.legend(fontsize=20)
plt.show()

We can see that the simulation results are better represented by the curves of Ref. [1] compared to the Langevin function. This was to be expected as the Langevin function is the magnetization curve of the real three dimensional system without dipole-dipole interaction. We can also see that the magnetization is smaller in the case of an external magnetic field perpendicular to the monolayer plane compared to the parallel case.

Feel free to experiment with different dipolar interaction parameters $\lambda$ up to around 4 and different area fractions $\phi$ up to around 0.4. For higher values the here used simple sampling method is not applicable as the particles form clusters of very high relaxation times exceeding normal simulation times by far. Therefore more advanced methods are necessary to increase the sampling performance.

It should also be noted that perhaps thereby one has to decrease the time step as for higher values of the dipolar interaction parameter and the area fraction the interaction between the particles is stronger.

References

[1] Tamás Kristóf and István Szalai. Magnetic properties in monolayers of a model polydisperse ferrofluid. Phys. Rev. E 72: 041105, 2005. DOI:10.1103/PhysRevE.72.041105