Precise extrapolation of the correlation function asymptotics in uniform tensor network states with application to the Bose-Hubbard and XXZ models

We analyze the problem of extracting the correlation length from infinite matrix product states (MPS) and corner transfer matrix (CTM) simulations. When the correlation length is calculated directly from the transfer matrix, it is typically significantly underestimated for finite bond dimensions used in numerical simulation. This is true even when one considers ground states at a distance from the critical point. We introduce extrapolation procedure to overcome this problem. To that end we quantify how much the dominant part of the MPS/CTM transfer matrix spectrum deviates from being continuous. The latter is necessary to capture the exact asymptotics of the correlation function where the exponential decay is typically modified by an additional algebraic term. By extrapolating such a refinement parameter to zero, we show that we are able to recover the exact value of the correlation length with high accuracy. In a generic setting, our method reduces the error by a factor of $\sim 100$ as compared to the results without extrapolation and a factor of $\sim 10$ as compared to simple extrapolation schemes using bond dimension. We test our approach in a number of solvable models both in 1d and 2d. Subsequently, we apply it to 1d XXZ spin-$\frac32$ and the Bose-Hubbard models in a massive regime in the vicinity of BKT critical point. We then fit the scaling form of the correlation length and extract the position of the critical point and obtain results comparable or better than those of other state-of-the-art numerical methods. Finally, we show how the algebraic part of the correlation function asymptotics can be directly recovered from the scaling of the dominant form factor within our approach. Our method provides the means for detailed studies of phase diagrams of quantum models in 1d and, through the finite correlation length scaling of projected entangled pair states, also in 2d.


I. INTRODUCTION
Tensor networks and related numerical renormalization group techniques allow to efficiently approximate systems of exponentially many degrees of freedom with manageable number of a few relevant ones, providing invaluable tools in the studies of strongly correlated systems. We are particularly interested in two such techniques. The first one is based on the matrix product state (MPS) representation of a many-body wave function [1]. It provides the underlying framework behind a family of state-of-the-art methods for approximating the low-energy states of local one-dimensional Hamiltonians [2][3][4], descendants of the seminal density matrix renormalization group (DMRG) algorithm [5,6]. The second one, closely related, is the corner transfer matrix (CTM) algorithm [7,8]. It is used to numerically solve classical systems in 2d and is a method of choice for contracting 2d quantum states described by the projected entangled pair states (PEPS) ansatz [9,10].
In this article we address the problem of precise extrapolation of the correlation length in such simulations. Correlation length is a fundamental quantity in the description of (quantum) many-body systems and their phase transitions. It provides valuable input into the nature of the phase being described, as well as informs about its boundaries. Apart from this general consideration, there are two immediate applications of the presented method that we mention below. Firstly, it is desirable to establish a reliable method to calculate the correlation length of PEPS using the CTM algorithm. It can be applied to characterize the critical behavior of 2d quantum systems within the PEPS approach and set up finite correlation length scaling both at zero [11,12] and at finite temperature [13]. A similar strategy can also be used in studies of classical systems in 3d [14]. Secondly, it allows to fit the scaling form of the correlation length in the vicinity of the Berezinskii-Kosterlitz-Thouless (BKT) critical point in one-dimensional Bose-Hubbard and XXZ type models. Besides precisely extracting the critical point position, the accurate knowledge of the scaling form is relevant, for instance, to understand the behavior of quantum fidelity in such systems [15].
Below we focus the discussion on MPS and notice that the argument could also be directly applied to CTM. It is well known that MPS with finite bond dimension is able to reproduce the ground state of the local gapped Hamiltonian up to an error which vanishes exponentially with the bond dimension [16]. Consequently, local observables can usually be simulated with similar precision. The correlation length, on the other hand, describes the tail of the correlation function. This tail is vanishing exponentially and as such there is no reason to expect that MPS would be able to capture it faithfully. More importantly, asymptotics of the connected correlation function calculated for MPS with finite bond dimension is purely arXiv:1801.08554v3 [cond-mat.str-el] 3 Dec 2018 exponential [1] while typically one expects that the exponential decay is modified by an additional algebraic term, as, e.g., in the Ornstein-Zernike formula for the correlation function in the context of the Ising-type models [17,18]. For those reasons, it is not straightforward to faithfully recover the asymptotics of the correlation function directly from MPS simulations with finite bond dimension. To highlight the problem, for the specific point in the XXZ spin-1 2 model which we discuss in detail later, MPS with bond dimension 4096 recovers the exact ground state energy with an error of the order of 10 −12 , but at the same time it is still underestimating the correlation length by a factor of 2.
In this article we propose an extrapolation scheme to overcome the above problem. For uniform MPS (CTM) which describes a translationally invariant system the correlation length is calculated from the ratio of the two largest eigenvalues of the site-to-site (column-to-column) transfer matrix. The spectrum of this transfer matrix is necessarily discrete for the finite bond dimension used in the numerical simulations. In order to recover the algebraic part of the correlation function asymptotics in Eq. (2), the spectrum would have to be continuous. In our approach, we look at the distance δ between the next dominant transfer matrix eigenvalues, e.g. the second and the third one, and employ it as a measure of deviation from the exact solution. One expects to recover δ = 0 in the limit of the exact representation of the ground state. For a given model, we calculate the MPS correlation length, as well as the refinement parameter δ, for a number of MPSs with increasing bond dimensions and subsequently extrapolate δ → 0 in order to recover the actual value of the correlation length. In order to benchmark our approach we analyze a number of models where the correlation length, or some related properties like the position of the critical point and its universality class, are known analytically. Based on this data we argue that the method proposed in this article is more reliable and produces more accurate results than the one that directly uses the bond dimension as a refinement parameter.
Correlation function asymptotics can be derived from the Euclidian path-integral representation of the ground state and the exact quantum transfer matrix (QTM). In Refs. [19][20][21] it was argued that the MPS transfer matrix can be understood as an approximation of the QTM obtained as a result of the renormalization procedure akin to Wilson's numerical renormalization group description of the impurity problem [22]. In this picture, the physical spin is interpreted as an impurity in -by construction translationally invariant -QTM, and the MPS transfer matrix retains only the degrees of freedom (along the virtual, imaginary time direction of the system) relevant for the description of correlations of such an impurity. In this article we further build on this intuition and observe how the form factors (i.e., matrix elements of the operator transfer matrix in a suitable basis, defined below) are effectively being renormalized. Most importantly, we argue that the exponent η of the algebraic part of the correlation function asymptotics in Eq. (2) is directly related to how the relevant dominant form factor decays as we approach the exact solution, δ → 0 (in a limit of infinite bond dimension).
We should finally contrast our approach with the finite entanglement scaling scheme where, for the 1d system in the vicinity of the critical point, finite MPS correlation length (a result of finite MPS bond dimension) is used similarly to the effective finite size of the system in order to postulate a scaling hypothesis. The position of the critical point and critical exponents can then be extracted by proper renormalization and collapse of the data obtained for different bond dimensions [23][24][25][26][27][28][29]. Our method enables obtaining the scaling form of the correlation length without assuming the scaling hypothesis and as such can be used to independently corroborate some of the results found with finite entanglement scaling. We remark that obtaining the correlation length outside the critical regimes is beyond standard finite entanglement scaling schemes as it requires knowledge of a nonuniversal scaling function. Furthermore, our extrapolations can also be applied beyond the scaling regime, or more generally, when it might be impossible to postulate a scaling ansatz. In this context, for instance, it should be possible to apply our method to obtain a more precise description of nonequilibrium dynamical properties, such as spreading of the correlations in the excited system. The two methods are equivalent exactly at the critical point as we note that our refinement parameter δ is directly proportional to the inverse of the MPS correlation length in this case.
More importantly, both methods have to be combined in studies of quantum systems in 2d. In that case both the limited bond dimension of PEPS and the finite bond dimension of CTM used to contract it contribute to errors. Setting up finite correlation length scaling [11][12][13] (or finite entanglement scaling, as it is often called in the context of 1d systems) requires extrapolating the error of the correlation length resulting from CTM to zero in a controlled way, which is the goal of this article. As such, we anticipate mapping out phase diagrams of quantum systems both in 1d and 2d as an important application of our method.
The rest of the article is organized as follows. In Sec. II we introduce relevant notation focusing on MPS and quantify the general arguments from the Introduction. In Sec. III we briefly summarize our method discussing, in particular, how different properties of a studied model can be used to further refine the measure of error used for extrapolation. We benchmark our approach in Sec. IV. We study XY and XXZ spin-1 2 models, where the exact value of the correlation length is known. Then we focus on XXZ spin-3 2 model and conclude with the Bose-Hubbard model with unit filling as a nontrivial implementation of our method. In Sec. V we test our approach in the context of 2d models and their PEPS description. We employ corner transfer matrix method to analyze exactly solvable statistical models: the classical 2d Ising model and 8-vertex model. Finally, we apply our technique to thermal quantum states in 2d, including interacting spinless fermions, where the correlation length is not known otherwise. In Sec. VI, we discuss how the exponent of the leading algebraic part of the correlation function naturally emerges within our approach and uncover its connection with the form factors. We conclude in Sec. VII. Appendix A illustrates the problems related to fitting the asymptotic directly, where inaccurate results are obtained away from the critical point. Paradoxically, it is possible to obtain much better results for the critical systems within such an approach. This approach is indeed widely used in the literature. Finally, in Appendix B we argue that methods that provide good extrapolation of energy per lattice site are not suited to work well for the correlation length.

II. NOTATION
In this article we focus on infinite, translationally invariant systems. Setting up the notation, uniform matrix product states take the form where |s = | . . . , s 1 , s 2 , s 3 , s 4 , . . . and A sn are D×D matrices with parameter D usually referred to as the MPS bond dimension. The MPS transfer matrix (TM) is defined in a standard way It is the key object in the calculation of the static correlation function. In order to calculate the expectation values related to some operator o it is also convenient to define the operator transfer matrix For larger unit cell consisting of L sites, i.e. if MPS is translationally invariant only when shifted by L lattice sites (i.e. due to spontaneous breaking of translational symmetry), those L sites are combined into one to calculate the transfer matrix in Eqs. (4,5).
In our approach we focus on the eigenvalues of the transfer matrix T A , with j = 0, 1, . . . , D 2 − 1 and |λ 0 | > |λ 1 | ≥ |λ 2 | ≥ . . .. Eq. (6) singles out (minus log of) the absolute value and the phase into j ≥ 0 and φ j ∈ (−π/L, π/L], respectively. We have introduced the period L, so that the correlation length is measured in the units of lattice spacing. Note that some information about phase is lost for L > 1. Proper normalization of the state |Ψ entails that λ 0 = 1. We additionally assume that the largest eigenvalue of the TM is unique. This ensures the state |Ψ in Eq. (3) is properly defined. It is well known that the correlation length ξ A associated with given normalized MPS is set by the second largest transfer matrix eigenvalue as and φ 1 captures the leading period of oscillations of the correlation function. More precisely, the connected correlation function of operators o and q at distance R can be expressed as where the form factors f oq j are defined as with right |ϕ j ) and left (ϕ j | eigenvectors of the transfer matrix normalized as (ϕ i |ϕ j ) = δ ij . Eq. (8) leads to purely exponential decay of asymptotics of the correlation function as in Eq. (1), dictated typically by the second largest eigenvalue of T A . More generally, this decay is established by the largest eigenvalues (j > 0) for which the corresponding form factors f oq j are nonzero. The asymptotics of the correlation function, however, typically contains an algebraic factor as in Eq. (2). For the MPS to recover the algebraic part R −η of the asymptotics, the spectrum of the transfer matrix would have to be continuous. This is clearly impossible for the finite bond dimension used in numerical simulations. This is a well-known fact in the studies of critical points. Here, the decay of the correlation function is purely algebraic. This fact is used as one of the arguments of why simulating such points with MPS is particularly challenging. As argued above, those issues are still present also far away from the critical point in the context of precise extrapolation of the correlation length.
The discrete dominant eigenvalues of the transfer matrix can at best approximate the continuous spectrum related to the exact quantum transfer matrix, as presented pictorially in Fig. 1. As suggested in that picture, it can be expected that the correlation length obtained as ξ A = 1/ 1 underestimates the exact value, as 1 would be localized inside the band and not on its edge. Consequently, one has to resort to extrapolation in order to recover the true value of the correlation length.
An alternative approach would be to fit the asymptotics in Eq. (2) for the intermediate values of the distance R. However, away from the critical point, this Illustration of the idea behind the scheme. The figure represents (the logarithm of) the dominant part of the transfer matrix spectrum in a generic situation. The blue line represents continuous band necessary to recover the algebraic part of the correlation function asymptotics in Eq. (2). In this case the exact correlation length is set by the gap between the bottom of the band and the origin, 0 = 0. The spectrum of the transfer matrix for finite bond dimension MPS, represented here by red marks, is necessarily discrete and as such can only approximate the continuous band. Consequently, 1/ 1 is typically underestimating the true value of the correlation length. We employ δ = 2 − 1 as a natural measure of how well the discrete spectrum is able to approximate the exact continuous one. By computing 1(D) and δ(D) for some number of MPSs with different bond dimensions D, we extract the correlation length by extrapolating δ → 0.
would require the ability to fit the correlation function asymptotics for distances between the physical correlation length and the length scale that results from the discreteness of the MPS transfer matrix. As we illustrate in the Appendix A, such scale separation might not be accessible in the MPS simulations. The above approach becomes viable at the critical point where there is no physical length scale that has to be respected.

III. SUMMARY OF THE APPROACH
In this article we employ consecutive largest TM eigenvalues to quantify the divergence from the continuous spectrum necessary to capture the algebraic part of the asymptotics. In the simplest case we use the distance between the third and second eigenvalue, i.e., as a refinement parameter that measures the deviation from the exact solution. If needed, the above simple measure can be further refined by taking into account the fact that some form factors may vanish, the transfer matrix can be degenerate and that the system might display some symmetries. The above situations are summarized below and discussed in the subsequent sections of the article. We calculate 1 (D) and δ(D) for a few MPSs obtained for different bond dimensions D, where we observe that the dependence is usually smooth and regular. This allows us to extrapolate δ → 0 in order to extract the true value of the correlation length with good precision. We compare this approach with the one where 1/D is used as a refinement parameter. We observe that 1 is significantly less regular as a function of D -especially away from the critical point -and the result of the extrapolation is less reliable. The refinement parameter defined in Eq. (10) proves to be a good starting point, being sufficient in many simpler cases. However, one of the advantages of quantifying the distance from the exact solution using intrinsic quantities calculated for a given MPS approximation -in contrast to external parameter such as the bond dimension -is that we can easily take into account additional information about the state to further refine it if necessary. This allows to uncover additional information that the MPS description is carrying as well as increase the precision.
Firstly, let us focus only on the part of the transfer matrix spectrum relevant for some particular correlation function C oq (R). To that end, we take into account only those transfer matrix eigenvalues for which the corresponding form factors f oq j are nonzero (within numerical precision), or dominant as compared to the other ones. We mark such eigenvalues with a tilde and additional superscript,˜ oq k . In such cases, we define the refinement parameter as We apply this definition throughout the article as the most reliable indicator of which TM eigenvalues are relevant and should be taken into account. Secondly, dominant TM eigenvalues are usually found in groups with well-defined complex phases corresponding to periods of oscillation of the correlation function. The nontrivial correspondence between those phases and the minima of the dispersion relation of the Hamiltonian for which a given MPS is the ground state is discussed in Ref. 19. In order to define the refinement parameter, we can focus only on part of TM spectrum with a given complex phase ϕ: We discuss this approach further in Sec. IV A, where we study the incommensurate phase of the XY model. Thirdly, there are models for which the dominant TM eigenvalues are either degenerate or are effectively becoming degenerate with the increasing bond dimension. In such a case it is necessary to define the refinement parameter as where the eigenvalues 1 , . . . , n−1 are (near) degenerate. We observe that behavior e.g. for the XXZ spin-1 2 model in Sec. IV B, where 4 dominant eigenvalues are near degenerate and we use n = 5 in the definition above. More generally, even without degeneracy, any definition of the refinement parameter δ = n − 1 with n > 1 should lead to the same extrapolated value of the correlation length. In practice, however, we observe that using smaller n allows for more precise results.
Finally, if the system has some local symmetry and MPS is implemented to take it into account, then the TM spectrum splits into groups with well-defined symmetry charge u. We can define the refinement parameter using only the eigenvalues belonging to one symmetry sector: This is equivalent to selecting a particular subset of correlators corresponding to a given symmetry sector. We use the above approach for the XXZ and Bose-Hubbard models in Secs. IV B-IV D, which all have U (1) symmetry.
It should be pointed out that the above features are not independent and can be used simultaneously. Ultimately, the reliability of the extrapolation hinges on the consistency of the data obtained for different bond dimensions and through application of different refinement features.
Let us now comment on the extrapolation model that we employ. By analyzing a number of exactly solvable systems we observe that very good results are obtained if one extrapolates by fitting the function where 1/ e is the extrapolated value of the correlation length and where we assume that the error of the inverse of the correlation function, − exact , is vanishing as a power law with δ. The exponent b is usually slightly smaller than 1 and in many cases linear fit, i.e. fixing b = 1 in Eq. (15) proves to be sufficient. It is also a good starting point, which can then be further tested by allowing b = 1 and checking if this significantly improves the quality of the fit [30]. We use 95% confidence bounds from the non-linear fit in order to estimate an error of extrapolation. We observe that it provides a sensible measure of the quality of the result.

IV. MATRIX PRODUCT STATES SIMULATIONS
In this section we benchmark our approach in a range of models of increasing difficulty. As we move to more difficult models, we illustrate different ways of introducing refinement parameters as briefly discussed in the previous section.

A. XY model
We start with the one-dimensional XY model with anisotropy parameter γ and magnetic field g. This model is exactly solvable and the asymptotic form of the connected correlation functions is long known [31]. We cite the relevant results below. The numerical results in this section were obtained using the variational uniform matrix product states algorithm (VUMPS) of Ref. 32, with one-site unit cell, which is based on the time-dependent variational principle approach [33,34]. All the states for different bond dimensions were converged with the norm of the energy gradient below 10 −12 . Similarly, the maximal change of the Schmidt values in the last iterations of the algorithm (which is another strict measure of convergence) was of the same order.

Ising model
We start with the Ising model by setting γ = 1 in Eq. (16) and note that the TM spectrum is real and positive in this case. We collect the numerical results in Fig. 2 and discuss them below.
In the paramagnetic phase, for g > 1, the correlation functions behave asymptotically as where the inverse of the correlation length 1/ξ = exact = ln g. Note the additional factor of two in the exponential part of C zz (R) which is halving the correlation length that appears there. The results for g = 1.01 are shown in Fig. 2(a). Several observations are in order. Without extrapolation, even for the relatively large bond dimension D = 512 for which the smallest Schmidt value of the MPS bipartition is of the order of 10 −14 , the relative error of the correlation length is still 2%.
The dependence of 1 on δ = 2 − 1 , Eq. (10), is close to linear. This behavior is already seen for the smallest bond dimensions presented on the plot. Linear regression allows for extrapolation of the true correlation length to within relative error below 0.1%, which can be made even smaller by neglecting the smallest D shown in the picture. It is interesting to note that the eigenvalues 1 and 2 , used here to calculate the distance δ above, contribute to different correlators C xx (R) and C yy (R), respectively. Nevertheless, this approach proves to work well in this model. We obtain a consistent result and similar accuracy when we take into account nonzero form factors, and consider only the part of the TM spectrum that contributes to C xx (R), i.e. using Eq. (11). The above can be contrasted with direct application of the bond dimension as a refinement parameter, where the first natural choice is δ = 1/D. In this case 1 is oscillating as a function of 1/D, making extrapolation significantly less reliable as it becomes arbitrary which points to choose for extrapolation. Linear regression for D = 32-512 recovers the correlation length with relative error 1%, more than an order of magnitude worse than our approach. Indeed, by exploring solvability of the model and in particular its Schmidt spectrum, it was argued in Ref. 20 that 1 should be approaching the exact value much slower than linearly in 1/D. This ex-  Distance from exact solution δ (c) γ = 1; g = 1.5; ξ = 1/ exact = 2.4663 . . . plains why linear regression is still underestimating the true value of ξ -a feature which for sufficiently large D is shared by all the models studied in this article.

Distance from exact solution δ
In Fig. 2(b) we focus on parts of TM spectrum contributing to different correlators: C xx (R), C yy (R) and C zz (R). In the case of a paramagnetic Ising model we observe that each eigenvalue has exactly one form factor which is nonzero. It is either f xx j , f yy j or f zz j . This allows to recover the fact that the correlation length associated with C zz (R) is halved as compared to the other two, in agreement with the exact result. This shows that such information is encoded, and can be directly extracted from the MPS TM. It is worth noticing that the dominant part of the TM spectrum is relatively sparse: all points in Fig. 2(b) were obtained using the information from up to 11 largest TM eigenvalues and this was enough to distinguish and extrapolate two correlation lengths differing by a factor of two, even for the largest D = 512 used there.
In Fig. 2(c), we show the results for g = 1.5, illustrating that the problem with precise extrapolation is present even far away from the critical point when the correlation length is of the order of few sites only. Even in this simple case, without resorting to extrapolation, it is virtually impossible (as the smallest Schmidt values are falling below numerical precision) to recover the true correlation length with relative error below 0.75%. On top of that the value of the relative error, say for fixed D, clearly depends on the distance from the critical point -compare with Fig. 2(a) for g = 1.01. This makes any fits that use correlation length obtained directly from MPS with fixed D, e.g., extracting critical exponents or the position of the critical point out of it, much less trustworthy. Proper extrapolation, as suggested in this article, allows to significantly mitigate this problem.
In Fig. 2(d), we show results for the Ising model in the ferromagnetic phase with g = 0.99. In the regime 0 < g < 1, the correlation functions behave asymptotically as with the inverse of the correlation length 1/ξ = exact = −2 ln g. All the observations made for the paramagnetic phase above fully apply here as well.
Finally, in Fig. 2(e) we show the validity of the general extrapolation ansatz introduced in Eq. (15). The error, that is the distance between 1 and the exact value, is vanishing as a power law of the refinement parameter proposed in this article with the exponent close to 1. In this simple model, however, we observe that linear regression -as discussed in the text above -is sufficient and a more general power law does not result in qualitative improvement of the results in this case. Part of the reason might be that even without extrapolation relative errors are already small here -at least as compared to other, more complicated models discussed below. More importantly, the exponent b is very close to 1 here. When 1/D is used as a refinement parameter, the error of does not follow a simple functional form, as can be seen in Fig. 2(f). While it can be locally approximated by a power law, it is evidently flattening, making it a poor ansatz for extrapolation. We note that the critical point is a single exception here, as observed in the context of finite entanglement scaling [24][25][26][27].

Incommensurate ferromagnetic phase
In this section we focus on the incommensurate ferromagnetic part of the phase diagram of the XY model, g 2 + γ 2 < 1, where the correlation function is not vanishing monotonically but has an additional oscillating term.
For g > 0 and 0 < γ < 1, the leading asymptotics of the correlation functions is In this case the asymptotic behavior may be additionally modified by an oscillating term.
In Fig. 3(a) we show the full TM spectrum on a complex plane. The dominant eigenvalues form groups with well-defined complex phases [19], 0 and ±ϕ XY , respectively. They correspond to frequency of oscillations of the correlation functions. We enlarge the ϕ XY branch in panel (b). It can be seen that the exact phase ϕ XY is well reproduced in the simulation, especially for the dominant eigenvalue. In panel (c) we show the results of extrapolation using linear fit. In this model, the simple distance from Eq. (10) used without any additional refinement results in not-to-smooth functional dependence. Nevertheless, linear regression still reproduces an exact value up to 1%. The data are significantly smoother if one focuses only on the part of the TM spectrum with a ϕ XY complex phase, or additionally takes into account nonzero form factors. The error resulting from linear regression is, however, still ∼ 1%. All those approaches yield much better results than a linear fit as a function of 1/D for which the error is ∼ 4%. This is clarified in panel (d) where we present the data on a log-log plot, showing that the dependence of relative error on δ is better described by a power law with the exponent that is close to 1 in our approach. Indeed, selecting points that have the same complex phase ϕ XY and applying a nonlinear fit allows to reduce the error of extrapolation further to ∼ 0.3%.  In this section we analyze the results for spin-1 2 XXZ model
The numerical results in this and in the following sections were obtained using the iDMRG algorithm [39] with a two-site unit cell incorporating U(1) symmetry [40,41], which in this case corresponds to the conservation of  We collect the results in Fig. 4, where we focus on the part of the TM spectrum corresponding to U (1) charge u = 1 with the refinement parameter δ defined in this sector according to Eq. (14). We present the splitting of the dominant part of the TM spectrum for a single D in the inset of panel (a). Figure 4 shows the data in a log-log scale to highlight power-law dependence of the error of on δ. For completeness, in panel (b) we show the dependence of the error on 1/D. While it is smooth, it does not seem to follow a clear functional form again making it not very useful for precise extrapolation.
We also run simulations with VUMPS with o onesite unit cell [42] and without U(1) symmetry where we used the refinement parameter δ = 5 − 1 , as defined in Eq. (13), to take into account near degeneracy of 1 , . . . , 4 -see inset of Fig. 4(a). We present those results to show how to deal with degeneracies when they become an issue for the simplest refinement parameter in Eq. (10). Those results are represented by circles in   Finally, we collect results of the actual nonlinear fits in Table I. The values are averaged over range of bond dimensions taken into account, where we use 8 to 13 points with largest D. Additionally, we also take into account dominant form factors. This is especially relevant for the u = 0 sector where there is near degeneracy of two dominant eigenvalues, which still has to be resolved (see inset of Fig. 4(a)). Apart from ∆ = 1.1, the exact value of the correlation length is recovered with the error well below 1%. It is ∼ 3% for ∆ = 1.1, however the correlation length is approaching 10 4 here and we extrapolate from MPS correlation lengths underestimating the exact one by almost a factor of 2. Obtaining such results from the simulation of a finite system is practically impossible, showing the effectiveness of the infinite uniform approach.
We finish this section with an observation that the TM spectra which we obtain from numerics in this model are real. There are both positive and negative eigenvalues for the one-site unit cell implementation (complex phases ϕ = 0, π) which corresponds to the monotonic and staggered part of the correlation function asymptotics -see Eq. (3.36) in Ref. 35. This information about the phase is lost in the two-site unit cell transfer matrix, where the spectrum is effectively squared and strictly positive, which shows some advantage of using as small a unit cell as possible. It is interesting to contrast this with the spectrum of the quantum transfer matrix, which is complex with the complex phase changing in a continuous way; see e.g., Ref. 37 for an in-depth discussion. In that case the minimal gap of the quantum transfer matrix is not dictating the actual correlation length as contributions of part of the QTM spectrum band are effectively canceling out. Apparently, MPS is capturing only the physically relevant part of the spectrum making it purely real in our case. A similar situation -though slightly more compli- . In panel (b) we collect correlation lengths obtained from our extrapolation procedure. Subsequently, we fit the form of log(ξ(∆)) given in Eq. (19). This allows us to recover the exact position of the critical point, ∆c = 1, with excellent accuracy. See text for details. cated -arises in the 8-vertex model, which we discuss below in the context of the CTM algorithm.
C. XXZ spin- 3 2 model In this section, we consider spin-3 2 XXZ model where S x,y,z m are standard spin-3 2 operators acting on site m. Similarly to the spin-1 2 case discussed in the previous section, the model has a critical region for −1 ≤ ∆ ≤ 1 with a BKT critical point at ∆ c = 1 separating the gapped phase for ∆ > 1 [43][44][45][46]. Again, we focus on the latter, where the correlation length scales as Even though the exact value of the critical point is known, the model cannot be solved exactly and the value of the correlation length is not known analytically. Additionally, it is particularly challenging to approximate ∆ the ground state due to very strong quantum fluctuations [47,48]. As such, the model provides a good test for numerical methods. We extrapolate the correlation lengths using our method and subsequently fit the scaling form of Eq. (19) (with higher-order corrections) in order to extract ∆ c and the parameter B. We collect the results in Fig. 5 and a few selected correlation lengths in Table II. We proceed similarly as in the previous section. We define the refinement parameter by taking into account information about the symmetry sector, Eq. (14), and extrapolate by fitting the general power law in Eq. (15) to 8-13 points with the largest bond dimension D used in the simulation. Consecutive D form a geometric series with a step 2 1/4 approaching 10000 for the most challenging points. We show some of such fits in Fig. 5(a). This also includes the critical point at ∆ = 1, where we obtain the extrapolated value of the inverse correlation length as e ∼ 10 −6 ± 7 · 10 −6 , which is 2 orders of magnitude smaller then the value of 1 obtained for the largest D = 5792 converged here. This shows both the limits and validity of our approach, as within the estimation error we are able to recover the exact value (zero) at the critical point.
We compare those results with the ones recently reported in Ref. 47, which were obtained by fitting the scaling ansatz capturing behavior of the energy gap in the finite system. Namely, ∆ c = 0.995 ± 0.004 and B = 0.50 ± 0.02 (∆ c = 0.989 ± 0.01 and B = 0.58 ± 0.04) for open (periodic) boundary conditions, where the systems up to 280 (72) spins were used. We are able to access the range of correlation lengths which are an order of magnitude larger then system sizes possible in state-of-the-art finite system simulations and also take into account subleading correction to the scaling -which become increasingly important at a distance from the critical point. As such, we expect our results to be more accurate, which can be seen in the precision with which we were able to localize the critical point.  TABLE III. Extrapolated correlation lengths in the Bose-Hubbard model with unit filling, n = 1. We show the correlation lengths for b † 0 bR and n0nR . Notice that the second one is approximately halved as compared to the first one. Additionally, we show maximal bond dimension D used in the simulation and MPS correlation length for that D.

D. Bose-Hubbard model
We conclude this part with the Bose-Hubbard model in one dimension, where b m are bosonic annihilation operators acting on site m and n m = b † m b m is the particle number operator. Below we set the energy scale by fixing the Coulomb repulsion U = 1 and consider a system with unit filling per lattice site, n m = 1. The model has a quantum phase transition between the gapped Mott insulator phase for J < J c , and the gapless superfluid phase for J > J c in the Berezinskii-Kosterlitz-Thouless universality class [49,50]. We focus on the gapped phase and proceed identically as in the previous sections.
In the iDMRG simulations we truncate the local Fock space at six particles, checking that this is enough to obtain converged results. We employ U(1) symmetry, which in this case corresponds to conservation of the total particle number, m n m . This model proves to be less challenging for MPS simulations than the XXZ spin-3 2 model from the previous section and all the results were obtained with bond dimension up to 5792. We collect the results of our extrapolation procedure in Fig. 6 and in Table III. We calculate both the correlation length associated with b † 0 b R and n 0 n R correlators. We observe that ξ nn is halved as compared to ξ b † b within the estimated extrapolation errors; see Table III. Subsequently, we focus on ξ b † b , which can be extrapolated with higher accuracy, and fit ln ξ b † b = a 1 +B(J c −J) −1/2 +a 2 (J c −J) 1/2 . We use ξ ∈ [10, 5000], for which the estimated extrapolation errors are well below 1%. We obtain the position of the critical point as J c = 0.3048(3) and nonuniversal constant B = 1.61(4). Additionally, a 1 = −1.34 (15) and a 2 = −3.52 (24).
For the collection of results on the critical point position obtained in a multitude of different studies, see Table  1  To conclude, we comment that the range of estimates of the critical point position, as collected in Ref. 50, illustrates very well the complexity of such studies. It emphasizes the necessity of using methods that are unbiased and able to precisely capture extreme values of the correlation length. The former is provided by MPSbased schemes, evidently seen by the consistency of the results obtained using those methods, see Ref. 50. The latter can be provided by working directly in the thermodynamic limit, which allows to avoid problems posed by strong finite-size effects and limitations on the possible system sizes in such simulations. On the other hand, a proper extrapolation scheme of the correlation length significantly lessens the systematic limitations caused by finite bond dimension always present in MPS simulations. Both allow us to expect excellent accuracy of the results presented here, especially as we keep in mind the quality of the data obtained for the XXZ model in the previous section.

V. CORNER TRANSFER MATRIX SIMULATIONS
In the corner transfer matrix methods the infinite environment of a given site (or sites in a unit cell) in a 2d tensor network is approximated by a combination of four corner C j and four top T j tensors of finite size, as depicted in Fig. 7(i). This allows to compute expectation values of local operators of interest as well as their correlation functions. To that end, we define column-tocolumn transfer matrix as shown in Fig. 7(ii).
We employ the general corner transfer matrix algorithm described in Ref. 52, suitable for nonsymmetric problems (i.e., when corners C j are not Hermitian) and its less expensive variant described in Ref. 53. As compared to the former one, it effectively avoids squaring small Schmidt values of the enlarged corners, which are later inverted in the algorithm. Furthermore, it reduces the leading cost of the algorithm. This allows to reach larger CTM bond dimensions, which we call D.

A. Classical 2d Ising model
The 2d classical Ising model can be exactly mapped on the 1d XY model [54,55], which we discussed in details in Sec. IV A and as such we are not going to repeat those results here.
The main point that we would like to make here is that the proposed extrapolation procedure gives essentially the same result for two algorithms: the CTM method for 2d classical model and the VUMPS algorithm for the corresponding XY model. The CTM method was used in its symmetric [7,8] as well as nonsymmetric [52,53] form. We conclude that the accuracy of the extrapolation method described in this article, and more generally renormalization of the exact quantum transfer matrix in CTM/MPS simulations [19][20][21], is mostly independent of the particular algorithm used to obtain it.

B. 8-vertex model
In this section we test our approach in a numerically significantly more challenging 8-vertex model. We use the standard formulation of the model, see, e.g. Refs. 37, 56, and 57, with local two-state variables living on the edges of a 2d square lattice. Each local variable is represented by an arrow pointing toward one of the two adjacent lattice vertices. Only configurations with an even number of arrows pointing out of any vertex are allowed and contribution of each configuration to the partition function is calculated as a product of Boltzmann weights for each vertex, assigned as in Fig. 7(iii). The weights are parametrized as  a > b, c, d). We focus on the vicinity of the continuous critical point, where the critical temperature T c is set by the relation a c = b c + c c + d c . The system is in the ferroelectric phase for T < T c and the disordered phase for T > T c . The expectation value of the vertical arrow direction serves as the order parameter in this model and the correlation length associated with this observable was calculated analytically in Ref. 37.
As this problem is in general not symmetric (for c = d), simple CTM implementation [7] cannot be used [57] and it is necessary to employ the most general CTM suitable for such problems [52,53]. See Ref. 57 for some recent CTM studies of symmetric, but not exactly solvable, modification of this model. We show the results of simulations in Fig. 8 In the ferroelectric phase, T = 0.998T c , it is necessary to focus on the part of the TM spectrum that contributes to the order parameter correlation function. We define δ as in Eq. (11). The correlation length is extrapolated up to an error below 0.1% using nonlinear fit, as shown in Fig. 8(a). To that end we used TM eigenvalues for which the form factors plotted in Fig. 8(b) are in the ∼ 10 −7 band. The corresponding error when δ = 1/D is used is ∼ 5%.
It is worth noting that the TM spectrum contains another, longer scale of length obtained from the two largest TM eigenvalues 1,2 . The order parameter form factors corresponding to those eigenvalues are zero up to the numerical precision, i.e. they are not visible in Fig. 8(b) (Note that form factors are defined as a product of two numbers of similar magnitude. This means that the values larger than 10 −20 are considered to be nonzero). This length scale corresponds to another band of QTM eigenvalues, which, however, does not contribute to order parameter correlation function as the corresponding form factors vanish due to symmetries of the model [58]. In order to break this symmetry we regard the 8-vertex model as a special case of a more general 16-vertex model, where vertices with three-in-one-out and three-out-one-in arrow configurations are considered together with the ones already depicted in Fig. 7(iii); see e.g. Ref. 59-61. 1,2 contributes to the correlation functions, and have nonzero form factors, for operators build from those additional vertices. The corresponding correlation length is represented by the lower dashed line in Fig. 8(a).
In the disordered phase, T = 1.002T c , the situation is much simpler as the dominant eigenvalues have nonzero order parameter form factors. In Fig. 8(c) we recover the true correlation length with error ∼ 0.1%, both by focusing on TM eigenvalues with form factors in the dominant band, ∼ 10 −7 in Fig. 8(d), and by taking δ = 2 − 1 . The second one works well even though the form factor for 2 (marked as stars in Fig. 8(d)) is orders of magnitude smaller than for 1 (circles). It is worth observing that the QTM has complex spec-trum, and the above physical correlation lengths are shorter than could have been expected from considering only the absolute value of the largest QTM eigenvalues [37]. The longest scales effectively cancel out due to a combination of continuously changing complex phase of QTM eigenvalues, together with the proper symmetry of the form factors (periodicity in the space of parameters quantifying the spectrum). On the other hand the spectrum of the column-to-column TM is real and positive. It directly describes the physical length scales, which (if necessary) can be distinguished with the help of nonzero form factors corresponding to proper operators.

C. 2d quantum states
The PEPS corresponding to the partition function of 2d classical models -analyzed in the previous sections -are given analytically. Here, we focus on 2d quantum systems, where PEPS is obtained as a result of a suitable variational procedure. We argue that our method can be applied in such case as well. In particular, we present results for a quantum Ising model on a square lattice, where we set g = 2.5. We also consider a system of interacting spinless fermions on a honeycomb lattice, at half filling, n = 1/2. We set V = 2. Here, c i is a fermionic annihilation operator on site i and n i = c † i c i is a fermion number operator.
We focus on states at finite temperature, where their PEPS description was obtained using variational tensor network renormalization (VTNR) [62][63][64]. We show the extrapolated correlation lengths in Fig. 9. Again, a clean functional behavior of 1 as a function of the refinement parameter δ = 2 − 1 introduced here is observed. Together with all the other results presented in this article, this allows us to conclude the superiority of our method also in such a case.

VI. FORM FACTORS RENORMALIZATION
In this section we focus on the behavior of the dominant form factors. We collect the data for all the models studied in this article in Fig. 10. For each D we plot the nonzero form factors and the corresponding largest TM eigenvalue. As explained below, we plot the form factors as a function of an error of an inverse correlation length, i.e. the distance between˜ oo 1 and the actual value of the correlation length inverse, exact , for some correlation function C oo (R). In each panel, D is increased from right to left. For the models where the exact value of the correlation length is not know analytically we use the result of our extrapolation procedure.
The main observation here is that the form factors are decreasing as a power law of an error defined above. Also, the exponent of that scaling matches the value of η in the correlation function asymptotics in Eq. (2). We base this observation on the models where η is known analytically. Indeed, for the Ising model [31] in the paramagnetic phase, Fig. 10(a), η = 1/2, 3/2 and 2 for C xx (R), C yy (R) and C zz (R) correlators, respectively. In the ferromagnetic phase, Fig. 10(b) we have η = 2, 3 and 2 for the above correlators. At the critical point, Fig. 10(c), η = 1/4, 9/4 and 2, respectively. The above values are in very good agreement with the results of fits in Fig. 10(ac).
It is worth noting that for C yy (R) the leading dependence, η = 1, is obtained from the sector with complex phase ϕ = 0, see Fig. 3(a). This is in agreement with the analytical result where the leading behavior is monotonic [31]. The oscillating part of C yy (R), with frequency ϕ = ϕ XY , appears only with larger powers of 1/R in the algebraic part. Indeed, this can also be seen in the scaling of the respective form factors in Fig. 10(d).
Before we discuss other models, we provide an argument for why such a relation is to be expected. Generally, starting with the quantum transfer matrix, the correlation function can be expressed as C oo (R) = d kf oo ( k)e − ( k)R , where e − ( k) are eigenvalues of QTM, parametrized by some set of continuous parameters k. f oo ( k)d k are the corresponding form factors. We assume that the relevant contributions in some frequency of oscillations of the correlation function can be collected as e −R/ξ ymax Here the asterisk means that we integrate over k contributing to given frequency of correlation function oscillations. In order to obtain the asymptotics, C oo (R) ∼ e −R/ξ R −η , the integrated correlatorf oo (y) would have to scale as y η−1 in the limit of small y. We can view the MPS TM as an approximation to the true QTM resulting from a renormalization group procedure which captures relevant degrees of freedom for the effective impurity problem along the virtual (imaginary time) direction of a true QTM [19][20][21]. This means that the dominant eigenvalue of the TM contributing to some C oo (R),˜ oo 1 , should represent contributions from the range of the smallest y ∈ [0, y 1 ], with˜ oo 1 − exact ∼ y 1 . The corresponding form factor is then obtained by averaging over the same range,f oo 1 y1 0f oo (y)dy. Collecting those results, together with the expectation that f oo (y) ∼ y η−1 , we obtaiñ The above argument is rather qualitative and should be understood as representing general intuition of what renormalization of the (virtual) degrees of freedom of QTM looks like in the numerical procedure leading to given MPS approximation with finite bond dimension. It is, however, hard to expect that a more formal derivation is possible, because the numerical algorithm is directly targeting variational energy as a figure of merit. As a result, it is not directly related to the full TM. Nevertheless, the results of this section give strong support to the above argument.
We can further test it in the other models considered in this article. The XXZ spin chain provides an interesting example. The usual correlation function asymptotics is modified by logarithmic corrections at the isotropic critical point. It is true for spin-1 2 and spin-3 2 models. The theoretical prediction reads C xx (R) = C zz (R) ∼ x =˜ oo 10. Nonzero form factors corresponding to dominant TM eigenvalues for all the models discussed in this article. We observe that they decay (as we increase the bond dimension) as a power law with the error of the inverse correlation length,˜ oo 1 , extracted directly from the transfer matrix. We observe that the exponent of the power law coincides with the corresponding exponent η of the correlation function asymptotics in Eq. (2). In each plot we show the fitted value of the exponent, as well as the value of η when the analytical result is available. See text for details. In each panel, the order of the legend corresponds to the respective position of different lines. Recall that f oo n is the form factor (for some Coo(R)) corresponding to n.˜ oo m are the eigenvalues for which the form factor is nonzero. Those nonzero form factors are marked asf oo n .
with η = 1 [45,66,67]. This relation is supported by numerical studies of finite systems [46,68]. We further corroborate this in the Appendix A using results of our iDMRG simulations. Those logarithmic corrections, which are stronger in the spin-3 2 case, should also manifest themselves in the scaling of the form factors. Similarly as above, in this case we expect to recover the following relatioñ For XXZ spin-3 2 model at the critical point, we test it in Fig. 10(e) and obtain excellent agreement with the above prediction. The situation, however, becomes less clear in the gapped phase. For ∆ = 1.01 shown in the loglog plot, Fig. 10(e), the scaling clearly deviates from the straight line, i.e. pure power law. If we assume that the vicinity of the critical point is still influencing the scaling (at least for this range of D) and allow for logarithmic corrections also in this case, then the data become consistent with η ≈ 2, for both C xx (R) and C zz (R).
Such claims can be supported by analyzing the XXZ spin-1 2 model where the situation is similar, although the effects of logarithmic corrections are weaker. In the massive phase, ∆ = 1.3, shown in the log-log plot, Fig. 10(f), the scaling deviates again from the clear power law. If we however assume that the logarithmic corrections are still relevant here, we should include them in the fit using Eq. (24). This allows to recover the exponent η ≈ 2 for both C xx (R) and C zz (R). This is in very good agreement with the analytical results, η = 2 [35,38]. The dotted line ∼ x 2 in the plot serves as a guide for the eye.
In the Bose-Hubbard model with unit filling the situation seems to be much simpler. In Fig. 10(g) we analyze the form factors for b † 0 b R and n 0 n R and the powerlaw fits allow us to obtain the values of η very close to the ones predicted theoretically [69][70][71][72][73]. For J = 0.304, very close to the critical point, we obtain η = 2.00 and η = 0.24 respectively. The exact values are 2 and 1 4 . In the Mott insulator phase, J = 0.26 < J c , we have η = 2.01 and η = 0.47, respectively, which should be compared with the expected 2 and 1 2 . Finally, in the 8-vertex model for J = 0.2, J = J = 0.1 simulated here, we analyze the form factors corresponding to the order parameter. From the power-law fits, see Fig. 10(h), we obtain η = 3.05 and η = 3.92 for ferroelectric (T = 0.998T c ) and disordered (T = 1.002T c ) phases respectively.
To conclude this section, a few remarks are in order. We should note that the power-law fits in Fig. 10 are susceptible to the value of the inverse of the true correlation length exact , which is especially relevant when the exact value is not known. As a result, errors of up to a couple of percent might be expected here. Nevertheless, we can contrast the indirect approach presented here with fitting the asymptotics of the correlation function directly. As we illustrate in Appendix A, the latter can give substandard results away from the critical point within MPS simulation. As shown there, even in the simple Ising model in the ferromagnetic phase, the values of both ξ and η cannot be extracted with high accuracy from fitting Eq. (2) directly. The indirect method discussed in this section reduces the error by an order of magnitude.
Finally, the data shown for all the models in this section are consistent and well described by the scaling in Eqs. (23) and (24). Additionally, they are in good agreement with the analytical values of η (when available). This allows us to gain deeper understanding about the information encoded in the MPS TM and its relation with the true QTM [19][20][21]. It also allows us to expect good accuracy also for the models where the value of η is not otherwise known.

VII. CONCLUSION
The main message of this article is that some care has to be taken when extracting long-distance properties from MPS and CTM simulations. The extrapolation procedure introduced here allows to extract the correlation length with much better accuracy than the widely used method based on the bond dimension. For onedimensional systems, in the vicinity of the critical point, it can be used in parallel with the finite entanglement scaling to corroborate the results of the latter. We note that, for systems in two dimensions, one has to combine the two approaches. In that case our method can be used as a prerequisite to control CTM contraction error and set up finite correlation length scaling [11][12][13].
We test our procedure, among others, in the gapped phase of the one-dimensional Bose-Hubbard model with unit filling. This allows us to fit the scaling form of the correlation function in the vicinity of the BKT critical point and extract, among others, the position of the critical point with high accuracy.
In this appendix, we illustrate problems related to fitting the correlation function asymptotics directly, which provides further evidence of superiority of the extrapolation scheme proposed in this article. In Fig. 11(a), we show C xx (R) in the ferromagnetic phase of the Ising model. We check that the results are converged in the bond dimension. We fit the asymptotic form in Eq. (2) for different ranges of the distance R to that data.
We expect three regimes of R. For relatively short distances 1 R ξ the behavior of the correlation function is still strongly influenced by the vicinity of the critical point, for which C xx (R) ∼ R −η1 , with η 1 = 1/4 in our example. In the other extreme limit, R 1/δ, where δ =˜ xx 2 −˜ xx 1 is the measure of deviation from the continuous spectrum, we end up with purely exponential behavior, e −R˜ xx 1 , i.e. the one in Eq. (1). This is an ultimate consequence of finite bond dimension used in the numerical simulations. Finally, let us see if one can recover the exact asymptotic, e −R/ξ R −η , in the intermediate regime ξ R 1/δ. The combination of required scale separation, numerical precision and limitations imposed by finite bond dimension makes this interesting limit very hard to attain in practice, especially when η = η 1 (here, η = 2).
A smooth transition between such three limits can indeed be recognized in Fig. 11(a). In our example ξ = 49.7 . . . and 1/δ ≈ 700 (D = 512), as can be read from Fig. 2(d). It is worth observing that the largest values of locally fitted correlation length, even though they are larger than the value given by the largest TM eigenvalue, are still almost 2% away from the exact one. Additionally, they strongly depend on the window of R's used in the fit. As a result, extracting it in this way is not very reliable. For comparison, our extrapolation procedure gives a result that is an order of magnitude more accurate, see Fig. 2(d). Similarly η ≈ 1.5, which is the largest local estimate in Fig. 11(a) is far from being precise. A significantly better estimate is obtained from scaling of the form factors in Fig. 10(b). We stress here that while the discussion above is quite general, the illustrative example of the Ising model is well known not to be challenging for MPS-based methods. As such, we could expect even more severe problems for more demanding systems.
Paradoxically, the situation at the critical point is much simpler, even though such points are generally harder to simulate with MPS. Above all, there is no physical length scale ξ in this case. As such, we can expect to recover the exact asymptotics for 1 R 1/δ, with finite D effects becoming relevant on larger distances only. This, in principle, makes direct extraction of the correlation function asymptotics much more straightforward in this case (there might still be problems related, e.g., with spontaneous symmetry breaking resulting from finite D, or other effects of inefficient description of critical points with MPS).
As an illustrative example, in Fig. 11(b) we show C xx (R) for the isotropic antiferromagnetic Heisenberg spin- 3 2 chain. The difference between the results for two bond dimensions shown there, as well as difference from C zz (R), are well below 10 −5 . In this case the scale 1/δ ≈ 3000, as can be read from Fig. 5(a). The refinement parameter δ introduced in this article is not suited for extrapolation of local quantities, such as energy. This is illustrated using the example of (a) Ising model and (b) XXZ spin-1 2 model, where we show the error of ground state energy ∆E converged for different D. In panel (a) D = 8-100 and in (b) D = 32-1448. This implies that the methods which do work well for the extrapolation of the energy would not work for the correlation length.
The above example, apart from corroborating older results, clearly shows that extracting the correlation function asymptotics from MPS simulations at the critical point is a viable method of obtaining the exponent η, unlike for the system away from criticality. Indeed, even logarithmic correction is clearly recovered in the example and has to be taken into account in the fit. Likewise, we should note that obtaining η in the critical systems from the fits to the form factors, as in Sec VI and Fig. 10(c,e,f) seems to yield comparable precision of the results. The latter method, however, proves to be superior away from the critical point.