We choose a cuboidal thin film permalloy sample measuring $120 \times 120 \times 10 \,\text{nm}^{3}$. The choice of a cuboid is important as it ensures that the finite difference method employed by OOMMF does not introduce errors due to irregular boundaries that cannot be discretised well. We choose the thin film geometry to be thin enough so that the variation of magnetisation dynamics along the out-of-film direction can be neglected. Material parameters based on permalloy are:
An external magnetic bias field with magnitude $80 \,\text{kA/m}$ is applied along the direction $e = (1, 0.715, 0)$. We choose the external magnetic field direction slightly off the sample diagonal in order to break the system’s symmetry and thus avoid degenerate eigenmodes. First, we initialize the system with a uniform out-of-plane magnetisation $m_{0} = (0, 0, 1)$. The system is allowed to relax for $5 \,\text{ns}$, which was found to be sufficient time to obtain a well-converged equilibrium magnetisation configuration. We refer to this stage of simulation as the relaxation stage, and its final relaxed magnetisation configuration is saved to serve as the initial configuration for the next dynamic stage. Because we want to use a well defined method that is supported by all simulation tools, we minimize the system’s energy by integrating the LLG equation with a large, quasistatic Gilbert damping $\alpha = 1$ for $5 \,\text{ns}$. In the next step (dynamic stage), a simulation is started using the equilibrium magnetisation configuration from the relaxation stage as the initial configuration. Now, the direction of an external magnetic field is altered to $e = (1, 0.7, 0)$. This simulation stage runs for $T = 20 \,\text{ns}$ while the (average and spatially resolved) magnetisation $M(t)$ is recorded every $\Delta t = 5 \,\text{ps}$. The Gilbert damping in this dynamic simulation stage is $\alpha = 0.008$.
Details of this standard problem specification can be found in Ref. 1.
Firstly, all required modules are imported.
import oommfc as oc
import discretisedfield as df
import micromagneticmodel as mm
Now, we specify all simulation parameters.
import numpy as np
lx = ly = 120e-9 # x and y dimensions of the sample(m)
lz = 10e-9 # sample thickness (m)
dx = dy = dz = 5e-9 # discretisation in x, y, and z directions (m)
Ms = 8e5 # saturation magnetisation (A/m)
A = 1.3e-11 # exchange energy constant (J/m)
H = 8e4 * np.array([0.81345856316858023, 0.58162287266553481, 0.0])
alpha = 0.008 # Gilbert damping
gamma0 = 2.211e5
Now, the system object can be created and mesh, magnetisation, hamiltonian, and dynamics are specified.
mesh = df.Mesh(p1=(0, 0, 0), p2=(lx, ly, lz), cell=(dx, dy, dz))
system = mm.System(name="stdprobfmr")
system.energy = mm.Exchange(A=A) + mm.Demag() + mm.Zeeman(H=H)
system.dynamics = mm.Precession(gamma0=gamma0) + mm.Damping(alpha=alpha)
system.m = df.Field(mesh, nvdim=3, value=(0, 0, 1), norm=Ms)
Finally, the system is relaxed.
md = oc.MinDriver()
md.drive(system)
Running OOMMF (ExeOOMMFRunner)[2023/10/23 16:04]... (0.4 s)
We can now load the relaxed state to the Field object and plot the $z$ slice of magnetisation.
system.m.sel("z").resample((10, 10)).mpl()
In the dynamic stage, we use the relaxed state from the relaxation stage.
# Change external magnetic field.
H = 8e4 * np.array([0.81923192051904048, 0.57346234436332832, 0.0])
system.energy.zeeman.H = H
Finally, we run the multiple stage simulation using TimeDriver
.
T = 20e-9
n = 4000
td = oc.TimeDriver()
td.drive(system, t=T, n=n)
Running OOMMF (ExeOOMMFRunner)[2023/10/23 16:04]... (47.8 s)
From the obtained vector field samples, we can compute the average of magnetisation $y$ component and plot its time evolution.
import matplotlib.pyplot as plt
t = system.table.data["t"].values
my = system.table.data["mx"].values
# Plot <my> time evolution.
plt.figure(figsize=(8, 6))
plt.plot(t, my)
plt.xlabel("t (ns)")
plt.ylabel("my average")
plt.grid()
From the $<m_{y}>$ time evolution, we can compute and plot its Fourier transform.
import scipy.fft
psd = np.log10(np.abs(scipy.fft.fft(my)) ** 2)
f_axis = scipy.fft.fftfreq(4000, d=20e-9 / 4000)
plt.plot(f_axis / 1e9, psd)
plt.xlim([6, 12])
plt.xlabel("f (GHz)")
plt.ylabel("Psa (a.u.)")
plt.grid()
[1] A. Baker, M. Beg, G. Ashton, M. Albert, D. Chernyshenko, W. Wang, S. Zhang, M.-A. Bisotti, M. Franchin, C.L. Hu, R. Stamps, T. Hesjedal, and H. Fangohr, J. Magn. Magn. Mater. 421, 428 (2017).