  We develop a new type of orthogonal polynomial, the modified discrete
Laguerre (MDL) polynomials, designed to accelerate the computation of bosonic
Matsubara sums in statistical physics. The MDL polynomials lead to a rapidly
convergent Gaussian "quadrature" scheme for Matsubara sums, and more generally
for any sum $F(0)/2 + F(h) + F(2h) + \cdots$ of exponentially decaying summands
$F(nh) = f(nh)e^{-nhs}$ where $hs>0$. We demonstrate this technique for
computation of finite-temperature Casimir forces arising from quantum field
theory, where evaluation of the summand $F$ requires expensive electromagnetic
simulations. A key advantage of our scheme, compared to previous methods, is
that the convergence rate is nearly independent of the spacing $h$
(proportional to the thermodynamic temperature). We also prove convergence for
any polynomially decaying $F$.
