Quantifying time in Monte Carlo simulations:
Application to relaxation processes and AC
susceptibilities of magnetic nanoparticles assemblies

A. Morjane, J.-G. Malherbe, Juan J Alonso, F. Vernay, V. Russier

► To cite this version:

A. Morjane, J.-G. Malherbe, Juan J Alonso, F. Vernay, V. Russier. Quantifying time in Monte Carlo simulations: Application to relaxation processes and AC susceptibilities of magnetic nanoparticles assemblies. Journal of Magnetism and Magnetic Materials, 2025, 630, pp.173443. 10.1016/j.jmmm.2025.173443. hal-05232688

HAL Id: hal-05232688

https://hal.science/hal-05232688v1

Submitted on 1 Sep 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Quantifying time in Monte Carlo simulations: application to relaxation
processes and AC susceptibilities of magnetic nanoparticles assemblies.

A. Morjane,¹,² J.-G. Malherbe,²,\* Juan J. Alonso,³,⁴,† F. Vernay,¹ and V. Russier²,‡

¹PROMES-CNRS (UPR-8521), Université de Perpignan Via Domitia,
Rambla de la thermodynamique, Tecnosud, 66100 Perpignan, France.

²ICMPE, UMR 7182 CNRS and UPE 2-8 rue Henri Dunant 94320 Thiais, France.

³Física Aplicada I, Universidad de Málaga, 29071 Málaga, Spain

⁴Instituto Carlos I de Física Teórica y Computacional,
Universidad de Málaga, 29071 Málaga, Spain

The study of the response of magnetic nanoparticles (MNP) assemblies to an external alternating magnetic field is of great interest for applications such as hyperthermia. The key quantity here is the complex susceptibility and its behavior in terms of temperature and frequency. From a theoretical point of view it can be obtained by Monte Carlo (MC) simulation with the time quantified Monte Carlo (TQMC) method if a physical time is associated with the MC step. Here we revisit this method by showing that the time unit can be derived from the MC stochastic process of the isolated particle. We first obtain a MC unit of time from the relaxation of the system at fixed temperature. Then this unit of time is used to compute complex susceptibilities. We show that it is now possible to match the TQMC results with actual experimental results regarding frequency dependent in phase susceptibilities and quantify the unit of time in seconds. Finally we show that the time unit obtained for the isolated particle remains valid when considering interacting particles such as the Heisenberg coupling or dipole dipole interactions.

<footnote>\* e-mail address: malherbe@u-pec.fr</footnote>

<footnote>† e-mail address: jjalonso@uma.es</footnote>

<footnote>‡ e-mail address: vincent.russier@cnrs.fr</footnote>
# I. INTRODUCTION

Magnetic nanoparticles (MNP) assemblies still arouse a great interest on the experimental and theoretical points of view due to the quite large field of applications¹. In particular the dynamic properties, either stemming from the blocking or Néel process or resulting from the application of an oscillating external field, are of great importance such as the induced heat transfer to the embedding medium giving rise to the magnetic hyperthermia with applications in cancer therapy²⁻⁴ or catalysis⁵ for instance. In the framework of mesoscopic scale models such as the family of effective one spin models, where the MNP are represented as homogeneously polarized single domain objects, the simulation of dynamic properties of magnetic nanoparticles assemblies can be performed through different methods or approaches. Many of them start from the stochastic LLG equation namely the deterministic Gilbert equation to which the thermal fluctuations are included through a fluctuating field, leading to the Fokker-Planck equation for the spin orientation probability distribution. On an other hand, the Monte Carlo simulation method (MC), besides the two-states scheme⁶,⁷ and the kinetic MC⁸, has been adapted in order to quantify the Monte Carlo step to an actual time step measured in seconds. After the dynamical Monte Carlo scheme used in⁹ where dynamics are introduced via a temperature independent constraint on the move amplitude of the moments, the actual pioneering works in this direction are the so-called time quantified Monte Carlo (TQMC) method by Nowak et al.¹⁰, and extended by Cheng et al.¹¹. As a general rule, the different applications of the TQMC concern ensembles of MNP fixed in position and either non interacting¹⁰,¹²⁻¹⁴, namely isolated MNP, or including short range Heisenberg¹¹ or dipolar interactions (DDI)¹⁵⁻¹⁷. The starting point is the mapping of the time step from the stochastic LLG equation to the one of the Monte Carlo step through a temperature dependent constraint on the amplitude of the moments displacements in the MC sampling assuming the heat-bath sampling¹⁰⁻¹². Melenev et al.¹⁴ used another method for the mapping of the TQMC scheme by introducing the number of MC steps corresponding to the relaxation time of the Néel process as a MC time scale, and focused on the superparamagnetic relaxation of non interacting NP. Billoni et al.¹³ used the TQMC in order to study the behavior of switching time of an isolated MNP in the Stoner-Wolffhart model. Stariolo et al.¹⁵ used the TQMC to calculate both the relaxation time and the autocorrelation function of an ensemble of NP on a 2-D square lattice, including DDI. More recently, this method was used for instance in¹⁶,¹⁷ for modeling the influence of DDI on the hyperthermia properties of MNP assemblies. Stricto-sensu, the TQMC leads to a mapping of the MC step on an actual time scale but this is related to the gyromagnetic factor $\gamma$ and the damping factor $\alpha$. However, the precise knowledge of $\alpha$ for magnetic nanoparticles is missing (see for instance¹⁸) and this makes difficult the correspondence between the TQMC step and the actual time scale. This may be compared to the usual correspondence made in the FC/ZFC protocol between the relaxation time at the blocking temperature and the measuring time where the prefactor of the Arrhenius form of the Néel relaxation time is generally estimated with a large uncertainty (10⁻⁸ to 10⁻¹⁰ s). Hence a convenient determination of the time scale through the TQMC is still missing as was mentioned in¹⁹. For this reason, in Refs.¹⁶,²⁰ the true time scale is not made explicit but instead stated as directly proportional to the MC step (MCS) once the temperature dependent constraint mentioned above is applied but in the framework of the Metropolis sampling. A similar recipe for performing dynamical MC from Metropolis sampling was already introduced in²¹⁻²³ imposing a temperature independent acceptance rate fixed in between 0.3 and 0.4. Similarly, in²⁴ the usual Metropolis sampling is used with the constraint of a constant acceptance ratio to perform MC simulations of the susceptibility in terms of the frequency $\omega$ with assumed proportionality between 1/ω and MCS.

In this paper, we address this question by demonstrating how a MC unit of time can be obtained from the details of the simulation of the isolated particles. We show that this unit of time remains valid when dealing
with interacting particles. Moreover we show how to extend the TQMC scheme to the Metropolis sampling instead of the heat-bath one which was used to map the MC time step to the LLG one$^{10,11}$. Therefore we explicit a time unit for both the heath bath and Metropolis samplings. To this aim, we propose an alternative way to use the TQMC. We start from the mapping deduced in$^{10,11}$ without specifying the values of the gyromagnetic factor nor the damping factor $\alpha$, seen as characteristics of the microscopic features of the system not taken into account in the mesoscopic scale models. Instead we leave a free parameter for the correspondence between the MC step and the actual time step. First we demonstrate that the TQMC step indeed behaves as a convenient time scale for both non interacting and interacting assemblies of MNP represented by dipolar hard spheres plus uniaxial anisotropy (the so-called effective one spin model).

Then we show from the Néel relaxation process of non interacting MNP that we recover for the relaxation time in terms of temperature $T$ and anisotropy barrier height in $k_B T$ unit, $\sigma = K_u/k_B T$ either the Brown$^{25}$ or the Bessais et al.$^{26}$ form by fitting a multiplicative constant. This is interpreted as a first step to quantify the TQMC scheme. Then we propose a mapping leading to a correspondence of the MC time with the true time scale through the behavior of the real part of the susceptibility with respect to the frequency of real systems which are standard quantities considered in experiments.

This paper is organized as follows. After presenting the model in Section II, the TQMC is revisited in Section III by defining a unit of time from the details of the simulation in the sense that it is sampling choice independent. In Section IV this time scale is tested first for non interacting case. Then we show how it is possible to map the simulation results with experimental ones in order to obtain the true time scale. Finally we demonstrate that this unit of time remains valid for interacting systems. The paper ends with a brief conclusion.

## II. MODEL

We model an assembly of $N$ spherical MNP in the single domain regime, free of super exchange interactions and characterized by a uniaxial magnetocrystalline anisotropy energy (MAE) and monodisperse in size. We place ourselves in the framework of the effective one spin model where each single domain MNP is assumed to be uniformly magnetized with a temperature independent saturation magnetization $M_s$. Hence, we consider a system of dipolar hard spheres of moment $\vec{\mu}_i = \mu_i \hat{\mu}_i$, with $\mu_i = M_s v_i$, interacting through the usual dipole dipole interaction (DDI), subjected to a one-body anisotropy energy, $K_i v_i (\hat{\mu}_i \cdot \hat{\mu}_i)^2$ and to an external field $\vec{H} = H_0 \hat{h}$ where $H_0$ may be time dependent. $K_i$, $\hat{\mu}_i$ and $v_i$ are the anisotropy constant, the easy axis and the volume of the particle $i$ respectively. The MNP ensemble is monodisperse with MNP diameter $d$ and moment $\mu$. The Hamiltonian of the system is given by

$$ \beta H = \frac{1}{2}\beta\epsilon_d \sum_{i \neq j} \frac{\hat{\mu}_i \hat{\mu}_j - 3(\hat{\mu}_i \hat{r}_{ij})(\hat{\mu}_j \hat{r}_{ij})}{(r_{ij}/d)^3} - \beta K_u v(d) \sum_i (\hat{\mu}_i \cdot \hat{\mu}_i)^2 + \frac{1}{2}\beta \sum_{i \neq j} v_{sr}(r_{ij}) - \beta \mu_0 H_0 \mu \sum_i \hat{h} \cdot \hat{\mu}_i \quad (1) $$
with $\epsilon_d = \frac{\mu_0 \mu^2}{4\pi d^3}$

where $\hat{r}_{ij}$ is the unit vector carried by the vector joining sites $i$ and $j$, $\beta = 1/(k_B T)$ is the inverse temperature. $v_{sr}(r)$ is the short range contribution, taken as the hard sphere potential, $v_{sr}(r > d) = 0$ and $v_{sr}(r < d) = \infty$. We introduce a reference temperature, say $T_0$ ($\beta_0 = 1/k_B T_0$) and the reduced temperature is chosen as $T^* = T/T_0$. In the following we restrict ourselves to an assembly of fixed MNP with a fixed distribution of easy axes and accordingly we can drop the short ranged contribution, $v_{sr}$. Moreover the MNP are located
on the nodes of a perfect lattice of either FCC for of BCC symmetry. Equation (1) is then rewritten as

$$ \beta H = \frac{1}{T^*} \left( \frac{\lambda_d}{2} \sum_{i \neq j} \frac{\hat{\mu}_i \hat{\mu}_j - 3(\hat{\mu}_i \hat{r}_{ij})(\hat{\mu}_j \hat{r}_{ij})}{(r_{ij}/d)^3} - \lambda_u \sum_i (\hat{n}_i \cdot \hat{\mu}_i)^2 - h \sum_i \hat{h} \cdot \hat{\mu}_i \right) \quad (2) \\ \lambda_u = \beta_0 K_u v(d); \quad \lambda_d = \beta_0 \epsilon_d; \quad h = \beta_0 \mu_0 \mu H_0 \equiv 2\lambda_u \frac{H_0}{H_K} $$

which introduces the DDI, MAE and Zeeman coupling constants, $\lambda_d$, $\lambda_u$ and $h$, where $H_K = 2K/\mu_0 M_s$ is the usual anisotropy field.

We will also consider the case of an effective one spin model where the particles interact through the usual Heisenberg interaction of coupling constant $J$, in unit $1/\beta_0$, restricted to first neighbour pairs $<ij>$ instead of the DDI above. Then, the hamiltonian is given by

$$ \beta H = \frac{1}{T^*} \left( J \sum_{<ij>} \hat{\mu}_i \cdot \hat{\mu}_j - \lambda_u \sum_i (\hat{n}_i \cdot \hat{\mu}_i)^2 - h \sum_i \hat{h} \cdot \hat{\mu}_i \right) \quad (3) $$

Here, we consider only the ferromagnetic Heisenberg model, $J < 0$, and the ordered phase at low temperature is either of ferromagnetic (FM) or spin-glass (SG) nature according to the value of $\lambda_u$^27 when the easy axes are randomly distributed. The simulation box is a cube with edge length $L$, the total number of dipoles is $N$ and the volume fraction is $\phi_v = N\pi d^3/(6L^3)$. We consider periodic boundary conditions by repeating the simulation cubic box identically in the 3 dimensions. The long range DDI interaction is treated through the Ewald summation technique^28,29.

### III. METHOD

We start from the TQMC scheme as introduced first in Ref.¹⁰ and then improved in Ref.¹¹, based on the mapping of the MC time to the actual time $t$ of the LLG equation. We refer the reader to these works for details on how this mapping is deduced. We just need to recall that it is based on the precise tuning of the aperture of the cone, $R$, in which the new moment at site $i$ say $\hat{\mu}_i^{(n)}$ is chosen from its preceding direction, $\hat{\mu}_i^{(o)}$, resulting in the following relation between the time step $\delta t$ and $R$

$$ \delta t = a \frac{R^2}{T^*} = a R_0^2 \quad (4) $$

where $a$ is a constant and we have introduced $R_0$, independent of $T^*$. The constant $a$ is a function of both the gyromagnetic factor $\gamma$ and the damping coefficient $\alpha$ and can be calculated, at least in principle once the value of $\alpha$ is known. Conversely, we consider that it includes the microscopic parameters, not explicitly taken into account in the effective one spin model and our alternative approach is to keep $a$ as a fitting parameter, instead of determining its value from known values of $\gamma$ and $\alpha$. Notice that in Ref.¹¹, the relation (4) is deduced without starting from the equilibrium condition, and includes the value the acceptance rate $\Gamma$ of the MC process used with the heat-bath acceptance rate to the first order in $\Delta E$ as opposed to the initial mapping performed in¹⁰. Therefore, at least in principle, equation (4) holds only for the heat-bath acceptance rate and for sufficiently small displacements so that the first order expansion of $\Gamma$ holds.
Now, from equation (4) we deduce that the time duration corresponding to $N_{MCS}$ Monte Carlo steps characterized by a constrained aperture cone $R^2 = R_0^2 T^*$ is $T^*$ independent and is given by

$$ \Delta t = a N_{MCS} R_0^2 \quad (5) $$

from which we introduce a unit of time for the TQMC scheme corresponding to $\Delta N$ Monte Carlo steps

$$ u_t = \Delta N R_0^2 / 2. \quad (6) $$

This definition holds a priori only when the heat-bath acceptance rate is used and $R_0$ is small enough. The factor 1/2 is introduced for convenience (see below equation (7)). In order to generalize to the Metropolis acceptance rate, we first notice the disymmetry of the latter with respect to $\Delta E = 0$ and second that for small displacements, where the heat-bath acceptance rate leads to $\Gamma \sim 0.5^{-}$, the Metropolis one leads to $\Gamma$ being closer to $\Gamma \sim 1^{-}$. Indeed if we assume that about half of the moment trials lead to $\Delta E < 0$ where $\Gamma_{met} = 1$ then in average and still for small values of $\beta \Delta E$ the Metropolis acceptance rate is about $\Gamma_{met} \approx (1 - \beta |\Delta E|/2)$ while $\Gamma_{hb} \approx (1 - \beta \Delta E/2)/2 \sim 0.5$, we can extend equation (6) to $u_{t,met} \approx 2u_{t,hb}$ when using the Metropolis acceptance rate instead of the heat-bath one. For not sufficiently small values of $R_0$, some deviations may occur (particularly for the Metropolis sampling). In such cases we may extend equation (6) to $u_t = \Gamma_{acc} \Delta N R_0^2$ where $\Gamma_{acc}$ is the value of the acceptance rate for a given temperature. For a MC path performed at fixed $T^*$ as a relaxation process, $\Gamma_{acc}$ is obviously a constant. In our simulations, we limit ourselves to reasonably small values of $R_0$ for which $\Gamma_{acc}$ is very weakly dependent of $T^*$ in the range of temperatures considered in the MC study$^{30}$. We thus propose to define an operational unit of time, both $T^*$ and sampling choice independent,

$$ u_t = \bar{\Gamma}_{acc} \Delta N R_0^2 \quad (7) $$

where $\bar{\Gamma}_{acc}$ is an average over the temperatures considered, excluding the high temperature range where the dynamics become reversible and the time unit is no longer relevant. Notice that in the framework of the heat-bath acceptance rate and when $R_0 \to 0$, equations (6) and (7) coincide. A note should be made about the precise determination of $\bar{\Gamma}_{acc}$. For isolated particles at $T^* = 1$ and Metropolis scheme, $u_t$ is increased by a factor 3.4 when going from $R_0 = 0.2$ to 0.4. This contrasts with the variation of $\Gamma_{acc}$ in the significant range of temperature, excluding the reversible region: we get ca. 0.6% and 1.3% from $T^* = 0.7$ to 1 at $R_0 = 0.2$ and 0.4 respectively. The error induced in $u_t$ by the variation of $\Gamma_{acc}$ is thus negligible. The unit of time $u_t$ varies essentially through the simulation parameters, $\Delta N$ and $R_0$.

From this point, we have to distinguish between the instantaneous MC time step $i$, which is nothing but the index of Monte Carlo step (as usual one MC step corresponds to a trial move for all the moments) and the MC time $t_{MC}$, which includes the unit of time as defined by either equation (6) or (7), namely $t_{MC} = u_t i$. In most of our simulations, especially for relaxation process, we divide the whole number of Monte Carlo steps, say $N_m$, into intervals of $\Delta N$ playing the role of a finite $\Delta t$ time interval. Then we consider sets of $\Delta N$ steps which will be labelled by $i_{MC} \in [1, N_m / \Delta N]$ and define the MC time as $t_{MC} = i_{MC} u_t$.

This approach is compatible with the one of Melenev$^{14}$, where instead of the introduction of the unit of time, (6, 7), the number of MC steps $N_{rel,T_0}$ corresponding to the relaxation time $\tau_0$ for a given temperature $T_0$ is used to define dimensionless time: $t^* = i_{MC} / N_{rel,T_0} = t / \tau_0$. Note that $t_{MC}$ and $t^*$ are dimensionless times. In the next section, we will explain how it is possible to map this simulation time with a time in seconds by using experimental data of specific systems.
The observables considered in this work are, besides the total energy, on the one hand the magnetization $M(t_{MC})$ projected on the direction of the external field $\hat{h}$ as the thermal average of the instantaneous $m(i)$ over $\Delta N$ MC steps and a large number (100 or 200) realizations of the systems as the studied properties are "time dependent" :

$$M(t_{MC}) = \langle m(i_{MC}) \rangle_{\Delta N} = \frac{1}{\Delta N} \sum_{i=i_{MC}\Delta N}^{(i_{MC}+1)\Delta N} m(i) \quad \text{with} \quad m(i) = \frac{1}{N} \sum_{k=1}^{N} \hat{\mu}_k \cdot \hat{h}. \qquad (8)$$

On the other hand the susceptibility $\chi(\omega)$ in terms of the frequency $\omega$, determined from the method introduced by Andersson et al.⁶. For the latter, the external reduced field is taken as $h(t) = h_0 \sin(\omega t)$ and we have

$$\chi' = \frac{2}{h_0 N_p} \sum_{i=1}^{N_p} \sin\left(\frac{2\pi}{N_p} i\right) m(i) \quad \chi'' = \frac{2}{h_0 N_p} \sum_{i=1}^{N_p} \cos\left(\frac{2\pi}{N_p} i\right) m(i) \qquad (9)$$

where $N_p$ is the number of MC steps corresponding to one period, i.e. $N_p.u_t = 2\pi/\omega$. The period corresponds to the measuring time $t_m = N_p u_t$. On a practical point of view, the whole MC run involves $N_m$ MC steps covering a number $p$ of periods, the averages involved in equation (9) are performed over the last $p/2$ periods. Equation (9) is written for $\Delta N = 1$ for the sake of simplicity, but can easily be translated to the more general case of $\Delta N > 1$ as

$$\chi' = \frac{2}{h_0 N'_p} \sum_{i=1}^{N'_p} \sin\left(\frac{2\pi}{N'_p} i\right) \langle m(i\Delta N) \rangle_{\Delta N} \quad \text{with} \quad N'_p = \frac{N_m}{\Delta N} \frac{1}{p} \qquad (10)$$

and similarly for $\chi''$. In any case, when translated to the actual MC time, the value of $\Delta N$ does not play any role. From the magnetization $m(i_{MC})$, we will mainly look at the relaxation either in the absence of an external field and starting for instance from a state with all the moments in a fixed direction, or in the presence of a static external field $h_0$ with the system evolving to the equilibrium value $m_{eq}(h_0)$.

## IV. RESULTS

The characterization of the time evolution of dynamic properties calculated in the framework of the TQMC is done first through the relaxation of the magnetization $M(t)$ (Unless indicated, in what follows $t$ stands for the MC time $t_{MC}$). For this purpose we consider the ensemble of NP with all the easy axes either aligned in the $\hat{z}$ direction, textured case, or randomly distributed and with all the moments aligned along the $\hat{z}$ axis with $\mu_z = 1$ as the initial state. Then we let the system evolve in the absence of external field towards the equilibrium state characterized by $M_z = 0$ according to the standard Néel process.

We first check that $u_t$ defined in equation (6) corresponding to the heat-bath acceptance rate actually plays the role of a unit of time for the TQMC Monte Carlo. To this aim, we consider a set of relaxation runs characterized by the same value of $u_t$ obtained with different couples of values of $(\Delta N_{MCS}, R_0^2)$ for a typical case of DDI interacting particles (fig.(1)). A very similar result is obtained for non interacting particles. As one sees, the relaxation curves corresponding to a given value of $u_t$ coincide. In the same way, we make coincide the $M(t)$ corresponding to different values of the Monte Carlo unit of time $u_t$ when displayed in terms of $t = i_{MC}u_t$. (fig.(2)). At this point only the heat-bath acceptance rate for reasonable values of $R_0$ is used according to $^{10}$. In order to test our suggestion allowing to use also the Metropolis acceptance rate according to the extension from equation (6) to (7) we perform magnetization relaxations by using both
acceptance rates. See the left panel of figure (3). As expected the Metropolis algorithm is more efficient than the heat-bath one. However, once again we make the corresponding $M(t)$ coincide by using the unit of time given in equation (7), see the right panel of figure (3). Another example is given in figure (4) in the case of NP interacting through DDI where we also include the variation of $\Gamma_{acc}$ along the whole relaxation path. The latter is clearly constant although the relaxation does not follow a single exponential decay due to the DDI. A similar behavior has been mentioned in³¹ for non interacting MNP at least for $\Gamma_{acc} \ge 0.40$.

We also note from figures (3,4) that one can extend the value taken by $R_0$ at least up to ~ 0.3 in order to satisfy the dynamics. It is noteworthy that the efficiency is proportional to $R_0^2$. As an illustration, in figure (4). the paths corresponding to $R_0 = 0.05$ and 0.3 use 1000 and 200 steps of $\Delta = 20$ MCS respectively although the first one is ~ 6 times shorter in real time.

In order to clearly identify the relaxation time, we let evolve $M(t)$ up to some value, say $t_1$ and look at the relaxation of $M(t)/M(t_1)$ in terms of $(t - t_1)$. Doing this, we exclude the very short time relaxation process. We thus obtain the relaxation time $\tau$ from the fit of $M(t)/M(t_1)$ on $\exp(-(t - t_1)/\tau)$. An example is shown in figure (5). It is worth mentioning the longitudinal (inter-well) relaxation time, $\tau_\|$ is much larger than the transverse (intra-well) one, $\tau_\perp$ (see figure (6)). As a result, the relaxation time for a random distribution of easy axes is quite close to $\tau_\|$.

In order to determine how the longitudinal relaxation time depends on $\lambda_u$ and $T^*$ for non interacting MNP, two sets of relaxation simulation runs are performed. First, we consider a set of values of $\lambda_u$ at constant $T^* = 1$, which remains to examine different values of $\sigma = \lambda_u/T^*$ at constant $T^*$ and second, a set of values of $T^*$ with the same value of $\sigma$, namely with $\lambda_u = T^*\sigma$. From the latter, we clearly evidence that $\tau$ behaves as $\tau = \exp(\sigma)f(\sigma)/T^*$ in agreement with Refs.²⁵,²⁶. For the function $f(\sigma)$ our simulations fit well the analytical results of either Brown²⁵ or Bessais et al.²⁶, both valid for $\sigma \ge 2$ (here we consider only $\sigma > 2.6$), as is shown on figure (7), namely

$$ \tau = b_0 \frac{\exp(\sigma)}{\sigma^{3/2}T^*} \quad \text{or} \quad \tau = b_1 \frac{\exp(\sigma)}{(1 + \sigma/4)^{5/2}T^*} \quad (11) $$

respectively. We must notice that in the range of values of $\sigma$ considered here, we cannot distinguish clearly between these two analytical forms of $\tau(\sigma, T^*)$. The relaxation time $\tau$ given in equation (11) is still unitless until an actual value in second is associated to the constant $b_0$ or $b_1$. This is not straightforward, as already noted in²⁶ for instance, and is to be compared to the estimation of $\tau_0 \sim 10^{-8} - 10^{-10}s$ generally used for the Arrhenius form. Following²⁵ as well as²⁶, the constants $b_0$ or $b_1$ should depend on the NP properties as $\sim \beta_0 v M_s \alpha/\gamma$ where $\alpha$ and $\gamma$ are the damping and the gyromagnetic factors respectively. In any case such a correspondence would need the precise knowledge of the NP properties. In principle this could be achieved from the knowledge of $K_u$ and the experimental relaxation time, but from our knowledge the latter is not available. We consider that a better way to map the actual time to the MC time is to focus on the behavior of $T_{max}(\chi')$ or $T_{max}(\chi'')$ in terms of $\omega$ where $T_{max}(\chi')$ ($T_{max}(\chi'')$) is the temperature of the maximum of the AC susceptibility $\chi'$ ($\chi''$).

Now in the case of non interacting NP, we compare our simulation results for the AC susceptibility $\chi(\omega)$ to the analytical result obtained from the Shliomis and Stepanov model³²,³³ in terms of the reduced variable $y = \omega\tau$ where only the longitudinal relaxation time is retained and $\tau = \tau_\|$, given for a random distribution
of easy axes by

$$
\begin{equation}
\begin{gathered}
\chi' = \frac{\chi_0 + \chi_1 y^2}{1 + y^2} \qquad \chi'' = \frac{y\chi_\parallel}{1 + y^2} \qquad (12) \\
\text{with} \quad \chi_0 = \frac{1}{3T^*}; \quad \chi_1 = \chi_0 \left( \frac{1}{\sigma} + \frac{1}{2\sigma^2} + \frac{5}{4\sigma^3} \right); \quad \chi_\parallel = \chi_0 \left( 1 - \frac{1}{\sigma} - \frac{1}{2\sigma^2} - \frac{5}{4\sigma^3} \right)
\end{gathered}
\end{equation}
$$

where the limit $\sigma \gg 1$, which is valid down to $\sigma \simeq 3.5^{33}$, has been taken. In order to compare our simulated susceptibility with equation (12), $y$ is taken with $\tau$ from equation (11) according to Ref. 26 with the constant $b_1$ fitted from our TQMC simulations of the relaxation described above and $\omega = 2\pi(N_p u_t)^{-1}$ with $u_t$ given by equation (7). The agreement that we get between our simulated $\chi$ and the analytic results Ref. 32,33, give us confidence on the way we introduce the time unit $u_t$ since the latter involves 3 decades of frequency. Notice that the nice agreement we get between our simulations and equation (12) (see figure (8)) holds only for sufficiently small values of the field amplitude $h$ in order to avoid non linear terms in the response to the AC field. As an heuristic criterion, we get from the maximum value of $\chi'$, $h\chi'_{max} \lesssim 0.15$.

From $\chi'(\omega, T^*)$ and $\chi''(\omega, T^*)$ we can deduce the behavior of $T'_{max}(t_m)$ and $T''_{max}(t_m)$ which are the locations of the maxima of $\chi'$ and $\chi''$ in terms of $T^*$ at fixed values of the measuring time $t_m = 2\pi/\omega$. On figure (9), we give $T'_{max}(t_m)$ and $T''_{max}(t_m)$ resulting from both the analytical result of equation (12)32 and the MC simulations in the case of the random distribution of easy axes. The curves displayed in figure (9) show the dependence of $T^*_{max}$ on the measuring time $t_m$ for several values of $\lambda_u$. Note the time, temperature and energies remain expressed in terms of $u_t$ and $T_0$.

To map the TQMC time on a true time we need to compare these curves to the equivalent ones obtained from real AC susceptibility experiments for a sample of non-interacting MNP with given anisotropy $K_u v(d)$. In order to do this, we first note that once a reference temperature $T_0$ is chosen, $T^*$ is translated into a real one and $\lambda_u$ gives $K_u v(d)$. Then from the $T'_{max}^{exp}$ found in experiments for a given $\omega^{exp}$ measured in $s^{-1}$ we may find the corresponding MC frequency $2\pi/(N_p u_t)$ for which $T'_{max}$ coincides with $T'_{max}^{exp}$ for the appropriate value of $\lambda_u$. This allow to quantify the TQMC time unit in such a way that $u_t$ represents a dimensioned time according to :

$$ u_t = \frac{2\pi}{\omega^{exp}N_p} \qquad (13) $$

We emphasize that this mapping of $u_t$ is dependent of the specific system as it should be for the resulting mapping of either $b_0$ or $b_1$ since the relaxation time is obviously system dependent.

To test our mapping procedure on an actual experimental system is not straightforward since the AC measurements on truly isolated MNP assemblies are only few$^{34,35}$. In Refs.$^{34,35}$ the DDI are made tunable through a non magnetic SiO$_2$ layer of variable thickness. The usual criterion to quantify the importance of the interactions between MNP from AC measurements is the value taken by the so-called Mydosh coefficient, namely the slope $\Phi = -\Delta T_m(\omega)/(T_m(\omega)\Delta\ln(\omega))$, ($\simeq -\partial\ln(T_m(\omega))/\partial\ln(\omega)$ if the variation of $T_m$ remains weak) with $\Phi \gtrsim 0.1$ or noticeably smaller in the non interacting or interacting cases respectively.

As an illustration, we look at the experimental result of Ref. [34] where $T_m(\chi')$ and $T_m(\chi'')$ are reported for a likely non interacting MNP assembly thanks to the SiO$_2$ coating of the particles, when the thickness of the later is of 17 nm for maghemite NP of 8 nm in diameter. Using the anisotropy constant $K_u = 3.0 \ 10^4 \text{Jm}^{-3}$ given in [34] and choosing $\lambda_u = 6$ we get $T_0 = 97K$.

We quantify $u_t$ from equation (13) for $f = 600 \text{ Hz}$: with $T_0$ this leads to map $T'_m \simeq 52 \text{ K}$ or $T'_m = 0.535$ which is done both by MC and the Shliomis-Stepanov$^{32}$ approximation, which is found to fit perfectly
well the value of $T_m^{*'}$ obtained by the simulation, see figure (8). We then get for the MC measuring time $N_p.u_t = 5.05 \times 10^5$. Finally, we translate to the MC measuring time corresponding to 10 Hz and deduce the temperature $T_m^{*'}(\chi'(f = 10Hz)) = 0.385$ ($T_m' = 37K$) which is reasonably close to the experimental one ($T_m' \simeq 44 K$).

Similarly, from susceptibility curves given in Ref.³⁵, we estimate the values of $T_m'$ for the NP of diameter $d = 15$ nm at ca 185 K and 155 K for $f = 10^3$ and 10 Hz respectively. Using $k_{DC}^{eff} \simeq 3 \times 10^4 Jm^{-3}$ from³⁵, we get $T_0 = 318$ K with the choice $\lambda_u = 12$, leading to $T_{m}^{*'} = 0.58$ for $f = 10^3$ Hz. We then obtain the MC measuring time $N_p.u_t = 2.30 \times 10^9$ to map this reduced temperature value and from our method, where the approximation of equation (12) is used, we obtain $T_{m}^{*'} = 0.472$ to represent the temperature peak of $\chi'$ at 10 Hz, leading to $T_m' = 150$ K close to the experimental value ($\sim 155$ K).

However, the comparison must be taken with care since $\Phi$ still takes a rather small value ($\sim 0.04$) for the samples of Ref.³⁴ and we thus cannot totally rule out a very weak interaction. Similarly, in Ref.³⁵, $\Phi$ remains smaller than $\sim 0.08$ and the AC blocking temperature $T_b(\omega)$ as determined from the peak of $\chi''$ depends on the SiO₂ layer thickness $\Delta$ up to $\Delta \simeq d$ where $d$ is the NP diameter. Moreover, the comparison between the experimental susceptibility and either the MC simulated one or the result given by³² should include the polydispersity even of reduced range, through a distribution of $\lambda_u$, expected to broaden the peak.

Now we demonstrate how the quantification of time made with isolated particles holds when the coupling between particles is not negligible with respect to the MAE.

Fig (10) show that starting from an ordered configuration of the dipoles along the $z$-axis, the relaxation of the system at different temperatures is no longer exponential for both a short range coupling or a long range DDI. Nevertheless, it is possible to use the properties of isolated particles to define a unit of time for the simulations of the interacting system. Taking as a unit of time the relaxation time $N_{rel,T_0}$ of the isolated particle with $R_0$ at $T_0^*$, it is possible to obtain the same relaxation curves using different dynamics (i.e. Metropolis sampling with different values of $R_0$) but preserving the dependence in terms of temperature. In Fig (10 left) the relaxation is plotted as a function of a reduced time $t^* = i_{MC}/N_{rel,T_0}$. A similar plot should be obtained using the unit of time $u_t = R_0^2 \bar{\Gamma}_{acc}$. Note that the dominant parameter in the unit of time is $R_0^2$. We checked that with Metropolis algorithm and moderate $R_0$, the same results are obtained taking $\Gamma_{acc}$ either of the simulation for an isolated particle with $R_0^2$ or of one for the interacting particles.

This holds also for $\chi'(\omega)$ as shown in fig (11) where the same curves where obtained with different dynamical processes by adjusting $N_p$ with respect to $N_{rel,T_0}$. This allows to save computational time because short frequencies should be studied by using quick relaxing dynamics (with large $R_0$). Note however that only moderate values of $R_0$ should be considered and consistency tests should be performed. As an example, in Fig (10 left), for $T^* = 2$ the same relaxation is obtained for $R_0 = 0.2$ and $R_0 = 0.4$ but not for $R_0 = 1$.

Clearly in the presence of non negligible long range DDI, we get a relaxation according to a stretched exponential behavior as already found in³⁶ (see figure (10 right)). Besides this relaxation process in the absence of external field, we have also studied the DDI and Heisenberg interaction on the dynamics of $\chi'$ and $\chi''$. As it is well known the interactions between NP strongly reduce the decrease of $T_m(\chi')$ as the measuring time $t_m$ is increased. This is the consequence of the collective behavior expected at low temperature leading to a freezing of the system with a power law dependence

$$t_m \propto (T_m' - T_c)^{-z\nu} \qquad (14)$$

with $T_c$ the super-paramagnetic (SPM)/super spin-glass (SSG) transition temperature between the SPM and
the expected SSG phases. The decrease of $T'_m$ together with the power law dependence of $t_m$ according to equation (14) are well reproduced in our simulations (see fig. (12)).

V. CONCLUDING REMARKS

The purpose of this work is mainly to revisit the TQMC scheme in order to provide a way to map the simulation scheme on an actual time scale. This is done first by showing that the precise control of the temperature dependent cone amplitude of the moments displacements leads to the definition of a MC time scale which can be obtained by different sets of the simulation parameters ($\Delta N, R_0$). Moreover we show that the scheme is not limited to the heat-bath sampling. The limitations of our TQMC development concern mainly the limit of $R_0$ which must be smaller than a threshold value obtainable from a trial run. In our calculations, we have got $R_0 < 0.4$ beyond which we have found deviation in the dynamics when compared to those obtained at smaller $R_0$ and the same value of $u_t$. Notice that, in the framework of the Metropolis sampling, this limits the acceptance rate to $\Gamma_{acc} \gtrsim 0.7$ (see figure (4)).

Then we focus on the frequency dependent susceptibilities, $\chi'(\omega)$ and $\chi''(\omega)$, as their frequency behavior is well known and experimental results for isolated MNP are available. We emphasize that our scheme holds also both for long range dipolar and short range Heisenberg interactions. This work gives ways to compare different theoretical or simulation results on the same models concerning dynamical properties of MNP. It also helps to bridge the gap between experimental and simulation results.

The first application of our TQMC will be to study models of disordered ensembles of MNP where a SG transition is obtained from both equilibrium MC simulations and experiments and to compare the transition temperatures obtained from the usual equilibrium MC simulations, with the ones obtained from TQMC simulations of AC susceptibilities and FC-ZFC dynamical experiments.

VI. ACKNOWLEDGEMENTS

We thank Dr. A. Tamion for providing us with the data relative to the experimental results of susceptibility of Ref.³⁴. This work was granted access to the HCP resources of CINES under allocations 2024-A0160906180 made by GENCI, CINES, France. A.M., F.V. and V.R. acknowledge the support of the French Agence Nationale de la Recherche (ANR) through the project ANR NanoHype under grant number ANR-21-CE09-0043. J.J.A. thanks SCBI at University of Málaga for additional computer time.
<figure><img src="image_2.png" /><figcaption>FIG. 1. Relaxation of the magnetization from $M_z = 1$ in absence of external field. FCC lattice with 500 NP and easy axes distribution textured along $\hat{z}$, $\lambda_u = 4$, $\lambda_d = 1$ and $T^* = 1$. Two realizations for each values of $u_t$ (see equation (6)) with $u_t$ going from 0.25 to 12.5 as indicated. These are obtained with pairs $(\Delta N, R_0) = (200, 0.05)$ and (800, 0.025); (200, 0.10) and (3200, 0.025); (2000, 0.05) and (8000, 0.025); (200, 0.2) and (800, 0.10); (2500, 0.1) and (10000, 0.025) for $u_t = 0.25$ to 12.5 respectively.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIG. 2. Relaxation of the magnetization starting from $M_z = 1$, displayed in terms of the MC time with $u_t$ from equation (6). Heat bath sampling. FCC lattice with 500 NP, easy axes distribution textured along $\hat{z}$ and $\lambda_u = 4$. Left: non interacting system; right: interacting system with $\lambda_d = 1$. $T^* = 1$.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. Relaxation of the magnetization starting from $M_z = 1$, calculated from either the heat-bath (hb) or Metropolis (met) sampling acceptance rates and two values for $R_0$ as indicated, without (left) or with (right) the extension of the definition of the MC time unit when going from equations (6) to (7). Ensembles of non interacting particles with $N = 500$, random distribution of easy axes and $\lambda_u = 8$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. Relaxation of interacting NP with dipolar interactions for 5 different realizations of the unit of time $u_t$ from equation (7) in the Metropolis sampling. $\lambda_d = 0.5$, $\lambda_u = 6$, $T^* = 1$. Flat curves on the upper part: (left y-axis) evolution of the acceptance rate; curves on the lower part (right y-axis) magnetization $M(t)$ starting from $M(t) = 1$, in log scale.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 5. Relaxation of $M_z(t)$ in terms of t (in MC time units) relative to the value of $M_z(t_1)$ for ensembles of 1720 non interacting NP at $T^* = 1$ and $\lambda_u = 10, 9, 8, 7, 6$ and 5 from right to left. $t_1 = 19$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Magnetization relaxation in absence of external field, starting from $M_z(t = 0) = 1$ for a non interacting system with a textured distribution of easy axes along either the $\hat{z}$-axis and $\lambda_u = 5$ (top curve) or the $\hat{x}$-axis and $\lambda_u = 5, 7$ and $10$ (rest of the curves from left to right) probing the longitudinal $\tau_\|$ and transverse relaxation times $\tau_\perp$ respectively. $T^* = 1.0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. 7. Relaxation time $\tau$ in terms of $\sigma$ at $T^* = 1$ for non interacting NPM with random (circles) or textured along $\hat{z}$ (triangles) distribution of easy axes. $\tau$ is expressed in $u_t$ units from equation (7). Solid and dashed lines are the fits on the analytical results taken from Refs.$^{26}$ and$^{25}$ respectively.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 8. Comparison between the susceptibility calculated from the approximation of$^{32}$ with $\tau$ taken from$^{26}$ (solid lines) and our TQMC simulations (symbols) in the absence of DDI and with random distribution of the easy axes. Left (right) panel displays curves for $\chi'$ ($\chi''$). The curves are indexed by the actual MC measuring time (see text). In each panel, the left hand side curve ($t_m = 5.05 \times 10^5$) is for $\lambda_u = 6$ and the other ones for $\lambda_u = 8$. The product $\chi'_{max} \cdot h_0$ ranges from 0.048 to 0.17 and is seen as a criterion for the LRT regime to be valid ($\chi'_{max} \cdot h_0 \lesssim 0.15$).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>FIG. 9. MC measuring time $t_m$ in terms of the temperature $T_m^*$ of the maximum of $\chi'(\omega = 2\pi/t_m)$. Colored lines/symbols stand from results obtained from the approximation of $^{32}$. Black symbols stand for TQMC results. Circles, triangles and squares are for $\lambda_y = 4, 6, 8$ respectively. Non interacting system with random distribution of easy axes.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 10. Magnetization relaxation for ensembles of interacting particles. Left : short range Heisenberg interactions for a BCC lattice of 2000 NP with randomly distributed easy axes. $J = 1$ and $\lambda_u = 6$. Two reduced temperatures $T_0^* = 1$ (upper curves) and $T^* = 2$ (lower curves) and different dynamics with $R_0 = 1, 0.4$ and $0.2$. Time is given as $t^* = i/N_{rel,T_0}$.</figcaption></figure>

Right : long range DDI for a FCC lattice of 864 NP with easy axes distribution textured along $\hat{z}$, $T^* = 0.6$ and $\lambda_u = 6$. and different values of $\lambda_d$ as indicated. Time is in $u_t$ units. The black thin lines are fits according to the stretched exponential : $M_z(t > t_1) \propto \exp(-ct^\beta)$. $\beta = 0$ for $\lambda_d = 0$ and $0.2$, $\beta = 0.485, 0.258$ for $\lambda_d = 0.4$ and $0.6$ respectively.
<figure><img src="image_12.png" /><figcaption>FIG. 11. In phase susceptibility in terms of temperature $T^*$ for interacting NP (same model as in Fig. (10-left)) obtained from TQMC simulations for several values of the dimensionless measuring times $t_m^* = N_p/N_{rel,T_0}$. Curve displayed with continuous lines are from bottom to top for $t_m^* = 40.305, 48.366, 60.450$ and $80.610$ using $R_0 = 0.4$. Crosses and purple line compare AC susceptibility for $t_m^* = 40.305$ obtained from dynamics with $R_0 = 0.2$ and $R_0 = 0.4$ respectively. The vertical line corresponds to the SPM/SSG transition temperature obtained from equilibrium MC simulations.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. MC measuring time $t_m$ in terms of the temperature $T_{max}^*$ of the maximum of $\chi'(\omega = 2\pi/t_m)$ for interacting NP with DDI, for $\lambda_u = 4$ and the indicated values of $\lambda_d$ compared to the non interacting case. Data are obtained from TQMC simulations. The blue line is the result from$^{32}$. The vertical red bar indicates the SPM/SSG transition temperature obtained from the fit of equation (14).</figcaption></figure>
¹ S. Bedanta and W. Kleemann. *J. Phys. D*, 42:013001, 2009.

² E.C. Abenojar, S. Wickramasinghe, J. Bas-Concepcion, and A.C.S. Samia. *Prog. Nat. Sci.: Mater. Int.*, 26:440, 2016.

³ Jesus G. Ovejero, Federico Spizzo, M. Puerto Morales, and Lucia Del Bianco. *Materials*, 14:6416, 2021.

⁴ H. Gavilán, K. Simeonidis, E. Myrovali, E. Mazarío, O. Chubykalo-Fesenko, R. Chantrell, Ll. Balcells, M. Angelakeris, M. P. Morales, and D. Serantes. *Nanoscale*, 13:15631–15646, 2021.

⁵ A. Bordet, L.-M. Lacroix, P.-F. Fazzini, J. Carrey, K. Soulantica, and B. Chaudret. *Angew. Chem. Int. Ed.*, 55:15894, 2016.

⁶ J.-O. Andersson, P. Svedlindh C. Djurberg, T. Jonsson, and P. Nordblad. *Phys. Rev. B*, 56:13983, 1997.

⁷ D. Brinis, A. Laggoun, D. Ledue, and R. Patte. *J. Appl. Phys.*, 115:173906, 2014.

⁸ R. P. Tan, J. Carrey, and M. Respaud. *Phys. Rev. B*, 90:214421, 2014.

⁹ J. García-Otero, M. Porto, J. Rivas, and A. Bunde. *Phys. Rev. Lett.*, 84:167, 2000.

¹⁰ U. Nowak, R. W. Chantrell, and E. C. Kennedy. *Phys. Rev. Lett.*, 84:163, 2000.

¹¹ X. Z. Cheng, M. B. A. Jalil, and H. K. Lee. *Phys. Rev. B*, 73:224438, 2006.

¹² O. Chubykalo, U. Nowak, R. Smirnov-Rueda, M. A. Wongsam, R. W. Chantrell, and J. M. Gonzalez. *Phys. Rev. B*, 67:064422, 2003.

¹³ O. Billoni and D.A. Stariolo. *J. Magn. Magn. Mat.*, 316:49, 2007.

¹⁴ P. V. Melenev, Yu. L. Raikher, V. V. Rusakov, and R. Perzynski. *Phys. Rev. B*, 86:104423, 2012.

¹⁵ D.A. Stariolo and O.V. Billoni. *J. Phys. D: Appl. Phys.*, 41:205010, 2008.

¹⁶ D. Serantes, D. Baldomir, C. Martinez-Boubeta, K. Simeonidis, M. Angelakeris, E. Natividad, M. Castro, A. Mediano, D.-X. Chen, A. Sanchez, LI. Balcells, and B. Martínez. *J. Appl. Phys.*, 108:073918, 2010.

¹⁷ R. Fu, Y. Yan, C. Roberts, Z. Liu, and Y. Chen. *Sci. Rep.*, 8:4704, 2018.

¹⁸ C. H. Back, R. Allenspach, W. Weber, S.S.P. Parkin, D. Weller, E.L. Garwin, and H.C. Siegmann. *Science*, 285:864, 1999.

¹⁹ O. Hovorka, J. Barker, G. Friedman, and R. W. Chantrell. *Phys. Rev. B*, 89:104410, 2014.

²⁰ D. Serantes, V. Vega, W. O. Rosa, V. M. Prida, B. Hernando, M. Pereiro, and D. Baldomir. *Phys. Rev. B*, 86:104431, 2012.

²¹ D.A. Dimitrov and G.M. Wysin. *Phys. Rev. B*, 96:9237, 1996.

²² L. Wang, J. Ding, H.Z. Kong, and Y. P. Feng. *Phys. Rev. B*, 64:214410, 2001.

²³ S. Lepadatu. *Journal of Applied Physics*, 130(16):163902, 2021.

²⁴ H.-M. Nguyen and P.-Y. Hsiao. *Appl. Phys. Lett.*, 95:222508, 2009.

²⁵ W. F. Brown. *Phys. Rev.*, 130:1677, 1963.

²⁶ L. Bessais, L. Ben Jaffel, and J.L. Dormann. *Phys. Rev. B*, 45:7805, 1992.

²⁷ M. Itakura. Frozen quasi-long-range order in the random anisotropy heisenberg magnet. *Phys. Rev. B*, 68:100405, 2003.

²⁸ M. P. Allen and D. J. Tildesley. *Computer Simulation of Liquids*. Oxford Science Publications, 1987.

²⁹ J.-J. Weis and D. Levesque. *Phys. Rev. E*, 48:3728, 1993.

³⁰ We note that in our systems of fixed NP with the moments as unique degree of freedom, by using the usual Metropolis sampling in the Monte Carlo simulation and imposing a constant acceptance rate through a self consistent adjustment of the displacement cone $R$, one obtains as a general rule $R \sim \sqrt{T}$ at least in the non reversible temperature range and sufficiently large values of the acceptance rate, which corroborate $\Gamma \sim \text{cste}$ with $R = R_0\sqrt{T}$.

³¹ J.C. Zapata and J. Restrepo. *Computation*, 9:124, 2021.

³² M.I. Shliomis and V.I. Stepanov. *J. Magn. Magn. Mat.*, 122:176, 1993.

³³ J.G. García-Palacios. *Adv. Chem. Phys.*, 112:1, 2000.
<footnote n="34">34 J. A. De-Toro, P. S. Normile, S. S. Lee, D. Salazar, J. L. Cheong, P. Muniz, J. M. Riveiro, M. Hillenkamp, F. Tournus, A. Tamion, and P. Nordblad. *J. Phys. Chem. C*, 117:10213, 2013.</footnote>

<footnote n="35">35 P. García-Acevedo, M. González-Gómez, A. Arnosa-Prieto, L. de Castro-Alves, Y. Piñeiro, and J. Rivas. *Adv. Sci.*, 10:2203397, 2023.</footnote>

<footnote n="36">36 O. Petracic, X. Chen, S. Bedanta, W. Kleemann, S. Sahoo, S. Cardoso, and P.P. Freitas. *J. Magn. Magn. Mat.*, 300:192, 2006.</footnote>