# Comparative study of magnetic properties of Mn3+ magnetic clusters in GaN
using classical and quantum mechanical approach

Y. K. Edathumkandy  Institute of Physics, Polish Academy of Sciences, PL
02-668 Warszawa, Poland  D. Sztenkiel<EMAIL_ADDRESS>Institute of
Physics, Polish Academy of Sciences, PL 02-668 Warszawa, Poland

(August 27, 2024)

###### Abstract

Currently, simulations of many-body quantum systems are known to be
computationally too demanding to be solved on classical computers. The main
problem is that the computation time and memory necessary for performing the
calculations usually grow exponentially with the number of particles $N$. An
efficient approach to simulate many-body quantum systems is the use of
classical approximation. However, it is known that at least at low
temperatures, the allowed spin fluctuations in this approach are overestimated
what results in enhanced thermal fluctuations. It is therefore timely and
important to assess the validity of the classical approximation. To this end,
in this work, we compare the results of numerical calculations of small Mn3+
magnetic clusters in GaN, where the Mn spins are treated classically with
those where they are treated quantum-mechanically (crystal field model). In
the first case, we solve the Landau-Lifshitz-Gilbert (LLG) equation that
describes the precessional dynamics of spins represented by classical vectors.
On the other hand, in the crystal field model, the state of Mn3+ ion ($d^{4}$
configuration with $S=2$, $L=2$) is characterized by the set of orbital and
spin quantum numbers $|m_{s},m_{L}>$. Particular attention is paid to use
numerical parameters that ensure the same single ion magnetic anisotropy in
both classical and quantum approximation. Finally, a detailed comparative
study of magnetization $\mathbf{M}(\mathbf{H},T)$ as a function of the
magnetic field $\mathbf{H}$, temperature $T$, number of ions in a given
cluster $N$ and the strength of super-exchange interaction $J$, obtained from
both approaches will be presented.

##  I Introduction

Simulations of quantum systems using classical computers bear the potential
for advancing scientific knowledge of properties of many materials and their
time evolution. However, such simulations may cause diverse problems.
Difficulties arise from the exponential growth of required resources
(computing time and computer memory) with the number of elements in the system
[1, 2, 3]. For example, the number of parameters required to describe the
state of $N$ interacting fermionic particles scales as $2^{N}$. Consequently,
the memory allotted to save the state increases exponentially with $N$
necessitating the use of vast gigabytes of data. Similarly, the computing time
required for performing the simulations is also subject to the same
exponential scaling with the number of particles, demanding, in some cases,
several years to complete the simulation. Even though approximation techniques
such as quantum Monte-Carlo methods [4] and renormalisation group algorithms
developed by Wilson [5] have been designed, they are not entirely successful
in the sufficient description of the state of the systems under investigation.
Fortunately, many-body quantum systems can be, in principle, simulated using
classical approximations. The Landau-Lifshitz-Gilbert (LLG) equation with
Langevin dynamics has been applied extensively to compute magnetic properties
of different materials and nanoparticles. With the significant improvement of
the computing power and efficiency of parallel algorithms, large-scale
atomistic spin model or micromagnetic simulations using the stochastic LLG
(sLLG) method gradually become possible. For example, ferromagnetic hysteresis
loops, spin waves, domain wall motion, and precessional or thermally assisted
magnetization switching processes have been investigated using this approach
[6, 7, 8, 9, 10, 11]. Classical approximations can reduce the overall
computational complexity, but they usually face difficulties in proper
reproduction of temperature-dependent properties of the system due to
neglecting of spin quantization. It is generally accepted that allowed spin
fluctuations in the classical approach are overestimated and an appropriate
re-scaling of simulation temperature is required to resolve this issue [6, 7].
It is important then to somehow assess the validity of the classical
approximation and its underlying assumptions. To accomplish this goal, we
compare magnetizations $\mathbf{M}(\mathbf{H},T)$ as a function of the
magnetic field ($\mathbf{H}$) and temperature ($T$) obtained by the sLLG
method with those obtained by the exact approach, that is by the quantum
mechanical crystal field model (CFM). Here we choose Ga1-xMnxN dilute magnetic
semiconductor, as this material has been extensively and thoroughly
investigated [12, 16, 17, 18, 19, 20, 21, 22, 14, 13, 15], and the
experimental findings has been remarkably well reproduced by the exact CFM
approach [16, 17, 19, 22]. Even though, there are many different cluster types
in GaN with randomly distributed Mn ions [23], in this article we consider
only a few exemplary ones. In particular, the simulations are restricted here
to very small magnetic clusters consisting of up to four Mn3+ ions coupled by
the ferromagnetic super-exchange interaction
$-J\mathbf{S}_{1}\cdot\mathbf{S}_{2}$, as shown in Fig. 1.

In this paper, we show that the comparison between classical and quantum
approaches improves with the increasing magnitude of the exchange coupling
parameter $J^{q}$ and it is the best for large number of ions $N$ and strong
magnitude of the exchange energy as compared to the thermal one, that is
$J^{q}|\mathbf{\hat{S}}_{i}|\cdot|\mathbf{\hat{S}}_{j}|\gg k_{B}T$. It is
shown, that it is possible to only partially explain the inconsistency between
classical and quantum results presented in temperature-dependent magnetization
curves by referring to the concept of quantum mechanical “stiffness”. On the
contrary, we find that the key factor is the value of effective spin of the
system. Depending on the relative strength between the thermal and the
exchange energy we deal with uncoupled ions characterized by $S=2$ or strongly
bound clusters with large total spin value $S_{T}=NS$. For weak coupling
between the spins
$J^{q}|\mathbf{\hat{S}}_{i}|\cdot|\mathbf{\hat{S}}_{j}|<k_{B}T$, the
investigated clusters consist of practically uncoupled spins. Then the
classical approximation underestimates the actual value of magnetization.
Contrary, in the strong coupling regime, where
$J^{q}|\mathbf{\hat{S}}_{i}|\cdot|\mathbf{\hat{S}}_{j}|\gg k_{B}T$, we deal
with one magnetic system characterized by high total spin value $S_{T}=NS$. In
this limit, $S_{T}$ increases with the number of ions $N$ in the cluster,
which improves the comparison between two approaches. These conclusions are
consistent with the common relation, that the classical and quantum results
differ for a small value of $S$ and both treatments merge for
$S\rightarrow\infty$.

Figure 1: Selected clusters composed of magnetic ions with a single exchange
constant. Lines indicate exchange bonds. Only clusters with up to four ions
are considered.

##  II MODEL

###  II.1 Crystal field model

When an ion is inserted in a solid, Coulomb interaction of its charge
distribution with the neighbouring charges in the crystal together with the
hybridization of relevant atomic orbitals should be taken into account. This
effect is called the crystal field interaction, which causes the quenching of
the orbital angular momentum and the development of the single-ion magnetic
anisotropy. By incorporating in the Hamiltonian additional terms corresponding
to the spin-orbit interaction and the influence of external magnetic field
$\mathbf{B}$, it is possible to simulate magnetic properties of transition
metal ions placed in a crystal environment. Such simulations, termed crystal
field model, were developed by Vallin and coworkers et al. to study the
magnetization of a Cr ion in II-VI dilute magnetic semiconductors [24]. Here
we use this model to obtain $\mathbf{M}(\mathbf{H},T)$ curves of magnetic
clusters composed of up to 4 Mn ions [22]. Detailed reviews on this
theoretical model are given in Ref. 25, 24, 26, 12, 16, 17, 22 and we follow
this approach. The Hamiltonian for an isolated Mn3+ impurity in GaN crystal is
of the form,

$\mathcal{H}^{q}(j)=\mathcal{H}^{q}_{CF}+\mathcal{H}^{q}_{JT}(j)+\mathcal{H}^{q}_{TR}+\mathcal{H}^{q}_{SO}+\mathcal{H}^{q}_{B}$
(1)

where
$\mathcal{H}^{q}_{CF}=-\frac{2}{3}B_{4}(\hat{O}_{4}^{0}-20\sqrt{2}\hat{O}_{4}^{3})$
is the cubic crystal field of tetrahedral symmetry,
$\mathcal{H}^{q}_{\textrm{JT}}=\tilde{B}_{2}^{0}\hat{\Theta}_{4}^{0}+\tilde{B}_{4}^{0}\hat{\Theta}_{4}^{2}$
describes the static Jahn-Teller distortion of tetragonal symmetry and
$\mathcal{H}^{q}_{TR}=B_{2}^{0}\hat{O}_{4}^{0}+B_{4}^{0}\hat{O}_{4}^{2}$
represents the trigonal distortion along the GaN hexagonal $\mathbb{c}$-axis.
This model incorporates also the spin orbit coupling
$\mathcal{H}^{q}_{S0}=\lambda\hat{\textbf{L}}\cdot\hat{\textbf{S}}$ and
interaction with the external magnetic field
$\mathcal{H}^{q}_{B}=\mu_{B}(g_{L}\hat{\textbf{L}}+g_{S}\hat{\textbf{S}})\textbf{B}$
with $g_{S}=2$, $g_{L}=1$ and $\mu_{B}$ representing the Bohr magneton. Here
$\hat{\textbf{S}}$ and $\hat{\textbf{L}}$ describe the spin and angular-
momentum operator, respectively. The symbols $\hat{\Theta}$ and $\hat{O}$
represent the equivalent Stevens operators along one of the cubic axis
$e_{\textrm{JT}}^{j}$ ($j=A,B,C$) and the trigonal axis
$\left[001\right]\parallel\mathbf{c}$-axis of GaN, respectively. These three
vectors $e_{\textrm{JT}}^{j}$ produce magnetic easy axes and are given by

$\begin{split}\scriptstyle\textbf{e}^{A}_{\textrm{JT}}=\left[\sqrt{\frac{2}{3}},0,\sqrt{\frac{1}{3}}\right],\textbf{e}^{B}_{\textrm{JT}}=\left[-\sqrt{\frac{1}{6}},-\sqrt{\frac{1}{2}},\sqrt{\frac{1}{3}}\right],\\\
\scriptstyle\textbf{e}^{C}_{\textrm{JT}}=\left[-\sqrt{\frac{1}{6}},\sqrt{\frac{1}{2}},\sqrt{\frac{1}{3}}\right]\end{split}$
(2)

In the case of single Mn3+ impurity in GaN with 3d4 configuration ($L=2$,
$S=2$), the matrix representation of the Hamiltonian is obtained using the
uncoupled spin and angular momentum basis states $|m_{L},m_{S}>$ with
$m_{S}=-2,-1,0,+1,+2$ and $m_{L}=-2,-1,0,+1,+2$. Thus the full matrix has
dimensions of ($25$ x $25$). By diagonalization of this matrix, one can obtain
the energy eigen-values $E$ and eigen-functions $|\phi>$. Then average
magnetic moment of Mn3+ ion ($\textbf{M}^{q}$) can be obtained from the
following formula,

$<\textbf{M}^{q}>=\frac{1}{Z}(Z_{A}<\textbf{m}^{A}>+Z_{B}<\textbf{m}^{B}>+Z_{C}<\textbf{m}^{C}>)$
(3)

where $Z_{j}=\sum_{i}exp(-E^{j}_{i}/k_{B}T)$ is the partition function at the
$j$-th JT center, $Z=Z_{A}+Z_{B}+Z_{C}$, and

$<\textbf{m}^{j}>=\frac{-\mu_{B}\sum_{i=1}^{25}<\phi_{i}^{j}|g_{L}\hat{\textbf{L}}+g_{S}\hat{\textbf{S}}|\phi_{i}^{j}>exp(-E^{j}_{i}/k_{B}T)}{\sum_{i=1}^{25}exp(-E^{j}_{i}/k_{B}T)}$
(4)

The symbols $E^{j}_{i}$ and $\phi_{i}^{j}$ correspond to $i$-th energy level
and eigen-state of Mn3+ ion being in the $j$-th JT center.

In this work, we extend the CFM model of single substitutional Mn3+ in GaN by
considering small clusters of Mn ions coupled by the ferromagnetic super-
exchange interaction
$\mathcal{H}^{q}_{Exchange}=-\sum_{\left[i,k\right]}J^{q}\hat{\textbf{S}}_{i}\cdot\hat{\textbf{S}}_{k}$,
where $\left[i,k\right]$ denotes all the nearest neighbour pairs within a
given cluster and $J^{q}$ is the value of exchange coupling constant. Since
these calculations are computationally demanding, the simulations are
restricted to clusters composed of up to four ions. The selected cluster types
are shown in Fig. 1. Only the ferromagnetic nearest neighbour (nn) interaction
is taken into account. Then, the relevant eigen-energy and eigen-functions of
a pair, triplet and quartet are found by the numerical diagonalization of full
($25^{2}$x$25^{2}$), ($25^{3}$x$25^{3}$) and ($25^{4}$x$25^{4}$) Hamiltonian
matrix respectively (c.f. Ref. 22 for details). Crystal field parameters are
listed in Tab. 1 and most of their numerical values are taken from Ref. 19.

$B_{4}$ | $B^{0}_{2}$ | $B^{0}_{4}$ | $\widetilde{B}^{0}_{2}$ | $\widetilde{B}^{0}_{4}$ | $\lambda_{TT}$ | $\lambda_{TE}$
---|---|---|---|---|---|---
11.44 | 3.19 | -0.423 | -5.85 | -1.17 | 5.5 | 11.5
Table 1: Parameters of the crystal field model of Mn3+ ion in GaN. All values
are in meV.

###  II.2 Classical model

The goal of this work is to compare the results of the classical and quantum
mechanical approach. Therefore the classical model should mimic the quantum
one as closely as possible. Here we use a classical atomic spin model, which
treats each atom as possessing a localised magnetic moment. The effective spin
Hamiltonian contains contributions from the exchange interactions,
magnetocrystalline anisotropies (Jahn-Teller and trigonal deformation) and
interaction with the external magnetic field. The same terms are present in
the quantum approach. Now, we neglect the spin-orbit interaction and the cubic
crystal field splitting (its contributions to the magnetic anisotropy is
negligible).

Therefore, the complete Hamiltonian in the classical approach is given by

$\mathcal{H}^{c}=\mathcal{H}^{c}_{B}+\mathcal{H}^{c}_{TR}+\mathcal{H}^{c}_{\textrm{JT}}+\mathcal{H}^{c}_{Exchange}$
(5)

where, $\mathcal{H}^{c}_{B}=-\mu_{S}\sum_{i}\textbf{m}_{i}\textbf{H}$ is the
Zeeman energy describing the interaction between spins and an external
magnetic field H. Now $\mathbb{m}$ is local spin moment direction, treated as
a classical vector with norm |$\mathbb{m}$|= 1, and $\mu_{S}=g\mu_{B}S$ is the
actual atomic moment, with $S=2$. The second and third term represent the
magnetic anisotropy energy $\mathcal{H}^{c}_{Anizo}$. Similarly as in the
quantum model, $\mathcal{H}^{c}_{Anizo}$ is expressed as a sum of uniaxial
$\mathcal{H}_{TR}$ and triaxial $\mathcal{H}_{\textrm{JT}}$ magnetocrystalline
anisotropy energy. These terms, in the classical approach, are defined as

$\mathcal{H}^{c}_{TR}=-\frac{1}{2}K_{TR}\sum_{i}\frac{1}{2}[m^{2}_{iz}-(m^{2}_{ix}+m^{2}_{iy}]$
(6)

$\mathcal{H}^{c}_{\textrm{JT}}=-\frac{1}{2}K_{\textrm{JT}}\sum_{i}\sum_{j=A,B,C}(\textbf{m}_{i}.\textbf{e}^{j}_{\textrm{JT}})^{4}$
(7)

$K_{TR}$, $K_{\textrm{JT}}$ are the anisotropy constants for the trigonal and
Jahn-Teller distortion, respectively. The last term represents a ferromagnetic
super-exchange interaction between spin moments

$\mathcal{H}^{c}_{Exchange}=-\sum_{\left[i,k\right]}J^{c}\textbf{m}_{i}\cdot\textbf{m}_{k}$
(8)

The exact value of $J^{c}$ (and $J^{q}$) is not known. Ab initio simulations
predict that the magnitudes of $J$ are rather high, between 15 [27] and 62 meV
[28]. In our recent work, we approximate the lower limit for the value of nn
exchange coupling between Mn ions in GaN, which result in $J\geq 8.6$ meV. All
values of $J$ are given assuming that the norm of Mn3+ spins $\mathbb{S}$ are
$|\mathbb{S}|=2$. Here, for comparison purposes, we choose a few different
values of the superexchange coupling parameter $J^{q}$, that is, 0.25, 0.5, 2
and 10 meV.

$J^{q}$ | 0.25 | 0.5 | 2 | 10
---|---|---|---|---
$J^{c}$ | 1.5 | 3 | 12 | 60
Table 2: Exchange coupling constants used for quantum ($J^{q}$) and classical
($J^{c}$) simulations. All values are in meV.

One should remember that in classical simulations, according to the standard
convention, we deal with magnetic spin moment $\mathbb{m}$ treated as a
classical vector with norm $|\mathbb{m}|=1$. On the other hand, the norm of
the quantized spin operator is $|\mathbb{\hat{S}}|=\sqrt{S(S+1)}$ what gives
$|\mathbb{\hat{S}}|=\sqrt{6}$ for S=2. Therefore the used values of $J^{c}$
should be six times larger than those applied in the crystal field model
approach ($J^{q}$), that is $J^{c}=S(S+1)J^{q}$. This ensures the same range
of exchange energy in both classical
($-J^{c}|\mathbb{m}_{1}|\cdot|\mathbb{m}_{2}$|) and quantum
($-J^{q}|\mathbb{\hat{S}}_{1}|\cdot|\mathbb{\hat{S}}_{2}$|) approximations.
The values of exchange constants are tabulated in Tab. 2.

Figure 2: (Color online) Motion of a magnetic moment (m) in an effective
magnetic field ($\textbf{H}_{eff}$).
$\textbf{m}\times(\textbf{m}\times\textbf{H}_{eff})$ and
$\textbf{m}\times\textbf{H}_{eff}$ represents damping and precessional terms,
respectively.

The dynamics of the magnetic spin moment treated as a classical vector in an
effective magnetic field $\textbf{H}_{eff}$ is governed by the Landau-
Lifshitz-Gilbert equation (Fig. 2) given by

$\frac{\partial\textbf{m}_{i}}{\partial
t}=-\frac{\gamma}{1+\alpha^{2}}[\textbf{m}_{i}\times\textbf{H}^{i}_{eff}+\alpha\textbf{m}_{i}\times(\textbf{m}_{i}\times\textbf{H}^{i}_{eff}]$
(9)

where $\alpha=1$ is the Gilbert damping constant, $\gamma$ is the gyromagnetic
ratio and the $\textbf{H}^{i}_{eff}$ is the net effective field acting on the
$i-th$ magnetic moment. This effective magnetic field is obtained by
calculating the spin Hamiltonian’s first derivative.

$\textbf{H}^{i}_{eff}=-\frac{1}{\mu_{s}}\frac{\partial\mathcal{H}}{\partial\textbf{m}_{i}}$
(10)

As we know, the LLG equation is strictly applicable for zero temperature
simulations. Thermal effects cause fluctuations in spin moments, and at
elevated temperatures, they can even overwhelm the exchange interactions
leading to a ferromagnetic-paramagnetic transition. To investigate the
temperature-dependent properties, we use Langevin dynamics by adding a
temperature-dependent term to $\textbf{H}_{eff}$. Then, the effective field is
written as

$\textbf{H}^{i}_{eff}=-\frac{1}{\mu_{s}}\frac{\partial\mathcal{H}}{\partial\textbf{m}_{i}}+\textbf{H}^{i}_{Th}$
(11)

The fundamental idea of Langevin dynamics is to represent the thermal
fluctuations by a Gaussian white noise term. Width of this Gaussian term
increases as the temperature increases, which means a more substantial thermal
effect on the system [6]. The thermal fluctuations are represented by a
Gaussian distribution $\bm{\Gamma}(t)$ in 3D space with a mean of zero and a
standard deviation of 1. At each time, an instantaneous thermal field acting
on each magnetic spin moment is given by

$\textbf{H}_{Th}=\bm{\Gamma(t)}\sqrt{\frac{2\alpha k_{B}T}{\gamma\mu_{S}\delta
t}}$ (12)

where $k_{B}$ is the Boltzman constant, $T$ is the simulation temperature, and
$\delta t=0.5\cdot 10^{-5}$ s is the integration time step. Further, we have
used $15\cdot 10^{8}$ equilibration and averaging steps for the simulation of
our system. The Langevin dynamics with the LLG equation is a stochastic
differential equation. One should note that the same results can be obtained
using the classical Monte-Carlo method [6].

##  III Numerical RESULTS

###  III.1 Adjustment of numerical parameters

Figure 3: (Color online) The single ion magnetic anisotropy energy
$E=-\mathbb{M}\cdot\mu_{0}\mathbb{H}$ at $T=2$ K and $\mu_{0}H=3$ T, where
$\mathbb{M}=\mathbb{M}(\mathbb{H},T)$. $\mathbb{H}$ can take all possible
values on the sphere surface with radius $\mu_{0}H$. The unit vector
describing the direction of given external magnetic field $\mathbb{H}$ is
represented using spherical coordinates $\theta,\phi$, namely
$\textbf{e}_{H}=[\mathrm{sin}\theta\mathrm{cos}\phi,\mathrm{sin}\theta\mathrm{sin}\phi,\mathrm{cos}\theta]$.
The magnetization $\mathbb{M}$ was calculated using (a) the quantum
$\mathbb{M}^{q}$ and (b) the classical $\mathbb{M}^{c}$ approach. The
difference in magnetic anisotropy energy for the quantum($\Delta E^{q}$) and
classical ($\Delta E^{c}$) simulations are 0.07 and 0.05 meV, respectively. No
contribution from the trigonal deformation is taken into account.

In order to perform the comparative study, we should use numerical parameters
that ensure the same single ion magnetic anisotropy in both classical and
quantum approximation, at least at low temperatures. Crystal field model
parameters are listed in Tab. 1. In order to adjust classical anisotropy
constants $K_{TR}$ and $K_{\textrm{JT}}$ we plot in Fig. 3 a single ion
magnetic anisotropy energy $-\mathbb{M}\cdot\mu_{0}\mathbb{H}$ at $T=2$ K and
$\mu_{0}H=3$ T, where $\mathbb{M}=\mathbb{M}(\mathbb{H},T)$. Only the Jahn-
Teller effect is included in these calculations with no contribution from
trigonal deformation. The applied external magnetic field vector $\mathbb{H}$
is represented using spherical coordinates. Such approach allows us to compare
quantum $-\mathbb{M}^{q}\cdot\mu_{0}\mathbb{H}$ (Fig. 3.a) and classical
$-\mathbb{M}^{c}\cdot\mu_{0}\mathbb{H}$ (Fig. 3.b) results, where ${M}^{q}$
corresponds to the magnetization obtained from CFM simulations, whereas
${M}^{c}$ describes the classical magnetization computed using the sLLG
equation. In Fig. 3, we observe very similar qualitative results obtained from
both approaches. The easy magnetic axes are along three different Jahn-Teller
directions $\mathbb{e}_{JT}$. The satisfactory fit of classical results to
quantum ones is obtained by taking $K_{\textrm{JT}}=0.75$ meV.

Figure 4: (Color online) Magnetization per one ion as a function of the
magnetic field of different magnetic clusters using (a) classical and (b)
quantum approaches at $T=2$ K. Solid and dotted lines denote magnetic easy
($\mathbb{H}\perp\mathbb{c}$) and hard axis ($\mathbb{H}$ $||$ $\mathbb{c}$)
respectively. The magnetization per ion increases with the number of ions $N$
in a given cluster due to the presence of ferromagnetic superexchange coupling
between Mn3+ ions.

Next, in Fig. 4 we plot the magnetization per one ion of different magnetic
clusters coupled by the ferromagnetic superexchange $J^{q}=0.5$ meV as a
function of the magnetic field $\mathbb{H}$ at $T=2$ K. As we can see from
Fig. 4, the magnetization per ion increases with the number of ions $N$ in a
given cluster. This increase in magnetization with $N$ is expected since a
ferromagnetic superexchange coupling exists between the Mn3+ ions. However,
the most prominent feature is the presence of a significant magnetic
anisotropy, with the easy axis perpendicular to the $\mathbb{c}$ axis of GaN
($\mathbb{H}\perp\mathbb{c}$), and the hard one parallel to $\mathbb{c}$
($\mathbb{H}$ $||$ $\mathbb{c}$). This magnetic anisotropy is mostly
controlled by trigonal deformation [12, 17, 19, 22]. In experimental
conditions, such deformation can be modified by the external electric field
through the inverse piezoelectromagnetic effect or by stretching the material
epitaxially [19]. By setting $K_{TR}=-1.24$ meV in classical simulations, we
reproduce quantum curves both qualitatively and quantitatively, as shown in
Fig. 5. Having adjusted classical parameters $K_{TR}$ and $K_{\textrm{JT}}$ at
$T=2$ K, we can perform comparative studies of magnetization obtained from the
quantum and the classical approach as a function of $T$.

Figure 5: (Color online) Comparison of magnetization per one ion as a function
of the magnetic field of different magnetic clusters at $T=2$ K with field
applied along the (a) magnetic easy ($\mathbb{H}\perp\mathbb{c}$) and (b) hard
axis ($\mathbb{H}$ $||$ $\mathbb{c}$). Solid and dotted lines denote quantum
and classical simulations respectively . Figure 6: (Color online)
Magnetization versus temperature of different magnetic clusters in an external
magnetic field of $\mu_{0}H=1$ T applied along the magnetic easy
($\mathbb{H}\perp\mathbb{c}$) and hard axis ($\mathbb{H}$ $||$ $\mathbb{c}$).
Solid and dotted lines denote quantum and classical simulations respectively.
The super-exchange coupling $J^{q}$ is 0.25 meV (a and b), 0.5 meV (c and d) ,
2 meV (e and f) and 10 meV (g and h). The dashed magenta lines indicate the
slope of the Curie law. Arrows in panels (a)-(d) indicate the value of
$T_{exch}=J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|/k_{B}$, that
is the temperature at which the magnitude of exchange energy equals the
thermal energy. In panels (e)-(h) $T_{exch}>100$ K.

###  III.2 Temperature-dependent magnetization

Studies on the magnetization of various clusters as a function of temperature
(up to 100 K) are shown in Fig. 6. Since the calculations are computationally
demanding, we restrict the upper limit of simulation temperatures to 100 K.
Both in classical and quantum simulations, a constant magnetic field of
$\mu_{0}H=1$ T is implemented. We present results for $\mathbb{H}$ applied
both along the magnetic easy (Fig. 6.a, c, e, g) and the hard axis (Fig. 6.b,
d, f, h). The magnetization is calculated for four different values of
$J^{q}$, namely 0.25, 0.5, 2, 10 meV. The $M(T)$ curves are presented in a
double logarithmic scale. As shown in Fig. 6, at the lowest temperatures
$T\lesssim 1$ K, the magnetization saturates. This effect is more pronounced
for larger clusters. At elevated temperatures, $T\gtrsim 20$ K one would
expect that ${M}(T)$ follows the Curie law (represented by dashed magenta
lines in Fig. 6). However, the $1/T$ dependency is obeyed only by singles, and
larger clusters with the magnitude of exchange energy
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|$ greater than the
thermal energy $k_{B}T$ ($0.1\leqslant T\leqslant 100$ K), as shown in Fig. 6
e-h. On the contrary, when
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|<k_{B}T$, the
individual spins creating a given cluster are partially decoupled. Then the
decrease in magnetization ${M}(T)$ with increasing $T$ is steeper for such
clusters than for the magnetization of single ion and ${M}(T)$ expected from
Curie law. At the highest temperatures investigated here, the magnetization of
larger clusters with low $J^{q}$ approaches the single ion results. This
situation is exemplified in Fig. 6 a-d, with $J^{q}=0.25$ and $0.5$ meV. The
detailed analysis of this effect, with an additional comparison with
experimental data, is also reported in Ref. 22.

Interestingly, we observe a small valley in the quantum $M^{q}(T)$ curves
centered at $T\approx 3$ K when the magnetic field is applied along the hard
axis ( $\mathbb{H}$ $||$ $\mathbb{c}$), as shown in Fig. 6 b, d, f, h. The
magnitude of this dip increases with the cluster’s size, and its center
slightly shifts to higher temperatures with $N$. It was only possible to
partially reproduce this effect using classical simulations. Nevertheless,
this phenomenon awaits its experimental verification.

Figure 7: (Color online) Comparison of finite temperature magnetization
studies using classical (${M}^{c}$) and quantum (${M}^{q}$) simulations. An
external field of $\mu_{0}H=1$ T is applied along magnetic easy
($\mathbb{H}\perp\mathbb{c}$) and hard axis ($\mathbb{H}$ $||$ $\mathbb{c}$)
with $J^{q}$= 0.25 meV (a and b), 0.5 meV (c and d) , 2 meV (e and f) and 10
meV (g and h). The $M^{c}/M^{q}$ ratios were calculated using the data
presented in Fig. 6.
$T_{exch}=J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|/k_{B}$
describes the temperature at which the magnitude of exchange energy equals the
thermal energy.

###  III.3 Comparative analysis of temperature-dependent magnetization
obtained using classical and quantum simulations

In the previous section, we can observe a discrepancy between results of
classical and quantum simulations presented in $\mathbb{M}(T)$ graphs. In
order to shed light on this disparity, we calculate the ratio of classical
magnetization to the quantum one $M^{c}/M^{q}$ for different values of $J^{q}$
at $\mu_{0}H=1$ T, as shown in Fig. 7. Ideally, this ratio should be exactly
one in the whole temperature range. However, the presented data indicate that
$M^{c}/M^{q}$ depends on the temperature, the value of super-exchange
coupling, and the size of the cluster $N$. First, we focus on results obtained
for the single non-interacting ion (black lines in Fig. 7). For the magnetic
easy axis ($\mathbb{H}\perp\mathbb{c}$), the ratio $M^{c}/M^{q}$ approaches
unity at very low temperatures $T\lesssim 2$ K (please note that we adjust the
classical parameters at $T=2$ K, as shown in Fig. 5) . For $\mathbb{H}$ $||$
$\mathbb{c}$, this ratio is slightly higher that one $M^{c}/M^{q}\simeq 1.2$
at $T\lesssim 2$ K, due to the fact that it was not possible to perfectly
reproduce the shape of $M^{q}(\mathbb{H})$ using classical simulations, as
seen in Fig. 5. As the temperature increases, $M^{c}/M^{q}$ ratio decreases
rather fast both for the easy and hard magnetization, and saturates at the
value of $M^{c}/M^{q}\simeq 0.75$ for $T\gtrsim 10$ K. Surprisingly, at very
high temperatures, where quantum and classical results should merge, the
$M^{c}/M^{q}$ ratio is low and significantly below unity. The reason for this
is explained in section III.4.

Now we turn to the data obtained for larger clusters and small values of
exchange coupling $J^{q}=0.25$ and $0.5$ meV. We find here that the
$M^{c}/M^{q}$ ratio clearly increases with the cluster size at low
temperatures, $T\leqslant 30$ K, as shown in Fig. 7 a-d. However $M^{c}/M^{q}$
still decreases with $T$. At high temperatures, where
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|<k_{B}T$, results for
larger clusters with $N\geq 2$ approach the single ion value, as expected. We
observe an improvement in the comparison of the results of classical and
quantum simulations with the increase of the exchange coupling $J^{q}$ to 2
meV. Finally, as shown in Fig. 7 g-h, $M^{c}/M^{q}$ is close to unity for the
tetrahedron quartet with $J^{q}=10$ meV in almost whole studied temperature
range.

To sum up, we see that the $M^{c}/M^{q}$ ratio depends strongly on the value
of super-exchange coupling $J$ and the size of the cluster $N$. The comparison
between results of classical and quantum simulations improves with increasing
magnitude of $J^{q}$. It seems, that for
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|<k_{B}T$ the
$M^{c}/M^{q}$ ratio is rather low. However, when
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|\gg k_{B}T$, the
$M^{c}/M^{q}$ ratio for $N\geq 2$ approaches unity, and it is close to one for
$N=4$ in almost whole studied temperature range.

###  III.4 Analysis

Figure 8: (Color online) Available number of states contained in the energy
interval ($E$, $E+\Delta E$), with $\Delta E=0.1$ meV for (a) singlet (b) pair
(c) closed triplet (d) tetrahedron quartet with $J^{q}$ = 0.5 meV and
$\mu_{0}H=1$ T. The magnetic field is applied along the hard axis $\mathbb{H}$
$||$ $\mathbb{c}$. The energy is presented relative to the energy of the
ground state $E_{GS}$.

In Fig. 7, we observed a significant discrepancy between the results of
quantum and classical simulations as a function of the temperature. As stated
previously, the above disparity stems from the classical nature of the sLLG
model [29, 7]. Generally, one can compute $M(T)$ using the classical approach
and reproduce rather well the Curie temperature $T_{\mathrm{C}}$, but not the
shape of experimental curves [7]. Even in the paramagnetic case studied here,
the classical treatment deviates significantly from results obtained using the
quantum treatment. Several methods have been proposed in the past to resolve
this issue [29, 7, 30, 31]. For example, Kuz’min introduced a simple
phenomenological equation to quantify the shape of spontaneous $M(T)$ [29].
Other approaches proceed by implementing temperature rescaling [7] or by
incorporating quantum Bose-Einstein statistics instead of classical Boltzmann
statistics [30]. The rescaling approach allows to modify the classical
thermodynamic quantities toward the quantum mechanical solution and include in
the simulations the underlying quantum effects in $ad$ $hoc$, and approximate
way [32]. However, these methods are primarily applicable to the bulk material
with a well defined critical temperature $T_{C}$.

To have a better insight into why the $M^{c}/M^{q}$ ratio is small for low
values of $J^{q}$, we calculated a quasi density of states (qDOS), which
describes the number of quantum states available within an energy range from
$E$ to $E+\Delta E$. The input data come from CFM simulations, and the results
are shown in Fig. 8. A simple explanation of the disparity between classical
and quantum simulation results refers to the concept of enhanced thermal spin
fluctuations in the classical limit owing to the effect of quantum mechanical
“stiffness” [7]. The physical interpretation is that in quantum mechanics,
only specific discrete energy levels exist, as shown in Fig. 8 a, for the case
of a single non-interacting Mn3+ ion. When the temperature $T$ is sufficiently
low, the thermal energy $k_{B}T$ becomes comparable with the energy difference
between the discrete quantum levels. The system populates mostly the ground
state with a low probability of occupying higher excited states. This
dramatically reduces the thermal fluctuations of the quantum system. On the
contrary, in the classical case, a continuous energy spectrum is observed, and
the allowed spin fluctuations are enhanced compared to the quantum case.
Hence, the classical approach fails in describing the magnetization at low
temperatures where $M^{c}<M^{q}$. On increasing temperature, the probability
of occupying higher excited states increases and the quantum mechanics
gradually merges with the classical one, bringing $M^{c}$ closer to $M^{q}$.
Therefore, one might expect that quantum simulations will give enhanced
results than the classical ones at lower temperatures, and the $M^{c}/M^{q}$
ratio should increase with $T$ when the quantum system approaches the
classical limit. However, a significant departure from this expected
characteristic is presented in Fig. 7. For low values of $J^{q}$, the
$M^{c}/M^{q}$ ratio acquires maximum (not minimum) values at lower
temperatures, and it decreases with $T$ (instead of increasing towards unity).
Finally, $M^{c}/M^{q}$ saturates at low level for $T\gtrsim 50$ K.

The figure 8 also illustrates a substantial variation of qDOS for different
cluster types. It is of a discrete type for the singlet, but acquires a quasi-
continuous form for the pair, closed triplet and tetrahedron quartet. This
gradual transition from the discrete to the quasi-continuous form of qDOS as
the cluster size increases is consistent with the data presented in Fig. 7
c,d. The correspondence between classical and quantum results is the best at
$T<20$ K for tetrahedron quartet (as compared with other clusters with $N<4$),
that is for the largest cluster whose quantum level structure, depicted in
Fig. 8 d, most closely resembles the classical case. The conclusion is that,
it is possible to only partially explain the inconsistency between classical
and quantum results presented in temperature-dependent magnetization curves by
referring to the concept of quantum mechanical “stiffness”.

Instead, one can invoke simple analytic expressions for magnetic
susceptibility in order to explain the observed differences between $M^{c}$
and $M^{q}$. It is well known that in the isotropic case, the magnetization of
a single ion can be calculated analytically, using both quantum mechanical and
classical considerations, by referring to the Brillouin and Langevin
functions, respectively. When $\mu_{B}\mu_{0}H/k_{B}T$ is small (low to
moderate $\mu_{0}H$, high $T$), it is possible to derive a quantum
$\chi^{q}=\frac{\mu_{0}g^{2}\mu_{B}^{2}S(S+1)}{3k_{B}T}$ and a classical
$\chi^{c}=\frac{\mu_{0}g^{2}\mu_{B}^{2}S^{2}}{3k_{B}T}$ approximations of
magnetic susceptibilities (the Curie law). It is of interest to compare the
two expressions. Clearly, the ratio $\frac{\chi^{c}}{\chi^{q}}=\frac{S}{S+1}$
depends on the value of spin $S$, and both approaches merge for
$S\rightarrow\infty$. In our opinion, the presented above relation explains
rather well the temperature-dependent data depicted in Fig. 7. For single non
interacting ions, $M^{c}/M^{q}$ saturates for high temperatures at
$M^{c}/M^{q}\simeq 0.75$. This saturation value agrees reasonably well with
the analytical level $\chi^{c}/\chi^{q}=2/3$ obtained for $S=2$. The observed
small discrepancy between these two magnitudes originates from the presence of
the magnetic anisotropy in our numerical calculations, whereas
$\chi^{c}/\chi^{q}$ was obtained for an isotropic case. Similarly, for larger
clusters and a small value of exchange coupling $J^{q}=0.25$ and $0.5$ meV
(Fig. 7 a-d), the individual spins are practically decoupled at $T\gtrsim 50$
K, and $M^{c}/M^{q}$ for $N\geqslant 2$ approaches the single ion value.
Interestingly, the increase of temperature does not lead to the merging of
classical and quantum results as $M^{c}/M^{q}$ have to stay close to the
fundamental single ion limit defined by $\chi^{c}/\chi^{q}=2/3$ with $S=2$.
However, in the case when the exchange energy is very strong $J^{q}=10$ meV,
that is $J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|\gg k_{B}T$, we
observe quite different characteristics (Fig. 7 g-h). At high temperatures,
the $M^{c}/M^{q}$ ratio increases with $N$, and approaches unity for $N=4$.
For $J^{q}=10$ meV, the Mn3+ ions are strongly bound in the whole temperature
range investigated here, $0.25\geqslant T\geqslant 100$ K, behaving as an
effective single ion with the total spin of $S_{T}=NS$. Then $M^{c}/M^{q}$
ratio at high $T$ follows roughly a simple relation
$\chi^{c}/\chi^{q}=NS/(NS+1)$, which is equal to 4/5, 6/7 and 8/9 for the
pair, triplet and quartet respectively. It is consistent with the simple fact,
that the classical and quantum simulations merge for $S\rightarrow\infty$.
However, in this paper, depending on the relative strength between the thermal
and the exchange energy, we deal with uncoupled ions characterized by $S=2$ or
strongly bound clusters with large total spin value $S_{T}=NS$.

The analysis presented above explains why using the classical approach, it is
very hard to properly model the shape of the experimental magnetization curves
of different ferromagnetic materials as a function of temperature. The problem
is the incorrect form of the classical $M(T)$ curve when compared to
experiment [7, 30, 33]. At low temperatures, well below $T_{C}$, the magnetic
ions are strongly bound creating percolating magnetic system with a high total
spin $S_{T}$. In this regime, the classical approach can give correct and
satisfactory results. However, by increasing $T$, but still below $T_{C}$,
spins become weakly coupled (a large effective cluster gradually breaks down
into small pieces with low $S$), and then the classical approximation
underestimates the actual value of the magnetization.

These conclusions are substantiated by the results shown in Fig. 7 e-f. For
$J^{q}=2$ meV the magnitude of the exchange energy is still slightly greater
than the thermal energy in the entire studied temperature range, that is
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|>k_{B}T$. However
$M^{c}/M^{q}$ for larger clusters decreases with increasing $T$ and gradually
approaches the single ion limit. Only in the limit of very strong coupling
with $J^{q}=10$ meV, where
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|>6k_{B}T$, the given
cluster acts as an effective single-ion spin with $S_{T}=NS$, as shown in Fig.
7 g-h.

##  IV CONCLUSIONS

We have performed a comparative study of magnetization
$\mathbb{M}(\mathbb{H},T)$ of Mn3+ ions ($S=2$, $L=2$) in GaN using classical
and quantum mechanical approaches. We have investigated small magnetic
clusters consisting of isolated ions, pairs, triplets, and quartets of Mn3+
ions coupled by the ferromagnetic super-exchange interaction. We have found
significant differences between the results of both approaches (classical and
quantum crystal field model). Especially, we found that the ratio of classical
magnetization to the quantum one $M^{c}/M^{q}$ strongly varies with
temperature and it depends on the choice of the exchange coupling constant
$J^{q}$ and the size of the cluster $N$. The $M^{c}/M^{q}$ ratio improves with
increasing of $J^{q}$ and approaches unity in the strong coupling limit
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|\gg k_{B}T$ and for
large value of $N$. Our results show that it is very complex task, to
reproduce the complicated quantum crystal field model CFM using a simple
classical approximation. We have found that the key factor is the magnitude of
the effective spin of the system, which depends on the relative strength
between the thermal and the exchange energy. For weak coupling between the
spins $J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|<k_{B}T$, the
investigated clusters consist of practically uncoupled spins characterized by
$S=2$. Then classical approximation underestimates the actual value of
magnetization. Contrary, in the strong coupling regime, characterized by
$J^{q}|\mathbb{\hat{S}}_{i}|\cdot|\mathbb{\hat{S}}_{j}|\gg k_{B}T$, we deal
with one magnetic system characterized by high total spin value $S_{T}=NS$. In
this limit, $S_{T}$ increases with the number of ions $N$ in the cluster, what
improves the comparison between $M^{c}$ and $M^{q}$. These conclusions are
consistent with the simple fact, that the classical and quantum simulations
differ for small value of $S$ and both treatments merge for
$S\rightarrow\infty$.

## Acknowledgments

We would like to thank M. Sawicki for proofreading of the manuscript and
valuable suggestions. The work is supported by the National Science Centre
(Poland) through project OPUS 2018/31/B/ST3/03438. This research was carried
out with the support of the Interdisciplinary Centre for Mathematical and
Computational Modelling (ICM) University of Warsaw under grant no GB77-6

## References

  * [1] K. L. Brown, W. J. Munro, and V. M. Kendon.  Using quantum computers for quantum simulation.  Entropy, 12(11):2268–2307, 2010.
  * [2] Y. Zhou, E. M. Stoudenmire, and X. Waintal.  What limits the simulation of quantum computers?  Phys. Rev. X, 10:041038, Nov 2020.
  * [3] D. A Meyer.  Quantum computing classical physics.  Phil. Trans. R. Soc. 360(1792):395–405, 2002.
  * Foulkes et al. [2001] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal.  Quantum monte carlo simulations of solids.  _Rev. Mod. Phys._ , 73(1):33, 2001.  doi: 10.1103/revmodphys.73.33.
  * Wilson [1971] K. G. Wilson.  Renormalization group and critical phenomena. I. Renormalization group and the Kadanoff scaling picture.  _Phys. Rev. B_ , 4(9):3174, 1971.  doi: 10.1103/physrevb.4.3174.
  * Evans et al. [2014] R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis, and R. W. Chantrell.  Atomistic spin model simulations of magnetic nanomaterials.  _Journal of Physics: Condensed Matter_ , 26(10):103202, 2014.  doi: 10.1088/0953-8984/26/10/103202.
  * Evans et al. [2015] R. F. L. Evans, U. Atxitia, and R. W. Chantrell.  Quantitative simulation of temperature-dependent magnetization dynamics and equilibrium properties of elemental ferromagnets.  _Phys. Rev. B_ , 91(14):144425, 2015.  doi: 10.1103/physrevb.91.144425.
  * Schlickeiser et al. [2014] F. Schlickeiser, U. Ritzmann, D. Hinzke, and U. Nowak.  Role of entropy in domain wall motion in thermal gradients.  _Phys. Rev. Lett._ , 113:097201, 2014.  doi: 10.1103/PhysRevLett.113.097201.
  * Ellis and Chantrell [2015] M. O. A. Ellis and R. W. Chantrell.  Switching times of nanoscale FePt: finite size effects on the linear reversal mechanism.  _Appl. Phys. Lett._ , 106(16):162407, 2015.  doi: 10.1063/1.4919051.
  * Bender et al. [2017] S. A. Bender, H. Skarsvåg, A. Brataas, and R. A. Duine.  Enhanced spin conductance of a thin-film insulating antiferromagnet.  _Phys. Rev. Lett._ , 119:056804, 2017.  doi: 10.1103/PhysRevLett.119.056804.
  * Hinzke and Nowak [2011] D. Hinzke and U. Nowak.  Domain wall motion by the magnonic spin seebeck effect.  _Phys. Rev. Lett._ , 107:027205, Jul 2011.  doi: 10.1103/PhysRevLett.107.027205.
  * Gosk et al. [2005] J. Gosk, M. Zajac, A. Wolos, M. Kaminska, A. Twardowski, I. Grzegory, M. Bockowski, and S. Porowski.  Magnetic anisotropy of bulk GaN: Mn single crystals codoped with Mg acceptors.  _Phys. Rev. B_ , 71(9):094432, 2005.  doi: 10.1103/physrevb.71.094432.
  * [13] S. Stefanowicz, G. Kunert, C. Simserides, J. A. Majewski, W. Stefanowicz, C. Kruse, S. Figge, Tian Li, R. Jakieła, K. N. Trohidou, A. Bonanni, D. Hommel, M. Sawicki, and T. Dietl.  Phase diagram and critical behavior of a random ferromagnet Ga1-xMnxN.  Phys. Rev. B, 88:081201(R), 2013.
  * [14] M. Grodzicki, P. Mazur, and A. Sabik.  Electronic properties of p-GaN co-doped with Mn by thermal process: Surface studies.  Surface Science, 689:121460, 2019.
  * [15] E. Sarigiannidou, F. Wilhelm, E. Monroy, R. M. Galera, E. Bellet-Amalric, A. Rogalev, J. Goulon, J. Cibert, and H. Mariette.  Intrinsic ferromagnetism in wurtzite (Ga,Mn)N semiconductor.  Phys. Rev. B, 74:041306, 2006.
  * Wolos et al. [2004] A. Wolos, A. Wysmolek, M. Kaminska, A. Twardowski, M. Bockowski, I. Grzegory, S. Porowski, and M. Potemski.  Neutral Mn acceptor in bulk GaN in high magnetic fields.  _Phys. Rev. B_ , 70(24):245202, 2004.  doi: 10.1063/1.1994087.
  * Stefanowicz et al. [2010] W. Stefanowicz, D. Sztenkiel, B. Faina, A. Grois, M. Rovezzi, T. Devillers, F. d’Acapito, A. Navarro-Quezada, T. Li, R. ł Jakieła, et al.  Structural and paramagnetic properties of dilute Ga1-xMnxN.  _Phys. Rev. B_ , 81(23):235210, 2010.  doi: 10.1103/physrevb.81.235210.
  * Bonanni et al. [2011] A. Bonanni, M. Sawicki, T. Devillers, W. Stefanowicz, B. Faina, Tian Li, T. E. Winkler, D. Sztenkiel, A. Navarro-Quezada, M. Rovezzi, R. Jakieła, A. Grois, M. Wegscheider, W. Jantsch, J. Suffczyński, F. d’Acapito, A. Meingast, G. Kothleitner, and T. Dietl.  Experimental probing of exchange interactions between localized spins in the dilute magnetic insulator (Ga,Mn)N.  _Phys. Rev. B_ , 84:035206, 2011.  doi: 10.1103/PhysRevB.84.035206.
  * Sztenkiel et al. [2016] D. Sztenkiel, M. Foltyn, G. P. Mazur, R. Adhikari, K. Kosiel, K. Gas, M. Zgirski, R. Kruszka, R. Jakiela, T. Li, et al.  Stretching magnetism with an electric field in a nitride semiconductor.  _Nat.Commun_ , 7(1):1–9, 2016.  doi: 10.1038/ncomms13232.
  * Gas et al. [2018] K. Gas, J. Z. Domagala, R. Jakiela, G. Kunert, P. Dluzewski, E Piskorska-Hommel, W. Paszkowicz, D. Sztenkiel, M. J. Winiarski, D. Kowalska, R. Szukiewicz, T. Baraniecki, A. Miszczuk, D. Hommel, and M. Sawicki.  Impact of substrate temperature on magnetic properties of plasma-assisted molecular beam epitaxy grown (Ga,Mn)N.  _J. of Alloys and Compounds_ , 747:946, 2018.  doi: https://doi.org/10.1016/j.jallcom.2018.03.056.
  * Gas et al. [2021] K. Gas, G. Kunert, P. Dluzewski, R. Jakiela, D. Hommel, and M. Sawicki.  Improved-sensitivity integral squid magnetometry of (Ga,Mn)N thin films in proximity to Mg-doped GaN.  _Journal of Alloys and Compounds_ , 868:159119, 2021.  ISSN 0925-8388.  doi: https://doi.org/10.1016/j.jallcom.2021.159119.
  * Sztenkiel et al. [2020] D. Sztenkiel, K. Gas, J.Z. Domagala, D. Hommel, and M. Sawicki.  Crystal field model simulations of magnetic response of pairs, triplets and quartets of Mn3+ ions in gan.  _New Journal of Physics_ , 22(12):123016, 2020\.  doi: 10.1088/1367-2630/abcce7.
  * [23] Y. Shapira and V. Bindilatti.  Magnetization-step studies of antiferromagnetic clusters and single ions: Exchange, anisotropy, and statistics.  J. Appl. Phys, 92(8):4155–4185, 2002.
  * Vallin and Watkins [1974] J. T. Vallin and G. D. Watkins.  Epr of Cr2+ in II-VI lattices.  _Phys. Rev. B_ , 9(5):2051, 1974.  doi: 10.1103/physrevb.9.2051.
  * Vallin et al. [1970] J. T. Vallin, G. A. Slack, S. Roberts, and A. E. Hughes.  Infrared absorption in some II-VI compounds doped with Cr.  _Phys. Rev. B_ , 2(11):4313, 1970.  doi: 10.1103/physrevb.2.4313.
  * Tracy et al. [2005] J. L. Tracy, G. Franzese, Ashlee Byrd, J. Garner, T. M. Pekarek, I. Miotkowski, and A. K. Ramdas.  Anisotropic magnetization of the III-VI diluted magnetic semiconductor In1-xMnxS in the mixed state.  _Phys. Rev. B_ , 72(16):165201, 2005.  doi: 10.1103/physrevb.72.165201.
  * Szwacki et al. [2011] N. G. Szwacki, J. A. Majewski, and T. Dietl.  Aggregation and magnetism of Cr, Mn, and Fe cations in GaN.  _Phys. Rev. B_ , 83(18):184417, 2011.  doi: 10.1103/PhysRevB.83.184417.
  * Sato et al. [2010] K. Sato, L. Bergqvist, J. Kudrnovskỳ, P. H. Dederichs, O. Eriksson, I. Turek, B. Sanyal, G. Bouzerar, H. Katayama-Yoshida, V. A. Dinh, et al.  First-principles theory of dilute magnetic semiconductors.  _Rev. Mod. Phys._ , 82(2):1633, 2010.  doi: 10.1103/revmodphys.82.1633.
  * Kuz’min [2005] M. D. Kuz’min.  Shape of temperature dependence of spontaneous magnetization of ferromagnets: Quantitative analysis.  _Phys. Rev. Lett._ , 94:107204, 2005.  doi: 10.1103/PhysRevLett.94.107204.
  * Bergqvist and Bergman [2018] L. Bergqvist and A. Bergman.  Realistic finite temperature simulations of magnetic systems using quantum statistics.  _Phys. Rev. Materials_ , 2:013802, 2018.  doi: 10.1103/PhysRevMaterials.2.013802.
  * Woo et al. [2015] C. H. Woo, Haohua Wen, A. A. Semenov, S. L. Dudarev, and Pui-Wai Ma.  Quantum heat bath for spin-lattice dynamics.  _Phys. Rev. B_ , 91:104306, 2015.  doi: 10.1103/PhysRevB.91.104306.
  * Körmann et al. [2010] F. Körmann, A. Dick, T. Hickel, and J. Neugebauer.  Rescaled Monte Carlo approach for magnetic systems: Ab initio thermodynamics of bcc iron.  _Phys. Rev. B_ , 81:134425, Apr 2010.  doi: 10.1103/PhysRevB.81.134425.
  * Körmann et al. [2011] F. Körmann, A. Dick, T. Hickel, and J. Neugebauer.  Role of spin quantization in determining the thermodynamic properties of magnetic transition metals.  _Phys. Rev. B_ , 83:165114, Apr 2011.  doi: 10.1103/PhysRevB.83.165114.