I'm trying to find an analytic form for the Hankel transform of a Sersic profile: $$ \begin{align} f_n(r) &= \frac{e^{-r^{1/n}}}{n \Gamma(2n)}\\ F_n(k) &= \int_0^\infty f(r) r J_0(k r) dr \end{align} $$
Background: This profile is commonly used in astronomy to model galaxy surface brightness profiles. $n=1/2$ is a Gaussian. $n=1$ is an exponential, which is a good fit for galaxies like the Milky Way. $n=4$ is called a DeVaucouleurs profile and is a good fit to many very large galaxies in clusters. Most galaxies are well fit by a Sersic profile in the range $n=1$ to $n=5$. When drawing these galaxies, we normally draw them in Fourier space, since we need to convolve the profile by a point-spread function, which is easiest to do via multiplication in Fourier space. Thus we generally use the Hankel transform of the profile to directly draw the $F(k)$ values for a given galaxy in our Fourier-space image.
The Hankel transform for n=4 (and higher) is fairly slow to compute, so I have been looking for an analytic formula that I might be able to use. The two values of $n$ for which I know the answer are: $$ \begin{align} F_{1/2}(k) &= e^{-k^2/4}\\ F_1(k) &= (1+k^2)^{-3/2} \end{align} $$ The other possibly helpful bit of information I have is the MacLaurin series in general $$ \begin{align} F_n(k) &= 1 - \frac{1}{4} \frac{\Gamma(4n)}{\Gamma(2n)} k^2 + \frac{1}{64}\frac{\Gamma(6n)}{\Gamma(2n)} k^4 - \frac{1}{2304}\frac{\Gamma(8n)}{\Gamma(2n)} k^6 + \dots \\ &= \sum_{i=0}^\infty \frac{(-1)^i \Gamma(2(i+1)n)}{2^{2i} \Gamma(i+1)^2 \Gamma(2n)} k^{2i} \end{align} $$ Unfortunately, when $n>1$, the coefficients in the MacLaurin series get big very quickly, so its not a useful approximation unless k is very, very small. I'm not sure formally about the radius of convergence, but practically, it diverges pretty quickly. For instance, the series for n=4 is $$ 1 - (6.5 \times 10^7) k^2 + (8.0 \times 10^{16}) k^4 - (7.1 \times 10^{26}) k^6 + (2.7 \times 10^{37}) k^8 - \dots $$ so it's useless above around $k=10^{-5}$.
Can anyone offer some guidance about how I might be able to find an analytic formula that has this MacLaurin series? Or alternatively some transformation I could make to this function that would render the series expansion convergent and therefore tractable as a fitting formula?
Thanks in advance for any tips you might be able to offer.
Note: I tried using a Pade approximant to the MacLaurin series, and it extends the range of accuracy by a factor of a few, but that's about it.