{'fraction_non_alphanumeric': 0.04134610954991477, 'fraction_numerical': 0.017160363016538442, 'mean_word_length': 4.351947731755424, 'pattern_counts': {'":': 0, '<': 16, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perturbative renormalization group, exact results and high temperature series to order 21 for the N-vector spin models on the square lattice

High temperature expansions for the susceptibility and the second correlation moment of the classical N-vector model (also known as the O(N) symmetric Heisenberg classical spin model or the as the lattice O(N) nonlinear sigma model) on the square lattice are extended from order beta^{14} to beta^{21} for arbitrary N. For the second field derivative of the susceptibility the series expansion is extended from order beta^{14} to beta^{17}. For -2<N<2, a numerical analysis of the series is performed in order to compare the critical exponents gamma(N), nu(N) and Delta(N) to exact (though nonrigorous) formulas and to compute the"dimensionless four point coupling constant"g_r(N). For N>2, we present a study of the analiticity properties of chi, xi etc. in the complex beta-plane and describe a method to estimate the parameters which characterize their low-temperature behaviors. We compare our series estimates to the predictions of the perturbative renormalization group theory, to exact (but nonrigorous or conjectured) formulas and to the results of the 1/N expansion, always finding a good agreement.


I. INTRODUCTION
We have extended the high temperature (HT) series expansion of the zero field susceptibility χ(N; β) and of the second correlation moment µ 2 (N; β) to order β 21 and the expansion of the second field derivative of the susceptibility χ 4 (N; β) to order β 17 for the Nvector model [1] (also known as the O(N) symmetric Heisenberg classical spin model) on all bipartite lattices in d = 2, 3, 4, 5, .. space dimensions [2][3][4]. The series coefficients have been determined by using the vertex renormalized linked cluster expansion (LCE) method [5][6][7] and have been expressed as explicit functions of the spin dimensionality N. This calculation pursues and improves our previous work [8,9] to a considerable extent: it summarizes into a convenient format a large body of information for an infinite set of universality classes and offers further insights into the properties of the N-vector model by enabling us to vary with continuity the crucial parameter N and to study how various physical quantities depend on N.
Strictly speaking, the N−vector model is defined only for positive integer N. Therefore it is possible to construct infinitely many "analytic interpolations" in the variable N of the HT coefficients and, as a consequence, of the physical quantities. We have performed the "natural" analytic interpolation by which the HT coefficients are rational functions of N. This is the most interesting interpolation because it coincides with that used in the 1/N expansion as well as in the usual renormalization group (RG) treatments and moreover it is unique in the sense of Carlson theorem [10].
Next interesting step, on which we will report elsewhere [4], is to compile tables of HT coefficients analytically interpolated both in N and in the space dimension d. The "natural" analytic interpolation, with respect to d, of the HT coefficients, which is polynomial in d and equivalent [11] to the one of the Fisher-Wilson method [12], is also unique in the sense above specified. We will thus be able to describe accurately the general (N, d) universality class, to achieve also a "view from HT" of presently inevitable RG approximation schemes such as the ǫ expansions (at the upper and at the lower critical dimension), the 1/N and the 1/d expansions and possibly to gain a more detailed knowledge of their limitations. As pointed out in Ref. [11,13] possible violations of the convexity of the free energy or of the Lee-Yang property might occur for noninteger values of N or for noninteger values of d, respectively, but, of course, they do not question the long known conceptual and analytical advantages of treating N or d as continuous parameters.
This paper is devoted to the study of the N-vector model on the square lattice and its main result reported in the appendices, is the tabulation of the HT series expansion coefficients. We have made an effort to keep our exposition selfcontained, in particular by reporting all known HT expansion coefficients of χ, µ 2 and χ 4 and not only the newly computed ones. It is worth noting that further sizable extensions of these series are not too difficult and now are in progress. Several interesting, but somewhat intricate, computational procedures that have made this laborious calculation (and its forthcoming extensions) possible will be illustrated elsewhere [4].
We shall also update, refine and extend our previous numerical analysis of O(β 14 ) series presented in Ref. [9], but our discussion will be more sketchy whenever there is some overlap with that reference.
The paper is organized as follows. In section II we introduce our notations and definitions.
In section III and IV we present a numerical study of our new HT series which are now long enough for a reliable assessment of the uncertainties of the analysis initiated in [8,9]. In section III we examine the range −2 ≤ N < 2 in which the N-vector model is known (or expected) to have an ordinary power law critical point at finite temperature and we estimate the critical exponents γ(N) of the susceptibility, ν(N) of the correlation length, as well as the gap exponent ∆(N) from χ 4 in order to compare them with exact (nonrigorous) formulas proposed some time ago [14][15][16]. We do not reanalyse here the N = 2 case (the Kosterlitz-Thouless model) which has been already studied in Ref. [17], but we only wish to call the readers attention on a few tiny numerical errors which have crept in the HT series coefficients at orders 17-20 reported in [17] and were due to an accidental contamination of a numerical file in the final stage of that calculation. In this paper we have corrected such errors, which of course could be annoying to those who wish to extend the computation and first have to make sure that they are able to reproduce correctly the existing data. We have also checked that, being so small, these errors were of no consequence at all either on the qualitative and on the quantitative results of the analysis in Refs. [17][18][19][20][21], which therefore does not need to be repeated until significantly longer series and/or better methods of analysis can really offer new insights, for example on the questions raised by Ref. [22].
In section IV we examine the set of models with N > 2, which are expected to behave quite differently [23]. In the last two decades their features have been extensively explored by various analytical and numerical techniques, with the main motivation that they are lattice regularizations of the field theoretic non-linear O(N)-symmetric σ-models, which share the crucial asymptotic freedom property with four dimensional gauge field theories, but are much easier to study. We update our previous survey of the nearby singularities of χ, µ 2 (and χ 4 ) in the complex inverse temperature β plane [9] and we still find no indication of any physical critical point at finite real β. On the contrary, we point out that the low temperature behavior appears to join smoothly onto the high temperature behavior so that several parameters which characterize the low temperature behavior can be computed in terms of HT series and full consistency is obtained with the predictions of the perturbative RG.
We end the paper by comparing our conclusions to some related recent works which, either by direct stochastic simulations or by analytic approximations such as the 1/N expansion, also test and confirm the predictions of RG.
In the appendices we report the closed form expressions for the HT series coefficients of χ, µ 2 and χ 4 as functions of the spin dimensionality N and their evaluation for a few specific values of N. Electronic files containing these data may be requested from the authors. The present tabulation extends and supersedes the one in Ref. [9] which, unfortunately, is marred by a few misprints.

II. DEFINITIONS AND NOTATIONS
The Hamiltonian H of the N-vector model is: where v(x) is a N-component classical spin of unit length at the lattice site x, and the sum extends to all nearest neighbor pairs of sites.
The HT expansion coefficients of all correlation functions C(N; β) = r f r (N)β r are simple rational functions of N, namely f r (N) = P r (N)/Q r (N) where P r (N) and Q r (N) are integer coefficient polynomials in the variable N. Therefore the same property is true for the expansion coefficients of the susceptibility of the second correlation moment of the second field derivative of the susceptibility It should be noticed that our definitions of χ and µ 2 differ by a factor 1 N and our definition of χ 4 differs by a factor 3 N (N +2) from that of Ref. [7] (apart from misprints, also adopted in Ref. [9] as far as χ 4 is concerned).
Like in any calculation of HT series, the correctness of the numerical results is a decisive issue. Our confidence on the validity of this work is based not only on the numerous direct and indirect internal tests passed by our code, but especially on the fact that the space dimension d and the spin dimensionality N enter simply as parameters which can be varied to produce, always by the same procedure, series in complete agreement with those already available (sometimes to a higher order) for specific values of N and d and in the spherical model limit (namely the limit as N → ∞ at fixedβ = β/N) [1]. More precisely let us observe that from the tabulation in the Appendix it appears thatã r (N) ≡ N r a r (N), s r (N) ≡ N r s r (N), andd r (N) ≡ N r d r (N) have a finite limit as N → ∞, so that the spherical model limit of our series coefficients can be immediately read. The HT coefficients of the susceptibility χ (s) (β) and the second correlation moment µ (s) 2 (β) of the spherical model can then be used to check in our tables the coefficients of the highest power of N in the numerator polynomials of a r (N) and s r (N) respectively. It should be noticed that χ 4 is O(1/N) for N → ∞, as it is expected because, in the spherical model limit, only the two-spin connected correlation functions are nonvanishing. However the quantityχ 4 (N, β) ≡ Nχ 4 (N, β) has a finite large N limit which will be denoted by χ (s) 4 (β). We remind the interested reader that, in Ref. [24], we have tabulated the HT coefficients for the spherical model susceptibility χ (s) on the square lattice through order 63, and that the HT coefficients for µ ). Similarly also the N = 0 limit, which corresponds to the self-avoiding walk (SAW) model can be easily obtained from our HT series after expressing all quantities in terms ofβ. Again, it is the quantityχ 4 (N,β) that has a finite limit for N = 0. The HT coefficients are then essentially given by the constant terms in the numerators of a r (N), s r (N), and d r (N), up to simple factors from the denominators.
More (partial) checks are trivially obtained by setting N = 1, 2 or 3 in our formulas and comparing the results to the corresponding available expansions cited below.
Of course, if a "complete set" of such checks were available, it could be used to reconstruct the whole computation.
It is interesting to recall that, more than two decades ago, HT series valid for all N and for a general lattice have been computed up to order β 8 ( β 9 for loosely packed lattices) [26]. Later on, in the case of a square lattice, the series were extended through β 11 [8] and then through β 14 [7,9]. On the other hand, very long expansions, on the square lattice, have been computed recently, both for the susceptibility and for the second correlation moment in the special cases N = 0 [27,28] (the self-avoiding walk (SAW) model), and N = 2 [17,19]( the Kosterlitz-Thouless model), by highly efficient algorithms, whose performance, however, does not excel for space dimensionality larger than two or which cannot be extended to other values of N. More precisely, for the susceptibility, the published series extend through orders β 43 (recently pushed to β 51 [29]) for N = 0 and β 20 for N = 2, and for the second correlation moment through β 27 and β 20 , respectively. The longest published expansions of χ 4 , valid for any N, presently extend only through β 14 [9]. The N = 1 case (the spin 1/2 Ising model), which is much simpler because the model is partially solved, should be considered separately: in this case the available series for χ and µ 2 extending to β 54 are tabulated in Ref. [30] while the series for χ 4 [31] extends only to β 17 . When our work was being completed, another calculation valid for any N was announced [32] for the nearest neighbor correlation function, χ and µ 2 (but not for χ 4 ) by the technique of group character expansion. This procedure seems to be efficient only in 2 space dimensions and to be presently feasible up to order 21 on the square lattice, to order 30 on the exagonal lattice and to order 15 on the triangular lattice [33]. It is reassuring that our general results, although obtained by a completely different procedure, agree throughout their common extent also with the specific cases N = 2, 3, 4, 8 tabulated in Ref. [33].

III. ANALYSIS OF THE HT SERIES FOR N < 2
We will now discuss some of the information that can be extracted from the series and update the analysis first presented in [9].
Let us recall that an exact expression for the critical exponent ν(N) within the range −2 ≤ N < 2 has been conjectured in Ref. [14] on the basis of an approximate RG analysis. Later on the same expression and an analogous one for the exponent η(N) were derived [15] by observing that a special O(N) spin model (assumed to be a faithful representative of this universality class) can be mapped into a soluble [16] loop gas model.
In order to test numerically the validity of (5) and (6), we have estimated γ(N) and ν(N) by forming first order inhomogeneous differential approximants (DA) [34] of the susceptibility χ and of the "second moment" correlation length squared ξ 2 = µ 2 /2dχ respectively. We also have computed the gap exponent ∆(N) from χ 4 . Our numerical procedure consists in averaging over all estimates from DA's in the class selected by the protocol of analysis of Ref. [35] which use at least 16 series coefficients in the case of χ and ξ, and at least 14 coefficients in the case of χ 4 . For each value of N, we first estimate β c (N) and γ(N) from χ, then we use β c (N) to bias the computation of ν(N) from ξ 2 and of ∆(N) from χ 4 /χ. As a measure of the uncertainties we have taken three times the rms deviation of the approximant estimates.
Alternatively, we have assumed the validity of the hyperscaling relation for −2 ≤ N < 2 and used also the series for χ 4 in the computation of ν(N) by resorting to the so called "critical point renormalization" [34,36]. In this case we have estimated ν(N) by examining the singularity at z = 1 of the series r h r (N)z r with coefficients h r (N) = d r (N)/t r (N) where χ 2 (N; β) = r t r (N)β r . Similarly the exponent ∆(N) has been determined from the series with coefficients l r (N) = d r (N)/a r (N) and the exponent γ(N) has been obtained in terms of the series for µ 2 and for ξ 2 . This procedure does not require the knowledge of β c (N), but only seventeen term series are available for the computation of ν(N) and ∆(N).
In the central region of the plot, approximately for −1 < N < 1.5, both numerical procedures we have followed yield very accurate estimates agreeing with the exact formulas within a small fraction of a percent. Near both ends of the interval −2 < N < 2 the agreement deteriorates because the series have to crossover to different singularity structures in order to exhibit either a gaussian behavior for N = −2 or a Kosterlitz-Thouless behavior for N = 2. However, since the exponent estimates always move in the right direction as the number of series coefficients is increased, we are confident that, in these border regions, we are simply facing a numerical approximation problem rather than a breakdown of the exact formulas (5), (6) and therefore we can conclude that their validity as well as the validity of the hyperscaling relation (7) are convincingly supported also by our HT series study.
We can also estimate with fair accuracy the (nonuniversal) critical amplitudes ofχ 4 for N = 0 and N = 1 which might be useful for reference and comparison with other numerical calculations. Let us recall that in the Ising model case the critical amplitude of the susceptibility c χ (1) has been computed exactly to be c χ (1) = 0.962581732.. and that also the amplitudes of the first few subleading confluent corrections to scaling are known [37]. Since the first confluent corrections are found to be analytic, and indeed it has been argued [38] that there are no irrelevant-variable corrections to scaling in the thermodynamic quantities of the two-dimensional Ising model, we expect that we can rely quite simply on near diagonal Padé approximants (PA's) of (β c − β) γ(1) χ to obtain an accurate estimate of c χ (1), even from not particularly long series. Indeed, from our O(β 21 ) expansion, we get c χ (1) = 0.96261(3) using as a bias the known values of γ(1) and β c (1). The critical amplitude ofχ 4 has not yet been evaluated exactly, but since the structure of the confluent corrections to scaling should be similar to that of χ also this amplitude should be accurately estimated. Our biased estimate of this quantity isĉ 4 (1) = 4.378(2) which compares well with the estimateĉ 4 (1) = 4.37(1) from the fifteen term series of Ref. [39]. Also in the N = 0 case ( now in terms of the variableβ) the structure of the confluent corrections is likely to be favorable since both the long expansion computed for χ in Ref. [28] and the results of extensive stochastic simulations [40] are consistent with an analytic dominant confluent correction. In this case we get c χ (0) = 1.0524(8) andĉ 4 (0) = 6.62 (2).
In terms of χ, ξ, andχ 4 we can also compute the "dimensionless renormalized four point coupling constant"ĝ r (N) as the value of at the critical point β c (N). If we assume that the inequality γ(N) + 2ν(N) − 2∆(N) ≥ 0, (rigorously proved to hold as an equality for N = 1), is also true for any −2 ≤ N < 2 , thenĝ r (N) is a bounded (nonnegative) universal amplitude combination whose vanishing is a sufficient condition for gaussian behavior at criticality, or, in lattice field theory language, for "triviality" of the continuum limit theory defined by the critical N−vector model [25]. Notice that our normalization ofĝ r (N) is the same as the one adopted in Ref. [7] and differs by a factor N +8 8πN (N +2) from the normalization traditionally chosen in the field theoretic renormalization group treatments [43].
For 0 ≤ N ≤ 2 we have estimatedĝ r (N) by evaluating both near diagonal PA's and first order inhomogeneous DA's of the series for 1/ĝ r (N, β) at the critical values β c (N). The two procedures yield results which are perfectly consistent within their numerical uncertainties. In Fig. 2 we have reported our estimates ofĝ r (N) for various values of 0 ≤ N < ∞ and compared our results to other computations in the literature.
It should be noticed that, as a reflection of the growing complexity of the critical singularity structure, the uncertainty of our numerical results is very low for N = 1 and relatively modest for N = 0, but it is much larger in the N = 2 case. However our new estimates appear to be generally more accurate than previous ones. The corresponding calculation for N > 2 will be discussed in next section.

IV. ANALYSIS OF THE HT SERIES FOR N > 2
In this range of values of N the general features of N-vector model are expected to change qualitatively: reliable, although nonrigorous (and sometimes questioned [46,47]), RG calculations at low temperature [23] indicate that the model is asymptotically free, namely that it becomes critical only at zero temperature. The asymptotic behaviors of the " second moment" correlation length, of the susceptibility and of χ 4 as β → ∞ are predicted to be are the first (renormalization scheme independent) low temperature perturbation expansion coefficients of the RG beta and gamma functions [23] and c ξ (N), c χ (N) and c 4 (N) are universal quantities which clearly cannot be computed in (low temperature) perturbation theory. The (nonuniversal) constants H 1 , H 2 , K 1 and K 2 appearing in (9) and (10) can be calculated in low temperature perturbation theory, and, on the square lattice [48,49], they come out rather small but not completely negligible in the range of values of β in which we shall be able to compute reliably χ and ξ . Numerical estimates for H 1 , H 2 , K 1 and K 2 can be found in Ref. [49] and for brevity are not reported here, although we use them in the calculations. Unfortunately, the analogous O( 1 β ) corrections have not yet been computed for χ 4 . As a consequence of (9), (10) and (11), for large β, whereĉ 4 (N) ≡ Nc 4 (N). Let us notice that the asymptotic formula (9) is valid both for the "second moment" correlation length ξ(N, β) (which here is employed exclusively), and for "exponential correlation length" ξ exp (N, β) but, a priori, with different multiplicative constants c ξ (N) and c exp ξ (N), respectively. However, it has been repeatedly noticed that ξ exp and ξ are numerically very close in the critical region [50]. This fact is confirmed by a recent analytic calculation for large N of the (universal) ratio c exp Moreover, for N = 3, this ratio has been measured [52] by a high precision MC method at β = 1.7 and 1.8, fully confirming the quantitative reliability of the 1/N expansion (14) down to very low values of N . Therefore, with very good approximation, we are justified in simply identifying c exp ξ (N) and c ξ (N) even for small N. These results are of direct interest here because the coefficient c exp ξ (N) can be computed exactly [53] by the thermodynamical Bethe Ansatz and its value, with our normalization conventions, is Let us now turn to series analysis. With sufficiently long series, like those analyzed here, even the simple plot of the HT coefficients of χ, versus their order, reported in Fig.3, is suggestive enough that the series for N > 2 must behave quite differently from those for N ≤ 2. Notice that in the plot the coefficients have been conveniently normalized, as indicated in the figure caption, in order to make the behaviors for different values of N easily comparable. In the N ≤ 2 case the coefficients are positive (and monotonically increasing) in agreement with the fact that the nearest singularity is located on the real positive β axis (and that the antiferromagnetic singularity located at −β c is much weaker). On the contrary, in the N > 2 case the coefficients do not remain positive as their order grows and they display what can be safely interpreted as the onset of an oscillatory trend. This feature is related via Darboux theorem [54] to the fact that, for N > 2, the nearest singularities of χ in the complex β plane become unphysical, as we have first pointed out some time ago [9]. This applies also to µ 2 , as well as to χ 4 with the only difference that for these quantities the oscillating behavior of the expansion coefficients should set in at higher orders. Our interpretation of these general features is impressively confirmed by a study [24] of the spherical model which, in spite of a priori legitimate mistrust about exchanging the large β and the large N limits, turns out to be a completely reliable guide to the qualitative behavior of the N−vector model even for not too large N ≥ 3, as it has been also argued some time ago [9,55].
As already mentioned above, an arbitrarily large number of HT expansion coefficients for χ (s) , µ (s) 2 and χ (s) 4 of the spherical model can be easily computed, for any lattice. In the case of the square lattice they exhibit regular cyclic alternations in sign of period 8 related to the symmetric quartet structure of the nearest unphysical singularities in the complex β plane. We have accurately mapped out in Ref. [24] the whole set of singularities, all of which are square root branch points. In the vicinity of β = ∞, this set has the characteristic structure dictated by asymptotic freedom (which was first discussed for the case of QCD in Ref. [56]), the analiticity domain of χ (s) , being a wedge with zero opening angle which contains the real β axis. It is quite likely that these features of the spherical model persist also down to all finite N ≥ 3 [57], although a complete study of this question is presently infeasible.
The transition from the N < 2 regime characterized by a conventional power law critical point, to the N > 2 asymptotically free regime characterized by unphysical singularities can be closely followed by locating the position of the nearest singularity β n in the first quadrant of the complex β plane as a function of N. For convenience in the graphical representation of the results, we shall use in what follows the scaled variableβ ≡ β/N and plot the estimates ofβ n versus x ≡ 1 − 1/N, rather than versus N. The trajectory of the singular pointβ n as a function of x in the complexβ plane can be traced out as described in Ref. [9] either by using PAs to locate the nearest singularity of the log derivative of χ or by directly computing DA's of χ. The results of both procedures agree perfectly within the numerical uncertainties.
In Fig.4 we have plotted the real and the imaginary part ofβ n as functions of x in the range 0 ≤ x ≤ 1. For −2 ≤ N ≤N < ∼ 2.2, the singularityβ n is still a real critical point, but for N >N it splits into a pair of complex conjugate singularities which move into the complex plane and, as N → ∞, tend to the limiting pointsβ ± ≃ 0.32162(1±i). In particular for N = 3 the nearest pair is located atβ n = 0.58(5) ± 0.14(5)i, while for N = 4 we havẽ β n = 0.55(4) ± 0.22 (5)i, for N = 6 we haveβ n = 0.50(2) ± 0.28 (2)i and for N = 10 we haveβ n = 0.44(2) ± 0.31(1)i. Although in general they may be weaker, the corresponding antiferromagnetic singularities will follow trajectories symmetrical with respect to the Im(β) axis so that, for all N >N the set of the nearest singularities will form a quartet with the same symmetry.
It is certainly conceivable [46] that, in contrast with the perturbative RG predictions, when the nearest singularities become complex, a further real critical singularity might appear so that, even for N ≥ 3, it would be still possible to relate the steep dependence on β of χ, ξ and χ 4 to a conventional finite temperature phase transition, but we have not been able to find any numerical indication of such a possibility for not too largeβ. More precisely neither Dlog PA's nor DA's exhibit any real and numerically stable singularity in their range of sensitivity. Another argument against the existence of critical points for finite values of β comes from the observation illustrated below, that, by a simple procedure, the high temperature behavior of ξ, χ and χ 4 can be smoothly extrapolated onto their low temperature behavior (9), (10) and (11) as predicted by the perturbative RG. This is feasible for any N, although the procedure is numerically very accurate only for N > 3, since, for N = 3, Im(β n ) is small and therefore the behavior of χ or ξ on the realβ axis is more strongly perturbed in the vicinity of Re(β n ). Both the often reported failure in observing asymptotic scaling by MC simulations of the N = 3 model at moderate values ofβ and the better successes for larger values of N find a completely plausible explanation in this picture. Of course, the results of our extrapolation scheme would be difficult to explain if the high temperature region were separated by a critical point from the low temperature region.
Let us now describe an approximation scheme which enables us to estimate low tem-perature perturbative parameters such as b 0 (N), as well as nonperturbative parameters like c ξ (N), c χ (N) etc. entering into the asymptotic formulas (9), (10) and (11), in terms of our HT series. Since ξ, χ etc. are exponentially fast varying quantities at large values of β, neither PA's nor DA's are well suited for a straightforward extrapolation of the HT series from small to (relatively) large β values. We should rather work with quantities which vary slowly enough to be well represented by PA's or DA's. Let us observe that, if χ has the asymptotic behavior (10), then for large enoughβ ). (16) The log derivative of χ, which is a slowly varying quantity, can be approximated by near diagonal PA's and then we can reliably extrapolate the quantity B χ (N,β). In practice, due to the finite extension of our series and to the intricate analytic structure of χ, we do not expect that this is a good approximation for largeβ and we rather make the reasonable (and successful) assumption that theβ independent parameter b 0 (N) is best approximated by evaluating B χ (N,β) at some finite real valueβ =β s where it is stationary or it shows the slowest variation whenβ is varied. Consistency of this approximation scheme requires that, as the number of HT coefficients used in the calculation is increased, the stationary value B χ (N,β s ) stabilizes and thatβ s → ∞. It can be checked that this actually happens in the N = ∞ case in which arbitrarily long HT expansions can be studied and also that our approximation scheme converges rapidly to the expected result.
A further check of the correctness of our procedures comes from the obvious remark that similar estimates of the same parameter b 0 (N) should be obtained starting either with the correlation length ξ and computing the quantity ). (17) or starting with χ 4 and computing the quantity In Fig.5 we have plotted B χ (N,β) versusβ for various values of N, in order to show that a stationary pointβ s (N) actually exists aroundβ ≃ 0.55 for any N, and that the size of the neighbourhood ofβ s (N) where B ξ (N,β) varies slowly withβ, grows with N. Notice that Re(β s (N)) < ∼ 0.55 for N > 4 and therefore on the border of the convergence region of the series or slightly outside it.
In Fig.6 we have plotted versus N our numerical estimates of b 0 (N) from the quantities B χ (N,β s (N)), B ξ (N,β s (N)) and B 4 (N,β s (N)) and have compared them to the expected value (12). Each point represents the average of the near diagonal PA's using at least 14 series coefficients for the quantities B χ (N,β s ) and B ξ (N,β s ), and at least 10 coefficients for B 4 (N,β s ). We have reported relative errors of 5% suggestive both of the scatter of the estimates obtained by the various PA's and of the systematic uncertainties of our extrapolation procedure.
In conclusion, it appears that from our high temperature data for ξ, χ and χ 4 , we have been able to extract completely consistent and correct estimates of the low temperature perturbation parameter b 0 (N) which characterizes the exponential asymptotic growth of these quantities, and in general that the deviation from the expected value (12) of b 0 (N) is never larger than 5% over a wide range of values of N.
In quite a similar way, assuming that b 0 (N) is given by (12), we can estimate the exponents of the power law prefactors in (9), (10) and (11). As it must be expected, the errors in this computation are somewhat higher, but they do not exceed 20 − 30%.
Let us now show that by a similar approximation procedure we can also estimate the constant c ξ (N). We have simply to compute the HT series of the slowly varying quantity obtained by dividing out the exponential factor in the asymptotic behavior (9) of ξ 2 . We then form near diagonal PA's to C ξ (N,β) and use them to evaluate the quantity at the valueβ s where it is stationary. In analogy with the previous computation this is taken to be the best approximation of c ξ (N) 2 . It is observed that also in this case the stationary values occur forβ ≃ 0.5.
Similarly, we can estimate c χ (N) by studying the HT series for the quantity ) .
Unfortunately no exact formula is known for c χ (N) , but we can compare our numerical estimates to the 1/N expansion through O(1/N) of c χ (N) which has been computed in Ref. [49] or to an analytic formula recently guessed in Ref. [61] with no other theoretical justification than a formal analogy with the exact formula (15) for c ξ (N).
An estimate of c 4 (N) could be obtained starting with the series for However the O( 1 Nβ ) corrections are not known, and moreover the 1/N expansion of c 4 (N) which has been computed in Ref. [44] it is practically useless, except for very large N, since the subleading term is quite large. Therefore we do not report our estimates for c 4 (N).
By the same method we have also directly estimated the universal quantity c r (N) which appears in the asymptotic expression of the ratio and have compared it to the 1/N expansion Let us notice that, for large N, the unphysical singularities of ξ 2 and χ tend to cancel in the ratio and that the 1/N correction in (26) is not very large.
We have reported our numerical estimates for c ξ (N), c χ (N) and c r (N) in Table 1 and in Fig.7 where they are compared to the exact or conjectured formulas and to their 1/N expansions. Like in the previous Fig.6, the error bars we have attached to our data points are fairly subjective in that they include a "statistical" contribution (describing the spread of the estimates from various approximants) which is not large in general, while the main part of the uncertainty comes from our estimate of the systematic error. As it appears from the Fig.7 and from Table 1, our central estimates for c ξ (N) and the exact formula agree within 1 − 2% on the whole range of N except for the lowest values of N. It should be observed that we have not our reported estimates for N = 3 since in this case the nearby unphysical singularities have a very small imaginary part Im(β n ) and there is a large spread in the stationary values of C ξ and C χ . This makes difficult to estimate unambiguously the values of c ξ and c χ . However, if we shift to only slightly larger values of N, such as N = 3.5, then Im(β n ) is already sufficiently large for our procedure to work appropriately and we can estimate c ξ (3.5) = 0.028 (8) to be compared to the exact value c exp ξ (3.5) = 0.0273 and, similarly, c χ (3.5) = 0.021(2), while the conjectured formula gives c exact χ (3.5) = 0.0266. In both cases the discrepancy is less than 20%. At N = 4 the exact value is c exact ξ (4) = 0.0416 and we find c ξ (4) = 0.039(1), which is off only by 6%. For larger N the agreement is much closer as it is shown in Table 1. Our estimates also agree well with the conjectured exact formulas for c χ (N) and c r (N). The discrepancy from these formulas or from their 1/N expansions does not exceed 5 − 10% for small values of N but it gets significantly smaller already for moderately large N. All numerical results are collected in Table 1. We believe that both the failure to reproduce accurately the N = 3 values of the parameters and the other general features of our approximations should not be surprising if we take into account the analytic structure in theβ complex plane of the quantities to be extrapolated and we consider that our computational method is the simplest and most direct possible and also that we are still working at moderate values ofβ where, for small values of N, the correlation length is not very large. Our approximation procedures should not however be suspected to be "ad hoc", since they were proposed and the first results [9] were published before either the exact formula (15) and the 1/N expansions became known.
We should also at this point recall our remark [24] (resumed and applied to fourteen term series in [62]) that the precision of these estimates might be significantly improved by performing a conformal transformation of the complexβ plane in order to remove at least the quartet of the nearest unphysical singularities before applying our approximation procedure. We think, however, that the success of our straightforward treatment over a wide range of values of N cannot be accidental and that simply getting a higher level of accuracy could hardly be more convincing of the validity of the RG picture of scaling and indicative of the purely numerical origin of the discrepancies for the lowest values of N. Therefore we shall not pursue here our old suggestion.
Finally, we can also estimateĝ r (N) for N ≥ 3, by forming PAs to the series expansion of 1/ĝ r (N, β) and evaluating them at their stationary points. We have reported in Fig.2 our estimates and have included for comparison the field theoretic estimate [63] for N = 3, and MC estimates [63,64,42] and other HT series estimates [44] for N = 3 and N = 4. It should also be noticed that our results are entirely consistent with the large N limit, in which we haveĝ The 1/N correction has been computed recently [44]. Also the accuracy of this calculation is satisfactory and the maximum error, for N ≥ 3, can be rated not to exceed 5%. Results and conclusions in complete agreement with ours are reached in the somewhat different analysis of the HT series presented in Ref. [33].

V. CONCLUSIONS
We have presented our estimates of the low temperature quantities b 0 (N), c ξ (N) etc. defined by (9), (10) and (11), obtained by a procedure which can essentially be seen as a simple improvement of the "matching method" introduced long ago in Ref. [65] and since used several times with more or less unconvincing results, due either to inadequate implementation and/or to incorrect supplementary assumptions. The initial paper [65] is an example of the former defect: the low temperature behavior was inadequately accounted for by one loop perturbation expansion and, on the HT side, too short series were used resulting into an unreliable matching. On the other hand Ref. [59] is an example of both shortcomings since the use of HT series (at that time extending to ten terms only) was supplemented with the (now appearing obviously incorrect) conjecture that χ and µ 2 have all positive HT coefficients. Indeed even if we made the weaker assumption that there are at most finitely many negative expansion coefficients this would clearly imply that the nearest singularity of (for example) χ is located on the real positive β axis. If also asymptotic freedom holds, then χ should be a regular analytic function in the whole finite complex β plane, contrary to the numerical evidence presented in the previous section.
We have tried to avoid the shortcomings of the previous approaches by the simplest possible treatment of sufficiently long HT series and by excluding unwarranted supplementary assumptions.
A brief review of some earlier references which are closely related to our analysis already appears in [9]. Here we shall mention only some later studies and address the reader to Ref. [9] for a long (but surely still incomplete) list of the abundant prior literature.
It is worthwhile to recall that recently, in the N = 3 case, a new method for extrapolating finite volume MC data to infinite volume [51] has been used to test the onset of the asymptotic behavior (9) by obtaining the second moment correlation length ξ up to values as large as 10 5 lattice units and agreement with (9) has been found within 4% atβ = 1. where ξ ≃ 10 5 . New MC data are now available [61] also for the susceptibility which yield Analogous results for the correlation length in the N = 3 case had been presented also in a computation [66] extending to ξ ≃ 15000. In an earlier high precision multigrid MC study devoted to the N = 4 model on lattices of size up to 256 2 [50], the asymptotic behavior of ξ and χ had been found to be perfectly compatible with (9) and (10) and the quantities c ξ (4) and c χ (4) had been estimated to be c ξ (4) = 0.0342 (20) and c χ (4) = 0.0329 (16). Moreover in that study the possibility of an ordinary critical pointβ c ≤ 1.25 was excluded, and it was stressed that the data could be compatible with a valueβ c ≥ 1.25 only assuming implausibly large values for the critical exponents. ( A far away power singularity with a large exponent is likely to be merely an effective representation of an exponential behavior.) Also the MC single cluster simulation of Ref. [67] for the N = 4 and N = 8 models gave good support to the asymptotic formulas (10) and (9) and produced estimates for c ξ completely consistent with (15).
Finally, on the side of the analytic approaches, we should mention the study of the scaling behavior in Refs. [68], whose results include a computation of the leading term of the 1/N expansion of c ξ , in complete agreement with the exact result (15), and of the same expansion for c χ .
In conclusion, we can summarize our main results as follows: a) By this and previous work [9] we have shown that our general N HT series are a useful tool also for obtaining high precision estimates of critical parameters in somewhat unconventional contexts, giving further support to qualitative and quantitative results obtained by entirely different approximation methods. b) In the −2 ≤ N < 2 vector models case we have confirmed, with high accuracy, the explicit formulas obtained by (semirigorous) model solving, for the critical exponents γ(N), ν(N) and ∆(N). We have also computed the "dimensionless renormalized four point coupling constant"ĝ r (N) for N = 0, 1, 2 in complete agreement with other estimates, but with higher accuracy. c) For the N > 3 vector models, we can somehow extrapolate the HT series to the border of (or beyond) their region of convergence reliably enough to reconstruct the quantitative features of low temperature behavior and we can obtain a set of (hardly accidental) consistency checks with the predictions of the perturbative RG, of exact solutions and of 1/N expansions with an accuracy practically uniform with respect to N. As shown by Table 1, our estimates of the parameter c ξ (N) agree well with the exact calculation by the Bethe Ansatz (under the assumption (15)). On the other hand our estimates for c χ (N), c r (N) and g r (N) are completely consistent with their 1/N expansions.
Of course we must say that, strictly speaking, purely numerical computations cannot validate the RG predictions: only complete proofs can settle the question, but they are still to come. Therefore, in principle, further discussion of this subject may still be considered healthy and welcome as long as it may stimulate either to design a rigorous justification of the generally accepted RG picture or to produce viable mechanisms for evading the expected asymptotic freedom regime while respecting the now established heuristic evidence. We believe, however, that the continuing accumulation of unambiguous, consistent and increasingly accurate numerical support for the RG predictions from a variety of independent approaches leaves little if any space for alternative pictures.

ACKNOWLEDGMENTS
This work has been partially supported by MURST.  N)) For the coefficients which follow it is typographically more convenient to set a r (N) = P r (N)/Q r (N) and to tabulate separately the numerator polynomial P r (N) and the denominator polynomial Q r (N).