QCD at finite chemical potential
with six time slices
Abstract
We investigate the Taylor expansion of the baryon number susceptibility, and hence, pressure, in a series in the baryon chemical potential () through a lattice simulation with light dynamical staggered quarks at a finer lattice cutoff . We determine the QCD cross over coupling at . We find the radius of convergence of the series at various temeperatures, and bound the location of the QCD critical point to be and . We also investigate the extrapolation of various susceptibilities and linkages to finite chemical potential.
pacs:
12.38.Aw, 11.15.Ha, 05.70.FhI Introduction
We report physics obtained in QCD with two light flavours of dynamical staggered quarks at finite temperature, , and chemical potential, , at lattice spacing . We investigate the physics at finite chemical potential cep using the method of Taylor expansions that was developed in pressure and used for QCD earlier in nt4 ; biswa . One of the quantities we investigate is the radius of convergence of the Taylor series, through which we estimate the QCD critical point. We also investigate the dependence on of various other quantities of physical interest. Finally, we investigate the linkage of quantum numbers and its dependence on and . Related earlier works are others .
That phase transitions are rounded off by finite size effects was discovered long back by van Hove. The most familiar aspects are seen when simulations are directly performed in the vicinity of the critical coupling. Quantities which would diverge in the thermodynamic (infinite volume) limit are finite. As a result, a lattice computation never sees a singularity, but infers its existence from some measures. Proofs of the existence usually involve testing extrapolations: such are the main remaining problems at finite temperature and vanishing chemical potential. A welldeveloped finite size scaling theory can be used to study the size, , dependence of such quantities and extract critical exponents. To date, the immensity of computational requirements has prevented full use of this theory for QCD.
The study of the effect of finite size rounding of critical points on series expansions is, to the contrary, in its infancy. The clearest fact about such effects is the following: since quantities which should diverge in the thermodynamic limit merely have finite peaks at finite , series expansions, strictly speaking, have finite radii of convergence for finite . In the limit the radius of convergence estimates the nearest critical point. The mechanism by which this limit is reached is straightforward. Were one to examine some estimator of the radius of convergence at order , , one would find that up to some the would approach a finite value. Such behaviour was exhibited for the series expansion in QCD in nt4 . For larger the would diverge, in accordance with van Hove’s theorem. With increasing one would find that approaches infinity. How the scaling of with codes for the critical exponent is at present unknown.
Another question that can be answered by the series expansion is where the singularity lies. In the case of QCD, the question is whether the finite radius of convergence of the series expansion in is due to a singularity at real . If it is, then the successive coefficients in the expansion are positive. At finite one must examine the signs of the series coefficients for . If these are positive, then the singularity which limits the expansion when is on the real axis. In nt4 this argument was used to justify the identification of the radius of convergence with the critical point. This argument is also implicit in misha06 .
In summary, the quark number susceptibilities are Taylor coefficients in the expansion of the pressure in powers of the chemical potential. From the series expansion for the pressure,
(1) 
we define the nonlinear susceptibilities (NLS) of the th order, . The second order susceptibility, also called the quark number susceptibility (QNS) gottlieb , has the expansion
(2) 
This series is expected to diverge at the QCD critical end point. We define the radius of convergence of this series as
(3) 
When successive estimators for are equal within statistical errors to the same value , we have identified the plateau in the radius of convergence. This corresponds to the critical point, provided the singularity in the series occurs at a real value of . In turn, this is the case when the coefficients from which the estimates are made are all positive.
In the next section we present the details of the simulation and the extraction of the critical coupling. This is followed by a section in which we report the main results, namely, the extraction of the QNS up to the eighth order. This results in five terms of the series for the pressure, and four terms of the series for the baryon number susceptibility. Using these we report our result for the radius of convergence of the series, and extract from this our best estimate of the critical point. In the section after this we discuss the extrapolation of physical quantities to large chemical potentials. This extrapolation throws more light on the nature of the QCD critical point.
Note that the series in eqs. (1,2) cannot be continued beyond even when all the terms are known exactly. The truncated series expansion fails even faster. As a result, it becomes difficult to extrapolate physical quantities to large values of . One way to use the series to better advantage is well known: the method of Padé approximants. The existing theory of Padé approximants baker is adapted to the case where each known series coefficient has infinite numerical accuracy. When coefficients are extracted through Monte Carlo estimates, and hence have statistical errors, new issues arise. We believe that it would be useful to extend the theory of Padé approximants in this direction. In the appendix we make a beginning which is adequate for the purpose of this paper.
Ii Simulations
5.39  2284  33  
5.40  22599  88  10099  48  919  35  
5.41  50584  197  14580  131  
5.415  17518  179  14044  158  
5.42  39649  164  35649  165  27974  140  
5.425  50589  189  47329  214  53563  267  
5.43  54619  218  41349  147  41869  202  
5.46  309  13  10719  13  1214  13  
5.54  969  7  
5.60  2891  4  
5.75  3626  4 
Operator  

1.611 (1)  1.611 (2)  
1.611 (1)  1.611 (2)  
0.031 (3)  0.295 (5)  
0.293 (9)  0.291 (9) 
The simulations were performed using the Ralgorithm for hybrid molecular dynamics. This uses a finite step size, , for the molecular dynamics. Our main data set is generated using and a total trajectory length of in MD time units. We performed tests of the accuracy and efficiency of the choices.
Most of our computations were performed with . This was found to be adequate for computations at . We checked our results at by running a long computation with and trajectory length of 3 MD time units. We found complete agreement between the runs with two different time steps. In Table 2 we show the comparison of bulk quantities computed in the two runs.
Changing the trajectory length from to at , and did not change the results for thermodynamic quantities within errors. However, near the longer trajectories were more effective at reducing the autocorrelation time. For example, we found that the longest autocorrelation time at was trajectories for , and it reduced to 36 for . As a result, the CPU time taken to produce decorrelated configurations is reduced by a factor of about 2.5 on taking the longer trajectories. At the effective speedup, computed in the same way, was a little under a factor of two. In the high temperature phase the autocorrelation times were very small, and there was little to be gained by using longer molecular dynamics trajectories. There were no changes in thermodynamic quantities on changing the trajectory length.
A range of was scanned to locate the bare coupling at the finite temperature crossover, . The crossover was located by the peak of the unrenormalized Polyakov loop susceptibility,
(4) 
Here is the gauge link in the time direction at the spatial site on the time slice . We shall show later that the fourth order quark number susceptibility also peaks at the same coupling. This is closely related to the inflection point of the second order susceptibility which is used by various groups hotqcd .
From the peak of we identify , where the uncertainty is due to resolution, and not a statistical uncertainty. There is a little finite volume shift in the position of the peak of at the smallest volume, but no such shift is observed in going from to (see Figure 2). While some volume dependence is visible in the peak of , with data from just these three volumes it is not possible to decide whether there is a crossover or a critical point at . However, it is possible to rule out a first order transition, since definitely drops with increasing , as shown in Figure 2.
Subsequently, runs were performed on lattices of size and on a grid of to determine the scale. The scale determination used the value of the plaquette to obtain the renormalized gauge coupling in the scheme. The errors in the scale setting involve the uncertainty in the location of the crossover coupling, , statistical errors in plaquette measurements, and scheme dependence estimated by evaluating the scale also in the E and V schemes. In the range of temperatures within 20% of , the largest errors came from the uncertainty in the determination of . Better results can only be obtained by using larger spatial volumes. At larger temperatures, the scheme dependence of the scale set the largest errors. These can be reduced by going to smaller lattice spacing, i.e., to larger . The scale setting using the crossover for is compatible within errors with that obtained earlier for a similar setting of scales for .
Iii Quark number susceptibilities
A quick reminder of our notation pressure ; nt4 is in order. A quark number susceptibility is obtained by taking a derivative of the pressure with respect to the chemical potential. In two flavour QCD there are two possible chemical potentials, and . If one takes derivatives with respect to and with , then the order of the quark number susceptibility is . Since the and quarks are degenerate, and indistinguishable at , we denote the susceptibilities by when and when . The susceptibilities are constructed from expectation values of a string of operators sandwiched between quark propagators. The operator is the operator with insertions of into a single fermion loop, and hence contributes only to . The operators are products , and may contribute to several of the . The construction of on the lattice is given in detail in nt4 . We shall discuss results for as well as the expectation values (since we discuss only the connected pieces of these expectation values, we have not used separate notation for that).
The quark number susceptibilities are obtained as expectation values of fermion loops with various operator insertions nt4 . These are evaluated as usual through stochastic estimators. The computations were optimized using the methods of nt4 . The need to use large number of stochastic vectors has been discussed in detail elsewhere sewm . We have taken 500 random vectors for each trace evaluation. With this we are able to control statistical errors on loops with up to six operator insertions. Even so, loops with larger numbers of insertions remain noisy. Thus, at on the largest volume, gives a signal at 53 and at 23, whereas and give signals at and respectively. For the two highest susceptibilities, this level of the signal is an improvement over the corresponding results with at equal . It was our experience at coarser lattice spacing that one needs lattices with larger spatial volumes to control loops with more insertions. We see this also at the current lattice spacings, at the smallest volumes, even loops with six insertions are hard to control.
In the following sections we will often compare results for and . These results are meaningful only if they are done holding other factors fixed. We shall therefore compare the new results obtained on lattices with those obtained earlier on with the same quark mass, .
iii.1 Second order
The lowest order quark number susceptibilities are shown in Figure 3. The diagonal susceptibility, seems to show significant dependence on , i.e., the lattice spacing. This is not a surprise; after all, even in the quenched theory a similar effect was seen quenched . The offdiagonal susceptibility, seems to scale better with the lattice spacing.
Note that takes contributions only from , whereas has contributions from this as well as . The results shown in Figure 3 indicate that the quarkline disconnected operator has, at best, marginal lattice spacing dependence. Most of the lattice spacing dependence seen in therefore comes from . This last expectation value is the response to a chemical potential on the isospin component and hence was called in some of our early papers. Both this and change rapidly near and the “inflection point”, i.e., the point at which the slope is maximum can be used as a corroborative measure of . Because of the numerical uncertainties in taking derivatives of noisy data, we will instead use the peak in the fourth order susceptibility. We discuss these next.
iii.2 Fourth order
Two of the fourthorder susceptibilities are shown in Figure 4. Both and peak at . This was already seen in earlier simulations with . Within the resolution of our measurements, we see that the peak in these quantities comes at exactly the same coupling as the peak in , at both and 6. Like the second order susceptibilities, these too have significant cutoff dependence. is much smaller than either of these susceptibilities and shows no special structure near .
The peak at can be resolved into a single operator expectation value, . This expectation value peaks at and falls off rapidly on both sides, as shown in Figure 5. Therefore it can serve as a good measure of the critical coupling. The expectation value , on the other hand shows a crossover near . One could construct yet another measure of the critical coupling from the point of steepest slope of this expectation value, or from its variance, the expectation value . This last quantity contributes to eighth order susceptibilities.
Using the fourth and second order quark number susceptibilities, one can form the first two terms of the series expansion of nt4 . From these coefficients can obtains the lowest order estimate of the radius of convergence of this series, . This is shown as a function of in Figure 6. Note that the large dependence on the lattice spacings seen in each of the susceptibilities almost cancel out in the estimate of the radius of convergence. The radius of convergence has smaller dependence on the lattice spacing.
iii.3 Sixth order
The sixth order NLS is shown in Figure 7. It has been pointed out earlier that has the form of a rounded step function, and that successively higher order NLS have the form of rounded derivatives of the step function structure . For example, the fourth order NLS has a peak. The sixth order NLS changes sign near and has a maximum and a minimum flanking the zero. This behaviour is clearly visible in Figure 7. This peculiar structure comes from the behaviour of the expectation value , also shown in Figure 7. Note that the measurement of is noisier than that of .
The quarkline connected operator expectation value at this order is . This is shown in Figure 8. Note that this has interesting structure below and that the structure is seen for both and 6. As a result, one cannot use or its variance for determining .
iii.4 Eighth order
The eighth order NLS are fairly noisy below and in the vicinity of . This is partly due to the fact that operators with multiple quark loops, such as , become noisier as the number of loops increases. However, singleloop operators such as are also not under sufficient control at these lattice volumes. We exhibit the expectation values and in Figure 9. Note that there could be structure in the former below , but this is currently obscured by noise. The latter has a single sharp peak at , as argued before, and shows that identified by the peaks in , , and are identical within the resolution of our study.
iii.5 Radius of convergence
The radius of convergence of the series expansion can be used to estimate the position of the critical end point of QCD as before. The radius of convergence gives the distance from the origin where the expansion ceases to hold. The corresponding singularity lies at real if the coefficients of the series expansion are all positive. As one comes down in from large values of , the series coefficients remain positive down to some lowest value of . At this temperature the radius of convergence is independent of the order at which it is evaluated and has value . Thus, our estimate of the position of the critical end point is
(5) 
with a lattice spacing of and a renormalized quark mass that corresponds to tuning , when the spatial size of the box is . In comparison, with a lattice spacing of at the same renormalized quark mass and the same spatial volume, it was found that remained unchanged whereas one had . Extrapolation of this result to the thermodynamic limit, on the coarse lattice yielded an estimate , which, although statistically compatible with the finite volume result, had a lower mean. It would be interesting to check how much the new estimate of the QCD critical end point is lowered upon taking the thermodynamic limit. However, this extrapolation lies outside the scope of the current work because of the CPU resources needed.
Iv Physics at finite
iv.1 Fluctuations and the quark number susceptibility
Baryon number fluctuations, by an amount from the expectation value in a grand canonical ensemble, have a spectrum
(6) 
It has been proposed that the susceptibilities be measured in eventtoevent fluctuations in heavyion collisions fluct . Indeed, the divergence of the width of the spectrum of fluctuations could be one signal for the detection of the QCD critical point in experiments critexpt . In view of this, it is interesting to estimate as a function of along the critical isotherm.
While the truncated series expansion can be used to estimate the radius of convergence, it cannot be used to extrapolate the susceptibility up to that point. As shown in Figure 10, the series expansion for taken to orders and fail to agree long before the radius of convergence is reached; nor do they show any divergence at . In order to extrapolate the QNS, one must therefore find more robust techniques.
The usual method is to convert the series to a Padé approximant (see mariapade for a previous application to QCD). There is extensive literature on the use of these methods when the series expansion is known exactly. In Appendix A we extend this theory to the case relevant to our study, i.e., when the series coefficients are known only through a Monte Carlo procedure, and hence are known with certain errors. The appendix examines error propagation in the Padé approximants, and sets out the basic methods to control these errors. For our purposes we use the Padé approximants labeled in the notation in Appendix A.
The Padé approximants and are shown in Figure 10. It is interesting to note that they diverge as is approached. While they disagree with the series expansions as the radius of convergence is approached, they remain consistent with each other except very close to the divergence. Note that the errors are large near the divergence. This seems unavoidable, since any error in the coefficients will be magnified near the pole. There are also large errors beyond the pole. It should be possible to control these in future work.
Note that in two flavour QCD one has , at all , as long as the isospin chemical potential remains zero. So there are two independent susceptibilities, and . In terms of the previously computed quantities, they are linkage ,
(7) 
It turns out that remains small within errors even at larger chemical potential, so that the behaviour exhibited in Figure 10 for is also almost quantitatively correct for after an overall normalization by a factor of . In particular, the divergence in is also reflected in . This has consequences which we deal with next.
iv.2 Linkage
Earlier works have introduced quantities which measure whether two quantum numbers vary together in thermodynamic fluctuations koch ; linkage . The most straightforward measure, called the linkage, utilizes diagonal and offdiagonal QNS in the form of the ratio
(8) 
for any two quantum numbers and . The linkage gives the thermal averaged amount of the quantum number excited per unit in a thermal fluctuation taking place in the grand canonical ensemble. In twoflavour QCD one may measure the linkage between U and D quantum numbers (conventionally +1 for quarks, 1 for antiquarks of the correct flavour, and zero otherwise). Also related is the linkage between the baryon number, , and the electrical charge, .
In Figure 11 we show the temperature dependence of . At this quantity should be , since the lightest excitation is a pion, and the two charged pions each give a contribution of , whereas the neutral pion gives a contribution of 0. In the hightemperature phase, when the lightest excitations are quark quasiparticles, the linkage should vanish. We see a rapid cross over between these two regimes, with a very small but nonzero value being reached at . We also exhibit the linkage . At this quantity is expected to vanish, since the lightest charged particle, the pion, has no baryon number. In an ideal quark gas, this linkage has value . One sees a rapid crossover between these two values in the vicinity of , exactly as for .
In the chiral limit, i.e., when the quark masses vanish, and a second order phase transition to occur, one would expect that becomes infinitely peaked at . As a result, one expects the diagonal susceptibilities to become infinitely sharp, and the linkages to jump abruptly across the transition. Some part of the rounding in the linkages is therefore due to the fact that the quark masses are finite. However, the rounding of the crossover in the linkages would be a direct demonstration that there is no abrupt change from the hadronic to the quark phase: one may use either description over a small range of temperatures near . This could have implications for the description of hadronization in a fireball, a process which, at the moment, has a very crude description in terms of the CooperFrye mechanism cooperfrye . However, a part of the rounding is also due to finite volume effects, and it is hard to disentangle the two in our computation. It would be an interesting future computation to understand quantitatively what part of this slow crossover is a finite volume effect and how much is the effect of a finite quark mass.
Since we have control over the higher order NLS, we can construct a Taylor series expansion for the linkage and examine its behaviour at finite chemical potential. For the analytic continuation of the linkage, we perform a jackknife analysis. In each jackknife bin the Padé approximant is evaluated at the chemical potential of interest. The mean of these values is used as the estimator for the continuation, and the 68% interval, evaluated nonparametrically, is quoted as the error bound. The results are shown in Figure 12 for several different temperatures.
At the highest temperature, i.e., , the linkage is close to zero at and remains zero for . At temperatures below , the linkage is nonzero at but evolves smoothly with the chemical potential. For we find a smooth increase with , the linkage decreasing with larger . This illustrates the important point that a finite radius of convergence for one susceptibility does not imply divergences in other quantities.
Interestingly, the linkage seems to fall marginally with increasing . This mild effect can be traced to the fact that the fourth order coefficient in the Taylor expansion of this linkage is small and negative. This results in a fall at large . It would be useful to check whether this persists at larger volumes, and whether higher order terms in the expansion turn this around and cause the linkage to increase. An interesting alternative possibility is that the fall in is physical, as is the rise in , and the two together imply the existence of a phase analogous to the quarkyonic phase at large largen . It is therefore of interest to check these results further. Unfortunately, both checks require massive investment of CPU resources, but are interesting enough that we hope to return to this soon.
V Summary
We have examined QCD with two flavours of dynamical staggered quarks at finite temperature with lattice spacing and bare quark masses tuned to . This quark mass is expected to correspond to , and hence our new results are directly comparable to the older results which were obtained on a coarser lattice with nt4 . Our simulations were performed on lattices with size , where is the spatial size of the box. We used the Ralgorithm with a step size in MD time units of . We have checked that decreasing this by an order of magnitude to does not change thermodynamic results (see Table 2). Similarly, we have checked that the physics results remain unchanged when the trajectory length is changed.
We identified the cross over at vanishing chemical potential through the Polyakov loop susceptibility, , (see Figure 2) and then cross checked this through two measures related to the QNS. One is the peaking of and (see Figure 4) which is related to the “inflection point” of the QNS. The other is the peaking of the operator , (see Figure 9) which is related to a similar inflection point in . These measures are consistent with each other within the accuracy of our computations. The scale setting using this identification of the cross over is consistent with the earlier scale setting using coarser lattice spacing nt4 .
We presented results for the NLS up to eighth order. There are clear lattice spacing effects, as expected. These are roughly consistent with earlier determinations of some of these quantities in quenched theory. While the lattice spacing artifacts for the NLS are very large, sometimes as much as 100%, the effect on the radius of convergence is much smaller (see Figure 6). Our estimate of the critical point of QCD is based on this radius of convergence. The critical point occurs when the radius of convergence identifies a singularity on the real axis, through the fact that the series coefficients are all positive. Caveats on this are presented in the introduction. Our estimate of the critical point using finite volume data is (see eq. 5)
This should be compared with our earlier estimate on the same lattice volume and same renormalized quark mass which gave . This is a change of about 26%, and is statistically significant. Extension of our results to larger volumes is outside the scope of this work. In simulations with , the estimate of dropped by about 16% on extrapolating to infinite volume.
The series expansion is a good tool for extracting the radius of convergence, and, through it, the critical point. However, as we show in Figure 10, it is a bad tool to extrapolate physical results to high . Padé approximants adjusted to give the same series expansion seem to perform better, even after taking into account the propagation of statistical errors. One sees the divergence of the susceptibility at the critical end point, something that the series expansion misses altogether.
We also examined the linkages and (see Figure 11). At vanishing chemical potential they show a rapid cross over from the values expected in the hadronic phase to those expected for the nearly ideal quarks. The rounding of this transition is closely related to the question of how sharp the hadronization transition can be in heavyion collisions. A discussion of the issue was presented in Section IV.2.
The measurements of the linkages were extended to finite chemical potential using Padé approximants. Unlike the QNS, they evolve smoothly through the critical point. Interestingly, on isotherms below , with increasing , the linkage between U and D quantum numbers changes towards the ideal quark gas, whereas the linkage between B and Q changes away from the quark gas. This could indicate the presence of a quarkyonic phase of QCD matter, although technicalities need to be sorted out before one can establish this.
The computations were performed over the last two years on the Cray X1 of the Indian Lattice Gauge Theory Initiative (ILGTI) at TIFR. We thank Ajay Salve for singlehandedly taking care of the machine during this extended period.
Appendix A Padé approximants
We follow Baker’s definition baker of a Padé approximant pade ; gaunt . The series expansion
(9) 
evaluated to order can be used to define the Padé approximant of order ,
(10) 
where and are polynomials in or order up to and respectively. From the matching condition it follows that . Introduce the notation—
(11) 
Then, writing out the matching condition order by order, one obtains the Padé approximants by solving first for the denominator
(12) 
(with the convention that for ) and then for the numerator
(13)  
The practical importance of Padé approximants arises from the fact that if the series has a radius of convergence as , then the series expansion is reliable only for , whereas the Padé approximants can be used for analytic continuation beyond this. Much of the standard theory of Padé approximants deals with the cases when the coefficient matrix in eq. (12) has vanishing determinant, and the information which can then be extracted.
Here we concentrate on a different problem— that of controlling errors when the series coefficients are obtained by a Monte Carlo program, and hence have a given statistical distribution. We found no discussion of this in the literature, although it is likely that sporadic attempts to answer related questions have been made in the past. These questions become important now that new developments in QCD at finite chemical potential lead us to analyze series coefficients obtained in a Monte Carlo process.
When the Padé coefficients are welldefined, the joint probability distribution of the series coefficients can be transformed into that of the coefficients of the Padé approximant using the usual Jacobian formula—
(14) 
Take the example of , where and , so that . Assume that and are drawn from independent Gaussian distributions of unit mean—
(15) 
Then the joint distribution of the Padé coefficients can be written down. The marginal distribution of the Padé coefficient , being the ratio of two Gaussian distributed numbers, is well known gaussrat and given by
(16)  
Note that , and hence is a finite number which depends only on . As a result, the marginal distribution of is not exponentially damped at infinity. The powerlaw damping comes from the factor . Clearly this distribution has a welldefined expectation value for , but the variance and higher cumulants do not exist. Thus, statistical measurements of are not subject to the central limit theorem.
A similar phenomenon occurs with any . Assume that the series coefficients are statistically independent and drawn from a Gaussian of unit mean, then the probability distribution of the Padé coefficients can be written as
(17) 
where is the Jacobian of the transformation given by for and , and is a quadratic form obtained by transforming the arguments of the Gaussians.
The Jacobian of this transformation is
(18) 
where is the derivative of with respect to . Multiplying the second row from the bottom by and subtracting that from the last row gives
(19) 
where we have used the relation to write .
The quadratic form in the argument of the exponent can be manipulated into a particularly useful form by completing the squares—
(20) 
where the real symmetric matrix and the vector can be easily written down. We do this next for the special case when all the are equal to .
Define the sequence of polynomials
(21) 
In terms of these, one writes
(22) 
In order to find the determinant of , we do the following row operations— starting from the top, add to each row times the next row. This reduces the determinant to a lower triangular form
(23) 
The solution of the equation can be obtained by the same operations. They yield the reduced equation
(24) 
where we have used the relation , to reduce the vector . This gives—
(25) 
Finally,
(26) 
Clearly, remains finite in the limit , so that does not damp the marginal distribution of . In fact, that damping comes from the factor of . As a result, has well defined mean but its variance is undetermined. Thus, estimators of evade the central limit theorem.
Nevertheless, the situation is pretty well under control. The appropriate question to ask of a distribution such as that in eq. (A) in the context of parameter estimation is not the value of the variance, but an appropriate measure of the variation in the estimate. One could quote the width at half maximum, or the limits such that 68% of the probability lies within these limits. Along with this one asks, if we make measurements of then how does such a measure of variation change with .
A numerical investigation shows that when , the modal value is . Since the distribution is skew, the modal value and the mean are different. The full width at half maximum is contained in , and this range contains 56.1% of the integral. The 68% probability interval is . Either of these ranges can be quoted as an estimate of the error in the modal value.
To answer the question about the distribution of means, we use the characteristic function. If is the distribution of , then the Fourier transform is called the characteristic function. Since is nonnegative and integrable, being a probability distribution, it is also square integrable, so that the characteristic function exists. The characteristic function of the mean of numbers, , is
(27) 
Fourier transforming this gives the distribution of . While this general method remains valid for the distributions , above, it does not seem possible to perform the Fourier transformations in closed form. So, instead of writing down an impenetrable formula for the distribution of means, we investigate useful subsidiary questions.
Define the skew of a distribution by
(28) 
where denotes the modal (most probable) value, and is the mean. The skew is nonzero for every skewed distribution, being positive if the distribution is skewed to the left and negative otherwise. For the distribution , we found . For the distribution of means of values, decreases. A Monte Carlo estimate for indicates that . This estimate was obtained using values of between 1 and 100. A similar result was obtained for a measure of skewness that compares the median and the mean in a manner analogous to eq. (28).
At the median of a distribution, the cumulative distribution becomes equal to 0.5. The errors on the median can be defined as the points at which the cumulative distribution is either 0.34 above or below. The distribution of the means of numbers narrows rapidly, and we find in a Monte Carlo estimate that both these intervals decrease as .
In conclusion, for the estimation of and confidence limits on the estimate, it matters little that the central limit theorem does not hold. The mean is well defined, and its difference with the mode and median scale with a factor of . The 68% confidence limits also scale as . We therefore quote the mean and 68% confidence limits on it as estimators for the Padé coefficients. These estimators are easy to incorporate into jackknife and bootstrap analyses. In parts of our analysis the estimators of the series coefficients are also not Gaussian distributed; even so, the nonparametric statistical analysis outlined here suffices.
References

(1)
A. Barducci et al., Phys. Lett. , B 231 (1989) 463;
M. A. Halasz et al., Phys. Rev. , D 58 (1998) 096007;
J. Berges and K. Rajagopal, Nucl. Phys. , B 538 (1999) 215.  (2) R. V. Gavai and S. Gupta, Phys. Rev. , D 68 (2003) 034506.
 (3) R. V. Gavai and S. Gupta, Phys. Rev. , D 71 (2005) 114014.
 (4) C. R. Allton et al., Phys. Rev. , D 71 (2005) 054508.

(5)
Z. Fodor and S. Katz, J. H. E. P. , 0203 (2002) 014;
C. R. Allton et al., Phys. Rev. , D 66 (2002) 074507;
M.P. Lombardo and M. d’Elia, Phys. Rev. , D 67 (2003) 014505;
Ph. de Forcrand and O. Philipsen, Nucl. Phys. , B 642 (2002) 290;
C. R. Allton, et al., Phys. Rev. , D 68 (2003) 014507;
Z. Fodor and S. Katz, J. H. E. P. , 0404 (2004) 050;
Ph. de Forcrand and O. Philipsen, J. H. E. P. 0701 (2007) 077;
C. Bernard et al., Phys. Rev. , D 77 (2008) 014503.  (6) M. A. Stephanov, Phys. Rev. , D 73 (2006) 094508.
 (7) S. A. Gottlieb et al., Phys. Rev. Lett. , 59 (1987) 2247.
 (8) G. A. Baker and P. GravesMorris, Encyclopedia of Mathematics: Padé Approximants, Vol 13, Part 1, AddisonWesley Publishing Company, Reading, Massachusetts, (1981).

(9)
Y. Aoki et al., Phys. Lett. , B 643 (2006) 46;
M. Cheng et al., Phys. Rev. , D 74 (2006) 054507;
C. Detar et al., PoS LAT2007 (2007) 179.  (10) R. V. Gavai and S. Gupta, Nucl. Phys. , A 785 (2007) 18.
 (11) R. V. Gavai and S. Gupta, Phys. Rev. , D 67 (2003) 034501.
 (12) R. V. Gavai and S. Gupta, Phys. Rev. , D 72 (2005) 054006.

(13)
M. Asakawa, U. Heinz and B. Muller, Phys. Rev. Lett. , 85 (2000) 2072;
S. Jeon and V. Koch, Phys. Rev. Lett. , 85 (2000) 2076. 
(14)
M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. , D 60 (1999) 114028;
M. A. Stephanov, K. Rajagopal and E. Shuryak, Phys. Rev. Lett. , 81 (1998) 4816;
G. S. F. Stephans, nuclex/0607030;
P. Sorensen (STAR Collaboration), nuclex/0701028.  (15) M. P. Lombardo, heplat/0509181.
 (16) V. Koch, A. Majumder and J. Randrup, Phys. Rev. Lett. , 95 (2005) 182301.
 (17) R. V. Gavai and S. Gupta, Phys. Rev. , D 73 (2006) 014004.
 (18) F. Cooper and G. Frye, Phys. Rev. , D 10 (1974) 186.
 (19) L. McLerran and R. D. Pisarski, Nucl. Phys. , A 796 (2007) 83.
 (20) H. Pade, Ann. de l’Ecole Norm. Sup. Serie, 9, Suppl. (1892) 3.
 (21) D. S. Gaunt and A. J. Guttmann, p. 181, Phase Transitions and Critical Phenomena, Vol. 3, eds. C. Domb and M. S. Green, Academic Press, London, 1974.
 (22) R. C. Geary, J. Royal Stat. Soc., 93 (1930) 442.