Evaluate $\int_{0}^\infty\frac{2}{({M-1})!}\sqrt{\frac{a}{x}}K_1\left(2\sqrt{\frac{a}{x}}\right)x^{M-1}\frac{e^{-\frac{x}{b}}}{b^M}\,dx$ How to evaluate this integral or at least rewrite it as special (non-elementary) function step by step ?
$$\int_{0}^\infty\frac{2}{({M-1})!}\sqrt{\frac{a}{x}}K_1\left(2\sqrt{\frac{a}{x}}\right)x^{M-1}\frac{e^{-\frac{x}{b}}}{b^M}\,dx$$
Where $M$ is a positive integer that is $M=1,2,3,4,...$ , $a,b$ are two positive number and ${K_1}\left( . \right)$ is the modified Bessel function of order 1.
If you could provide some reference, I would be pretty much graceful !
Thank you for your enthusiasm !

It seems useful to rearrange the integral:
\begin{align}
 I&=\int_{0}^\infty\frac{2}{({M-1})!}\sqrt{\frac{a}{x}}K_1\left(2\sqrt{\frac{a}{x}}\right)x^{M-1}\frac{e^{-\frac{x}{b}}}{b^M}\,dx\\
 &=\frac{2b^{-M}\sqrt{a}}{({M-1})!}\int_{0}^\infty K_1\left(2\sqrt{\frac{a}{x}}\right)x^{M-3/2}e^{-\frac{x}{b}}\,dx\\
 &=\frac{2}{({M-1})!}\left( \frac{a}{b} \right)^M\int_{0}^\infty K_1\left(\frac{2}{\sqrt s}\right)s^{M-3/2}e^{-\frac{a}{b}s}\,ds
\end{align}
which expresses that the integral is a function of $a/b$.
To obtain the Meijer function expression given by @MariuszIwaniuk in the comment, one may use the Mellin-Barnes representation for the modified Bessel function
\begin{equation}
 K_{\nu}\left(z\right)=\frac{(\frac{1}{2}z)^{\nu}}{4\pi i}\int_{c-i\infty}^{c+i\infty}\Gamma\left(t\right)\Gamma\left(t-\nu\right)(\tfrac{1}{2}z)^{-2t}\,dt
\end{equation}
valid for $c>\max⁡(\Re⁡\nu,0),|\arg⁡ z|<\pi/2$. Here, with $\nu=1,z=2/\sqrt s,c>1$, we obtain, by changing the order of integration
\begin{align}
 I&=\frac{2}{({M-1})!}\left( \frac{a}{b} \right)^M\int_{0}^\infty K_1\left(\frac{2}{\sqrt s}\right)s^{M-3/2}e^{-\frac{a}{b}s}\,ds\\
 &=\frac{1}{({M-1})!}\left( \frac{a}{b} \right)^M\frac{1}{2\pi i}
 \int_{c-i\infty}^{c+i\infty}\Gamma\left(t\right)\Gamma\left(t-1\right)\,dt\int_{0}^\infty s^{t+M-2}e^{-\frac{a}{b}s}\,ds
\end{align}
The integral over $s$ gives directly a Gamma function,
\begin{equation}
 I=\frac{1}{({M-1})!}\left( \frac{a}{b} \right)\frac{1}{2\pi i}
 \int_{c-i\infty}^{c+i\infty}\Gamma\left(t\right)\Gamma\left(t-1\right)\Gamma(t+M-1)\left(\frac ba\right)^t\,dt
\end{equation}
From the definition of the Meijer function
\begin{equation}
{G^{m,n}_{p,q}}\left(z;\begin{array}{c}
a_1,\ldots,a_p\\
b_1,\ldots,b_q
\end{array} \right)=\frac{1}{2\pi\mathrm{i}}\int_{L%
}{\frac{\prod\limits_{\ell=1}^{m}\Gamma\left(b_{\ell}-t\right
)\prod\limits_{\ell=1}^{n}\Gamma\left(1-a_{\ell}+t\right)}{\left(\prod\limits_%
{\ell=m}^{q-1}\Gamma\left(1-b_{\ell+1}+t\right)\prod\limits_{\ell=n}^{p-1}%
\Gamma\left(a_{\ell+1}-t\right)\right)}}z^{t}\,dt
\end{equation}
we can use $z=b/a$, $a_1=1,a_2=2,a_3=2-M$, with $m=q=0, n=p=3$ and an integral along the vertical axis from $c-i\infty$ to $c+i\infty$, as the real parts of the the poles of $\Gamma\left(1-a_{\ell}+t\right)$ are all less than $1$. We deduce
\begin{equation}
 I=\frac{1}{({M-1})!}\left( \frac{a}{b} \right)G^{0,3}_{3,0}\left(\frac{b}{a};\begin{array}{c}
1,2,2-M\\
\textrm{---}
\end{array} \right)
\end{equation}
The proposed expression can be retrieved by combining the identity
\begin{equation}
 I=\frac{1}{({M-1})!}\left( \frac{a}{b} \right)G^{3,0}_{0,3}\left(\frac{a}{b};\begin{array}{c}
 \textrm{---}\\
0,-1,M-1
\end{array} \right)
\end{equation}
and the second identity
\begin{equation}
 I=\frac{1}{({M-1})!}G^{3,0}_{0,3}\left(\frac{a}{b};\begin{array}{c}
 \textrm{---}\\
0,1,M
\end{array} \right)
\end{equation}
