Dimensionality crossover in critical behaviour of ultrathin ferromagnetic films

Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

Contents lists available at ScienceDirect

Journal of Magnetism and Magnetic Materials journal homepage: www.elsevier.com/locate/jmmm

Dimensionality crossover in critical behaviour of ultrathin ferromagnetic films Pavel V. Prudnikov n, Vladimir V. Prudnikov, Maria A. Menshikova, Natalia I. Piskunova Department of Theoretical Physics, Omsk State University, Omsk 644077, Russia

art ic l e i nf o

a b s t r a c t

Article history: Received 20 October 2014 Received in revised form 12 March 2015 Accepted 25 March 2015 Available online 27 March 2015

We propose the model which takes account of magnetocrystalline anisotropy effects in thin magnetic films. The dimensionality crossover from two-dimensional monolayer to three-dimensional system in multilayer magnetic films is studied using a Monte Carlo technique. Finite-size scaling is applied for the determination of the critical characteristics as a function of film thickness. The transition to intermediate planar phase is discussed. & 2015 Elsevier B.V. All rights reserved.

Keywords: Ultrathin magnetic films Magnetocrystalline anisotropy Dimensionality crossover Monte Carlo simulation Critical behaviour

1. Introduction The behaviour of ultrathin magnetic films has become of great technological importance due to the applications in magnetic storage devices [1]. It has been theoretically suggested that the highest areal density advantage for heat assisted magnetic recording (HAMR) can be achieved when the maximum heating temperature is closed to or exceed the Curie temperature [2]. HAMR is the most appropriate technology to achieve magnetic recording densities 1 TBit/in2 [3,4]. In this respect, it is important to understand the temperature evolution of magnetization in thin magnetic films, especially at temperatures close to or above the Curie temperature. Magnetic order in thin ferromagnetic films is very complex due to a strong influence of the shape and the magnetocrystalline anisotropies of the sample. In the past decade, a considerable amount of experimental results on different aspects of magnetism in ultrathin films has appeared [5]. Nevertheless it is difficult to reach general conclusions even in seemingly basic things such as the kind of magnetic order at low temperatures. In view of this complexity, theoretical work on simplified models and computer simulations is essential for rationalizing and guiding new experimental work. In the vicinity of the critical temperature Tc, the thermodynamic observables associated with statistical models display n

Corresponding author. E-mail address: [email protected] (P.V. Prudnikov).

http://dx.doi.org/10.1016/j.jmmm.2015.03.075 0304-8853/& 2015 Elsevier B.V. All rights reserved.

universal characteristics, which may be parametrized in terms of critical exponents. These quantities tend to zero or infinity at the transition and depend only on the spatial dimensionality of the system, the range of the interactions and the number of components of the order parameter. The dimensionality aspects of magnetic and structural phase transitions represent one of the key problems of ultrathin film [6,7]. For magnetic systems, the spin dimensionality as well as the spatial extension determines the universality class, giving rise to a great number of ordering phenomena on different length scales. Furthermore, there are transition regions not represented by any universality class with corresponding critical exponents, but representing something in between. For example, the dependence of the critical exponents of thin magnetic films with thickness exhibit such a transition, in which the exponents are continuously changes with increasing film thickness of the layers, from two-dimensional (2D) Ising (β ¼ 0.125) to three-dimensional (3D) Heisenberg (β ¼ 0.364) [8] behaviour. The crossover from 2D to 3D critical exponents has been observed in thin films of Ni on W(110) [9] and of Co on Cu(111) [10] as their thickness is increased. Thus, studying the dimensional crossover of a system as its thickness is increased often provides a significant degree of insight. It is well known now the fundamental role of competing interactions in the emerging features of low-dimensional systems. Among a wide number of numerical and theoretical investigations on equilibrium and dynamical properties of several model Hamiltonian of low-dimensional magnets, Heisenberg-like models are one of the most widely used to approach real magnetic

78

P.V. Prudnikov et al. / Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

materials. In fact, anisotropic versions of the Heisenberg model describe several compounds as K2NiF4 [11], BaCo2(AsO4)2 [12], CoCl2-GIC [13] and Rb2CrCl4 [14]. In early experiments on Fe/Cu(100) films [15,16], a spin reorientation transition (SRT) from a region with perpendicular magnetization to one with in-plane magnetization was observed. SRT transitions in films have actively been studied both theoretically [17,18] and experimentally [19–21], recently. In this paper we investigate the magnetic phase transition in anisotropic Heisenberg thin films using extensive Monte Carlo (MC) simulations. We have found the dimensional crossover of magnetization m and susceptibility χm from 2D to 3D like with increasing film thickness. Estimated values of the critical exponents for different thicknesses demonstrate crossover from 2D Ising universality class to 3D Heisenberg through 3D Ising class. The paper is organized in the following manner. The anisotropic Heisenberg model and the Monte Carlo simulation scheme are discussed in the next section (Section 2). The numerical results for universal dimensionality behaviour of thin films with increasing film thickness and for SRT transition are given in Section 3. The paper ends with concluding remarks and summary in Section 4.

i, j

(1)

where Si = (Six, Siy, Siz ) is a unit vector in the direction of the classical magnetic moment at lattice site i, the sum is extended over nearest-neighbor pairs on the cubic lattice, J being the exchange constant, and Δ characterizes the amount of anisotropy. Thus, Δ ¼0 corresponds to the isotropic Heisenberg case, Δ ¼1 – the Ising case. The effective anisotropy constant Δ (N) as a function of film thickness N was chosen from experimental studies of the Curie temperature TC for thin films of Ni(111)/W(110) [9] with different thicknesses of Ni film. Microscopic nature of anisotropy in films of Fe, Co, Ni and its dependence of film thicknesses N is determined by influence of crystalline field of substrate surface, magnetic single-ion anisotropy and dipole–dipole interaction of magnetic moments of atoms in film and their concurrence. Therefore, the calculation of anisotropy effects in magnetic films is very complicated task. It was proposed that Δ is proportional to the critical temperature for different film thicknesses. In approximation procedure the fact that Ni films with a large number of layers demonstrate bulk critical properties of 3D isotropic Heisenberg magnets was used [23,24]. The simulations were carried out for films with sizes L × L × N with the use of periodic and free boundary conditions for the inplane and out-plane directions, respectively. L × L represents the number of spins in each layer of the film and N is the number of layers. We considered films with L¼ 32, 48, 64 and N ranging from a monolayer to 32 layers. Temperature T of system is changed in interval [0.01; 5.01] J/kB with step ΔTstep = 0.02. As a starting configuration we always used a completely ordered ferromagnetic state. We measured the magnetization

1/2 ⎡ ⎛ Ns ⎞2 ⎛ Ns ⎞2 ⎛ Ns ⎞2⎤ ⎢⎜ x⎟ z⎟ ⎥ y⎟ ⎜ ⎜ ⎢ ⎜∑ Si ⎟ + ⎜∑ Si ⎟ + ⎜∑ Si ⎟ ⎥ ⎠ ⎝ i ⎠ ⎝ i ⎠⎦ ⎣⎝ i

, (2)

out-plane magnetization

1 Ns

mz =

Ns

∑ Siz

, (3)

i

the in-plane magnetization

1 Ns

m∥ =

1/2 ⎡ ⎛ Ns ⎞2 ⎛ Ns ⎞2⎤ ⎢⎜ x⎟ y⎟ ⎥ ⎜ ⎢ ⎜∑ Si ⎟ + ⎜∑ Si ⎟ ⎥ ⎠ ⎝ i ⎠⎦ ⎣⎝ i

, (4)

and an orientational order parameter [25,26]

nhα − nvα nhα + nvα

Oα =

,

(5)

where Ns = NL2 is a total number of spins in film, angle brackets denote the statistical averaging, α ∈ {x, y, z} , nh and nv are the number of horizontal and vertical pairs of nearest neighbor spins with antialigned {x, y} components, respectively,

⎧ ⎪

∑ ⎨1 − sgn[S α (rx , ry ), S α (rx ⎪

r

We have performed Monte Carlo simulations of macroscopic spin system with crystalline structure of ultrathin film which is described by the anisotropic Heisenberg model [22] with Hamiltonian

H = − J ∑ [(1 − Δ (N))(Six S jx + Siy S jy ) + Siz S zj ],

m=

nhα =

2. Model and methods

1 Ns

nvα =

⎩ ⎧ ⎪

⎫ , + 1, ry )] ⎬ ⎪ ⎭

∑ ⎨1 − sgn[S α (rx , ry ), S α (rx , ry ⎪

r

⎩

⎪

⎫ . + 1)] ⎬ ⎪ ⎭ ⎪

We define the magnetic susceptibility

χm ∼

[〈m2〉]

−

[〈m〉]2

(6)

χm as (7)

and orientational susceptibility χO as

χO ∼ [〈Oz2 〉] − [〈Oz 〉]2

(8)

Temperature dependence of the susceptibility has been calculated for different lattice sizes to estimate the critical temperature Tc. The position of the susceptibility maximum allowed us to determine range of values of the critical temperature. The spin configurations of the films are updated using the Swendsen–Wang cluster algorithm [27]. The spin-flip algorithm for O(n) models was proposed in [28]. The spin system is divided into clusters. The bonds between sites in cluster are created with probability 1 − exp[ − 2J (Si r)(S j r)/T ] only if the condition (Si r)(S j r) > 0 is true, where r is a random unit vector. After that each cluster is flipped with probability 1/2. The temperature dependencies of the Binder cumulant U4 (T )

U4 =

1⎛ [〈m4 〉] ⎞ ⎜3 − ⎟ 2⎝ [〈m2〉]2 ⎠

(9)

were calculated to clarify the critical temperatures of second-order phase transition in samples. The scaling dependence of the cumulant

U4 (L, T ) = u [L1/ ν (T − Tc )].

(10)

makes it possible to determine the critical temperature Tc from the coordinate of the intersection points of the curves specifying the temperature dependence U4 (L, T ) for different L. We consider finite size scaling form for film geometry [29] to find how m and χ scale with the size L and thickness N of the systems and use this to extract the effective critical exponents from our results. The basic finite-size scaling ansatz [30] rests on

P.V. Prudnikov et al. / Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

an assumption that only a single correlation length ξ is needed to describe the critical properties of thin films. Hence, the empirical scaling forms for m and χ at some fixed N's can be written as

˜ (L1/ ντ , N) 〈m (T , N) 〉 = L−β / νm χ (T , N) = Lγ / νχ˜ (L1/ ντ , N)

(11)

where γ, β, and ν are the effective critical exponents associated with χ, m, and ξ, respectively. It should be clarified that for finite N the presence of non-zero anisotropy in films always leads to the asymptotic critical behavior of the Ising-like type with d ¼2. However, for large N this asymptotic critical region is unobservable and only effective exponents can be detected in modeling and experiments. In this paper the crossover of effective exponents is studied for different values of N. For monolayer system with N ¼ 1, the effective exponents are ˜ are the critical exponents for the 2D system. The functions χ˜ and m scaling functions for a given N and τ, where τ = T /Tc − 1 is the reduced temperature. The effective exponent 1/ν can be extracted from the derivative of the cumulant U4 with respect to L at Tc owing to its variation with system size as L1/ν

dU4 ∼ L1/ ν. dT

(12)

The scaling dependence of the cumulant allowed us to calculate values of the critical exponent ν. We detected two types of phase transitions in films on basis of calculated temperature dependencies of the magnetization m(T), m∥ (T ), mz(T) and susceptibilities χm (T ), χO (T ). The transition from ferromagnetic phase to paramagnetic phase was detected for all thicknesses of film. This type of transition is characterized by peak of magnetic susceptibility χm (T ) which is presented in Fig. 1 for films with N ¼2 and N ¼ 10. We also calculated the heat capacity C (T ) ∼ [〈E 2 (T ) 〉] − [〈E (T ) 〉]2. The obtained results are presented in Fig. 1a for thin films with thickness N ¼ 2. Temperature dependence of heat capacity demonstrates the singular critical behaviour at the same temperature as the magnetic susceptibility. However, measurements show that the singularity of heat capacity is smoothed over with increasing thickness of Heisenberg film. Though in vicinity of critical temperature mz(T) tends to zero, but critical behaviour of films is determined by mz fluctuations. Monte Carlo data which confirm this fact are presented in Fig. 2. This phenomenon corresponds with results of studies of 2D anisotropic Heisenberg model in [22]. The observable dimensional crossover in magnetic characteristics with increasing of film thickness will

79

be discussed in Section 3.1. For the case of Ising films, a first study of crossover from d ¼2 to d ¼3 with increasing N has been given in [31] by Binder. For films with N = 9–22 we detected the spin reorientation phase transition from an out-of-plane magnetization ordering to an in-plane configuration. It is characterized by behaviour of the magnetization presented in Fig. 1b. For comparison, we demonstrate the temperature dependencies of m(T), m∥ (T ), mz(T) and susceptibilities χm (T ), χO (T ) in the same figure. The curve for susceptibility χm (T ) is characterized by two peaks in Fig. 1b. The first peak at T = 0.91(4) is pronounced and the fast decrease of the outof-plane magnetization is occurred near this temperature. The susceptibility χO (T ) has a maximum at the same temperature, that is why first peak corresponds to spin reorientation transition from ferromagnetic phase to planar ferromagnetic phase. Spin reorientation phase transition will be discussed in Section 3.2.

3. Results of Monte Carlo modeling 3.1. Dimensionality effects in critical behaviour of ultrathin films The temperature dependencies of total magnetization m (T , N) and susceptibility χm (T , N) are presented in Fig. 3 for different thicknesses of thin film. These curves were obtained by averaging over 3000 samples for each size of film. We divided all obtained data against N into four groups. The typical behaviour of m (T , N) and χm (T , N) for each group is demonstrated in Fig. 3. The measurement error is within the symbol. For first group of films with N ≤ 5 (Fig. 3a) the values of universal critical characteristics are essentially the same as for 2D model and from this fact it can be implied that these ultrathin films fall into the 2D universality class. We calculated the temperature dependencies of ln(L)/ln[U4/(T − TC )] for different temperatures above critical and for different linear sizes L of film (Fig. 4). These dependencies allow us to calculate the effective value of critical exponent ν. This value was calculated using linear approximation of linear size L to infinity and linear approximation of temperature to critical value [33]. Note that if Eq. (11) correctly encapsulates the nature of magnetic critical behaviour in films, we can extract the values of effective exponents β /ν and γ /ν from the slopes of the log–log plots of m or χ against L at Tc. The values of exponent β can be obtained from temperature dependence of the magnetization near the critical point m ∼ (Tc − T ) β .

Fig. 1. Temperature dependencies of the total magnetization m, the in-plane magnetization m∥, the out-plane magnetization mz, χm (T ) , χO (T ) and the heat capacity C(T) for films with N ¼ 2 (a) and with N ¼ 10 (b). The error bars are smaller then the size of the symbols.

80

P.V. Prudnikov et al. / Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

Then, we consider the f (m) = L β / νm as a scaling function of

− T )/Tc . We calculated this scaling function for thin films with thickness from 2 to 5 monolayers with L ¼32, 48, 64. The scaling functions obtained for different L collapsed into a single curve with the correct values of critical temperature and effective critical exponents. The results are shown in Fig. 5. The obtained values of exponents β, ν and γ for films with thicknesses N ≤ 5 are for N ¼2: β = 0.126(8), ν = 1.010(17), γ = 1.816(69); for N ¼ 3: β = 0.128(8), ν = 1.011(27), γ = 1.770(94); for N ¼4: β = 0.126(9), ν = 0.986(21), γ = 1.713(112); for N ¼5: L1/ ν (Tc

Fig. 2. Monte Carlo data for the out-plane magnetization mz plotted vs. time [in units of Monte Carlo steps per spin (MCS/s)].

β = 0.129(9), ν = 0.972(59), γ = 1.609(150). These values are very close to exact values of the exponents β2DIsing = 1/8, ν2DIsing = 1, γ2DIsing = 7/4 (Fig. 6) for the 2D Ising model. Calculated value of

γ = 1.816(69) for N ¼ 2 is in a good agreement with experimental value γ = 1.75(2) which was measured for Fe/W(110) bilayers in [32]. We performed the calculation of the effective dimension deff based on the hyperscaling relation deff = γ /ν + 2β /ν . We have found that in case with N = 2–5 deff has a value equal to 2 within

Fig. 4. Temperature dependence of ln(L)/ln[U4/(T − TC )].

Fig. 3. Temperature dependencies of total magnetization m (T , N) and susceptibility χ (T , N) for thin films with N ¼ 1 (full symbol), N ¼6 (open symbol) (a); N ¼ 8 (full symbol), N ¼ 10 (open symbol) (b); N ¼ 15 (full symbol), N ¼ 17 (open symbol) (c); N ¼ 26 (full symbol), N ¼31 (open symbol) (d), where L ¼ 32 – (), L ¼48 – (■ ), L ¼ 64 – (▴ ). The error bars are smaller then the size of the symbols.

P.V. Prudnikov et al. / Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

81

Table 1 Critical temperatures Tc and TSRT , critical exponents ν, β, γ for different N.

Fig. 5. Dependence f (m) = L β / νm as a function of L1/ ν (Tc − T )/Tc .

N

Tc

TSRT

Δ

ν

β

γ

2 3 4 5 6 7 8 9 10 12 13 15 16 17 21 22 23 26 30 31

1.03(8) 1.15(11) 1.25(5) 1.31(10) 1.35(4) 1.39(4) 1.45(8) 1.49(10) 1.57(8) 1.63(11) 1.87(12) 3.91(20) 4.15(14) 4.09(15) 1.61(15) 1.55(14) 1.46(8) 1.43(10) 1.43(6) 1.43(5)

– – – – – – 0 0.65(2) 0.91(4) 1.43(6) 1.87(7) 2.71(7) 2.89(8) 2.87(15) 1.17(3) 0.77(3) 0 – – –

0.581 0.636 0.687 0.734 0.777 0.816 0.851 0.882 0.909 0.952 0.967 0.986 0.990 0.990 0.950 0.930 0.906 0.811 0.630 0.575

1.010(17) 1.011(27) 0.986(21) 0.972(59) 0.974(62) 0.916(43) 0.981(28) 0.861(39) 0.839(32) 0.720(29) 0.703(26) 0.632(21) 0.634(28) 0.658(27) 0.722(36) 0.706(26) 0.713(25) 0.723(18) 0.752(36) 0.759(40)

0.126(8) 0.128(8) 0.126(9) 0.129(9) 0.170(11) 0.184(5) 0.286(11) 0.294(14) 0.310(13) 0.299(11) 0.314(5) 0.324(7) 0.329(8) 0.328(6) 0.344(2) 0.352(6) 0.351(3) 0.358(1) 0.370(2) 0.368(2)

1.816(69) 1.770(94) 1.713(112) 1.609(150) 1.803(76) 1.609(81) 1.586(52) 1.426(61) 1.225(61) 1.215(27) 1.195(37) 1.141(33) 1.182(53) 1.132(42) 1.206(43) 1.264(75) 1.294(84) 1.396(121) 1.410(86) 1.414(77)

Fig. 7. Dependence of magnetization from reduced temperature τ = (Tc − T )/Tc for films with N ¼ 15 and with N¼ 17.

Fig. 6. Dimensionality crossover in values of the critical exponents γ, ν and β as functions of the film thickness N.

error bars. This confirms the 2D type of universality for ultrathin films at the first group. We calculated values of the critical exponents β (N), γ (N) and ν (N) (Fig. 6) for different thicknesses N of the magnetic film. As a result, the dimensional crossover from two-dimensional Ising model to three-dimensional Heisenberg model with increasing film thickness was found from these dependencies. So, for thicknesses with 6 ≤ N ≤ 12 the values of exponent β move to 3D values, as expected. This second group of films forms a crossover range. For films with thicknesses from 13 to 21 the values of critical exponents correspond to value β for 3D Ising model. It is the third group of films. Films with large thicknesses N ≳ 22 demonstrate the critical

behaviour typical for 3D Heisenberg model and form the fourth group. In the total, the thickness dependence of critical exponents demonstrates the crossover from 2D Ising-like behaviour to 3D Ising-like and then to 3D Heisenberg-like behaviour. Critical temperatures for each group of films are presented in Table 1. The anisotropy constant Δ is close to 1 for thin films with thicknesses from 15 to 17. This case corresponds to the Ising model. For these films we were able to select two temperature intervals corresponding to the 2D and the 3D Ising critical behaviour. In the first temperature interval near the critical point obtained values of the exponent β (L = 64) = 0.124(5) for N ¼ 15 (Fig. 7) and β (L = 64) = 0.127(7) for N ¼17 are close to value β ¼0.125 for 2D Ising model. In the second temperature interval the values of critical exponent are β (L = 64) = 0.319(3) for N ¼ 15 and β (L = 64) = 0.328(6) for N ¼17. They correspond to the 3D Ising critical behaviour with β = 0.3249(6) [34]. This is because, the system demonstrates the 3D critical behaviour when the correlation length is increased with approaching to the critical temperature in all directions before it reaches a film boundary, but after this, the correlation length continues to increase only inplane direction of the film and system shows the critical behaviour of the 2D system.

82

P.V. Prudnikov et al. / Journal of Magnetism and Magnetic Materials 387 (2015) 77–82

films the temperature of SRT has been determined from measurements of peak positions for the orientational susceptibility χO (T ). For these systems we detected that the spin-reorientation phase transition is changed by the second-order phase transition from out-plane ferromagnetic phase to paramagnetic phase with increasing of temperature. The calculated values of critical exponents demonstrate the dimensional crossover from values typical for critical behaviour of the 2D Ising model to characteristic values for critical behaviour of the 3D Heisenberg model. Results of this phenomenon are in good agreement with results of experimental investigations of the magnetic properties of Ni and Co films as a function of thickness in [9,10].

Acknowledgments Fig. 8. Critical values of the anisotropy parameter ΔC (N) as a function of the film thickness N.

3.2. Spin reorientation transitions The thickness dependence of critical temperatures is presented in Table 1 for film with thicknesses from N ¼2 to N ¼31. The notation Tc is used for temperatures of the second-order phase transition from ferromagnetic phase to paramagnetic, and TSRT is used for the spin reorientation transition from ferromagnetic phase to planar ferromagnetic phase. In the experimental [21] and theoretical [17] investigations devoted to the study of single-layer magnetic materials, it is predicted that the spin reorientation transition is a weak first-order transition. On basis of the obtained thickness dependence of TC(N) we calculated the critical values of the anisotropy parameter ΔC (N) for different thicknesses of film which are given in Fig. 8 and denoted by strikethrough squares. The comparison of these values of ΔC (N) with the meanfield values of Δ (N) in Eq. (1) which are given in Fig. 8 by solid curve shows correspondence with each other. Difference of their values is caused by contribution of critical fluctuations to TC(N) and so to ΔC (N). The large differences between values of ΔC (N) and Δ (N) are observed in crossover range with N = 14–19 which are determined by sequence of two transitions for these N to the planar XY-like phase and to the 3D Ising-like phase.

4. Conclusions We have studied the magnetic properties of Heisenberg-like ultrathin films with different thicknesses from monolayer case with N ¼1 to case with N ¼ 32 layers using extensive Monte Carlo simulations. We have detected the dimensional crossover of magnetic characteristics of films from two-dimensional to threedimensional behaviour with increasing film thickness. We have examined the critical regime of these systems in details and extracted effective critical exponents based on a finite-size scaling method for films. For ultrathin films with thicknesses N = 9–22 the spin-reorientation phase transition was discovered from an out-of-plane magnetization ordering to an in-plane configuration. For these

The reported study was supported by Russian Scientific Fund through Project no. 14-12-00562. The simulations were supported by the Supercomputing Center of Lomonosov Moscow State University and Joint Supercomputer Center of the Russian Academy of Sciences.

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34]

C. Chappert, A. Fert, F. Nguyen van Dau, Nat. Mater. 6 (2007) 813. A. Lyberatos, K.Y. Guslienko, J. Appl. Phys. 94 (2003) 1119. R. E Rottmayer, S. Batra, D. Buechel, et al., IEEE Trans. Magn. 42 (2006) 2417. M.A. Seigler, W.A. Challener, E. Gage, et al., IEEE Trans. Magn. 44 (2008) 119. C.A.F. Vaz, J.A.C. Bland, G. Lauhoff, Rep. Prog. Phys. 71 (2008) 056501. M. Pärnaste, M. Marcellini, E. Holmström, N. Bock, J. Fransson, O. Eriksson, B. Hjörvarsson, J. Phys.: Condens. Matter 19 (2007) 246213. E.T. Papaioannou, V. Kapaklis, A. Taroni, M. Marcellini, B. Hjörvarsson, J. Phys.: Condens. Matter 22 (2010) 236004. K. Chen, A.M. Ferrenberg, D.P. Landau, Phys. Rev. B 48 (1993) 3249. Y. Li, K. Baberschke, Phys. Rev. Lett. 68 (1992) 1208. F. Huang, M.T. Kief, G.J. Mankey, R.F. Willis, Phys. Rev. B 49 (1994) 3962. R.J. Birgeneau, J. Skalyo, Jr. Shirane, G. Shirane, J. Appl. Phys. 41 (1970) 1303. L.P. Regnault, P. Burlet, J. Rossat-Mignod, Physica 86 (1977) 660. D.G. Wiesler, H. Zabel, Phys. Rev. B 36 (1987) 7303. L.M. Castro, J.A. Plascak, A.S.T. Pires, J. Magn. Magn. Mater. 248 (2002) 62. D.P. Pappas, K.P. Kamper, H. Hopster, Phys. Rev. Lett. 64 (1990) 3179. R. Allenspach, A. Bischof, Phys. Rev. Lett. 69 (1992) 3385. M. Carubelli, O.V. Billoni, S.A. Pighín, et al., Phys. Rev. B 77 (2008) 134417. Rogelio Díaz-Méndez, R. Mulet, J. Magn. Magn. Mater. 324 (2012) 128. W. Cai, S. Muraishi, J. Shi, Y. Nakamura, W. Liu, R. Yu, Appl. Phys. A 109 (2012) 69. C. Bordel, J. Juraszek, D.W. Cooke, et al., Phys. Rev. Lett. 109 (2012) 117201. A.B. MacIsaac, K. Dé Bell, J.P. Whitehead, Phys. Rev. Lett. 80 (1998) 616. K. Binder, D.P. Landau, Phys. Rev. B 13 (1976) 1140. F. Huang, M.T. Kief, G.J. Mankey, R.F. Willis, Phys. Rev. B 49 (1994) 3962. A. Ney, A. Scherz, P. Poulopoulos, et al., Phys. Rev. B 65 (2001) 024411. I. Booth, A.B. MacIsaac, J.P. Whitehead, K. Dé Bell, Phys. Rev. Lett. 75 (1995) 950. M.C. Ambrose, R.L. Stamps, Phys. Rev. B 87 (2013) 184417. R.H. Swendsen, J.S. Wang, Phys. Rev. Lett. 58 (1987) 86. U. Wolff, Phys. Rev. Lett. 62 (1989) 361. O. Dillmann, W. Janke, M. Muller, K. Binder, J. Chem. Phys. 114 (2001) 5853. V. Privman, in: V. Privman (Ed.), Finite Size Scaling and Numerical Simulation of Statistical Systems, World Scientific, Singapore, 1990. K. Binder, Thin Solid Films 20 (1974) 367. M.J. Dunlavy, D. Venus, Phys. Rev. B. 69 (2004) 094411. P.V. Prudnikov, M.A. Medvedeva, Prog. Theor. Phys. 127 (N3) (2012) 369. W. Janke, D.A. Johnston, R. Villanova, Physica A 281 (2000) 207.