PHYSICAL REVIEW B VOLUME 60, NUMBER 16 15 OCTOBER 1999-II Essential finite-size effect in the two-dimensional XY model S. G. Chung Department of Physics, Western Michigan University, Kalamazoo, Michigan 49008 Received 1 February 1999 The thermodynamics of the two-dimensional 2D XY model is formulated by a transfer-matrix method and analyzed by a density-matrix renormalization group. The finite-size scaling and the beta function of the model are studied by the Roomany-Wyld renormalization-group theory. It is found that the 2D XY model has an essential finite-size effect and the Berezinskii-Kosterlitz-Thouless transition with the critical temperature TBKT 0.892 appears in a finite system of 2000­3000 spins as a massless to massive transition with the effective critical temperature Tc 1.07 0.01. S0163-1829 99 03340-8 The two-dimensional 2D XY model describes classical compared to 0.98 of canonical MC simulations. The critical planar spins with nearest-neighbor interactions. The Hamil- temperature 0.894 thus obtained is very close to the seem- tonian is ingly exact value 0.892. In the above, the thermodynamics limit is discussed. In practice, however, when certain physical systems are consid- HXY cos i j . 1 ered by the 2D XY model, the strong finite-size effect is ij essential. In fact, Bramwell and Holdsworth BH Ref. 8 This model has been intensively studied as a model which studied this issue in the context of layered Heisenberg mag- undergoes the Berezinskii-Kosterlitz-Thouless BKT nets with planar anisotropy, which are well regarded as transition,1 the binding-unbinding transition of a vortex- quasi-2D XY systems.9 In these layered magnets, there is a antivortex pair that coexists with spin waves below the tran- very sharp crossover from a 3D behavior to a 2D fluctuation sition temperature. The remarkable point of the BKT theory when lowering temperature. Moreover, the intralayer ex- is that it predicts an essential singularity,2 an exponential change coupling J and interlayer one J have a typical ratio growth of the correlation length and other thermodynamic J/J 103 to 104, and therefore only fluctuations of length quantities near the transition in contrast to the power-law scale less than the order of (J/J )1/2 30­ 100 are two- behavior in a second-order transition. This prediction has dimensional. These layered magnets are thus well described been tested numerically and reported confirmed again and as a finite-size, 103­ 104 spins or sites, 2D XY model. Using again.3 However, there are still some legitimate questions. In the linearized RG equations for finite-size scaling,10 BH fact, there are three or four unknown fitting parameters in found an effective transition temperature Tc 1.080 0.004 the theory. Moreover, it was pointed out that critical region for 1024 spins and a power-law behavior of magnetization to which the renormalization-group RG equations confi- with critical exponent 0.23, explaining some experimen- dently apply is very narrow, (T TBKT)/TBKT 10 2.4 How- tal findings.9 A remarkable point is that 0.23 is universal, ever, due to a strong finite-size effect, ln L L is the system irrespective of lattice types, spin values, and degree of planar size dependence of quantities, approaching such a narrow anisotropy of layered magnets, and claimed to be a signature temperature region by standard Monte Carlo MC simula- of 2D XY behavior. While experiments are nicely explained, tions is almost impossible.3 the BH theory may be questionable in the following points. Three methods have been proposed for handling this First, it relies on the BKT theory that the correlation length finite-size problem. One is Olsson's self-consistent boundary behaves as2 condition method.5 This method can virtually eliminate finite-size fluctuations. Another is the matching method of exp , 2 Hasenbusch and Pinn.6 This method is based on the exact c T TBKT 1/2 solution available for the BCSOS body-centered solid-on- solid model. Comparing the block-spin RG flow of the XY where c is a constant of order 2. However, strictly speaking,4 model with that of the BCSOS model at long distance the it is correct only in a narrow temperature window, (T two models belong to the same universality class , the TBKT)/TBKT 10 2, and its use outside this region needs a method introduces systematic errors which decay as L 2, justification. The effective temperature Tc is defined in BH thereby overcoming the strong finite-size problem. These as the temperature when L. Second, the ``shifted'' tem- two methods remarkably agree on the critical temperature, perature giving TBKT 0.892. The third and recent attempt is a short- 2 time dynamic approach due to Zheng, Schulz, and Trimper.7 T* L TBKT Since the critical exponents at equilibrium enter into the 4c lnL 2 3 short-time dynamic scaling, and the nonequilibrium spatial plays an important role in BH. This is the temperature where correlation length is small in the short-time dynamic evolu- the renormalized spin-wave stiffness in the finite system tion, the method can handle lower temperatures, up to 0.94 takes the universal value 2/ . This is determined by Monte 0163-1829/99/60 16 /11761 4 /$15.00 PRB 60 11 761 ©1999 The American Physical Society 11 762 S. G. CHUNG PRB 60 Carlo as the temperature where the magnetization becomes a where 0 is the lowest energy. The last equality holds for spin-wave estimate at T TBKT , N M 1. The free energy per site is given by M 1 1/16 f T 0 2L2 . 4 M . 12 The argument seems less precise and again uses the BKT The correlation function can be evaluated likewise,11 theory beyond its confirmed validity. In this paper, we study the finite-size effect in the 2D XY 0 r 0 0 r 0 , 13 model by a transfer matrix TM method11 with the help of where the density-matrix RG DMRG Ref. 12 technology and 0 denotes the eigenvector for the eigenvalue 0 . The problem is thus reduced to the TM eigenvalue problem Roomany-Wyld RG theory.13 Our approach does not involve 10 . any uncontrolled parameters, nor rely on the BKT theory. It To analyze the TM equation 10 , we use DMRG recently demonstrates that the finite-size XY model with 2000­3000 developed for 1d quantum systems.12 In fact, the TM equa- spins behaves precisely like a system with the transition tem- tion contains the Hamiltonian in a Boltzmann factor, so it is perature Tc 1.07 0.01 separating the massless and massive essentially a quantum M-body problem. The TM-DMRG phases in agreement with the Tc 1.080 0.004 in the BH method is now a powerful tool for studying the thermody- theory for 1024 spins. namic properties of 1d quantum systems.14 In the latter case, To calculate the partition function the Suzuki-Trotter decomposition of the Boltzmann factor N leads to a non-Hermitian TM equation, requiring a little bit Z d of numerical innovations.11,15 In the present case, the TM i exp HXY , 5 i 1 operator K is real and symmetric, and less troublesome. On the other hand, unlike the quantum 1d cases where the TM where i ( i1 , i2 ,..., iM), consider a ribbon geometry; M equations are discrete from the outset, our TM equation is an sites in the width direction and N sites along the long, circu- integral equation, and needs a descretization. For that pur- lar direction. With a negligible error at the edges, the Boltz- pose, we use a Gaussian integration formula, transforming mann factor in Eq. 5 can be written as Eq. 10 to N p K i , i 1 6 M w i 1 iw jwk¯K i jk¯ :i j k ¯ n i jk¯ i jk¯ 1 with the TM operator exp n n i j k ¯ , 14 M 1 where n(ijk¯) means n( i j k¯) and i (xi 1), K i , i 1 exp etc., with w j 1 2 cos ij i 1j i and xi being the Gaussian abscissas and weight factors.16 If we rescale as cos ij ij 1 cos i 1j i 1j 1 Mwiwkwk¯ 1/2 i jk¯ i jk¯ , 15 cos then the eigenvector is Euclid normalized and the TM equa- i j 1 i 1 j 1 . 7 tion becomes Now write OS ii W ij:i j OS jj W jk:j k N N 1 i jk¯ d 1 d i 1 N 1 8 i 1 i 1 OS kk ¯ i jk¯ and expand as exp n i j k ¯ , 16 where W(i j:i j ) denotes the Boltzmann factor on the right- 1 N 1 n* N 1 n 1 , 9 hand side of Eq. 7 with i 1 i and j 1 j , and n OS(ii ) (wiwi )1/2, etc. where n is an orthonormal complete set. We choose The DMRG procedure to get the lowest eigenvalues n is n such that the TM equation as follows. Step i , consider three sites, that is the p3 p3 eigenvalue problem. Note that, unlike the 1d quantum cases 2 where the degree of freedom at each site is small, the degree d K , n exp n n . 10 0 of freedom here, p, needs to be large, so we have to use the three-site algorithm whose accuracy was discussed before.17 With these manipulations and using the TM equation repeat- Step ii , construct the density matrix out of the ground state edly, we arrive at 0 , Z exp nN exp 0N , 11 i j :i j 0* i j k 0 ijk , 17 n k PRB 60 ESSENTIAL FINITE-SIZE EFFECT IN THE TWO- . . . 11 763 FIG. 3. Beta function vs temperature for L 45, 47, 49, and 51 with L L 2 from top to bottom. Again p 12 and m 60. To- FIG. 1. Internal energy vs temperature. Open circles are data gether with Fig. 2, the effective critical temperature is determined to from the equation of motion method for the kinetic XY model and be Tc 1.07 0.01. diamonds are Monte Carlo data available up to T 1.5 Ref. 18 . Solid circles are from the TM-DMRG. O, and then a RG transformation T2*W T3 . The seven-site system is now represented by T3*T3 . Repeating the proce- diagonalize it, and retain the largest m states. Denote the m dure, one can get a larger and larger system. eigenvectors by X1 ,X2 ,...,Xm and eigenvalues by One needs to keep track of, in the course of RG transfor- 1 , 2 ,..., m . So (X1X2¯Xm) O(ij: ) constitutes a p2 mations, the spin operator (cos ,sin ) at each site to calcu- m transformation matrix for the change of basis late various expectation values such as the correlation func- . The renormalized TM operator is tion 13 . For example, in the RG operation *W*W T2 in steps i and ii , the angle the caret denotes an operator O* :i j OS ii in the middle site is a one-particle operator, diagonal in the ii jj original angle representation, and therefore represented as W i j:i j OS j j W jk: j k O i j: 1 i j: i j ii j jj . 19 T2 k : k . 18 Thus the RG transformation new is Symbolically *W*W T2 , where the solid circles represent OS. Step iii , go back to step i with the renormalized op- erator T O :i j 2*T2 for the five sites, and diagonalize the resulting ii j jj O i j: new . 20 pm2 pm2 eigenvalue problem to get the new ground state ii jj 0 . Step iv , repeat step ii to construct the reduced den- Let us first check the accuracy of our TM-DMRG method. sity matrix, diagonalize it to get a new transformation matrix For this purpose, we calculate the internal energy per site E T2 T 0 /M . 21 It is now well established that the ground-state energy per site 0 /M in the thermodynamic limit can be easily obtained by measuring the energy increase in successive steps of the DMRG procedure. In the present case, the typical system size to get the convergence is 30. We have started the calcu- lation with p m 12 and checked the complete conver- gence up to p m 24. In the T 0.1 case, we have checked the convergence by calculating up to the p m 32 case. The result is given in Fig. 1, where the kinetic energy T/2 is added to compare with the latest simulation on the kinetic XY model.18 The agreement is perfect near T 1, but both at low T 0.1 and high T 1.3 the TM-DMRG result is about 10% larger than the simulation result. Since the finite-size effect is FIG. 2. L gap(L) vs temperature for L 37, 43, 49, 55, and 61 negligibly small at higher T, the discrepancy is a puzzle at from bottom to top. The result is for p 12 and m 60. The cases the moment. However, similar discrepancies have been ob- p 12, m 66 and p 16, m 60 give less than 1% correction to served for the specific heat among MC calculations,3 so the the result. problem could reside in simulations. 11 764 S. G. CHUNG PRB 60 We now calculate not only the smallest 0 but also the temperature is determined to be Tc 1.07 0.01, in agree- next smallest 1 . In fact, one can show easily that the cor- ment with the Tc 1.080 0.004 in the BH theory8 for 1024 relation function 13 behaves for r 1 like spins. To conclude, the thermodynamics of the 2D XY model is exp 1 0 r . 22 formulated by a transfer-matrix method and analyzed by a Thus gap(L) 1 0 measures the gap energy of the sys- density-matrix renormalization group. The finite-size scaling tem at size L alias M in Eq. 21 . Then due to the and the beta function of the model are studied by the Roomany-Wyld RG theory,13 the effective critical tempera- Roomany-Wyld renormalization-group theory. It is found ture Tc of the finite-size 2D XY model can be located by that the 2D XY model has an essential finite-size effect and calculating L gap(L) versus T. Figure 2 shows the result the Berezinskii-Kosterlitz-Thouless transition with the criti- for L 37, 43, 49, 55, and 61 using p 12 and m 60. cal temperature TBKT 0.892 manifests itself in a finite sys- The cases p 16, m 60 and p 12, m 66 give less than tem of 2000­3000 spins as a massless to massive transition 1% correction to the result. We have also calculated the with the effective critical temperature T Roomany-Wyld approximate for the beta function c 1.07 0.01. After the completion of this work, I was informed of a 1 ln G new finite-size-scaling MC method which gives the critical RW L /GL /ln L/L temperature T LL T D BKT 0.893, very close to the seemingly exact LDL /GLGL 1/2 , 23 value 0.892, and demonstrates that the critical exponents where GL gap(L) and DL GL / T. Figure 3 shows the 1/4 and 0.48.19 result for L 45, 47, 49, and 51 with L L 2. Again p 12 and m 60. Small drift of the zero point of the beta This work was partially supported by the NSF under function towards higher T Grant No. DMR 980009N and utilized the SGI/CRAY Ori- c for larger L may be due to an accumulation of the truncation (p m m) error in DMRG. gin 2000 at the National Center for Supercomputing Appli- Based on the results in Figs. 2 and 3, the effective critical cations at the University of Illinois at Urbana-Champaign. 1 Z. L. Berezinskii, Zh. Eksp. Teer. Fiz. 59, 907 1970 Sov. Phys. 10 J. M. Kosterlitz, J. Phys. C 7, 1046 1974 ; K. Binder, Monte JETP 32, 493 1971 ; J. M. Kosterlitz and D. J. Thouless, J. Carlo Methods, Topics in Current Physics Vol. 7 Springer, Ber- Phys. C 6, 1181 1973 . lin, 1986 . 2 C. Itzykson and J.-M. Douffe, Statistical Field Theory Cam- 11 S. G. Chung, D. R. A. Johlen, and J. Kastrup, Phys. Rev. B 48, bridge University Press, New York, 1989 , Vol. 1. 5049 1993 ; S. G. Chung, J. Phys.: Condens. Matter 9, L219 3 J. Tobochnik and G. V. Chester, Phys. Rev. B 20, 3761 1979 ; R. 1997 . Gupta and C. F. Baillie, ibid. 45, 2883 1992 ; J. R. Lee and S. 12 S. R. White, Phys. Rev. B 48, 10 345 1993 . Teitel, ibid. 46, 3247 1992 ; N. Schultka and E. Manousakis, 13 H. H. Roomany and H. W. Wyld, Phys. Rev. D 21, 3341 1980 . ibid. 49, 12 071 1994 ; M. Campostrini, A. Pelisetto, P. Rossi, 14 X. Wang and T. Xiang, Phys. Rev. B 56, 5061 1997 ; N. Shibata, and E. Vicari, ibid. 54, 7301 1996 . J. Phys. Soc. Jpn. 66, 2221 1997 . 4 J. M. Greif, D. L. Goodstein, and A. F. Silva-Moreira, Phys. Rev. 15 B. Ammon, M. Troyer, T. M. Rice, and N. Shibata, Phys. Rev. B 25, 6838 1982 ; J. L. Cardy, ibid. 26, 6311 1982 ; P. Min- Lett. 82, 3855 1999 . nhagen, Rev. Mod. Phys. 59, 1001 1987 . 16 Handbook of Mathematical Functions, edited by M. Abramowitz 5 P. Olsson, Phys. Rev. Lett. 73, 3339 1994 . and I. A. Stegun, Natl. Bur. Stand. U.S. Applied Mathematics 6 M. Hasenbusch and K. Pinn, J. Phys. A 30, 63 1997 . Series No. 55 U.S. GPO, Washington, D.C., 1972 , Chap. 25. 7 B. Zheng, M. Schulz, and S. Trimper, Phys. Rev. E 59, 1351 17 S. G. Chung, J. Phys.: Condens. Matter 9, L619 1997 ; S. G. 1999 . Chung, Current Topics in Physics, edited by Y. M. Cho, J. B. 8 S. T. Bramwell and P. C. W. Holdsworth, J. Phys.: Condens. Hong, and C. N. Yang World Scientific, Hong Kong, 1998 , Matter 5, L53 1993 . Vol. 1, p. 295. 9 L. P. Regnault and J. Rossat-Mignod, in Magnetic Properties of 18 X. Leoncini, A. D. Verga, and S. Ruffo, Phys. Rev. E 57, 6377 Layered Transition Metal Compounds, edited by L. J. de Jongh 1998 . Kluwer, Deventer, 1990 . 19 J.-K. Kim, Phys. Lett. A 223, 261 1996 .