Critical behavior and entanglement of the random transverse-field Ising model
between one and two dimensions
We consider disordered ladders of the transverse-field Ising model and study their critical properties and entanglement entropy for varying width, , by numerical application of the strong disorder renormalization group method. We demonstrate that the critical properties of the ladders for any finite are controlled by the infinite disorder fixed point of the random chain and the correction to scaling exponents contain information about the two-dimensional model. We calculate sample dependent pseudo-critical points and study the shift of the mean values as well as scaling of the width of the distributions and show that both are characterized by the same exponent, . We also study scaling of the critical magnetization, investigate critical dynamical scaling as well as the behavior of the critical entanglement entropy. Analyzing the -dependence of the results we have obtained accurate estimates for the critical exponents of the two-dimensional model: , and .
In nature there are materials, which are in a way between two integer dimensions, such as they are built from -dimensional layers having a finite width, . Examples are thin filmsthinfilms (), magnetic multilayersmajkrzak91 () or ladders of quantum spinsdagotto96 (). One interesting question for such multilayer systems is the properties of critical fluctuations, when the linear extent of the layers, , goes to infinity. If the system is classical having thermal fluctuations, finite-size scaling theoryfisherme (); barber () can be applied. One basic observation of this theory is that for any finite the critical behavior is controlled by the fixed point of the -dimensional system, but the scaling functions in terms of the variable, , involve also the critical exponents of the -dimensional system. For example the critical points, , measured at a finite width, , approach the true -dimensional critical point, , as
where the shift exponent, , generally corresponds to the correlation-length exponent, , in the -dimensional system.
In a quantum system having a quantum critical point at zero temperature, , by varying a control parameter, , the dimensional cross-over is a more subtle problem. If the -dimensional critical quantum system is isomorphic with a -dimensional classical systemkogut (), then results of finite-size scaling can be transferred to the quantum system, too. This is the case, e.g. for the quantum critical point of the -dimensional transverse-field Ising model, which is equivalent to the critical point of the classical -dimensional Ising model. However, the situation is more complicated for antiferromagnetic models with continuous symmetry, such as for Heisenberg antiferromagnetic spin ladders. In this case the form of low-energy excitations could sensitively depend on the value of : if the ladder contains even number of legs there is a gap, whereas for odd number of legs the system is gaplessdagotto96 (). In the following for quantum systems we restrict ourselves to models with a discrete symmetry, such as to the transverse-field Ising model.
In disordered systems, in which besides deterministic (thermal or quantum) fluctuations there are also disorder fluctuations in a sample of finite width one can define and measure a sample-dependent pseudo-critical point, (or ), and study its distributiondomany (). In particular one concerns the shift of the mean value, , and the scaling of the width of the distribution, . In this case besides the shift exponent, , which is defined analogously to Eq.(1) one should determine the width exponent, , too, which is defined by the scaling relation:
According to renormalization group theoryaharony () the finite-size scaling behavior of random classical systems depends on the relevance or irrelevance of the disorderharris (). If the disorder represents an irrelevant perturbation at the pure system’s fixed point, which happens if the correlation length exponent in the pure system satisfies , than for the disordered system we have and and the thermodynamic quantities at the fixed point are self-averaging. On the contrary for relevant disorder, which happens for , there is a new conventional random fixed point with a correlation-length exponent, ccfs (), and we have . In this fixed point there is a lack of self-averaging. These predictions, which have been debated for some timepsz (), were checked later for various modelsdomany (); aharony (); paz2 (); mai04 (); PS2005 ().
For quantum systems quenched disorder is perfectly correlated in the (imaginary) time direction, therefore generally it has a more profound effect at a quantum critical pointqsg (). In some cases the critical properties of the random model are controlled by a so called infinite disorder fixed pointfisher (), in which the disorder fluctuations play a completely dominant rǒle over quantum fluctuations. This happens, among others for the random transverse-field Ising model, as shown by analytical resultsfisher () in and numerical resultspich98 (); motrunich00 (); lin00 (); karevski01 () in . Finite-size scaling has been tested for the model and a new scenario is observedilrm (). The finite-size transition points, denoted by in a system of length, , are shown to be characterized by two different exponents, . This means, that asymptotically , which is just the opposite limit as known for irrelevant disorder.
In the present paper we go to the two-dimensional problem and study the finite-size scaling properties of ladders of random transverse-field Ising models. For this investigations we use a numerical implementation of the so called strong disorder renormalization group methodim (). As in this method is expected to be asymptotically exact in large scales. In the numerical implementation of the method we have used efficient computer algorithms and in this way we could treat ladders with a large number of sites: we went up to lengths for legs and used random samples. Our aim with these investigations is threefold. First, we want to clarify the form of finite-size scaling valid for this random quantum model. Second, using the appropriate form of the scaling Ansatz we want to calculate estimates for the critical exponents of the 2d model. Previous studiespich98 (); motrunich00 (); lin00 (); karevski01 (); vojta09 () in this respect have quite large error bars and we want to increase the accuracy of the estimates considerably. Our third aim is to calculate also the entanglement entropyAmicoetal08 () in the ladder geometry and study its cross-over behavior between onerefael () and two dimensionslin07 (); yu07 ().
The structure of the rest of the paper is the following. The model and the method of the calculation is presented in Sec. II. In Sec. III finite-size transition points are calculated and their distribution (shift and width) is analyzed. In Sec. IV we present calculations at the critical point about the magnetization and the dynamical scaling behavior. Results about the entanglement entropy are presented in Sec.V. Our paper is closed by a discussion.
Ii Model and method
ii.1 Random transverse-field Ising ladder
We consider the random transverse-field Ising model in a ladder geometry in which the sites, and , are taken from a strip of the square lattice of length, , and width, . We use periodic boundary conditions in both directions. The model is defined by the Hamiltonian:
in terms of the Pauli-matrices, . Here the first sum runs over nearest neighbor sites and the couplings and the transverse fields are independent random numbers, which are taken from the distributions, and , respectively. For concreteness we use box-like distributions: , for and , for ; , for and , for . We consider the system at and use as the quantum control parameter.
In the thermodynamic limit, , the system in Eq.(3) displays a paramagnetic phase, for , and a ferromagnetic phase, for . In between there is a random quantum critical point at and we are going to study its properties for various widths, .
ii.2 Strong disorder renormalization group method
The model is studied by the strong disorder renormalization group methodim (), which has been introduced by Ma, Dasgupta and Humdh () and later developed by D. Fisherfisher () and others. In this method the largest local term in the Hamiltonian (either a coupling or a transverse field) is successively eliminated and at the same time new terms are generated between remaining sites. If the largest term is a coupling, say connecting sites and , ( being the energy-scale at the given RG step), then after renormalization the two sites form a spin cluster with an effective moment , where in the starting situation each spin has unit moment, . The spin cluster is put in an effective transverse field of strength: , which is obtained in second order perturbation calculation. On the other hand, if the largest local term is a transverse-field, say , then site is eliminated and new couplings are generated between each pairs of spins, which are nearest neighbors to . If say and are nearest neighbor spins to , than the new coupling connecting them is given by: , also in second order perturbation calculation. If the sites and are already connected by a coupling, , than for the renormalized coupling we take . This last step is justified if the renormalized couplings have a very broad distribution, which is indeed the case at infinite disorder fixed points. The renormalization is repeated: at each step one more site is eliminated and the energy scale is continuously lowered. For a finite system the renormalization is stopped at the last site, where we keep the energy-scale, , and the total moment, , as well as the structure of the clusters.
ii.3 Known exact results in the chain geometry
The renormalization has special characters in the chain geometry, i.e. with . In this case the topology of the system stays invariant under renormalization and the couplings and the transverse fields are dual variables. From this follows that at the quantum critical point the couplings and the transverse fields are decimated symmetrically, thus the critical point is located at pfeuty (). The RG equations for the distribution function of the couplings and that of the transverse fields can be written in closed form as an integro-differential equation, which has been solved analytically both at the critical pointfisher () and in the off-critical region, in the so-called Griffiths-phasei02 (). Here we list the main results.
The energy-scale, , and the length-scale, , are related as:
with an exponent: . (Here can be the size of a finite system and is a reference energy scale.) The average spin-spin correlation function is defined as , where denotes the ground-state average and stands for the averaging over quenched disorder. In the vicinity of the critical point has an exponential decay:
in which the correlation length, , is divergent at the critical point as:
with . At the critical point the average correlations have a power-law decay:
with a decay exponent:
The average cluster moment, , is related to the energy-scale, as:
with . The average cluster moment can be expressed also with the size as , where the fractal dimension of the cluster is expressed by the other exponents as:
Iii Finite-size critical points
iii.1 Results in the chain geometry
In the chain geometry finite-size critical points are studied in Ref.ilrm (), in which they are located by different methods, which all are based on the free-fermion mapping of the problemlsm (). The finite-size critical points are shown to satisfy the micro-canonical condition:
from which follows that the distribution of is Gaussian with zero mean and with a mean deviation of . Consequently the width-exponent of the distribution is given by:
On the other hand the shift-exponent is given by:
although in some cases (c.f. for periodic boundary conditions) the prefactor of the scaling function can be vanishing.
iii.2 The doubling method
In the ladder geometry, i.e. for , the free-fermionic mapping is no longer valid, therefore new methods have to be utilized to locate pseudo-critical points. Here we used the doubling method combined with the strong disorder renormalization group.
In the doubling procedurePS2005 () for a given random sample () of length and width , we construct a replicated sample () of length and width by gluing two copies of () together and study the ratio of the magnetizations: , which are calculated by the strong disorder renormalization group method. In Fig.1 we illustrate the dependence of the total magnetic moments, and , for a given sample of a -leg ladder. The corresponding ratio of the magnetizations, , is shown in the upper inset of this figure. It is seen, that in the ordered phase: this ratio approaches . On the other hand in the disordered phase: the magnetizations approach their minimal values, which in the SDRG method can be and , respectively, thus we have . In between there is a sudden change in the value of this ratio, which can be used to define a sample-dependent pseudo-critical point, .
There is another possibility, if we consider the ratio of the two gaps: , which are also calculated by the strong disorder renormalization group method. In Fig.1 we show the two log-gaps, and , for the same sample as before and the corresponding ratio, , is put in the upper inset of this figure. It is seen that this ratio in the ordered phase, , approaches and in the disordered phase: , goes to . In between this ratio has a quick variation and we can fix the point where to define a sample-dependent pseudo-critical point, .
iii.3 Numerical results
iii.3.1 Comparison of the two definitions
In the doubling method we have calculated pseudo-critical points by using both ratios. We have observed, that for a given sample calculated from the ratio of the gaps is always somewhat smaller, than , which is obtained from the ratio of the magnetizations. This is illustrated in the upper inset of Fig.1 for a given sample. We have also calculated for several realizations the ratio of the two pseudo-critical points, , which is shown in the lower inset of Fig.1 as a function of for the leg ladder for various lengths, and . The relative difference between the two pseudo-critical points is indeed vary small, it is of the order of and this is decreasing with increasing and . In the following we restrict ourselves to those pseudo-critical points, which are calculated from the ratio of the magnetization and which have a relative precision of for each sample.
iii.3.2 Distribution of finite-size critical points
We have calculated the distribution of pseudo-critical points for ladders with a fixed number of legs, , for varying lengths, , with . Indeed, for the largest values of the relation, is well satisfied. For the leg ladder the distribution of the values for various lengths are shown in Fig.2, which are obtained for realizations for each cases. As seen in this figure the width of the distribution is decreasing with increasing and there is only a weak shift of the position of the maximum. The distributions somewhat deviate from Gaussians, they are asymmetric, as can be seen in the log-lin plot in the inset of Fig.2. With increasing , however, the skewness of the distribution is decreasing, which is in agreement with the expectation, that in the limit we get back the corresponding results for chains.
iii.3.3 ”True” critical points for ladders
For a fixed value of the number of legs, , we have calculated the mean value of the pseudo-critical points. We have observed that the -dependence of becomes weaker and weaker with increasing , which is in agreement with the fact, that the system approaches more and more the chain geometry. Due to this one can obtain accurate estimates in the thermodynamic limit for the ”true” critical points of ladders, which are listed in Table 1 for different number of legs. Here the errors are merely due to disorder fluctuations since for the finite length effects are negligible.
Here we also list our estimate for the chain, , which agrees within the error of the calculation with the exact result: and (see Sec.III.3.4).
These data approach the critical point in the 2d system, , see Fig.3. Here the corrections for large- are expected to have a power-law form, and analogously to Eq.(1), it contains the shift exponent, , of the 2d system.
Estimates for the effective (-dependent) values of the shift exponent are obtained from the ratio of the second and the first finite differences:
which are calculated at the central point of five-point fits. The effective exponents are given in the upper inset of Fig.3, which are extrapolated as , thus we obtain:
iii.3.4 Scaling of the width of the distribution
We have measured the standard deviation of the distribution of the pseudo-critical points, , for ladders with legs and with a varying length, . This quantity is expected to scale with the length as:
where the scaling function, , behaves for small arguments as: . From this follows, that for finite widths we have:
with a prefactor, , which behaves for large as , with an exponent . We have checked this scenario by analyzing the data for . First, for a fixed we have fitted a function , with a free parameter, . We have found that for each widths, , the exponent agrees with , within a few thousands of error, as illustrated in Fig.4.
In the next step we have fixed the value of and estimated the limiting value of for large , which is denoted by . These limiting values are presented in Table 1, which are analyzed for large . As seen in the upper inset of Fig. 4 in a log-log plot the values are asymptotically on a straight line. We have calculated effective, -dependent exponents: , which are presented in the lower inset of Fig.4 as a function of . These effective exponents have a weak -dependence and we estimate its limiting value as . With this we have for the width exponent in :
Closing this section we try to decide in a direct way about the relation between the two exponents, and . For this we form the scaled difference: (see Eq.(14)), which scales as , as well as the scaled standard deviation: , which scales as , and form their ratio, . As seen in the lower inset of Fig.3 this ratio approaches a finite value which can be estimated as . Thus we can conclude that at the infinite disorder fixed point of the random transverse-field Ising model the shift and the width exponents are equal and they correspond to the correlation length exponent of the model.
Iv Scaling at the critical point
Having estimates for the critical points of random ladders with legs, , we have calculated scaling of the magnetization at the critical point as well as the critical dynamical scaling. These calculations are made for lengths up to and for realizations.
We have calculated the average total magnetic moment at the critical point, , for varying lengths, , which is expected to scale as:
with a scaling function, which behaves for small arguments as: . Then, for a finite width, , we have:
with a prefactor, which for large behaves as: , with .
The scaling Ansatz in Eq.(21) is checked in Fig.5. Then, we have calculated the limiting value of , which is denoted by and which is presented as a function of in a log-log plot in the lower inset of Fig.5. Effective, -dependent exponents are calculated, which are extrapolated in the upper inset of Fig.5 giving . Thus the fractal dimension in is and from Eq.(10) we obtain for the magnetization scaling dimension:
iv.2 Dynamical scaling
At an infinite disorder fixed point there is a special form of dynamical scaling, as given in Eq.(4). The energy scale of a sample at the end of the renormalization can be defined either by the value of the last decimated (log) coupling or by the last decimated (log) transverse field . The distribution of as well as are shown in Fig.6 in upper and in the lower panel, respectively, for the -leg ladder for various values of the length, . An appropriate scaling collapse of the date is observed in terms of the scaling variables, and , with , as illustrated in the insets.
In order to have a more quantitative picture about dynamical scaling we consider the mean value: and the standard deviation, and similarly, and . All these quantities are expected to scale in the same way, for example with we have:
with . For a finite width, , we have then:
with for large with .
We have checked that the scaling form in Eq.(24) is indeed satisfied for all values of and than calculated the limiting value of , which is denoted by . As illustrated in Fig.7 the scaling functions of the typical energy-scales have only a very weak dependence, and we estimate (not shown) a small exponent: . Thus we have for the exponent in the model:
V Entanglement entropy
In the ladder geometry we consider a block, , which contains all the legs and has a length, . Consequently the block has two parallel lines of width, , at which it has contact with the rest of the system, . The entanglement of with is quantified by the von-Neumann entropyAmicoetal08 ():
in terms of the reduced density matrix: , where denotes a pure state (in our case the ground state) of the complete system.
At the critical point of a random quantum system the properties of which are controlled by an infinite disorder fixed point the asymptotic behavior of the entropy in the large limit can be obtained by the strong disorder renormalization group method. For the random transverse-field Ising model entanglement between and are given by such renormalized spin clusters, which contain sites both in and in , and the cluster is eliminated at some point of the renormalizationrefael (); lin07 (); yu07 (). Due to the very broad distribution of the effective couplings and transverse fields, the cluster at the energy scale of its decimation is in a so called GHZ entangled state of the form: . Each such cluster contribute by an amount of to the entanglement entropy, thus calculation of the entropy is equivalent to a cluster counting problem, which is illustrated in Fig.8.
In the chain geometry the asymptotic behavior of the entropy at the critical point is obtained from the analytical solution of the RG equations asrefael ():
where is a non-universal constant, which depends on the form of the disorder, whereas the prefactor of the logarithm, , which is also called as the effective central charge, is universal and given by: . This result is checked numerically in Ref.[il08, ]. In the two-dimensional case, which is expected to hold for , there are somewhat conflicting numerical results at the critical point. Lin et al.lin07 () have observed a double-logarithmic multiplicative factor to the area-law: whereas later Yu et al.yu07 () argued to have only a subleading logarithmic term to the area law: .
Here we study numerically the critical ladder systems with various number of legs and try to identify the cross-over between one- and two dimensions. To illustrate the dependence of the entanglement entropy in Fig.9 we show as a function of for different number of legs for . (We have checked, that the asymptotic results does not change for .) The central parts of the curves are very well linear having approximately the same slope, which is consistent with the exact result for the chain geometry. Thus we conclude that the effective central charge, , does not depend on the number of legs.
In the next step we fix , calculate the non-universal term: and take its limit, , for large (but still with ). As illustrated in the upper inset of Fig.9 the -dependent correction term is approximated as and the asymptotic non-universal terms, , are shown for different number of legs in the lower inset of Fig.9. One can see that starting with the chain, , first is decreasing, has a minimum around and then starts to increase. This increase for large is approximately linear, we have fitted: . This linear increase is compatible with the area law, which should hold for non-critical systems and for large blocks. Our analysis can be used to clarify the one- to two-dimensional cross-over of the entropy in the limit . However, our data can not be used to make predictions further, for , i.e. for the two-dimensional case. For this one should analyze the occurrent and possibly very weak dependence of the prefactor of the linear term of , which however can not be done with our data, which are only up to .
In this paper we have studied the critical properties and the entanglement entropy of random transverse-field Ising models in the ladder geometry by the strong disorder renormalization group method. In our numerical calculation we went up to w=20 legs and with a length up to L=4096 for realizations. In principle the sizes of the systems could have been increased further, but it was not necessary. With we have already reached the limit where no further systematic finite-size effects are seen. On the other hand for larger values of we would have obtained too large errors in calculating quantities, such as through two-point fit.
First, we have calculated sample dependent finite-size critical points, which are obtained by the doubling procedure and the strong disorder renormalization group method. We have analyzed the shift of the mean value of the transition points and the width of the distribution as a function of the number of legs, , and estimated the exponents of the model, and , respectively. These are found to be identical and given by the correlation-length exponent of the model. Consequently the scaling behavior of the pseudo-critical points of the random transverse-field Ising model is in the same form as for classical and conventional random fixed pointsaharony (). In this respect there is a difference with the modelilrm (), in which . For this latter model probably the free-fermionic character could be the reason for the different scaling properties. Our estimate for the correlation length exponent, , is clearly larger than the possible limiting value of , which has been observed in the model and in some other random systemspsz ().
Scaling at the critical point for different quantities is analyzed in a similar way, what we summarize here as follows. Let us consider a physical observable, , which at the critical point has the mean value, . This quantity scales with the critical exponent of the model, , as:
where the scaling function, , for small arguments behaves as:
where is the critical exponent in the model. Consequently for a finite , but for , we have
with and . In general we measure the scaling function for different widths, estimate the exponent and calculate the critical exponent in as: . Since the exponents in are exactly known and the correction term, , is comparatively small we have obtained quite accurate exponents in . In the following we compare the estimates for the different critical exponents in the infinite disorder fixed point, which are listed in Table 2.
Here besides different numerical strong disorder renormalization group results there are also Monte Carlo simulations, both for the random transverse-field Ising model and for the random contact process. This latter model is expected to belong to the same universality classhiv (), at least for strong enough disorder. It is seen in Table 2 that our estimates fit to the trend of the previous results and generally have a somewhat smaller error.
We have also studied the scaling behavior of the entanglement entropy in the ladder geometry. For a fixed width, , the entropy is found to grow logarithmically with the length of the block, , and the prefactor is found independent of . On the other hand the independent term of the entropy is found to have a linear dependence, at least for large enough , which corresponds to the are law for this systems.
The investigations presented in this work can be naturally continued for larger and larger widths and approaching the case, , which corresponds to the two-dimensional model. However, with increasing the numerical computation becomes more and more costly. The reason for this is the fact that the connected clusters in the strong disorder renormalization group method are typically of size , which for large becomes fully connected after decimating a small percent of the transverse fields. The number of further renormalization steps grows in a naïve approach as , so that by this method one can not go further than or in . To treat larger systems improved algorithms are necessary. Studies in this direction are in progress.
Acknowledgements.This work has been supported by the Hungarian National Research Fund under grant No OTKA K62588, K75324 and K77629 and by a German-Hungarian exchange program (DFG-MTA). We are grateful to P. Szépfalusy, H. Rieger and Y-C. Lin for useful discussions.
- (1) L. B. Freund and S. Suresh, Thin film materials, (Cambridge University Press, Cambridge, 2004).
- (2) C. F. Majkrzak, J. Kwo, M. Hong, Y. Yafet, D. Gibbs, C. L. Chien, and J. Bohr, Adv. Phys. 40, 99 (1991).
- (3) E. Dagotto and T. M. Rice, Science 271, 618 (1996).
- (4) M. E. Fisher and M. N. Barber, Phys. Rev. Lett. 28, 1516 (1972).
- (5) M. N. Barber, Phase Transitions and Critical Phenomena Vol. 8 [eds. C. Domb and J. L. Lebowitz] 146 (Academic Press, London, 1983).
- (6) J. Kogut, Rev. Mod. Phys. 51, 659 (1979).
- (7) S. Wiseman and E. Domany, Phys. Rev. Lett. 81 (1998) 22; Phys Rev E 58 (1998) 2938.
- (8) A. Aharony, A.B. Harris and S. Wiseman, Phys. Rev. Lett. 81 (1998) 252.
- (9) A. B. Harris, J. Phys. C 7 , 1671 ( 1974)
- (10) J. T. Chayes et al., Phys. Rev. Lett. 57, 299 (1986).
- (11) F. Pázmándi, R.T. Scalettar, and G.T. Zimányi, Phys. Rev. Lett. 79, 5130 (1997).
- (12) K. Bernardet, F. Pazmandi and G. Batrouni Phys. Rev. Lett. 84 4477 (2000).
- (13) M. T. Mercaldo, J-Ch. Anglès d’Auriac, and F. Iglói Phys. Rev. E 69, 056112 (2004).
- (14) C. Monthus and T. Garel, Eur. Phys. J. B 48, 393-403 (2005).
- (15) For reviews, see: H. Rieger and A. P Young, in Complex Behavior of Glassy Systems, ed. M. Rubi and C. Perez-Vicente, Lecture Notes in Physics 492, p. 256, Springer-Verlag, Heidelberg, 1997; R. N. Bhatt, in Spin glasses and random fields A. P. Young Ed., World Scientific (Singapore, 1998).
- (16) D.S. Fisher, Phys. Rev. Lett. 69, 534 (1992); Phys. Rev. B 51, 6411 (1995).
- (17) C. Pich, A.P. Young, H. Rieger and N. Kawashima, Phys. Rev. Lett. 81 5916 (1998).
- (18) O. Motrunich, S.-C. Mau, D.A. Huse and D.S. Fisher, Phys. Rev. B61 1160 (2000).
- (19) Y.-C. Lin, N. Kawashima, F. Iglói and H. Rieger, Progress in Theor. Phys. 138, (Suppl.) 479 (2000).
- (20) D. Karevski, Y-C. Lin, H. Rieger, N. Kawashima and F. Iglói, Eur. Phys. J. B 20 267-276 (2001).
- (21) F. Iglói, Y-C. Lin, H. Rieger, and C. Monthus, Phys. Rev. B76, 064421 (2007)
- (22) F. Iglói and C. Monthus, Physics Reports 412, 277, (2005).
- (23) T. Vojta, A. Farquhar and J. Mast, Phys. Rev. E79, 011111 (2009).
- (24) L. Amico, R. Fazio, A. Osterloh and V. Vedral, Rev. Mod. Phys. 80 517 (2008).
- (25) G. Refael and J. E. Moore, Phys. Rev. Lett. 93, 260602 (2004).
- (26) Y-C. Lin, F. Iglói and H. Rieger, Phys. Rev. Lett. 99:147202 (2007).
- (27) R. Yu, H. Saleur and S. Haas, Phys. Rev. B77:140402 (2008).
- (28) S.K. Ma, C. Dasgupta and C.-K. Hu, Phys. Rev. Lett. 43, 1434 (1979); C. Dasgupta and S.K. Ma, Phys. Rev. B22, 1305 (1980).
- (29) P. Pfeuty, Phys. Lett. 72A, 245 (1979).
- (30) F. Iglói, Phys. Rev. B 65, 064416 (2002).
- (31) E. Lieb, T. Schultz and D. Mattis, Annals of Phys. 16, 407 (1961).
- (32) F. Iglói and Y-C. Lin, J. Stat. Mech. P06004 (2008).
- (33) J. Hooyberghs, F. Iglói and C. Vanderzande, Phys. Rev. Lett. 90 100601, (2003); Phys. Rev. E 69, 066140 (2004).