{'fraction_non_alphanumeric': 0.027822415477203263, 'fraction_numerical': 0.011374736112423985, 'mean_word_length': 5.137690968864823, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Micromagnetic Monte Carlo method with variable magnetization length based on the Landau-Lifshitz-Bloch equation for computation of large-scale thermodynamic equilibrium states

An efficient method for computing thermodynamic equilibrium states at the micromagnetic length scale is introduced, using the Markov chain Monte Carlo method. Trial moves include not only rotations of vectors, but also a change in their magnetization length. The method is parameterized using the longitudinal susceptibility, reproduces the same Maxwell-Boltzmann distribution as the stochastic Landau-Lifshitz-Bloch equation, and is applicable both below and above the Curie temperature. The algorithm is fully parallel, can be executed on graphical processing units, and efficiently includes the long range dipolar interaction. This method is generally useful for computing finite-temperature relaxation states both for uniform and non-uniform temperature profiles, and can be considered as complementary to zero-temperature micromagnetic energy minimization solvers, with comparable computation time. Compared to the dynamic approach it is shown the micromagnetic Monte Carlo method is up to almost 20 times faster. Moreover, unlike quasi-zero temperature approaches which do not take into account the magnetization length distribution and stochasticity, the method is better suited for structures with unbroken symmetry around the applied field axis, granular films, and at higher temperatures and fields. In particular, applications to finite-temperature hysteresis loop modelling, chiral magnetic thin films, granular magnetic media, and artificial spin ices are discussed.


I. Introduction
Monte Carlo computations are widely used in many fields of research, including statistical physics and atomistic spin lattice modelling. For the Ising and Heisenberg Hamiltonian spin lattice models, the Metropolis Monte Carlo [1] algorithm has proved popular.
This allows computation not only of magnetic parameters temperature dependences and phase transition temperatures, but also hysteresis loops [2]. In micromagnetics modelling a common need is computation of relaxed magnetization states. For zero-temperature models this may be accomplished efficiently using energy minimization solvers, such as a steepest descent method [3]. Finite-temperature micromagnetic models include the stochastic Landau-Lifshitz-Gilbert equation (sLLG) [4], as well as the Landau-Lifshitz-Bloch (LLB) equation [5] and its stochastic forms (sLLB) [6,7]. Whilst the sLLG equation is a reasonable approximation at low temperatures, it fails to take into account the magnetization length distribution, which has been shown to play an important role in magnetization reversal [8], including linear and elliptical reversal mechanisms [9], particularly for temperatures approaching the Curie point. Thus, whilst the sLLB equation may be used to compute finite-temperature relaxation states, this dynamic approach to relaxation is very inefficient and requires lengthy computation times.
Alternative methods, such as a Monte Carlo approach, are required to efficiently compute relaxation states of thermodynamic equilibrium, which should ideally be applicable across the entire temperature range, both below and above the Curie temperature; for this reason correct reproduction of the magnetization length probability distribution at the micromagnetic length scale is essential.
A previous work extended the Monte Carlo method to the micromagnetic length scale, the Curie temperature, the distribution of the free energy does not follow a Boltzmann distribution, but a Maxwell-Boltzmann type distribution [7]. Moreover, for the purposes of computing thermodynamic equilibrium states at the micromagnetic length scale the hybrid Monte Carlo approach is inefficient, since an ensemble of conjugate momenta must be evolved over many iterations (typically 100 or more), before a new micromagnetic magnetization configuration may be accepted or rejected.
Here we describe a micromagnetic Monte Carlo (MMC) method based on the Markov Chain approach, thus similar to the Metropolis Monte Carlo method [1] commonly used for atomistic spin lattice simulations, but including trials moves for magnetization length change.
The resulting method not only correctly reproduces the expected Maxwell-Boltzmann distribution for the free energy, but generates new micromagnetic magnetization configurations every iteration. The algorithm has been implemented in Boris [17], and is publically available and open-source [18]. In Section II the algorithm and related theory are described and tested.
In Section III the inclusion of the long range demagnetizing interaction in the parallel MMC implementation is described, and tested by computation of hysteresis loops and comparison with established methods. In Section IV the application of MMC to chiral magnetic films with Dzyaloshinsky-Moriya interaction (DMI) [19,20] is discussed. Finally, in Section V the application of MMC to the study of artificial spin ices is addressed, before summarising possible future developments in the concluding remarks.

II. Micromagnetic Markov Chain Monte Carlo Method
For atomistic spin lattice models, where the spins have fixed length, at thermodynamic equilibrium the internal energy E follows a Boltzmann distribution, Z is the partition function, kB is the Boltzmann constant, and T is the temperature. At the micromagnetic length scale however, where the magnetization length can vary due to thermodynamic averaging over many atomistic spins, the distribution of the free energy follows a Maxwell-Boltzmann type distribution [7], given in Equation (1).
required thermodynamic equilibrium distribution. A sufficient condition for convergence is that of detailed balance: is the transition probability from state A to state B in the Markov chain. Thus we obtain the following ratio: The Markov chain transition probability is given by: (2) 6 Here ) (m f l is the magnetization length probability distribution, with Zl a renormalization factor, given by: The longitudinal term contribution in Equation (7), where V is the micromagnetic computational cellsize volume (e.g. V = h 3 for a cubic cell with cellsize h), is the same which appears in the sLLB or LLB equation, and which gives rise to a longitudinal field. It is parametrized using the relative longitudinal susceptibility given as [5]: Here is the Langevin function, µ is the atomic moment, TC is the Curie temperature, and me = Me / Ms0 is the normalized temperature-dependent equilibrium magnetization length given as [5]: The remaining terms which contribute to the free energy ) (m F in Equation (6), include all the usual micromagnetic terms, and these may be obtained directly from the corresponding energy density terms by multiplying with the computational cellsize volume V. There are a few terms which require special attention however, namely exchange interactions (direct and DMI) and demagnetizing interaction. First we treat the direct exchange interaction, and discuss the remaining terms later. The micromagnetic exchange free energy for a given magnetization vector i is given below, where A is the temperature-dependent exchange stiffness: It is important to note the multiplicative factor of 2 here: the micromagnetic exchange energy density (see e.g. Ref. [17]) expression is derived from the atomistic Heisenberg direct exchange Hamiltonian contribution which includes compensation for double-counting of spins. In order to calculate the free energy contribution of a single magnetization vector however, just as with the atomistic Monte Carlo method, the factor of 2 must be restored.
Single-site energy terms are obtained directly from the corresponding micromagnetic energy density expressions [17], for example the uniaxial magnetocrystalline anisotropy is included as: Here A e is the symmetry axis direction, K1 and K2 are the 2 nd and 4 th order uniaxial anisotropy constants respectively. Cubic magnetocrystalline anisotropy is included as: Here K1 and K2 are the 4 th and 6 th order cubic anisotropy constants respectively, i, i, and i are direction cosines of the magnetization vector.
Thus the MMC method consists of the following steps: 1) for a vector i perform a rotation trial move in a cone of given solid angle. 2) For the same vector i perform a magnetization length change trial move. 3) Compute the total free energy change, including the longitudinal free energy change using the energy term appearing in the exponent of Equation (7). 4) Accept the compound trial move with probability given in Equation (5). This procedure is done exactly once for each micromagnetic vector, which completes an MMC iteration. The algorithm is fully parallel, implemented with the red-black checkerboard decomposition scheme discussed previously for the atomistic Monte Carlo algorithm [21,22], and is thus suitable for computations on graphical processing units (GPUs). The final question is how the solid angle for the rotation trial moves should be chosen, and how the magnetization length should be changed. With time-quantized Monte Carlo (TQMC) [23,24], the cone solid angle is set such that real-time processes may be reproduced, for example the effect of field sweep rate on hysteresis loops [25], or thermally activated grain reversal time [26]. This may also be possible with MMC, however this topic is outside the scope of the current work. Instead, the cone solid angle is adaptively set such that the acceptance rate is kept at an optimal value of ~0.5 [22], which allows rapid thermalization across the entire applicable temperature range. In practice this needs to be capped to a maximum value, and we use   0.03 which is reached close to TC, and also for T > TC this latter constant value is used. It should be noted that since the algorithm generates a magnetization length distribution governed by Equation (7), the probability of generating a vector with |M| > MS0 is practically zero.
We now verify the MMC method correctly reproduces the target probability obtained by fitting numerically computed distributions, as shown in Figure 1, with Equation . This shows the MMC method correctly reproduces the magnetization length probability distribution across the entire simulated temperature range. Finally, we verify the transverse distribution in Equation (6) is reproduced correctly.
For this, the angular deviation from mean direction probability distribution is computed. Thus for each micromagnetic magnetization configuration the mean direction is obtained, and for each vector in the ensemble the polar angle  to this mean direction is found. Using 100 bins for , between the minimum and maximum  values, a histogram is obtained. This is then averaged over >50,000 ensembles. Example results are shown in Figure 1(d)-(f), both for MMC and sLLB, showing excellent agreement. Note, the histograms thus obtained are weighted by the bin solid angle differential surface area on the unit sphere, hence the computed probability tends to zero as  tends to zero or to  radians.

III. Hysteresis Loop Modelling
One important application for MMC is finite-temperature hysteresis loop modelling.
Here we concentrate on the static hysteresis limit; using kinetic Monte Carlo [28,29] it is possible to simulate the effect of field sweep rate. Alternatively TQMC may be employed, and time quantization may also be extended to MMC, however this is left for a future work. For computation of hysteresis loops, in general it is essential the demagnetizing interaction is included, allowing for shape anisotropy effects. Such long range interactions are notoriously difficult to include in a parallel Monte Carlo algorithm. Since each spin interacts with all other spins, an exact computation of the interaction energy is not compatible with domain decomposition methods. For the atomistic spin lattice case the dipole-dipole interaction was previously approximately included using a complex stream processing domain decomposition method [30]. Here we discuss a simpler alternative approach, which only requires a single demagnetizing field computation per MMC iteration, and may be achieved using the usual efficient FFT-based convolution method [17,31]. The demagnetizing field may be obtained using the discrete convolution sum shown below, where N is the demagnetizing tensor: The free energy contribution of magnetization vector i is given by Fd,i, where as for the exchange interaction the value is compensated for the factor of 1/2 arising from double counting of magnetization vectors in the usual demagnetizing energy density formula. At the start of an MMC iteration the demagnetizing field is fully computed, and for each trial move the energy change Fd,i is obtained from Equation (14) as Once a trial move is accepted, the demagnetizing field is not immediately updated. Thus whilst this approach is also inevitably an approximation, it may be made more accurate by having more than one demagnetizing field update per MMC iteration; for example since a red-black checkerboard decomposition scheme is used, the demagnetizing field can be updated before each red and black parallel passes. In turn, the red and black passes can be further decomposed to include more demagnetizing field updates. This approach was numerically tested extensively, also against a serial algorithm implementation where Fd,i was computed exactly.
Results comparing the parallel and serial MMC implementations are shown in Appendix B, with no statistical difference found between them. This method works since locally the demagnetizing interaction is relatively weak, and the change in demagnetizing field from one iteration to another is small. Thus, in thermodynamic equilibrium, locally the trial move acceptance probability is largely driven by the exchange interaction, with |Fex,i| being over 1 to 2 orders of magnitude larger than |Fd,i|, even for ultrathin films with perpendicular magnetization (see Appendix B). The effect of the demagnetizing field becomes apparent only over many iterations and on a large spatial scale, appearing as a bias in the acceptance probability. First, we take a simple example to verify the shape anisotropy effect is correctly reproduced, by comparing the MMC method with the sLLB equation solution. Easy and hard axis hysteresis loops in a 400 × 200 × 5 nm Ni80Fe20 ellipse, with 5 nm cubic discretization cellsize, are computed at 300 K and 600 K, with results in Figure 3 showing excellent agreement. Note, for all the hysteresis loops shown in this work only the increasing field sweep has been simulated, with the decreasing field sweep completed through inverse symmetry. For sLLB, at each field step the magnetization configuration solution is accepted when the average normalized torque value falls below 10 -4 (smaller values lead to excessive simulation times).
With the MMC method we need to ensure enough iterations are used to obtain correct thermodynamic equilibrium states. From atomistic Monte Carlo simulations of hysteresis loops, it is known the switching field estimation becomes increasingly more accurate as the number of Monte Carlo iterations increases [2]. One could run the MMC algorithm with a fixed number of iterations for each field step. Here a more efficient approach is taken; a chunk with a fixed number of iterations (2000)  As a further test we show the hysteresis loops obtained for a 400 × 200 × 100 nm Ni80Fe20 ellipsoid at 300 K, also discretized using a 5 nm cubic cellsize, computed along the easy, medium, and hard axes. Here MMC is compared with results obtained using the steepest descent (SDesc) energy minimizer [3], where field steps are considered solved when the maximum normalized torque value falls below 10 -5 . This is used as a quasi-zero temperature method, where the material parameters are simply rescaled for the required temperature, but otherwise the SDesc method does not include stochasticity or a magnetization length distribution. A very good agreement is observed for the medium and hard axes, with hysteresis loop behaviour largely dictated by the shape anisotropy effect. A reasonable agreement is also obtained for the easy axis hysteresis loop, however here stochasticity also plays a part, resulting in the magnetization configuration switching to a vortex state close to zero field sooner for MMC compared to SDesc; the zero-field vortex state is shown in the inset to Figure 4  scaling law [32] is used, with l being the anisotropy term order. The simulation space consists of a 800 nm  800 nm  20 nm mesh with periodic boundary conditions in the plane, and cubic discretization cell with 2 nm cellsize (total of 1.6 million cells included in simulation). The grains are generated using three-dimensional Voronoi tessellation with average 10 nm spacing between seed points and with phase separation (i.e. grains interact through the demagnetizing field only). A typical granular structure generated is shown in the inset to Figure 5, with the resulting hysteresis loops obtained from MMC and temperature-rescaled SDesc also shown for T = 300 K. The field was applied along the in-plane horizontal direction (along cubic easy axis). A good agreement is obtained between the two methods, confirming the demagnetizing interaction is correctly evaluated in the implemented MMC method. Granular films are important for magnetic storage media, including heat-assisted magnetic recording [33]. Specialized methods exist for computing hysteresis loops in granular media, including a thermal activation model [34] and micromagnetic kinetic theory [35]. The MMC approach presented here however, is a generalpurpose method for computing thermodynamic equilibrium states; the resulting collection of states are the same as those that would be obtained by integrating the sLLB equation given sufficient simulation time. For relaxation problems the MMC method is over 2 to 3 orders of magnitude faster compared to sLLB. We also note for the granular film problem of Figure 5, the MMC method is slightly faster (typical computation time ~1 h on a modern GPU) than the SDesc method, which requires many iterations to reach the normalized torque convergence threshold of 10 -5 .

IV. Magnetic Thin Films with DMI
The free energy contribution of the interfacial DMI term is given below, which is included in addition to the direct exchange contribution of Equation (10).
Here Hysteresis loops at 300 K, 350 K, and 400 K are shown in Figure 6, where the field is applied perpendicular to the plane. The hysteresis loops are typical of experimental results [37]. A labyrinth domain structure is observed at zero field in all cases, resulting in sheared hysteresis loops. Inspecting the increasing field sweep, for T = 300 K the labyrinth domain structure is suddenly formed at a negative field, nucleated through thermal activation, as indicated in the inset to Figure 6(a). As the field strength is increased the labyrinth domain structure is gradually reduced, eventually forming a small number of skyrmions; further increasing the field results in thermal annihilation of skyrmions, leading to saturation. The magnetization reversal process depends on the sample temperature, not only due to increased effect of stochasticity, but more importantly due to the temperature dependence of the DMI constant and effective anisotropy.
Thus at higher temperatures, instead of nucleating a labyrinth domain structure, first skyrmions are nucleated, which then grown into a labyrinth structure as shown in Figure 6(b),(c). Also at higher temperatures, increasing the field towards saturation results in a dense array of skyrmions formed out of the labyrinth structure, approximately arranged into a hexagonal latticesee inset in Figure 6(b). The thermodynamic equilibrium states arising at each field step for the hysteresis loops in Figure 6 are detailed in the Supplemental Information.  Finally, we also compute the hysteresis loop with a non-uniform temperature profile, in particular a Gaussian temperature profile as in Ref. [38], with 300 nm width, reaching a maximum of 400 K at the centre, from 300 K at the extremities. Such non-uniform temperature profiles are encountered in ultrafast laser-induced skyrmion nucleation studies [37][38][39][40]. Whilst the MMC method evidently does not reproduce the dynamics, it may be used to analyse the states resulting on a long time-scale, such as mixed labyrinth domain and skyrmion states.
Results are shown in Figure 7. In contrast to the simulation with a uniform temperature, where a labyrinth domain structure is nucleated, instead a number of skyrmions are nucleated inside the Gaussian profile at a larger negative field, which then grow into a labyrinth domain structure. On the positive side of the hysteresis loop, first a number of skyrmions are formed from the labyrinth structure inside the Gaussian profile, resulting in a mixed state as noted in a recent experimental study [37].
The MMC method may thus be used to study thermodynamic equilibrium states also for materials with DMI, which includes hysteresis loop modelling, but also relaxed states in magnetic nanostructures, for example as observed in a recent study on Co/Pt multi-layered nano-dots [41]. For such multi-layered structures the demagnetizing field used for the MMC method may also be computed using the multi-layered convolution algorithm [31]. Analysis of skyrmion arrangements and temperature-driven diffusive motion in magnetic nano-structures was previously done using a quasi-particle Monte Carlo model [42]. The MMC method may also be used for such studies to analyse the probability distribution of skyrmion positions, having the advantage of treating the energy terms on an equal footing with micromagnetic approaches, with full spatial resolution. Thus in particular, the MMC method naturally takes into account skyrmion deformations, thermal nucleation and annihilation of skyrmions.

V. Artificial Spin Ices
Artificial spin ices (ASI) is an active and growing field of research; for recent reviews see Refs. [43,44]. Here we briefly point out the applicability of the MMC method to the study of ASI. Monte Carlo methods have been successfully applied to the study of ASI, relying on a point dipole approximation for the magnetic moments of islands comprising ASI arrays [45,46]. On the one hand this allows computation of the blocking temperature [47], TB; annealing with T > TB (normally TB < TC) allows the magnetic islands to settle into a different overall energy state, as the energy barriers for flipping their moment directions are overcome, dictated by magnetic frustration due to the dipolar interaction. Use of the dipole approximation Monte Carlo method also allows study of statistics of the vertex population types (vertex configurations grouped by energy state). The MMC method presented here may also be used, having the advantage of not relying on a dipole approximation for the magnetic islands moments. This allows for example resolving the domain configuration of magnetic islands, which may not be in a single domain state, particularly for larger magnetic islands. Effects due to magnetic islands shape can in principle be taken into account, and also connected ASI arrays may be studied, where the exchange interaction contributes to frustration. Thus the MMC approach is a more general method compared to the point dipole Monte Carlo approach; the computational complexity is greater, however efficient use may be made of GPUs which, as we show below, allows simulations of ASI arrays of dimensions approaching those used in experimental studies. In Figure 8 we show results for a square ASI array, with dimensions of 6 µm × 6 µm × 25 nm. The Ni80Fe20 islands are stadium-shaped and of dimensions 200 × 40 × 25 nm, with lattice spacing of 250 nm. The simulation space was discretized using a 5 nm cubic cellsize, with in-plane periodic boundary conditions, thus resulting in 7.2 million simulation cells. An example region of the simulated ASI array is shown in Figure 8(a). Here we compute the blocking temperature TB, as well as the vertex population statisticsfor a square ASI, as is well known, there are 4 vertex types as given e.g. in Ref. [48]; for brevity the definitions are not repeated here. The simplest approach to computing TB is by analysing field-cooled (FC) and zero-field-cooled (ZFC) curvessee e.g. Ref. [47]. We also use this method with MMC.
Thus, starting from a large temperature near TC, a small field of 10 kA/m, but large enough to switch the magnetic islands along the field direction, is applied -FCand the magnetization length in thermodynamic equilibrium along the applied field direction (horizontal direction in Figure 8) is recorded. The same is repeated in zero field -ZFCand the temperature differential of the FC and ZFC difference curve is plotted in Figure 8(c); the point of minimum is the blocking temperature, obtained as TB = 656 K.
Next, an alternative method of computing TB, which has the advantage of obtaining statistical information on vertex populations, is used. The ASI array is saturated in a large field (100 kA/m) at room temperature, along the horizontal direction, then relaxed at zero field; this starting state, showing all the horizontal islands pointing towards the right (red color) is indicated in Figure 8(a). Next, an annealing temperature is set and thermodynamic equilibrium achieved, before relaxing the ASI array at room temperature again where the vertex types are counted. This is repeated with increasing annealing temperature, and results are shown in Thus the MMC method may also be used to study ASI, including not only square arrays, but also kagome lattices [49] and magnetic quasicrystals [50].

VI. Conclusions
In hysteresis loops, and grain reversal times. Finally, the method introduced here was applied to ferromagnetic materials, however a future work will investigate extensions to two-sublattice models, allowing applications to ferrimagnetic and antiferromagnetic materials [51]. The SDesc and MMC methods generally agree well at low temperatures far below TC, especially when magnetization reversal occurs largely through coherent rotation. An important exception to this is for processes where magnetization reversal is principally driven by thermally activated reversal, and we have given an example in the main text for chiral magnetic thin films. For such problems quasi-zero temperature methods (SDesc) are inadequate. Here we further show also for high temperatures the SDesc method is increasingly inaccurate, necessitating the use of MMC. Figure 9 shows the Ni80Fe20 ellipsoid problem of Figure 4, with hysteresis loops given along the easy axis for 300 K, 500 K, 600 K and 650 K (TC = 680 K). As the temperature increases the MMC and SDesc solutions increasingly diverge. The difference between the solutions at high fields is due to lack of stochasticity in the SDesc method, which becomes increasingly important at higher temperatures, as the average angular deviation from the mean direction increases (see e.g. Figure 1(d)-(f)). Thus the SDesc solutions reach the magnetization saturation value at low fields, failing to reproduce a realistic high field behaviourcomparisons of high field behaviour for the SDesc and MMC solutions are shown in Figure 10. For the MMC solutions however large fields are required at high temperatures to fully saturate the magnetization along the applied field direction, requiring narrowing of the angular deviation from the mean direction probability distribution. It should be noted, for the sLLB equation

Appendix A -Comparison of MMC and quasi-zero temperature solutions at high temperatures and fields
form of Ref. [7] for which the current MMC method reproduces the same thermodynamic equilibrium propertiesthe un-weighted average magnetization does not coincide with the equilibrium magnetization, as noted in Ref. [7] and as may easily be verified numerically.
Conceptually this is problematic, since in atomistic modelling the two values are identical. This problem has been addressed recently [52], with an alternative form of the sLLB equation obtained from the LLB equation in Ref. [5]. It may be possible to extend the MMC method to this alternative sLLB formulation, however this is left for future work.

Appendix B -Free energy contributions to MMC trial moves
The free energy contributions to MMC trial moves are computed a) for the Ni80Fe20 ellipsoid problem of Figure 4, and b) for the ultrathin Co/Pt film (here 1 nm Co thickness) of Figure 6, in order to illustrate the relative importance of the different interactions to the trial move acceptance rate. Results are shown in Figure 11 as a function of temperature, with a) the ellipsoid initialized at zero field with the vortex state shown in the insets to Figure 4, Figure 4, initialized at zero field with the vortex state shown in the insets to Figure   4, and (b) the ultrathin Co/Pt film (here 1 nm Co thickness) of Figure 6.
Finally, we compare the parallel MMC algorithm, where the demagnetizing field is update once per iteration, with the serial MMC version where the demagnetizing field is updated after every accepted move. The latter thus exactly takes into account the long range demagnetizing interaction. Magnetization switching events are computed as a function of number of iterations for the ellipse problem of Figure 3. Starting from the remanence state with magnetization pointing along the -x direction, the coercive field value of H = 12 kA/m is applied, and the magnetization along the field recorded over 10,000 iterations. Individual switching events are shown in Figure 12(a), as well as averages over 500 switching events for both algorithms. The switching probability is strongly influenced by the demagnetizing field, which gives rise to an effective uniaxial anisotropy. The probability of crossing the x axis is plotted in Figure 12(b) as a function of iteration number. These results show there is no significant difference between the two approaches, again justifying the use of Equation (14) for including the demagnetizing interaction in the parallel MMC algorithm.

Appendix C -MMC and sLLB computation times comparison
Comparisons of computation times between the MMC and sLLB methods are shown here. These are given in TABLE I for the ellipse and ellipsoid problems of Figure 3 and Figure   4 respectively For MMC the demagnetizing field is computed once per iteration as explained in the main text. For sLLB with the RK4 method, which consists of 4 sub-steps per time-step iteration, normally the demagnetizing field is computed once per sub-step. In order to improve the computation time a recently proposed method of speeding up explicit evaluation methods using demagnetizing field polynomial extrapolation [53] was used, thus also employing one demagnetizing field evaluation per iteration. As expected, computation times for MMC are significantly smaller, with speedup factors of almost 20 obtained as shown in TABLE I.  (Fig. 3) and ellipsoid (Fig. 4) problems. These are computed on the CPU, respectively GPU as indicated, with speedup factor defined as time for sLLB divided by time for MMC computation completion.

Data Availability
The data that support the findings of this study are available from the corresponding author upon reasonable request, and the micromagnetic Monte Carlo algorithm code is available at Ref. [18].

Supplemental Information
Videos of relaxed magnetization image sequences along the increasing field sweep are shown for each of Figure 6(a), (b), (c) and Figure 7 hysteresis loops.