# Multicanonical Sampling of the Space of States of -Vector Models

###### Abstract

Problems of temperature behavior of specific heat are solved by the entropy simulation method for Ising models on a simple square lattice and a square spin ice (SSI) lattice with nearest neighbor interaction, models of hexagonal lattices with short-range (SR) dipole interaction, as well as with long-range (LR) dipole interaction and free boundary conditions, and models of spin quasilattices with finite interaction radius. It is established that systems of a finite number of Ising spins with LR dipole interaction can have unusual thermodynamic properties characterized by several specific-heat peaks in the absence of an external magnetic field. For a parallel multicanonical sampling method, optimal schemes are found empirically for partitioning the space of states into energy bands for Ising and SSI models, methods of concatenation and renormalization of histograms are discussed, and a flatness criterion of histograms is proposed. It is established that there is no phase transition in a model with nearest neighbor interaction on a hexagonal lattice, while the temperature behavior of specific heat exhibits singularity in the same model, in case of LR interaction. A spin quasilattice is found that exhibits a nonzero value of residual entropy.

## I Introduction

Any equilibrium statistical properties of complex systems consisting of interacting particles can be calculated if the value of energy of each of states is known. The knowledge of the shape of the interaction energy distribution in the space of these states, , i.e., the distribution of the density of states (DOS) , or the so-called energy landscape Wales2003 (), could allow one to precisely calculate thermodynamic averages. However, an exponential increase in the number of configurations and a complex growth of the number of allowed values of energy as a function of makes this problem extremely complicated. The equilibrium thermodynamics of a system of interacting magnetic moments can be analyzed by various approximate numerical methods. One of the most promising and actively developed methods that allows the sampling of the space of states and calculate the is the method, proposed by Wang and Landau PhysRevLett.86.2050 (); Landau2004 (), which is called the WangâLandau (WL) method in western scientific literature. This method allows one to dynamically construct the by conditional random walks over the entire energy landscape by forming the statistical weights for degenerate or quasidegenerate (i.e., differing from each other by infinitely small interaction energy) states.

The wide spectrum of applied and fundamental problems, general function, immense field of applications, elegance, and obvious simplicity of the WL method have been the reason for the widest variety of its applications in statistical physics, biophysics, and other fields starting from spin systems PhysRevLett.86.2050 (); PhysRevE.64.056101 (); doi:10.1142/S0129183102003243 (); ref1 (); PhysRevLett.90.120201 (); PhysRevE.70.066128 (); Brown2005 (); PhysRevE.71.046705 (); PhysRevE.72.056710 (), quantum systems PhysRevE.82.046703 (), atomic clusters PhysRevE.89.013311 (); Calvo2003 (), dipole PhysRevB.72.214203 () and spin PhysRevLett.110.210603 () glasses, liquid crystals PhysRevE.72.036702 (), fluids Desgranges2009 (), XY model Ngo2010 (), the BlumeâCapel model PhysRevE.92.022134 (), the Potts model Caparica2015447 (), biomolecules PhysRevE.90.042715 (), protein folding Rathore2002 (), polymer films TSAI2006 (), as well as in many other fields of science, for example, in solving optimization problems ArgollodeMenezes2003428 (), development of combinatorial number theory 0305-4470-36-24-304 (), and others.

The relatively fast convergence of the WL method makes it excellent for calculating the DOS; however, one should note that there are some questions on the accuracy of calculations, which have been pointed out by many authors (see, for example, PhysRevE.72.025701 ()). Certainly, the WL algorithm is one of the most interesting and unusual achievements of the last decade in Monte Carlo (MC) simulation methods PhysRevE.75.046701 ().

The method is based on the algorithm for calculating the DOS â the relative number of possible states (configurations) for the energy level . Thus, one can calculate any thermodynamic quantities that describe configurations arising during random walks over the space of states, including the free energy, in a wide range of temperatures.

The use of histograms for obtaining significantly increased the level of simulation. This progress made it possible to obtain results for models of complex systems under test and allowed one to make significant progress in the numerical simulation and description of phase transitions. In contrast to the well-known and widely used MC algorithms, such as the Metropolis algorithm Metropolis1953 (), the SwendsenâWang cluster MC simulation PhysRevLett.58.86 (), the parallel tempering or replica exchange MC algorithm earl2005parallel (), the Wolff algorithm PhysRevLett.62.361 (), and others, the WL algorithm allows one to obtain information on the canonical distribution for given temperature . The values of for the same system may range from several units to several orders of magnitude. For example, in the Ising model, there are only two ground states on a square lattice of spins; i.e., ; however, in this case .

The results of investigation of the efficiency and the convergence of the method Xu2015 (); Belardinelli2008 (); Maerzke2012 (), as well as ways for its improvement PhysRevE.72.025701 (); Bornn2013749 (); Bauer2010 (), are regularly published in the scientific literature. There also exist other constraints that have not yet been completely investigated, for example, a constraint on the speed of construction of flat histograms (see, for example, PhysRevLett.92.097201 ()). The renormalization scheme of the WL method was substantiated in PhysRevE.75.046701 (). Among other problems formulated in PhysRevE.75.046701 (); Belardinelli2007 (), which are very important for the practical application of the method, we can highlight the following. When a histogram can be considered as flat? What is the relationship between the value of the modification factor and the calculation error? Does there exist any universal method for controlling the convergence rate of the WL algorithm?

The energy space may have a very rough structure, strong irregularities, steps, and forbidden gaps and discontinuities, whose neglect may have a significant effect on the calculation error in thermodynamic quantities. Such landscapes can be observed, for example, in spin glasses and spin ice, which, owing to their internal structure or lattice topology, can be characterized by the absence or even unattainability of equilibrium due, in particular, to the macroscopic degeneracy of the ground state. Therefore, the application of the WL algorithm to the calculation of equilibrium thermodynamic properties of spin systems that can pass to the state of spin glass or spin ice would be very interesting. Note that, in the physics of phase transitions, the problem of transition from the paramagnetic state to the spin glass state has not yet been completely analyzed.

In the present study, we consider in detail serial and parallel WL algorithms, try to answer the questions posed above, and produce a solution to the problem of calculation of the specific heat of the square Ising model by the WL method to substantiate the efficiency of the algorithm, a solution to the problem of the DOS of SSI with a finite number of spins by a parallel WL algorithm, a solution for fully connected models of hexagonal spin lattices with long-range (LR) dipoleâ dipole interaction, and a solution for spin quasilattices with finite interaction radius.

## Ii Serial Method of Calculation

Just as the Metropolis method landau2000guide (), the WL method belongs to the group of MC methods. These methods are based on conditional random walks over the space of states. However, the algorithms used in these methods differ by the approach to sampling, i.e., by the method of constructing a sampled population, a part of the universe of states, that is covered by a single experiment or by a series of serial or parallel experiments. The WL method uses an equiprobable sampling to construct the DOS. The probability distribution of energy levels, or the DOS, is represented as a histogram , see Fig. 1.

In addition (see, for example, PhysRevE.84.065702 (); kalyan2016joint ()), there exists a possibility to calculate a multidimensional distribution of the DOS of any thermodynamic quantity measured in a numerical experiment and to obtain its mean value

(1) |

The histogram is updated at every step irrespective of whether the configuration is accepted or rejected. This approach allows one to significantly reduce the operation time of the algorithm.

The operation principle and the details of implementation of a serial WL algorithm are described in the fundamental works PhysRevLett.86.2050 (); Landau2004 (), see also silant2011 (); shchure2014 (). In addition to , a histogram is formed that serves as an indicator that all possible energy levels of the system are visited uniformly. The following initial values are defined: , , , the modification factor is ; the effect of the accuracy of this factor on the result and the convergence rate will be discussed below.

The operation of the algorithm consists in step-by-step generation of a chain of states of the system. At every MC step, a candidate for a new configuration is suggested that differs from the previous one by one flipped spin. With regard to the probability , either the new configuration is accepted (), or the old one is returned :

(2) |

The updating probability

(3) |

depends on the probability distribution of energy states obtained at previous iterations of the algorithm, where and are energies corresponding to the old and new configurations.

The histograms are updated irrespective of whether the configuration is accepted or rejected according to the rule

(4) |

In practice, in order to achieve the accuracy limits guaranteed by the computer, it is desirable to use instead of . Then the updating condition of is expressed as

(5) |

The sampling is performed until the histogram becomes flat with certain accuracy. At this moment, the WL step is terminated, and new values of and , are set. Thus, the algorithms starts a new sampling cycle with higher updating accuracy of . The accuracy of flatness is determined by the maximum deviation of each element of the histogram from its mean value. Usually, it is 80%. According to the results presented in PhysRevLett.86.2050 (), an increase in the threshold value deteriorates the convergence of the algorithm. No increase in the accuracy of the result is observed. The answer to the question of how the uniformity of the distribution of the auxiliary histogram affects the accuracy of calculations requires additional investigation. However, one can assert that an increase in the uniformity of distribution leads to a decrease in the efficiency of the algorithm.

A new distribution is formed on the basis of the previous one. The modification factor f simultaneously serves as a criterion of completion and an indicator of computing speed. The algorithm persists while , i.e., until a certain minimum value is attained. Taking into account that every MC step increases by a factor f, by varying , we actually vary the accuracy of due changing the number of full WL cycles. This fact determines a balance between the accuracy of and the efficiency of the algorithm.

## Iii Parallel Sampling

The parallel WL method is a good symbiosis of two MC algorithms: the serial WL method and the replica exchange method, which is sometimes called a parallel tempering method earl2005parallel (). In the classical WL algorithm, a single serial process produces a random walk over the whole space of states from to . In the parallel computation scheme PhysRevLett.110.210603 (); vogel2014scalable (), the energy space is divided into overlapping intervals (energy windows or bands). The value of overlap may vary depending on the conditions of the problem PhysRevLett.110.210603 (); vogel2014scalable () or on the form of the Gibbs distribution. Each energy interval is independently calculated by processes.

During operation, each process performs random walks in its unique replica of the system and has individual independent distributions and and its own value of .

During operation, processes and that occur in adjacent windows and have their own respective histograms and exchange configurations and with energies and , respectively, with probability

(6) |

Note that a random process in one window chooses for exchange a random process in one of two adjacent windows. The exchange of configurations between processes in windows that are not adjacent is not admitted because, as a result of such an exchange, energy would fall outside the limits of the energy window with greater probability.

Upon completing the exchange procedure, according to the rule (4), one should update the histograms of each process involved in the exchange. To reduce the calculation error, when a given value of the flatness criterion of is reached in all processes in one energy window, is averaged over all processes of this window, and each process takes the average value vogel2014scalable (). All values of the histogram are set equal to zero.

### iii.1 Energy Intervals

To increase the efficiency of the method, one can divide the space of states from to into equal intervals with a given value of overlapping. Such energy âbandsâ can be sampled independently. The boundaries are defined as follows:

(7) |

(8) |

(9) |

where - is the width of the energy band, and are the lower and upper boundaries of the th interval, where the numbering of starts from , and is the overlapping level of the intervals. For each model, the size of the window is chosen empirically to obtain optimal computing time.

A nonoptimal value of may reduce the convergence rate of the algorithm but does not affect the accuracy of calculations (except for the case of ). A too large value of leads to a significant increase in the size of the window; as a consequence, a large number of states should be visited. For a small value of , the probability of failure during configuration exchange noticeably increases. One of energies to be exchanged may fall outside the energy window of the process that accepts a configuration. The optimal value of for the SSI model is .

In PhysRevLett.110.210603 (), the authors additionally proposed to increase dynamically and reduce as the energy level increases. This approach is justified in the case of systems with rough energy landscape. Near sharp jumps in , one should reduce the values of and , while, in regions with smoother landscape, one should increase the size of the window.

### iii.2 Falling Outside the Energy Window

During operation of the algorithm, it is critically important to ensure that the values of energy obtained satisfy the condition

(10) |

Otherwise, problems associated with the return of the process to its own energy window may arise. Before starting the algorithm, one should balance the system, i.e., make MC steps until each random process accepts a replica whose energy falls in its own energy window. For a known configuration of , one should start balancing from this configuration and make subsequent MC steps taking only those configurations that increase energy.

When condition (10) is violated, one should cancel the configurations () during sampling. The exchange of and configurations between processes and can occur only provided the following conditions are satisfied:

(11) |

Otherwise the exchange is cancelled.

### iii.3 Concatenation and Renormalization of Histograms

A preliminary result of the parallel WL method is a set of histograms bounded by appropriate intervals , which should be concatenated into a resulting histogram . The histograms may significantly differ in height in different energy windows.

The process of concatenation consists in determining a gluing point followed by the normalization of each th part to guarantee the continuity of . An example of a normalized set is shown in Fig. 1. Two histograms are concatenated in a region where their growth rates (angles of inclination) coincide, i.e., where the values of derivatives are maximally close to each other.

To study the models given below, we apply simple two-point differentiation:

(12) |

This approach is effective in number series with small variance, see Fig. 2b. However, it cannot be applied to the analysis of systems with âdenseâ and âcoarseâ energy landscape in view of the large dispersion of the histogram values, see Fig. 2a. In such a case, one should apply smoothing techniques before differentiation.

The authors of vogel2014scalable () propose using a five-point approximation with given sequence length when the distribution over states is smooth and has hardly any forbidden gaps. The scope of approximation should be large enough to obtain sufficiently smooth derivatives. Examples of differentiation with scopes of 1 and 10 are shown in Figs. 2c and 2d, respectively. In case of an unsatisfactory result of matching the DOS regions, one can apply additional line-smoothing techniques.

To calculate specific heat, one can use a relative DOS; therefore, we calculated the specific heat of a magnetic system on the basis of an unnormalized histogram as a derivative of the internal energy:

(13) |

If one needs to calculate other thermodynamic quantities, for example, entropy, one should obtain true values of the degeneracy multiplicity of levels. In the Ising ferromagnetic model, this can be done easily: if the first element in the histogram is , then all the other elements can be normalized by . In the Ising model, the ground state has multiplicity two; the same is observed in the SSI model irrespective of the dipoleâdipole interaction radius.

Since the parallel method of multicanonical sampling needs further development, all the resulted presented below were obtained by the serial WL method.

## Iv Results and Discussion

### iv.1 Critical Phenomena and Phase Transitions in -Vector Models

The divergence of thermodynamic quantities and the stability or instability of thermodynamic behavior for certain parameters (for example, the divergence of specific heat and magnetic susceptibility for given values of the external magnetic field and temperature) are usually called âcritical phenomenaâ Ma1976 (). They accompany a second-order phase transition that occurs at certain critical temperature . It is assumed that, in this transition, two phases are transformed into each other without energy dissipation.

The continuous character of thermodynamic functions describing a second-order phase transition allows one to investigate the behavior of a system near . It was even established that the behavior near the phase-transition temperature can be described by power laws whose exponents depend on a very small number of parameters such as the dimension and the number of degrees of freedom of the system Fisher1974 (); Wilson1979 (); Fisher1998 () (although the symmetry of the lattice and the interaction radius may change the character of critical behavior Fisher1974 (); Ma1976 (); Jose1977 (); Kosterlitz1978 ()). In fact, this means that the critical behavior of most systems can be described by two Hamiltonians involving microscopic pair interactions Vaz2008 ().

1. The Potts model of states in which every th spin can be in one of Q possible discrete states (orientations) . If two neighboring spins have the same orientation, then they make a contribution to the total energy of the configuration; otherwise, they make no contribution:

(14) |

where is the Kronecker delta denotes summation over all interacting neighbors, and is the exchange constant.

2. The -vector model, which is characterized by the fact that a spin can be in a continuum of states:

(15) |

where is an -dimensional unit vector at site that interacts with vector situated at cite .

In this work, to study the thermodynamics of phase transitions in complex spin systems by the WL method considered above, we sampled spaces of states of - and -vector models. The latter include the SSI model, the hexagonal spin ice (HSI) model, and plane quasilattices of spin ice, the so-called âspin snow.â

The Ising and SSI models were investigated by short-range (SR) interaction; i.e., the interaction between spins was taken into account only in the first coordination sphere. The hexagonal model was calculated with both SR and LR interactions; the model with LR interaction took into account interactions between all spins.

Below we will show that, depending on the radius and character of interaction, the same -vector model may exhibit a significant variation in the behavior of specific heat near , which is accompanied by the emergence of several maxima. We also show that, in the case of HSI in the limit of an infinite number of particles, the shape of a sample with free boundary conditions does not affect the temperature behavior of specific heat near irrespective of the interaction radius.

### iv.2 The Ising Model on a Square Lattice

To check the accuracy of the WL algorithm, we compared its results with the exact solution for the simplest and well-studied Ising model on a simple square lattice.

The Hamiltonian of the model is

(16) |

where and take values .

The size of the square lattice was spins with periodic boundary conditions. The flatness criterion of the histogram was . We used 24 cycles of WL sampling; i.e., at the last iteration, we had . During operation, we performed altogether MC steps. The flatness of the histogram was checked every steps.

Figure 3 presents a comparison of the temperature dependence of specific heat obtained by the WL method and the exact analytic solution for Ferdinand1969 (). The temperature of the peak was for the WL method and for the exact solution. The values of specific heat at the peak were and in dimensionless units, respectively. The relative deviation of specific heat from the exact solution was less than . Figure 3 presents the data of one run of the algorithm. The deviation between independent runs of the algorithm was less than the point size in the figure.

### iv.3 Artificial Arrays of Superspin Ice

The nanoarchitecture of artificial superspin (macrospin) ice represents an array of single-domain ferromagnetic nanoparticles of given geometry arranged on a substrate or in the bulk of a nonmagnetic material. The nanoparticles are usually made of thin-film magnetic materials, for example, permalloy or cobalt. The volume of a nanoisland is on the order of , and the magnetization is on the order of . The shape anisotropy due to the geometry of the nanoisland leads to the parallel ordering of atomic magnetic moments along the long axis, while the behavior of the magnetic moment of the nanoisland is Ising-like (superspin or macrospin), as demonstrated in a number of experimental works wang2006artificial (); PhysRevB.77.094418 (); Lederman1993 (); PhysRevB.80.140409 ().

In Nefedev2010 (); Ivanov2011 (); Nefedev2011 (), the authors presented the results of numerical simulation of magnetic force microscopy (MFM) experiments on the basis of which one can conclude that the MFM images of nanoparticle arrays (see, for example, RevModPhys.85.1473 () do not contain non-singledomain states of nanoislands and that the deviation from the uniform magnetization due to the presence of edge effects can be neglected even for square nanoparticles. For rectangular thin-film islands with significant shape anisotropy, the edge effects are still weaker. This give grounds for the application of the Ising spin model with dipole interaction to describe the equilibrium thermodynamics of dipole ice.

The model of artificial dipole ice was presented in our recent paper Shevchenko2017 (). The Hamiltonian of dipole interaction is

(17) |

and the energy of dipoleâdipole interaction for an pair of spins is

(18) |

where is a dimensional coefficient, is the total magnetic moment of the island, and is the lattice parameter. Each nanoisland is considered as an Ising-like magnetic dipole.

All the solutions of the problems on the temperature behavior of specific heat presented in this paper have been obtained within equilibrium thermodynamics, i.e., under the assumption that the energy barriers controlled by the shape anisotropy or by other kinds of anisotropies are overcome.

#### iv.3.1 Square spin ice

Using a parallel WL method, we obtained the DOS for the SSI model (Fig. 4) with a size of unit cells (9940 dipoles) and SR interaction (four neighbors) with free boundary conditions. The energy space was partitioned into 80 equal intervals with 80% overlapping. The algorithm was run on 400 computer cores, with five replicas per interval.

The results obtained by the WL method are in good quantitative and qualitative agreement with the results obtained by the Metropolis algorithm to within experimental error (Fig. 5). The peak temperature is . However, this value is different from the result obtained in silva2012thermodynamics (). The difference in temperature is attributed to the different number of interacting nearest neighbors used in the model. The dipoleâdipole interaction is of antiferromagnetic character; therefore, it leads to a decrease in the temperature of phase transition with increasing number of interacting neighbors in SSI.

#### iv.3.2 Hexagonal spin ice

We investigated two series of samples consisting of different numbers of particles formed into arrays of hexagonal spin ice of square and hexagonal (âfrom the centerâ) shape with different interaction radii. In the SR case, the interaction only with four (or less, on the boundaries) nearest neighbors was taken into account. In the LR case, each spin interacted with each other. The thermodynamics of the samples was calculated in 24 steps of a serial WL algorithm.

The model of spin ice in the form of a square sample of hexagonal lattice with free boundary conditions, dipole interaction (17), and cells, i.e., magnetic moments, is schematically shown in Fig. 6. The temperature behavior of the specific heat of (), (), and () models is demonstrated in Figs. 7 and 8.

In case of LR interaction (Fig. 7), the temperature behavior of specific heat exhibits an anomalous character, and two temperature peaks are observed; i.e., the critical behavior of the system of dipoles is substantially changed near . The first peak increases with increasing number of particles, while the second peak decreases. Figure 8 shows that the specific heat peak in models with SR interaction decreases with increasing number of particles.

The model of spin ice for a hexagonal sample of hexagonal lattice with free boundary conditions and dipole interaction (17) is schematically demonstrated in Fig. 9. The temperature behavior of the specific heat of samples of this shape is shown in Figs. 10 and 11.

The values of specific heat peaks for square and hexagonal samples of hexagonal lattice are shown in Fig. 12 as a function of the size of the system. We constructed approximations by the following equations:

1) SR (square) ,

2) SR (hexagonal) ,

3) LR (square) second peak ,

4) LR (hexagonal) second peak ,

5) LR (square) first peak ,

6) LR (hexagonal) first peak .

Figure 12 implies that, as , one can neglect the effect of boundaries on the critical phenomena near . The specific heat peaks converge in square and hexagonal samples. In the SR case, the height of the specific heat peak tends to a value of ; in the LR case, the specific heat has a singularity in a low-temperature region, and the height of the second peak tends to a value of as the size of the system increases. This implies that there is no phase transition in the model with SR, whereas, in the model with LR, a phase transition occurs. The question of how many neighbors should one take into account in order not to lose the main thermodynamic properties of the system requires additional analysis.

#### iv.3.3 Spin snow

Magnetic nanoparticles provide a truly inexhaustible variety of many-body systems for both experimenters and theoreticians, because the possibilities of designing nanoparticle arrays of arbitrary architecture are actually infinite. A free choice of the geometry of an array, as well as the geometry of its nanocomponents, allows the design of any lattices that have no analogs in nature. We found an analogy with an infinite variety of snowflakes of ordinary water ice that form ordinary snow. The present section is titled âspin snow,â because by âsnowâ one usually means a set of ice crystals. In this section, as an example, we demonstrate the possibilities of simulation of arbitrary samples of spin snow with a given number of neighbors, i.e., with arbitrary interaction radius.

In particular, we analyzed the quasilattices shown in Fig. 13 by a serial WL method. The calculations were performed with regard to different coordination spheres, depending on which the number of neighbors in various samples ranged from 14 to 57. It turns out that such lattices exhibit very interesting thermodynamic properties.

For example, in a sample consisting of 400 spins (Fig. 13a), the radius of the coordination sphere was taken equal to four lattice constants. Therefore, the number of neighbors ranged from 18 to 20 depending on the location of the spins. It is noteworthy that the temperature behavior of specific heat exhibits three finite peaks. The graph of specific heat of a sample consisting of 108 spins (Fig. 13b) has two characteristic global peaks. The radius of the coordination sphere was taken equal to six lattice constants, and the number of neighbors ranged from 37 to 57. We observed a nonzero value of residual entropy .

The temperature behavior of the specific heat of a spin snow lattice consisting of 432 spins (Fig. 13c) has two peaks. The interaction between dipoles was restricted to four coordination spheres, which amounted to from 14 to 16 neighbors per particle.

## V Conclusions

In this study, we have considered serial and parallel implementations of the WL algorithm. The scheme of numerical calculation of the DOS proposed here involves random walks over states, and the transition probability between states depends on the ratio of the occurrence frequency of a previous state to the occurrence frequency of a new state. Thus, the WL method suggests more frequent visit of energy levels with smaller value of . The method shows good convergence of results to the exact solution of the twodimensional Ising model, as well as convergence to the results of independent simulation of the SSI model by the Metropolis method.

The divergence of specific heat for in Fig. 5 is due to an insufficient number of Metropolis steps.

A parallel WL algorithm has many variable parameters, such as , , , , , and the size and overlapping of energy windows, which affect the speed and accuracy of calculations. However, currently there does not exist a universal algorithm for choosing these parameters. Fine tuning of the algorithm is performed depending on the complexity of the energy landscape and the type of the system considered.

The anomalous behavior of specific heat observed in the arrays should be investigated in detail to find out the reasons behind the variation of critical behavior near the critical temperature of phase transition. Here two directions of research are possible that are aimed at the theoretical description of the second-order phase transition.

The first direction can obviously be related to the fact that the classical idea of the second-order phase transition suggests mutual transformation between two phases, which is usually explained in terms of percolation theory, i.e., in terms of the emergence of a âpercolation cluster.â The presence of one or several peaks near can be attributed to the existence of several transitions between possibly coexisting phases and the creation and subsequent transformation of a percolation cluster under LR dipole interaction. The question of the concept of percolation cluster in systems with alternating-sign LR interaction also remains open.

The second direction of research can be associated with the determination of the critical behavior, which is no longer described within universal and simple power laws; i.e., the hypothesis of the universality of the critical behavior of second-order phase transitions within one of -vector models requires verification.

## Acknowledgments

We are grateful to Far Eastern Federal University for providing supercomputer facilities and to Kseniya Valerâevna Shapovalova for preparing interesting samples of spin snow. We are also grateful to Prof. Yutaka Okabe for valuable and helpful discussions. This work was supported by a grant from the President of the Russian Federation for young scientists and graduate students within the program for the development of the priority direction âStrategic Information Technologies, Including the Creation of Supercomputers and Software Development,â project nos. SP-946.2015.5 and SP-1675.2015.5, and under state task âMagnetic Properties Multiscale Structure of Nanomaterialsâ (task no. 3.7383.2017/B4).

## References

- (1) T. F. Middleton and D. J. Wales, J. Chem. Phys. 118, 4583 (2003).
- (2) F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001)
- (3) D. P. Landau, S. H. Tsai and M. Exler, Am. J. Phys. 72, 1294 (2004).
- (4) F. Wang and D. P. Landau, Phys. Rev. E 64, 056101 (2001).
- (5) B. J. Schulz, B. Kurt and M. Müller, Int. J. Mod. Phys. C 13, 477 (2002).
- (6) A. Proykova and D. Stauffer, Open Phys. 3, 209 (2005).
- (7) M. Troyer, S. Wessel and F. Alet, Phys. Rev. Lett. 90, 120201 (2003).
- (8) A. Malakis, A. Peratzakis and N. G. Fytas, Phys. Rev. E 70, 066128 (2004).
- (9) G. Brown and T. C. Schulthess, J. Appl. Phys. 97, 10E303 (2005).
- (10) B. J. Schulz, K. Binder and M. Müller, Phys. Rev. E 71, 046705 (2005).
- (11) S. Reynal and H. T. Diep, Phys. Rev. E 72, 056710 (2005).
- (12) F. Calvo, Phys. Rev. E 82, 046703 (2010).
- (13) Y. L. Xie, P. Chu, Y. L. Wang et al., Phys. Rev. E 89, 013311 (2014).
- (14) F. Calvo and P. Parneix, J. Chem. Phys. 119, 256 (2003).
- (15) J. Snider and C. Y. Clare, Phys. Rev. B 72, 214203 (2005).
- (16) T. Vogel, Y. W. Li, T. Wüst et al, Phys. Rev. Lett. 110, 210603 (2013).
- (17) D. Jayasri, V. S. S. Sastry and K. P. N. Murthy, Phys. Rev. E 72, 036702 (2005).
- (18) C. Desgranges and J. Delhommelle, J. Chem. Phys. 130, 244109 (2009).
- (19) V. T. Ngo, D. T. Hoang and H. T. Diep, Phys. Rev. E 82, 041123 (2010).
- (20) W. Kwak, J. Jeong, J. Lee et al., Phys. Rev. E 92, 022134 (2015).
- (21) AA Caparica, S. A. Leão and C. J. DaSilva, Physica A: Statistical Mechanics and its Applications 438, 447 (2015).
- (22) J. Liu, B. Song, Y. Yao et al., Phys. Rev. E 90, 042715 (2014).
- (23) N. Rathore and J. J. de Pablo, J. Chem. Phys. 116, 7225 (2002).
- (24) S. H. Tsai, F. Wang and D. P. Landau, Brazilian journal of physics 36, 635 (2006).
- (25) M. A. De Menezes and A. R. Lima, Physica A: Statistical Mechanics and its Applications 323, 428 (2003).
- (26) V. Mustonen and R. Rajesh, J. Phys. A: Math. Gen. 36, 6651 (2003).
- (27) C. Zhou and R. N. Bhatt, Phys. Rev. E 72, 025701 (2005).
- (28) R. E. Belardinelli and V. D. Pereyra, Phys. Rev. E 75, 046701 (2007).
- (29) N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth et al., J. Chem. Phys. 21, 1087 (1953).
- (30) R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58, 86 (1987).
- (31) D. J. Earl and M. W. Deem, Phys. Chem. Chem. Phys. 7, 3910 (2005).
- (32) U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
- (33) L. Bornn, P. E. Jacob and P. Del Moral, J. Comput. Graph. Stat. 22, 749 (2013).
- (34) B. Bauer, E. Gull, S. Trebst et al., Journal of Statistical Mechanics: Theory and Experiment 2010, No. P01020 (2010).
- (35) S. Xu, X. Zhou, Y. Jiang et al., Sci. China Phys., Mech. Astron., 58, 1 (2015).
- (36) R. E. Belardinelli, S. Manzi and V. D. Pereyra, Phys. Rev. E 78, 067701 (2008).
- (37) K. A. Maerzke, L. Gai, P. T. Cummings et al., J. Chem. Phys. 137, 204105 (2012).
- (38) P. Dayal, S. Trebst, S. Wessel et al., Phys. Rev. Lett. 92, 097201 (2004).
- (39) R. E. Belardinelli and V. D. Pereyra, J. Chem. Phys., 127, 184105 (2007).
- (40) R. F. Wang, C. Nisoli, R. S. Freitas et al., Nature 439, 303 (2006).
- (41) Y. Qi, T. Brintlinger and J. Cumings, Phys. Rev. B 77, 094418 (2008).
- (42) M. Lederman, G. A. Gibson and S. Schultz, J. Appl. Phys. 73, 6961 (1993).
- (43) G. Möller and R. Moessner, Phys. Rev. B 80, 140409 (2009).
- (44) K. V. Nefedev, Y. P. Ivanov and A. A. Peretyatko, Methods and Tools of Parallel Programming Multicomputers: Second Russia-Taiwan Symposium, MTPP 2010, Vladivostok, Russia, May 16-19, 2010, Revised Selected Papers 6083, 260 (2010).
- (45) Y. P. Ivanov, K. V. Nefedev, A. I. Iljin et al., Journal of Physics: Conference Series 266, 012117 (2011).
- (46) K. V. Nefedev, Y. P. Ivanov, A. A. Peretyatko et al., Solid State Phenomena 168, 325 (2011).
- (47) C. Nisoli, R. Moessner and P. Schiffe, Rev. Mod. Phys. 85, 1473 (2013).
- (48) D. P Landau and K. Binder, Cambridge Univ. Press, 384 (2000).
- (49) G. Brown, Kh. Odbadrakh, D. M. Nicholson et al., Phys. Rev. E 84, 065702 (2011).
- (50) M. S. Kalyan, R. Bharath, V. S. S. Sastry et al., J. Stat. Phys. 163, 197 (2016).
- (51) I. A. Silantâeva and P. N. Vorontsov-Velâyaminov, Vychisl. Metody Programm. 12, 397 (2011).
- (52) L. N. Shchur, Mekh., Upravl. Inform. 6 (6), 160 (2014).
- (53) T. Vogel, Y. W. Li, T. Wüst et al., Phys. Revi. E 90, 023302 (2014).
- (54) R. C. Silva, F. S. Nascimento, L. A. S. Mól et al., New Journal of Physics 14, 015008 (2012).
- (55) C.A.F. Vaz, J.A.C. Bland, G. Lauhoff, Reports on Progress in Physics, 71, 056501 (2008).
- (56) M.E. Fisher, Rev. Mod. Phys. 46 597 (1974).
- (57) K.G. Wilson, Sci. Am. 241 140 (1979).
- (58) M.E. Fisher, Rev. Mod. Phys. 70 653 (1998).
- (59) E.E. Stanley, Rev. Mod. Phys. 71 S358 (1999)
- (60) S.K. Ma, Modern Theory of Critical Phenomena (Reading, MA: Benjamin), (1976).
- (61) J.V. Jose, L.P. Kadanoff, S. Kirkpatrick and D.R. Nelson, Phys. Rev. B 16 1217 (1977).
- (62) J.M. Kosterlitz and D.J. Thouless, Progress in Low Temperature Physics vol. VII B ed. D.F. Brewer (Amsterdam: North-Holland) p. 371 (1978).
- (63) Y. Shevchenko, A. Makarov, K. Nefedev, Phys. Let. A, 381 428 (2017).
- (64) A.E. Ferdinand and M.E. Fisher, Phys. Rev. B 185, (832) 1969.