A guide of the parameters used in the algorithm.

## Abstract

The chapter refers to a modification of the so-called adding probability used in cluster Monte Carlo algorithms. The modification is based on the fact that in real systems, different properties can influence its clusterization. Finally, an additional factor related to property disorder was introduced into the adding probability, which leads to more effective free energy minimization during MC iteration. As a measure of the disorder, we proposed to use a local information entropy. The proposed approach was tested and compared with the classical methods, showing its high efficiency in simulations of multiphase magnetic systems where magnetic anisotropy was used as the property influencing the system clusterization.

### Keywords

- Monte Carlo simulations
- cluster Monte Carlo methods
- magnetic simulations
- entropy methods
- ultra-high coercive alloys

## 1. Introduction

The problem of simulation magnetization processes of multiphase magnetic materials seems to be important, regarding the tendency for applications of high-efficient permanent magnets with reduced or without rare earth elements. Recently, we reported unique hard magnetic properties of Tb-Fe-B-Nb bulk alloys, i.e., coercive field over 7 T at room temperature, attributed to a specific microstructure of dendrite-like Tb_{2}Fe_{14}B grains [1, 2]. In this system, Tb and Fe magnetic moments are coupled antiferromagnetically, which is responsible for relatively low magnetic remanence (μ_{0}*M*_{r} ≈ 0.3 T) and in consequence |*BH*|_{max} (about 13 kJ/m^{3}). However, the Fe-Nb-B-Tb bulk alloys can be considered as a material with extremely high resistance to the external magnetic field and can be a source of magnetic anisotropy in powders as well as bulk spring-exchange composites containing magnetically soft and ultra-high coercive phases. For this reason, the ability of simulating of such systems is very useful in the process of examining and designing of spring-exchange composites in the pre-lab phase.

Monte Carlo (MC) algorithms are now widely used to clarify various physical phenomena, as well as to investigate their potential application in modern technology. Among many simulation methods, the Metropolis MC (MMC) approach [3, 4] is especially attractive in statistical physics for the determination of system physical quantities in thermodynamic equilibrium. The MMC algorithm realizes an ergodic stochastic process, ensuring the fulfillment of the detailed balance condition. One of the bright examples of the application of the MMC algorithm is the Ising model of spins located on the nodes of some lattice. Indeed, using this method one can study a course of magnetic ordering and its dependence of temperature and details of interactions between the spins.

The MMC method utilizes the single-spin-flip procedure to change the spin configuration; however, in many cases (e.g., simulations of magnetization processes) a more effective algorithm is needed. The simplest approach relies on the generation of a cluster of uniformly oriented spins and their subsequent flip to reach new state of the system. The main question is how to determine the cluster and how to establish a rule of its acceptance, simultaneously satisfying the detailed balance condition. Classical approaches, based on the Kasteleyn-Fortuin theorem [5, 6], were proposed by Swendsen and Wang (SW) [7] as well as by Wolff [8] who assumed a specific cluster-building procedure controlled by the so-called adding probability. It is known that the cluster Monte Carlo methods (CMC) are very efficient in the analysis of critical phenomena, e.g., transformation from ferromagnetic to paramagnetic phase [9, 10]. In contrast, their application for studying magnetization processes of systems far below the Curie point produces artificial results, which can be demonstrated for the systems containing magnetically different phases (e.g., hard and soft) as well as geometrical irregularities. The cluster-building algorithms implemented within SW and Wolf approaches are steered by the exchange interactions and the system temperature, but they are not sensitive to other features, potentially affecting the clusterization of spins.

In order to broaden the applications of the CMC methods to simulations of real magnetic composites, we proposed a new method based on some modification of the SW/Wolff adding probability and a particular Metropolis-like algorithm, ensuring the principle of detailed balance [11]. The idea is based on the fact that some kind of regions of the system, characterized by a local disorder of selected system property, constitutes natural barriers for the extension of clusters. In the case of magnetic multiphase composites, spatial distribution of the magnetic anisotropy can be considered as the property affecting the cluster formation.

In the chapter, the disorder-based CMC algorithm is introduced and discussed in a context of classical CMC methods. We show that the new simulation procedure is efficient leading to physically reliable results, especially for multiphase magnetic composites.

## 2. Local disorder-based CMC method

Let the considered physical system be characterized by a discrete spectrum of microscopic energy states (microstates) labeled by *α* is proportional to the Boltzmann factor: *F* can be calculated using the following formula [12, 13]:

where ^{100} states—the number that makes the summation occurring in Eq. (1) impossible.

In order to estimate the average

where probabilities *α* to *β* state) is independent of the states preceding *α*. Moreover, It is also assumed that the detailed balance condition,

The transition probability

As an example, let us consider Ising model of

where *single-spin-flip algorithm*) [12, 13]. This procedure results in uniform distribution of the selection probabilities, i.e., each new state participates in simulations with probability *critical slowing down*) of the Metropolis algorithm is the reason for the difficulties in the generation of statistically independent spin configurations that are needed for the calculation of the estimator

The solution of the critical slowing down problem was proposed by Swendsen-Wang [7] and later by Wolff [8]. The approach developed by Wolf (*cluster-flipping algorithm*) based on the generation of the uniformly oriented spin cluster and its subsequent flipping [12]. In contrast to single-spin-flip algorithm, this procedure easily destroys domains of correlated spins and allows the system to walk through the configuration space. The Wolff algorithm is recognized to be more effective than the Swendsen-Wang one [16, 17, 18]. The selection of a spin cluster starts from randomly chosen spin to which the neighbors occupying the same spin state are added with the probability

Despite the great achievements, the Wolff model is not able to correctly simulate magnetic phenomena that occur far below the critical point in real magnetic materials, including those that are composed of various magnetic phases as well as those containing geometrical irregularities. To be more precise, one can consider, as an example, remagnetization of the system build of two magnetically hard and soft ferromagnetic spheres coupled by a narrow bridge. Let us assume that in the initial state, the magnetization of the system is collinear with the direction of the external magnetic field and then the magnetic field is switched in the opposite direction. What happens to the system is that the magnetization of the soft sphere will follow the change of the magnetic, and then the similar behavior of the hard sphere is expected. Unfortunately, for systems with strong spin-spin coupling, the two-step behavior of the considered system cannot be modeled using Wolff clusterization algorithm. Indeed, every attempt to build a cluster within the two-sphere magnetic system results in that the hard and soft spheres belong to the same cluster independently on magnetic anisotropy and geometry of the system. So, even if the two-step remagnetization process is energetically preferred, the simulated magnetization curve consists of one step related to the common spin rotation. Taking into account the problems encountered during modeling the remagnetization of magnetically inhomogeneous system, we propose a modification of the Wolff algorithm that relies on an assumption that some regions of the system, characterized by a disorder of selected system property, can serve barriers for an extension of magnetic clusters.

We began from the introduction a distribution of system property

Using the Stirling formula, the previous equation takes the form

where *local information entropy* [19]:

So the number of “microstates” corresponding to the

The probability of particular local distribution

Now we can introduce the perception of *local disorder* of the

which gives

We assume that local disorder of

The above equation is the clue of our method which expresses the decrease of the classical adding probability by the disorder-based factor. A set of property, for which the disorder is determined (by the local information entropy), should be chosen dependently on a specific application. The *α* parameter expresses some enhancement or weakness of the influence of the entropy on the system clusterization. It should be defined regarding the specific problem (see next paragraphs).

## 3. Simulation procedure

The main simulation procedure consists of a series of Monte Carlo steps, and it is based on the classical Metropolis algorithm applied to the so-called spin continuous approach. However, the main difference lies in the cluster-building procedure which is executed with a small probability *Pcl*, as shown in Figure 1.

First of all, a random node *i* of the system is chosen, and then it is decided, regarding the probability *Pcl*, whether to analyze the selected node or build a cluster, starting from the node as a seed. In the first case, the algorithm goes to the typical Metropolis procedure and spin continuous method, i.e., the spin direction is randomly modified by angle ±*θ*, and then the energy difference ∆E between the new and the old configuration is calculated. The energy of the system is computed in the frame of the 3D Heisenberg model:

where *Jij* is the exchange parameter, *Si* is the spin vector on site *i*, *Ki* is the anisotropy constant (per site), *ni* is the easy magnetization axis, *g* is the Lande factor, μ_{B} is the Bohr magneton, μ_{0} is the vacuum permeability, *Hi* is the magnetic field on site *i*, *D* is the dipolar constant, *eij* is the directional versor between the *i*th and *j*th nodes, and *rij* is the distance between the *i*th and *j*th nodes. The change is accepted with the Metropolis acceptance probability. In the second case, the algorithm is similar; however, it is necessary to find a cluster around the selected node and, if it exists, carry out a coherent rotation of all spins belonging to the cluster. The procedure of the cluster building is the key and important point. We used the Wolff algorithm but with the modified adding probability

where *E*_{coupling} is the exchange interaction energy between the neighboring spins.

The local information entropy is computed in the defined sphere around each node. In general, the choice of properties (used in the entropy calculation) will depend on the problem being considered; however, for magnetic systems, the natural limit of cluster growth is the change in the value and direction of magnetic anisotropy *K*. Therefore, an optimal feature is a set of three components [*Kx*, *Ky*, *Kz*], whereby the nonmagnetic nodes do not participate in the entropy calculations. In addition, at the beginning of each cluster-building procedure, the *α* coefficient is drawn that will weaken (*α* < 1) or strengthen (*α* > 1) the influence of a local property disorder on the system clusterization. Moreover, it is recommended that from time to time (typically about 20% cases), it completely ignores the impact of disorder and builds a cluster based on the standard Wolff method.

The presented algorithm gives some freedom of the cluster building for which the thermodynamic balance is fulfilled. Indeed, even if the cluster rotation slightly disturbs the balance, the remaining single-spin MCM iterations restore it again. The only condition is that the *P*_{cl} and *θ* values are relatively low (it can be determined experimentally for cases when the results do not depend on the parameters). Finally, one MC step consists of *N*_{iter} iterations presented in Figure 1, and the *N*_{step} steps are taken to obtain the magnetization of the system <*Sz*> (the average spin in the external field direction).

The most important parameters of the algorithm, their possible values as well as our recommendations from a practical point of view are summarized in Table 1.

Parameter | Description | Recommendation and notes |
---|---|---|

θ | The angle at which the direction of the spin or cluster can be changed in a single iteration | Note, it affects the energy changes in a single step. A large θ value will freeze the system, while the small will make it very loose as well as more iterations will be needed to obtain the same change |

P_{cl} | Probability of cluster analysis instead of spin in a single iteration | If many independent clusters in the system are expected, then the value of the parameter should be raised to give everyone a chance to be analyzed. However, high values may destroy the thermodynamic equilibrium |

S_{range} | The range of local information entropy | The range of local entropy should be selected depending on the geometry of the system. The larger range increases the chance to separation between the magnetic grains and a small common area |

N_{iter} | The number of iterations contained in one MC step | Typically, we suggest N_{iter} = 3 N where N is the number of all spins in the system |

N_{step} | The number of MC steps used to calculate the average spin of the system <S_{Z}> | We propose N_{step} = 400 |

α | Modifier of the impact of local information entropy | For α < 1 and for α > 1, the S_{loc} impact will be depressed and strengthened, respectively. The best approach is to draw an α value before each cluster search |

P_{Wolff} | The probability of searching for a cluster based only on the Wolff method | In some cases, you can build a cluster ignoring the impact of local entropy. Typically, we suggest P_{Wolff} = 0.2 |

The algorithm can be, in some aspects, paralleled [20, 21, 22, 23]. The main limitations of the parallelization process are due to three reasons. First of all, each step in the Monte Carlo procedure should be based on the system modified in the previous step. In particular, the decision to accept the new spin direction depends on the current direction of the neighboring spins due to the exchange energy. Consequently, a situation in which different threads are testing a new configuration for the two neighboring spins at the same time should be refused. Similarly, simultaneous analysis of the whole cluster and individual spins inside it (as well as spins interacting with it) is not allowed. The second thing to keep in mind is to ensure that each state of the system can be selected with the same probability. This means that none of the actions can interfere with the probability of choosing a spin for analysis. For example, the spin cannot be temporarily omitted in the draw as well as the analyzed cluster must always be found in the real time and cannot be taken from a fixed database. Finally, it should be taken into account that a change of the spin direction inside a cluster, which is constructed based on the

Despite all the above restrictions, there are three time-consuming problems in the algorithm which can be parallelized: the main MC loop, the calculation of cluster energy as well as the cluster-building procedure around the chosen spin. Our experience shows that the parallelization of the algorithm accelerates the computation time more than 10 times (depending on the probability of cluster analysis and the system size) and maintains the correctness of results. Figure 2 shows a comparison of computation time (needed for 400 MC steps) for different *P*_{cl}, system size, and its magnetic properties. On this graph, the independent quantity is a number of threads applied for the simulations (HP ProLiant DL580 G7, 4x Intel Xeon 8C X7560 2.27 GHz).

Generally, this aspect of calculations is very complex. Nevertheless, one can conclude that the parallelization is more effective for the systems with small number of the cluster-building trials.

## 4. Application of the method to multiphase magnetic systems

In order to show an efficiency of the disorder-based cluster MC algorithm and its comparison with the classical Wolff method, the two systems containing magnetically hard and soft phases were analyzed. The 3D system space consists of 50 × 50 × 50 (1,25,000) nodes; the spins are arranged in shapes of joined spheres with hard and soft magnetic properties. The bonding between the spheres, further called “bridges,” is different for the two cases, and it equals 1 or 7 nodes, receptively. Figure 3 depicts an example of the system with one-bridge coupling between the hard and soft spheres. The parameters of the magnetic phases and simulation procedure are listed in Table 2.

Parameter | Soft phase | Hard phase |
---|---|---|

Exchange coupling J | 1.5e–2 eV | 1.5e–2 eV |

Anisotropy constant K | 0 | 1e–3 eV |

Simulation procedure | ||

Dipolar constant D | 1.8e–7 eV | |

Thermal energy k_{B}T | 1e–4, 1 e–5, 1e–6 eV | |

α parameter | Random value from 0 to 25 | |

Cluster analysis probability P_{cl} | 10^{−3} | |

Range of the entropy S_{range} | ±3 nodes | |

Spin change angle θ | π/100 | |

System size (in one direction) n | 50 | |

Number of iteration in one MC step | 3n^{3} = 3,75,000 | |

Number of MC steps for magnetization averaging | 400 |

The difference between classical and disorder-based approaches lies in the definition of adding probability. It can be demonstrated by the cluster-building successes when the procedure starts from the same cluster “seed.” Figures 4 and 5 show the comparison for the one- and seven-bridge systems, respectively. The cluster seed is placed in the center, and the white color is attributed to 100 cluster-building successes per 100 trials.

As shown, using the classical Wolff algorithm, all spins in the system belong to the same cluster independently on their magnetic properties. This means that during simulations, a possibility of closing the clusters inside the hard or soft phase is not available. The modification of adding probability causes that the cluster-building procedure can distinguish the hard and soft phase with relatively high probability. However, with lower probability, the cluster contains whole spins in the system (like in the Wolff algorithm) which is related to a random value of the α parameter. The modification of adding probability causes that the configurations necessary for modeling of real magnetization processes are available during MC iterations. In other words, if magnetization processes require separate behavior of the hard and soft phase, the simulation procedure will test such a possibility.

As a final test, the so-called reverse magnetization curves were simulated using the classical as well as modified adding probability. Initially, the system was saturated in the field direction (all spins are directed up), and then the field was switched off. During calculations the field was increased in the opposite direction. Magnetization, determined as average spin value in the field direction, as a function of the external magnetic field for all examined cases is shown in Figure 6.

The difference between *P*_{add} and _{B}*T* above 1e–6 eV, this process consists of a subsequent change of the “soft” spins and next, coherent flip of the “hard” spins. Such a scenario is fully confirmed by the spin configurations depicting the reverse magnetization process for the seven-bridge system at k_{B}*T* = 1e–4 eV (see Figure 7).

Regarding the fact that in some cases also for *Padd* the two-step behavior occurs, the main question is which curve is physically correct. For this reason, it is worth to analyze the curves in fields when the soft phase changed spin direction using

It is clear that the introduced modification of adding probability results in more effective finding of the free energy minimum and, therefore, produces physically reliable results which allows modeling multiphase magnetic systems.

## 5. Conclusions

As it was shown, in some conditions (low temperature and/or strong exchange coupling) the classical CMC algorithm can produce incorrect results in the application of multiphase magnetic systems. The problem lies in the adding probability used in the cluster-building procedures. This problem can be overcome by the proposed modification of *Padd* accompanied by the specific Metropolis-like algorithm. The key point is to recognize a property that influences the system clusterization and, next, introduction of the local information entropy of this property into the adding probability.

In summary, one can conclude that the method proposed for modeling of magnetic systems allows sampling spin configurations practically inaccessible by original Wolff clusterization procedure. An interesting future of our method is the acceleration of the energy minimization process.

## Acknowledgments

This work was supported by the National Science Centre in Poland by grant no. 2015/19/B/ST8/02636 (A. Chrobak, G. Ziółkowski). D. Chrobak acknowledges support from the National Science Centre of Poland, grant no. 2016/21/B/ST8/02737.