2. Scaling Rules
A typical MC Metropolis algorithm refers to single spins (magnetic moments) placed in the so-called nodes with defined distances r. This well-known algorithm lies in a random walk over the system nodes and change of the spin direction with the acceptance probability equal to unity if it decreases the energy of the system. If the energy increases, the spin-flip can be accepted with the Metropolis probability exp(−ΔE/kBT) where ΔE = Enew − Eold is the difference in energy between the new and old spin configurations and kBT is the thermal energy.
Our approach involves the magnetic energy of the spin system (further called system energy) calculated in the frame of the continuous 3D Heisenberg model [23]:(1) E=−∑i,jJijSi⋅Sj−∑iKi(S^i⋅n^i)2−gμBμ0∑iHi⋅Si+D∑i,jSi⋅Sj−3(Si⋅eij)(Sj⋅eij)rij3 where Jij is the exchange parameter, Si is the spin vector on site i, Ki is the anisotropy constant (per site), ni^ is the versor of the easy magnetization axis, Si^ is the versor of the spin, 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.
In order to effectively simulate magnetization processes of ferromagnetic materials, it is necessary to use the disorder-based cluster MC method that we also consider appropriate for multi-phase magnetic systems. This approach is based on the Wolff clusterization technique [17] applied to a spin continuous system (described in detail in [22]). Here we present only a general idea based on modification of the so-called adding probability (adding a spin to a cluster) by an additional factor attributed to a local configuration (information) entropy of the selected system’s property (the magnetic anisotropy in our case). Finally, the adding probability takes the form (2) Pijadd=(1−exp(−EijcouplingkBT))exp(−αSiloc) where Eijcoupling is the direct exchange coupling energy between the spins attributed to nodes i and j, Siloc is the local configuration entropy of anisotropy (calculated in the defined sphere around the i-th node) and α is the factor responsible for strengthening and weakening of the entropy impact on the adding probability.
It turns out that the entropy-based modification of adding probability supports and significantly streamlines simulations of magnetization processes in multiphase magnetic systems [24,25]. In order to make a profit from this method, as a tool for designing magnetic materials, it is necessary to analyze relatively large-scale systems, which has a key meaning for both the qualitative and quantitative analyses. Here, we propose the scaling rules, which enable the modeling of large-scale mesoscopic systems, with an acceptable level of computing resource consumption.
Let us consider an initial 3D system of nodes, each occupied by a single spin represented by the vector S, spaced apart by r (refer to Figure 1a). Our re-scaled system is now composed of a set of super-cells characterized by spin vector S’ and distances between them r’ (Figure 1b). The r’/r ratio defines an enlargement factor n which has a key meaning for further considerations. Each super-cell reflects the magnetic moment S’ of n3 spins located in the box of dimensions nr × nr × nr.
In this way, the number of nodes (representing the system) can be reduced by the n3 ratio. In order to carry out MC simulations, it is necessary to keep (i) an equivalence between energies of the original and re-scaled systems when the direction of spin S or S’ changes and (ii) detailed thermodynamic balance during the iteration steps.
The first condition can be fulfilled by analyzing Equation (1) and taking into account that the acceptance and adding probabilities include a change of the system energy related to the spin’s direction changes. Generally, we assume that the “spin” S’, attributed to the super-cell in the re-scaled system (3D case and ferromagnetic coupling), is equal to n3S. In addition, all linear distances between the nodes are re-scaled by taking r’ = nr. The energy change calculated for the re-scaled system can be described using similar relations to the single-spin system but with S’, r’ and some apparent quantities like J’, K’ and D’.
For the change of the exchange energy, the apparent exchange integral parameter Jij’ (energy difference between parallel and anti-parallel super-cells) should be defined in the following way:(3) ΔEij′=−2Jij′S1′S2′=−2Jij′S1S2n6=−2JijS1S2n2
The last term expresses exchange coupling between the neighboring super-cells and includes only the spins on the common wall, as shown in Figure 1b. From this equation, one can determine (4) J′=Jn−4
Equation (1) contains the so-called dipolar energy term, for which the scaling rules should account for the dependence on the relative position between interacting spins. For a simplification of this long-range energy term, we found that the value of the D parameter can remain unchanged for the re-scaled system. The conducted (however, not shown here) comparison between energies calculated for the group of single spins S and the corresponding magnetic super-cells (with re-calculated spins S’ and distances r’), reveals a small error of order 1%. Apart from that, for hard magnetic phases, these kinds of interactions have, generally, marginal meaning due to the high value of magnetic anisotropy energy. Nevertheless, for soft magnetic phases, the error must be taken into consideration.
In addition, the anisotropy constant should be re-scaled taking into account the n3 and n2 factors for anisotropy in the volume and on the surface, respectively. Table 1 summarizes the proposed scaling rules for all parameters appearing in the expression of system energy Equation (1).
The proposed scaling rules should also ensure a detailed thermodynamic balance, similar to the single-spin system. In the MC type of simulation, such a balance is realized through the proper choice of the adding probability, which is a kind of competition between energy change (caused by MC step) and thermal energy kBT. For the re-scaled system, the MC step means a collective change of all spins in a chosen super-cell. It is worth writing an explicit ΔE’ with Equation (1) using the scaling rules (5) ΔEi′=−∑jJ′ijΔS′i⋅S′j−K′i(ΔSi′^⋅ni^)2−gμBμ0Hi⋅ΔS′i+D∑jΔS′i⋅S′j−3(ΔS′i⋅eij)(S′j⋅eij)r′ij3=−n2∑jJijΔSi⋅Sj−n3Ki(ΔSi^⋅ni^)2−n3gμBμ0Hi⋅ΔSi+n3D∑jΔSi⋅Sj−3(ΔSi⋅eij)(Sj⋅eij)rij3
Consequently, (6) ΔEi′=n3Δei where (7) Δei=−n−1∑jJijΔSi⋅Sj−Ki(ΔSi^⋅ni^)2−gμBμ0Hi⋅ΔSi+D∑jΔSi⋅Sj−3(ΔSi⋅eij)(Sj⋅eij)rij3 is the energy change of, let us say, a reduced “e-system”, similar to the original one but with the exchange energy divided by the n parameter. It is interesting that the energy change of the super-cell ΔE’ is equivalent to the energy change of the single spin in n3 e-systems simultaneously. Through this conclusion, the re-scaled system can be simulated with the thermodynamic balance determined for the reduced one, which is more effective, taking into account the fact that exp(−ΔE’/kBT) << exp(−Δe/kBT). It should be emphasized that the most important assumption is the coherent rotation of all spins in the super-cells. This causes a restriction in the temperature range, which has to be significantly lower than the Curie point. Therefore, the value kBT can be considered a kind of “annealing temperature”, enabling effective relaxation of the system and leading to minimization of its free energy.
In summary, the algorithm used and modus operandi are briefly described in the following points:
Step 1. Choose the 3D system size, fitting it to computing resources.
Step 2. Generate a structure of the system and assign the magnetic quantities (Jij, Si, Ki and distance between the nodes rij) to each i-th node. In this stage, the nodes should be considered atoms with a localized magnetic moment and the parameters should reflect their exchange coupling, magnetic moment and anisotropy.
Step 3. Define the n parameter as the re-scaling factor.
Step 4. Assign to each i-th node the re-scaled magnetic quantities (Jij’, Si’, Ki’ and distance between the nodes rij’) using the relations listed in Table 1.
Step 5. Define temperature T (or better kBT in relation to Jij) and external magnetic field H.
Step 6. Relax the system by carrying out Nrelax simulation steps (see below).
Step 7. Calculate average values (in our case, average magnetization) by carrying out Navr simulation steps.
Step 8. Change the field H and go to 6.
The simulation step is based on the cluster MC method, described in detail in [22]. Figure 2 schematically shows the single iteration step including the cluster building procedure as the modification of the Metropolis algorithm.
Focusing on the re-scaling rules, the most important aspect is the determination of the acceptance probability, which depends on the competition between the system energy change (calculated using Equation (5)) and the thermal energy (kBT). As was demonstrated (Equations (5)–(7)), this ratio can be calculated in the following two ways:Calculate the energy change (ΔE’) using the re-scaled quantities and take exp(−n−3ΔE’/kBT) as the acceptance probability.
Calculate the energy change (Δe) using the unscaled distances rij, spins Si, anisotropy constants Ki and Jij/n as the exchange integral parameter, taking exp(−Δe/kBT) as the acceptance probability.
The question may arise of how a limit of the system enlargement can be defined by the n parameter. As the super-cell should be magnetically uniform, its size should not exceed a critical diameter for which magnetization reversal takes place by coherent rotation. More precisely, based on the system parameters (J, S) one can determine a coherence radius Rcoh which is the maximum size of a uniformly magnetized particle [26]. For various magnetic materials, it ranges from more than 10 nm to even more than 1 µm for magnetically soft and hard magnets, respectively [27].
An important aspect of the proposed approach is a direct relation between low-level calculations and the finite element Monte Carlo magnetic simulations. The quantities describing the re-scaled system, e.g., spin vector, exchange integral parameter and magnetic anisotropy attributed to the super-cell, are derived from the unscaled one. In a situation when the nodes in the unscaled system represent atoms in a crystal structure, the mentioned quantities can be determined using, for example, density functional theory (DFT) based calculations.