MATHEMATICS OF COMPUTATION VOLUME 44, NUMBER 169 JANUARY 1985, PAGES 177-190 Gaussian Quadrature Involving Einstein and Fermi Functions With an Application to Summation of Series* By Walter Gautschi and Gradimir V. Milovanovic Abstract. Polynomials trk() = trk(; d\), k = 0,1,2,..., are constructed which are orthogo- nal with respect to the weight distributions d\(t) = (t/(e' - l))r dt and d\(t) = (I/O' + V)Ydt, r = 1,2, on (0, oo). Moment-related methods being inadequate, a discre- tized Stieltjes procedure is used to generate the coefficients ak, ßk in the recursion formula "¡t+i(0 = (' - «*K(0 - /V*-i(0. * - 0,1,2,..., ir0(f) = 1, «■_!(<) - 0. The discretiza- tion is effected by the Gauss-Laguerre and a composite Fejér quadrature rule, respectively. Numerical values of ak, ßk, as well as associated error constants, are provided for 0 < k ^ 39. These allow the construction of Gaussian quadrature formulae, including error terms, with up to 40 points. Examples of n-point formulae, n = 5(5)40, are provided in the supplements section at the end of this issue. Such quadrature formulae may prove useful in solid state physics calculations and can also be applied to sum slowly convergent series. 1. Introduction. We are interested in Gaussian quadrature on [0, oo] relative to the weight functions er(t) = (t/(e' - l))r and <pr(f) = (l/(e' + l))r. These functions arise, for example, in solid state physics and are referred to, when r = 1, as Einstein and Fermi functions, respectively. Integrals with respect to the measure dX(t) = er(t) dt, r = 1 and 2, are widely used in phonon statistics and lattice specific heats [7, §10], [1, §2.4], and occur also in the study of radiative recombination processes [9, §9.2]. Specifically, the average energy of a quantum harmonic oscillator of frequency « at temperature T (representing thermal vibrations of crystal lattice atoms) is given by — h<¿ (1.1) u = exp(hu/kT) - 1 ' where h = h/2m (h is the Planck constant) and k is the Boltzmann constant. Therefore, U = kTex(hu/kT), where ex is Einstein's function. Letting g(u) denote the phonon density of states function, and integrating (1.1) over the entire frequency range, the total energy of thermal vibration of the crystal lattice becomes ReceivedApril 11, 1983; revisedNovember14, 1983and April 3, 1984. 1980Mathematics Subject Classification. Primary33A65, 65D32;Secondary 65A05, 65B10. * The work of the first author was sponsored in part by the National Science Foundation under grant MCS-7927158A1. ©1985 American Mathematical Society 0025-5718/85 $1.00 + $.25 per page 177