# Can spin-lattice dynamics model hysteresis loops?

G. dos Santos  CONICET, Mendoza, 5500, Argentina  Facultad de IngenierÃa,
Universidad de Mendoza, Mendoza, 5500 Argentina<EMAIL_ADDRESS>F.
Romá  Departamento de Física, Universidad Nacional de San Luis, Instituto de
Física Aplicada (INFAP), Consejo Nacional de Investigaciones Científicas y
Técnicas (CONICET), Chacabuco 917, D5700BWS San Luis, Argentina  R. Aparicio
CONICET, Mendoza, 5500, Argentina  Facultad de Ingeniería, Universidad de
Mendoza, Mendoza, 5500 Argentina  J. Tranchida  CEA, DES/IRESNE/DEC, 13018
Saint Paul Lès Durance, France  E.M. Bringa  CONICET, Mendoza, 5500, Argentina
Facultad de IngenierÃa, Universidad de Mendoza, Mendoza, 5500 Argentina
Centro de Nanotecnología Aplicada, Facultad de Ciencias, Universidad Mayor,
Santiago, Chile 8580745

(August 27, 2024)

###### Abstract

We use spin-lattice dynamics (SLD) simulations to study the possibility of
modeling the low-temperature magnetic hysteresis behavior of a bulk
ferromagnetic material. The temporal evolution of the magnetic and mechanical
degrees of freedom is obtained through a set of two coupled Langevin
equations. Hysteresis loops are calculated for different angles between the
external field and the magnetocrystalline anisotropy axes. The influence of
several relevant parameters is studied, including field frequency, magnetic
damping, magnetic anisotropy (magnitude and type), magnetic exchange, and
system size. The role played by a moving lattice is also discussed. At very
low temperatures we find that the exchange, and lattice dynamics barely affect
the loops, while the field frequency and magnetic damping have a large effect
on it. The influence of anisotropy magnitude and symmetry is found to follow
the expected behavior. We show that a careful choice of simulations parameters
allows for an excellent agreement between SLD and the paradigmatic Stoner-
Wohlfarth model, opening the possibility for more accurate inclusion of
lattice defects and thermal effects.

††preprint: APS/123-QED

keywords: hysteresis, spin-lattice dynamics, ferromagnetic, Fe

##  I  Introduction

For basic research and applications, magnetic studies are of fundamental
relevance as they are an important tool for revealing otherwise hidden
structural, thermodynamic and physicochemical properties of a material.
Hysteresis loops, for instance, are usually measured to characterize the
dynamic behavior of bulk samples that exhibit microstructures like grain
boundaries, defects, etc, and multiple domains [1, 2]. Furthermore, such
measurements are also performed for low-dimensional magnets such as magnetic
nanoparticles (NPs) or thin films [3]. The power generated by a sample subject
to an alternating magnetic field can be directly determined from the shape of
a hysteresis loop. In fact, the Specific Absorption Rate (SAR), defined as the
absorbed energy per unit of mass, is proportional to the area of this curve
[4]. In magnetic hyperthermia [5, 6, 7], for instance, SAR provides the
heating efficiency of a specific type of magnetic NP.

There are different classical models that are used to describe magnetic
hysteresis [1, 8]. Perhaps, the most popular of them is that of Stoner-
Wohlfarth (SW) [9, 10]. This model is relatively simple to evaluate
numerically, allowing easy comparison with experimental results. However, the
model includes an important number of approximations. Among them, it is
considered that the material behaves as a single magnetic domain. This is
known as the macrospin approximation, which implies that below a certain
length scale, typically on the order a few nanometers, the internal structure
of the material is neglected. In the SW model the Hamiltonian is the sum of a
Zeeman term and an uniaxial anisotropy term. Hysteretic behavior arises as a
consequence of a simple zero-temperature dynamic protocol: as the external
field varies with time, the magnetization can either change continuously while
the system remains at an energy minimum, or it can jump abruptly to another
value when that minimum becomes unstable. In addition, although in the
original model thermal fluctuations were not considered, further approaches
have attempted to incorporate temperature effects [11]. Fortunately,
micromagnetic simulations allow to overcome many of these limitations [12,
13]. Thus, it is possible to study a system made up of a large number of
interacting macrospins with a Hamiltonian including any complex energetic
term. In particular, the equilibrium or nonequilibrium magnetization dynamics
of such model can be calculated by numerically solving the phenomenological
stochastic Landau-Lifshitz-Gilbert (sLLG) equation [14]. From first principles
it can be demonstrated that this same theoretical framework is valid to carry
out atomistic spin dynamic (ASD) simulations [15, 16]. This simulation method
has been successfully applied to the modeling of hysteresis loops [17, 18, 19,
20]. Another alternative to include temperature effects in hysteresis loops is
to use Metropolis Monte Carlo (MC) simulations for a fixed lattice [21, 22].

More recently, there has been increased interest in studying more realistic
models that take into account both spin and lattice degrees of freedom, and
the coupling between them [23, 24, 25, 26, 27]. These spin-lattice dynamics
(SLD) simulations can use a Langevin bath for the lattice and sLLG for the
spins [28]. Although this method is very powerful, its implementation can be
computationally demanding and, as the complexity of a system increases, the
number of parameters that characterize it also increases, and therefore the
practical use of this numerical scheme may not be feasible.

In this work we use SLD calculations with sLLG spin dynamics to simulate the
hysteresis behavior of a ferromagnetic material. We use physical parameters
that are typical for bulk bcc iron. We focus on determining optimal simulation
parameters that allow to obtain reliable results at low temperatures. Also, we
explore how the hysteresis loops depend on the different physical properties
that characterize these magnetic systems: anisotropy, exchange, damping, and
lattice vibrations. The paradigmatic SW model is taken as a reference to
achieve this goal. The paper is organized as follows. In Sec. II, the
theoretical models used and the simulation details are presented. In Sec. III,
the main results of this study are shown. Finally, a summary and the
conclusions of the work are drawn in Sec. IV. In order to simplify several
parametric studies, as a first approach we run spin dynamics simulations with
the positions of the atoms fixed at their ideal equilibrium values. The
influence of the spin-lattice coupling is explored and discussed in the last
part of Sec. III.

##  II  Methods

###  II.1  SLD Model

Let us consider an ensemble of $N$ atoms each endowed with a classical
magnetic moment. The Hamiltonian of the model is

$\mathcal{H}=\sum_{i=1}^{N}\frac{\left|\bm{p}_{i}\right|^{2}}{2m_{i}}+\sum_{i,j,i\neq
j}^{N}V(r_{ij})+\mathcal{H}_{\rm mag}.$ (1)

The first term in Eq. (1) accounts for the kinetic energy of the $N$ atoms of
the system where $\bm{p}_{i}$ and $m_{i}$ represent, respectively, the linear
momentum and mass of the $i$-th atom ($m_{i}=55.845$ u for iron). The second
term is the interatomic potential describing the interactions between pairs of
atoms separated by a distance $r_{ij}$. As usual in molecular dynamics
simulations of metals, we use a classical embedded atom model potential which
describes well a broad spectrum of iron properties [29]. The interatomic
cutoff distance for this potential was set to $0.57$ nm.

The coupling among the spin and lattice degrees of freedom is provided through
the last term in Eq. (1), a magnetic Hamiltonian defined as

$\mathcal{H}_{\rm
mag}=-\mu_{0}\mu\bm{H}\cdot\sum_{i=0}^{N}\bm{s}_{i}-\sum_{i,j,i\neq
j}^{N}J(r_{ij})\bm{s}_{i}\cdot\bm{s}_{j}+\mathcal{H}_{\rm ani}.$ (2)

Here, $\bm{s}_{i}$ is a classical unitary vector representing the spin of
$i$-th atom. The first term in Eq. (2) is the Zeeman energy, the interaction
of each spin with an external uniform magnetic field $\bm{H}$, where
$\mu=2.2\mu_{B}$ is the atomic magnetic moment for iron ($\mu_{B}$ is the Bohr
magneton) and $\mu_{0}$ is the vacuum permeability constant. Note that the
Zeeman energy within LAMMPS is defined slightly differently than in Eq. 2. The
second term is just a Heisenberg Hamiltonian describing the interaction
between spins, where $J(r_{ij})$ is an interatomic distance-dependent exchange
coupling which is defined as the following Bethe-Slater curve [30, 31],

$J(r_{ij})=4\alpha\left(\frac{r_{ij}}{\delta}\right)^{2}\left[1-\gamma\left(\frac{r_{ij}}{\delta}\right)^{2}\right]e^{-\left(\frac{r_{ij}}{\delta}\right)^{2}}\Theta(R_{c}-r_{ij}),$
(3)

being $\Theta(R_{c}-r_{ij})$ the Heaviside step function and $R_{c}$ the
cutoff distance. The coefficients in Eq. (3) must be fitted to ab-initio data.
In the present work, we have set $\alpha=25.498$ meV, $\gamma=0.281$, and
$\delta=0.1999$ nm which fits very well the ab-initio data from Pajda et. al
[32]. In addition, we choose $R_{c}=0.35$ nm.

Finally, the last term in Eq. (2) is responsible for computing the magneto-
crystalline anisotropy. In this work we consider one of two, either uniaxial
$\mathcal{H}_{\rm uni}$ or cubic $\mathcal{H}_{\rm cub}$ anisotropy. The
corresponding expressions are given in Eqs. (4) and (5):

$\mathcal{H}_{\rm uni}=-K_{1}\sum_{i=1}^{N}(\bm{s}_{i}\cdot\bm{n})^{2}$ (4)

and

$\begin{split}\mathcal{H}_{\rm
cub}=\sum_{i=1}^{N}K_{1}[(\bm{s}_{i}\cdot\bm{n}_{1})^{2}(\bm{s}_{i}\cdot\bm{n}_{2})^{2}+(\bm{s}_{i}\cdot\bm{n}_{2})^{2}(\bm{s}_{i}\cdot\bm{n}_{3})^{2}+\\\
(\bm{s}_{i}\cdot\bm{n}_{1})^{2}(\bm{s}_{i}\cdot\bm{n}_{3})^{2}]-K_{2}(\bm{s}_{i}\cdot\bm{n}_{1})^{2}(\bm{s}_{i}\cdot\bm{n}_{2})^{2}(\bm{s}_{i}\cdot\bm{n}_{3})^{2}.\end{split}$
(5)

Here, the unit vectors $\bm{n}_{1}$, $\bm{n}_{2}$, and $\bm{n}_{3}$ lie along
the three crystallographic directions $[100]$, $[010]$, and $[001]$,
respectively. In Eq. (4), $\bm{n}$ is also a unit vector that in general could
point along any of these axes. The first term in the cubic anisotropy energy
above is defined with different sign within LAMMPS. $K_{1}$ and $K_{2}$ are
the magneto crystalline anisotropy constants which we set $K_{1}=35$
$\mu$eV/atom and $K_{2V}=3.6$ $\mu$eV/atom (equivalents to volumetric
anisotropies $K_{1V}=470$ kJ/m3 and $K_{2V}=46$ kJ/m3). In order to make
comparisons, we have used the same value of $K_{1}$ in both anisotropy
equations. Since this constant is positive (and also $K_{2}>0$), then the easy
axes of magnetization in Eqs. (4) and (5) are given by the unit vector defined
above. Note that the values of $K_{1}$ and $K_{2}$ are ten times larger than
those usually used to model bulk bcc iron [3]. However, such larger anisotropy
magnitude has been considered for Fe NP, and we also consider bulk values for
selected runs. As we will discuss later, using a large anisotropy in the
simulations helps to quickly stabilize the magnetization of the system. In
addition, to speed up our calculations we have neglected long-range dipolar
interactions in the magnetic Hamiltonian. Since we simulate small systems (see
below), this approximation should not affect the validity of our results.

###  II.2  Simulation details

We perform SLD simulations using the SPIN package of the software LAMMPS [33,
34]. The evolution of the system is described by two coupled Lanvegin
equations, one for the spin and another one for the lattice degrees of
freedom. In particular, the former is the well known sLLG equation introduced
by Brown [14].

Each Langevin equation has a damping term and a random force (or torque) which
are connected through the “fluctuation-dissipation” theorem. We have chosen
the damping constants equal to $\lambda_{s}=0.5$ (for the spin degrees of
freedom) and $\lambda_{L}=1.0$ s-1 (for the lattice). We use separate Langevin
thermostats for the lattice and spin subsystems, but both are set to the same
temperature $T=10$ K. For a more detailed description of these SLD simulations
see Ref. [28].

The hysteresis loops were calculated applying an alternating magnetic field
$\bf{H}$, and averaging the curves over up to ten different cycles. The field
is varied discretely: remains constant during certain simulation time, and
then jumps to reach the value given by a function $H_{max}cos(2\pi ft)$.
$H_{max}$ is set to be larger than the expected saturation field, and jumps do
not have the same magnitude along the field, given that the simulations time
at each field value is kept constant.

The simulation time step was set to $0.1$ fs and each cycle took around $2-8$
ns, giving MHz frequencies. In particular, we have simulated field frequencies
of $f_{0}$, $f_{0}/2$, $f_{0}/4$, and $f_{0}/8$, with $f_{0}=500$ MHz. These
values are equivalent to sweeping rates of approximately $2.2\times 10^{9}$
$T/s$, $1.1\times 10^{9}$ $T/s$, $0.55\times 10^{9}$ $T/s$ and $0.275\times
10^{9}$ $T/s$ respectively. We note that, in general, in experiments
hysteresis loops are obtained using fields that change almost continuously,
and the measurement of the magnetic moment of a sample can take up to several
microseconds. Those time scales are well beyond the feasibility of SLD, or of
other simulation methods like ASD since typically, simulations only reach
nanosecond scale with sweep rates similar to ours [17, 19, 18]. Nevertheless,
as we show in the next section, the hysteresis loops quickly converge to a
limiting curve as the frequency is decreased.

Simulations were run for different angles, $\phi=0^{\circ}$, $45^{\circ}$, and
$90^{\circ}$, between the external field $\bf{H}$ and the easy anisotropy
axes. In practice, this was done by changing the uniaxial anisotropy axis
directions while keeping $\bf H$ aligned in the $[001]$ direction ($z$ axis).
For uniaxial anisotropy, for instance, for $\phi=90^{\circ}$ or $45^{\circ}$
the unit vector $\bm{n}$ was oriented along the $[100]$ or $[101]$ direction.
Most simulations where run with magnetic uniaxial anisotropy according to Eq.
4.

In all cases, we use cubic samples of dimensions $(10\times 10\times
10)a_{0}^{3}$ ($a_{0}=0.286$ nm is the lattice parameter of bcc Fe), with
periodic boundary conditions in all directions. The resulting system contains
$2000$ atoms. To show that this size is large enough to capture the main
hysteresis properties of the model, we also run a few simulations using a
larger system with $(15\times 15\times 15)a_{0}^{3}$ cells (6750 atoms), as
shown in the Supplemental Material (SM).

Finally, we have to mention that the total magnetization $M$ as well as all
components $M_{x}$, $M_{y}$ and $M_{z}$, will be expressed normalized to the
ideal Fe bulk saturation magnetization $M_{s}$. Saturation magnetization is
given by the maximum magnetic moment per unit volume. For bulk iron, in the
volume of a bcc unit cell ($a_{0}^{3}$) there are 2 spins with magnetic moment
$\mu=2.2\mu_{B}$ and, therefore $M_{s}=\frac{2\times
2.2\mu_{B}}{a_{0}^{3}}\approx 1730$ kA/m. In this way, $M=1.0$ means $M=1730$
kA/m.

##  III  Results and discussion

In this section, we first analyze how the computational parameters, such as
the simulation run time and the field frequency, affect the final outcome of
the calculations. The aim is to establish optimal parameters that allow
minimizing the computational cost without affecting the reliability of the
calculations. Next, we focus on studying how the hysteresis loops can change
when different physical properties like anisotropy, exchange interaction, and
damping parameter, are modified. In a first stage, we analyze these effects
without coupling the lattice vibrations to the the spins dynamics, i.e. we
consider for most cases spin dynamics simulations with the atoms fixed at
their ideal-lattice positions. The influence of coupling the lattice
vibrations is analyzed on the final part of this section, where we have run
full SLD simulations.

###  III.1 Effect of simulation time

We start analyzing the effect of the simulation run time $t_{\rm sim}$. This
quantity represents the interval during which the external field $\bm{H}$ is
held constant before increasing or decreasing it by a given amount. It is
basically the total time of simulation used to calculate each individual
discrete point of a hysteresis loop. After changing the external field, the
magnetic moments must relax to a new stable configuration and therefore
$t_{\rm sim}$ must be large enough to allow this process to occur. It is
important, then, to determine an optimal value of this quantity that allows to
achieve this goal.

In principle, one can imagine an individual spin precessing around the field
direction, with a frequency which increases with the field magnitude. For Fe,
the resulting Larmor frequency gives a period of 36 ps for a 1T field, and
several precession periods are needed to describe magnetization evolution.
However, for damped dynamics, the frequency remains the same but the spin will
spiral down towards the field direction, allowing shorter simulations [15].
Figure 1 shows the time evolution of the components of the normalized
magnetization along the three Cartesian axes, $M_{x}$, $M_{y}$, $M_{z}$,
corresponding to directions $[100]$, $[010]$, and $[001]$, respectively,
during the simulation of an entire loop for the case of $\phi=90^{\circ}$.
Unless otherwise stated, for this and subsequent calculations we use uniaxial
anisotropy. As it can be noticed by inspection of this figure, the applied
field (green dashed curve and right axis) is maintained at a constant value
for 90 ps before increasing or decreasing it to the next value, so $t_{\rm
sim}=90$ ps. We are interested here in the stabilization of the magnetization
along the field direction, $M_{z}$. During the first moments after the field
value is varied a fluctuation period is observed, specially in the region
around the switching field, before $M_{z}$ reaches a stable value
approximately after 45 ps. This is better appreciated in Fig. 1(b), where a
zoom in the region around the first magnetization switching is displayed. In
this figure, we have also included the average value of $M_{z}$, $\langle
M_{z}\rangle$, which is obtained in our simulations from the last 30 ps of
each step, where a well stabilized magnetization is observed. We note that
$t_{\rm sim}=90$ ps is similar to simulation times used in micromagnetic and
ASD calculations of hysteresis loops [35, 13, 36, 17, 19, 18]. We have checked
that longer simulation times do not significantly affect the final results, as
it is shown in Fig. S3 in the SM. It is crucial that our magnetization
dynamics is well described near the switching field values. For fields of 0.5
T, the precession period is $\sim 72$ ps, and one would need a few ns of
simulation time for low damping. However, thanks to the high damping values
discussed below, shorter simulation times can be employed.

Figure 1: Time evolution of the components of the magnetization throughout an
entire cycle (a) and half a cycle (b). The field value at each time is also
included (green dashed line and right axis). Note that each field value is
kept constant for 90 ps ($9\times 10^{5}$ steps) of simulation time. In (b),
the average value of $M_{z}$, $\langle M_{z}\rangle$ (blue dotted line),
obtained from the last 30 ps of simulation, is also added for comparison.
Figure 2: Top images: resulting hysteresis loops obtained at the converged
frequency for $\phi=0^{\circ}$ (top left) and $\phi=90^{\circ}$ (top right).
Bottom images: snapshots of a fraction of the system, showing typical spin
configurations at different stages of the hysteresis loop for the case of
$\phi=0^{\circ}$ (bottom left) and $\phi=90^{\circ}$ (bottom right). Each
snapshot, (a), (b), (c) and (d), correspond to the points marked on the
hysteresis curves (upper images). In the snapshots, the atoms are represented
as red spheres and the spins are colored according to their orientation along
the field direction ($S_{z}$). Each cubic cell represents a fraction of the
system with a volume of (8 Å)3, extracted from the center of the original
system. Figure 3: Histograms of individual spins orientation along the three
axes. The histograms correspond to the state of the system described in the
snapshot (c) of Fig. 2 for the case of $\phi=90^{\circ}$.

Using $t_{sim}=90$ ps and frequency $f_{0}/4$, we calculate the hysteresis
loops for the cases of $\phi=0^{\circ}$ and $\phi=90^{\circ}$. These curves,
along with snapshots showing typical spin configurations at different stages
of the process, are shown in Fig. 2. As we can see, the system behaves
qualitatively according to what the SW model predicts, i.e. all spins are
roughly in sync like a single macrospin [9, 10]. This is expected for low-
temperature simulations, but deviations would occur at higher temperatures.
Still, Fig. 3 shows that there is a roughly Gaussian distribution of spin
values, an appreciable deviation from the simpler macrospin assumption. The
spin orientation histograms in this figure correspond to the state of the
system described by the snapshot (c) of the $\phi=90^{\circ}$ case of Fig. 2
(bottom right).

Another observation is the behavior of the magnetization components around
$H\approx 0$ (being $H$ the modulus of $\bf H$) for the case of
$\phi=90^{\circ}$. In Fig. 1 we can see that when this field takes values
close to zero, $M_{z}\approx 0$ (and also $M_{y}\approx 0$) while $M_{x}$
takes values close to saturation. This means that, at this stage, all spins
are aligned approximately along the $x$ direction (the easy axis of
magnetization for the uniaxial anisotropy), showing that the simulation
reproduces well the expected behavior for $\phi=90^{\circ}$. The snapshots in
Fig. 2 (right panels) confirm this interpretation.

Since these tests confirm that our simulations are consistent with some
general physical properties of the system, from now on we set $t_{\rm sim}=90$
ps.

###  III.2 Hysteresis loops and the effect of field frequency

As it was argued in Section II, using SLD simulations it is not possible to
calculate hysteresis loops at typical very low experimental frequencies, due
to intrinsic time limitations of atomistic simulation approaches. Instead, we
are limited to working with high frequencies in the order of MHz. However,
high frequency hysteresis loops appear to converge to a limit curve as the
frequency is decreased and, as the following analysis shows, this limit loop
should not be so different from those measured experimentally or calculated
theoretically. Following this line of argumentation, we calculate the
hysteresis loops for different frequencies. In Fig. 4, we show how the area of
the hysteresis loops depends on the field frequency for the cases with
$\phi=0^{\circ}$ and $\phi=45^{\circ}$. Note that we exclude the case
$\phi=90^{\circ}$ from this analysis, since for this angle the curves do not
present a hysteretic behavior. The curves in Fig. 4 exhibit a typical
crossover from a high dynamic regime (for which the area of the loops is
large) to a low dynamic regime (a slow variation of the area) [37].

Given that the simulation time for a given field value is always the same and,
in principle, long enough to cover spin relaxation for high damping, higher
frequency would only mean larger changes in the discrete field values. The
larger the field jumps, the larger the change across the energy landscape, and
the system might not adapt fast enough within the 90 ps of simulated time,
decreasing spin switching probability and generating wider loops. The
magnitude of the jumps around the region of the switching field are $\sim 0.3$
T for $f_{0}$ and $\sim 0.06$ T for $f_{0}/4$. Therefore, the stabilization of
the frequency for $\phi=0^{\circ}$ occurs when the field jumps are about one
order of magnitude lower than the value of the converged coercive field.

For $\phi=45^{\circ}$, we found convergence at a frequency of $f_{0}/2$ while
for the $\phi=0^{\circ}$ case it is not until a frequency of $f_{0}/4$ that
the loop area reaches a stable value. At these frequencies the area of the
loops coincides (within the statistical errors) with the predictions of the SW
model. Figure S2 in the SM shows the corresponding hysteresis curves
calculated at several different frequencies for the case of $\phi=0^{\circ}$,
evidencing convergence as frequency goes down. This is similar to results
showing convergence as loop time increases in ASD simulations [38] for
simulation times similar to our SLD simulations. This analysis shows that,
although the convergence frequencies found in this study are much higher than
those typicallly used in experiments, they are low enough to allow the
magnetization to equilibrate with the field direction at each field step. This
means that, at this rates, the ballistic regime is avoided and, therefore, our
simulated hysteresis loops should not present large discrepancies with
experiments, as discussed for example by Westmoreland et. al [19].

Figure 4: Loop area vs field frequency for the cases with $\phi=0^{\circ}$ and
$\phi=45^{\circ}$. The values of the area are normalized by the area of the SW
loop.

In Fig. 5 we compare the hysteresis loops calculated for $\phi=0^{\circ}$,
$\phi=45^{\circ}$, and $\phi=90^{\circ}$ at a single frequency $f_{0}/4$, with
those of the SW model [3]. As we can see, a good agreement is observed in all
cases. Nevertheless, it is important to stress that we do not expect that the
simulation curves to fit exactly to the theoretical ones. Unlike the SW model
(for which the macrospin approximation is considered at $T=0$), the SDL
simulations are carried out at low but finite temperature for many magnetic
atomic moments which, along a hysteresis cycle, do not rotate coherently in a
perfect way (see Fig. 3). For example, for the $\phi=0^{\circ}$ case, we do
not obtain a rectangular curve because thermal fluctuations and also the
existence of many degrees of freedom, make it easier for the magnetization to
change its orientation at a external field value slightly lower than predicted
by the SW model.

Figure 5: Hysteresis curves calculated with SLD simulations (full lines) at
different angles $\phi$ compared to the theoretical ones according to the SW
model (dashed lines). The SLD curves correspond to the average curve obtained
from many different cycles. In all cases a field frequencies of $f=f_{0}/4$
was used. Figure 6: (a) Loops for different values of uniaxial anisotropy
$K_{1}$. (b) Loop area normalized to the area of the loop with $K_{1}^{*}=35$
$\mu$eV/atom, versus $K_{1}/K_{1}^{*}$. The dashed line is the prediction of
the SW model.

###  III.3 Effect of anisotropy

According to the SW model, the coercive field and, therefore, the area of the
hysteresis loop when the external field is applied along the easy axis
($\phi=0^{\circ}$), is proportional to the anisotropy constant $K_{1}$ in Eq.
(4). We have explored this dependency in our SLD simulations by reducing and
increasing the original anisotropy constant value (throughout this test called
$K_{1}^{*}=35$ $\mu$eV/atom) by a factor of 10. Figure 6 (a) shows that the
resulting simulated loops qualitatively exhibit the expected behavior, that
is, the loop broadens (narrows) according to the increase (reduction) of the
anisotropy constant. This is further confirmed in Fig. 6 (b) where it can be
seen that the loop area follows an approximately linear behavior as a function
of $K_{1}/K_{1}^{*}$, deviating slightly from the SW model prediction only for
small values of this constant. These deviations arise from the combination of
two effects. On one hand, this is due to the fact that simulated loops are
“rounded” near the point of spin reversal, deviating from perfect rectangular
loops as explained above. The discrepancy with the SW model area become
increasingly important as the loop narrows. On the other hand, as it was
mentioned in the methods section, large anisotropy helps stabilizing the
magnetization of the system [18]. For the case of low anisotropy
($K_{1}=0.1K_{1}^{*}$) there is poor stabilization of $M_{z}$ during the 90 ps
of simulation of each field step (specially near the region of magnetization
switching), as seen in Figure S6. Much longer simulation times would be
required for low anisotropy values. Given that our smallest anisotropy is the
Fe bulk anisotropy, we could compare with some experimental results. MOKE
measurements of Fe(100) appear to give an slightly lower coercive field than
our simulation estimate [39], although those measurements were obtained at
higher temperatures and for multi-domain samples.

We have also analyzed the influence of anisotropy symmetry. Loops calculated
with SLD simulations for systems with uniaxial and cubic anisotropies, Eqs. 4
and 5, are shown in Fig. 7 for the case $\phi=0^{\circ}$. We compare these
curves with the corresponding one for the SW model with uniaxial anisotropy.
The reason is that at $T=0$ K both simulation curves should coincide with the
latter. Our results show that for the case of $\phi=0^{\circ}$ the hysteresis
loop does not depend on the type of anisotropy. Therefore, the small
differences between the curves shown in Fig. 7 can be attributed to the effect
of temperature, which is more pronounced in the case of cubic anisotropy. The
area of the loop is slightly smaller for cubic anisotropy, similar to the
results reported in ref. [35]. Usov and Peschany [40] considered the SW model
for uniaxial and cubic anisotropies, for randomly oriented nanoparticles. They
report that the maximum normalized coercive field in both cases is 1, in
agreement with our result. For a random collection of anisotropy orientations,
the resulting coercive field for cubic anisotropy was $\sim 0.68$ of the one
for uniaxial.

Figure 7: Hysteresis loops for uniaxial (red full line) and cubic (blue
circles) anisotropies for the $\phi=0^{\circ}$ case and frequency
$\nicefrac{{f_{0}}}{{4}}$. The results are compared with the SW prediction for
uniaxial anisotropy (red dashed line).

It is more interesting to analyze the difference between uniaxial and cubic
anisotropies, when the field is applied in the $[111]$ direction. For this
field orientation, we consider a SW model whose Hamiltonian has an arbitrary
anisotropy term. We start by applying an external field strong enough such
that there is a single energy minimum for this Hamiltonian, for which the
magnetization $\bm{M}$ points in the direction of $\bm{H}$. Then, we decrease
the intensity of the field using small jumps. At each step we use a steepest
descent method to try to escape from the energy minimum, so that if it becomes
unstable, a new minimum can be reached. Using this simple algorithm, we have
simulated the zero temperature dynamics of the SW model with both uniaxial and
cubic anisotropies. As shown in Fig. 8, in this case the corresponding
hysteresis loops for the SW model are very different. However, as expected,
the loops calculated with the SLD simulations agrees very well with the curves
for the SW model for both, cubic and uniaxial anisotropies. In the case of
cubic anisotropy, while the remanent magnetization for both curves (SW and
SLD) is $\pm 1/\sqrt{3}$, there is a qualitative difference between them for
field values greater than the coercive one: a small hysteresis behavior
present in the loop of the SW model is not well reproduced in the SLD
simulation. The reason for this is that the energy barrier separating these
states is of approximately $0.3$ $\mu$eV/atom, several orders of magnitude
smaller than $K_{B}T\sim 8.6\times 10^{-3}$ eV, where $K_{B}$ is the Boltzmann
constant. Only a simulation at much lower temperature would be able to
reproduce this behavior. The cubic case has a lower coercive field, which
could be somewhat expected from results for a collection of NP with random
orientation [40].

Figure 8: Hysteresis loops for an external field applied along the [111]
direction ($\phi\approx 54.7^{\circ}$) and frequency
$\nicefrac{{f_{0}}}{{4}}$, for cubic and uniaxial anisotropy. The loops
obtained with SLD are compared to the predictions of an extension of the SW
model that incorporates cubic anisotropy (red short-dotted line), and to the
SW with uniaxial anisotropy (green dot-dashed line).

###  III.4 Effect of exchange interactions

Exchange energy does not play a role in SW, and within that framework there
should be no effect from the exchange interaction on the simulated hysteresis
loops. To test this within the SLD framework, we have run simulations with
different values of the function $J(r_{ij})$. In Fig. 9, we compare the curves
obtained using the original exchange $J(r_{ij})^{*}$, with the ones obtained
by increasing or reducing this interaction by a factor of 10. As it can be
seen, loops are nearly unaffected by these changes. Increasing the magnitude
of the exchange interaction is expected to give higher values of magnetization
and broader loops [41], but for low temperatures (as in this case) this effect
is expected to be very small. Larger differences should be noticed at higher
temperatures where thermal fluctuations start to play a major role in the
magnetization behavior. Nevertheless, even at this low temperatures, a few
subtle effects can be noticed: (i) for the simulation with
$J(r_{ij})=0.1\times J(r_{ij})^{*}$ a small but consistent reduction of
saturation magnetization is observed, and (ii) for $J(r_{ij})=10\times
J(r_{ij})^{*}$ the loop becomes more rectangular. This last effect is also
observed in other theoretical approaches based on the SW model which include
exchange interactions through a mean field approach as in Refs. [42, 43].
Lower/higher exchange will lead to lower/higher Curie temperature and spin
ordering will be modified. For the same value of external field, the histogram
of spin values is narrower and with a higher mean value for the higher
exchange, facilitating the macrospin behavior, leading to a larger saturation
field and helping with the sudden spin flip which causes a more square loop.

Figure 9: Simulated hysteresis curves for the case of $\phi=0^{\circ}$
considering different magnitudes of the exchange coupling.

###  III.5 Effect of damping parameter

The value of the Gilbert damping parameter in the sLLG equation is usually
chosen to be in the range $0.01-1$. Experimental estimates for Fe are around
0.001. However, the possibility of increasing typical spin damping values for
numerical convenience is often explored. In SLD simulations a value of
$\lambda_{s}=0.1$ was used in ref. [26]. ASD [44] and micromagnetic [35, 13]
simulations use values in the range $0.1-1$.

Figure 10: (a) Hysteresis curves for the case of $\phi=90^{\circ}$ with a
field frequency of $f_{0}$ and with damping parameters of $\lambda_{s}=0.1$
(dashed line) and $\lambda_{s}=0.5$ (full line). The plot shows several cycles
for each case. (b) Time evolution of the components of the magnetization
throughout the simulation of half a cycle for the cases considered in (a). The
field value at each time is also included (green dashed line and right axis).
Note that each field value is kept constant for 90 ps ($9\times 10^{5}$ steps)
of simulation time.

In general, as the damping increases, there is a faster energy loss, and the
relaxation time required for the system to reach a steady state decreases,
with spins aligned with the preferred orientation determined by the external
field and the anisotropy. Therefore, it is possible to use higher frequencies
to calculate the hysteresis loops, which become narrower as damping grows. In
the same spirit, in most of the simulations presented here we use a larger
damping value, $\lambda_{s}=0.5$. Using lower damping is possible but would
require lower frequencies as discussed below. The parameter $\lambda_{s}$ is
related to the optimal value of the simulation run time $t_{\rm sim}$.

In micromagnetic simulations, it has been proposed to use both very high
damping and high loop sweep rates (SR), in order to achieve computational
efficiency [13]. The rationale behind this proposal was as follows. The
magnitude of the coercive field is a function of the measurement time, and the
attempt frequency for spin flips. This measurement time can be related to SR,
and the flip frequency can be assumed proportional to damping. Therefore,
loops calculated with the same value of SR$/\lambda_{s}$ would have the same
coercive field which gives the loop width. Based on this, they were able to
use high SR in their simulations to obtain hysteresis loops that closely match
those obtained at very low SR, by employing extremely high damping values, and
large time steps. Recently, this methodology was successfully applied to
simulate hysteresis of iron oxide magnetic nanoparticles with application to
hyperthermia [36].

In Fig. 10 we use $f_{0}$ for the case with $\phi=90^{\circ}$. For
$\lambda_{s}=0.5$ the hysteresis loop in (a) is reasonably close to the SW
prediction, and the magnetization versus time in (b) indicates reasonable
convergence. However, for $\lambda_{s}=0.1$, the loop in (a) is poorly
defined, as expected from the lack of convergence towards stable magnetization
values shown at the top of panel (b). The possibility of reducing the
simulation time in micromagnetic simulations by increasing both, the SR and
the damping, leaving constant the ratio SR$/\lambda_{s}$, would be equivalent,
in our case, to keep constant the ratio $f/\lambda_{s}$. However, SLD
simulations at high damping might depart from the desired dynamics and should
be treated with care. In our simulations, by choosing $\lambda_{s}=0.5$, we
were able to set $t_{\rm sim}=90$ ps and use $f_{0}/4$ as discussed at the
beginning of this section. For lower values of damping, smaller frequencies
would be required, significantly increasing the computational costs. Following
the previous scaling, instance, for $\lambda_{s}=0.1$, around 400 ps ($4\times
10^{6}$ steps) are required to stabilize magnetization at each field value, in
contrasts with the 90 ps ($9\times 10^{5}$ steps) that are required for
$\lambda_{s}=0.5$. We have run some simulations for such a low damping and
nanosecond steps, and verified that this indeed the case, explaining the lack
of convergence in Fig. 10(a).

Figure 11: Comparison of the hysteresis loops obtained with (lattice moving)
and without (lattice frozen) the coupling to the lattice dynamics. All the
figures correspond the case of $\phi=0^{\circ}$ at different field
frequencies, $f_{0}$ (top), $f_{0}/2$ (middle) and $f_{0}/4$ (bottom). The
inset on the top figure shows histograms of the spins orientation along the
field direction ($S_{z}$) and correspond to the state of the system just
before the spin flip, marked with circles in the respective curves.

###  III.6 Effect of coupling the lattice dynamics

All the results presented up to this point correspond to frozen-lattice
simulations, meaning that the lattice vibrations are not included in the
dynamics. We now analyze the effect of full SDL calculations actually coupling
the spin and lattice degrees of freedom, which we refer as moving-lattice
simulations. In Fig. 11 we compare the simulated hysteresis loops for the case
of $\phi=0^{\circ}$ obtained with (moving-lattice) and without (frozen-
lattice) spin-lattice coupling at different frequencies. The former gives
lower coercivity and narrower loops in comparison with the frozen-lattice
case. This effect is more pronounced for higher frequencies. This is
consistent with the fact that in a system that incorporates new degrees of
freedom (in this case those of the lattice), the energy barrier that must be
overcome to reverse the magnetization must decrease. In other words, when a
reversal field is applied, a moving lattice provides additional paths for the
magnetization relaxation process. However this effect is small at a frequency
$f_{0}/4$ where stable magnetization is easily achieved for the chosen sweep
rate, showing that lattice dynamics contributions are not very important at
this low temperatures ($T=10$ K), as expected. Larger differences are expected
as the temperature increases. As it was previously shown for iron NP at
equilibrium, these contributions would become significant at temperatures
higher than 300 K [26, 45]. The equilibrium magnetization for a moving lattice
is always smaller than that of a frozen lattice for all temperatures [26]. In
addition, in the present approach the anisotropy contributions do not depend
on the lattice. Larger differences between moving and frozen lattice
simulations might be observed for Hamiltonians that include anisotropy terms
which are sensitive to lattice details, like orientation, vibrations, and
atomic volume [46, 47].

##  IV  Summary and Conclusions

We use Spin-Lattice Dynamic (SLD) simulations to calculate the hysteresis
loops of a bulk ferromagnetic system at a low temperature. SLD can include
temperature effects and lattice defects in a relatively simple way,
complementing spin-dynamic simulations and micromagnetic simulations.

In order to validate our SLD results, we compare to the Stoner-Wolfhart (SW)
model, which is widely used in experimental and simulation studies. In the SW
model, there is no lattice, no temperature, and no exchange interactions. It
assumes uniaxial anysotropy and Zeeman energy contributions, and spins in the
volume are assumed to behave like a single macrospin. We present comparisons
for different angles between the external field and the anisotropy axis. We
analyze the effects of several parameters like simulation run time, field
frequency, damping, exchange, anisotropy and lattice-vibrations. We find that
loops obtained with SLD agree very well with those calculated with SW,
provided that several parameters are chosen with care.

Given the short timestep required to integrate the atomic degrees of freedom,
usually $\sim 0.1$ fs - 1 fs, one needs relatively long simulations to achieve
a stable magnetization at a given value of the applied magnetic field.
Typically, this would involve several precession periods. Using the Larmor
frequency for Fe gives a period of $\sim 36$ ps for a field of $1$T. In order
to speed-up spin-dynamics convergence, a Gilbert damping $\lambda_{s}=1.0$ is
used in most ASD simulations [17, 19, 18, 20] and values much larger than 1
have been used in micromagnetic simulations of hysteresis loops [13]. In our
simulations, we find that for a given damping, as frequency decreases, the
loop area also decreases until reaching a constant value. For bulk iron, and a
Gilbert damping of $0.5$, we find that frequencies of $125$ MHz or lower
provide such a constant loop area, for simulation times of $90$ ps for each
field value. Values between 30 and 200 ps (depending on the system under study
and the Gilbert damping employed) are also used for ASD simulations [17, 18,
20]. For the case where the anisotropy axis and the field direction are
perpendicular, $\phi=90^{\circ}$, a converged loop can be found at even higher
frequencies. Agreement between SLD and SW suggests that loops using this
limiting frequency would not be much different from the experimental or
theoretical estimates. This result partially solves the problem related to the
computational costs in atomistic simulations of hysteresis loops, allowing to
reproduce experimental hysteresis loops, running SLD simulations at much
higher frequencies than those typically used in experiments.

From our simulations at a temperature of $10$ K, spins behave almost like a
macrospin, but there is some expected spread in magnetization values,
associated with thermal effects and exchange interactions. The exchange
energy, which is much larger than Zeeman or anisotropy energies, does not
significantly affect hysteresis loops at low temperature. The coercive field
depends linearly on the magnitude of the uniaxial anisotropy, as in the SW
model. Cubic anisotropy and uniaxial anisotropy produce nearly the same loops
for the external field aligned with an easy anisotropy axis. However, there
are important differences for the case of $\phi\approx 54.7^{\circ}$
misalignment, where the applied field lies along the [111] direction. For this
case, cubic anisotropy leads to a lower coercivity value.

When the atoms vibrations are coupled to spin dynamics, we find that lattice
dynamics contributions depend on the loop frequency, but they are very small,
as expected for such a low temperature. SLD simulations of hysteresis loops at
higher temperature, or including anisotropy terms sensitive to lattice details
[46, 47], are expected to behave differently than the simulations presented
here.

Future SLD simulations can be applied to hysteresis loops for systems of
technological interest, including nanoparticles [48] and bulk systems with
defects [49, 50, 51, 20, 45], providing information which can inform
micromagnetic and atomistic spin dynamics simulations. Large computer cluster
allow for simulations of increasingly larger systems by distributing spatial
scales, but time scale cannot be split amongst different parallel processes.
Therefore, new approaches for accelerated dynamics [52], bypassing this time
scale problem, will be required for more efficient simulations of hysteresis
loops.

###### Acknowledgements.

GDS and EMB thanks support from a SIIP-UNCUYO-2022-2023 grant, and from PICTO-
UUMM-2019-00048. We thank computer run time in the cluster TOKO
(toko.uncu.edu.ar). The authors thank support from IPAC-2019 grant from
Sistema Nacional de Computación de Alto Desempeño (SNCAD), for run time in the
cluster Dirac http://dirac.df.uba.ar/. FR acknowledges financial support from
CONICET (Argentina) under Project No. PIP 112-202001-01294-CO and Universidad
Nacional de San Luis (Argentina) under Project No. PROICO 03-2220.

## References

  * Bertotti [1998] G. Bertotti,  _Hysteresis in Magnetism_ (Academic Press, Elsevier, London, 1998).
  * Mayergoyz [2003] I. Mayergoyz,  _Mathematical Models of Hysteresis and their Applications_ (Academic Press, Elsevier, London, 2003).
  * Cullity and Graham [2011] B. D. Cullity and C. D. Graham,  _Introduction to magnetic materials_ (John Wiley & Sons, New Jersey, 2011).
  * Mehdaoui  _et al._ [2011] B. Mehdaoui, A. Meffre, J. Carrey, S. Lachaize, L.-M. Lacroix, M. Gougeon, B. Chaudret, and M. Respaud, Optimal size of nanoparticles for magnetic hyperthermia: a combined theoretical and experimental study, Advanced Functional Materials 21, 4573 (2011).
  * Shaterabadi  _et al._ [2018] Z. Shaterabadi, G. Nabiyouni, and M. Soleymani, Physics responsible for heating efficiency and self-controlled temperature rise of magnetic nanoparticles in magnetic hyperthermia therapy, Progress in biophysics and molecular biology 133, 9 (2018).
  * Pankhurst  _et al._ [2003] Q. A. Pankhurst, J. Connolly, S. K. Jones, and J. Dobson, Applications of magnetic nanoparticles in biomedicine, Journal of physics D: Applied physics 36, R167 (2003).
  * Hedayatnasab  _et al._ [2017] Z. Hedayatnasab, F. Abnisa, and W. M. A. W. Daud, Review on magnetic nanoparticles for magnetic nanofluid hyperthermia application, Materials & Design 123, 174 (2017).
  * Liorzou  _et al._ [2000] F. Liorzou, B. Phelps, and D. Atherton, Macroscopic models of magnetization, IEEE Transactions on Magnetics 36, 418 (2000).
  * Stoner and Wohlfarth [1948] E. C. Stoner and E. Wohlfarth, A mechanism of magnetic hysteresis in heterogeneous alloys, Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical Sciences 240, 599 (1948).
  * Tannous and Gieraltowski [2008] C. Tannous and J. Gieraltowski, The stoner–wohlfarth model of ferromagnetism, European journal of physics 29, 475 (2008).
  * Lanci and Kent [2003] L. Lanci and D. V. Kent, Introduction of thermal activation in forward modeling of hysteresis loops for single-domain magnetic particles and implications for the interpretation of the day diagram, Journal of Geophysical Research: Solid Earth 108 (2003).
  * Bertotti  _et al._ [2008] G. Bertotti, I. Mayergoyz, and C. Serpico,  _Nonlinear magnetization dynamics in nanosystems_ (Elsevier, Oxford, 2008).
  * Behbahani  _et al._ [2020] R. Behbahani, M. L. Plumer, and I. Saika-Voivod, Coarse-graining in micromagnetic simulations of dynamic hysteresis loops, Journal of Physics: Condensed Matter 32, 35LT01 (2020).
  * Brown [1963] W. F. Brown, Thermal fluctuations of a single-domain particle, Physical Review 130, 1677 (1963).
  * Evans  _et al._ [2014] R. F. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. Ellis, and R. W. Chantrell, Atomistic spin model simulations of magnetic nanomaterials, Journal of Physics: Condensed Matter 26, 103202 (2014).
  * Eriksson  _et al._ [2017] O. Eriksson, A. Bergman, L. Bergqvist, and J. Hellsvik,  _Atomistic Spin Dynamics_ (Oxford University Press, Oxford, 2017).
  * Böttcher  _et al._ [2011] D. Böttcher, A. Ernst, and J. Henk, Atomistic magnetization dynamics in nanostructures based on first principles calculations: application to co nanoislands on cu (111), Journal of Physics: Condensed Matter 23, 296003 (2011).
  * Alzate-Cardona  _et al._ [2019] J. Alzate-Cardona, D. Sabogal-Suárez, R. Evans, and E. Restrepo-Parra, Optimal phase space sampling for monte carlo simulations of heisenberg spin systems, Journal of Physics: Condensed Matter 31, 095802 (2019).
  * Westmoreland  _et al._ [2020] S. C. Westmoreland, C. Skelland, T. Shoji, M. Yano, A. Kato, M. Ito, G. Hrkac, T. Schrefl, R. F. Evans, and R. W. Chantrell, Atomistic simulations of $\alpha$-fe/nd2fe14b magnetic core/shell nanocomposites with enhanced energy product for high temperature permanent magnet applications, Journal of Applied Physics 127, 133901 (2020).
  * Jenkins  _et al._ [2021] S. Jenkins, R. W. Chantrell, and R. F. Evans, Atomistic origin of the athermal training effect in granular irmn/cofe bilayers, Physical Review B 103, 104419 (2021).
  * Tajiri  _et al._ [2018] T. Tajiri, H. Deguchi, M. Mito, K. Konishi, S. Miyahara, and A. Kohno, Effect of size on the magnetic properties and crystal structure of magnetically frustrated dymn 2 o 5 nanoparticles, Physical Review B 98, 064409 (2018).
  * Essajai  _et al._ [2019] R. Essajai, Y. Benhouria, A. Rachadi, M. Qjani, A. Mzerd, and N. Hassanain, Shape-dependent structural and magnetic properties of fe nanoparticles studied through simulation methods, RSC advances 9, 22057 (2019).
  * Shimada  _et al._ [2015] T. Shimada, K. Ouchi, I. Ikeda, Y. Ishii, and T. Kitamura, Magnetic instability criterion for spin–lattice systems, Computational Materials Science 97, 216 (2015).
  * Ma  _et al._ [2016] P.-W. Ma, S. Dudarev, and C. Woo, Spilady: A parallel cpu and gpu code for spin–lattice magnetic molecular dynamics simulations, Computer Physics Communications 207, 350 (2016).
  * Wu  _et al._ [2018] X. Wu, Z. Liu, and T. Luo, Magnon and phonon dispersion, lifetime, and thermal conductivity of iron from spin-lattice dynamics simulations, Journal of Applied Physics 123, 085109 (2018).
  * Dos Santos  _et al._ [2020] G. Dos Santos, R. Aparicio, D. Linares, E. Miranda, J. Tranchida, G. Pastor, and E. Bringa, Size-and temperature-dependent magnetization of iron nanoclusters, Physical Review B 102, 184426 (2020).
  * dos Santos  _et al._ [2021] G. dos Santos, R. Meyer, R. Aparicio, J. Tranchida, E. M. Bringa, and H. M. Urbassek, Spin-lattice dynamics of surface vs core magnetization in fe nanoparticles, Applied Physics Letters 119, 012404 (2021).
  * Tranchida  _et al._ [2018] J. Tranchida, S. Plimpton, P. Thibaudeau, and A. P. Thompson, Massively parallel symplectic algorithm for coupled magnetic spin dynamics and molecular dynamics, Journal of Computational Physics 372, 406 (2018).
  * Chamati  _et al._ [2006] H. Chamati, N. Papanicolaou, Y. Mishin, and D. Papaconstantopoulos, Embedded-atom potential for fe and its application to self-diffusion on fe (1 0 0), Surface Science 600, 1793 (2006).
  * Kaneyoshi [1992] T. Kaneyoshi,  _Introduction to amorphous magnets_ (World Scientific Publishing Company, Singapore, 1992).
  * Yosida  _et al._ [1996] K. Yosida, D. C. Mattis, and K. Yosida,  _THEORY OF MAGNETISM.: Edition en anglais_ , Vol. 122 (Springer Science & Business Media, Berlin, 1996).
  * Pajda  _et al._ [2001] M. Pajda, J. Kudrnovskỳ, I. Turek, V. Drchal, and P. Bruno, Ab initio calculations of exchange interactions, spin-wave stiffness constants, and curie temperatures of fe, co, and ni, Physical Review B 64, 174402 (2001).
  * Plimpton [1995] S. Plimpton, Fast parallel algorithms for short-range molecular dynamics, Journal of Computational Physics 117, 1 (1995).
  * [34] http://lammps.sandia.gov/doc/fix_nve_spin.html.
  * Aurélio and Vejpravova [2020] D. Aurélio and J. Vejpravova, Understanding magnetization dynamics of a magnetic nanoparticle with a disordered shell using micromagnetic simulations, Nanomaterials 10, 1149 (2020).
  * Behbahani  _et al._ [2021] R. Behbahani, M. L. Plumer, and I. Saika-Voivod, Multiscale modelling of magnetostatic effects on magnetic nanoparticles with application to hyperthermia, Journal of Physics: Condensed Matter 33, 215801 (2021).
  * Dolz  _et al._ [2020] M. I. Dolz, S. D. Calderón Rivero, H. Pastoriza, and F. Romá, Magnetic hysteresis behavior of granular manganite la$0.67$ca$0.33$mno$3$ nanotubes, Physcial Review B 101, 174425 (2020).
  * Chureemart  _et al._ [2017] P. Chureemart, R. Evans, R. Chantrell, P.-W. Huang, K. Wang, G. Ju, and J. Chureemart, Hybrid design for advanced magnetic recording media: Combining exchange-coupled composite media with coupled granular continuous media, Physical Review Applied 8, 024016 (2017).
  * Bansmann  _et al._ [1992] J. Bansmann, M. Getzlaff, C. Westphal, and G. Schönhense, Surface hysteresis curves of fe (110) and fe (100) crystals in ultrahigh vacuum—evidence of adsorbate influences, Journal of magnetism and magnetic materials 117, 38 (1992).
  * Usov and Peschany [1997] N. Usov and S. Peschany, Theoretical hysteresis loops for single-domain particles with cubic anisotropy, Journal of magnetism and magnetic materials 174, 247 (1997).
  * Kachkachi  _et al._ [2000] H. Kachkachi, A. Ezzir, M. Nogues, and E. Tronc, Surface effects in nanoparticles: application to maghemite-fe o, The European Physical Journal B 14, 681 (2000).
  * Atherton and Beattie [1990] D. L. Atherton and J. Beattie, A mean field stoner-wohlfarth hysteresis model, IEEE transactions on magnetics 26, 3059 (1990).
  * Zhang  _et al._ [2003] H.-w. Zhang, S.-y. Zhang, B.-g. Shen, and H. Kronmüller, The magnetization behavior of nanocrystalline permanent magnets based on the stoner–wohlfarth model, Journal of magnetism and magnetic materials 260, 352 (2003).
  * Evans [2020] R. F. Evans, Atomistic spin dynamics, in  _Handbook of Materials Modeling: Applications: Current and Emerging Materials_ , edited by W. Andreoni and S. Yip (Springer, Cham, 2020) pp. 427–448.
  * Meyer  _et al._ [2022] R. Meyer, G. dos Santos, R. Aparicio, E. M. Bringa, and H. M. Urbassek, Influence of vacancies on the temperature-dependent magnetism of bulk fe: A spin-lattice dynamics approach, Computational Condensed Matter , e00662 (2022).
  * Nieves  _et al._ [2021] P. Nieves, J. Tranchida, S. Arapan, and D. Legut, Spin-lattice model for cubic crystals, Physical Review B 103, 094437 (2021).
  * Dednam  _et al._ [2022] W. Dednam, C. Sabater, A. E. Botha, E. B. Lombardi, J. Fernández-Rossier, and M. J. Caturla, Spin-lattice dynamics simulation of the einstein–de haas effect, Computational Materials Science 209, 111359 (2022).
  * Lak  _et al._ [2021] A. Lak, S. Disch, and P. Bender, Embracing defects and disorder in magnetic nanoparticles, Advanced Science 8, 2002682 (2021).
  * Palanisamy  _et al._ [2021] D. Palanisamy, A. Kovács, O. Hegde, R. E. Dunin-Borkowski, D. Raabe, T. Hickel, and B. Gault, Influence of crystalline defects on magnetic nanodomains in a rare-earth-free magnetocrystalline anisotropic alloy, Physical Review Materials 5, 064403 (2021).
  * Han  _et al._ [2021] L. Han, Z. Rao, I. R. Souza Filho, F. Maccari, Y. Wei, G. Wu, A. Ahmadian, X. Zhou, O. Gutfleisch, D. Ponge, _et al._ , Ultrastrong and ductile soft magnetic high-entropy alloys via coherent ordered nanoprecipitates, Advanced Materials 33, 2102139 (2021).
  * Bersweiler  _et al._ [2021] M. Bersweiler, E. P. Sinaga, I. Peral, N. Adachi, P. Bender, N.-J. Steinke, E. P. Gilbert, Y. Todaka, A. Michels, and Y. Oba, Revealing defect-induced spin disorder in nanocrystalline ni, Physical Review Materials 5, 044409 (2021).
  * Plimpton  _et al._ [2020] S. J. Plimpton, D. Perez, and A. F. Voter, Parallel algorithms for hyperdynamics and local hyperdynamics, The Journal of Chemical Physics 153, 054116 (2020).

## Supplementary Material

This supplementary material presents extra figures intended to expand the
discussion and provide extra evidence and support the points discussed in the
main text.

###  IV.1 Average curve

Figure S1 presents an example of how the resulting hysteresis loops where
obtained. The resulting curve was obtained by averaging several individual
cycles.

Figure S1: 4 different hysteresis loops for the $\phi=0^{\circ}$ case at
frequency $f=f_{0}$ plotted along the average curve.

###  IV.2 Effect of frequency and loop convergence

Figure S2 shows the hysteresis loops at different frequencies for the case
where the anisotropy axis and the direction of the external field form an
angle of $\phi=0^{\circ}$. The plot evidence that the loops narrow when the
frequency is decreased until the convergence frequency is found and the
resulting loop coincides with the Stoner Wohlfarth prediction (green curve).

Figure S2: (Color online) Hysteresis loops for the $\phi=0^{\circ}$ case at
different field frequencies $f_{0}$, $f_{0}/2$, $f_{0}/4$ and $f_{0}/8$,
compared to the rectangular loop from the SW model (green curve).

###  IV.3 Increased simulation time

As it is argued in the main text, a simulation time of $t_{sim}=90$ ps for
each field value is sufficient to produce reliable hysteresis loops. Figure S3
shows that increasing this time by a factor of approximately 2 has no
appreciable effect in the shape of the resulting curves.

Figure S3: Hysteresis loops for the $\phi=90^{\circ}$ case and frequency
$f=f_{0}/2$ with different simulations times ($t_{\rm sim}$) for each point of
the curve

###  IV.4 Size effects

Figure S4 compares two individual hysteresis loops obtained for a bcc lattice
in a cubic box of volume $V=(10a_{0})^{3}$ (the system under study) and a
larger system with $V=(15a_{0})^{3}$. The simulations correspond to the case
of $\phi=0^{\circ}$ and frequency $f_{0}/4$.

Figure S5 compares the fluctuations of $M_{z}$, the component of the
magnetization along the flied direction, during several field steps for three
different system volumes.

Both figures show that simulating a system with $V=(10a_{0})^{3}$ is
sufficient to avoid large size-effects.

Figure S4: Hysteresis cycles for the $\phi=0^{\circ}$ case and frequency
$f_{0}/4$ for different system sizes. Figure S5: Comparison of the
fluctuations of the component of the magnetization along the applied field
direction $M_{z}$ during the time evolution of a simulation throughout half a
cycle for different system sizes. The simulations correspond in all cases to
the case of $\phi=0^{\circ}$ and frequency $f_{0}/4$. The field value at each
time is also included (green dashed line and right axis).

###  IV.5 Anisotropy and magnetization fluctuations

As it was argued in the Methods section of the main text, we have used mainly
uniaxial anisotropy with an anisotropy constant equal to $K_{1}^{*}=35$
$\mu$eV/atom, ten times larger than the Fe bulk value. The main reason for
this is that, using a large anisotropy in the simulations helps to quickly
stabilize the magnetization of the system for the $\phi=0^{\circ}$ case, and
therefore, less simulation time for each field value is required. Figure S6
shows this effect by comparing the time-evolution of the components of the
magnetization for two different anisotropy intensities.

Figure S6: Comparison of the fluctuations of the components of the
magnetization for the simulations with $\phi=0^{\circ}$, frequency $f_{0}/4$
and different anisotropy constants, during the same number of applied field
steps. The Gilbert damping in these simulations is $\lambda_{s}=0.5$.