Author: Marijan Beg, Ryan Pepper
Date: 11 May 2016
This problem is to calculate the single domain limit of a cubic magnetic particle. This is the size $L$ of equal energy for the so-called flower state (which one may also call a splayed state or a modified single-domain state) on the one hand, and the vortex or curling state on the other hand.
Geometry:
A cube with edge length, $L$, expressed in units of the intrinsic length scale, $l_\text{ex} = \sqrt{A/K_\text{m}}$, where $K_\text{m}$ is a magnetostatic energy density, $K_\text{m} = \frac{1}{2}\mu_{0}M_\text{s}^{2}$.
Material parameters:
More details about the standard problem 3 can be found in Ref. 1.
In [1]:
!rm -rf standard_problem3/ # Delete old result files (if any).
Firstly, we import all necessary modules.
In [2]:
import sys
sys.path.append('../')
from sim import Sim
from atlases import BoxAtlas
from meshes import RectangularMesh
from energies.exchange import UniformExchange
from energies.demag import Demag
from energies.zeeman import FixedZeeman
from energies.anisotropy import UniaxialAnisotropy
The following two functions are used for initialising the system's magnetisation [1].
In [3]:
# Function for initiaising the flower state.
def m_init_flower(pos):
x, y, z = pos[0]/1e-9, pos[1]/1e-9, pos[2]/1e-9
mx = 0
my = 2*z - 1
mz = -2*y + 1
norm_squared = mx**2 + my**2 + mz**2
if norm_squared <= 0.05:
return (1, 0, 0)
else:
return (mx, my, mz)
# Function for initialising the vortex state.
def m_init_vortex(pos):
x, y, z = pos[0]/1e-9, pos[1]/1e-9, pos[2]/1e-9
mx = 0
my = np.sin(np.pi/2 * (x-0.5))
mz = np.cos(np.pi/2 * (x-0.5))
return (mx, my, mz)
The following function is used for convenience. It takes two arguments:
It returns the relaxed simulation object.
In [4]:
import numpy as np
def relaxed_state(L, m_init):
mu0 = 4*np.pi*1e-7 # magnetic constant (H/m)
N = 16 # discretisation in one dimension
cubesize = 100e-9 # cude edge length (m)
cellsize = cubesize/N # discretisation in all three dimensions.
lex = cubesize/L # exchange length.
Km = 1e6 # magnetostatic energy density (J/m**3)
Ms = np.sqrt(2*Km/mu0) # magnetisation saturation (A/m)
A = 0.5 * mu0 * Ms**2 * lex**2 # exchange energy constant
K1 = 0.1*Km # Uniaxial anisotropy constant
axis = (0, 0, 1) # Uniaxial anisotropy easy-axis
cmin = (0, 0, 0) # Minimum sample coordinate.
cmax = (cubesize, cubesize, cubesize) # Maximum sample coordinate.
d = (cellsize, cellsize, cellsize) # Discretisation.
atlas = BoxAtlas(cmin, cmax) # Create an atlas object.
mesh = RectangularMesh(atlas, d) # Create a mesh object.
sim = Sim(mesh, Ms, name='standard_problem3') # Create a simulation object.
sim.add(UniformExchange(A)) # Add exchange energy.
sim.add(Demag()) # Add demagnetisation energy.
sim.add(UniaxialAnisotropy(K1, axis)) # Add uniaxial anisotropy energy.
sim.set_m(m_init) # Initialise the system.
sim.relax() # Relax the magnetisation.
return sim
Vortex state:
In [5]:
sim_vortex = relaxed_state(8, m_init_vortex)
print 'The relaxed state energy is {} J'.format(sim_vortex.total_energy())
# Plot the magnetisation in the sample slice.
%matplotlib inline
sim_vortex.m.plot_slice('y', 50e-9, xsize=6)
Flower state:
In [6]:
sim_flower = relaxed_state(8, m_init_flower)
print 'The relaxed state energy is {} J'.format(sim_flower.total_energy())
# Plot the magnetisation in the sample slice.
sim_flower.m.plot_slice('z', 50e-9, xsize=6)
Now, we can plot the energies of both vortex and flower states as a function of cube edge length. This will give us an idea where the state transition occurrs.
In [7]:
L_array = np.linspace(8, 9, 11) # values of L for which the system is relaxed.
vortex_energies = []
flower_energies = []
for L in L_array:
sim_vortex = relaxed_state(L, m_init_vortex)
sim_flower = relaxed_state(L, m_init_flower)
vortex_energies.append(sim_vortex.total_energy())
flower_energies.append(sim_flower.total_energy())
# Plot the energy dependences.
import matplotlib.pyplot as plt
plt.plot(L_array, vortex_energies, 'o-', label='vortex')
plt.plot(L_array, flower_energies, 'o-', label='flower')
plt.xlabel('L (lex)')
plt.ylabel('E')
plt.grid()
plt.legend()
Out[7]:
We now know that the energy crossing occurrs between $8l_\text{ex}$ and $9l_\text{ex}$, so a bisection algorithm can be used to find the exact crossing.
In [8]:
from scipy.optimize import bisect
def energy_difference(L):
sim_vortex = relaxed_state(L, m_init_vortex)
sim_flower = relaxed_state(L, m_init_flower)
return sim_vortex.total_energy() - sim_flower.total_energy()
cross_section = bisect(energy_difference, 8, 9, xtol=0.1)
print 'The transition between vortex and flower states occurs at {}*lex'.format(cross_section)