An entropic simulational study of the spin- 1 Baxter-Wu model in a crystal field
L. N. Jorge, P. H. L. Martins, C. J. DaSilva, L. S. Ferreira, A. A. Caparica
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug An entropic simulational study of the spin- Baxter-Wu model in a crystal field
L. N. Jorge
Instituto Federal de Educa¸c˜ao, Ciˆencia e Tecnologia do Mato Grosso Campus C´aceres - Prof. Oleg´ario Baldo,CEP 78201-380, C´aceres, Mato Grosso. Brazil ∗ P. H. L. Martins
Instituto de F´ısica, Universidade Federal de Mato Grosso, CEP 78060-900, Cuiab´a, Mato Grosso, Brazil.
Claudio J. DaSilva
Instituto Federal de Educa¸c˜ao, Ciˆencia e Tecnologia de Goi´as, Rua 76, Centro, Goiˆania - GO, Brazil.
L. S. Ferreira and A. A. Caparica
Instituto de F´ısica, Universidade Federal de Goi´as,C.P. 131 CEP 74001-970, Goiˆania, Goi´as, Brazil.
We investigate the critical behavior of the two-dimensional spin-1 Baxter-Wu model in a crystalfield using entropic sampling simulations with the joint density of states. We obtain the temperature-crystal field phase diagram, which includes a tetracritical line ending at a pentacritical point. Afinite-size scaling analysis of the maximum of the specific heat, while changing the crystal fieldanisotropy, is used to obtain a precise location of the pentacritical point. Our results give thecritical temperature and crystal field as T pc = 0 . D pc = 1 . I. INTRODUCTION
The spin-1 Baxter-Wu (BW) model in a crystal field[1–3] is a generalization of the original spin − BW model[4–7], which includes a crystal field anisotropic term D , inaddition to the three-spin interaction. The Hamiltonianof the model considered here is H = − J X h i,j,k i s i s j s k + D N X i s i , (1)where the spin variables are located at the sites of a tri-angular lattice and assume the values s i = ± , J > D is the anisotropy due to the crystallinefield. The first sum extends over all triangular faces whilethe second runs over all lattice sites. Note that, for thespin- case, the second term in Eq. 1 is just a constant.The Hamiltonian in Eq. 1 resembles the Blume-Capel(BC) model[8, 9] case, which is described by a phase dia-gram that presents a line of continuous phase transitionsmeeting a discontinuous one at a tricritical point.In the thermodynamic limit, at zero temperature andfor D <
2, four ordered states are present in the groundstate, one ferromagnetic, (+ + +), and other three ferri-magnetic phases: (+ − − ), ( − − +) and ( − + − ), so thatthe spin-1 / D >
2, thestate that minimizes the energy is the microstate (000),with all spins s i = 0. For D = 2, one has a first-orderphase transition with all those five phases coexisting. ∗ [email protected] The first proposal to study the BW model with in-teraction anisotropy was made by Kinzel et al. [1], ana-lyzing the correlation length behavior, the authors con-cluded that a continuous transition only occurs for thepure case, s = 1 / D c = 1 . T c = 1 . et al. [3] by applying finite-size scaling and conformal invariance, analyzed both thespin-1 and spin-3 / D c = 0 . T c = 1 . D = 0 was studied using entropic simu-lations, resulting in the observation that the BW modelpresents a mixture of phase transitions[10], related to atetracritical point. That was settled after the analysisof the energy probability, the configurations in the criti-cal region, and by separating the density of states in twoparts, sorting out the ferromagnetic and ferrimagnetic ar-rangements. Besides, this previous study analyzed howthe choice of lattice sizes affects the results on the ther-modynamic limit for the critical exponents and criticaltemperature, and it is clear now that choosing multipleof three lattice sizes or not is equivalent.Entropic simulations[11–14] are excellent to studyphase transitions and critical phenomena. The resultsobtained by this technique have revealed important char-acteristics regarding different models[10, 15–17]. The es-timation of the joint density of states[14, 17, 18] brings arange of additional information, besides making it possi-ble to obtain thermodynamic properties for any temper-ature and coupling constants, allowing the constructionof the phase diagram.In the present work, we completely characterize thephase diagram and determine the position of the multi-critical point by studying the model using the joint den-sity of states. We also analyze the possibility that theline, said to be of second-order, might be a line of co-existence of four configurations ending in a pentacriticalpoint. Finally, we perform a study of the specific heatin the region of the first-order transition to understandwhat causes this transition.This work is organized as follows: in the next section,we describe the computational details. In section III, wepresent the results and the phase diagram, and we stillreport an anomaly that we found out in the specific heat.The last section comprises the conclusions and some fur-ther remarks. II. ENTROPIC SAMPLING SIMULATIONSWITH THE JOINT DENSITY OF STATES
The study of the effects of spin anisotropy in a mag-netic system has arisen interest in many situations[19–21]. An appropriate approach to tackle these problems isto construct a joint density of states[22], with a secondparameter that characterizes the model. To investigatedimers, some authors have used a joint density of stateswith three-variables dependence[23], g ( N , U, N ), where N is the number of the energetic dimers, U is the energybetween interacting dimers, and N is the total numberof dimers.In the present case, entropic simulations with a jointdensity of states can be performed by defining E = X h i,j,k i s i s j s k , (2)and E = X i s i . (3)Hence the total energy assumes the form E = − JE + DE . The probability of flipping from a configurationwith ( E , E ), to a new one, with ( E ′ , E ′ ), is given by p [( E , E ) → ( E ′ , E ′ )] = min (cid:18) g ( E , E ) g ( E ′ , E ′ ) , (cid:19) . (4)After each Monte Carlo step (MCS), the joint den-sity of states and the histogram are updated asln g ( E , E ) → ln g ( E , E ) + ln f , where f is the mod-ification factor, initially set as f = e = 2 . ... and H ( E , E ) → H ( E , E ) + 1, respectively[12]. It is im-portant to mention that each MCS consists of flipping L spins, whether the attempt is accepted or not. This procedure continues until the histogram reaches its flat-ness – i.e., each term of the histogram is at least 80%of its average over all the possible values of E and E .Then the modification factor is updated to f i +1 = √ f i and the histogram is reset. Furthermore, to halt the sim-ulations, we use the parameter ǫ = | T c ( t ) − T c (0) | [13].To calculate ǫ , we use the following procedure: begin-ning from the modification factor f , when the flatnesscriterion is satisfied, the temperature of the peak of thespecific heat T c (0) is calculated, using the current den-sity of states at the end of each Wang-Landau (WL) level,and whenever the flatness of the histogram is verified, itsvalue is updated. If ǫ remains less than 10 − and the his-togram satisfies the flatness condition for the same levelof WL, then the microcanonical averages and the den-sity of states are saved and the simulation is finished.For all lattice sizes, simulations started from a single rununtil the Wang-Landau level f , since up to that pointthe current density of states is not yet biased and canproceed to any final result. This procedure saves about60% of CPU time, in comparison with the case of start-ing from the first Wang-Landau level. Ref. 16 appliedthis approach with great success.According to the canonical ensemble, the partitionfunction can be written as Z ( T, D ) = X E ,E g ( E , E ) e ( JE − DE ) /k B T , (5)where g ( E , E ) is the joint density of states. So, foreach lattice size, we set up g ( E , E ) in order to obtainthe thermodynamic quantities for any temperature andany value of the crystal field D .Since the function g ( E , E ) does not depend on tem-perature or the crystal field, we can calculate any thermo-dynamic quantity without performing a new simulationrun.The canonical averages of any thermodynamic quan-tity A can be calculated as h A i T,D = P E ,E h A i E ,E P ( E , E , J, D ) P E ,E P ( E , E , J, D ) , (6)where h A i E ,E is the microcanonical average accumu-lated during each simulation run and P ( E , E , J, D ) = g ( E , E ) e ( JE − DE ) /k B T , (7)is the canonical distribution.Performing entropic simulations with a joint densityof states requires huge computational effort, so we con-sidered here only small lattice sizes, yielding excellent re-sults even though. Ref. 24 demonstrated that when usingthe order parameter as the total magnetization, one canconsider systems with non-multiple of three lattice sizeswithout loss of generality of the model [10]. Thus we ransimulations for L = 8, 10, 14, and 16 with n = 24, 20,20, and 16 independent runs, respectively. χ T L=8L=10L=14L=16
FIG. 1. Magnetic susceptibility for D = 1 . III. PHASE DIAGRAMA. Determination of the pentacritical point
In this work, we determine the critical line of the D × T phase diagram using the maximum of the magnetic sus-ceptibility for the largest lattice size taken in our simu-lations, since, as Fig.1 shows, all maxima are very closeto each other, around the same temperature, thereby afinite-size scaling behavior is not visible to the eyes. Infact, finite-size scaling effects for the critical tempera-ture are very subtle when one uses lattice sizes that arenot multiples of three[10, 16, 24]. Thus, the differencebetween the critical temperature at the thermodynamiclimit and those obtained for the finite sizes is small, mak-ing the critical line observed in the diagram very close tothat expected for an infinite system. Hence, for each D value there is a temperature T c that corresponds to themaximum of the magnetic susceptibility.The zero crystal field case, D = 0, was addressed inreferences 25 and 10. In the first work, the authors per-formed MC simulations using the Metropolis algorithmand found a critical temperature of T c = 1 . T c =1 . T = 1 . -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 (E +E )L -2 FIG. 2. Energy probability distribution as a function of( E + E ) /L for different values of D , for L = 16. Thetemperatures where the two peaks occur are shown at theright side of the graph. inverse of the microcanonical temperature β = dS/dE ,and a discontinuity in the order parameter around thetransition temperature. For the ferrimagnetic ones, acontinuous transition, with only one peak in P(E), onlyone inflection point in β , and the order parameter witha typical second-order transition behavior. The analysisof the spin configurations showed the coexistence of fer-romagnetic, ferrimagnetic, and paramagnetic clusters indifferent configurations, revealing a tetracritical behav-ior.To find out if the same behavior persists for D = 0, weanalyze the energy probability distribution for other val-ues of the crystal field, shown in Fig 2. The temperaturesat which the two peaks have the same height are shownon the right. The first curve represents the energy prob-ability for D = − .
0, where we see two peaks. However,this double peak characteristic can be a finite-size effectand, a more detailed study in this region is still miss-ing. For positive values of the field, the finite-size effectwill not remove the two peaks, since we have for the case D = 0 . D pc = 1.68288(62) D L -2 c v D L=8L=10L=14L=16
FIG. 3. Size dependence of the locations of crystal field func-tion at the extrema in the specific heat for lattice sizes L = 8,10, 14, and 16. The inset shows behavior of the maximumheight of the specific heat while increasing the crystal field.Its maximum value is achieved at the transition to the para-magnetic phase. behavior the specific heat maximum, C v max ( D ). The evo-lution of the specific heat maximum versus D results in apeak, which is associated to the beginning of a first orderphase transition. This behavior is presented in the insetof Fig 3. The critical crystal field is obtained from thescaling law[21] D c ( L ) = D ∞ c + bL − d . (8)Figure 3 shows the best fit for d = 2, which gives D c =1 . T c ( L ) = T ∞ c + a L − d , (9)where T ∞ c is the critical temperature at the thermody-namic limit and a is a quantity-dependent constant, en-abling the determination of T c by extrapolating L → ∞ ( L − d = 0) from the linear fitting given by the tempera-ture of the maxima of the specific heat and the suscepti-bility.Figure 4 shows the linear fitting of the temperatureof the maxima of specific heat and susceptibility against L − , for the critical crystal field, which converges to T c when L − →
0. The value for the temperature of thepentacritical point was T c = 0 . T pc = 0.98030(10) T L -2 χ C v U M T L=8L=10L=14L=16
FIG. 4. Temperature of the maxima of specific heat and sus-ceptibility versus L − . The error bars are smaller than sym-bols. The inset shows the cumulant of the magnetization forthe simulated lattice sizes. The crossing of the cumulantstakes place in a region very close to that estimated for thepentacritical point. D c T c Costa[2] 1 . . . . . . ferro-ferri phases from disorder. The dashed line repre-sents the points of coexistence of the five phases: the fourfundamental states and the (0 0 0). The point D = 2 isa fivefold one, where the five states are equally probable.For D > s i = 0 remains, sinceit corresponds to the minimum of the free energy. Thelocation of the pentacritical point is described in TableI, together with some results from previous works, forcomparison. We have obtained D c = 1 . T c = 0 . FIG. 5. Phase diagram for the Baxter-Wu spin-1 model withcrystal field anisotropy in the D versus T plane. The penta-critical point separates the first order phase transition regionand criticality. transitions has been complemented previous works andallows authors to review their results. To the best of ourknowledge, this duality of the BW model was unknownuntil now. B. Anomaly of the specific heat for large crystalfield in the first order region
By analyzing the specific heat, we detected a un-usual behavior for large values of the crystal field, above D = 1 . T = 0 .
80, but it becomes more evident asthe value of the crystal field increases (in fact, this sec-ond peak does not increase, but the first peak decreasesas one approaches D = 2).In Figs. 6a and 6b, the anomaly in C v and the densityof spins are shown for D = 1 . C v and the density of spins for D = 1 . T = 0 .
80. To understandwhich transition is occurring for each peak in the specificheat, we analyzed the density of spins ± n = 1 − n ± , goes inthe opposite direction, drastically increasing in the regionclose to the pseudo-critical temperature. In fact, thisincrease in the number of spins 0 after the pseudo-criticaltemperature manifests itself in a very moderate way forvalues of D >
1, and gradually increases while increasing
D=1.990 (a) c v T L=16L=14L=100.00.20.40.60.81.00.3 0.4 0.5 0.6 0.7 0.8 0.9
D=1.990 (b) n T L=16L=14L=10
FIG. 6. (a) Double peak structure observed for specific heatat D = 1 . D = 1 . T = 0 . the crystal field.This anomaly of the specific heat can be related tothe presence of vacancies, or the Schottky defect, in acrystalline lattice [29, 30], if the non-magnetic spins areinterpreted as vacancies. The peak that presents finite-size effects, in general, is related to the order-disordertransition of the system. However, it can be seen in Figs.6b and 7b, that after the transition, the number of spins ± T → ∞ , the number of non-magnetic spins goes to 1 / ± /
3. It can alsobe noted that the crossing of the ± L = 16 and for crystal field D > . L = 8) the anomaly is already visible and fora smaller value of the crystal field D > . D=1.997 (a) c v T L=16L=14L=100.00.20.40.60.81.00.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
D=1.997 (b) n T L=16L=14L=10
FIG. 7. (a) Double peak structure observed for specific heatat D = 1 . D = 1 . T = 0 . IV. CONCLUSIONS
In this work, we have constructed the temperature-crystal field phase diagram of the Baxter-Wu model us-ing entropic sampling simulations. The simulations wereperformed with a joint density of states, g ( E , E ), sothat with just one simulation run one calculates ther-modynamic quantities for any temperature and crystalfield values. To determine the point at which we haveonly discontinuous transitions, we observe the evolutionof the maximum of the specific heat as a function of thecrystal field. A peak is expected to appear for a cer-tain value of the crystal field, from which the region ofdiscontinuous transitions begins. This point was esti-mated as D pc = 1 . T pc = 0 . D c = 1 . T c = 1 . D c = 0 . T c = 1 . . s i = 0, which canbe associated with the Schottky defect in a crystallinelattice. V. ACKNOWLEDGMENTS
We acknowledge the computer resources provided byLCC-UFG and IF-UFMT. L. S. F. acknowledges the sup-port by CAPES, and L. N. Jorge acknowledges the sup-port by FAPEG. [1] W. Kinzel, E. Domany, and A. Aharony,J. Phys. A: Math. Gen. , L417 (1981).[2] M. L. M. Costa, J. C. Xavier, and J. A. Plascak,Phys. Rev. B , 104103 (2004).[3] D. A. Dias, J. C. Xavier, and J. A. Plascak,Phys. Rev. E , 012103 (2017).[4] D. Wood and H. Griffiths, J. Phys. C. , L253 (1972).[5] R. J. Baxter and F. Wu, Phys. Rev. Lett. , 1294 (1973).[6] R. Baxter and F. Wu, Aust. J. Phys. , 357 (1974).[7] R. Baxter, Aust. J. Phys. , 369 (1974).[8] M. Blume, Phys. Rev. , 517 (1966).[9] H. W. Capel, Physica (Amsterdam) , 966 (1966).[10] L. N. Jorge, L. S. Ferreira, and A. A. Caparica, PhysicaA (Amsterdam) , 123417 (2020).[11] F. Wang and D. Landau, Phys. Rev. E , 056101 (2001).[12] A. A. Caparica and A. G. Cunha-Netto,Phys. Rev. E , 046702 (2012).[13] A. A. Caparica, Phys. Rev. E , 043301 (2014).[14] L. Ferreira, L. Jorge, S. Le˜ao, and A. Caparica, J. Com-put. Phys. , 130 (2018).[15] A. A. Caparica and C. J. DaSilva,Braz. J. Phys. , 713 (2015). [16] L. Jorge, L. Ferreira, and A. Caparica, Phys. Rev. E , 032141 (2019).[17] L. S. Ferreira, L. N. Jorge, C. J. DaSilva, and A. A.Caparica, “Thoroughly analysis of the phase diagram forthe Bell-Lavis model: An entropic simulational study,”(2020), arXiv:2008.01028 [cond-mat.stat-mech].[18] C. Zhou, T. C. Schulthess, S. Torbr¨ugge, and D. P. Lan-dau, Phys. Rev. Lett. , 120201 (2006).[19] J. Plascak and D. Landau, Phys. Rev. E , 015103(2003).[20] L. Bahmad, A. Benyoussef, and A. El Kenz, Phys. Rev.B , 094412 (2007).[21] J. Zierenberg, N. G. Fytas, and W. Janke, Phys. Rev. E , 032126 (2015).[22] W. Kwak, J. Jeong, J. Lee, and D. Kim, Phys. Rev. E , 022134 (2015).[23] L. S. Ferreira, ´A. A. Caparica, L. N. Jorge, and M. A.Neto, Chem. Phys. , 119 (2019).[24] L. N. Jorge, L. S. Ferreira, S. A. Le˜ao, and A. A. Ca-parica, Braz. J. Phys. , 556 (2016).[25] M. Costa and J. Plascak, Braz. J. Phys. , 419 (2004).[26] C. M. Care, J. Phys. A: Math. Gen. , 1481 (1993). [27] M. E. Fisher and A. N. Berker, Phys. Rev. B , 2507(1982).[28] A. Yamagata, J. Phys. A: Math. Gen. , 519 (1993). [29] C. Kittel, Introduction to solid state physics (Wiley,2005).[30] A. W, N. Ashcroft, N. Mermin, N. Mermin, and B. P.Company,