How to solve this set of definite integrals? I'm looking to solve the integrals
$$
I_n=\int_0^\infty \frac{x^ne^{-2ax}}{\sqrt{x^4+1}}dx.
$$
for $a>0$.
This should be reduced to finding $I_0$ as numerically I find that
$$
I_n=\left(-2\right)^{-n}\frac{d^nI_0}{da^n}.
$$
Mathematica gives a result in terms of the MeijerG function $G^{51}_{15}\left(x|^{\frac{1}{4}}_{-\frac{1}{2},-\frac{1}{4},-\frac{1}{4},0,\frac{1}{4}}\right)$ but I cannot manage to prove it. I have looked for integral representations of this and hypergeometric functions in Gradshteyn but couldn't find anything useful.
Any help is much welcomed.
Thanks
 A: One way to get the result is to use the Mellin transform of the integral considered as a function of $a$:
\begin{align}
 I_0&=\int_0^\infty \frac{e^{-2ax}}{\sqrt{x^4+1}}\,dx\\
 \mathcal{M}\left[ I_0\right]&=\int_0^\infty a^{s-1}\,da\int_0^\infty \frac{e^{-2ax}}{\sqrt{x^4+1}}\,dx\\
 &=\int_0^\infty \frac{1}{\sqrt{x^4+1}}\,dx \int_0^\infty a^{s-1}e^{-2ax}\,da\\
 &=2^{-s}\Gamma(s)\int_0^\infty \frac{x^{-s}}{\sqrt{x^4+1}}\,dx
\end{align}
which is valid for $\Re(s)>0$. Now, changing $x=y^{1/4}$ in the integral, one obtains
\begin{align}
 \mathcal{M}\left[ I_n\right]&=2^{-s-2}\Gamma\left(s\right)\int_0^\infty\frac{y^{\frac{-s-3}{4}}}{\sqrt{y+1}}\,dy\\
 &=2^{-s-2}\Gamma(s)B\left(\frac{-s+1}{4},\frac{s+1}{4}\right)\\
 &=\frac{2^{-s}}{4\sqrt{\pi}}\Gamma\left(s\right)\Gamma\left(\frac{-s+1}{4}\right)\Gamma\left(\frac{s+1}{4}\right)
\end{align}
(The integral is the Mellin transform of $1/\sqrt{y+1}$ taken at $(n-s+1)/4$). This result is valid for $0<\Re(s)<1$. Taking the inverse transform, 
\begin{equation}
 I_0=\frac{1}{4\sqrt{\pi}}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(2a)^{-s}\Gamma\left(s\right)\Gamma\left(\frac{-s+1}{4}\right)\Gamma\left(\frac{s+1}{4}\right)\,ds
\end{equation}
with $0<\sigma<1$. Changing $s=4t$, one can express
\begin{equation}
 I_n=\frac{1}{\sqrt{\pi}}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(16a^4)^{-t}\Gamma\left(4t\right)\Gamma\left(t+\frac{1}{4}\right)\Gamma\left(\frac{1}{4}-t\right)\,dt
\end{equation} 
Expanding $\Gamma(4t)$ using the duplication formula:
\begin{equation}
 I_0=\frac{\sqrt{2}}{8\pi^2}\frac{1}{2i\pi}\int_{\sigma-i\infty}^{\sigma+i\infty}(\frac{a^4}{16})^{-t}\Gamma\left(t\right)\Gamma^2\left( t+\frac{1}{4} \right)\Gamma\left( t+\frac{1}{2} \right)\Gamma\left( t+\frac{3}{4} \right)\Gamma\left(\frac{1}{4}-t\right)\,dt
\end{equation} 
With $b_1=\frac{1}{4},a_1=1,a_2=\frac{3}{4},a_3=\frac{3}{4},a_4=\frac{1}{2},a_5=\frac{1}{4}$ and $\sigma=1/8$, we can express the result using the integral representation of the Meijer function DLMF:
\begin{equation}
 I_0=\frac{\sqrt{2}}{8\pi^2}{G^{1,5}_{5,1}}\left(\left.\frac{16}{a^4}\right|{1,\frac{3}{4},\frac{3}{4},\frac{1}{2},\frac{1}{4}\atop \frac{1}{4}}\right)
\end{equation} 
which, using these identities, can be transformed into
\begin{align}
 I_0&=\frac{\sqrt{2}}{8\pi^2}{G^{5,1}_{1,5}}\left(\left.\frac{a^4}{16}\right|{\frac{3}{4}\atop 0,\frac{1}{4},\frac{1}{4},\frac{1}{2},\frac{3}{4}}\right)\\
 &=\frac{\sqrt{2}}{32\pi^2}a^2{G^{5,1}_{1,5}}\left(\left.\frac{a^4}{16}\right|{\frac{1}{4}\atop -\frac{1}{2},-\frac{1}{4},-\frac{1}{4},0,\frac{1}{4}}\right)
\end{align} 
A: It's a general result that an integral over $[0, \infty)$ of two Meijer G-function of rational powers of $x$ is again a G-function:
$$\int_0^\infty \frac {x^n e^{-2 a x}} {\sqrt {x^4 + 1}} dx =
\frac 1 {\sqrt \pi} \int_0^\infty x^n
 G_{0, 1}^{1, 0} \left( 2 a x \middle| {- \atop 0} \right)
 G_{1, 1}^{1, 1} \left( x^4 \middle| {\frac 1 2 \atop 0} \right) dx = \\
\frac {2^{n - 3/2}} {\pi^2 a^{n + 1}}
 G_{5, 1}^{1, 5} \left( \frac {16} {a^4} \middle|
  {\frac {-n} 4, \frac {-n + 1} 4, \frac {-n + 2} 4, \frac {-n + 3} 4, \frac 1 2
   \atop 0} \right).$$
The resulting G-function is not reducible to a sum of hypergeometric functions because one of the four coefficients $(-n + i)/4$ is half-integral for any integral $n$.
