Self-Adaptive Acceptance Rate-Driven Markov Chain Monte Carlo Method Applied to the Study of Magnetic Nanoparticles

Abstract

A standard canonical Markov Chain Monte Carlo method implemented with a single-macrospin movement Metropolis dynamics was conducted to study the hysteretic properties of an ensemble of independent and non-interacting magnetic nanoparticles with uniaxial magneto-crystalline anisotropy randomly distributed. In our model, the acceptance-rate algorithm allows accepting new updates at a constant rate by means of a self-adaptive mechanism of the amplitude of Néel rotation of magnetic moments. The influence of this proposal upon the magnetic properties of our system is explored by analyzing the behavior of the magnetization versus field isotherms for a wide range of acceptance rates. Our results allows reproduction of the occurrence of both blocked and superparamagnetic states for high and low acceptance-rate values respectively, from which a link with the measurement time is inferred. Finally, the interplay between acceptance rate with temperature in hysteresis curves and the time evolution of the saturation processes is also presented and discussed.

1. Introduction

The theoretical study of magnetic nanoparticle systems dates to the pioneering work of E. C. Stoner and E. P. Wohlfarth. (1948) [1], L. Néel (1949) [2] and W. J. Brown (1963) [3]. These works set the starting point for current developments and applications in the field of magnetic fluids, which include magnetic resonance imaging, magnetic hyperthermia for cancer treatment, among others. [4,5,6,7].

Due to the mathematical complexity of systems composed of many particles, it is necessary to implement numerical simulations carried out by computer, through algorithms and simulation methods to recreate their behaviors. For magnetic nanoparticle systems, the stochastic differential Landau–Lifshitz–Gilbert (LLG) [8,9] equation or the respective Fokker–Planck (FP) [10] equation, are usually integrated to reproduce the movement of magnetic moments and the appropriate probability distribution. On the other hand, some authors prefer to use Monte Carlo (MC) simulations based on Metropolis–Hastings (MH) dynamics for this purpose [11,12]. Monte Carlo methods, as is well established, can be based on sampling of discrete events or on Markov chains. This latter is known as Markov chain Monte Carlo (MCMC), from which the MH algorithm is the most popular MCMC method to generate Markov chains according to a certain proposal probability distribution. In a classical physical system of magnetic moments in contact with a thermal reservoir, such a distribution is given by the Maxwell-Boltzmann statistics. The MCMC method, which uses the Bayesian inversion approach, has been demonstrated to be a powerful tool to estimate unknown observables according to a prior knowledge as it can be found in several reported works [13,14,15,16,17].

Any of these methods is really feasible; however, the MC technique differs from LG or FP in that it lacks a time variable in real units, which may in principle generate limitations when one wishes to tackle kinetic or dynamic properties.

In the search for MC algorithms that allow the inclusion and reproduction of dynamic quantities, some interesting proposals have emerged; among them is the one by P. V. Melenev et al. (2012) [18]. In their work, authors propose that the evolution of the magnetic moments is carried out using the conventional Metropolis–Hastings algorithm, additionally, special rules are imposed on the rotations that they can perform to recreate metastable states that allow access to magnetic properties such as hysteresis. They propose, for example, that the Monte Carlo step (MCS) takes the role of time (t) and that the rotations (of a random nature) are obtained from angles taking values between 0 and $\delta \theta$, being the angular parameter $\delta \theta$ chosen ad hoc. Calibrating appropriately MCS and $\delta \theta$ they obtain results comparable with the LLG (or FP) equation.

Another idea that also stands out is the one put forward by D. A. Dimitrov and G. M. Wysin (1996) [19]. They use a model similar to Melenev but they force the algorithm to accept and reject magnetic moment movements at a certain constant rate, which they call the acceptance rate $\mathsf{\Gamma }$. This control over rotations forces to $\delta \theta$ to be adjusted. The authors state that in this way it is possible to sample the phase space at a uniform rate, simulating dynamic properties, and be assertive that the Monte Carlo step can be viewed as a time variable. To test these ideas, they obtain Zero Field Cooling (ZFC) and Field Cooling (FC) curves, and calculate the blocking temperature for a system of Cobalt nanoparticles. As a result, they show that their thoughts are like-minded with the experiment; moreover, they obtain a transformation of MCS to t.

In summary, there is an alternative and promising method that can be directly comparable with LLG, FP and experimental results, which we believe should be thoroughly studied. Therefore, this paper aims to build on this research to investigate the role of the acceptance rate and how actually affects the properties of a magnetic nanoparticles ensemble. We choose the magnetization and study its response in the presence of magnetic fields at constant temperature. Curves are simulated and computed for different values of $\mathsf{\Gamma }$. The influence of this quantity on the behavior of the system is analyzed and we conclude that the acceptance rate plays an important role in the relaxation processes.

Finally, we show that the model and computational method used can recreate the dynamics of magnetic moments under Néel rotations (magnetic moment rotates internally with respect to the mono-crystalline lattice, see Section 2.2). The method can be extended to include Brownian motion present in realistic magnetic fluids. This can be done by implementing translations and rotations of the particles in the presence of viscous media. This is to develop a computational work focused on the calculation of the Specific Loss Power (SLP) for biomedical applications in alternative cancer treatments using the magnetic hyperthermia technique.

2. Materials and Methods

2.1. Magnetic Nanoparticle Model

A physical system composed of a set of spherical nanoparticles of equal size, sufficiently spaced so that the interactions between them are negligible and uniformly distributed in a solid matrix in such a way that translations and rotations are forbidden, is considered. Each particle (see Figure 1) is assumed to be a single-magnetic-domain with uniaxial magneto-crystalline anisotropy; as a result, its properties can be characterized through a magnetic moment $\overrightarrow{\mu }$ and an easy axis $\widehat{n}$.

The magnitude of the magnetic moment $\mu$ is equal to the product of $M_S$ (saturation magnetization of the single-domain) with $\mathsf{\Omega }$ (nanoparticle volume). It is supposed that thermal fluctuations do not change $M_S$ and neither dilate the particle. Likewise, the direction of the easy axis is not affected by the thermal bath nor by internal or external interactions, which means that its initial orientation remains unchanged at all times.

With respect to the foregoing, the simplest Hamiltonian governing the behavior of a magnetic moment in the presence of an external magnetic field $\overrightarrow{H}$ is:

$$\mathcal{H}=-K_{eff}\mathsf{\Omega }{\left(\frac{\overrightarrow{\mu }·\widehat{n}}{\mu }\right)}^2-{\mu }_0\overrightarrow{\mu }·\overrightarrow{H},$$

(1)

where the first term is the anisotropy potential energy that describes the interaction between the magnetic moment and the easy magnetization axis, with $K_{eff}$ the effective magnetic anisotropy constant. Surface and shape contribution to the anisotropy are neglected for simplicity. The second term is the Zeeman potential energy and expresses the coupling between the magnetic moment and the field. Temperature is included using the Metropolis algorithm as indicated in Section 2.3.

2.2. Néel Relaxation

Due to the existence of uniaxial magnetic anisotropy, it is evident from Equation (1) that the magnetic moment has two stable and energetically equal orientations (one of them parallel to the easy axis and the other anti-parallel). These two orientations are separated by an energy barrier equal to $K_{eff}\mathsf{\Omega }$ and for a given absolute temperature T, there is a probability that the moment spontaneously changes from one direction to the other because of thermal fluctuations. The average time for this change to occur is known as the Néel relaxation time [20] and it is given by the expression:

$${\tau }_N={\tau }_{0}exp\left(\frac{K_{eff}\mathsf{\Omega }}{k_{B}T}\right),$$

(2)

with ${\tau }_0$ a characteristic time of the magnetic solid that takes typical values of ${10}^{-9}$ s (or less) [21] and Boltzmann’s constant $k_B$. This equation is only valid for zero or very weak external magnetic fields compared to the anisotropy interaction.

Suppose now that we want to measure the magnetization of a nanoparticle and that the measurement takes a time ${\tau }_m$ to be completed. If the measurement time is greater than the Néel time (i.e., ${\tau }_m>{\tau }_N$) the result is that the magnetic moment oscillates many times during measurement and therefore the average magnetization is zero (see Figure 2a). Conversely, if this time is smaller than the relaxation time (${\tau }_m<{\tau }_N$) then the magnetic moment will not oscillate, and it will remain in a blocked state causing the magnetization of the nanoparticle remains well defined (see Figure 2b). The first situation is called superparamagnetic state since the magnetization behaves like that of a paramagnet with macroscopic magnetic moment. The second one is known as a blocked state because the moment remains aligned along a given orientation.

2.3. Metropolis Algorithm Driven by Acceptance Rate

The Metropolis–Hastings algorithm is used to carry out magnetic moment movements and updates along the phase space of the system. For this algorithm, a single Monte Carlo step (MCS) is defined as the elapsed time to visit all the particles (N) in the system and to attempt a new magnetic moment configuration (${\overrightarrow{\mu }}^{\prime }$) from a previous one ($\overrightarrow{\mu }$) for every single magnetic moment, by forming in this way the Markov chain. In this scheme, a particle is chosen at random, and the direction of its magnetic moment is perturbed to access a new configuration. If the energy difference between microstates is $\mathsf{\Delta }E\le 0$ then the movement is accepted, conversely, if $\mathsf{\Delta }E>0$ the movement is accepted only with a Boltzmann transition probability given by $W=exp(-\mathsf{\Delta }E/k_{B}T)$ by comparison with a random number $r\in [0,1)$. If $r\le W$ the motion is accepted, otherwise is definitively rejected and a new random particle is again selected. The goal of this algorithm, accordingly to a canonical ensemble in contact with a thermal bath as is known in statistical physics, is to make the magnetic moments obey a Boltzmann probability distribution given by:

$$\rho (\mathbb{X})=\frac{1}{Z}{e}^{-E(\mathbb{X})/k_{B}T},$$

(3)

with E being the eigenvalue of the Hamiltonian of the system (see Equation (1)), Z is the partition function, which serves as normalizing factor, and the vector $\mathbb{X}$ in phase space stands symbolically for the set of variables describing the degrees of freedom, e.g., $\mathbb{X}=(\overrightarrow{{\mu }_1},\overrightarrow{{\mu }_2},\cdots ,\overrightarrow{{\mu }_N})$. If any of the magnetic moments change from $\overrightarrow{\mu }\to {\overrightarrow{\mu }}^{\prime }$, a new microstate $\mathbb{Y}$ is generated. Thus, the probability density given by Equation (3) describes the statistical weight with which the configuration or microstate $\mathbb{X}$ occurs in thermal equilibrium [22].

Bayes theorem states that the posterior distribution is given by [15,23]:

$$\pi (\overrightarrow{\mu }|{\overrightarrow{\mu }}^{\prime })=\frac{\pi ({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu }){\pi }_0(\overrightarrow{\mu })}{\int \pi ({\overrightarrow{\mu }}^{\prime },\overrightarrow{\mu }){\pi }_0(\overrightarrow{\mu })d\overrightarrow{\mu }}.$$

(4)

Here, $\pi$ is a Lebesgue density, so ${\pi }_0(\overrightarrow{\mu })$ is the prior density. In our case, let $Q({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu })$ or $Q(\mathbb{Y}|\mathbb{X})$ be the proposal or instrumental density, i.e., a Markov transition density describing how to go from state $\mathbb{X}$ to $\mathbb{Y}$.

To sample phase space, we start from a current or initial state for the magnetic moment $\overrightarrow{\mu }$ with any given random orientation. The trial movement is executed so that the new orientation is close to the old one, with ${\overrightarrow{\mu }}^{\prime }$ selected inside a cone whose main axis is along $\overrightarrow{\mu }$ (see Figure 3). This update is achieved by means of a double random rotation matrix $R(\theta ,\varphi )$ (visit Appendix A.1) involving a polar angle $\theta \in [0,\delta \theta ]$ and an azimuthal one $\varphi \in [0,2\pi )$, being $\delta \theta$ the cone aperture. Angles $\theta$ and $\varphi$ are chosen randomly from uniform distributions $U(0,1)$. Thus, ${\overrightarrow{\mu }}^{\prime }=R(\theta ,\varphi )\overrightarrow{\mu }$ and the way as the angles are chosen, makes the Markov matrix Q to be symmetric, i.e., $Q({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu })=Q(\overrightarrow{\mu }|{\overrightarrow{\mu }}^{\prime })$.

A sketch of the algorithm is given in Algorithm 1, where $\alpha$ is called the acceptance/rejection probability. With this election, the principle of detailed balance is fulfilled and the ratio of transition probabilities $\pi ({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu })/\pi (\overrightarrow{\mu }|{\overrightarrow{\mu }}^{\prime })$ depends only on the energy change $\mathsf{\Delta }E$.

Algorithm 1 Metropolis–Hastings Algorithm.

Given the current state $\overrightarrow{\mu }$, then:

Generate $\theta \sim U(0,\delta \theta )$ and $\varphi \sim U(0,2\pi )$;      ▹ Random Uniform Distribution.  Propose the new candidate ${\overrightarrow{\mu }}^{\prime }\sim Q({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu })$;            ▹ See Appendix A.1.  Calculate $\alpha =min\{1,\frac{\rho ({\overrightarrow{\mu }}^{\prime })}{\rho (\overrightarrow{\mu })}\frac{Q(\overrightarrow{\mu }|{\overrightarrow{\mu }}^{\prime })}{Q({\overrightarrow{\mu }}^{\prime }|\overrightarrow{\mu })}\}$;       ▹ Acceptance/Rejection Probability.  Draw $r\sim U(0,1)$;                 ▹ Random Uniform Distribution.     if ($r\le \alpha$) then   Accept ${\overrightarrow{\mu }}^{\prime }$, set $\overrightarrow{\mu }={\overrightarrow{\mu }}^{\prime }$ and define the flag $\gamma =1$;  else   Reject ${\overrightarrow{\mu }}^{\prime }$, $\overrightarrow{\mu }$ does not change and do $\gamma =0$;  end if

Return $\overrightarrow{\mu }$ and $\gamma$.

In discrete event-based Monte Carlo simulations on classical magnetic systems, new orientations ${\overrightarrow{\mu }}^{\prime }$ of the magnetic moments are chosen completely at random, i.e., $\theta \in [0,2\pi ]$ instead of $\theta \in [0,\delta \theta ]$, and independent on the initial state $\overrightarrow{\mu }$. The consequence of this is that the system always exhibits a superparamagnetic behavior regardless of the temperature value, as Dimitrov and Wysin [19] warned in their paper, since the system can rapidly escape from metastable states responsible for magnetic hysteresis. Please note that $\delta \theta$ is a parameter that has not been specified so far, and it is the one responsible for controlling the convergence rate of the algorithm and how the exploration of the phase space is performed. If it is too small most of the moves will be accepted and vice versa.

The undesired effect of having the system always in a superparamagnetic regime can be overcome precisely through a proper handling of the $\delta \theta$ parameter. To do so we try to reproduce a dynamic similar to that of the LLG framework (where the system always evolves and explores likely microstates), by stating that $\delta \theta$ must be modified in a self-adaptive manner such that the phase space is sampled at a constant rate. This is achieved by including an additional acceptance rate ${\mathsf{\Gamma }}_{\theta }$, which is calculated by counting the number of accepted movements in a large enough number of MC steps $N_{MC}$. ${\mathsf{\Gamma }}_{\theta }$ is thus calculated and $\delta \theta$ is updated every single Monte Carlo step (MCS). Thus, the cone aperture $\delta \theta$ is adjusted so that statistically ${\mathsf{\Gamma }}_{\theta }$ remains constant as much as possible within certain tolerance range during the simulation of any curve. For such purposes, it is imposed that all the particles are selected and attempted to be perturbed so that a specific $\delta \theta$ influences equally the behavior of the set (see Algorithm 2).

Algorithm 2 Main Algorithm.

Set the initial conditions: magnetic field, temperature, magnetic moments orientations, the easy axes orientations and so on.

$N_{acc}=0$;                           ▹ Total Accepted Movements.    for$i\le N$do   Run Metropolis algorithm and return $\gamma$ value;               ▹ See Algorithm 1.   $N_{acc}=N_{acc}+\gamma$;  end for   Compute ${\mathsf{\Gamma }}_{\theta }=\frac{N_{acc}}{N}\times 100\%$;                     ▹ Acceptation Rate.  Update $\delta \theta$;                              ▹ See Algorithm 3.  Compute the physical properties.

The self-adaptation process of $\delta \theta$ is as follows (see Algorithm 3): an initial value for ${\mathsf{\Gamma }}_{\theta }$ is chosen, namely a ${\mathsf{\Gamma }}_0$ corresponding to the target acceptance rate we pretend to achieve, and then the acceptance rate is calculated after a Monte Carlo step. If ${\mathsf{\Gamma }}_{\theta }>{\mathsf{\Gamma }}_0+2\%$, which means more accepted movements (this happens when $\delta \theta$ is small), then $\delta \theta$ is increased by $20\%$ taking care that the maximum value of $\pi$ is not exceeded. On the other hand, if ${\mathsf{\Gamma }}_{\theta }<{\mathsf{\Gamma }}_0-2\%$, which means less accepted movements (this happens when $\delta \theta$ is big), then $\delta \theta$ is reduced by $20\%$. We can do this because Q is not uniquely determined and some arbitrariness in the explicit choice of it remains.

Algorithm 3 Cone Aperture Update.

if (${\mathsf{\Gamma }}_{\theta }<{\mathsf{\Gamma }}_0-2$) then   $\delta \theta =0.8\delta \theta$;              ▹ Reduce the cone aperture by 20%.  else if (${\mathsf{\Gamma }}_{\theta }>{\mathsf{\Gamma }}_0+2$) then   $\delta \theta =1.2\delta \theta$;             ▹ Increase the cone aperture by 20%.  else   $\delta \theta$ does not change;  end if   $\delta \theta =min(\pi ,\delta \theta )$.           ▹ The cone aperture cannot exceed $\pi$.

With the election of such percentages, i.e., with $2\%$ for the confidence interval of ${\mathsf{\Gamma }}_{\theta }$ and $20\%$ for the increase/decrease of the cone aperture $\delta \theta$, we managed to guarantee ${\mathsf{\Gamma }}_{\theta }\approx const$. For other percentages, in our tests, a stability in the acceptance rate was not achieved. In Appendix A.2 we show the results of some tests and diagnostics over this algorithm, including autocorrelation plots, effective sample size and thin size.

3. Results and Discussion

In our simulation we used a total number of $N={10}^3$ particles with magnetic anisotropy constant $K_{eff}=10$ kJm${}^{-3}$ and radius of $R=7$ nm. A saturation magnetization per particle $M_S=446$ kAm${}^{-1}$ corresponding to magnetite [24] and $M_0=N{M}_S$ for the ensemble were employed. Additionally, for all experiences the easy axes were randomly oriented. Numerical experiments are explained and discussed below.

3.1. Magnetization, Acceptance Rate and Cone Aperture vs. Magnetic Field

This section explores the dependence of magnetization both as function of a time-dependent field and the respective MCS dependence when a constant magnetic field is applied. At the beginning of each time-dependent field experience, all the magnetic moments are randomly oriented, and they are thermalized until saturation. To accomplish this initial state, 2000 MCS are required. For curves such as those shown in Figure 4, 200 points are used and a magnetic field in the z-direction is chosen according to the rule:

$$H_i=H_{0}cos\left(\frac{\pi }{200N_{MC}}i\right),$$

(5)

with $H_0\approx 40$ kAm${}^{-1}$ (or 500 Oe) the field amplitude and the time variable $i\in [0,200]N_{MC}$ given in MCS units. This particular choice is made since the magnetic hysteresis cycles are symmetrical and it is enough to calculate the upper branch, so a half oscillation in the field is only necessary, whereas the lower branch is obtained by simultaneously reflecting the first one with respect to the H (or x) and M (or y) axes. At each i-th time value, a constant field is assumed, and the evolution is performed using the Metropolis algorithm. Half of the total MCS, namely $N_{MC}/2$, are used to thermalize the system at a given new field value and the other half is used for data collection and averages calculation purposes. A total number of $N_{MC}=100$ MCS was employed at every point per experiment.

Figure 4 shows the results for (a) the reduced magnetization $M/M_0$, (b) the acceptance rate and (c) the cone aperture as a function of the reduced magnetic field $H/H_0$ at $T=100$ K. Each point on the curve of magnetization is the result of averaging the magnetic moments of the particles in the field-direction during data acquisition. The same averaging process is performed for ${\mathsf{\Gamma }}_{\theta }$ and $\delta \theta$ (see Section 2.3). All curves are in addition the result of configurational averages performed over 100 independent numerical experiments.

As is observed in Figure 4a, the system exhibits hysteresis for values of ${\mathsf{\Gamma }}_{\theta }$ above $20\%$, whereas for $10\%$ of acceptance the system behaves as a superparamagnet. As the acceptance rate increases above $20\%$ both the coercivity and the squareness of the loop increase. It is remarkable the degree of control on the acceptance rate ${\mathsf{\Gamma }}_{\theta }$ above $20\%$ as can be observed in Figure 4b, which is constant for all the range of values of the external field. Only for $10\%$ of acceptance a noticeable variation is observed only at low fields. In Section 3.2 we discuss this feature where both the temperature and the fact of having the system in a superparamagnetic state are key factors to understand such phenomenology.

As concerns to Figure 4c, for higher values of ${\mathsf{\Gamma }}_{\theta }$, the angular parameter of the cone aperture $\delta \theta$ is smaller. This fact means that the system remains confined at metastable states as long as the movements are strongly bounded. This is clearly demonstrated by looking at the cycle for ${\mathsf{\Gamma }}_{\theta }=90\%$ where the hysteresis is the biggest one and the coercive force is also the strongest one, which is a clear indicative of the high degree of metastability. Likewise, the peaks observed for $\delta \theta$ are the ones where the switching fields take place. For such fields, the cone aperture becomes greater as a physical mechanism to allow magnetization reversal. In addition, remanence (at zero field) is due to the coupling of the magnetic moments along that component of the easy axes favoring the direction of the applied field. On the other hand, for low ${\mathsf{\Gamma }}_{\theta }$ values, the aperture cone is greater, which means that between step and step, the magnetic moments can rotate more freely, performing broader trajectories in phase space and therefore they can easily escape from metastable states so the magnetic moments remain in an alternating regime giving rise to a zero average magnetization (see Section 2.2). As a result, no hysteresis is observed, and the superparamagnetic state is finally established.

The progressive disappearance of coercivity as the acceptance rate decreases, is consistent with a continuous transition from a blocked state to a superparamagnetic one. In this sense, blocked states are characterized by small values of the cone aperture $\delta \theta$, which is in agreement with what is expected for a so-called “blocked” state. In contrast, at zero coercivity (${\mathsf{\Gamma }}_{\theta }=10\%$), the widest cone apertures are consistent to what is expected for a “superparamagnetic” regime. Such a gradual disappearance of coercivity coincides with typical observations of the behavior of $M(H)$ curves as a function of the measurement time ${\tau }_m$, in which the state of the nanoparticle (superparamagnetic or blocked) depends on the comparison between the measurement time and the Néel relaxation time ${\tau }_N$ (see Section 2.2). Thus, at first glance, for a given system, our acceptance rate ${\mathsf{\Gamma }}_{\theta }$, which is indeed a frequency of acceptance, seems to be related to the reciprocal of the measurement time, i.e., ${\mathsf{\Gamma }}_{\theta }\sim 1/{\tau }_m$. This fact is endorsed by the Néel-Arrhenius relationship observed between the remanence ($M_r$) and $1/{\mathsf{\Gamma }}_{\theta }$ for ${\mathsf{\Gamma }}_{\theta }\ge 30\%$, i.e., the expression $M_r(t)=M_0{e}^{-t/{\tau }_N}$ is fulfilled by taking the variable $t=1/{\mathsf{\Gamma }}_{\theta }$. We want to recall at this point that the maximum number of MCS was kept constant. As we show below, very small values of ${\mathsf{\Gamma }}_{\theta }$ imply a high rate of rejections, and therefore phase space sampling is not so efficient, and averages are not reliable.

Another interpretation could be in turn one for which the ${\mathsf{\Gamma }}_{\theta }$ parameter goes proportional to ${\tau }_N$. i.e., ${\mathsf{\Gamma }}_{\theta }\sim {\tau }_N$ and the measurement time ${\tau }_m$ refers to the number of MCS used for data acquisition. Thus, if the MCS number is kept constant, as in our case, the interplay between blocked and superparamagnetic states can be achieved by tuning ${\tau }_N$ or equivalently ${\mathsf{\Gamma }}_{\theta }$. Thus, high values of ${\mathsf{\Gamma }}_{\theta }$ give rise to the compliance of ${\tau }_m\ll {\tau }_N$ favoring the appearance of blocked states, as is effectively observed. Conversely, small values of ${\mathsf{\Gamma }}_{\theta }$ will favor the inequality ${\tau }_m\gg {\tau }_N$ and superparamagnetic states will be observed. In both scenarios, the role played by ${\mathsf{\Gamma }}_{\theta }$ recovers the $M(H)$ behavior observed in the superparamagnetic theory.

Concerning the use of a constant field and the time dependence of the observables, results are summarized in Figure 5. In this case, the time dependence in MCS units is obtained for (a) the reduced magnetization, (b) the acceptance rate and (c) the cone aperture angle at $T=100$ K. In this case, the time evolution of the magnetization process was obtained at a constant field of 200 kAm${}^{-1}$ (or 2500 Oe) starting from an initial configuration of magnetic moments randomly distributed. Curves are in turn the configurational average over 100 experiments statistically independent.

Figure 5a shows the effect of the acceptance rate in the magnetization process under an applied field. Concretely, for small values of ${\mathsf{\Gamma }}_{\theta }$, the system can easily reach the saturation state within a small range of MCS. This behavior can be ascribed to the fact that more microstates compatible with those where the alignment of the magnetic moments with the field are generated as the cone aperture is increasingly wider. In particular, for ${\mathsf{\Gamma }}_{\theta }$ up to $40\%$, magnetization curves seem to overlap, and discrepancies in the magnetization mechanism, are only observable for greater percentages. Moreover, the system progressively becomes magnetically harder making the saturation state more difficult to reach as the acceptance rate increases.

Figure 5b shows the time evolution of the ${\mathsf{\Gamma }}_{\theta }$ parameter. In all cases, the system seems to initiate its dynamics at ${\mathsf{\Gamma }}_0\approx 50\%$, related to the random initial configuration of the magnetic moments. Once the dynamics starts, ${\mathsf{\Gamma }}_{\theta }$ evolves toward the expected values ${\mathsf{\Gamma }}_0$, initially determined, reaching steady states within the first 20 MCS. Analogously, the time evolution of the cone aperture angle $\delta \theta$ is shown in Figure 5c. Here, an initial value of $\delta \theta =45$° was defined. In this case, convergence is also achieved but in a different fashion, passing through a well-defined peak, dealing with the particular mechanism the system follows in phase space during the magnetization process for every ${\mathsf{\Gamma }}_{\theta }$ value considered.

3.2. Superparamagnetic Behavior and Temperature Dependence

Figure 6 summarizes the circumstances under which both the blocked behavior (that one characterized by open hysteresis loops drawn with solid lines) and the superparamagnetic state (dashed lines with no hysteresis) take place at different temperatures. Results are shown in Figure 6a–c for ${\mathsf{\Gamma }}_{\theta }$ percentages of $10\%$, $50\%$ and $90\%$, respectively. In Figure 6d,e the corresponding behaviors of the acceptance rate and the cone aperture angle are respectively plotted.

Temperature is a key factor in determining the behavior of acceptance rate and the cone aperture. This is clearly demonstrated in Figure 6d,e. As temperature increases along with its respective fluctuations, new microstates, which are not energetically favorable at low temperatures, can now be achieved. An interesting feature occurs along the continuous transition at some blocking temperature $T_B$, which is dependent on the acceptance rate, for which the system passes from a blocked state with a coercive force different from zero to a superparamagnetic one. At this last stage, the high fluctuations in the orientations of the magnetic moments at low fields, make that ${\mathsf{\Gamma }}_{\theta }$ exhibits a peak and concomitantly $\delta \theta$ reaches the maximum value (180°), which makes difficult, under the temperature and field conditions, to equilibrate the desired value of ${\mathsf{\Gamma }}_{\theta }$. The appearance of the peak for ${\mathsf{\Gamma }}_{\theta }$, above the target value ${\mathsf{\Gamma }}_0$, is due to the increment in the new movement acceptance because of those new microstates achieved during the self-regulation process of $\delta \theta$. Blocked states require comparatively small $\delta \theta$ values contrary to those of the superparamagnetic state (see the comparison between continuous and dashed lines in Figure 6e).

More specifically, for low fields close to zero, the orientations energetically favorable are those dictated by the easy anisotropy axes, which are doubly degenerated. Thus, thermal fluctuations are the ones responsible for the moments to alternate not only along such directions but also in between, giving rise to the excess of acceptance rate observed. In consequence an average magnetization close to zero is obtained.

In contrast to the low-field scenario, at high fields (positive or negative) the most likely and privileged orientations are those satisfying the alignment criterion between the magnetic moments and the applied field. Thus, orientations energetically not favorable, although thermally probable, represent a smaller population than those corresponding to zero field. This is the reason an excess in the acceptance rate is not observed. Additionally, we want to stress that our results also show that the superparamagnetic state is achieved at different blocking temperatures depending on ${\mathsf{\Gamma }}_{\theta }$. This fact leads us to conclude that the acceptance rate must be related to the measurement time ${\tau }_m$ involved in the following expression for the blocking temperature (see Section 2.2):

$$T_B=\frac{K_{eff}\mathsf{\Omega }}{k_{B}ln({\tau }_m/{\tau }_0)}.$$

(6)

To validate the above reasoning, Figure 7 shows the $M(H)$ curves for ${\mathsf{\Gamma }}_{\theta }=50\%$ and for some selected temperatures. As observed, some superparamagnetic states are possible to reproduce with constant acceptance rate, i.e., sampling of the phase space occurs at constant speed, except for the one at the highest temperature (400 K). On this basis we can point out that when temperature is high enough the Boltzmann distribution makes any orientation to null fields highly probable, and the acceptance rate increases. If temperature increases indefinitely, all the microstates become equiprobable for any applied field, and the acceptance rate is expected to increase up to $100\%$. Such a limit case is inferred from the Boltzmann probability distribution $P(E)\propto exp(-E/k_{B}T)$ for $T\to \infty$.

4. Conclusions

In this work, we have implemented a novel algorithm, which allows reproducing both the blocked and superparamagnetic states of a system of independent magnetic nanoparticles with uniaxial magneto-crystalline anisotropy randomly distributed. The method presented is based on the Markov chain Monte Carlo method with Metropolis–Hastings algorithm for which the magnetic moments movements are proposed to be accepted at a constant rate ${\mathsf{\Gamma }}_{\theta }$ as phase space is sampled. Consequently, the aperture of the rotations in the updates of the magnetic moments must be self-regulated. Isotherms of $M(H)$ curves show that a constant acceptance rate makes cone aperture of the rotations of the magnetic moments must lie below certain upper bounds. The amplitude of such an aperture is the responsible one for the occurrence of either blocked or superparamagnetic states.

For high values of ${\mathsf{\Gamma }}_{\theta }$, much more microstates must be accepted, so the upper bound for $\delta \theta$ must decrease to satisfy the constant acceptance rate condition. In this case, exploration of the phase space is slow, and it takes time for the system to find states of relaxation. In contrast, for small values of ${\mathsf{\Gamma }}_{\theta }$, much more microstates are rejected, so the upper bound for $\delta \theta$ must increase. In this case, exploration of the phase space is faster, and the system relaxes more easily. Concomitantly, temperature plays a key role in these processes since it helps to make more likely energetically unfavorable events. This causes an extra excess in the acceptance rate and the cone aperture must be readjusted to equilibrate such an unbalance. Additionally, our results allow also to show, from the set of isotherms in the $M(H)$ curves, that the election of a predefined acceptance rate can give rise to different blocking temperatures. This fact leads us to conclude that the acceptance rate must be related to the measurement time.

Finally, a ${\mathsf{\Gamma }}_{\theta }$ value of $10\%$ implies that most of the movements of the magnetic moments are rejected so the exploration of the phase space to find representative microstates is not efficient. In other words, importance sampling is incomplete to guarantee reliable averages of observables. For this reason, we do not recommend using such small values of ${\mathsf{\Gamma }}_{\theta }$.