Critical properties of the HintermannMerlini model
Abstract
Many critical properties of the HintermannMerlini model are known exactly through the mapping to the eightvertex model. Wu [J. Phys. C 8, 2262 (1975)] calculated the spontaneous magnetizations of the model on two sublattices by relating them to the conjectured spontaneous magnetization and polarization of the eightvertex model, respectively. The latter conjecture remains unproved. In this paper, we numerically study the critical properties of the model by means of a finitesize scaling analysis based on transfer matrix calculations and Monte Carlo simulations. All analytic predictions for the model are confirmed by our numerical results. The central charge is found for critical manifold investigated. In addition, some unpredicted geometry properties of the model are studied. Fractal dimensions of the largest Ising clusters on two sublattices are determined. The fractal dimension of the largest Ising cluster on the sublattice A takes fixed value , while that for sublattice B varies continuously with the parameters of the model.
pacs:
05.50.+q, 64.60.Cn, 64.60.Fr, 75.10.HkI Introduction
The exact solutionsIsing (); YangIsing () of the twodimensional (2D) Ising model significantly promote the research of phase transitions and critical phenomena. After that, the Ising model becomes one of the most famous lattice models in statistical physics. Another famous Ising system, known as the BaxterWu modelbw (), is defined on the triangular lattice with pure threespin interactions. The model was firstly proposed by Wood and Griffithswood () and exactly solved by Baxter and Wubw () by relating the model with the coloring problem on the honeycomb lattice. The solution gives the critical exponents () and (), which are exactly the same as those of the 4state Potts modelpotts (); wfypotts (), which can be derived via the Coulomb gas theorydenNijs1 (); denNijs2 (); cg (). This means that the BaxterWu model belongs to the universality class of the 4state Potts model. The critical properties of the 4state Potts model are modified by logarithmic correctionslog2 (); log () due to the second temperature field, which is marginally irrelevantNBRS (); log2 (). With the two leading temperature fields simultaneously vanishing, the leading critical singularities of the BaxterWu modeldengbw () do not have logarithmic factors. Deng et al. generalize the BaxterWu model in dengbw () where the spins are allowed to be states ( can be larger than 2) and the up and downtriangles can have different coupling constants. Both generalizations lead to discontinuous phase transitions.
In 1972, Hintermann and Merlinihintermann () considered an Ising system on the UnionJack lattice (as shown in Fig. 1)
(1) 
where the sum takes over all the square unit cells and () are the coupling constants. The model is similar to the BaxterWu model in the sense that both are Ising systems with pure threespin interactions. However the critical properties of the HM model are much more complicated. For the ferromagnetic case (), the model has fourfold degenerate ground states, as shown in Fig. 2. The nature of the phase transition breaking the symmetry of the order parameter cannot be determined by the dimensionality of the system and the symmetry properties of the ground statesSuzuki (). The exponents may vary with some tuning parameter without changing the symmetry of the order parameter. Such behavior has been found, e.g., in the AshkinTeller model (AT) AT (); ATMC (), the eightvertex model 8vertex (); exactbook (), the 2D model in a fourfold anisotropic fieldXYh (), and the ferromagnetic Ising model with antiferromagnetic nextnearest neighbor couplings on square latticej1j2henk (); j1j2sandvik (); j1j2guo ().
Through mapping to the eightvertex model8vertex (); exactbook (), the free energy of the HM model has been found exactlyhintermann (). The critical manifold and critical exponent varying with the ratio of couplings were obtained. Based on the conjectured spontaneous magnetization 8vmag () and the spontaneous polarization 8vpolar () of the eightvertex model, Wu calculated the spontaneous magnetizations of the model on two sublattices. The results show that the two magnetizations possess different critical exponentswubw (). The spontaneous magnetization of the eightvertex model was derived later 8vmagproof (), however, the spontaneous polarization remains a conjecture. It is thus very necessary to verify the results numerically.
In current paper, we study the critical behavior of the ferromagnetic HM model with and numerically. The numerical procedure includes transfer matrix calculations and Monte Carlo simulations. The critical properties we studied includes not only the verification of the analytic predictions, but also the fractal structure of spin clusters, which is studied for the first time.
The paper is organized in the following way: In Sec. II, we summarize the theoretical results of the model. In Sec. III, we introduce the transfer matrix method and present the associated numerical results. In Sec. IV, we describe the algorithm used in Monte Carlo simulations and give the associated numerical results. In Sec. V, we define Ising clusters on two sublattices and numerically determine the fractal dimensions of the corresponding largest clusters, respectively. We summarize in Sec. VI.
Ii Exact solutions
In this section, we summarize the analytical results for the HM modelhintermann (); wubw ().
Assign arrows between the nearestneighboring Ising spins of sublattice A (Fig. 3): if the two spins beside the arrow are the same, the arrow is rightward (or upward), otherwise it is leftward (or downward). There is a twotoone correspondence between the Ising spin configurations and the arrow configurations.
After taking sum over the spin (in the center of the unit cell), the Boltzmann weights of the vertices are
(2) 
Thus the model is mapped to the symmetric eightvertex model8vertex (); exactbook (), which has been exactly solved by Baxter. For the ferromagnetic case considered in current paper, the critical manifold is given by
(3) 
which leads to for the uniform case ().
The singularity of the free energy density is governed by
(4) 
when is not an integer (The case that = an integer is not considered in current paper). is the critical temperature and is given by wubw ()
(6)  
For the ferromagnetic case that we consider, is determined by (6). Since and , the critical exponent is thus found
(7) 
For the uniform case, (7) gives .
A very interesting critical property of the HM model is that the spontaneous magnetizations and of the sublattice A and B possess different critical exponents:
(8) 
The critical exponents and were obtained by Wuwubw ()
(9) 
(It should be noted that the author reversely wrote the two critical exponents in Eq. (16) of wubw ().) According to the scaling law, this leads to two magnetic exponents and
(10) 
Here has fixed value , while varies continuously with the parameters of the model.
In following sections, we will numerically study the critical properties of the model. The critical points and critical exponents will be numerically verified. Especially, Wu’s result of () is based on the unproved conjecture of the spontaneous polarization8vpolar () of the eightvertex model, thus numerical verification is very necessary. Our numerical studies are focused on the subspace that and . By defining , critical manifold and exponents are expressed as functions of . It is worthy to notice that there is a symmetry for the transformation . At a special point and its dual , , the analytic results give . This is the 4state Potts point.
Iii Transfer matrix calculations
As shown in Fig. 4, we define the rowtorow transfer matrix
(11)  
where and are the states of two neighboring rows, respectively. Here periodic boundary conditions are applied. For a system with rows, the partition sum is found to be
(12) 
with the periodic boundary conditions applied. In the limit , the free energy density is determined by the leading eigenvalue of
(13) 
For the HM model, the dimension of the matrix is . To numerically calculate the eigenvalues of , we used the sparse matrix technique, which sharply reduces the requirement of computer memory for storing the matrix elements. (For details of this technique, see henk1982 (); henk1989 (); henk1993 (); Qian2004 ().) We are able to calculate the eigenvalues of with up to 22. In our calculations, we restrict the system size to even values, because the ordered configurations, as shown in Fig. 2, do not fit well in odd systems.
The critical properties can be revealed by calculating three scaled gaps :
(14) 
where are three subleading eigenvalues of the matrix , respectively. According to the conformal invariance theoryconformal (), the scaled gap is related with a correlation length
(15) 
where governs the decay of a correlation function . According to the finitesize scalingfss (), the gap in the vicinity of a critical point scales as
(16)  
where is the corresponding scaling dimension, which is related to an exponent according to the conformal invarianceconformal (). is an irrelevant field, and is the corresponding irrelevant exponent. ,, and are unknown constants.
We focus on the the energyenergy correlation and two magnetic correlations and . Generally speaking, the magnetic correlation function is defined as . For the HM model, since the spontaneous magnetizations on the two sublattices behave differently, we look two types of magnetic correlation functions: with and on the A sublattice and on the B sublattice.
Let be the largest eigenvalue in the subspace that breaks the spin updown symmetry, which means that the associated eigenvector satisfies
(17) 
where is the operator flipping spins. Thus, the scaled gap is identified as and the corresponding correlation is .
Let and be the second and the third largest eigenvalue in the subspace keeping the spin updown symmetry. It is not a priori clear which of the two corresponding gaps is the thermal one or the magnetic one . Based on the magnitudes of and , we identify and as and , respectively.
We then numerically solve the finitesize scaling equation
(18) 
for , respectively. The solution satisfies
(19) 
Here is an unknown constant. Since and , converges to the critical point with increasing system sizes. (19) is used to determine the critical point in our numerical procedure. For the uniform case, we obtain , which is in good agreement with the exact solution . We have also estimated the critical points of the cases with , and 5, where and . Both and have been used to estimate the critical points, we list the best estimations in Table 1. Our numerical estimations of are consistent with the theoretical results in a high accuracy.
Exactly at the critical point , (16) reduces to
(20) 
which is used to determine the scaling dimensions and . For the uniform case , it gives , and , which are consistent with the analytic results (Table 1). We applied the same procedure to . The numerical estimations and the associated analytic results are consistent, see Table 1.
We also calculate the scaled gap in the vicinity of critical point. The scaling behavior of seems different from that of or , as shown in Fig. 7 for the uniform case. This is the result of in Eq. (16). The location of the minimum of for a given can be denoted as , which converges to the critical point when . This property can also be used to estimate , e.g., as was done in henkcubic (). But we don’t do this in current paper. By calculating at the estimated , we obtain the scaling dimension according to (20). The results are listed in Table 1.
The finitesizescaling behavior of the free energy density at the critical point determines the central charge according to fren1 (); fren2 ()
(21) 
Fitting the data of the free energy density according to (21), we obtain the central charge for the uniform case. For the other cases, the numerical estimations of the central charge are listed in Table 1. The central charge for all takes the fixed value as expected for transitions breaking fourfold symmetry.
Furthermore, we consider the 4state Potts point . A finitesize scaling analysis based on transfer matrix calculations is performed at this point. The numerical estimations obtained are listed in Table 1. No logarithmic corrections are found, as in the BaxterWu model.
1  2  3  4  5  
T  0.4406867935  0.3046889317  0.2406059125  0.225147108  0.2017629641  0.1751991102  
TM  0.44068679(5)  0.3046889(1)  0.2406059(1)  0.2251472(2)  0.2017629(1)  0.1751991(2)  
T  4/3  1.39668184  1.47604048  3/2  1.53960311  1.58921160  
MC  1.332(2)  1.397(2)  1.478(3)  1.500(1)  1.540(3)  1.590(2)  
T  2/3  0.60331816  0.52395952  1/2  0.46039689  0.41078840  
TM  0.66666(1)  0.603318(1)  0.52396(1)  0.50000(1)  0.46040(2)  0.41079(1)  
T  15/8  15/8  15/8  15/8  15/8  15/8  
MC  1.874(1)  1.874(2)  1.876(2)  1.874(2)  1.877(3)  1.874(2)  
T  1/8  1/8  1/8  1/8  1/8  1/8  
TM  0.125000(1)  0.125000(1)  0.125000(1)  0.12499(1)  0.12500(1)  0.12498(3)  
T  11/6  1.84917046  1.86901012  15/8  1.88490078  1.89730290  
MC  1.833(1)  1.848(2)  1.870(2)  1.875(1)  1.886(2)  1.896(4)  
T  1/6  0.15082954  0.13098988  1/8  0.11509922  0.10269710  
TM  0.166666(1)  0.150829(1)  0.130990(1)  0.12500(1)  0.115099(1)  0.102696(2)  
T  1  1  1  1  1  1  
MC  1.000000(1)  1.000000(1)  0.99998(3)  1.0000(1)  1.0000(2)  0.9999(1) 
Iv Monte Carlo simulations
The BaxterWu model has been simulated by the Metropolis algorithmbwMC (), the Wang Landau algorithmwanglandau (), and a cluster algorithmhenkbw (); dengbw (). The cluster algorithm is similar to the SwendsenWang algorithmsw () for the Potts model. In this algorithm, the triangular lattice is divided as three triangular sublattices. By randomly freezing the spins on one of the sublattices, the other two sublattices compose a honeycomb lattice with pair interactions, then a SwendsenWang type algorithm can be formed to update the spins on this honeycomb lattice. In current paper, we suitably modify this cluster algorithm to simulate the HM model.
The algorithm is divided into three steps:

Step 1, update the spins on sublattice A.
Let the spins on sublattice B be unchanged (frozen), the interactions of the spins reduce to pair interactions on sublattice A.
A. bonds. An vertical (horizontal) edge on sublattice A belongs to two triangles, the interactions of the left (upper) and right (down) triangles can be denoted and respectively, which may take the values , or . Place a bond on this edge with probability if both the products of spins in the left (upper) triangle and in the right (down) triangle are 1, if only the product in the left (upper) triangle is 1, if only the product in the right (down) triangle is 1, otherwise.
B. clusters. A cluster is defined as a group of sites on sublattice A connected through the bonds.
C. update spins. Independently flip the spins of each cluster with probability 1/2.

Step 2, update the spins on sublattice B and half of the spins on sublattice A.
This step is very similar to step 1, but we freeze only half of the spins on sublattice A, which are labeled (Fig. 8). The other spins (on sublattice B+) are updated, whose interactions also reduce to pair interactions when spins are frozen.

Step 3, update the spins on sublattice B and the other half of the spins on sublattice A.
In this step, the spins labeled (in Fig. 8) are frozen, other spins (on sublattice B+) are updated.
In a complete sweep, all the spins are updated twice.
In the simulations, the sampled variables include the specific heat , the magnetizations on sublattice A, and on sublattice B.
The specific heat is calculated from the fluctuation of energy density
(22) 
where is the linear size of the system, and means the ensemble average. The magnetizations are defined as
(23)  
(24) 
where is the number of total sites on sublattice A or B.
In order to demonstrate our numerical procedure based on the Monte Carlo simulations, we take the uniform case () as an example. The cluster algorithm is very efficient, which easily allows us to do meaningful simulations for system with linear size up to . All the simulations are performed at the critical point . After equilibrating the system, samples were taken for each system size. Using the LevenbergMarquardt algorithm, we fit the data of according to the finitesize scaling formulafss ()
(25) 
where are the leading correctiontoscaling term with the leading irrelevant exponent. is the dimensionality of the lattice. , , and are unknown parameters. The fitting yields , which is in good agreement with the exact result .
The finitesize scaling behaviors of the magnetizations and are
(26)  
(27) 
respectively. The loglog plot of and versus are shown in Fig. 9. The fits yield , and , where and are in good agreement with the analytic predictions and , respectively.
We also simulate the , and cases. The numerical estimations of and the corresponding analytic predictions are listed in Table 1. They are in good agreement. In these fits, no logarithmic corrections to scaling are found.
V Fractal structure of the model
In addition to calculate the critical exponents , , and , we also investigate the geometry properties of the model.
The configurations of the HM model are represented by Ising spins. Thus we can define “Ising clusters” for the model as in the Ising modelIsingCluster (). To do so, each sublattice is viewed as a square lattice. For two nearestneighboring Ising spins on a sublattice, they are considered to be in the same cluster if they have the same sign. Here, in saying “nearestneighboring”, the neighborhood of the site is the same to that of a site in a square lattice. This definition is applied to the two sublattices respectively, thus there are two types of Ising clusters. For clarity, we call the Ising cluster on sublattice A as “A cluster”, while the Ising cluster on sublattice B as “B cluster”.
The simulations are also performed at the critical points. The size of the largest A cluster and the largest B cluster are denoted and , respectively. They satisfy the finitesize scaling
(28)  
(29) 
which means that the largest A cluster and the largest B cluster are fractals, with and the fractal dimensions, respectively. Figure 10 is an illustration of and versus system size at the uniform case. Fitting the data according to (28) and (29), we find and .
The same procedure were applied to and . The results obtained are listed in Table. 2. The fractal dimension of the A cluster has a fixed value , while that of the B cluster varies with the ratio .
1  2  3  4  5  

1.888(1)  1.889(1)  1.890(3)  1.889(2)  1.888(1)  1.888(2)  
1.925(1)  1.931(1)  1.939(1)  1.941(1)  1.945(2)  1.951(2) 
Vi Conclusion and discussions
In conclusion, we have numerically studied the critical properties of the HM model by means of finitesize scaling analysis based on the transfer matrix calculations and the Monte Carlo simulations. For the critical points and the critical exponents (), (), and (), our numerical estimations are in good agreement with the corresponding analytic predictionshintermann (); wubw (). The analytic prediction for is based on the conjectured spontaneous polarization of the eightvertex model, which remains unproved. Our numerical results verified the correctness of the prediction.
In addition the central charge of the model is found to be . This is consistent with the fact that , and vary continuously with the parameters, which is related to the fourfold degeneracy of the ground states of the model. Usually, logarithmic correctionslog (); lvAT (); j1j2guo () show up in the 4state Potts point for the model with 4fold degenerate ground states due to the second temperature field. However, we do not see such logarithmic corrections in the 4state Potts point of the HM model. This is the same as the BaxterâWu model, in which the amplitude of the marginally irrelevant operator is zero.
Furthermore, two unpredicted critical exponents are found to describing geometry properties of the model. We define clusters based on the Isingspin configurations on the two sublattices of the model. The fractal dimension of the largest cluster on sublattice A takes fixed value , while the fractal dimension of the largest cluster on sublattice B varies continuously with the parameter of the model.
Acknowledgment
Ding thanks F. Y. Wu for valuable discussion in terms of understanding the analytical results of HM model. This work is supported by the National Science Foundation of China (NSFC) under Grant No. 11205005 (Ding) and 11175018 (Guo).
References
 (1) L. Onsager, Phys. Rev. 65, 117 (1944).
 (2) C. N. Yang, Phys. Rev. 85, 808 (1952).
 (3) R. J. Baxter and F. Y. Wu, Phys. Rev. Lett. 31, 1294 (1973).
 (4) D. W. Wood and H. P. Griffiths, J. Phys. C 5, L253 (1972).
 (5) R. B. Potts, Proc. Camb. Phys. Soc. 48, 106 (1952).
 (6) F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
 (7) M. P. M. den Nijs, J. Phys. A 12, 1857 (1979).
 (8) M. P. M. den Nijs, Phys. Rev. B 27, 1674 (1983).
 (9) B. Nienhuis, J. Stat. Phys. 34, 731 (1984).
 (10) B. Nienhuis, A. N. Berker, E. K. Riedel, M. Schick, Phys. Rev. Lett. 43, 737 (1979).
 (11) J. Salas and A. D. Sokal, J. Stat. Phys. 88, 567 (1997).
 (12) M. Nauenberg, D. J. Scalapino, Phys. Rev. Lett. 44, 837 (1980).
 (13) Y. Deng, W.A. Guo, J. R. Heringa, H. W. J. Blöte, and B. Neinhuis, Nucl. Phys. B 827[FS], 406 (2010).
 (14) A. Hintermann and D. Merlini, Phys. Lett. A, 41, 208 (1972).
 (15) M. Suzuki, Prog. Theor. Phys. 51, 1992 (1974).
 (16) J. Ashkin and E. Teller, Phys. Rev. 64, 178 (1943).
 (17) S. Wiseman and E. Domany, Phys. Rev. E 48, 4080 (1993).
 (18) J.P. Lv, Y. Deng, and Q.H. Chen, Phys. Rev. E 84, 021125(2011).
 (19) M. P. Nightingale and H. W. J. Blöte, Physica A 251, 211 (1998).
 (20) S. Jin, A. Sen, and A. W. Sandvik, Phys. Rev. Lett. 108, 045702 (2012).
 (21) S. Jin, A. Sen, W.A. Guo, and A. W. Sandvik, Phys. Rev. B 87, 144406(2013).
 (22) R. J. Baxter, Phys. Rev. Lett. 26, 832 (1971).
 (23) R. J. Baxter, Exactly Solved Models in Statistical Mechanics (Academic, London, 1982).
 (24) J. V. José, L. P. Kadanoff, S. Kirkpatrick, and D. R. Nelson, Phys. Rev. B 16, 1217 (1977).
 (25) M. N. Barber and R. J. Baxter, J. Phys. C 6, 2913 (1973).
 (26) R. J. Baxter and S. B. Kelland, J. Phys. C 7, L403 (1974).
 (27) F. Y. Wu, J. Phys. C 8, 2262 (1975).
 (28) R. J. Baxter, J. Stat. Phys. 15, 485 (1976).
 (29) H. W. J. Blöte and M. P. Nightingale, Physica A (Amsterdam) 112, 405 (1982).
 (30) H. W. J. Blöte and B. Nienhuis, J. Phys. A 22, 1415 (1989).
 (31) H. W. J. Blöte and M. P. Nightingale, Phys. Rev. B 47, 15046 (1993).
 (32) X. F. Qian, M. Wegewijs, and H. W. J. Blöte, Phys. Rev. E 69, 036127 (2004).
 (33) J. L. Cardy, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz (Academic Press, London, 1987), Vol. 11, p. 55, and references therein.
 (34) For reviews, see e.g. M. P. Nightingale in FiniteSize Scaling and Numerical Simulation of Statistical Systems, ed. V. Privman (World Scientific, Singapore 1990), and M. N. Barber in Phase Transitions and Critical Phenomena, Vol. 8, eds. C. Domb and J. L. Lebowitz (Academic, New York, 1983).
 (35) H. W. J. Blöte, M. P. Nightingale, Physica A (Amsterdam) 129, 1 (1984).
 (36) H. W. J. Blöte, J. L. Cardy, and M. P. Nightingale, Phys. Rev. Lett. 56, 742 (1986).
 (37) I. Affleck, Phys. Rev. Lett. 56, 746 (1986).
 (38) L. N. Shchur and W. Janke, arXiv:1007.1838.
 (39) N. Schreiber and J. Adler, J .Phys. A 38, 7253 (2005).
 (40) H. W. J. Blöte, J. R. Heringa, and E. Luijten, Comp. Phys. Comm. 147, 58 (2002).
 (41) R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. 58, 86 (1987).
 (42) W. Janke and A. M. J. Schakel, Phys. Rev. E 71, 036703 (2005).