Critical behavior of the classical spin-1 Ising model: a combined low-temperature series expansion and Metropolis Monte Carlo analysis
CCritical behavior of the classical spin-1 Ising model: a combined low-temperatureseries expansion and Metropolis Monte Carlo analysis
Amir Taheridehkordi and Roberto Zivieri ∗ Department of Physics and Physical Oceanography, Memorial University of Newfoundland,St. John’s, Newfoundland & Labrador A1B 3X7, Canada Department of Mathematical and Computer Sciences,Physical Sciences and Earth Sciences, University of Messina, Messina, Italy (Dated: July 20, 2020)In this paper, we theoretically study the critical properties of the classical spin-1 Ising modelusing two approaches: 1) the analytical low-temperature series expansion and 2) the numericalMetropolis Monte Carlo technique. Within this analysis, we discuss the critical behavior of one-,two- and three-dimensional systems modeled by the first-neighbor spin-1 Ising model for differenttypes of exchange interactions. The comparison of the results obtained according the MetropolisMonte Carlo simulations allows us to highlight the limits of the widely used mean-field theoryapproach. We also show, via a simple transformation, that for the special case where the bilinearand bicubic terms are set equal to zero in the Hamiltonian the partition function of the spin-1 Isingmodel can be reduced to that of the spin-1/2 Ising model with temperature dependent externalfield and temperature independent exchange interaction times an exponential factor depending onthe other terms of the Hamiltonian and confirm this result numerically by using the MetropolisMonte Carlo simulation. Finally, we investigate the dependence of the critical temperature on thestrength of long-range interactions included in the Ising Hamiltonian comparing it with that of thefirst-neighbor spin-1/2 Ising model.
I. INTRODUCTION
The classical spin-1/2 Ising model with nearest neigh-bor interactions for a lattice with N sites was sug-gested by Lenz in 1920’s is defined by the followingHamiltonian H spin1 / = − J (cid:88)
1, thesystem will assume another configuration, which musthave higher free energy with respect to ground stateconfiguration C in order to undertake a phase transi-tion at T (cid:54) = 0. We denote these two possible con-figurations by C = (+ + ... + + − + + ... + +) and C (cid:48) = (+ + ... + + 0 + + ... + +) where − and 0 meanspin down and spin-less respectively. The free energies, F and F (cid:48) , associated with these configurations are F = E + 4 J + 4 L + 2 H − k B T ln N, (5) F (cid:48) = E + 2 J + 4 L + H − k B T ln N. (6) In both cases, in the thermodynamic limit and for T (cid:54) = 0we get F , F (cid:48) < F . Thus, the 1D version of spin-1 Isingmodel does not exhibit any phase transition at non-zerotemperatures. However, for higher dimensions the sys-tem described by the spin-1 Ising model exhibits a phasetransition and, due to its enlarged parameter space, amuch richer variety of critical behavior with respect tothe spin-1/2 counterpart . Since the study of the 2D and3D spin-1 Ising model in its general form using either nu-merical or analytical tools becomes very interesting butat the same time very complicated, its critical behav-ior has only been studied in a few special cases in theliterature. In one case just J and H are assumed dif-ferent from zero and the system described by this spin-1Ising model has been solved by applying low-temperatureseries expansion . In another case the spin-1 Isingmodel with non-vanishing J , D , and K has been inves-tigated by using mean-field approximation . The latterone has been studied for K = 0 by other methods likeseries expansion , renormalization group theory , andMonte Carlo simulation .The paper is organized as follows. In section 2 weshortly introduce the analytical and numerical methodsapplied to the spin-1 Ising model. In Section 3 we firstlystudy the critical behavior of the classical spin-1 Isingmodel with nearest neighbor interaction in 2D square and3D cubic lattices, respectively comparing some criticalfeatures with those of the spin-1/2 Ising model, then lastpart of section 3 deals with the long-range spin-1 Isingmodel and the dependence of the critical temperature onthe strength of long-range interactions. Finally, section4 is devoted to the conclusions. II. METHODS
In this section we briefly discuss the analytical andnumerical approaches we have used to analyze the clas-sical spin-1 Ising model. In the first and second sub-sections we describe the mean-field theory and the low-temperature series expansion methods as representativeanalytical methods. The third subsection is a short de-scription of the Metropolis Monte Carlo simulation whichis a strong numerical tool that enables us to investigatethe critical properties of the system.
A. Mean-field theory
We use a systematic way of deriving the mean-fieldtheory for some Hamiltonian H in arbitrary dimensionand coordinate number z (cf. where Mean-field theoriesare studied). We begin from the Bogoliubov inequality F ≤ Φ = F + < H − H > , (7)where F is the true free energy of the system, H is atrial Hamiltonian depending on some variational param-eters which we will introduce, F is the correspondingfree energy, and < ... > denotes an average taken inthe ensemble defined by H . The mean-field free energy, F MF , is then defined by minimizing Φ with respect to thevariational parameters. In order to see how this methodworks let us consider the most general form of the clas-sical spin-1 Ising model given by (2). We introduce thetrial Hamiltonian as the following H = − H N (cid:88) i =1 s i − D N (cid:88) i =1 s i . (8)Hence, there are two variational parameters, H and D , which will be determined by minimizing the func-tional Φ. Assuming that the lattice is translationallyinvariant it is straightforward to find the partition func-tion and the mean-field free energy. Consequently onecan easily compute the magnetization, M = < s i > de-fined in dimensionless units as the thermal average of thespin variable and the thermal average of the square of thespin variable τ = < s i > as follows M = 2 e β ( D + LzM + Kzτ ) sinh[ β ( JzM + H + Lzτ )]1 + 2 e β ( D + LzM + Kzτ ) cosh[ β ( JzM + H + Lzτ )] , (9) τ = 2 e β ( D + LzM + Kzτ ) cosh[ β ( JzM + H + Lzτ )]1 + 2 e β ( D + LzM + Kzτ ) cosh[ β ( JzM + H + Lzτ )] . (10)Solving two self consistent equations (9) and (10) si-multaneously, one can find the mean-field magnetizationas a function of the temperature.We remind that, although the use of mean-field theorygives us some valuable information about the behavior ofthe system and specifically it allows us to reproduce thephase diagram of the model in a relatively simple andqualitative way, from a quantitative point of view theresults are not generally exact for the case of 2D and 3Dlattices so that in this case we need to use more preciseanalytical and numerical methods allowing us to have amore exact knowledge of the critical behavior of thesesystems also quantitatively. B. low-temperature series expansion method
Another possible way to investigate the spin-1 Isingmodel is to use low-temperature series expansion. Theidea is to start from a completely ordered configuration,i.e. the ground state, and then flip spins one by one, andtake all the configurations into account to compute thepartition function, Z as the following: Z = e − E kBT (1 + ∞ (cid:88) n =1 ∆ Z ( n ) N ) , (11)where E is ground state energy, and ∆ Z ( n ) N is the sum ofBoltzmann factors with energy that is measured with re-spect to the ground state energy when n spins are flipped Table I. The counts and the Boltzmann weights contributingto the to the low-temperature series expansion of the partitionfunction of a square lattice classical spin-1 Ising model fordifferent numbners of flipped spins ( N f ). N f Count Boltzmann weight1
N x y z N x yz uw N x y z N x y z uw N x y z u w N ( N − / x y z N ( N − x y z uw N ( N − / x y z u w N x y z N x y z uw N x y z uw N x y z u w N x y z u w N x y z u w N x y z N x y z uw N x y z uw N x y z u w N x y z u w N x y z u w N ( N − x y z N ( N − x y z uw N ( N − x y z uw N ( N − x y z u w N ( N − x y z u w N ( N − x y z u w N ( N − N + 62) / x y z N ( N − N + 62) / x y z uw N ( N − N + 62) / x y z u w N ( N − N + 62) / x y z u w N x y z u w N x y z u w starting from the ground state configuration. Two factorscontribute to the Boltzmann factor ∆ Z ( n ) N , namely thenumber of ways of flipping n spins with a specific Boltz-mann weight, or counts, and the corresponding Boltz-mann weights . We consider the spin-1 Ising model de-scribed by the Hamiltonian (2) for a square lattice as anexample. For the sake of simplicity, we introduce thefollowing parameters x = e − βJ , y = e − βH , z = e − βL , u = e − βD , w = e − βK , (12)where β = k B T . Then we need to calculate the counts fordifferent configurations. The computation of the countsis more complicated with respect to that of the spin-1/2Ising model because in the spin-1 Ising model each spincan choose among three possible states. Table I showsBoltzmann weights and their counts for some configura-tions in which a few spins are flipped. Starting from thesecalculations summarized in Table I we finally obtain thelow-temperature series expansion of the partition func-tion for the spin-1 Ising model : Z squarespin1 = e − βE (cid:18) N x yz uw + 2 N x y z u w + N x y z − N x y z u w + 4 N x y z uw +6 N x y z u w − N x y z u w +2 N x y z − N x y z uw +6 N x y z u w + 31 N x y z u w + N x y z u w + 12 N x y z u w +18 N x y z u w + O ( x ) (cid:19) (13)We will follow the same procedure by finding the Boltz-mann weights and the associated counts for the 3D cubiclattice in one of the special cases outlined in section 3. C. Metropolis Monte Carlo method
Monte Carlo method is a powerful numerical toolwidely used to evaluate discrete spin models like Isingmodels to investigate the behavior of the associated ther-modynamic functions of the models. It is also very popu-lar to study continuous spin systems like XY and Heisen-berg model, fluids, polymers, disordered materials, andlattice gauge theories (cf. where Monte Carlo Meth-ods are studied). In this paper we use Metropolis MonteCarlo simulation. The algorithm can be summarized inthree steps:1. Set up of the lattice sites. To do it, for examplein the case of 2D square lattice, we define a 2D arraywith l x × l x spin, which is called spin[ i ][ j ], where i and j determine a specific lattice site.2. Initialization of the system. We use a functionnamed init ( l x , J , K , D , L , H ) to set the initial state ofthe system, which in principal can be chosen arbitrarily.In the simulations we choose a completely ordered stateas the initial state with maximum magnetization. In thisfunction, l x is the number of spins in each direction, and J , K , D , L , and H are the values of the coefficients in theHamiltonian of the spin-1 Ising model given by equation(2).3. Use of a main loop in the main program to updatethe system many times. The function mc( T ) takes thetemperature T and uses Metropolis update. This func-tion at first chooses randomly one of the spins in thelattice and flips that spin. For instance, if the randomlychosen spin is +1 it is flipped to 0 or − P (initial → final) = (cid:40) , if E final < E initial e ( − β ( E final − E initial )) , otherwise By updating the system a sufficient number of times, iteventually reaches the equilibrium state at any temper-ature. Finally, it is possible to determine the thermody-namic functions such as magnetization, and susceptibilityusing following formulas : M = 1 N N (cid:88) i =1 s i , (14) χ = 1 k B T ( < M > − < M > ) . (15) III. RESULTS AND DISCUSSIONS: SPECIALCASES OF SPIN-1 ISING HAMILTONIAN FOR2D SQUARE AND 3D CUBIC LATTICES
Owing to the previous arguments, in principle we knowhow to calculate the partition function associated with(2) expressing the Hamiltonian of the spin-1 Ising modeland consequently we can characterize thermodynamicallythe model. However, since the general form of the Hamil-tonian is very complex and the number of parametersappearing in the parameter space is high, it is not possi-ble to study analytically and/or numerically in an exactway the critical behavior of the full Hamiltonian. It isthus useful to understand better the critical behavior ofthe spin-1 Ising model focusing our attention on somespecial cases with a lower number of parameters that arenumerically solvable. These cases will be discussed in thefollowing subsections. In the last subsection we will in-troduce the long-range spin-1 Ising model in analogy withthe well-known spin-1/2 model to find out how its criticaltemperature depends on the magnitude of the long-rangeinteraction.
A. Case 1
In this subsection, we restrict ourselves to the specificcase in which all the coefficients in (2) are zero except J and H , so that the Hamiltonian is given by H (1)Ising1 = − J (cid:88)
7) + N ( N − N + 2)6 − N ( N − − N (cid:19) x y + 21 N x y + 77 N x y + (cid:18) N ( N −
17) + 12 N ( N − N ( N − N ( N −
20) + 6 N ( N − (cid:19) x y + O ( x ) (cid:21) , (22) M (1)cube =1 − x − x + 5 x − x + 108 x − x − x − x + 1602 x + O ( x ) , (23) χ (1)cube = β [ x + 12 x − x + 189 x − x + 129 x + 336 x + 1232 x − x + O ( x )] . (24)Equations (20) and (21) obtained from combinatorialcombinations are in agreement with the results foundin using finite lattice method for a 2D square lattice.As shown by Enting, Guttmann and Jensenin at leastthe first 60 terms of low-temperature series expansion ofthermodynamic functions are needed to have a physicallyconsistent result. The critical temperature has been fi-nally approximated as followsexp( − Jk B T (1 , square) c,SE ) = x (1 , square) c,SE = 0 . ± . . (25) Figure 1. Magnetization M vs. reduced temperature t = k B TJ obtained by using Metropolis Monte Carlo simulation for asquare lattice with 40 ×
40 sites for the case 1. The error barsare smaller than size of the markers.Figure 2. Binder Cumulant vs. t = k B TJ for different sizelattices for Case 1. The error bars are smaller than size of themarkers. Also the calculation of the critical exponents β and γ associated to M and χ , respectively lead to the conclu-sion that the model belongs to the same universality classas spin-1/2 Ising model. Now, we apply the MetropolisMonte Carlo simulation on a 2D square lattice and exam-ine these conclusions. Figure 1 shows the magnetizationversus the reduced temperature, t = k B TJ for a squarelattice with 40 ×
40 spins. It shows that in the hightemperature regime, system is in a disordered phase andthe magnetization is zero, and for a reduced temperature t = k B TJ (cid:39) . U = 1 − < M > < M > , (26)is an observational tool to estimate critical points. Itturns out that the intersection of U − T curves, for net-works with different number of sites, gives the criticaltemperature of the lattice with a good accuracy . Fig-ure 2 shows how we can use the Binder cumulant to deter-mine the critical temperature considering three differentsize lattices. In this figure the Binder cumulants for three2D square lattices with 5 ×
5, 10 ×
10, and 15 ×
15 sitesare displayed. The intersection of the curves correspondsto the critical point given by k B T (1 , square) c, MC J = 1 . ± . . (27)This result is in the accordance with the critical tem-perature found from series expansion method given by(25).In order to complete our discussion about this specificcase we shall find some of the critical exponents of themodel using the data we have obtained from simulationby the so called finite lattice method . Let us firstlyevaluate the β critical exponent. In order to calculatethe β exponent what is usually done is to plot y = ln M vs. x = ln l x . The slope of this graph is the β exponent.Likewise the slope of y = ln χ vs. x = ln l x gives the γ exponent. We eventually find β = 0 . ± . , (28) γ = 1 . ± . . (29)These observations suggest that spin-1 Ising modelgoverned by Hamiltonian (16) belongs to the same uni-versality class of the spin-1/2 Ising model in agreementwith series expansion method. In order to have an ap-proximation of the critical temperature for a 3D cubiclattice we use again the Binder cumulant. The criticaltemperature is: k B T (1 , cube) c, MC J = 3 . ± . . (30) B. Case 2
In this subsection, we consider the Hamiltonian of thespin-1 Ising model given by (2) and we assume J = L =0. Therefore, the Hamiltonian of the system is: H (2)Ising1 = − K (cid:88)
0, and H → D max = − Kz . (40)Since we know the exact critical temperature of the spin-1/2 Ising model for the 2D square lattice, we can easilyfind the critical temperature for case 2 for such a lat-tice. The critical temperature of the 2D square latticespin-1/2 Ising model with Hamiltonian (38) equation orequivalently for the spin-1 Ising model described by (31)is T (2 , square) c = K k B ln(1 + √
2) (41)We now perform Metropolis Monte Carlo simulation for2D square lattice to compare the simulation results withthe analytical calculation with special regard to the crit-ical temperature. We set K = 1, and use the suitablevalue for D , and change H , from +1 to − T . In Figure 4 the results for a 40 ×
40 2D lat-tice are shown. The magnetization M is plotted versusthe external field, H for different values of temperature.As it can be seen, for very low T , there are two jumpsin the curve. One jump corresponds to a positive value Figure 4. M − H curves for case 2 with D = − .
8, and K = 1 at different temperatures: k B T = 0 . k B T = 0 . k B T = 0 . k B T = 1 .
5, for a 40 ×
40 square lattice obtainedfrom Metropolis Monte Carlo simulation. The error bars aresmaller than size of the markers. of H , in which the magnetization jumps from 1 to zero.Another jump occurs at a negative value of H , and inthis case magnetization has a sudden variation from zeroto −
1. This result is conceptually is in accordance withour analytical solution. In fact, we have two first orderphase transitions corresponding to the two values of H ,where coefficient Rβ is zero. Moreover, the critical tem-perature can be estimated by plotting the susceptibilityversus temperature; it turns out that k B T (2 , square) c, MC K = 0 . ± .
01 (42)which agrees with the expression of the analytical criticaltemperature given by (41). The qualitative behavior ofthe 3D lattice system is very similar to the square latticeas we expect and the critical phase transition occurs at k B T (2 , cube) c, MC K = 1 . ± . . (43) C. Case 3
In this section we consider the Hamiltonian (2) with K = 0: H (3)Ising1 = − J (cid:88)
67. On the other hand, at very low-temperatureswe have a completely ordered lattice with all sites spinup or down, i.e. | M | ≈
1. We emphasize that, for L = 0,the model reduces to the Blum-Capel model and corre-sponding mean-field expressions for free energy, magne-tization, and τ can be simply found from (45), (46) and(47). In particular, for H = 0 we get M = 2 e βD sinh [ βzJM ]1 + 2 e βD cosh [ βzJM ] . (48)Equation (48) enables us to calculate mean-field approx-imation for transition temperature, T ; it’s enough toexpand the expression on the right hand of (48) up tothe first order of M . As M → βzJM ] → βzJM and cosh[ βzJM ] →
1, thus zJβ = 1 + 12 e − β D , (49)where β = k B T . Equation (49) determines the curve ofphase transition in phase diagram of Blum-Capel modelat H = 0 plane, however up to now we do not know thetype of phase transition occurring at T . Using (46) we write down the first few terms of series expansion of thezero external field free energy around M = 0: F MF ( H = 0) (cid:39) a + a M + a M , (50)with a = − N k B T ln (1 + 2 e βD ) , (51) a = zJ − βzJ e βD e βD ] , (52) a = zJ
24 ( βzJ ) e βD e βD [ 6 e βD e βD − . (53)The essential condition for having a critical phase tran-sition according to Landau theory is a = 0 , a > . (54)Hence, to have a critical phase transition we get accord-ing to mean-field theory D > − zJ ln 43 . (55)In addition, at the tricritical (tc) point where the threephases predicted by the classical spin-1 Ising model be-come critical simultaneously we must have a = 0 , a = 0 . (56)Then, the tc point is determined by β (3)tc , MF zJ = 3 , D (3)tc , MF = − zJ ln 43 . (57)Mean-field solution in general is not the exact solu-tion but only approximated because it neglects the effectof dimensionality. The results of mean-field calculationsbecome more precise when the dimensionality of the sys-tem becomes larger and the mean-field predictions, likefor example the critical exponents, become exact if thedimensionality of the system is equal or higher than theupper critical dimension, d up , which is given by d up = ( γ + 2 β ) ν . (58)Mean-field usually gives good predictions for the phasediagrams of the 3D systems . In this respect, the inter-esting example is a 3D system like a cubic lattice at thetc point. At the tc point we have the following criticalexponents β tc = 14 , ν tc = 12 , γ tc = 1 , (59)(58) and (59) lead to d up = 3 . (60)It means that at the tc point, mean-field theory pro-vides a very good description of the model in 3D lat-tices in terms of critical exponents. So for the spin-1 3D Figure 5. Absolute magnetization M vs. reduced temperature t = k B T /J for different values of D for a 2D square latticewith 20 ×
20 sites described by the Blum-Capel Hamiltonianin zero external field according to Metropolis Monte Carlosimulation for
D/J = 1 . D/J = − . D/J = − . D/J = − .
97 and
D/J = − .
10. The critical phase transi-tion becomes first-order at D (3 , square)tc , MC /J = − . ± .
01, and t (3 , square)tc , MC = 0 . ± .
01. The error bars are smaller than sizeof the markers.
Ising model, which is described by the Hamiltonian (44),around the tc point, for zero external field and L = 0 thecritical exponents are the ones given by (59).Now we use the Metropolis Monte Carlo technique toinvestigate the behavior of 2D square and 3D cubic lat-tices defined by the Hamiltonian (44) and compare theresults with those of mean-field theory. Figure 5 showsthe absolute magnetization versus reduced temperaturefor different values of D obtained by Metropolis MonteCarlo simulation for a 2D square lattice. It agrees qual-itatively with mean-field but obviously leads to differentvalues for the tc point: D (3 , square)tc , MC J = − . ± . ,k B T (3 , square)tc , MC J = 0 . ± . . (61)We have already seen that mean-field approximationsuggests that the spin-1 Ising model governed by (44)does not exhibit any phase transitions when L isnonzero.Interestingly, Metropolis Monte Carlo simula-tion that is a more accurate method confirms this re-sult. For instance, Figure 6 proves the non-existence ofthe phase transition for a 2D square lattice with 20 × L is non-zero. For a 3Dcubic lattice the system behavior is similar and tc point Figure 6. Magnetization as a function of L plots obtained fromMetropolis Monte Carlo simulation for the spin-1 Ising modeldefined by Hamiltonian (44) for a 20 ×
20 square lattice with H = D = J = 0 for k B T /L = 1 . k B T /L = 2 . Bottom : k B T /L = 10 .
0, which indicates that no phase transition takesplace. Error bars are smaller than size of the markers. is given by D (3 , cube)tc , MC J = − . ± . ,T (3 , cube)tc , MC J = 1 . ± . . (62) D. Case 4
As another specific case we consider the following spin-1 Ising model H (4)Ising1 = − J (cid:88)
20 2D spin squarelattice governed by the Hamiltonian (63).
K/J k B T tc /J D tc /J term proportional to J expresses an exchange interactionbetween two neighbor sites if and only if both are occu-pied and the amount and sign of this interaction dependon spin variables of these two occupied sites. The secondterm of H lg proportional to K is the interaction energybetween a pair of filled neighbors regardless of their spinsand the last term proportional to D is a spin independenteffect of some external field with occupied sites with D playing the role of this field. Now we are ready to provethe equivalence of H (4)Ising1 given by (63) and H lg givenby (64). To do it we start from (64) and impose thefollowing transformation r i = t i s i (65)Equation (65) illustrates that r i can be +1, −
1, or 0.Obviously r i only has two possible values: +1, or 0. Soin the two last terms of (64) we can substitute t i with r i .Thus H lg in terms of new spin variable r i is H lg = − J (cid:88)
1) components .The model was originally inspired by the experimentalobservation that the continuous superfluid transition inHe with He impurity becomes a first order transitioninto normal and superfluid phase separation above somecritical He concentration. Blume, Emery and Griffithshave found the mean-field solution and have determinedthe approximated phase diagram and the tc of the model.Since we have found that Metropolis Monte Carlo re-sults qualitatively agree with the mean-field solution, e.g.phase diagrams obtained from the two methods are sim-ilar, in Table II we present the Metropolis Monte Carlosimulation results for the tc temperature T tc and the tccrystal field of strength D tc for a 2D square lattice. These results show that with increasing K there is an increaseof T tc and a negative increase of D tc . E. Case 5: long-range Ising spin-1 model
Long-range interaction that are typical of statisticalmechanical systems may affect the critical behavior ofthe corresponding models . Regarding this, in this sec-tion we deal with the 2D spin-1 Ising model with long-range spin interactions. We investigate the effect of thisfurther interaction using Metropolis Monte Carlo simu-lation comparing the results with the ones of the corre-sponding 2D spin-1/2 Ising model in the presence of thesame interaction. In analogy with long-range spin-1/2Ising model we define the long-range Hamiltonianfor spin-1 Ising model as follows H (5)lr = − (cid:88) ij Jr d + σij s i s j , (67)where s i = 0, 1, or − d is the lattice dimensionality, σ is the phenomenological parameter which determines theinteraction strength, r ij is the distance between a coupleof spins labeled by the indices i and j , and the subscript lr denotes long range. Explicitly, in the Metropolis algo-rithm if the flipped spin is at position ( x, y ) we have r ij = (cid:113) ( x − i ) + ( y − j ) . (68)In this analysis we limit ourselves to ferromagnetic ma-terials, i.e. J >
0. Notice that for spin-1/2 Ising modelthe Hamiltonian (67) has the same expression but s i = 1,or −
1. On the basis of this numerical simulation we caninvestigate the dependence of the critical temperature on σ . We do the Metropolis Monte Carlo simulation for dif-ferent values of parameter σ , considering the long-rangeinteraction and assuming that each spin interacts withother spins which their distance is equal or less than someradius R . In other words the summation in the Hamil-tonian (67) is carried out over all spins, which are in thecircle of radius R around the flipped spin in each stepof Metropolis algorithm. It means that in the algorithm r ij < R . Figure 7 shows the result for R = 7. As weexpect the critical temperature of the model decreaseswith parameter σ . The dependence of critical tempera-ture on parameter σ , basing upon the results of numericalmodelling of the least squares for spin-1/2 is given by k B T lr c, spin − / J = 2 . . e − . σ . (69)Similarly for spin-1 model we get: k B T lr c, spin − J = 1 . . e − . σ . (70)So in both cases we have T c = a + be − cσ , (71)where a is the critical temperature of the short-rangemodel and c ≈ . Figure 7. t c = k B T c J as a function of σ for the long-rangeinteraction Ising model Top : spin-1/2,
Bottom : spin-1 for asquare lattice with 40 ×
40 sites described by Hamiltonian (67)with R = 7 obtained from Metropolis Monte Carlo simulation.Error bars are smaller than size of the markers. IV. CONCLUSION
In the present work we studied the classical spin-1 Isingmodel using different analytical and numerical methodssuch as mean-field theory, series expansions and MonteCarlo simulation to investigate some critical properties ofthe model like critical temperature and critical exponentsfor 1D chain, 2D square lattice, and 3D cubic lattice. Wehave found that, albeit some similarities with the criticalbehavior of the classical spin-1/2 Ising model, becauseof the presence in the Hamiltonian of the spin-1 modelof more terms, i.e. more different types of interactionsbetween spin pairs, the critical properties of this modelare much richer and more variegated with respect to theones of the corresponding spin-1/2 Ising model.We have used mean-field theory that represents astrong mathematical tool to study the physics of themodel in some special cases. We have found that, for 3Dlattices near the tricritical point, the critical properties ofthe model can be described by mean-field theory. In par-ticular, we have found that the critical exponents aroundthe tricritical point calculated via the mean-field approx-imation are confirmed by Monte Carlo simulations. Onthe other hand, the mean-field results for 2D lattices areonly qualitatively but not quantitatively correct as high-lighted by Monte Carlo simulation. The simulation re-sults obtained for 2D square and 3D cubic lattices canbe easily extended to other types of lattices.We have shown that, for a special case of the spin-1Ising model where the bilinear and the bicubic terms areset equal to zero, it is possible to write the correspondingpartition function in arbitrary dimensions as the one ofthe spin-1/2 Ising model in agreement with our MonteCarlo simulation. Finally, we have investigated the long-range spin-1 Ising model Hamiltonian by including inthe Hamiltonian a long-range interaction term in anal-ogy with what was carried out for spin-1/2 Ising modeldetermining the dependence of the critical temperatureof the two models on the strength of this interaction.
ACKNOWLEDGEMENTS
This work was partially supported by National Groupof Mathematical Physics (GNFM-INdAM) and IstitutoNazionale di Alta Matematica F. Severi. ∗ [email protected] E. Ising, Zeitschrift f¨ur Physik , 253 (1925). S. G. Brush, Rev. Mod. Phys. , 883 (1967). R. Peierls and M. Born, Proceedings of the CambridgePhilosophical Society , 477 (1936). L. Onsager, Phys. Rev. , 117 (1944). R. J. Baxter, Journal of Statistical Physics , 1164(2012). C. N. Yang, Phys. Rev. , 808 (1952). J. M. Yeomans,
Statistical Mechanics of Phase Transitions (Oxford University Press, Oxford, 1993). M. E. Newman,
Monte Carlo Methods in Statistical Physics (Oxford University Press, New York, 1999). A. F. Sonsin, M. R. Cortes, D. R. Nunes, J. V. Gomes,and R. S. Costa, Journal of Physics: Conference Series , 012057 (2015). T. Kaneyoshi and J. Mielnicki, Journal of Physics: Con-densed Matter , 8773 (1990). A. Benyoussef, T. Biaz, M. Saber, and M. Touzani, Jour-nal of Physics C: Solid State Physics , 5349 (1987). T. Kaneyoshi, Journal of the Physical Society of Japan ,933 (1987). T. Kaneyoshi, Journal of Physics C: Solid State Physics , L679 (1988). N. Boccara, A. Elkenz, and M. Saber, Journal of Physics:Condensed Matter , 5721 (1989). T. Kaneyoshi, Journal of Physics C: Solid State Physics , L557 (1986). J. W. Tucker, Journal of Physics: Condensed Matter ,485 (1989). E. Juriinov and M. Juriin, Physica A: Statistical Mechanicsand its Applications , 554 (2016). M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev.A , 1071 (1971). K. Huang and P. Beale,
Statistical Mechanics (Wiley, NewYork, 1988). I. G. Enting, A. J. Guttmann, and I. Jensen, Journal ofPhysics A: Mathematical and General , 6987 (1994). P. F. Fox and A. J. Guttmann, Journal of Physics C: SolidState Physics , 913 (1973). D. M. Saul, M. Wortis, and D. Stauffer, Phys. Rev. B ,4964 (1974). J. Yeomans and M. E. Fisher, Phys. Rev. B (1981),10.1103/PhysRevB.24.2825. A. K. Jain, Phys. Rev. B (1980), 10.1103/Phys-RevB.22.445. C. M. Care, Journal of Physics A: Mathematical and Gen-eral , 1481 (1993). K. Binder and D. Heermann,
Monte Carlo Simulation inStatistical Physics (Springer, Berlin, 1997). W. Selke, Journal of Statistical Mechanics: Theory andExperiment , P04008 (2007). R. B. Griffiths, Physica , 689 (1967). F. Wu, Chinese Journal of Physics , 153 (1978). J. Ortega and W. Rheinboldt,
Iterative Solution of Non-linear Equations in Several Variables (SIAM PublicationsClassics in Appl. Math, Philadelphia, 2000). J. Als-Nielsen and R. J. Birgeneau, American Journal ofPhysics - AMER J PHYS , 554 (1977). A. Campa, T. Dauxois, and S. Ruffo, Physics Reports , 57 (2009). I. Balog, G. Tarjus, and M. Tissier, Journal of StatisticalMechanics: Theory and Experiment , P10017 (2014). M. E. Fisher, S.-k. Ma, and B. G. Nickel, Phys. Rev. Lett.29