Q:

2D Fourier transform of $1/(1+|\mathbf{r}|^6)$

In some physics literature I've been studying I ran into the Fourier transform of the function
$$f(\mathbf{r}) = \frac{1}{1+|\mathbf{r}|^6},$$
where $\mathbf{r}$ is a two dimensional vector. These papers (#1, see p. 21; and #2, see p. 3) claim this equals
$$\mathcal{F}[f(\mathbf{r})] = \left(\frac{2\pi^2}{3} \right) \left(\frac{e^{-k/2}}{k} \right) \left[e^{-k/2} - 2\sin(\pi/6 - \sqrt{3}k/2) \right],$$
where $k = |\mathbf{k}| = |(k_1,k_2)|$ with these being the Fourier variables corresponding to $\mathbf{r} = (x,y)$.
Unfortunately, I cannot find any more details on the derivation of this result. Understanding it would hopefully help me understand some more of the physics. I've tried carrying out the Fourier integral myself, but to no avail (although I'm far from an expert on the topic). The FourierTransform function of Mathematica is also of disappointingly little help, despite it seeming like such a manageable result.
I would appreciate any suggestions and/or pointers to transform rules that I may be missing that would make life easier.

A:

I) It is true that the 2D Fourier transform is
$$ \begin{align} \hat{h}({\bf k})~&=~\iint_{\mathbb{R}^2} \! \frac{\mathrm{d}^2r}{2\pi} \exp\left(\pm i {\bf r}\cdot{\bf k} \right) h(r) \\
~&=~\iint_{\mathbb{R}^2} \! \frac{\mathrm{d}^2r}{2\pi} \cos\left({\bf r}\cdot{\bf k} \right)h(r) \\
~&=~\int_{0}^{\infty} \! \mathrm{d}r~r J_{0}(kr) h(r) \tag{1} \end{align}$$
for suitable functions $h$ that are independent of the polar variable $\theta$, cf. integral (9.1.18) in Ref. 3, Wikipedia, and above comment by Hans Engler. However, we shall take a different route in this answer.
II) As a warmup to OP's question, consider the following function$^1$
$$g({\bf r})~=~\frac{1}{r^2+m^2}, \qquad {\rm Re}(m)~>~0. \tag{2}$$
The 2D Fourier transform reads
$$ \hat{g}({\bf k})
~=~\iint_{\mathbb{R}^2} \! \frac{\mathrm{d}^2r}{2\pi} \frac{\cos\left({\bf r}\cdot{\bf k} \right)}{r^2+m^2}
~=~\left. \int_{\mathbb{R}}\! \frac{\mathrm{d}x}{2\pi} \cos\left(xk \right) \int_{\mathbb{R}}\! \frac{\mathrm{d}y}{y^2+E^2} \right|_{E=\sqrt{x^2+m^2}}$$
$$~=~\left. \int_{\mathbb{R}}\! \frac{\mathrm{d}x}{2\pi} \frac{\cos\left(xk \right)}{2iE} \int_{\mathbb{R}}\! \mathrm{d}y\left( \frac{1}{y-iE} -\frac{1}{y+iE}  \right)\right|_{E=\sqrt{x^2+m^2}}$$
$$~=~\left. \int_{\mathbb{R}}\! \frac{\mathrm{d}x}{2\pi} \frac{\cos\left(xk \right)}{2iE} \left[   {\rm Ln}(y-iE) -{\rm Ln}(y+iE)   \right]_{y=-\infty}^{y=\infty}\right|_{E=\sqrt{x^2+m^2}}$$
$$~=~\int_{\mathbb{R}}\! \frac{\mathrm{d}x}{2} \frac{\cos\left(xk \right)}{\sqrt{x^2+m^2}}~=~K_0(m k)~\sim~\sqrt{\frac{\pi}{2km}}e^{-km}\quad\text{for}\quad k~\to~ \infty, \tag{3}$$
where we in the last equality used integral (9.6.25) in Ref. 3. Here $K_0$ is the modified Bessel function of 2nd kind. 
III) Now let us return to OP's question. From the partial fraction decomposition
$$f({\bf r})~=~\frac{1}{r^6+1}
~=~\frac{1}{3}\sum_{n=-1}^1 \frac{\omega^{2n}}{r^2+\omega^{2n}},\qquad \omega~:=~\exp\left[\frac{i\pi}{3}\right]~=~\frac{1+i\sqrt{3}}{2},\tag{4}$$
the 2D Fourier transform reads
$$ \hat{f}({\bf k})
~=~\iint_{\mathbb{R}^2} \! \frac{\mathrm{d}^2r}{2\pi} \frac{\cos\left({\bf r}\cdot{\bf k} \right)}{r^6+1}
~=~\frac{1}{3}\sum_{n=-1}^1 \omega^{2n}K_0(\omega^n k)  $$
$$\sim~\frac{1}{3}\sqrt{\frac{\pi}{2k}}\left(e^{-k}+2e^{-k/2}\sin\left(\frac{k\sqrt{3}}{2}\right)\right)\quad\text{for}\quad k~\to~ \infty, \tag{5}$$
which is our main result.
IV) Finally in Table 1 below we check numerically the asymptotic expansion (5) (third column) compared to the exact formula (fourth column).
In the second column we have calculated the characteristic sine factor from OP's Refs. 1 & 2. Note that Refs. 1 & 2 get the wrong sign, which can not be explained away by a different normalization convention.
$\downarrow$ Table 1. 
$$
\begin{array}{ccccccc} 
\text{Wave number } k&& \frac{2}{3}\sqrt{\frac{\pi}{2k}}e^{-k/2}\sin\left(\frac{k\sqrt{3}}{2}-\frac{\pi}{6}\right)&& \frac{2}{3}\sqrt{\frac{\pi}{2k}}e^{-k/2}\sin\left(\frac{k\sqrt{3}}{2}\right)   &&\text{Exact} ~\hat{f}(k)  \\ \\
4 &&1.1 \cdot 10^{-2} && -1.8 \cdot 10^{-2}&& -1.3 \cdot 10^{-2}  \\
11 &&4.2 \cdot 10^{-4} && -1.0 \cdot 10^{-4} && -9.2 \cdot 10^{-5}  \\
15 &&-1.2 \cdot 10^{-5} && 4.9 \cdot 10^{-5} && 4.8 \cdot 10^{-5}  \\
22 &&-9.4 \cdot 10^{-7} && 6.0 \cdot 10^{-6}&& 5.8 \cdot 10^{-6}  \\ 
\end{array}
$$
Conclusion: Taking into account possible different normalizations of the 2D Fourier transform, the main formula (5) does not agree with OP's Refs. 1 & 2, not even as an asymptotic large $k$ expansion. 
References:

X. Antoine & R. Duboscq, Computer Physics Communications 193 (2015) 95; eq. (5.43).
C.-H. Hsueh, T.-C. Lin, T.-L. Horng, & W.C. Wu,
Phys. Rev. A 86 (2012) 013619; p. 3.
Abramowitz & Stegun, Handbook of Mathematical Functions.

--
$^1$ Formulas (2) and (3) are well-known in physics as the Euclidean 2D propagator/Greens function for a massive field with mass $m$, although the roles of the Fourier variables ${\bf r}\leftrightarrow  {\bf k}$ are reversed. 

