Critical point determination from probability distribution functions in the three dimensional Ising model
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Critical point determination from probability distribution functions in the threedimensional Ising model
Francisco Sastre Departamento de Ingenier´ıa F´ısica, Divisi´on de Ciencias e Ingenier´ıas,Campus Le´on de la Universidad de Guanajuato,AP E-143, CP 37150, Le´on, Guanajuato, M´exico
In this work we propose a new numerical method to evaluate the critical point, the susceptibilitycritical exponent and the correlation length critical exponent of the three dimensional Ising modelwithout external field using an algorithm that evaluates directly the derivative of the logarithm ofthe probability distribution function with respect to the magnetisation. Using standard finite-sizescaling theory we found that correction-to-scaling effects are not present within this approach. Ourresults are in good agreement with previous reported values for the three dimensional Ising model.
PACS numbers:
I. INTRODUCTION
The Ising model has a great importance in statisticalmechanics since a great variety of techniques and meth-ods, analytical and numerical, have been formulated firston this model. There are several numerical algorithmsthat can be used to study the critical behavior in spinsystems, we can mention three types: • Those that require adjusting a control parameter,like the standard Monte Carlo or the Wolff algo-rithms [1]. • Algorithms that not required to fine tune any pa-rameter, examples of those kind are the InvasionCluster Algorithm [2], algorithms based on Self Or-ganisation [3] or the Locally Cluster Algorithm [4]. • Algorithms that evaluate the Density of States inthe microcanonical ensemble [5, 6].In this work we propose a new methodology based onthe algorithm proposed by Sastre et al. for the evaluationof effective temperatures in out of equilibrium Ising-likesystems [7]. The algorithm is a canonical generalisationof the work proposed by H¨uller and Pleimling [6] for theevaluation of the density of states in the two and threedimensional Ising model. It is important to point outthat variations of this algorithm have been successfullyimplemented in fluids with discrete potential interaction.The microcanonical version was used in [8] to evaluatethermodynamic properties in the supercritical region andin [9] for the evaluation of the Hight Temperature Ex-pansion coefficients of the Helmholtz free energy. Thecanonical version [10] was used for the evaluation of thecritical temperature and the correlation length criticalexponent in the Square-Well fluid with interaction rangeof 1 . ν , and the susceptibility, γ , in thethree dimensional Ising model in an efficient way.This article is organized as follows: in section II weexplain the basic definitions for the Ising model and thealgorithm used in this work, in section III we apply themethod to evaluate the critical point and the critical ex-ponents. We made our concluding remarks in the sec-tion IV. II. PROBABILITY DISTRIBUTIONS FOR THEISING MODEL
For a better understanding of the algorithm used inthis work we will review first how the algorithm proposedby H¨uller and Pleimling in the microcanonical ensembleworks. The algorithm uses a variation of the transitionvariable method [11–13] in order to evaluate the entropyas function of the energy and the magnetisation.The hamiltonian for the Ising model on a three di-mensional cubic lattice without external field and nearestneighbors interaction is1 k B T H = − β − X h i,j i σ i σ j , (1)where σ i = ± i th site, β = k B T /J isthe control parameter and
J > h i, j i indicatesthat the summation runs over all nearest neighbors pairson the lattice. If we consider a system with periodicboundary conditions and N = L spins, the magnetisa-tion, M = P i σ i , and the energy will be bounded. More-over, when a spin is flipped in the system we observe that∆ M = ± E/J = ± η , η = 0 , ± , ± , ± E µ , M k )is the number of microstates that share the same mag-netisation M k and energy E µ , for simplicity we will useΩ µ,k for this number. When the systems is in a givenmacrostate Ω µ,k and we flip a spin at random we canreach a new macrostate Ω ν,l , with E ν = E µ + 4 ηJ and M k = M l ±
2. The probability of reach Ω ν,l starting fromΩ µ,k will be given by P ( m )( µ,k ) , ( ν,l ) = V ( µ,k ) , ( ν,l ) N Ω µ,k , (2)the superscript ( m ) denotes that we are working in themicrocanonical ensemble and V ( µ,k ) , ( ν,l ) indicates howmany ways the system can reach Ω ν,l starting from Ω µ,k ,this quantity is purely geometric. We can also obtain thereverse probability with P ( m )( ν,l ) , ( µ,k ) = V ( ν,l ) , ( µ,k ) N Ω ν,l . (3)As any spin flip can be reversed, the relation V ( µ,k ) , ( ν,l ) = V ( ν,l ) , ( µ,k ) must be satisfied. For example, from the basestate with E µ = − N J and M k = N we can reach thestate with E ν = − J ( N −
4) and M l = N − N ways,then V ( µ,k ) , ( ν,l ) = N and, as Ω µ,k = 1, P ( m )( µ,k ) , ( ν,l ) = 1. Inthe reverse process we have Ω ν,l = N then P ( m )( ν,l ) , ( µ,k ) =1 /N . Combining equations (2) and (3) we obtain theimportant microcanonical relation P ( m )( µ,k ) , ( ν,l ) P ( m )( ν,l ) , ( µ,k ) = Ω ν,l Ω µ,k . (4)The last equation can be used to evaluate the micro-canonical derivatives 1 T = ∂S∂E , (5)and − BT = ∂S∂M , (6)where B is the magnetic external field. The change onthe entropy can be obtained using the following approx-imation∆ S = k B ln P ( m )( µ,k ) , ( ν,l ) P ( m )( ν,l ) , ( µ,k ) ≈ ∆ E T − ∆ M BT . (7)In the simulation the probabilities can be estimated withthe rate of attempts T ( µ,k ) , ( ν,l ) to go from Ω µ,k to Ω ν,l .The rate is given by the relation T ( µ,k ) , ( ν,l ) = z ( µ,k ) , ( ν,l ) z ( µ,k ) , (8)where z ( µ,k ) is the number of times that the systemspends in a macrostate Ω µ,k , and z ( µ,k ) , ( ν,l ) is the num-ber of times that the system attempts to change fromΩ µ,k to Ω ν,l . In this method, once that we fix the ranges[ E min , E max ] and [ M min , M max ], where the simulationwill be confined, the quantities z ( µ,k ) , ( ν,l ) and z ( µ,k ) canbe estimated in the following way: 1. With E min ≤ E µ ≤ E max and M min ≤ M k ≤ M max as initial state, a spin is chosen at randomand z ( µ,k ) is always incremented by 1.2. We evaluate the new values E ν and M l that thesystem would take if the chosen spin is flipped.3. If E min ≤ E ν ≤ E max and M min ≤ M l ≤ M max the quantity z ( µ,k ) , ( ν,l ) is incremented by 1.4. The spin flip attempt is accepted with probabilitymin (1 , z ( ν,l ) , ( µ,k ) z ( µ,k ) z ( µ,k ) , ( ν,l ) z ( ν,l ) ). This condition assures thatall macrostates are visited with equal probability,independently of their degeneracy.The values z ( µ,k ) , ( ν,l ) and z ( µ,k ) can be initialized withany positive integer, a safe option is 1, and after alarge number of spin flip attempts we will observe that T ( µ,k ) , ( ν,l ) → P ( µ,k ) , ( ν,l ) .M¨uller an Pleimling used this method for the determi-nation of the microcanonically defined spontaneous mag-netisation and the order parameter critical exponent forthe Ising model in two and three dimensions. This algo-rithm is highly efficient for evaluate the ratios Ω ν,l / Ω µ,k since it gives the freedom of restricting the calculationsto a chosen range in the energy and magnetisation. Addi-tional details of the method can be found in the originalwork [6].The microcanonical algorithm counts all the attemptsto change from a given macrostate to another, as long asthe final macrostate is an allowed one, while the gener-alisation proposed in [7] adds an additional condition tothe attempts count. The additional condition includes a”heat bath”, then we will need to incorporate an extrafactor in the ratio of probabilities P ( µ,k ) , ( ν,l ) P ( ν,l ) , ( µ,k ) = P ( m )( µ,k ) , ( ν,l ) P ( m )( ν,l ) , ( µ,k ) e − ∆ EkBT = Ω ν,l Ω µ,k , e − ∆ EkBT (9)here the absence of the superscript indicates that we areno longer in the microcanonical ensemble. CombiningEquations (7) and (9) we getln( P ( µ,k ) , ( ν,l ) ) − ln( P ( ν,l ) , ( µ,k ) ) ≈ − ∆ M Bk B T , (10)where we can drop the subscript µ and ν , since the righthand side of the equation depends only on ∆ M . In thesimulation the probabilities now can be estimated withthe rate of attempts T kl to go from a macrostate withmagnetisation M k (level k ) to a macrostate with mag-netisation M l (level l ). The quantity T kl will be givennow by the relation T kl = z kl z k , (11)where z kl is the number of times that the system at-tempts to change from level k to level l and z k is thenumber of times that the system spends in level k . Forthe estimation of the z kl and z k values, the detailed stepsare now:1. With M min ≤ M k ≤ M max as initial state, a spinis chosen at random and z k is always incrementedby 1.2. If the possible state M l , that would be reachedif the chosen spin is flipped, is allowed we eval-uate ∆ E between states M k and M l , and thequantity z kl is incremented by 1 with probabilitymin(1 , e − ∆ E/k B T ).3. The spin flip attempt is accepted with probabilitymin (1 , z lk z k z kl z l ).The values z kl and z k are initialized to 1 and after alarge number of spin flip attempts we will observe that T kl → P kl .Now we obtain the derivative of ln P , instead of thederivatives of S , with the following approach ∂ ln P∂M (cid:12)(cid:12)(cid:12)(cid:12) k ≈ ln P kl − ln P lk ≈ ln( T kl /T lk ) . (12)We will use the following function in our simulations g ( m ) = L − d ∂ ln P∂m , (13)where m = M/L d , since it can be obtained directly fromthe transition rates with the approximation g ( m k ) ≈
12 (ln[( T k,k +1 ) / ( T k +1 ,k )] − ln[( T k,k − ) / ( T k − ,k )]) , (14)where T k,k ± k are the transition rates from level k to itsadjacent levels and T k ± ,k are the transition rates fromthe adjacent levels to level k .We verified that this method is compatible with theansatz proposed by Tsypin and Bl¨ote[14] for the threedimensional Ising model at the critical point P ( m ) ∼ exp − (cid:18) mm (cid:19) − ! a (cid:18) mm (cid:19) + c ! , (15)where a, c and m are size depending fitting parameters.For L = 12 we performed three independent simulationsat β c = 0 . m . Our results are show in Figure 1 alongwith the curve given by Eq. (15), using a = 0 . c =0 . m = 0 . III. RESULTS
Now we will explain how to obtain the critical temper-ature and the correlation length critical exponent. It isa well known fact that there is a change in the probabil-ity distribution of the order parameter above and below a -0.4 -0.2 0 0.2 0.4 m -0.02-0.0100.010.02 g ( m ) |m|<0.694|m|<0.347|m|<0.174 Tsypin and Blöte m -0.04-0.02 g FIG. 1: (Color online) Graph of g as function of m forthree different range simulations. We observe that the re-striction in range does not affects the computed value of g .The continuous line reproduces Eq. (15) with the parameters a = 0 . , c = 0 .
859 and m = 0 . L = 12. There is areally good agreement with our results, except in the extremesof the curves, as shown in the inset. certain value of β , denoted as β c ( L ), that depends on thesystem size. For β values below β c ( L ) there is just onepeek at m = 0, that means that the curve g ( m ) crossesthe horizontal axis with a negative slope. For β valuesabove β c ( L ) there are two peeks at m = ± m sta in theprobability distribution, while g ( m ) crosses the horizon-tal axis in three places, at m = ± m sta with a negativeslope and at m = 0 with a positive slope. In Figure 2 wecan see clearly the symmetry breaking for three dimen-sional Ising model with linear size L = 8.As we can see there is a change in the sign of the slope ∂g/∂m around m = 0 as function of β , then we can find β c ( L ) restricting our simulations around m = 0. Weperformed simulations restricting the intervals to | m | . . g ( m ) is linear around m = 0 and the slopecan be easily obtained from a linear fit to the simulationdata. For every system size we evaluate the slope forseveral values of β , in Figure 3 we are illustrating howthe evaluation of β c ( L ) is performed for the case L = 12.Here we used a linear fit to the curves of the slope as -0.5 0 0.5 m -0.2-0.100.10.2 g ( m ) β < β c β > β c FIG. 2: (Color online) Symmetry breaking for the three di-mensional Ising model with L = 8. For β < β c the curvecross the horizontal axis at m = 0 with a negative slope. For β > β c we have three crossings, at m = 0 and at m = ± m sta and we observe that the slope at m = 0 is now positive. function of β in order to solve (cid:16) ∂g∂m (cid:17) β = β c = 0.From here we can use the Finite size scaling ansatz forthe critical temperature, given by the relation β c ( L ) ≈ β c + AL − /ν , (16)where β c is the critical control parameter for the infinitesystem, A is a non universal parameter and ν is the crit-ical exponent for the correlation length [16]. In principleEq. (16) is valid for sufficiently large L values. For smallsystems the last term changes to L − /ν (1 + BL − ω ), herethe parameter ω is the scale correction exponent, whosereported value for the three dimensional Ising model is ω ≈ .
81 [17].The simulations were carried out in systems with lin-ear sizes L = 8, 10, 12, 14, 16, 20 and 24, using 5 N × spin flip attempts and 120 independent runs for every setof parameters. With these values we obtain reliable datawith less CPU time compared with standard canonicalsimulations, since most of the spin flips are discardedwhen the system falls outside of the restricted range.However, all attempts, successful or not, are used forthe evaluation of the transition rates. We evaluated T c and ν performing a non-linear curve fitting to Eq. (16),in Figure 4 we show the evaluation of the critical point β -0.008-0.0040.00.0040.008 ∂ g___ ∂ m FIG. 3: Slopes of the function g ( m ) at m = 0 for the threedimensional Ising model with L = 12 as function of the controlparameter β . The dashed line is a linear fit to the simulationdata. along with the ν critical exponent. We must empha-sized that in our analysis the scale correction exponentis absent, which is a great advantage in numerical sim-ulations that study critical phenomena. This feature isalso observed in the evaluation of the critical temperaturefor the square-well fluid using an equivalent method [10].We think that the absence of scale corrections are relatedto the fact that in this method we evaluate the criticalpoint analyzing the behavior of the probability distribu-tion of the order parameter around M = 0, while in mosttraditional methods the critical point is evaluated ana-lyzing the behavior around the peaks of the probabilitydistribution.The results for the critical point is β c = 0 . ν = 0 . γ .We used the scaling ansatz for the probability distribu-tion function at the critical point proposed in Ref. [20] P ( m, L ) ≈ exp ( − A + A x + A x + . . . ) , (17) L −1 β c ( L ) FIG. 4: Evaluation of the critical temperature and the corre-lation length critical exponent for the three dimensional Isingmodel. The dashed line is a non-linear curve fit to Eq. (16).The results from the fit are β c = 0 . A = − . /ν = 1 . β c ν Source0.22165(65) 0.6301(88) This work0.22165452(8) 0.63020(12) Butera and Comi [18]0.221655(2) 0.6299(2) Deng and Bl¨ote [19]0.221654(2) 0.6308(4) Lundow and Campbell [17] where x = mL β/ν << β is the order parametercritical exponent. As we are restricting our simulationsto the range | m | ≤ . ∂ ln P∂m = 2 A mL β/ν + O ( x ) , (18)that we can combine with Eq. (13) to obtain the desiredscaling relation ∂g∂m ∼ L − γ/ν . (19)For the three dimensional case we performed simula-tions at our estimated critical point β c = 0 . L = 8, 10, 12, 14, 16, 20 and 24. In this case we
10 20 L ∂ g ___∂ m T=T c ||| FIG. 5: Evaluation of the critical exponent γ/ν for three di-mensional Ising model. We are showing a log-log graph ofthe derivatives ∂g/∂m as function of the linear size L . Thedashed line is the linear fit to Equation (19). From the fit weobtain the slope γ/ν = 1 . used N × spin flip attempts and 120 independent runsfor every L value. The critical exponent γ was obtainedfrom a linear fit to Eq. (19) as shown in Figure 5. Fromthis fit we obtain γ/ν = 1 . γ/ν = 1 . IV. CONCLUSIONS
We have presented a new method for the evaluation ofthe critical temperature and the critical exponents for thecorrelation length ν and the susceptibility γ on the threedimensional Ising model. Using the derivatives of theprobability distribution function for the magnetisationand small system sizes we obtain reliable results that arein good agreement with the reported values. The methodcan be used in a restricted range on the magnetisationand this feature reduces the computational time in thesimulations. One additional advantage of the method isthat scale corrections are not present, at least in the threedimensional Ising model. In future works we will studyif this advantage is present in other systems. V. ACKNOWLEDGEMENTS
This research was supported by Universidad de Gua-najuato (M´exico) under Proyecto DAIP 879/2016 andCONACyT (M´exico)(grant CB-2017-2018-A1-S-30736- F-2164).
References [1] U. Wolff, Phys. Rev. Lett. , 361 (1989).[2] J. Machta, Y. S. Choi, A. Lucke, T. Schweizer, and L. V.Chayes, Phys. Rev. Lett. , 2792 (1995).[3] Physica A: Statistical Mechanics and its Applications , 171 (1999), ISSN 0378-4371.[4] E. Faraggi and D. T. Robb, Phys. Rev. B , 134416(2008).[5] F. Wang and D. P. Landau, Phys. Rev. Lett. , 2050(2001).[6] A. H¨uller and M. Pleimling, International Journal ofModern Physics C , 947 (2002).[7] F. Sastre, I. Dornic, and H. Chat´e, Phys. Rev. Lett. ,267205 (2003).[8] F. Sastre, A. L. Benavides, J. Torres-Arenas, and A. Gil-Villegas, Phys. Rev. E , 033303 (2015).[9] F. Sastre, E. Moreno-Hilario, M. G. Sotelo-Serna, andA. Gil-Villegas, Molecular Physics , 351 (2018).[10] F. Sastre, Molecular Physics , e1593534 (2020).[11] The European Physical Journal B - Condensed Matter and Complex Systems (1998), ISSN 1434-6028.[12] The European Physical Journal B - Condensed Matterand Complex Systems (1998), ISSN 1434-6028.[13] M. Kastner, M. Promberger, and J. D. Mu˜noz, Phys.Rev. E , 7422 (2000).[14] M. M. Tsypin and H. W. J. Bl¨ote, Phys. Rev. E , 73(2000).[15] I. A. Campbell and P. H. Lundow, Phys. Rev. B ,014411 (2011).[16] A. M. Ferrenberg and D. P. Landau, Phys. Rev. B ,5081 (1991).[17] P. H. Lundow and I. A. Campbell, Phys. Rev. B ,024414 (2010).[18] P. Butera and M. Comi, Phys. Rev. B , 144431 (2002).[19] Y. Deng and H. W. J. Bl¨ote, Phys. Rev. E , 036125(2003).[20] K. Kaski, K. Binder, and J. D. Gunton, Phys. Rev. B29