PHYSICAL REVIEW B 77, 155413 共2008兲
Rapidly converging methods for the location of quantum critical points from finite-size data M. Roncaglia,1 L. Campos Venuti,2 and C. Degli Esposti Boschi3 1Max-Planck-Institut
für Quantenoptik, Hans-Kopfermann-Strasse 1, D-85748 Garching, Germany ISI, Villa Gualino, viale Settimio Severo 65, I-10133 Torino, Italy 3 CNR, Unità di Ricerca CNISM and Dipartimento di Fisica dell’Università di Bologna, viale Berti-Pichat 6/2, I-40127 Bologna, Italy 共Received 21 January 2008; revised manuscript received 12 March 2008; published 9 April 2008兲 2Fondazione
We analyze in detail, beyond the usual scaling hypothesis, the finite-size convergence of static quantities toward the thermodynamic limit. In this way, we are able to obtain sequences of pseudo-critical points, which display a faster convergence rate as compared to currently used methods. The approaches are valid in any spatial dimension and for any value of the dynamic exponent. We demonstrate the effectiveness of our methods both analytically, on the basis of the one dimensional XY model, and numerically, considering c = 1 transitions occurring in nonintegrable spin models. In particular, we show that these general methods are able to precisely locate the onset of the Berezinskii–Kosterlitz–Thouless transition making only use of ground-state properties on relatively small systems. DOI: 10.1103/PhysRevB.77.155413
PACS number共s兲: 05.70.Fh, 64.60.an, 75.40.Mg
I. INTRODUCTION
In the study of physical properties of phase transitions, a basic prerequisite is a reliable method to locate the critical point, whenever the latter is not known a priori from symmetry or duality arguments. Typically, in numerical or even experimental studies on finite samples, one obtains a sequence of pseudocritical points 共in the sense specified below兲 to be extrapolated to the true critical point in the thermodynamic limit 共TL兲. The extrapolation may be done with some polynomial fit in the inverse size of the system or, better, exploiting some fitting function derived on the basis of a scaling Ansatz or through the renormalization group 共RG兲. The point is especially relevant in the context of quantum phase transitions1 共QPTs兲 in lattice systems where the exponential growth of the dimension of the Hilbert space with the number of sites is a strong limitation on the accessible sizes with the current computational power and algorithms. One of the most used algorithms is still the Lanczos method for the virtually exact extraction of the low-lying energy levels; in the most favorable case of spin-1/2 models, one cannot go beyond some tens of sites. This limit can be moved to maybe a few thousands of sites by using the so-called density matrix renormalization group2 共DMRG兲 that has become the method of choice for one dimensional 共1D兲 problems due to its high level of accuracy. Nonetheless, if one considers two or even three dimensional systems, the situation is much worse: with the Lanczos algorithm, the largest lattices have only a few sites of linear extension and the DMRG is not particularly efficient. At present, the only other choice is quantum Monte Carlo 共QMC兲 共see, for instance, Ref. 3兲 that, however, suffers from a sign problem in the case of fermionic or frustrated systems and does not reach the level of accuracy of the DMRG. Very recently, there have been attempts to exploit both DMRG-like features and the QMC sampling tricks to design hybrid methods4,5 that are, however, still under verification. It is generally believed that a sequence of pseudocritical points, for example, the loci of maxima of finite-size susceptibilities, converges to the critical point as a power law L− 1098-0121/2008/77共15兲/155413共9兲
with a so-called shift exponent given by the inverse of the correlation length exponent . Hence, generally speaking, the larger is the , the slower is the convergence. This difficulty reaches its maximum for Berezinskii–Kosterlitz–Thouless 共BKT兲 transitions, in which the correlation length diverges with an essential singularity or, loosely speaking, “ = ⬁.” However, already in the seminal paper by Fisher and Barber,6 it was pointed out that the relation = 1 / is not always valid and depends, among other factors, on the boundary conditions. The most used method to locate quantum critical points in d = 1 by means of finite-size data is the so-called phenomenological renormalization group 共PRG兲, which is reviewed, for instance, in Ref. 7. Another convenient approach, which is the finite-size crossing method 共FSCM兲, was recently proposed in Ref. 8. The aim of this paper is to improve both of them by means of criteria that produce sequences of pseudocritical points that converge more rapidly. We will show that in our sequences, the shift exponent will have the form = 0 + 1 / , where 0 and/or 1 are larger than the corresponding values in the usual methods and, therefore, allow for a better convergence. The paper is organized as follows. In Sec. II, we illustrate the general arguments leading to the enhanced sequences, both in the framework of the FSCM and of the PRG 共Sec. II D兲. Special cases as the BKT transition 共Sec. II C兲 and that of logarithmic divergences 共Sec. II B兲 are separately discussed. In Sec. III, we illustrate the usefulness of the methods on the hand of analytic and numerical tests. In Sec. III A, we treat the XY spin-1/2 chain by using a series of exact calculations reported in the Appendix. Then, we move to two cases of spin chains for which no exact solution is available: in Sec. III B, we consider a spin-1 model with anisotropies in a parameter range that gives rise to a large value of , and in Sec. III C, we study the spin-1/2 model with next-to-nearest neighbor interactions that is known to undergo a BKT transition. In this case, we find a value for the critical coupling in agreement with the accepted one, which was found by using a model-specific investigation of the excited states.9 Section IV is devoted to conclusions.
155413-1
©2008 The American Physical Society
PHYSICAL REVIEW B 77, 155413 共2008兲
RONCAGLIA, CAMPOS VENUTI, AND DEGLI ESPOSTI BOSCHI II. DERIVATION OF RAPIDLY CONVERGING SEQUENCES
We consider systems in d spatial dimensions of linear size L and periodic boundary conditions 共PBCs兲. Let the transition be driven by a linear parameter g such that the Hamiltonian is H共g兲 = H0 + gV. Dealing with QPT, we consider the case of strictly zero temperature, T = 0, even if the arguments presented below can be simply extended to the finite-temperature case, replacing the parameter g with T 共and without using the dimensional crossover rule used below兲. The free energy density reduces to the ground-state 共GS兲 energy density, which, close to the critical point gc, shows a singularity in the second 共or higher兲 derivatives with respect to g as follows: 1 具H共g兲典 = e共g兲 = ereg共g兲 + esing关共g兲兴, Ld where ⬇ 兩g − gc兩− ⬅ t− is the correlation length. Note that as a consequence of the scaling hypothesis, the singular part of the energy esing is a universal quantity that depends only on , which is the relevant length scale close to the critical point. Hence, esing may be considered quite in general an even function of 共g − gc兲 that vanishes at the critical point. On the other hand, the bulk energy density at the critical point behaves as 共Privman–Fisher hypothesis兲 e共gc,L兲 = e⬁共gc兲 − L
−共d+兲
F共gc兲,
共1兲
where F共g兲 is a sort of a Casimir-type term that may depend on the actual geometry of the lattice. Note that this hypothesis has to be properly changed if one or more of the spatial dimensions are of infinite extent. Moreover, Eq. 共1兲 has been written in analogy with Eq. 共11.29兲 of Ref. 10 by using the dimensional crossover rule according to which the partition function and the thermodynamic 共static兲 properties of a d-dimensional quantum system are equivalent to those of a 共d + 兲-dimensional classical counterpart,1,10 where is the dynamic exponent.1 Then, for the implementation of our methods, we need to know by some other means the value of relating the energy gap ⌬ and the correlation length : ⌬ ⬀ −. Typically, but not always, energy and momentum in the continuum limit at the critical point satisfy a linear dispersion relation, E = vk, for small k so that = 1 and a relativistic effective field theory can be used to describe the universal features of the transition. In d = 1 the scale invariance at the critical point is often sufficient to also imply conformal invariance 共see Chap. 2 of Ref. 11兲, thanks to which several exact results can be obtained by using the powerful predictions of conformal field theories 共CFTs兲. For example, by mapping the space-time complex plane onto a cylinder whose circumference represents the finite chain of length L, we can identify F共gc兲 = cv共gc兲 / 6 in Eq. 共1兲, where c is the central charge of the theory. In the RG sense, moving away from criticality corresponds to perturbing the CFT with a relevant operator that destroys conformal invariance. However, this is not the only effect of varying the microscopic parameter g out of gc: in general, also the speed of elemen-
tary excitations gets renormalized in the unperturbed CFT part. For this reason, we say that v 共henceforth F兲 depends on g in the vicinity of gc. Scaling and dimensional arguments imply that in the thermodynamic off-critical regime L Ⰷ , the singular part of the energy behaves as esing ⬇ t2−␣ ,
共2兲
where ␣ = 2 − 共d + 兲. For a second order phase transition, ␣ ⬍ 1. After introducing the scaling variable z = 共L / 兲1/ = tL1/, the finite-size scaling 共FSS兲 theory asserts7 that in a system of length L, esing = C0L−共d+兲⌽e共z兲 + ¯ ,
共3兲
where ⌽e共z兲 is a universal function that, in the off-critical regime z Ⰷ 1, must behave as ⌽e共z兲 ⬇ z共d+兲 in order to recover Eq. 共2兲. Instead, for L Ⰶ we are in the critical regime and ⌽e共z兲 behaves as an analytic function that vanishes for z → 0. Here, we assume that the leading term in ⌽e共z兲 is quadratic in z 共see below why it cannot be linear兲, but the following arguments are easily generalizable to higher integer powers. The constant term C0⌽e共0兲 / Ld+ is already adsorbed in the nonuniversal part of the energy density at gc, as shown in Eq. 共1兲. Differentiating e共g兲 with respect to g gives the mean value b = 具V典 / Ld, whose singular part bsing behaves as bsing ⬇ sgn共g − gc兲t1−␣ = sgn共g − gc兲t共d+兲−1 .
共4兲
Considering FSS for the combination of Eqs. 共1兲 and 共3兲 and then differentiating, we find b共g,L兲 = b⬁,reg共g兲 + sgn共g − gc兲C0L1/−共d+兲⌽e⬘共z兲 − L−共d+兲关F⬘共gc兲 + F⬙共gc兲共g − gc兲 + ¯兴 + O共L−共d++⑀兲兲,
共5兲
where the subscript “⬁, reg” hereafter means regular in the TL. In order to write down the expression above, we used g = L1/z and assumed that the powers neglected in the last term are just larger than 共d + 兲. To illustrate this point, we could consider the irrelevant operator with the smallest scaling dimension dirr. At first order in perturbation theory with the renormalized coupling girr共L兲 = girr共0兲Ld+−dirr, the corrections to the GS energy density are of the form Cirrgirr共0兲L−dirr so that ⑀ = dirr − d − ⬎ 0. Note also that the amplitude Cirr can 2 兲 terms have to be included. For these and vanish and O共girr more details, we leave the reader to Ref. 12. Generically, we admit corrections with ⑀ ⱖ 0 that may come either from irrelevant operator in the continuum theory or from lattice effects. The case ⑀ = 0 corresponds to marginal perturbations and typically leads to logarithmic corrections. Notice that now the leading term of ⌽e⬘共z兲 is linear in z in the critical region. If we had admitted a linear term in ⌽e共z兲, then ⌽e⬘共0兲 ⫽ 0 and a finite jump discontinuity at finite L would be present in b共g , L兲. We can also calculate a = 具H0典 / Ld = e共g兲 − g关e共g兲 / g兴, which yields to a singular part that is similar to Eq. 共4兲 but with a changed sign. In fact, for gc ⫽ 0, the leading singular parts of a and b must cancel in the sum that gives back the
155413-2
PHYSICAL REVIEW B 77, 155413 共2008兲
RAPIDLY CONVERGING METHODS FOR THE LOCATION…
Lb共g,L兲 L2 b共g,L兲 . = Le共g,L兲 L2 e共g,L兲
energy e共g兲 = a共g兲 + gb共g兲, which does not contain that singularity. In particular, the scaling is a共g,L兲 = a⬁,reg共g兲 − sgn共g − gc兲gC0L1/−共d+兲⌽e⬘共z兲
Calling = 共g − gc兲 = tsgn共g − gc兲 and rewriting only the essential terms in the scaling Ansätze, we have the simplified forms
+ C0L−共d+兲⌽e共z兲 − L−共d+兲关F共g兲 − gF⬘共gc兲 + O共g − gc兲兴 + ¯ . The FSCM8 identifies the critical point with the limit of the sequence gLⴱ of single crossing points as follows: b共gLⴱ ,L兲 = b共gLⴱ ,L⬘兲,
e共g,L兲 = e⬁,reg共g兲 − L−共d+兲F共g兲 + D1L−共d++⑀兲 , b共g,L兲 = b⬁,reg共g兲 + L−共d+兲关2C0L2/ − F⬘共g兲兴 + D2L−共d++⑀兲 . Now, by putting these two relations in 共共8兲兲, we obtain
where L⬘ = L + ␦L. Applying this criterion to Eq. 共5兲, we obtain 共for ␦L Ⰶ L兲 gLⴱ − gc = −
共d + 兲F⬘共gc兲 L−2/, 2 2C0 − 共d + 兲 v
册
冋
⫽
2 . d+
This equation defines the shift exponent FSCM = 2 / and may converge very slowly when Ⰷ 1. The extremely difficult case is the BKT transition where formally = ⬁, but this latter situation must be treated in a different way 共see Sec. II C兲. Now, consider the quantity ⌫共g , L , ␥兲 = ␥e共g , L兲 − b共g , L兲 and suppose to be able to tune ␥ exactly at
␥ⴱ =
F⬘共gc兲 . F共gc兲
共6兲
It easily seen that ⌫共gc , L , ␥ⴱ兲 does not contain the Casimirtype term responsible for the critical point shift. In fact, the scaling of ⌫ is ⌫共g,L, ␥兲 = ⌫⬁,reg共g, ␥兲 − sgn共g − gc兲C0L1/−共d+兲⌽e⬘共z兲 − L−共d+兲关␥F共g兲 − F⬘共gc兲 + O共g − gc兲兴
gLⴱ − gc =
b共gLⴱ ,L + ␦L兲 − b共gLⴱ ,L兲 e共gLⴱ ,L + ␦L兲 − e共gLⴱ ,L兲
.
Notice that putting ␥ = 0 is equivalent to the FSCM applied to b共g , L兲 and the denominator has a definite sign about the critical point. Now, we find gLⴱ requiring that ␥共gLⴱ , L − ␦L兲 = ␥共gLⴱ , L兲, i.e., − −
b共gLⴱ ,L − ␦L兲 e共gLⴱ ,L − ␦L兲
or in the continuum version
=
b共gLⴱ ,L e共gLⴱ ,L
+ ␦L兲 − + ␦L兲 −
b共gLⴱ ,L兲 , e共gLⴱ ,L兲
共7兲
共9兲
A. Homogeneity condition
In the previous section, we have shown how to obtain a sequence of pseudocritical points with an improved shift exponent fast. Here, we provide another yet simpler equation for the determination of the critical point. The resulting pseudocritical sequence is characterized by the same shift exponent fast. However, in this case, we are able to prove convergence toward the critical point even in the extreme case of a BKT transition 共see Sec. II C兲. The idea is to require that at the critical point the L-dependent part of b is dominated by the Casimir-type term with power 共d + 兲 关see Eq. 共5兲兴. This condition is translated into the requirement that the L-derivative of b is a homogeneous function of degree 共d + + 1兲, i.e., 共d + + 1兲Lb共gLⴱ ,L兲 + LL2 b共gLⴱ ,L兲 = 0.
When ␥ is equal to ␥ⴱ given in Eq. 共6兲, the critical point found by this crossing method is approached as gLⴱ − gc ⬇ L−fast, with a shift exponent fast = 2 / + ⑀ = FSCM + ⑀ 关again this holds true provided that ⫽ 2 / 共d + 兲兴. The additional term ⑀ allows, in general, for a better convergence of the sequence gLⴱ . A possible algorithm for numerically finding the critical point in such a way is the following. If gLⴱ is at a crossing point of ⌫共g , L , ␥兲, then we have
b共gLⴱ ,L兲 e共gLⴱ ,L兲
D2F共gc兲 − D1F⬘共gc兲 ⑀共d + + ⑀兲 − L fast , 4C0F共gc兲 共d + − 2 兲
which gives a shift exponent fast as anticipated. The main result of this section is the crossing criterion 关Eq. 共7兲兴 that identifies the rapidly converging sequence 关Eq. 共9兲兴 to the critical value.
+ O共L−共d++⑀兲兲.
␥共gLⴱ ,L兲 =
共8兲
共10兲
Consequently, the corresponding sequence of pseudocritical points gLⴱ scales as gLⴱ − gc =
⑀D2共d + + ⑀兲 − L fast . 2C0 2 共d + − 2 兲
共11兲
Equation 共10兲 represents the homogeneity condition method 共HCM兲 that we propose for the efficient location of critical points. We stress here that b共g , L兲 being a GS property, this criterion does not require knowledge of excited states as in the case of the PRG method. This is a point in favor to the HCM since excited states are typically assessed with less numerical accuracy. In addition, the HCM is superior to the PRG in that it produces a faster converging sequence 共see Sec. II D兲. B. Case = 2 Õ (d + ) with logarithmic divergences
For completeness, we consider the case = d+2 that was excluded in the previous treatment. In this situation, the Ansatz requires the inclusion of logarithmic corrections as follows:
155413-3
PHYSICAL REVIEW B 77, 155413 共2008兲
RONCAGLIA, CAMPOS VENUTI, AND DEGLI ESPOSTI BOSCHI
b共g,L兲 = b⬁,reg共g兲 + C1 ln L − L−共d+兲F⬘共g兲 + D2L−共d++⑀兲 . The calculation of the critical point with the FSCM gives in this case gLⴱ − gc = −
共d + 兲F⬘共gc兲 −共d+兲 L , C1
whereas the HCM Eq. 共10兲 yields gLⴱ − gc = −
⑀D2共d + + ⑀兲 −共d++⑀兲 L . C1共d + 兲
These results are compatible with the exact calculations of the XY model 共see Sec. III A兲. Note that, in general, we should perform a similar calculation for
=
p , d+
p 僆 N + 1.
Namely, when the pth derivative of the free energy diverges logarithmically.
ence in L in order to drop F⬘共g兲. Formally, in the region L / Ⰶ 1, we can write down the condition
L关L3Lb共gLⴱ ,L兲兴 = 0.
It is worth noting that the latter condition is equal to the HCM 关Eq. 共10兲兴 when d = = 1. Treating L as a continuous variable, one can read off the shift exponent for the sequence BKT = ⑀ / 共n − 1 + 兲; if ⑀ = 0, then the gLⴱ that turns out to be fast corrections to scaling are also governed by 共another兲 marginal operator and we expect gLⴱ − gc ⬃ 共ln L兲−共m+1兲/共n−1+兲, where m is a positive integer 关m = 4 from Eq. 共22兲 in Ref. 12兴. Dealing with numerical simulations, it is very important to specify how one implements the finite-size differences in L. In fact, there are several finite-difference expressions used in the literature to express the derivatives and, here, the requirement is that they all reproduce Eq. 共13兲 in the limit L → ⬁. For example, if one takes a uniform step ␦L, then the following symmetric expression can be built: b⬘共g,L兲 ⬅
C. Scaling Ansatz for the Berezinskii–Kosterlitz–Thouless case
For d = 1 at the BKT, the correlation function in the TL behaves like ⬇ exp共at−兲. In the typical example of the classical two dimensional XY model, it is known that = 1 / 2. Instead, for the quantum Heisenberg model with frustration 共which we will consider in Sec. III C兲 Haldane suggested = 1.9 We also set = 1 because the 共effective兲 dimensionality in the BKT scenario is 2. The singular part of the finite-size energy density now is conveniently expressed in terms of y ⬅ L / ,7 so that e共g,L兲 = e⬁,reg共g兲 + L−2关C0⌽e共y兲 − F共g兲兴 + O共L−共2+⑀兲兲, where ⌽e共y兲 is a universal function that, in the off-critical regime y Ⰷ 1, must behave as ⌽e共y兲 ⬇ y 2. Again, in the quasicritical regime y Ⰶ 1 at any finite L, the energy density and its derivatives must be analytic in 共g − gc兲. The value ⌽e共0兲 can be absorbed in F共g兲 and it can be directly checked that the first contribution has to be at least quadratic in t = 兩g − gc兩 because, otherwise, a finite-size discontinuity in b共g , L兲 would be generated. For y Ⰶ 1, we adopt the following Ansatz 共justified from perturbed conformal field theory12兲: e共g,L兲 = e⬁,reg共g兲 + L 关K共at −2
−
− ln L兲
−n/
− F共g兲兴 + ¯ ,
where K is a constant and n is an integer larger than 1 关n = 3 from Eq. 共22兲 in Ref. 12兴. Hence,
冋
b⬙共g,L兲 ⬅
b共g,L + ␦L兲 − b共g,L − ␦L兲 , 2␦L
b共g,L + ␦L兲 − 2b共g,L兲 + b共g,L − ␦L兲 , 共␦L兲2
and then the precise condition to cancel the term O共L−2兲 becomes L3b⬙共gLⴱ ,L兲 + 关3L2 − 共␦L兲2兴b⬘共gLⴱ ,L兲 = 0.
冉
冊
−共n/兲−1
共14兲
In the limit of large L, we recover Eq. 共13兲 as required. Clearly, the correct discretization prescription must be identified not only for the BKT; when d and/or are not one, the suitable variant of Eq. 共14兲 for finite L has to be adopted with b⬙ weighted by L共d++1兲 and b⬘ weighted by a polynomial in L of degree 共d + 兲, whose coefficients depend on 共␦L兲. D. Phenomenological renormalization group revisited
The PRG method identifies the critical point with the limit of the sequence gLⴱ of crossing points satisfying L⌬共gLⴱ ,L兲 − 共L + ␦L兲⌬共gLⴱ ,L + ␦L兲 = 0.
共15兲
Here, ⌬共g , L兲 is the finite-size energy gap of the spectrum for which we may adopt the following form: ⌬共g,L兲 = Ld关e1共g,L兲 − e共g,L兲兴 = L−⌽共z兲 + L−共+⑀兲C , 共16兲
b共g,L兲 = b⬁,reg共g兲 + L−2 nKa−n/tn−1sgn共g − gc兲 ln L t ⫻ 1− a
共13兲
册
where ⌽共z兲 ⯝ ⌽0 + z⌽1 + z ⌽2 + . . . and C is a prefactor depending on the excited state 兩典 that we consider. The standard PRG approach in d = 1 relies on Eq. 共15兲 with = 1. More generally, a first test to identify is done by plotting the usual scaled gaps L⌬ and see if at the critical point they settle to a constant or not. If they do not, one is led to search for a better value of ⫽ 1 and solve Eq. 共15兲 with the correct value of the dynamic exponent. Note that the corrections come from the irrelevant contributions and the Casimir-type term for the gap is actually the constant term in ⌽共z兲. For 2
− F⬘共g兲 + O共L−共2+⑀兲兲, 共12兲
Now, we want to get rid of all the O共L−2兲 contributions that “hinder” the location of finite-size pseudocritical points. Hence, we first differentiate with respect to L to eliminate b⬁,reg共g兲, then multiply by L3 to isolate the term in square brackets in Eq. 共12兲, and, finally, set to zero a further differ-
155413-4
PHYSICAL REVIEW B 77, 155413 共2008兲
RAPIDLY CONVERGING METHODS FOR THE LOCATION…
quantum systems in d = 1, CFT 共 = 1兲 predicts ⌽0 = 2dv共g兲, where d is the scaling dimension of the operator ˆ that generates the excited state 兩典. In addition, ⌽1 , ⌽2 , . . ., can be computed in the framework of perturbed CFT, in the regime z Ⰶ 1 when a relevant operator 共g − gc兲R with scaling dimension x = 2 − 1 / is added to the critical field theory. For instance, from Eqs. 共7兲 and 共10兲 of Ref. 12 we have ⌽1 = b / 共2兲x, where b is the structure constant that appears as prefactor of the three point correlation function 具R共rជ1兲ˆ 共rជ2兲ˆ 共rជ3兲典. Now, in general, the PRG method gives 兩gLⴱ − gc兩 ⬇ L−PRG, where PRG = 1 / + ⑀.11 However, it may happen that ⌽1 = 0 in the scaling function 共examples are discussed in Sec. III B and at the end of Sec. III A for a specific excited state of the XY model兲. The shift exponent in such a case decreases to PRG = 1 / + ⑀ / 2. Nonetheless, there is a way to improve this behavior. In fact, the “extremum” 共instead of the zeros兲 of the quantity in the left hand side of Eq. 共15兲 with the Ansatz 关Eq. 共16兲兴 is exactly located at the critical point z = 0. In order to appreciate a shift from criticality, we have to include higher orders of = 共g − gc兲 in the irrelevant nonscaling term, i.e., L⌬共g,L兲 = ⌽0 + L2/2⌽2 + 共C共0兲 + C共0兲兲L−⑀ . In this case, the convergence is 兩gLⴱ − gc兩 ⬇ L−extr, where extr = 2 / + ⑀, which is better than the usual PRG not only because of the double exponent but also thanks to the coefficient in front: it is proportional to C共0兲, which is usually a small quantity. It is also worth noticing that in this case, extr = fast. III. TESTING THE METHODS A. XY model in transverse field
To analytically check our methods, we consider the 1D spin-1/2 XY model given by L
H=−兺 j=1
冋
册
共1 − 兲 y y 共1 + 兲 x x j j+1 + j j+1 + hzj . 共17兲 2 2
Throughout the paper, we will consider L even and PBC. This model can be solved exactly11,13 by means of a Jordan– Wigner transformation from spins to spinless fermions followed by a Bogoliubov transformation to arrive at a Hamiltonian of free quasiparticles. The number of original fermions N = 兺ic†i ci is not a conserved quantity, but its parity P ⬅ exp共iN兲 corresponds to a rotation around the z axis, and, therefore, is conserved. One should beware of a delicate issue concerning boundary conditions. Starting with PBC in Eq. 共17兲, the fermionic Hamiltonian turns to have PBC in the sector of odd parity P = −1. Instead, in the even parity sector P = 1—which comprises the ground state—antiperiodic boundary conditions must be used. In this sector, the model becomes
冉
H = 兺 ⌳共k兲 †k k − k
冊
1 , 2
共18兲
where k are Bogoliubov quasiparticles, k ranges in the first Brillouin zone, and the dispersion relation is
⌳共k兲 = 2冑2 + h2 + 共1 − 2兲cos2 k + 2h cos k.
共19兲
For ⫽ 0, this model displays an Ising transition at h = 1 with exponents = = 1, which means that the scaling variable is z = L共h − 1兲. To test our ideas, we need to calculate, for finite L, the GS energy eL共h , 兲 and the average potential bL共h , 兲 = −具zj 典 = eL / h in the quasicritical region given by z Ⰶ 1. According to Eq. 共18兲, the GS is given by eL共h , 兲 = −1 / 2L兺kn⌳共kn兲, where kn = 共2n + 1兲 / L , q = 0 , . . . , L − 1. We then expand the argument of the sum up to the desired order in z. For our purposes, we need eL共h , 兲 up to O共z2兲 and bL共h , 兲 up to O共z兲. The resulting sums are then evaluated with the aid of the Euler–Maclaurin formula.14 The results and the details are given in the Appendix. From Eqs. 共A2兲 and 共A3兲, one sees explicitly that the terms of order L−3 are absent both from eL and bL so that according to our definition Eq. 共5兲, we find ⑀ = 2. In passing, we notice that from Eq. 共A2兲, the Casimir-type term is −L−2兩兩 / 6 consistent with the CFT formula −L−2cv / 6. In fact, the central charge is c = 1 / 2 and, from the dispersion relation 关Eq. 共19兲兴, the spin velocity turns out to be v = 2兩兩. Now, by using Eqs. 共A2兲 and 共A3兲, we have all the elements to analytically derive the sequences of pseudocritical points. As far as the FSCM is concerned, the pseudocritical points are obtained by imposing bL共hLⴱ 兲 = bL+2共hLⴱ 兲 or, formally, LbL共hLⴱ 兲 = 0. Up to leading order in L, the solution is hLⴱ = 1 +
2 , 6L2
as already obtained in Ref. 8 for the Ising model 共兩兩 = 1兲. This explicitly shows that FSCM = 2, which is consistent with our prediction FSCM = 2 / 共note that = 1兲. For what concerns the “balancing trick” discussed in Sec. II, we can show that the solution of Eq. 共6兲 is given by ␥ⴱ = −1 / 22. The pseudocritical points are then given by imposing L⌫共hLⴱ , L , ␥ⴱ兲 = 0. In this case, at leading order, the solution is hLⴱ = 1 −
74共24 + 2 − 3兲 −4 L , 7202
which means fast = 4, which is once again consistent with our prediction fast = 2 / + ⑀ 共note that ⑀ = 2兲. From a numerical point of view, it is more profitable to use the HCM. From 3Lb + LL2 b = 0, we obtain hLⴱ = 1 +
74共22 − 3兲 −4 L 7202
with the same exponent fast = 4. Now, we proceed to discuss the PRG method for which knowledge of the lowest gap is required. The first excited state belongs to the sector with odd parity. Correspondingly, the finite-size gap is given, besides a constant term, by the difference between two Riemann sums where the sampling is taken over odd and even 共in units of / L兲 wave numbers. In analogy with Ref. 15, the lowest gap can be eventually written as ⌬L = 2共h − 1兲 + 关T共L兲 − 2T共L / 2兲兴, where T共L兲 = 共1 / 2兲兺2L−1 j=0 ⌳共j / L兲. Again, we refer to the Appendix for the details. By using the form of the gap 关Eq. 共A4兲兴, we can
155413-5
PHYSICAL REVIEW B 77, 155413 共2008兲
RONCAGLIA, CAMPOS VENUTI, AND DEGLI ESPOSTI BOSCHI
evaluate the various terms of the scaling function in Eq. 共16兲 as follows:
1
兩兩 3关3共h − 1兲 − 2共2 + h兲2 + 84兴 −2 L + ¯. +z+ 2 192兩兩3
0.9
⌽共z兲 =
共20兲
*
DL 0.8
Hence, ⌽0 = 兩兩 / 2 = 2vd, which is consistent with Ref. 15 for 兩兩 = 1 and with the known result d = 1 / 8 for the first excited state of the c = 1 / 2 minimal model.11 Moreover, since ⌽1 ⫽ 0, according to the general analysis in Sec. II D, there is no advantage in searching the minima or maxima in of the left hand side of Eq. 共15兲. Looking for the zeros by imposing the PRG equation L关L⌬L共hLⴱ 兲兴 = 0, yields hLⴱ = 1 +
3共42 − 3兲 −3 L + O共L−5兲. 48兩兩
This result for generic ⫽ 0 explicitly shows that the PRG shift exponent is PRG = 3 and proves a conjecture put forward by Hamer and Barber15 for the Ising model. From Eq. 共20兲 together with Eq. 共16兲, we read ⑀ = 2 as previously, so that the calculated exponent is consistent with the prediction PRG = 1 / + ⑀. The explicit calculation of the various shift exponents confirms that the HCM given by Eq. 共10兲 is superior to the standard PRG method. However, the things are different if we considered the excitation obtained by creating a welldefined quasiparticle with the smallest momentum k = 2 / L, instead of the very first excitation gap. In this case, the energy gap is given by Eq. 共19兲. By using the method of the extrema of the PRG quantity in the left hand side of Eq. 共15兲, we would have obtained a convergence to the critical point as L−4 while with the standard PRG only as L−2. B. c = 1 non-Berezenskii–Kosterlitz–Thouless transition
We choose the spin-1 Jz − D model on a chain 共d = 1兲 with PBC as follows: L
H = J 兺 Sជ j · Sជ j+1 + 共Jz − J兲Szj Szj+1 + D共Szj 兲2
共21兲
j=1
because the transition from the Haldane phase to the phase at large D is described by a c = 1 CFT 共 = 1兲 with continuously varying exponents. The Hamiltonian 关Eq. 共21兲兴 has been used to describe the magnetic properties of different quasi-1D compounds 共see Ref. 16 for a brief account兲. Here, by fixing J = 1 and Jz = 0.5, for which it has been already estimated = 2.38,17 we wish to test the methods described above in a severe case in which a 1 / scaling would give sublinear convergence in L. Indeed, for the FSCM, we expect FSCM = 0.84. The previous estimate by using the log-log plots of the finite-size gaps was Dc = 0.65,17 while in Ref. 18, it is found Dc = 0.635 by using the method of level crossing with antiperiodic boundary conditions, which is, however, specific to this transition. The first irrelevant operator allowed by the lattice symmetries has scaling dimension K + 2, where K = 2 − 1 so that ⑀ = min共K , ⑀lattice兲, where with ⑀lattice, we denote the smallest exponent of corrections arising
FSCM fast PRG PRG (minima) HCM
0.7
0.6 0
0.05
1/L
0.1
0.15
FIG. 1. 共Color online兲 The sequences of pseudocritical points obtained with the FSCM, the fast-convergent method of Eq. 共7兲, and the HCM described by Eq. 共14兲. Moreover, the sequences have been calculated with PRG methods, both standard and in our improved version. Continuous lines are algebraic best fits to the data.
from lattice contributions at finite L that are not captured in the framework of the 共relativistic兲 continuum theory. As discussed in Refs. 16 and 19, this is a case in which the linear term ⌽1 in the scaling function of the PRG vanishes because b = 0 for the sine-Gordon model, which is the effective field theory that describes the surroundings of the c = 1 line. So, it is convenient to also use our improved version of the PRG, as discussed in Sec. II D. The expected shift exponent is extr = 2 / + ⑀. In order to have an idea of the range of values of L to be used, let us imagine that ⑀ = K so that fast = extr = 2 + 1 / = 2.42. With L ⬎ 10, this exponent leads to variations in the pseudocritical points DLⴱ smaller than O共10−3兲. Hence, we prefer to illustrate the method with virtually exact numerical data obtained with the Lanczos algorithm by using L = 8 , 10, 12, 14, 16. The reason is that the DMRG would give rather accurate values for the energies but the estimates for b = 具共Sz兲2典 could not be sufficiently precise to appreciate the variations in DL obtained from the crossings. So, we use the DMRG only to extend the data to L = 18, 20 with 37 opz =0 timized states. In any case, the GS belongs to the Stot sector. = 0.647 with The FSCM with b = 具共Szj 兲2典 yields DFSCM c FSCM = 0.79 in reasonable agreement with the CFT expectations. The high-precision procedure based on Eq. 共7兲 yields Dfast c = 0.633 with fast = 7.6. The value of the shift exponent is definitely larger than what is expected, which could be due to the vanishing of the coefficient of the first irrelevant contribution with scaling dimensions K + 2 ⬍ ⑀lattice. In any case, from Fig. 1, one can clearly appreciate that the sequence DLFSCM converges more slowly than DLfast. As far as the PRG is concerned, with the standard procedure of finding the ze= 0.640 with PRG = 1.50. Inros of Eq. 共15兲, we get DPRG c stead, by using the improved method estimation by looking for the extremal value of the quantity on the left hand side of Eq. 共15兲, we find a rapidly converging sequence that, however, also oscillates between 0.636 and 0.635. Finally, a small oscillation in the fourth decimal place of about 0.6305 is also seen in the sequence DLhom obtained through the
155413-6
PHYSICAL REVIEW B 77, 155413 共2008兲
RAPIDLY CONVERGING METHODS FOR THE LOCATION…
homogeneity criterion 关Eq. 共14兲兴. As expected, the latter sequence converges in a fashion similar to DLfast. All of these results are summarized in Fig. 1. The data show that subtracting out the terms that induce a slow convergence of the pseudocritical points, one remains with quantities that inherit a residual 共possibly oscillating兲 L dependence from the specific lattice model and that could be very difficult to account for. Other factors that affect the extrapolation of the critical points at this level of accuracy are the sampling ␦D and the trade-off between computational accuracy and the maximum available size. To summarize, we observe that our improved methods yield very fast convergence to the critical point that we estimate to be Dc = 0.633⫾ 0.02, which is consistent with Ref. 18.
12 16 20 24 28
0.02
3
2
L b’’+(3L -16)b’
0.04
0 -0.02 -0.04 0.2
0.25
0.4
0.3
0.35
J 2 /J 1
0.4
0.45
0.5
0.35
C. c = 1 Berezenskii–Kosterlitz–Thouless transition
Now, we consider the spin-1/2 Heisenberg model with frustration due to next-to-nearest neighbor interaction as follows:
J2,L
L
*
0.3
0.25
H = 兺 J1S j · Sជ j+1 + J2S j · Sជ j+2 . j=1
The model is equivalent to a two-leg zigzag ladder with L / 2 rungs. The best estimate of the critical point was made by Okamoto and Nomura9 by using a model-specific crossing method; with exact diagonalizations up to L = 24, they determined J2c / J1 = 0.2411⫾ 0.0001. The model is gapless for J2 ⬍ J2c and has a doubly degenerate GS in the TL for J2 z ⬎ J2c. Again, the GS has Stot = 0. Without exploiting a priori information about the BKT character of transition 共if not the value of the dynamic exponent = 1兲, we tested the homogeneity criterion 关Eq. 共14兲兴 by using ladders with PBC of up to L / 2 = 16 rungs, an effective ␦L = 4, and 1024 DMRG states that ensure an accuracy of O共10−7兲 on the values of b = 具Sជ j · Sជ j+2典. The results reported in Fig. 2 are encouraging: while with the FSCM we would get no crossings at all, the zeros of Eq. 共14兲 yield a sequence of points converging to J2 ⬃ 0.25. The main problem comes from the left side of the transition where the truncation DMRG error for L = 28 and 32 induces some oscillations on the plotted quantity. We content ourselves with linear fits in 1 / L. If we exclude the point with L = 12, the fit is better even if we find J2c / J1 = 0.2553⫾ 0.0008; by selecting all the available points, instead, the fit is visibly worse but the extrapolated value is J2c / J1 = 0.242⫾ 0.006, which is in agreement with Ref. 9. As above, apart from the details of the extrapolation procedure, we see that the homogeneity criterion provides a viable procedure to locate the critical point in a BKT transition, where almost all existing generic methods fail. We remark that this analysis is solely based on b共g , L兲, namely, an observable evaluated on the GS, without invoking further assumptions on the nature of the excitations. IV. CONCLUSIONS
By making use only of finite-size quantities related to the ground state, we show how to generate sequences of pseud-
0.2 0
0.02
0.04
1/L
0.06
0.08
FIG. 2. 共Color online兲 Upper panel: Searching for the zeros of L3b⬙ + 共3L2 − 16兲b⬘ according to Eq. 共14兲 with J1 = 1 and ␦L = 4. Lower panel: Extrapolations of the pseudocritical points J2,L 共cyan line for all the points and blue line with L = 12 excluded兲; the black cross marks the accepted value 0.2411 共see text for further details兲.
ocritical points that converge very fast to the infinite-size critical point. The convergence is of the form L− with a shift exponent . In this paper, we propose a homogeneity condition method 共HCM兲, which is faster than the standard PRG in locating the critical points. Moreover, its validity is more general as it can be applied without modification to the difficult case of a Berezinskii–Kosterlitz–Thouless transition. The homogeneity method requires only the knowledge of b共g , L兲 = 具V典 / Ld, that is, the expectation value of the term that drives the transition. We also presented an improvement to the PRG method, which allows us, under certain conditions, to obtain pseudocritical sequences characterized by the same shift exponent as that for the HCM. However, this modification, relying on a particular form of the gap scaling function, is not valid in general. It holds true, for instance, for the sine-Gordon model that underlies a variety of transitions in 共1 + 1兲 dimension. The formulations of the approaches are sufficiently general to be applied in any spatial dimensionality. Even if we are primarily interested in quantum phase transitions ideally at zero temperature, in principle, the methods can be extended to problems of finite-temperature statistical mechanics. At variance with other accelerating methods found in the literature, e.g., the van den Broeck–Schwartz or the Bulirsch–Stoer ones 共reviewed in Chap. 9 of Ref. 11兲, the procedures presented here rely on the scaling behavior of
155413-7
PHYSICAL REVIEW B 77, 155413 共2008兲
RONCAGLIA, CAMPOS VENUTI, AND DEGLI ESPOSTI BOSCHI
thermodynamical quantities expected from physical and renormalization arguments. The validity of the methods has been tested with analytical calculations on the one dimensional XY model in transverse field and numerically on a nontrivial spin-1 chain with anisotropy. As an extreme case, we have shown that the homogeneity condition method provides a satisfactory location of the critical point also in the case of Berezinskii– Kosterlitz–Thouless transitions. These confirmations motivate us to consider systems in higher spatial dimensionality. In these cases, the numerical data are restricted to smaller system-size and the need for fast-converging pseudo-critical sequences is a prerequisite for the precise location of the critical points.
in 共a , b兲. Here, Bn are the Bernoulli numbers and the remainder Rm depends on f 共2m兲 on 共a , b兲. Some care must be taken when the function f diverges at the border of the Brillouin zone; in this case, one must keep a and b away from the borders. Moreover, sending m to infinity in Eq. 共A1兲, some sums must be regularized by using a Borel summation technique. The final result for the ground-state energy is eL共h, 兲 = e⬁共1, 兲 −
冋
+ 共h − 1兲 b⬁共1, 兲 − −
ACKNOWLEDGMENTS
This work was partially supported by the Italian MIUR through the PRIN Grant No. 2005021773. M.R. acknowledges support from the EU 共SCALA兲. Numerical calculations were performed on a cluster of machines made available by the Theoretical Group of the Bologna Section of the INFN.
L−1
kn =
where ␥C = 0.577 216. . . is the Euler–Mascheroni constant and the thermodynamic values are given by e⬁共1, 兲 = −
+
冉冊 z L
N
bL共h, 兲 = b⬁共1, 兲 −
j=0
冕
a
册
+兺 r=0
共2r兲!
ln共L兲 + ␥C + ln共8兩兩/兲 − 1 + O共z2兲. 兩兩
L−1
␦
␦2r
.
−2 73共3 − 2␥2兲 −4 L − L 12兩兩 2880兩兩3
T共L兲 =
f共x兲dx + 关f共a兲 + f共b兲兴 2 m
− 共h − 1兲
1 2 兩 ⌳兩h=1共kn兲 + O共z3兲 · 2 h
b
冑1 − 2
Finally, the sum
The resulting sums can be computed by using the Euler– Maclaurin formula 共see, e.g., Ref. 14兲
␦ 兺 f共a + j␦兲 =
2 arctan共冑1 − 2/兩兩兲
共A3兲
1 z ⌳h=1共kn兲 + 兩h⌳兩h=1共kn兲 eL共h, 兲 = − 兺 2L n=0 L 2
册
2兩兩 2 arctan共冑1 − 2/兩兩兲 + , 冑1 − 2
By using similar procedures, we obtain the following expression for average potential bL共h , 兲:
2n + 1 , L
冋
冋
b⬁共1, 兲 = −
where ⌳共k兲 is given by Eq. 共19兲. We need to investigate the above sum in the quasicritical region z = L / Ⰶ 1 so that it is sufficient to expand eL共h , 兲 in powers of 共h − 1兲 = z / L as follows: L−1
册
共A2兲
Here, we indicate how to compute the mean energy per site eL and the average potential bL, as required in Sec. III A. Consider, for example, the energy sum 1 兺 ⌳共kn兲, 2L n=0
−2 L 12兩兩
共h − 1兲2 ln共L兲 + ␥C + ln共8兩兩/兲 − 1 + O共z3兲, 2 兩兩
APPENDIX
eL共h, 兲 = −
兩兩 −2 73 共42 − 3兲 −4 L − L 6 360 4兩兩
冉 冊
2 j 1 ⌳ 兺 L j=0 L
for the evaluation of the finite-size gap in the PRG method, can be treated along similar lines. The final result for the gap is ⌬L = 共h − 1兲 +
B2r关f 共2r−1兲共b兲 − f 共2r−1兲共a兲兴 + Rm 共A1兲
valid for a function f with at least 2m continuous derivatives
155413-8
+
共22 + h − 1兲 −1 L 4兩兩
3关3共h − 1兲 − 2共2 + h兲2 + 84兴 −3 L + O共L−5兲. 192兩兩3 共A4兲
PHYSICAL REVIEW B 77, 155413 共2008兲
RAPIDLY CONVERGING METHODS FOR THE LOCATION… Sachdev, Quantum Phase Transitions 共Cambridge University, Cambridge, 1999兲. 2 U. Schollwöck, Rev. Mod. Phys. 77, 259 共2005兲. 3 H. G. Evertz, Adv. Phys. 52, 1 共2003兲. 4 A. W. Sandvik, arXiv:0710.3362 共unpublished兲. 5 N. Schuch, M. M. Wolf, F. Verstraete, and J. I. Cirac, Phys. Rev. Lett. 100, 040501 共2008兲. 6 M. E. Fisher and M. N. Barber, Phys. Rev. Lett. 28, 1516 共1972兲. 7 M. N. Barber, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. L. Lebowitz 共Academic, New York, 1983兲, Vol. 8. 8 L. Campos Venuti, C. Degli Esposti Boschi, M. Roncaglia, and A. Scaramucci, Phys. Rev. A 73, 010303共R兲 共2006兲. 9 K. Okamoto and K. Nomura, Phys. Lett. A 169, 433 共1992兲. 10 J. G. Brankov, D. M. Danchev, and N. S. Tonchev, Theory of Critical Phenomena in Finite-Size Systems 共World Scientific, 1 S.
Singapore, 2000兲. M. Henkel, Conformal Invariance and Critical Phenomena 共Springer, New York, 1999兲. 12 J. L. Cardy, J. Phys. A 19, L1093 共1986兲. 13 E. Lieb, T. Schultz, and D. Mattis, Ann. Phys. 共N.Y.兲 16, 407 共1961兲. 14 K. E. Atkinson, An Introduction to Numerical Analysis, 2nd ed. 共Wiley, New York, 1989兲. 15 C. J. Hamer and M. N. Barber, J. Phys. A 14, 241 共1981兲. 16 C. Degli Esposti Boschi and F. Ortolani, Eur. Phys. J. B 41, 503 共2004兲. 17 C. Degli Esposti Boschi, E. Ercolessi, F. Ortolani, and M. Roncaglia, Eur. Phys. J. B 35, 465 共2003兲. 18 W. Chen, K. Hida, and B. C. Sanctuary, Phys. Rev. B 67, 104401 共2003兲. 19 A. Kitazawa, J. Phys. A 30, L285 共1997兲. 11
155413-9
