Scaling corrections to the ground state energy of the spin-½ isotropic anti-ferromagnetic Heisenberg chain

Solutions to the Bethe Ansatz equations for the ground state of spin-½ isotropic anti-ferromagnetic periodic Heisenberg chains to length L = 221 are obtained by combining Lagrange interpolation with Newton–Raphson iteration. The long chain lengths allow many powers of a renormalization group running coupling constant to be included in fits to the ground state energy and make possible the confirmation of the convergence of the leading logarithmic term. The amplitude of this term is consistent with that expected on the basis of conformal field theory and the connection of the discrete spin-½ system to the continuum Wess–Zumino–Witten model. This resolves a decades old discrepancy based on analysis of shorter chains. An analytical improvement to the Hulthén wave-vector distribution is also provided.


Introduction
This paper presents the derivation and asymptotic analysis of the ground state energy of large (even) length L periodic chains of s=½ spins anti-ferromagnetically coupled as defined by the Hamiltonian where i labels both the sites and distance along the chain and the components of s   i are the Pauli spin matrices. Eigenstates of (1) can be labelled by total S and S z ; for the ground state S=S z =0. The motivation for these calculations was an attempt to resolve a long standing apparent failure of universality first noted by Affleck et al [1]. More evidence for failure was provided by Nomura [2] in an analysis of spin-½ chains to length L=16 384.
The history of this model problem dates to 1931 when Bethe [3] (an English translation appears in [4]) conjectured that all eigenstates of (1) could be found based on the solutions of nonlinear algebraic equations that he described. For the ground state these are equations for L/2 distinct real wave-vectors, each being associated with one of the s z = −½ spins on the chain. Subsequently Hulthén [5], starting from the Bethe Ansatz, obtained the ground state energy E L in the limit L→∞. Hulthén's result is The development of conformal field theory (CFT) in the 1980s led to the realization that E L would approach the Hulthén value Lε ∞ as 1/L with a coefficient constrained by the 'conformal anomaly' of the field theory in the universality class of the spin-½ anti-ferromagnet. Avdeev and Dörfel [6], building on earlier work by de Vega and Woynarovich [7], concluded and confirmed this value by numerical evaluation of E L for chains to length L=256. Hamer [8] came to the same conclusion independently.
Affleck et al [1] argue that the spin-½ anti-ferromagnet is in the universality class of the k=1 Wess-Zumino-Witten nonlinear σ model and use this to identify a marginal operator that will lead to corrections to E L beyond the conformal anomaly term. Arguments by Cardy [9] specify what form these corrections must take and Affleck et al conclude 1 ln In deriving (6), Nomura incorporated Cardy's conclusion that asymptotically the C correction term is proportional to g 3 with g a running coupling that satisfies the renormalization group equation with the ellipsis indicating unknown terms. Nomura attributed the difference between (5) and (6) to an invalid reliance by Woynarovich and Eckle on the Euler-Maclaurin sum formula resulting in uncontrolled errors. Of course, even if we accept (5) as unreliable, the apparent failure of universality remains because of the (4) and (6) difference that lies well outside Nomura's error assignment. But this assignment can be faulted because it does not include sensitivity to additional terms in β(g) or variation in the constant of integration of (7). The major goal of the present paper is to provide a realistic appraisal of such effects by first generating E L data for much longer chains. The outline of the remainder of the paper is as follows. Section 2 expands on the above historical summary by providing formulas useful for subsequent analysis. Particular attention is paid to the approximation necessary to derive (3) from the Hulthén [5] wave-vector distribution. Section 3 is a derivation of an analytic improvement on the Hulthén wave-vector distribution. The approximation that led to (3) can be applied to this improved distribution and leads to a result of the form of (4) but with C=(3/8)ln 2 (2)-a roughly factor 2 underestimate from the CFT prediction. Section 4 describes how Lagrange interpolation is combined with Newton-Raphson (NR) iteration for an efficient generator for E L for large L. Results to a maximum L=2 21 are given. Section 5 describes the analysis of the data in section 4 with the conclusion that the CFT result (4) is unambiguously confirmed.

E L alternatives from the Hulthén wave-vector distribution
The ground state is specified by a distinct real wave-vector k n , 0<k n <2π, for each of the L/2 overturned spins from a ferromagnetically aligned state. A convenient change of variable is to the 'rapidity' l = ( ) k cot 2 ; n n these variables for the ground state satisfy the Bethe ansatz equations (BAE) where å m and similar sums hereafter are understood to range over all L/2 listed half-integer values. The arccot function is its principal value; i.e. 0arccot(x)π for −∞ x ∞. The symmetry l l = --L n n 2 allows reduction to L/4 ((L+2)/4) independent variables if L/2 is even (odd). The BAE (8) with its discrete solutions l n can also be interpreted as defining a continuous function m m with n taking on half-integer values whenever l is one of the solutions l n of (8). This point of view is useful in a number of instances and, for example, allows one to identify n(∞)=0 and n(−∞)=L/2 as the limits on the range of possible n. The energy of the ground state is The NR solution of (11), dl = + - 1 1 converges quadratically but because it requires numerical matrix inversion it is practical in the elementary form given here only to lengths of several thousand. By reparameterizing l in some appropriate basis, solutions to lengths L≈2×10 6 have been obtained as described in section 4.
A useful initialization is provided by the Hulthén solution which is based on the approximation   L term. Indeed, the expected energy based on CFT is that given in (4). A formula equivalent to (10) for E L was derived by de Vega and Woynarovich [7]. The key identity, restricted here to the isotropic chain (1) and in the present notation, is where l n in the integral in (19) is understood to be the inverse of l ( ) n from (9). The proof of (19) follows trivially on making the change of integration variable from n to l using l n d d calculated from (9); explicitly, the last term in (19) is ò ò p pl l pl l l l l p pl where the integrals evaluated for the final equality are those that already appeared in (16) and (17). Substituting (20) into (19) yields (10) and completes the proof. Note that (20) relies only on the functional form of (9) and not that the l m in the sum in (9) satisfy the BAE. Thus any approximation to l n in (19) (and simultaneously to the l m in (9)) has the same effect on E L as does the approximation when applied directly to (10). In particular (19) reproduces the O(1) error found for the Hulthén approximation ¢ E L in (18). On the other hand, different E L estimates using (19) are possible if one drops the constraint that the analytic continuation from a discrete l n list to a continuous l n function be via (9). If the approximate l n is a known analytic function such as the Hulthén l ( ) n 0 in (14), it would seem more natural to define (19) by the condition that the same analytic l n appear in both sum and integral. As an example, if in (19), with the constraining (9) removed, we replace l n by the Hulthén l ( )  a result first obtained by Avdeev and Dörfel [6] and Hamer [8]. The leading correction in (21) agrees with the CFT value in (4) which illustrates the dramatic improvement to an o(1/L) error that has been achieved. This is a surprise since the removal of the constraint (9) has eliminated the justification for (19) to be the formula for E L . It is of course possible that the dramatic error reduction has just been accidental-which motivates further exploration in the following by observing what changes are induced in (21) when the modified (19) is applied to an analytical improvement on the Hulthén solution (14).

An improved wave-vector distribution
Before proceeding with the calculation of such an improvement, it is worth introducing a change of variable l n  n n in which we write the exact l n as The unit spacing between n sets the scale for distinguishing n n as either small or large. In the case that n n is small we find from (22) that  The numerical evidence discussed below suggests n | | n is bounded by 0.03 for all n and L and so in the following we use only the linear in n n versions of (23) and (24).
Any iteration of (11) to improve on the Hulthén l ( ) n 0 ideally starts with an exact evaluation of the sum contributing to ( ) h n 0 in (12). Fortunately this can be done analytically by complex contour and residue methods similar to that described in the appendix. The details involve no new concepts but it is useful to first recast the arccot function in the sum in (12) as an integral as was done in the first equality in (16). One can confirm that the singularities in this case at l m =   where (With the substitution y=ln(1+e x )=arcsinh(e x (1+e x /2)/(1+e x )), −∞<x<∞, integrals such as (25) can be done numerically as an equally weighted sum on a uniform grid with exponential convergence both with respect to the large | | x cutoffs and the (inverse) grid spacing.) To make the right hand side of (11) tractable we employ the Hulthén integral approximation (15). The necessary integral for the sum inˆ( the last equality being the linearized form of (23). The (integral approximated) contribution from dl   As an example of the accuracy of (33), comparisons of n ( ) n 1 with the exact n n obtained by NR iteration are shown in figure 1 for L=2 10 =1024 and L=2 21 =2097 152. Extrapolation from discrete to continuous n for 3 The resulting (29) has the structure of the Fourier transform equation that leads to the Hulthén solution (14). The unknown n ( ) the NR solution is by the use of (9). A single NR iteration starting from the Hulthén solution is graphically almost indistinguishable from the final (multiple iteration) result shown in figure 1 so that the error in n ( ) n 1 is inferred to be dominated by the Hulthén integral approximation made in its derivation. The qualitative agreement between n ( ) n 1 and the exact n n for the lengths L shown in figure 1 suggests we might use (33) to guess how n n will approach its L→∞ limit.
We consider two cases. In the limit L→∞ with n=O(1), we simplify (33) by first noting The Γ functions involving -( ) g y n can be expanded in an asymptotic series in 1/ln(L) and the remaining integrals  The first (L independent integral) term in (35) is 0.019 181, 0.012 768, 0.01 0168, K for = ¼ n , , while the remaining term suggests ν n (L)−ν n (∞) will be O(1/ln 2 (L)) for n=O (1). A second asymptotic result is the slope of n ( ) n 1 at the symmetry point n=L/4. We have ò n l p y y y . The equality in (37) is the result of an explicit calculation of the sum minus integral difference by the complex contour method described in the appendix. For the leading terms in (37) in the limit L→∞ we can set ¢ = ( ) y L cosh 2 1 and all ln(Γ(K)) terms containing the ratio ¢ t t then vanish because of the anti-symmetry of their imaginary parts under ¢ « y yinterchange. Further reduction analogous to that in (34) followed by asymptotic expansion of the Γ functions gives the energy correction ò ò

E L by NR with Lagrange interpolation
The NR procedure for finding the solutions l n of the BAE (8) as described in the discussion of equations (11)-(13) is impractical for chains of length L greater than several thousand. To deal with long chains we restrict the number of variables treated by NR to a limited number M of l n m >0, m=1, 2, K, M, and use a Lagrange interpolation scheme that accurately and efficiently determines all remaining l n in terms of the basis variables l . n m In our implementation of the Lagrange interpolation we assume l n is a smooth function l ( )   It remains to discuss the choice of the special n=n m . It is essential that there be a minimum number of n=n m = ¼ , , 1 with unit gaps since no accurate interpolation is possible near the end-points of the l interval. Following this are gaps of increasing size between n until a final n=n M near L/4. Because we are interpolating in x it is plausible that the spacing in the corresponding l = ( ) x n n 0 m m should be nearly uniform. A possible formula to achieve this is a modification of the Hulthén l ( ) ( ) n n 0 from (14) which we define as a r c t a ne x p 2 1 , 4 1 m where the last term has been added to provide extra flexibility and ¢ K K , are parameters dependent on L. If we define a crossover integer m x as the largest (integer) m for which which is an approximation to the defining equation for the crossover m x and ¢ = - for mm x in (42). These two equations can be reduced to where now the first equality determines K after which ¢ K follows trivially. For example, at L=2 21 with M=92 and m x =60, K =32.7464 and ¢ K =2.3858. This is a special case of the parameter choice which we have used and have found that for 2 10 L2 21 yields an absolute error in the scaled energy deviation that is less than 10 -41 . The observed trend suggests this bound will improve slightly for L>2 21 . The oscillatory behaviour of dl n is typical with the decay in amplitude towards large x the (deliberate) result of the large M-m x choice (44). The decay in amplitude towards small x is not easy to modify given the constraint of the functional form (41) but it is appropriate as it partially compensates for the growth in the number of l n terms between adjacent l n m at small x and leads to a substantial suppression of the dE n oscillation amplitude there. The further suppression of these oscillations by about 2 orders of magnitude to a final dE 1 2 (the rightmost point in figure 2) is also typical.
This explicit cancellation in the generation of d / E 1 2 implies a loss of significant digits but a more important loss of significance is the cancellation between terms in the sum for l n in (39). In total up to 15 digits can be lost which is in addition to the 12 digits lost in cancellation in the scaled energy deviation 1+E (corr) in (45). Significant cancellation between terms in (40) can also occur so in view of this, for all L>10 6 , a conservative extended precision arithmetic of 90 digits was used; somewhat less for shorter chains.   1 m given by (33) (see the note following (26) for an efficient integration scheme). Whatever initialization is chosen, this must be carried to a precision comparable to that used for a general iteration term-otherwise roundoff noise can lead to failure to converge.
The energies from the NR/Lagrange calculation for 2 10 L2 21 given below are exact to the 40 digits listed. They confirm the Nomura values for L2 14 . The last three columns are the fit deviations δ=E (corr) | fit −E (corr) of formulas described in the next section and referenced by their equation numbers. 5 Error determination is based on comparison with a more accurate calculation using expansion in a 204 Chebyshev polynomial basis set. The Lagrange interpolation scheme described here is easier to implement, much less susceptible to round-off error and instability and, for the same basis size, essentially equivalent in accuracy. The analysis of the energies (47) described in the following section unambiguously confirms the CFT result (4).

Energy scaling analysis
The logarithmic corrections to scaling beyond the conformal anomaly term are contained in E (corr) as defined in (45) and displayed in (47). A fit which includes analytic corrections, incorporates the expected CFT asymptotic scaling but treats C as a fitting parameter rather than fixed at the CFT value 3/8. The coupling g in (48) satisfies the renormalization group equation which allows for additional terms not present explicitly in (7). Integration of (49) yields ò a a = - where L 0 is a (model dependent) integration constant. When α(g)=0, (48) combined with (50) with L 0 =0.5653 is the Nomura [2] fit 6 . For long chains the analytic correction b/L 2 in (48) is unimportant and E (corr) /g(L) 3 , taken as a proxy for C, is shown as squares in figure 3 using the Nomura data for 2560L16 384. This is supplemented as crosses using the additional data (47). The trend with increasing L clearly shows Nomura's conclusion C≈0.365 is no longer tenable. As the first and most important modification of Nomura, allow L 0 as a free parameter keeping α(g)=0. Determine a sequence of C, b and L 0 from 3-point fits of E (corr) | fit_1 (48) to E (corr) data (47) at L=2 p , 2 p-1 and 2 p-2 , 12p21. The results for C for integer p are shown as the lowest diamond sequence in figure 3. The overlapping curve is a (negative) deviation from 3/8 that is proportional to y 2/3 =1/ln 2 (L) with an amplitude such that it passes through the point for L=2 21 . This agreement between the C deviation from 3/8 and a pure 1/ln 2 (L) power suggests that we introduce a constant α(g)=c 0 into β(g) in (49). The resulting added integral term≈c 0 g in (50) is a relative O(g 2 ) correction to the leading 1/g which in turn is an implicit relative O(1/ln 2 (L)) correction to g 3 in the 3-parameter E (corr) | fit_1 in (48). To check this determine a new sequence of C, b, c 0 and L 0 from 4-point fits to E (corr) (47) at L=2 p , K 2 p-3 , 13p 21. The results for this new C sequence are shown as the highest cross sequence in figure 3 together with a (positive) deviation curve proportional to y 4/3 =1/ln 4 (L). The new agreement between the C deviation and a pure power suggests the process we have started be continued, first with an α(g)=c 0 +c 2 g 2 and then with α(g)=c 0 +c 2 g 2 +c 4 g 4 , resulting in the two remaining cross sequences shown in figure 3. All these fits are plausibly consistent with a deviation pattern of alternating sign and of magnitude proportional to y 2 n/3 , n=1, 2,K . Although the fits become unstable for n much beyond four, the sequences already shown in figure 3 are impressive evidence that the CFT result (4) is correct and applies to the spin-½ isotropic Heisenberg chain.
As an alternative analysis consider keeping α(g)=0 throughout but changing E (corr) | fit_1 to include explicit corrections to g 3 of the form g 5+2n which are the leading corrections, deduced by power counting, that would be induced by the powers g 2n in α(g). Define   where the first column entries labelled by L are N=6 diamond estimates 10 5 (C -3/8). The nth following column entries are the corresponding scaled C ∞ deviations where I take 1/ln(L) as a reasonable proxy for g(L) and use the fitting function C L =C ∞ +c 1 /ln 8 (L)+K+c n /ln 6+2n (L). The Neville table using the N=6 cross estimates is similar. There is nothing in the extrapolated C ∞ values to suggest any systematic deviation from the CFT C=3/8 but they are disappointing-being at best an order of magnitude improvement on the N=6 fit in figure 3. To see whether more substantial improvement is possible I turn to higher order fits that utilize more of the large L raw data from (47).  The coefficients in (53) have been too severely truncated for (53) to be used to determine the fit accuracy but the deviations δ=E (corr) | fit -E (corr) calculated directly from the original multiple precision fit are recorded as column δ(53) in (47). Each zero in that column marks a datum that has been used in the fit to generate (53). The deviation of the C coefficient in (53) from 3/8 is fairly typical of such high order fits and illustrates that, even with data to L=2 21 , obtaining C to an accuracy better than 1 part in 10 5 does not seem possible. An intriguing observation is the ≈1.19 amplitude of the 1/L 2 term in (53); this is reasonably close to the 7π 2 /60=1.151K Avdeev and Dörfel [6] and Hamer [8] result given in (21) which if exact would be an even greater enigma than suggested in remarks following (21). Another observation regarding (53) is the rapid growth of the coefficients of g n with increasing n. While this might indicate the exact g series is asymptotic rather than convergent, a more likely reason is that no term in the g series has the exactly correct L dependence. This would imply the higher order terms must in part fit the lower order term errors which necessarily requires a growth in coefficients to compensate the decay in g n magnitude with increasing n. One can speculate that this error is intrinsic to the single variable approach-a proper renormalization group treatment probably requires a Wilson structure of coupled equations in many variables. A more accessible numerical question is whether the CFT amplitude C=3/8 can be imposed without a dramatic reduction in the quality of a fit such as (53). The answer is in the affirmative; an example is in which a new term has been introduced into β(g) in lieu of a variable C. The fit deviations of (54) are shown as δ  The deviations δ(55), shown as the last column in (47)  The deviations δ(56) are about 3% smaller than the δ(55)-an inconsequential difference. The effect of other changes such as including a 1/L 4 term or using only data from the largest chain lengths in the fits are also small. Provided the CFT amplitude C=3/8 is accepted, I estimate from (54)-(56) and other fits that = ( ) ( ) L 0.48364 2 . 57 0 It is interesting to note that with this value for L 0 rather than the Nomura L 0 =0.5653, the ratio E (corr) /g(L) 3 approaches 3/8 from above rather than below. While it is clear that no significance can be attached to any single high order term in the fits, the global quality of the fits as seen in the deviation columns in (47) is good and I estimate that for L>2 21 , the error in the fits (55) or (56) is bounded by ≈10 -13 with the maximum error at L≈10 21 .

Conclusions
Numerical solutions to the BAE for the ground state of spin-½ isotropic anti-ferromagnetic periodic Heisenberg chains to length L=2 21 have been obtained. An analysis of the ground state energy shows unambiguously that the amplitude of the leading logarithmic term is consistent with that expected on the basis of CFT and the