Similar to my answer here, we can use a generalization of Ramanujan's master theorem called the method of brackets to derive some of the results in Claude Leibovici's answer.
Assume that $b \ge0 $, $a >2b$, and $m \in \mathbb{Z} >-1$.
The exponential function has the series representation $$e^{- ax} = \sum_{j=0}^{\infty} \frac{(-1)^{j}}{j!}(ax)^{j} = \sum_{j=0}^{\infty} \phi{_j}(ax)^{j}. $$
And the square of the Bessel function of the first kind of order zero has the series representation $$J_{0}(bx)^{2} = {_2F_1} \left(\frac{1}{2};1,1;-(bz)^{2} \right) = \sum_{k=0}^{\infty} \phi_{k} \frac{\Gamma \left(k+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(k+1)^{2}} \, (bx)^{2k}.$$
Therefore, the integral $\int_{0}^{\infty}e^{-ax} x^{m} J_{0}(bx)^{2} \, \mathrm dx$ has the corresponding bracket series $$\sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \phi_{j,k} \, a^{j} \, b^{2k} \, \frac{\Gamma \left(k+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(k+1)^{2}} \, \langle j+2k+m+1 \rangle.$$
Since there are $2$ indices but only $1$ bracket, there are two cases to examine.
If we let $k$ be a free parameter, then the bracket vanishes if $j=-2k-m-1$, and we have the contribution $$\begin{align} S_{1} &= \sum_{k=0}^{\infty} \frac{(-1)^{k}}{k!} a^{-2k-m-1} \, b^{2k} \frac{\Gamma \left(k+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(k+1)^{2}} \, \Gamma \left(
2k+m+1 \right) \\ &= \sum_{k=0}^{\infty} \frac{1}{k!} a^{-(m+1)} \, \left( -\frac{b^{2}}{a^{2}} \right)^{k} \, \frac{\Gamma \left(k+ \frac{1}{2} \right)}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma(k+1)^{2}} \, \frac{2^{2k+m}}{\Gamma \left(\frac{1}{2} \right)} \, \Gamma \left(k+ \frac{m}{2} + \frac{1}{2} \right) \Gamma \left(k+ \frac{m}{2}+1 \right) \\ &= a^{-(m+1)} \, \frac{2^{m}}{\Gamma \left(\frac{1}{2} \right)} \, \Gamma \left(\frac{m}{2}+ \frac{1}{2} \right) \Gamma \left( \frac{m}{2}+1 \right) \, {_3F_2} \left(\frac{1}{2}, \frac{m+1}{2}, \frac{m+2}{2}; 1,1; -\frac{4b^{2}}{a^{2}} \right) \\ &= a^{-(m+1)} \, 2^{m} \, \Gamma \left(2 \left(\frac{m}{2}+ \frac{1}{2} \right) \right) 2^{-2(m/2+1/2)+1} \, {_3F_2} \left(\frac{1}{2}, \frac{m+1}{2}, \frac{m+2}{2}; 1,1; -\frac{4b^{2}}{a^{2}} \right) \\ &= a^{-(m+1)} \, \Gamma(m+1) \, {_3F_2} \left(\frac{1}{2}, \frac{m+1}{2}, \frac{m+2}{2}; 1,1; -\frac{4b^{2}}{a^{2}} \right). \end{align}$$
On the second and fourth lines I used Legendre's duplication formula.
If we let $j$ be a free parameter, then the bracket vanishes if $k=- \frac{j}{2} - \frac{m}{2} - \frac{1}{2}$, and we have the contribution $$ S_{2} = \frac{1}{2} \sum_{j=0}^{\infty}\frac{(-1)^{j}}{j!} \, a^{j} \, b^{-j-m-1} \, \frac{\color{red}{\Gamma \left(- \frac{j+m}{2} \right)}}{\Gamma \left(\frac{1}{2} \right)} \frac{1}{\Gamma \left(\frac{-j-m+1}{2} \right)^{2}} \, \Gamma \left(\frac{j+m+1}{2} \right). $$
Since not all terms of this series are finite, we discard it.
Therefore, we have $$\int_{0}^{\infty}e^{-ax} x^{m} J_{0}(bx)^{2} \, \mathrm dx = S_{1} = \Gamma(m+1) \, a^{-(m+1)}\, {_3F_2} \left(\frac{1}{2}, \frac{m+1}{2}, \frac{m+2}{2}; 1,1; -\frac{4b^{2}}{a^{2}} \right)$$ for $b \ge 0$, $a>2b$, and $m \in \mathbb{Z} >-1$.
Now assume that $b \ge 0$, $a>2b$, and $m \in \mathbb{Z} > -2$.
The product $J_{0}(bx) J_{1}(bx)$ has the series representation $$ \begin{align} J_{0}(bx)J_{1}(bx) &= \frac{bx}{2} \, {_2F_3} \left( \frac{3}{2}, 1; 1, 2, 2; (bx)^{2} \right) \\ &= \frac{bx}{2} \, {_1F_2} \left( \frac{3}{2}; 2, 2; (bx)^{2} \right) \\ &= \frac{bx}{2} \sum_{k=0}^{\infty} \phi_{k} \frac{\Gamma \left(k+ \frac{3}{2} \right)}{\Gamma \left(\frac{3}{2} \right)} \frac{1}{\Gamma \left(k+2 \right)^{2}} (bx)^{2k}. \end{align}$$
Therefore, the integral $\int_{0}^{\infty} e^{-ax} x^{m} J_{0}(bx)J_{1}(bx) \, \mathrm dx$ has the corresponding bracket series $$ \frac{b}{2} \sum_{j=0}^{\infty} \sum_{k=0}^{\infty} \phi_{j,k} \, a^{j} \, b^{2k} \, \frac{\Gamma \left(k+ \frac{3}{2} \right)}{\Gamma \left(\frac{3}{2} \right)} \frac{1}{\Gamma(k+2)^{2}} \, \langle j+2k+m+2 \rangle.$$
There are again two cases to examine.
If we let $j$ be a free parameter, then the bracket vanishes if $ j= -2k-m-2$, and we have the contribution $$ \begin{align} S_{3} &= \frac{b}{2} \sum_{k=0}^{\infty} \frac{(-1)^{k}}{k!} \, a^{-2k-m-2} \, b^{2k} \, \frac{\Gamma \left(k+ \frac{3}{2} \right)}{\Gamma \left(\frac{3}{2} \right)} \frac{1}{\Gamma(k+2)^{2}} \Gamma \left(2k+m+2 \right) \\ &= \small \frac{b}{2} \sum_{k=0}^{\infty} \frac{1}{k!} a^{-(m+2)} \left(-\frac{b^{2}}{a^{2}} \right)^{k} \, \frac{\Gamma \left(k+ \frac{3}{2} \right)}{\Gamma \left(\frac{3}{2} \right)} \frac{1}{\Gamma(k+2)^{2}} \frac{2^{2k+m+1}}{\Gamma \left(\frac{1}{2} \right)} \, \Gamma \left(k+ \frac{m}{2} + 1\right) \Gamma \left(k+ \frac{m}{2}+\frac{3}{2} \right) \\ &= \frac{b}{2} \frac{a^{-(m+2)} 2^{m+1}}{\Gamma\left(\frac{1}{2} \right)} \, \Gamma \left(\frac{m}{2}+1 \right) \Gamma \left(\frac{m}{2} + \frac{3}{2} \right) {_3F_2} \left(\frac{3}{2}, \frac{m+2}{2}; \frac{m+3}{2};2,2; - \frac{4b^{2}}{a^{2}} \right) \\ &= \frac{b}{2} \, \frac{a^{-(m+2)} 2^{m+1}}{\Gamma \left(\frac{1}{2} \right)} \, \Gamma \left(\frac{1}{2} \right) 2^{-m-1} \, \Gamma(m+2) \, {_3F_2} \left(\frac{3}{2}, \frac{m+2}{2}; \frac{m+3}{2};2,2; - \frac{4b^{2}}{a^{2}} \right) \\ &= \frac{b}{2} \, a^{-(m+2)} \, \Gamma(m+2) \, {_3F_2} \left(\frac{3}{2}, \frac{m+2}{2}; \frac{m+3}{2};2,2; - \frac{4b^{2}}{a^{2}} \right). \end{align}$$
If we let $k$ be a free parameter, we get a series like before that doesn't have all finite terms.
Therefore, we have $$\small \int_{0}^{\infty} e^{-ax} x^{m} J_{0}(bx) J_{1}(bx) \, \mathrm dx = S_{3} = \frac{b}{2} \, a^{-(m+2)} \, \Gamma(m+2) \, {_3F_2} \left(\frac{3}{2}, \frac{m+2}{2}, \frac{m+3}{2};2,2; - \frac{4b^{2}}{a^{2}} \right)$$ for $b \ge 0$, $a > 2b$, and $m \in \mathbb{Z} > -2$.
Since the Laplace transform defines an analytic function in the half-plane that it converges absolutely, analytic continuation extends both results to $\Re(a) >0$.
