Dynamic phase transitions under the presence of quenched randomness
Abstract
We present an extensive study of the effects of quenched disorder on the dynamic phase transitions of kinetic spin models in two dimensions. We undertake a numerical experiment performing Monte Carlo simulations of the squarelattice randombond Ising and BlumeCapel models under a periodically oscillating magnetic field. For the case of the BlumeCapel model we analyze the universality principles of the dynamic disorderedinduced continuous transition at the lowtemperature regime of the phase diagram. A detailed finitesize scaling analysis indicates that both nonequilibrium phase transitions belong to the universality class of the corresponding equilibrium random Ising model.
pacs:
64.60.an, 64.60.De, 64.60.Cn, 05.70.Jk, 05.70.LnI Introduction
In the last sixty years our understanding of equilibrium critical phenomena has developed to a point where wellestablished results are available for a wide variety of systems. In particular, the origin and/or the difference between equilibrium universality classes is by now well understood. This observation also partially holds for systems under the presence of quenched disorder. However, far less is known for the physical mechanisms underlying the nonequilibrium phase transitions of manybody interacting systems that are far from equilibrium and clearly a general classification of nonequilibrium phase transitions into nonequilibrium universality classes is missing.
We know today that when a ferromagnetic system, below its Curie temperature, is exposed to a timedependent oscillating magnetic field, it may exhibit a fascinating dynamic magnetic behavior Tome (). In a typical ferromagnetic system being subjected to an oscillating magnetic field, there occurs a competition between the time scales of the appliedfield period and the metastable lifetime, , of the system. When the period of the external field is selected to be smaller than , the timedependent magnetization tends to oscillate around a nonzero value, which corresponds to the dynamically ordered phase. In this region, the timedependent magnetization is not capable of following the external field instantaneously. However, for larger values of the period of the external field, the system is given enough time to follow the external field, and in this case the timedependent magnetization oscillates around its zero value, indicating a dynamically disordered phase. When the period of the external field becomes comparable to , a dynamic phase transition takes place between the dynamically ordered and disordered phases.
Throughout the years, there have been several theoretical Lo (); Zimmer (); Acharyya1 (); Chakrabarti (); Acharyya2 (); Acharyya3 (); Buendia1 (); Buendia2 (); Fujisaka (); Jang1 (); Jang2 (); Shi (); Punya (); Riego (); Keskin1 (); Keskin2 (); Robb1 (); Deviren (); Yuksel1 (); Yuksel2 (); Vatansever1 () and experimental studies He (); Robb (); Suen (); Berger (); Riego1 () dealing with dynamic phase transitions, as well as with the hysteresis properties of magnetic materials. The main conclusion emerging is that both the amplitude and the period of the timedependent magnetic field play a key role in dynamical critical phenomena (in addition to the usual temperature parameter). Furthermore, the characterization of universality classes in spin models driven by a timedependent oscillating magnetic field has also attracted a lot of interest lately Sides1 (); Sides2 (); Korniss (); Buendia3 (); Park (); Park2 (); Tauscher (); Buendia4 (); Vatansever2 (); Vatansever_Fytas (). Some of the main results are listed below:

The critical exponents of the kinetic Ising model were found to be compatible to those of the equilibrium Ising model at both two (2D) and three dimensions (3D) Sides1 (); Sides2 (); Korniss (); Vatansever2 (); Park ().

Buendía and Rikvold using soft Glauber dynamics estimated the critical exponents of the 2D Ising model and provided strong evidence that the characteristics of the dynamic phase transition are universal with respect to the choice of the stochastic dynamics Buendia3 ().

The role of surfaces at nonequilibrium phase transitions in Ising models has been elucidated by Park and Pleimling: The nonequilibrium surface exponents were found to be different than equilibrium critical surface ones Park2 ().

Experimental evidence by Riego et al. Riego1 () and numerical results by Buendía and Rikvold Buendia4 () verified that the equivalence of the dynamic phase transition to an equilibrium phase transition is limited to the area near the critical period and for zero bias.

Numerical simulations by Vatansever and Fytas showed that the nonequilibrium phase transition of the spin1 BlumeCapel model belongs to the universality class of the equilibrium Ising counterpart (at both 2D and 3D) Vatansever_Fytas (). General and very useful features of the dynamic phase transition of the BlumeCapel model can also be found in Refs. Buendia1 (); Keskin1 (); Keskin2 (); Deviren (); ShiWei (); Acharyya4 ().
The above results in 2D and 3D kinetic Ising and BlumeCapel models establish a mapping between the universality principles of the equilibrium and dynamic phase transitions of spin1/2 and spin1 models. They also provide additional support in favor of an earlier investigation of a GinzburgLandau model with a periodically changing field Fujisaka (), as well as with the symmetrybased arguments of Grinstein et al. in nonequilibrium critical phenomena Grinstein ().
Motivated by the current literature, in the present work we attempt to shed some light on the effect of quenched disorder on dynamic phase transitions. To the best of our knowledge, with the exception of a few meanfield and effectivefield theory treatments of the problem Gupinar12 (); Gupinar12b (); Akinci12 (); Vatansever13 (); Vatansever13b (); Vatansever15 (), there exists no dedicated (numerical) work. However, what we have mainly learned from the previous studies on the topic is that the dynamic character of a typical magnetic system driven by a timedependent magnetic field sensitively depends on the amount of disorder, accounting for reentrant phenomena and dynamic tricritical points Akinci12 (). In the current work we use as testcase platforms for our numerical experiment the Ising and BlumeCapel models on the square lattice under a timedependent magnetic field, diffusing disorder in the ferromagnetic exchange interactions. For the case of the BlumeCapel model we focus on the disorderedinduced continuous dynamic transition at the lowtemperature regime of the phase diagram. In a nutshell, our results indicate that the dynamic phase transitions of both the randombond Ising and BlumeCapel models belong to the universality class of the equilibrium random Ising model.
The outline of the remainder parts of the paper is as follows: In Sec. II we introduce the disordered versions of the Ising and BlumeCapel models and in Sec. III the thermodynamic observables necessary for the application of the finitesize scaling analysis. The details our simulation protocol are given in Sec. IV and the numerical results and discussion in Sec. V. Finally, Sec. VI contains a summary of our conclusions.
Ii Models
We consider the squarelattice randombond Ising and BlumeCapel (BC) models under a timedependent oscillating magnetic field, described by the following Hamiltonians
(1) 
and
(2) 
In the above Eqs. (1) and (2) indicates summation over nearest neighbors and the spin variable takes on the values for the Ising and for the BC model, respectively. The couplings denote the random ferromagnetic exchange interactions, drawn from a bimodal distribution of the form
(3) 
Following Refs. Malakis1 (); Malakis2 (); Malakis3 (), we choose and , so that defines the disorder strength; for the pure systems are recovered. A clear advantage of using the bimodal distribution (3) is that the critical temperature of the random Ising model is exactly known from duality relations as a function of the disorder strength via Fisch (); Kinzel ()
(4) 
For the case of the BlumeCapel Hamiltonian (2) denotes the crystalfield coupling that controls the density of vacancies (). For vacancies are suppressed and the model becomes equivalent to the Ising model. Finally, the term corresponds to a spatially uniform periodically oscillating magnetic field, so that all lattice sites are exposed to a squarewave magnetic field with amplitude and half period Korniss (); Buendia3 (); Park ().
A brief description of the BlumeCapel model’s phase diagram together with some necessary pinpoints of the current literature with respect to the effect of disorder on its critical behavior may be useful here: The phase diagram of the equilibrium pure BlumeCapel model in the crystalfield – temperature plane consists of a boundary that separates the ferromagnetic from the paramagnetic phase. The ferromagnetic phase is characterized by an ordered alignment of spins. On the other hand, the paramagnetic phase can be either a completely disordered arrangement at high temperature or a spin gas in a spin dominated environment for low temperatures and high crystal fields. At high temperatures and low crystal fields, the ferromagneticparamagnetic transition is a continuous phase transition in the Ising universality class, whereas at low temperatures and high crystal fields the transition is of firstorder character Capel (); Blume (). The model is thus a classic and paradigmatic example of a system with a tricritical point Lawrie (), where the two segments of the phase boundary meet. A detailed reproduction of the phase diagram of the model can be found in Ref. Zierenberg () and an accurate estimation of the location of the tricritical point has been given in Ref. Kwak (): . A lot of work has been also devoted in understanding the effects of quenched bond randomness on the universality aspects of the BlumeCapel model, especially in two dimensions, where any infinitesimal amount of disorder drives the firstorder transition at the lowtemperature regime to a continuous transition. Quantitative phase diagrams of the randombond BlumeCapel model at equilibrium have been constructed in Refs. Malakis1 (); Malakis2 () and, more recently, a dedicated numerical study at the firstorder transition regime revealed that the induced under disorder continuous transition belongs to the universality class of the random Ising model with logarithmic corrections Fytas17 ().
Iii Observables
In order to determine the universality aspects of the kinetic randombond Ising and BlumeCapel models, we shall consider the halfperiod dependencies of various thermodynamic observables. The main quantity of interest is the periodaveraged magnetization
(5) 
where the integration is performed over one cycle of the oscillating field. Given that for finite systems in the dynamically ordered phase the probability density of becomes bimodal, one has to measure the average norm of in order to capture symmetry breaking, so that defines the dynamic order parameter of the system. In the above Eq. (5), is the timedependent magnetization per site
(6) 
where defines the total number of spins and the linear dimension of the lattice.
To characterize and quantify the transition using finitesize scaling arguments we must also define quantities analogous to the susceptibility in equilibrium systems. The scaled variance of the dynamic order parameter
(7) 
has been suggested as a proxy for the nonequilibrium susceptibility, also theoretically justified via fluctuationdissipation relations Robb1 (). Similarly, one may also measure the scaled variance of the periodaveraged energy
(8) 
so that can be considered as the corresponding heat capacity. Here denotes the cycleaveraged energy corresponding to the cooperative part of the Hamiltonians (1) and (2). With the help of the dynamic order parameter we may define the fourthorder Binder cumulant Sides1 (); Sides2 ()
(9) 
a very useful observable for the characterization of universality classes Binder81 ().
Iv Simulation Details
We performed Monte Carlo simulations on square lattices with periodic boundary conditions using the singlesite update Metropolis algorithm Metropolis (); Binder (); Newman (). This approach, together with the stochastic Glauber dynamics Glauber:63 (), consists the standard recipe in kinetic Monte Carlo simulations Buendia3 (). Let us briefly outline below the steps of our computer algorithm:

A lattice site is selected randomly among the options.

The spin variable located at the selected site is flipped, keeping the other spins in the system fixed.

The energy change originating from this spin flip operation is calculated using the Hamiltonians of Eqs. (1) of (2) as follows: , where denotes the energy of the system after the trial switch of the selected spin and corresponds to the total energy of the system with the old spin configuration. The probability to accept the proposed spin update is given by:
(10) 
If the energy is lowered, the spin flip is always accepted.

If the energy is increased, a random number is generated, such that : If this number is less than or equal to the calculated Metropolis transition probability the selected spin is flipped. Otherwise, the old spin configuration remains unchanged.
Using the above scheme we simulated system sizes within the range . For each system size independent realizations of the disorder have been generated – see Fig. 1 for characteristic illustrations of disorder averages and their relative variance – and for each random sample the following simulation protocol has been used: The first periods of the external field have been discarded during the thermalization process and numerical data were collected and analyzed during the following periods of the field. The time unit in our simulations was one Monte Carlo step per site (MCSS) and error bars have been estimated using the jackknife method Newman (). To give a flavor of the actual CPU time of our computations we note that the simulation times needed for a single disorder realization of the kinetic Ising model on a single node of a Dual Intel Xeon E52690 V4 processor were 6 hours for and 11 days for . The analogous CPU times for the kinetic BlumeCapel model were 3 hours and 9 days for and , respectively. For the Ising model we fixed the value of the disorder strength to , whereas for the BlumeCapel model we focused on the value in the originally firstorder regime selecting now following Refs. Malakis2 (); Malakis3 (). Appropriate choices of the magneticfield strength, , and the temperature, and , ensured that the system lies in the multidroplet regime Park (). Here, Fisch (); Kinzel () and Malakis2 (); Malakis3 () are the equilibrium critical temperatures of the Ising and BlumeCapel models for the particular choices of and considered in this work.
For the fitting process on the numerical data we restricted ourselves to data with . As usual, to determine an acceptable we employed the standard test of goodness of fit Press (). Specifically, the value of our test is the probability of finding an value which is even larger than the one actually found from our data. We considered a fit as being fair only if .
V Results and discussion
As a starting point let us describe shortly the mechanism underlying the dynamical ordering in kinetic ferromagnets (here, under the presence of quenched randomness), as exemplified in Figs. 2  4 for the Ising model, and Figs. 5  8 for the BlumeCapel model. In both cases, results for a single realization of the disorder are shown over a system size of .
Figure 2 presents the time evolution of the magnetization and Fig. 3 the period dependencies of the dynamic order parameter of the kinetic randombond Ising model. Several comments are in order: For rapidly varying fields, Fig. 2(a), the magnetization does not have enough time to switch during a single half period and remains nearly constant for many successive field cycles, as also illustrated by the black line in Fig. 3. On the other hand, for slowly varying fields, Fig. 2(c), the magnetization follows the field, switching every half period, so that , as also shown by the blue line in Fig. 3. In other words, whereas in the dynamically disordered phase the ferromagnet is able to reverse its magnetization before the field changes again, in the dynamically ordered phase this is not possible and therefore the timedependent magnetization oscillates around a finite value. The competition between the magnetic field and the metastable state is captured by the halfperiod parameter (or by the normalized parameter , with being the metastable lifetime Park ()). Obviously, plays the role of the temperature in the equilibrium system. Now, the transition between the two regimes is characterized by strong fluctuations in , see Fig. 2(b) and the evolution of the red line in Fig. 3. This behavior is indicative of a dynamic phase transition and occurs for values of the half period close to the critical one (otherwise when ). Of course, since the value MCSS used for this illustration is slightly above , see also Fig. 10, the observed behavior includes as well some nonvanishing finitesize effects.
Some additional spatial aspects of the transition scenarios described above via the configurations of the local order parameter are shown in Fig. 4. Below , see panel (a), the majority of spins spend most of their time in the state, i.e., in the metastable phase during the first half period, and in the stable equilibrium phase during the second half period, except for equilibrium fluctuations. Thus most of the and the system is now in the dynamically ordered phase. On the other hand, when the period of the external field is selected to be bigger than the relaxation time of the system, above , see panel (c), the system follows the field in every half period, with some phase lag, and at all sites . The system lies in the dynamically disordered phase. Near and the expected dynamic phase transition, there are large clusters of both and values, within a sea of , as shown in Fig. 4(b).
Although the discussion above concentrated on the Ising case, analogous description and relevant conclusions may be drawn also for the dynamical ordering of the disorderedinduced continuous transition of the BlumeCapel model, as depicted in Figs. 5  7. Note that in this case the critical halfperiod of the system has been estimated to be (see also Fig. 11 below). However, we should underline here that for the case of the BlumeCapel model the value of the local order parameter does not distinguish between random distributions of and and regions of . To bring out this distinction, we present in Fig. 8 similar snapshots of the dynamic quadrupole moment over a full cycle of the external field, , where . In the spin1 BlumeCapel model the density of the vacancies is controlled by the crystalfield coupling and, thus, the value of the dynamic quadrupole moment changes depending on Vatansever_Fytas (). We point out that in Fig. 8, except for the red areas, the regions enclosed by finite values demonstrate the role played by the crystalfield coupling in the BlumeCapel model.
To further explore the nature of the dynamic phase transitions encountered in the above disordered kinetic models we performed a finitesize scaling analysis using the observables outlined in Sec. III. Previous studies in the field indicated that although scaling laws and finitesize scaling are tools that have been designed for the study of equilibrium phase transitions, they can be successfully applied as well to far from equilibrium systems Sides1 (); Sides2 (); Korniss (); Buendia3 (); Park ().
As an illustrative example for the case of the kinetic randombond Ising model and for a system size of we present in the main panel of Fig. 9 the finitesize behavior of the dynamic order parameter and in the lower inset the emerging dynamic susceptibility [see Eq.(7)]. The dynamic order parameter goes from a finite value to zero values as the half period increases showing a sharp change around the value of the half period that can be mapped to the respective peak in the plot of the dynamic susceptibility. The location of the maxima in may be used to define suitable pseudocritical half periods, denoted hereafter as . The corresponding maxima may be analogously denoted as . We also measured the energy and its scaled variance, the heat capacity [see Eq.(8)]. The upper inset of Fig. 9 shows the halfperiod dependency of the energy of the same system and the relevant heat capacity. In this case the maxima may be denoted as .
We start the presentation of our finitesize scaling analysis with Ising case. In the main panel of Fig. 10(a) we present the size evolution of the peaks of the dynamic susceptibility in a loglog scale. The solid line is a fit of the form Ferrenberg ()
(11) 
providing an estimate for the magnetic exponent ratio , in excellent agreement to the Ising value . The shift behavior of the corresponding peak locations is plotted in the inset of Fig. 10(a) as a function of . The solid line shows a fit of the usual shift form Fisher (); Privman (); Binder92 ()
(12) 
where defines the critical half period of the system and is the critical exponent of the correlation length. The obtained values and are listed also in the panel and, in particular, the value of the critical exponent appears to be in very good agreement to the value of the 2D equilibrium Ising model. This finding strongly supports the claim that the kinetic Ising model under the presence of quenched bond randomness shares the universality class of its corresponding equilibrium counterpart. Ideally, we would also like to observe the double logarithmic scaling behavior of the maxima of the heat capacity . Indeed, as it is shown in the main panel of Fig. 10(b), the data for the maxima of the heat capacity are adequately described by a fit of the form
(13) 
as predicted by Ref. Dotsenko () for the random Ising universality class. As a comparison, we plot the same data with respect to the simple logarithm of the system size in the corresponding inset. It is obvious that a fit , as shown by the solid line, does not capture the full scaling behavior.
So, where do we stand at this point: We have shown that the universality class of the dynamic phase transition encountered in an Ising model under the presence of quenched bond disorder is equivalent to that of its equilibrium counterpart with the inclusion of logarithmic corrections in the scaling of the heat capacity. We turn now our discussion to the dynamic phase transition of the BlumeCapel model with bond disorder. As mentioned previously in Sec. II, only very recently the claims of universality violation in the equilibrium randombond BlumeCapel model have been dispelled and it was shown that the induced under disorder continuous transition belongs to the universality class of the random Ising model Fytas17 (). We therefore expect, or at least hope, that the results presented in the current work will also be relevant to this reignited problem, yet from a nonequilibrium perspective.
The scaling aspects of the dynamic phase transition of the kinetic randombond BlumeCapel model at are shown in Fig. 11, following fully the presentation and analysis style of Fig. 10 and excluding the data for that suffer from strong finitesize effects. In this case an estimate is obtained for the magnetic exponent ratio , again compatible within errors to the Ising value . From the shift behavior of the corresponding pseudocritical half periods [inset of Fig. 11(a)] the critical halfperiod and the correlationlength exponent are estimated to be and , respectively. Again the estimate of supports the scenario presented above in Fig. 10 for the criticality in the dynamic phase transition of the randombond kinetic Ising model. Last but not least, in Fig. 11(b) the maxima of the heat capacity are plotted versus (main panel) and (inset) and as in the Ising case are much better described by the double logarithmic fit (13).
An alternative test of universality comes from the study of the fourthorder Binder cumulant defined in Eq. (9) for the case of the dynamic order parameter. In Fig. 12 we present our numerical data of for the kinetic randombond Ising (main panel) and BlumeCapel (inset) models. In both panels the vertical dashed line marks the critical halfperiod value of the system and the horizontal dotted line the universal value of the 2D equilibrium Ising model Salas (). Certainly, the crossing point is expected to depend on the lattice size (as it also shown in the figure) and the term universal is valid for given lattice shapes, boundary conditions, and isotropic interactions Selke_2 (); Selke_3 (). However, the data shown in Fig. 12 support, at least qualitatively, another instance of equilibrium Ising universality, since in both panels the crossing point is consistent to the value . We should note here that Hasenbusch et al. presented very strong evidence that the critical Binder cumulant of the equilibrium 2D randomly sitediluted Ising model maintains its puresystem value Hasenbusch08 (). In this respect, a dedicated study along the lines of Ref. Hasenbusch08 () for an accurate estimation of in the kinetic randombond Ising and BlumeCapel models would be welcome, but certainly goes beyond the scope of the current work.
Vi Conclusions
In the present work we investigated the effect of quenched disorder on the dynamic phase transition of kinetic spin models in two dimensions. In particular, we considered the squarelattice Ising and BlumeCapel models under a periodically oscillating magnetic field, the latter at its lowtemperature regime where the pure equilibrium system exhibits a firstorder phase transition. Using Monte Carlo simulations and finitesize scaling techniques we have been able to probe with good accuracy the values of the critical exponent and the magnetic exponent ratio , both of which were found to be compatible to those of the equilibrium Ising ferromagnet. An additional study of the scaling behavior of the heatcapacity revealed the double logarithmic divergence expected for the universality class of the random Ising model. To conclude, although universality is a cornerstone in the theory of critical phenomena, it stands on a less solid foundation for the case of nonequilibrium systems and for systems subject to quenched disorder. In the current work we have studied two systems where both of the above complications merge, yet arriving to the simplest scenario. We hope that our work will stimulate further research in the field of nonequilibrium critical phenomena at both numerical and analytical directions.
Acknowledgements.
The authors would like to thank P. A. Rikvold and W. Selke for many useful comments on the manuscript. The numerical calculations reported in this paper were performed at TÜBİTAK ULAKBIM (Turkish agency), High Performance and Grid Computing Center (TRUBA Resources).References
 T. Tomé and M.J. de Oliveira, Phys. Rev. A 41, 4251 (1990).
 W.S. Lo and R.A. Pelcovits, Phys. Rev. A 42, 7471 (1990).
 M.F. Zimmer, Phys. Rev. E 47, 3950 (1993).
 M. Acharyya and B.K. Chakrabarti, Phys. Rev. B 52, 6550 (1995).
 B.K. Chakrabarti and M. Acharyya, Rev. Mod. Phys. 71, 847 (1999).
 M. Acharyya, Phys. Rev. E 56, 1234 (1997).
 M. Acharyya, Phys. Rev. E 69, 027105 (2004).
 G.M. Buendía and E. Machado, Phys. Rev. E 58, 1260 (1998).
 G.M. Buendía and E. Machado, Phys. Rev. B 61, 14686 (2000).
 H. Fujisaka, H. Tutu, and P.A. Rikvold, Phys. Rev. B 63, 036109 (2001).
 H. Jang, M.J. Grimson, and C.K. Hall, Phys. Rev. E 68, 046115 (2003).
 H. Jang, M.J. Grimson, and C.K. Hall, Phys. Rev. B 67, 094411 (2003).
 X. Shi, G. Wei, and L. Li, Phys. Lett. A 372, 5922 (2008).
 A. Punya, R. Yimnirun, P. Laoratanakul, and Y. Laosiritaworn, Physica B 405, 3482 (2010).
 P. Riego and A. Berger, Phys. Rev. E 91, 062141 (2015).
 M. Keskin, O. Canko, and U. Temizer, Phys. Rev. E 72, 036125 (2005).
 M. Keskin, O. Canko, and Ü. Temizer, J. Exp. Theor. Phys. 104, 936 (2007).
 D.T. Robb, P.A. Rikvold, A. Berger, and M.A. Novotny, Phys. Rev. E 76, 021124 (2007).
 B. Deviren and M. Keskin, J. Magn. Magn. Mater. 324, 1051 (2012).
 Y. Yüksel, E. Vatansever, and H. Polat, J. Phys.: Condens. Matter 24, 436004 (2012).
 Y. Yüksel, E. Vatansever, U. Akinci, and H. Polat, Phys. Rev. E 85, 051123 (2012).
 E. Vatansever, Phys. Lett. A 381, 1535 (2017).
 Y.L. He and G.C. Wang, Phys. Rev. Lett. 70, 2336 (1993).
 D.T. Robb, Y.H. Xu, O. Hellwig, J. McCord, A. Berger, M.A. Novotny, and P.A. Rikvold, Phys. Rev. B 78, 134422 (2008).
 J.S. Suen and J.L. Erskine, Phys. Rev. Lett. 78, 3567 (1997).
 A. Berger, O. Idigoras, and P. Vavassori, Phys. Rev. Lett. 111, 190602 (2013).
 P. Riego, P. Vavassori, and A. Berger, Phys. Rev. Lett. 118, 117202 (2017).
 S.W. Sides, P.A. Rikvold, and M.A. Novotny, Phys. Rev. Lett. 81, 834 (1998).
 S.W. Sides, P.A. Rikvold, and M.A. Novotny, Phys. Rev. E 59, 2710 (1999).
 G. Korniss, C.J. White, P.A. Rikvold, and M.A. Novotny, Phys. Rev. E 63, 016120 (2000).
 G.M. Buendía and P.A. Rikvold, Phys. Rev. E 78, 051108 (2008).
 H. Park and M. Pleimling, Phys. Rev. E 87, 032145 (2013).
 H. Park and M. Pleimling Phys. Rev. Lett. 109, 175703 (2012).
 K. Tauscher and M. Pleimling, Phys. Rev. E 89, 022121 (2014).
 G.M. Buendía and P.A. Rikvold, Phys. Rev. B 96, 134306 (2017).
 E. Vatansever, arXiv:1706.03351.
 E. Vatansever and N.G. Fytas, Phys. Rev. E 97, 012122 (2018).
 X. Shi and G. Wei, Phys. Scr. 89, 075805 (2014).
 M. Acharyya and A. Halder, J. Magn. Magn. Mater. 426, 53 (2017).
 G. Grinstein, C. Jayaprakash, and Y. He, Phys. Rev. Lett. 55, 2527 (1985).
 G. Gulpinar and E. Vatansever, J. Stat. Phys. 146, 787 (2012).
 G. Gulpinar, E. Vatansever, and M. Agartioglu, Physica A 391, 3574 (2012).
 U. Akinci, Y. Yüksel, E. Vatansever, and H. Polat, Physica A 391, 5810 (2012).
 E. Vatansever, U. Akinci, and H. Polat, J. Magn. Magn. Mater. 344, 89 (2013).
 E. Vatansever, U. Akinci, Y. Yüksel, and H. Polat, J. Magn. Magn. Mater. 329, 14 (2013).
 E. Vatansever and H. Polat, Phys. Lett. A 379, 1568 (2015).
 A. Malakis, A.N. Berker, I.A. Hadjiagapiou, and N.G. Fytas, Phys. Rev. E 79, 011125 (2009).
 A. Malakis, A.N. Berker, I.A. Hadjiagapiou, N.G. Fytas, and T. Papakonstantinou, Phys. Rev. E 81, 041113 (2010).
 A. Malakis, A.N. Berker, N.G. Fytas, and T. Papakonstantinou, Phys. Rev. E 85, 061106 (2012).
 R. Fisch, J. Stat. Phys. 18, 111 (1978).
 W. Kinzel and E. Domany, Phys. Rev. B 23, 3421 (1981).
 H.W. Capel, Physica (Amsterdam) 32, 966 (1966).
 M. Blume, Phys. Rev. 141, 517 (1966).
 I.D. Lawrie and S. Sarbach, in: C. Domb, J.L. Lebowitz (Eds.), Phase Transitions and Critical Phenomena, Vol. 9 (Academic Press, London, 1984).
 J. Zierenberg, N.G. Fytas, M. Weigel, W. Janke, and A. Malakis, Eur. Phys. J. Special Topics 226, 789 (2017).
 W. Kwak, J. Jeong, J. Lee, and D.H. Kim, Phys. Rev. E 92, 022134 (2015).
 N.G. Fytas, J. Zierenberg, P.E. Theodorakis, M. Weigel, W. Janke, and A. Malakis, Phys. Rev. E 97, 040102(R) (2018).
 K. Binder, Z. Phys. B: Condens. Matter 43, 119 (1981); Phys. Rev. Lett. 47, 693 (1981).
 N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953).
 D.P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, Cambridge, U.K., 2000).
 M.E.J. Newman and G.T. Barkema, Monte Carlo Methods in Statistical Physics (Oxford University Press, New York, 1999).
 R.J. Glauber, J. Math. Phys. 4 , 294 (1963).
 A. Aharony and A. B. Harris, Phys. Rev. Lett. 77, 3700 (1996).
 S. Wiseman and E. Domany, Phys. Rev. Lett. 81, 22 (1998); Phys. Rev. E 58 2938 (1998).
 W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical Recipes in C, 2nd ed. (Cambridge University Press, Cambridge, 1992).
 A.M. Ferrenberg and D.P. Landau, Phys. Rev. B 44, 5081 (1991).
 M.E. Fisher, Critical Phenomena, edited by M.S. Green (Academic, London, 1971).
 V. Privman, Finite Size Scaling and Numerical Simulation of Statistical Systems (World Scientific, Singapore, 1990).
 K. Binder, Computational Methods in Field Theory, edited by C.B. Lang and H. Gausterer (Springer, Berlin, 1992).
 A.E. Ferdinand and M.E. Fisher, Phys. Rev. 185, 832 (1969).
 V.S. Dotsenko and V.S. Dotsenko, Sov. Phys. JETP Lett. 33, 37 (1981).
 J. Salas and A.D. Sokal, J. Stat. Phys. 98, 551 (2000).
 W. Selke, J. Stat. Mech. P04008 (2007).
 W. Selke and L.N. Shchur, Phys. Rev. E 80, 042104 (2009).
 M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, Phys. Rev. E 78, 011110 (2008).