Monte Carlo simulations of dynamic phase transitions in ferromagnetic thin-films

0.57 MB·English

by  Bahadır Ozan Aktaş

Monte Carlo simulations of dynamic phase transitions in ferromagnetic thin-films
Bahadır Ozan Akta¸s∗
Dokuz Eylul University, Graduate School of Natural and Applied Sciences, TR-35160 Izmir, Turkey
Erol Vatansever and Hamza Polat
Department of Physics, Dokuz Eylu¨l University, TR-35160 Izmir, Turkey
(Dated: February 12, 2015)
5
1 By means of detailed Monte Carlo (MC) simulations, we have presented dynamic phase transi-
0 tion (DPT) properties of ferromagnetic thin-films. Thermal variations of surface, bulk and total
2 dynamical order parameters (DOP) for a film and total order parameter for the films with differ-
ent thicknesses have been examined. Opposite regimes of the critical value of reduced exchange
b
e interaction (surface to bulk ratio) Rc at which the critical temperature becomes independent of
F film thickness L has been also taken into consideration. The average magnetizations of each layer
is reversed in these regimes. Based on the results, we have confirmed that the system represents
1
a crossover behavior in between ordinary to extraordinary transition in the presence of surface
1
exchangeenhancement.

] Keywords: Monte Carlo simulations, Magnetic thin-film, Surface magnetism, Surface enhancement phe-
h nomenon
c
e
m I. INTRODUCTION inEuTe(111)films decreasesignificantly differently from
- thoseofbulklayers[1]. Violbarbosaandcoworkersfound
t
a Magnetic properties of free surfaces drastically differ that the formation of the blocks of layers with robust
st from the bulk material, because the free surface breaks magnetic structure whereas the interblock interactions
. the translational symmetry (i.e. surface atoms are em- are relatively weak in fcc-Fe on Cu(001) film [3]. More-
t
a bedded in an environment of lower symmetry than that over,theenhancedsurfacemagnetismhasbeenthefocus
m of the inner atoms and consequently the exchange con- of such systems. For instance, Gd film has been inves-
- stants between atoms in the surface region may differ tigated experimentally. The thickness-dependent spin-
d
from the bulk value). The surface enhancement phe- polarized electronic structure of strained ultrathin and
n
o nomenon in finite magnetic materials has attracted con- thin films of Gd has been investigated by Waldfried et
c siderableamountofinterestforbothexperimentalists[1– al. [17]. They found that the surface magnetic struc-
[ 6] and theorists [7–16]. turedominatesthemagneticorderingoftheultrathinGd

Applied oscillating magnetic fields, depending on the films. Withdecreasingthicknesssomebulkbandsexhibit
3
v competition between the two time scales,namely the os- increasingly more passive magnetic behavior. Skomski
2 cillation period P of the external perturbation and the and coworkers also found that Gd films exhibit a mag-
4 relaxation time τ of the sample, a dynamic symmetry netic surface transition which occurs at about above the
4 breaking may take place causing a DPT. There are two bulk Curie temperature [18].
6
cases due to the competition between these time scales: In sense of dynamic phase transitions, a great deal of
0
P < τ and P > τ. In the first case, the system can- theoreticaleffortshasalsobeendevotedtotheinvestiga-
.
1 not relax within a complete cycle of the magnetic field tiononsuchsystems. Thedetailsofsurfaceenhancement
0 oscillation, hence the instantaneous magnetization M(t) phenomenon for the films were subjected to an exter-
5
oscillates in time around a nonzero value corresponding naloscillatoryfieldhavebeenintensivelypropoundedby
1
todynamicallyorderedstate(i.e. dynamicferromagnetic Akta¸s et al. by using effective field theory (EFT) [19].
:
v phase). In other case, M(t) can follow the external field Thegeneraltrendoffrequencydispersionbelongstocrit-
Xi withsomedelay,andthesystemexhibitsadynamicpara- ical temperature coordinate of the special point for dif-
magnetic behavior. The relaxation time τ can be con- ferent frequency and amplitude values has been demon-
r
a trolled by supplied energy with several different ways: stratedintheirwork. Nonequilibriumphasetransitionin
The agency of an adjustable parameter such as the field the kinetic Ising model on a two-layer square lattice has
amplitude, strength itself, the type of the exchange in- beenexaminedbyCankoetal. [20]. Dynamicphasedia-
teractions, and the temperature. The DPT point can gramshavebeen constructedin the plane of the reduced
be controlledby tuning mentioned competing factors to- temperature versus the amplitude. Similarly, dynamic
gether with the time period of external field. magnetic behavior of a mixed Ising system on a bilayer
Experimentalpointofview,Schierleandcoworkersob- squarelattice hasbeeninvestigatedbyErta¸sandKeskin
served that the magnetizations of the outermost layers [21]. They presentedthe dynamic phase diagramsin the
reducedtemperatureandmagneticfieldamplitudeplane
and the effects of interlayer coupling interaction on the
critical behavior of the system have been investigated in
∗ [email protected] their work.
2
In recent series of works by Pleimling and coworkers, In order to simulate the system, we employ the
surface criticality at a DPT and surface phase diagram Metropolis MC simulation algorithm [25, 26] to Eq. (1)
of the three-dimensional kinetic Ising model has been onanN×N×LsimplecubiclatticewhereN =70andwe
elucidated. In the first one of these studies, Park and apply periodic (free) boundaryconditions in direction(s)
Pleimling found that the nonequilibrium surface expo- parallel (perpendicular) to film plane. We have studied
nents do not coincide with those of the equilibrium crit- ultrathin-films with thickness L = 3,4,5 together with
ical surface [22]. In addition, in three space dimensions, a relatively thicker thin-film L = 20 to observe average
the surface phase diagram of the nonequilibrium system magnetizations of each layer for selected some system
differs markedly from that of the equilibrium system. parameters. For simplicity, the exchange couplings are
The values of the critical exponents have been deter- restricted to the ferromagnetic case.
mined through finite-size scaling by Park and Pleimling Configurations were generated by selecting the sites
in their followup investigation [23]. Their results have in sequence through the lattice and making single-spin-
showedthatthe studiednonequilibriumphase transition flip attempts, which were accepted or rejected according
belongstotheuniversalityclassoftheequilibriumthree- to the Metropolis algorithm, and N ×N ×L sites are
dimensional Ising model. The surface phase diagram of visited at each time step (a time step is defined as an
thethree-dimensionalkineticIsingmodelbelowtheequi- MC step per site or simply MCS). Data were generated
libriumcriticalpoint subjectedto a periodically oscillat- over 50 independent sample realizations by running the
ing magnetic field has also been presented by Taucher simulations for 50000 MCS per site after discarding the
and Pleimling [24]. They presented that surface phase first25000steps. Thisamountoftransientstepsisfound
diagramthatinpartsstronglyresemblesthecorrespond- to be sufficient for thermalization for the whole range of
ing equilibrium phase diagram, with an ordinary transi- the parameter sets. Throughout the analysis, oscillation
tion, an extraordinary transition, and a surface transi- period of the external field is kept fixed as P =100.
tion. Thesethreelinesmeetataspecialtransitionpoint. Ourprogramcalculatestheinstantaneousvaluesofthe
Forweaksurfacecouplings,however,thesurfacedoesnot bulk and surface magnetizations M and M , and the
s b
order. total magnetization M at time t. These quantities are
T
In this regard, our task in the present work is to shed defined as
some light on the DPT properties, -especially the evolu-
tionofcrossoverpointwithfieldamplitude-offerromag- 1 Ns 1 Nb
netic thin-films in the presence of ac driving fields. In Ms(t)= si, Mb(t)= sj,
N N
this present paper, the DPT properties in the presence s Xi=1 b Xj=1 (2)
of external oscillatory field of the system are studied by N M (t)+N M (t)
s s b b
M (t)=
MC simulation. The layout of the work is as follows: T
N +N
s b
Section 2 describes the model and the MC simulation
scheme, the numerical results are reported in Section 3, where N and N denote the number of spins in the sur-
s b
the paper ends with concluding remarks in Section 4. face and bulk layers, respectively. From the instanta-
neous magnetizations, we obtain the related order pa-
rameters as follows [27]:
II. SIMULATION
1 1
Q = M (t)dt, Q = M (t)dt,
We consider a ferromagnetic thin film with thickness s P I s b P I b
(3)
L described by spin-1/2 Hamiltonian 1
Q = M (t)dt
T P I T
H=− J s s −h(t) s (1)
ij i j i
Xhiji Xi Using Eq. (1), we calculate the total energy per spin
where s =±1is a two-statespinvariable,andJ is the 1
i ij E = Hdt (4)
nearest neighbor interaction energy. The summation in tot P(N +N )I
s b
the first term is taken over only the nearest neighbor in-
teractions whereas the summation in the second term is Consequently, the specific heat is defined as
carried out over the all lattice sites. In the second term,
h(t)=h0sin(ωt)representstheoscillatingmagneticfield, dEtot
C = . (5)
where h0 and ω are the amplitude and the angular fre- dT
quency of the applied field, respectively. The period of
the oscillating magnetic field is given by P = 2π/ω. If Wealsonotethatthevalueofthebulkexchangeinterac-
the lattice sites i and j belong to one of the two sur- tion J is fixed to unity, and we also use the normalized
b
faces of the film we have J = J , otherwise J = J , surfacetobulkratioofexchangeinteractionsR=J /J ,
ij s ij b s b
where Js and Jb denote the ferromagnetic surface and as well as the reduced field amplitude H0 = h0/Jb, and
bulk exchange interactions, respectively. reduced temperature Θ=k T/J .
B b
3
III. RESULTS AND DISCUSSION neighborbonds,andaDPT canbe observedatlowtem-
peratures. When we fixed rest of the parameters except
Based on the results, we mainly focus on the effect of reduced exchange and compared the strength of the or-
external oscillatory field amplitude in surface enhance- der parametersfor the films with different thicknesses in
ment phenomenon. First, in order to obtain a general two different value of R (namely R=0.25 and 2.75), we
insight on DPT characteristics, we plot Fig. (1) and (2) see the hierarchically sequence reversalin Fig. (2) (from
respectively. In Fig. (1) the thermal variations of sur- (a)to(d)and(b)to(e)and(c)to(f)). Inbetweenthese
face, bulk and total order parameters for the films with twovalue, there shouldbe a criticalpoint ofthe reduced
thickness L = 3 are shown. Following this, total dy- exchangeatwhichallthelayersseemtooscillateinphase
namical order parameters for the films with three differ- independently from the thickness.
ent thicknesses are shown in Fig. (2). We restrict our

discussions for two value of reduced exchange R = 0.25
and 2.75. Both in Fig. (1) and (2), the transition point 1.0 L = 3
4
increases with reduced exchange interaction at constant
> 5
values of the other system parameters. Hence, relatively (t) H0 = 0.0 H0 = 0.0
morethermalagitationisneededtomakethesystemdy- M T R=0.25 R=2.75
namically disordered for more interaction. Moreover, at <0.5
a constanttemperature, surfaces areweakly ordereddue
to the scarcity of dipole-dipole interaction per site for
R=0.25value. For R=2.75,surfacesdominate against (a) (d)
bulk due to the surface enhancement. From Fig. (1), 0.0
we see that both surfaces and bulk layers of magnetic 1.0
thin-films exhibit a phase transition at a certain critical
temperature independently from the value of R.
T
Q
0.5 H0 = 0.25 H0 = 0.25
P = 100 P = 100
1.0 Surface R=0.25 R=2.75
Bulk
Total
0.8 L=3 0.0 (b) (e)

P = 100
R=0.25 R=2.75 1.0
Q0.6 H0=0.25
0.4 T
Q
0.5 H0 = 0.5 H0 = 0.5
0.2 P = 100 P = 100
R=0.25 R=2.75
(a) (b)
0.0 (c) (f)
0 2 4 6 8 0 2 4 6 8 0.0
0 2 4 6 8 0 2 4 6 8
FIG. 1. (Color online) Surface, bulkand total orderparame- FIG. 2. (Color online) Average magnetization in static case
tersforafilmL=3layersfortwodifferentregimeofreduced (a)and(d),andthetotalorderparametersforthefilmswith
exchange interaction. three different thicknesses L=3,4,5. and for three different
valuesof H0.
From Fig. (2), one can easily see the effect of external
field amplitude H0 for fixed value of the system parame- In order to make the aforementioned phenomenon
tersfromthepanel(a)to(c)and(d)to(f). Criticaltem- more clear, we plot the dynamical order parameters of
perature exhibits a decreasing behavior with increasing each layers for a film with L = 20 in Fig. (3). For
H0 asaconsequenceofthewell-knownfollowingphysical this purpose, we choose a constant temperature value at
mechanism: For small amplitude values, the energy sup- which the system is well-below the transition point and
plied by the external oscillatory field cannot break the the thermal fluctuations can be ruled out. The magneti-
ferromagnetic energy induced order due to the nearest- zation M(t) cannot follow the external field h(t) (τ >P
neighbor exchange coupling through the system at low case) for each selected field amplitude values H0, conse-
temperatures. Hence, a DPT cannot be observed un- quently the dynamic ferromagnetism is enhanced. Re-
less a relatively large amount of thermal energy is sup- duced exchange varies from R=1.0 to 2.0 including the
plied to the system. As the field amplitude increases, critical value of itself. So, we have qualitatively different
it becomes dominant against the ferromagnetic nearest- two regimes: R < R and R > R . Below R the in-
c c c
4
(kBTc/J −R) plane, we evaluated the thermal variation
1.04
of specific heat for a given set of system parameters. A
(a)
typicalexampleisshowninFig. (4)forthefilmswithdif-
> 0.96
ferentthicknesses as L=3,4,5. The temperature values
)
t
M(k0.88 R = 1.0 cthoerretrsapnosnidtiionng tteomthpeermatauxriems.aForfosmpetchifiecphaenaetlc(uar)vteos a(cre)
< 1.1 1.6
1.2 1.7 H = 0.0 and (d) to (f), field amplitude H0 changes. Also, from
0.80 1.3 1.8 =3.25 (a) to (d), (b) to (e) and (c) to (f) reduced exchange
1.4 1.9
1.5 2.0 has two different value. DPT points shift towards the
0.72
lower temperature with increasing field values as well
1.04
as they shift towards higher temperature with increas-
(b)
ing reduced exchange. The related detailed story has
0.96
been explained above. Similarly, the effect of reduced
Q k0.88 exchange on specific heat peaks is easy to understand:
Thestrongerdipole-dipoleinteractionmakemorecontri-
0.80 H0 = 0.25 bution to total energy. Hence, more thermal agitation is
P = 100 needed to make the system dynamically disordered. In
=3.25
0.72 Fig. (4), the crossoverbehavior can be seen easily when
1.04 R changed from R < Rc to R > Rc (namely, from (a)
(c) to (d), (b) to (e) and (c) to (f)). Below Rc, thicker film
0.96 has more nearest-neighbor interaction per site, this cre-
ates more contribution to energy. Consequently, both
k
Q 0.88 strength of the peak and corresponding critical temper-
ature are relatively higher than the others.
H0 = 0.5
0.80
P = 100 In order to obtain a general overview of the non-
=3.25 equilibrium phase diagramin (k T /J −R) planes. For
0.72 B c b
0 5 10 15 20 this purpose,inFig. (5)weplotthe criticaltemperature
k versus R with selected film thickness values L = 3,4,5
FIG. 3. (Color online) Average magnetizations profiles for andforthreeselectedvaluesoffieldamplitudeH0. Since,
a film with L = 20 layers. The number accompanying each the temperature values at which the specific heat curves
curve denotes several values of reduced exchange varies from exhibit a sharp maximum correspond to the transition
1.0 to 2.0 including Rc. Each layerlabelled by k. temperatureofthinfilm,criticaltemperaturevalueshave
beenobtainedbyexaminingthethermalvariationofspe-
cific heat curves (a selected set has been given in Fig.
ner layers are highly ordered compare to surfaces. This (4)). Fig. (5) represents a characteristic phenomenon
canbebrieflyexplainedasfollows: Inmiddle ofthe film, peculiar to thin film systems. Namely, due to the exis-
there are relatively more neighboring per magnetic sites tence of reduced surfaces, there exists a special value of
which causes locally larger magnetic interaction. So it surfacetobulkratioofexchangeinteractionsRcatwhich
becomesmoredifficulttofollowtheexternalfieldforany the transition temperature of the film becomes indepen-
spin. Surface spins are embedded in an environment of dent of thickness L. Simply we can say that the curves
lower symmetry than that of the inner atoms. The ex- with different film thicknesses intersect each other. This
changeconstantbetweenatomsinthesurfaceregionmay fullysupportsarecentstudy[19]whereintheframework
differ from the bulk one. This regime corresponds to a ofanEFT the existence of aspecialtransitionpointwas
surface type of magnetic ordering. The opposite of the predicted. The critical temperature value of crossover
above scenarios can be considered also. Above Rc, in point for H0 = 0.0 (static case) is in a good agreement
the inner layers, although there are more neighboring, with previous studies [19, 22, 28]. However, variation of
there are far fewer exchange constant per magnetic sites Rc asafunctionofH0 isveryslowaccordingtoFig. (5),
whichcausesrelativelysmallermagneticinteractionthan and we see that the location of R barely deviates from
c
that ofthe surface one. Rc plays the main roleto obtain its equilibrium value with increasing H0. This deviation
the frontier of this crossover. The free surface cannot has been reported before by Yu¨ksel [12]. The discussion
break the translational symmetry since magnetic prop- on existence of this kind of deviation is an academic is-
erties of the free surfaces exactly overlap with the bulk sue and this may be due to insufficient data and cannot
oneatR . Wecansaymoregenerallythatthedeficiency be located accurately with such an approach with less
c
of the interaction per surface spin can be compensated effort. This was also reported in a MC simulation treat-
byincreasingthemodifiedexchangeinteractionstrength. ment of surface critical phenomena by Hasenbusch [29].
Moreover,the effect of external field can be also seen by From panel (a) to (c) the effect of field amplitude H0 on
following the panels from (a) to (c). the transition characteristics of the film can be seen and
In obtaining the critical frontiers depicted in it is also straightforward as stated before: Greater the
5

2.5 8
(a) L = 3 (d)
B
2.0 4 6 kBTC/Jb
5
C1.5 RH0= 0=. 50.0 HR=0 2=. 50.0 4 L = 3 H0 = 0.0
4
1.0
2 5
0.5 RC
(a)
0

0.0 8
2.5 H0 = 0.25
(b) H0 = 0.25 (e) H0 = 0.25 6 P = 100
2.0
P = 100 P = 100
R=0.5 R=2.5 4
C1.5

1.0 2
RC
(b)
0.5 0

8
0.0 H0 = 0.5
2.5 (c) (f) 6 P = 100
H0 = 0.5 H0 = 0.5
2.0
P = 100 P = 100 4
R=0.5 R=2.5
C1.5
2
1.0 RC (c)
0
0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0
R
0.0
0 2 4 6 8 0 2 4 6 8 FIG. 5. (Color online) The critical temperature versus re-
ducedexchangeinteractioninphaseplanesforthefilmswith
different L=3,4,5 thicknesses.
FIG. 4. (Color online) The specific heat versus temperature
for the films with different L=3,4,5 thicknesses.
marizedasfollows: We firstinvestigatethe thermalvari-
ations of related order parameters for the films with dif-
amplitude H0 means more energy transferredto the sys- ferent thicknesses. The effect of the field amplitude and
tem over half cycle by the external oscillatory field and
reduced exchange on the previous standard arguments
this makestransitionfromdynamicallyorderedto disor-
for static case has been propounded. In the vicinity
dered phase more easy. For R < R , we have ordinary
c of the dynamic ferromagnetic-paramagnetic phase tran-
transition behavior where the bulk magnetism is domi-
sition temperature, specific heat curves exhibit a sharp
nantagainstthe surface magnetismwhereasfor R>R ,
c peak which becomes more apparent for sufficiently high
the surface may exhibit enhanced magnetic behavior in
reduced exchange values (R > R regime) values. The
c
comparisonwith bulk. This is calledextraordinarytran-
thinnerfilmsintheabsenceofenhancedsurfaces(R<R
c
sition. Moreover, as shown in Fig. (5), for R < R ,
c regime) with high field amplitudes exhibit a weak peak
thicker films have greater transition temperatures while
in relatively lower temperatures.
for R > R , the transition temperature of the film de-
c According to our findings, an increment of the field
creases with increasing thickness. The results indicate
amplitude causes a decreasing in corresponding temper-
that the well-known surface enhancement properties of ature coordinateof the crossoverpoint in (k T /J −R)
B c b
thesystemmaychangeitscharacteristicsinthepresence
planes. Criticalvalueofsurfacetobulkratioofexchange
of an external oscillatory field.
interactionsR atwhichthetransitiontemperatureisin-
c
dependent of film thickness is not apparently responsive
to varying field amplitude values, but exhibits slow vari-
IV. CONCLUSION
ation as a function of H0. We confirmed the general
trend was generated by using EFT calculations before
In conclusion, we have applied MC simulations to [19]. Hence, we can say that the evolution of a crossover
studytheDPTcharacteristicsinthinferromagneticfilms is not from the limitation of EFT.
in the presence of oscillating magnetic fields. The fore- We hope that this study will shed light on further in-
most results obtained from simulation data can be sum- vestigationsofthedynamicnatureofcriticalphenomena
6
in pure crystalline ferromagnetic thin films and will be Performance and Grid Computing Center (TRUBA Re-
beneficial from both theoreticaland experimentalpoints sources) and this study was completed at Dokuz Eyll
of view. University,Graduate Schoolof NaturalandApplied Sci-
ences. One of the authors (B.O.A.) would like to thank
theTurkishEducationalFoundation(TEV)forfullschol-
ACKNOWLEDGEMENTS arship.
The numerical calculations in this paper were per-
formedatTU¨BI˙TAKULAKBI˙M(Turkishagency),High
[1] E. Schierle, E. Weschke, A. Gottberg, W. S¨ollinger, W. [12] Y. Yu¨ksel, Phys. Lett. A 377 (2013) 2494.
Heiss, G. Springholz, and G. Kaindl, Phys. Rev. Lett. [13] U¨. Akıncı,JMMM 368 (2014) 36.
101, (2008) 267202. [14] U¨. Akıncı,Thin Solid Films 550 (2014) 602.
[2] M. Ahlberg, M. Marcellini, A. Taroni, G. Andersson, [15] Y. Yu¨ksel, U¨. Akıncı,Physica B 433 (2014) 96.
M. Wolff, and B. Hjo¨rvarsson, Phys. Rev. B 81, (2010)
[16] Y. Yu¨ksel, Physica A 396 (2014) 9.
214429.
[17] C. Waldfried, T. McAvoy, D. Welipitiya, P. A. Dowben,
[3] C. E. ViolBarbosa, H. L. Meyerheim, E. Jal, J.-M. Ton-
E. Vescovo,Europhys. Lett. 42 (1998) 685.
nerre, M. Przybylski, L. M. Sandratskii, F. Yildiz, U.
[18] R.Skomski,C.Waldfried,P.A.Dowbeny,J.Phys.: Con-
Staub,andJ.Kirschner,Phys.Rev.B85,(2012)184414.
dens. Matter 10 (1998) 5833.
[4] A.Berger,O.Idigoras,andP.Vavassori,Phys.Rev.Lett. [19] B. O. Akta¸s, U¨. Akıncı, H. Polat, Thin Solid Films 562
111, (2013) 190602.
(2014) 680.
[5] B.-Y.Wang,J.-Y.Hong,K.-H.O.Yang,Y.-L.Chan,D.-
[20] O. Canko, E. Kantar, M. Keskin, Physica A 388 (2009)
H. Wei, H.-J. Lin, and M.-T. Lin, Phys. Rev. Lett. 110,
28.
(2013) 117203.
[21] M. Erta¸s, M. Keskin,Chin. Phys.B 22 (2013) 120507.
[6] O. Yal¸cın, S¸. U¨nlu¨er, S. Kazan, M. O¨zdemir, Y. O¨ner.
[22] H. Park, M. Pleimling, Phys. Rev. Lett. 109, (2012)
JMMM 373 (2015) 144.
175703.
[7] M. Saber, A. Ainane, F. Dujardin, B. Stb, Journal of
[23] H. Park, M. Pleimling, Phys. Rev.E 87, (2013) 032145.
Non-Crystalline Solids 250 (1999) 735.
[24] K. Taucher, M. Pleimling, Phys. Rev. E 89, (2014)
[8] A. Oubelkacem, A.Ainanea, J. J. deMiguel, J. Ricardo
022121.
deSousa, M. Saber, Physica A 358 (2005) 160.
[25] K. Binder, Monte Carlo Methods in Statistical Physics,
[9] S. Tuleja, J. Kecer, and V. Ilkovi, Phys. Stat. Sol. (b)
Springer, Berlin, 1979.
243, (2006) 1352.
[26] M.E.J.Newman,G.T.Barkema,MonteCarloMethods
[10] A. Zaim, M. Kerouad, Y. EL Amraoui, D. Baldomir,
in Statistical Physics, Oxford UniversityPress, 2001.
JMMM 316 (2007) e306.
[27] T.Tome,M.J.deOliveira,Phys.Rev.A41,(1990)4251.
[11] U¨.Akıncı, JMMM 329 (2013) 178. [28] T. Kaneyoshi, J. Phys. Condens. Matter 3 (1991) 4497.
[29] M. Hasenbusch,Phys.Rev. B 84, (2011) 134405.