Calculating the continued fraction coefficients of a sub-diagonal Padé approximant at arbitrary order

Abstract: The inspiral of two compact objects in gravitational wave astronomy is described by a post-Newtonian expansion in powers of $(v/c)$. In most cases, it is believed that the post-Newtonian expansion is asymptotically divergent. A standard technique for accelerating the convergence of a power series is to re-sum the series by means of a rational polynomial called a Pad\'e approximation. If we liken this approximation to a matrix, the best convergence is achieved by staying close to a diagonal Pad\'e approximation. This broadly presents two subsets of the approximation : a super-diagonal approximation $P^M_N$ and a sub-diagonal approximation $P_M^N$, where $M = N+\epsilon$, and $\epsilon$ takes the values of 0 or 1. Left as rational polynomials, the coefficients in both the numerator and denominator need to be re-calculated as the order of the initial power series approximation is increased. However, the sub-diagonal Pad\'e approximant is computationally advantageous as it can be expressed in terms of a Gauss-like continued fraction. Once in this form, each coefficient in the continued fraction is uniquely determined at each order. This means that as we increase the order of approximation of the original power series, we now have only one new additional coefficient to calculate in the continued fraction. While it is possible to provide explicit expressions for the continued fraction coefficients, they rapidly become unwieldy at high orders of approximation. It is also possible to numerically calculate the coefficients by means of ratios of Hankel determinants. However, these determinants can be ill-conditioned and lead to numerical instabilities. In this article, we present a method for calculating the continued fraction coefficients at arbitrary orders of approximation.