2.4. Example: Hysteresis loop for Stoner-Wohlfarth particle

This example is very similar to Example: Simple hysteresis loop but computes the hysteresis loop of a smaller ellipsoidal magnetic object. This allows to compare the results with the analytical solution given by the Stoner-Wohlfarth model. We use an ellipsoid whose x,y,z semi-axes have lengths 9 nm, 3 nm and 3 nm, respectively. (The mesh is contained in ellipsoid.nmesh.h5 and produced with `Netgen`_ from ellipsoid.geo):

../_images/ellipsoid_mesh.png

To compute the hysteresis loop for the ellipsoid, we use the script ellipsoid.py:

import nmag
from nmag import SI, at

#create simulation object
sim = nmag.Simulation()

# define magnetic material
Py = nmag.MagMaterial(name="Py",
                      Ms=SI(1e6,"A/m"),
                      exchange_coupling=SI(13.0e-12, "J/m"))

# load mesh: the mesh dimensions are scaled by 0.5 nm
sim.load_mesh("ellipsoid.nmesh.h5",
              [("ellipsoid", Py)],
              unit_length=SI(1e-9,"m"))

# set initial magnetisation
sim.set_m([1.,1.,0.])

Hs = nmag.vector_set(direction=[1.,1.,0.],
                     norm_list=[1.0, 0.995, [], -1.0,
                                -0.995, -0.990, [], 1.0],
                     units=1e6*SI('A/m'))

# loop over the applied fields Hs
sim.hysteresis(Hs, save=[('averages', at('convergence'))])

We apply external magnetic fields in [110] direction (i.e. 45 degrees between the x and the y-axis) to this system, with strengths in the range of 1000 kA/m down to -1000 kA/m in steps of 5 kA/m.

The save parameter is used to tell the hysteresis command what data to save, and how often. Here, we are only interested in saving the spatially averaged magnetisation values for every stage (i.e. meta-stable equilibrium before the applied field is changed).

2.4.1. Plotting the hysteresis loop

To extract the data needed for plotting the hysteresis loop we proceed as explained in the previous example Example: Simple hysteresis loop. We use the ncol command and extract the data into a text file named plot.dat:

$ ncol ellipsoid H_ext_0 H_ext_1 H_ext_2 m_Py_0 m_Py_1 m_Py_2 > plot.dat

We then use `Gnuplot`_ to plot the loop:

$ gnuplot make_plot.gnu

The gnuplot script make_plot.gnu is:

set term postscript eps enhanced color
set out 'hysteresis.eps'
set xlabel 'Applied field (kA/m)'
set ylabel 'M / Ms'
versor_x = 1/sqrt(2)
versor_y = 1/sqrt(2)
versor_z = 0.0
scalar_prod(x1,x2,x3) = x1*versor_x + x2*versor_y + x3*versor_z

set mxtics 5            # minor tics and grid
set ytics 1
set mytics 5
set grid xtics ytics mxtics mytics lt -1 lw 0.5, lt 0
plot [-1050:1050] [-1.2:1.2] \
  'plot.dat' u (scalar_prod($1,$2,$3)/1000):(scalar_prod($4,$5,$6)) t 'Stoner-Wohlfarth' w lp 4

Note that within the gnuplot file, we project the magnetisation data in the [1,1,0] direction because the applied field was acting in this direction. We obtain this hysteresis loop:

../_images/hysteresis1.png

The coercive field, which is located somewhere between 165 and 170 kA/m, can now be compared with the analytically known result for this particular system. To compute it, we need the demagnetizing factors Nx, Ny, Nz of the particle along the main axes. Since we deal with a prolate ellipsoid where two of the axes have the same dimension (y and z in this case), it is sufficient to compute the factor along the longest axis (x axis). The other two are easily derived from the relation Nx + Ny + Nz = 1. The expression to compute Nx is

\[N_x = \frac{1}{m^2-1} \cdot \left[ \frac{m}{2\sqrt{m^2-1}} \cdot \ln\left( \frac{m+\sqrt{m^2-1}}{m-\sqrt{m^2-1}} \right) - 1 \right]\]

where we call the length of the x semi-axis a, the length of the y (or z) semi-axis c, and take m to be the ratio m = a/c. Here, the value of Nx is therefore 0.1087, so we have Ny = Nz = 0.4456. With these values the shape anisotropy is easily computed according to the expression:

\[H_a = M_s \cdot \Delta N = M_s \cdot \left(N_z-N_x\right)\]

This gives Ha = 337 kA/m in the case of Ms = 1000 kA/m. The final step is to compute the coercive field hc using this analytical (Stoner-Wohlfarth) result:

\[h_c = \frac{H_c}{H_a} = \sin \theta_0 \cdot \cos \theta_0\]

Here, theta_0 is the angle between the easy-axis of the particle (x-axis in our case) and the direction of the applied field. Substituting theta_0 = 45 (degrees) in the formula, we obtain hc = 0.5, that is Hc = 0.5 * Ha = 168 kA/m. As we have seen before, the simulated hysteresis loop gives a value between 165 and 170 kA/m, which is in agreement with the analytical solution.

Note that this simulation is relatively slow due to a number of constraints: to get good Stoner-Wolfarth behaviour, we need to describe the shape of the ellipsoid well, and thus need a small edgelength when we generate the mesh. We further need uniform behaviour of the magnetisation, which limits the overall size of the ellipsoid. A general property of micromagnetic simulations is that the associated differential equations get stiffer if the edge lengths (or more generally: distances between neighbouring degrees of freedom) become smaller. Stiffer systems of differential equations are harder to intergrate, and thus take more time.