{'fraction_non_alphanumeric': 0.040707772470690406, 'fraction_numerical': 0.018671298306556665, 'mean_word_length': 4.246298405466971, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical Approximation to the Thermodynamic Integrals

We approximate boson thermodynamic integrals as polynomials in two variables chosen to give the correct limiting expansions and to smoothly interpolate into other regimes. With 10 free parameters, an accuracy of better than 0.009\% is achieved for the pressure, internal energy density and the number density. We also revisit the fermion case, originally addressed by Eggleton, Faulkner and Flannery (1973), and substantially improve the accuracy of their fits.


INTRODUCTION
The thermodynamic functions of fermions and bosons are expressed as integrals which cannot be evaluated analytically except in limiting cases of degeneracy or relativity. Numerical integration is not entirely straightforward, however, since the form of the integrand varies widely under different conditions. Thus the numerical method chosen must be carefully tested and optimized for the parameters at hand and different methods may be needed in different regimes. Furthermore, the computational time required can be an issue if the integrals need to be evaluated many times. For fermions, these problems can be overcome by using an ingenious numerical approach described by Eggleton, Faulkner, & Flannery (1973), hereafter referred to as EFF. They approximate the integrals by means of a polynomial in two variables. The polynomial form is chosen so that it yields the correct behavior in the four limiting situations for which the series can be derived analytically.
The coefficients are then optimized to accurately interpolate into regimes which are inaccessible to analytical analysis. We have found that a simple modification of the approach of EFF, which involves a single additional parameter, can significantly improve the accuracy for a given number of coefficients and, in addition, permits the approximation to merge exactly into the correct asymptotic behavior. This modification is superior to the one discussed by Pols et al. (1995). We discuss the fermion case briefly in §3.
The main purpose of the present paper is to analyze the thermodynamic integrals for bosons; any condensate contribution would, of course, have to be handled separately. While we follow the general approach of EFF, bosons require a different approximation scheme. This is discussed in §4, in which we also present results for various orders of approximation.
We begin by defining the integrals at hand in §2. The fermion and boson cases are discussed in §3 and 4, respectively, and brief concluding remarks are given in §5.

DEFINITIONS
For a relativistic system of particles, the standard form for the pressure is where g is a degeneracy factor, the energy E k = √ k 2 c 2 + m 2 c 4 , µ is the chemical potential and the upper (lower) sign refers to fermions (bosons).
The units employed are such that Boltzmann's constant k B = 1, so that the temperature is measured in units of energy. The angular integration is trivial for an infinite system and it is convenient to perform the standard integration by parts. Expressing the result in terms of the dimensionless variables of EFF and a dimensionless integration variable l = k/mc, we write the pressure in dimensionless form where E = √ l 2 + 1 − 1 and d = g 2π 2 mc 2 mc h 3 carries the dimensions of P . For fermions the value of ψ, or equivalently µ, is unrestricted. However for bosons, since the occupation probability has to be a positive quantity, ψt = µ mc 2 − 1 ≤ 0; thus ψ is always negative. For a gas containing antibosons (b) in thermal equilibrium with bosons (b), so that µb = −µ b , the requirement ψbt ≤ 0 yields the further restriction ψ b t ≡ ψt ≥ −2.
The remaining thermodynamic functions are easily obtained. We write the entropy density S in dimensionless form as The number density N written in dimensionless form is The energy density E written in dimensionless form is For e, the leading power series behavior of ρ obviously differs from that of ts, p and ρψt. We therefore work with the internal energy density 3 APPROXIMATION TO FERMION

Formalism
The approach we use follows that of EFF. We determine an approximation for the pressure in terms of the degeneracy parameter ψ and the temperature parameter t. Then, by eqs. (4-6), the remaining thermodynamic variables can be found by differentiation. The EFF scheme requires knowledge of the expansions of the pressure in the limiting cases. In the non-degenerate case (ND), ψ ≪ −1, the denominator in equation (3) can be expanded in powers of exp(ψ + 1/t) (Chandrasekhar 1958), leading to where, as before, the upper sign refers to fermions and the lower to bosons.
Here K 2 is a modified Bessel function which may be expanded (Abramowitz & Stegun 1965) in the extremely relativistic (ER, t ≫ 1) and non-relativistic (NR, t ≪ 1) limits. In the extremely degenerate case (ED) for fermions, in which ψ ≫ 1, the Sommerfeld expansion (Chandrasekhar 1958) may be employed and the extremely relativistic (ER, ψt ≫ 1), and non-relativistic (NR, ψt ≪ 1) limits obtained. Then the four limiting cases for the fermion pressure are with the corresponding limiting cases for the number density and internal energies easily obtained from equation (9) by differentiation.
The key to the EFF scheme is to find functions f (ψ) and g(ψ, t) such that equation (9) can be expressed as where a mn , b mn , c mn and d mn are coefficients. This is possible provided that ( which will be used to interpolate into all regions of ψ and t.
We choose for the derivative which differs from EFF by leaving a as a free parameter to be fitted, rather than taking it to be unity as they did. This simple change gives a significant improvement in the fit, because it allows greater freedom in matching the limiting expression for the entropy in the extremely degenerate case. In the EFF scheme the relative errors in the entropy in the degenerate limit remain substantial even when the order of the approximation is very large.
Integration yields For the non-degenerate case, f ≃ 4ae ψ−2 , which is small; while for extreme degeneracy, f ≃ 1 4 aψ 2 , which is large. We define g = t √ 1 + f , as in EFF, since we have found no advantage in introducing additional parameters here.
The leading terms of equation (9) yield values for the four "corner" coefficients and these are collected in Table 1. In addition, the leading contribution to the entropy density for the extremely degenerate case supplies constraints on the p M −1,0 and p M −1,N coefficients. These are also listed in Table 1.
We choose the internal energy density and number density to be thermodynamically consistent with equation (12), so by differentiation we find (

Results
We determined the polynomial coefficients p mn by a least squares fit to data for the pressure, which was obtained by accurate numerical integration using the same (f, g) grid as in EFF. This was done in order to facilitate comparisons with EFF. The accuracies of the number density are of the same order as those for the pressure, so this procedure yields acceptable results for the internal energy and number density as well.
Table 2 also shows that treating a as a free parameter significantly improves the fit as long as both M and N are greater than 1. For example, for the M = N = 3(2) case the improvement amounts to a factor of 8 (4). For fermions the number of coefficients p mn to be fitted is (N + 1)(M + 1).
The N = M = 2 case, whose coefficients are listed in Table 3, shows RMS deviations of 0.03%, which may be sufficient for many purposes. However, an order of magnitude improvement in the accuracy can be obtained by going to the N = M = 3 approximation, which is displayed in Table 4.
In many situations, it is preferable to have approximations that smoothly merge onto the exact limits. This prevents discontinuities from occurring when the exact results are used instead of the approximate scheme in extreme situations. We performed additional fits in which the 6 constraints from Table   1 were incorporated. For the M = N = 3 case (with 10+1 free parameters), the deviations, denoted by a dagger in Table 2, are about 4 times larger than in the unconstrained fit. Nevertheless, the overall accuracy is still acceptable for most numerical work, and if more accuracy is desired, one could increase M or N. Note that for the case M = N = 3, the improvement in the pressure error gained by fitting a amounts to a factor of 9. We list the coefficients for this case in Table 5.

Limiting Cases
The approximation to boson integrals follows by analogy to that for fermions.
However, the degenerate case is now denoted by |ψ| ≪ 1 since ψ must always be negative. Furthermore, a boson-antiboson gas in equilibrium has the additional restriction that ψt ≥ −2, which would preclude the existence of the non-degenerate, extremely relativistic case, although it is nevertheless useful to consider this formal limit below.
We first determine the expansions in the various limiting cases.

Non-degenerate
The expression for the pressure in the non-degenerate (ψ ≪ −1) case was given in equation (8) and from this the internal energy density and density are easily obtained.
For the non-relativistic case, in which t ≪ 1, the asymptotic series for the modified Bessel functions yields ND NR (18) As we have remarked, we consider the extremely relativistic (t ≫ 1) case even though it is not physical when one is dealing with a boson-antiboson gas in thermal equilibrium. Making the small argument expansion of the modified Bessel functions, one obtains

Extremely Degenerate
In the non-relativistic case (|ψ| ≪ 1, t ≪ 1), we use y = E/t as the integration variable in equation (3) to find The integrals can be written in terms of a function Φ which is defined by Erdélyi et al. (1953) as where ζ(z) is the Riemann zeta function. For the series expansion to be valid, |ψ| < 2π and s = 1, 2, 3 . . .. This function is related to the polylogarithm function employed by Haber & Weldon (1982a, b); Φ(e ψ , s, 1) = e −ψ Li s (e ψ ).

Numerical Approximation
We now follow the strategy employed for fermions by EFF. The preceding formulae suggest that for the non-relativistic (relativistic) case, we need a power series in t (t −1 ). Thus, a suitable choice is simply g = t. In the nondegenerate limit, powers of e ψ are required. For the degenerate case, powers of √ −ψ are indicated by equation (22); the situation is less clear for equation (23), which has a complicated structure, but it would seem reasonable to use a similar approach. We are therefore led to consider a variable h such that This does not uniquely specify the derivative and, after some experimentation, we chose which has the desired limiting behavior. The introduction of the parameter a allows greatly improved fits, as we have also found in the fermion case (see §3). This property is related to the limiting behavior of the number density and entropy in the degenerate limit. Integration then yields For small values of h and ψ, h ≃ √ −2aψ, while for large values h ≃ √ ae 1−ψ ; note that h ≥ 0 for all ψ ≤ 0.
We then write the polynomial expansion for the pressure where the coefficients of the polynomial are denoted by p mn . (Note that in the case h = 0, the interpretation h m = δ 0m is to be made.) For the ED case, where h and |ψ| are small, p should not contain a term proportional to h, which would yield a leading divergent 1/ √ −ψ contribution to the number density. Therefore, we must impose the condition p 1n = (M + 1)p 0n so that the contribution to equation (27) which is linear in h vanishes. Equation (27) then has the correct behavior in the four limiting cases and should be suitable for smoothly interpolating for all values of ψ and t. Since it is desirable that the internal energy density and number density be thermodynamically consistent with the pressure, we differentiate equation (27) to obtain The entropy density can then be obtained from s = (p + u)/t − ρψ. We note that the condition p 1n = (M + 1)p 0n leads to explicit cancellation of the leading order contributions to the entropy in the degenerate limit.
Note that equation (27) for bosons is actually the same as equation (12) for fermions as may be seen by setting f = 1/h and using the boson parameterization g = t. It is simply more convenient to use these forms so that infinite values of h for bosons or f for fermions are only encountered for infinitely large |µ|. In the boson case, it is necessary to be able to specify the case ψ = 0, which is just h = 0.
The leading terms for the limiting expressions for p and ρ, given in §4.1, yield values for selected coefficients p mn and these are collected in Table 6. For reference, we note that ζ(3/2) = 2.61238, ζ(5/2) = 1.34149, ζ(3) = 1.20206, and ζ(4) = π 4 /90. The fitted values for these coefficients should approach the values in Table 6 with increasing dimensions M and N, as indeed we find.
For the case M = 2 we have two different requirements on the coefficients, derived from the limiting cases for p and ρ. For p 20 these expressions agree for a = 0.9486 or a = 0.9999, while for p 2N they agree if a = 0.3854 or a = 0.9005. Using a = 0.91, for example, the expressions agree to ≤ 0.25% in both cases. On the other hand, for a = 1 the expressions for p 20 agree to 0.0003%, but the ones for p N 0 agree to only 2.7%. This large error explains why it is advantageous to employ a as an additional fitting parameter.

Results
The polynomial coefficients p mn in equations (27) and (28) were obtained by a least squares fit to "exact" results obtained by numerical Gaussian integration (the accuracy was at least six decimal places). The results were generated for a (t, ψ) grid consisting of t values of .05, .1, .5, 1, 5, 10 and 100 and, for each t, seven equally spaced values of ψ in the range −2 ≤ ψt ≤ 0.
The overall quality of the results obtained for various dimensionalities M and N is indicated in Table 7 in terms of the maximum modulus and the root mean square deviations, as for fermions. If the coefficients p mn are obtained by fitting just the pressure data, then a comparable fit is obtained for the internal energy. However the fit for the number density can be significantly improved by separately optimising the p mn coefficients. Therefore the numbers in Table 7 were obtained by simultaneously fitting the pressure and number density data, thus sacrificing some accuracy in p to achieve improvement in ρ (both the sacrifice and improvement are typically a factor of ∼ 4 in the RMS deviations). In addition to minimizing with respect to the (N + 1)M coefficients p mn , the results were minimized with respect to the parameter a of equation (26). Although a turned out to be fairly close to 1, there is a substantial improvement gained by allowing it to vary. For the cases listed a was in the range 0.7 − 1.1.
The simplest approximation that is sensible, M = 2 N = 1, produces results accurate to better than 5% with only 4+1 parameters, while the M = N = 2 case with 6 + 1 parameters achieves accuracies of < ∼ 1 2 %. These cases are sufficiently simple that they could be employed in hand calculations to yield reasonable estimates. We therefore tabulate the coefficients in Tables 8   and 9, respectively. Increasing the number of parameters obviously improves the accuracy and we would judge that the M = 3 N = 4 case is sufficient for most purposes. Here the MM deviations are < 0.02% and the RMS deviations are ∼ 0.003%. The 15 + 1 parameters are given in Table 10.
As in the case of fermions, it might be preferable, or even necessary, to have an approximation which merges exactly into the asymptotic limiting cases. Thus, it is natural to employ the values of Table 6 as constraints and fit the remaining parameters. For M = 3 and N = 4 this leaves 9+1 free parameters. The deviations, denoted by a dagger in Table 7, are, of course, larger than in the unconstrained case, but only by a factor of < ∼ 2 (the degradation was slightly worse for fermions). The coefficients for this case, tabulated in Table 11, are similar to those of Table 10 since the errors do not differ greatly.

CONCLUDING REMARKS
We have generated polynomial approximations to the thermodynamic integrals which embody and merge exactly into the analytic limiting cases and accurately interpolate over the entire range of temperature and chemical potential. The boson case had not been treated previously and we have improved upon the parameterization of EFF and the modification of Pols et al. (1995) for fermions. With 10 (11) free parameters, and 20 (16) terms, we achieve accuracies of better than 0.008% (0.03%) for the boson (fermion) pressure. Some simpler, but less accurate, approximations were also discussed.
In astrophysics, it is often the case that particles and antiparticles are in thermal equilibrium. Then, the net number density is N = N + − N − where N + (N − ) is the number density of particles (antiparticles), and the total pressure is P = P + + P − . In thermal equilibrium, we must have µ + = −µ − . The problem usually posed is such that one must determine the thermodynamic quantities given N and T . An obvious approach would be to employ the EFF scheme and to establish the limits of the pressure in analogy to equation (9) in the fermion case. Unfortunately, the non-degenerate, non-relativistic case has no single limit in terms of power laws because of pair formation. This is true in both the fermion and boson cases. Thus, the EFF scheme cannot be directly applied to pairs in equilibrium.
It is necessary, therefore, to treat pairs in terms of the particles and antiparticles separately. Thus, we can use the results of §3 for the number densities N + (f + , T ) and N − (f − , T ) of fermions and antifermions, respectively.
Equation (15) gives ψ ± in terms of f ± . The pair of simultaneous equations is then solved for f + and f − , as opposed to the inversion of a single equation N = N (f, T ) for f . For bosons, one simply replaces f ± by h ± in equation Table 1. Fermion coefficients, p mn , derived from the leading terms in the limiting expressions for p and s. The notation 5.9(−4), for example, indicates 5.9 × 10 −4 . † Employing the constraints of Table 1.   Table 6. Boson coefficients, p mn , derived from the leading terms in the limiting expressions for P and ρ.  Table 6.   Table 11. Boson coefficients p mn for M = 3, N = 4 including the constraints of Table 6; a = 1.040.