EDIT:
I think the OP was referring to the recurrence relation $$\frac{d}{dx} \left(x^{-n}j_{n}(x) \right)=-x^{-n}j_{n+1}(x). \tag{1} $$
From $(1)$ it follows immediately that $$\begin{align}\int_{0}^{\infty}x^{-n} j_{n+1}(x) \, dx &= -x^{-n} j_{n}(x) \Bigg|^{\infty}_{0} \\ &=-x^{n} \sqrt{\frac{\pi}{2x}} \, J_{n+\frac{1}{2}}(x) \Bigg|^{\infty}_{0} \\ &= \lim_{x \downarrow 0} \, x^{-n} \sqrt{\frac{\pi}{2x}} \, J_{n+\frac{1}{2}}(x) \tag{2} \\ &=\lim_{x \downarrow 0} \, x^{-n} \sqrt{\frac{\pi}{2x}}\left(\frac{1}{\Gamma\left(n + \frac{3}{2}\right)} \left(\frac{x}{2}\right)^{n+ \frac{1}{2}} +\mathcal{O}\left(x^{n+ \frac{5}{2}} \right)\right) \tag{3} \\ &= \sqrt{\pi} \, 2^{-n-1} \, \frac{1}{\Gamma \left(n+ \frac{3}{2}\right)} \\ &= \sqrt{\pi} \, 2^{-n-1} \, \frac{2^{n+1}}{\sqrt{\pi} (2n+1)!!} \tag{4}\\ &= \frac{1}{(2n+1)!!}. \end{align}$$
$(1)$: http://dlmf.nist.gov/10.51#E3
$(2)$: https://en.wikipedia.org/wiki/Bessel_function#Asymptotic_forms
$(3)$: https://en.wikipedia.org/wiki/Bessel_function#Bessel_functions_of_the_first_kind:_J.CE.B1
$(4)$: http://mathworld.wolfram.com/DoubleFactorial.html (2)
This might not be the particular approach you were seeking, but it can be evaluated using a property of the Mellin transform.
Ramanujan's master theorem states that if $f(x)$ has an expansion of the form $$ f(x) = \sum_{k=0}^{\infty} \frac{\phi(k)}{k!} (-x)^{k} ,$$
then
$$ \int_{0}^{\infty} x^{s-1} f(x) \, dx = \Gamma(s) \phi(-s) $$
for the values of $s$ for which the integral converges.
The hypergeometric representation of the Bessel function of the first kind of order $\alpha$ is $$ J_{\alpha}(x) = \frac{(\frac{x}{2})^{\alpha}}{\Gamma(\alpha +1)} \ {}_0F_{1} \left(\alpha +1; - \frac{x^{2}}{4} \right) = \frac{(\frac{x}{2})^{\alpha}}{\Gamma(\alpha +1)} \sum_{k=0}^{\infty} \frac{\Gamma(\alpha + 1)}{\Gamma(\alpha +1+k)} \frac{(-\frac{x^{2}}{4})^{k}}{k!} .$$
So for $s+n+1 >0$ and $s <2$, we have
$$ \begin{align} \int_{0}^{\infty} x^{s-1} j_{n+1}(x) \, dx &= \int_{0}^{\infty} x^{s-1} \sqrt{\frac{\pi}{2x}} \, J_{n+\frac{3}{2}} (x) \, dx \\ &= \int_{0}^{\infty} x^{s-1} \sqrt{\frac{\pi}{2x}} \frac{(\frac{x}{2})^{n + \frac{3}{2}}}{\Gamma(n +\frac{5}{2})} \, {}_0F_{1} \left(n + \frac{5}{2}; - \frac{x^{2}}{4} \right) \, dx \\ &= \frac{\sqrt{\pi}}{2^{n+2}} \frac{1}{\Gamma(n+\frac{5}{2})} \int_{0}^{\infty} x^{s+n} \, {}_0F_{1} \left(n + \frac{5}{2}; - \frac{x^{2}}{4} \right) \, dx \\ &= \frac{\sqrt{\pi}}{2^{n+2}} \frac{1}{\Gamma(n+\frac{5}{2})} \int_{0}^{\infty} (2 \sqrt{u})^{s+n} \, {}_0F_{1} \left(n + \frac{5}{2}; - u \right) \, \frac{du}{\sqrt{u}} \\ &= \frac{ \sqrt{\pi} \ 2^{s-2}}{\Gamma(n+\frac{5}{2})} \int_{0}^{\infty} u^{\frac{s+n+1}{2}-1} \, {}_0F_{1} \left(n + \frac{5}{2}; - u \right) \, du \\ &= \frac{ \sqrt{\pi} \ 2^{s-2}}{\Gamma(n+\frac{5}{2})} \, \Gamma \left(\frac{s+n+1}{2} \right) \frac{\Gamma(n+ \frac{5}{2})}{\Gamma(n + \frac{5}{2} - \frac{s+n+1}{2})} \\ &= \sqrt{\pi} \ 2^{s-2} \, \Gamma \left(\frac{s+n+1}{2} \right) \frac{1}{\Gamma(\frac{4-s+n}{2})} . \end{align}$$
Now let $s= -n+1$.
Then
$$ \begin{align} \int_{0}^{\infty} x^{-n} j_{n+1}(x) \ dx &= \sqrt{\pi} \ 2^{-n-1} \Gamma (1) \, \frac{1}{\Gamma \left(n + \frac{3}{2} \right)} \\ &= \sqrt{\pi} \ 2^{-n-1} \frac{2^{n+1}}{\sqrt{\pi} (2n+1)!!} \tag{1} \\ &= \frac{1}{(2n+1)!!} \end{align}$$
$(1)$: http://mathworld.wolfram.com/DoubleFactorial.html (2)
The answer can be represented as an integral which appears at the bottom of this answer.
The first step is to transform the operator $D_x$ to one $D_w$ via the transformations $x=\sqrt{u}$ and then $w=u-a^2.$ We find that
$$D_x=(x^2-a^2)\frac{1}{x}\frac{d}{dx}\frac{1}{x}\frac{d}{dx} + \frac{2}{x} \frac{d}{dx} \to 4D_w,\, D_w=w\,\frac{d^2}{dw^2}+\frac{d}{dw} .$$
The $D_w$ operator is related to the Laguerre polynomials which allows us to make significant progress. The OP's problem then becomes
$$\exp{(i\lambda D_x)}\, j_0(x) \to \exp{(c D_w)}\, j_0(\sqrt{w+a^2}) \,\,\,, (c=4i\lambda). $$
Expand the sinc function, using the abbreviation $y=a^2$,
$$ j_0(\sqrt{w+y}) = \sum_{k=0}^\infty (-1)^k\frac{(w+y)^k}{(2k+1)!}=
\sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)!} \sum_{m=0}^{k} \binom{k}{m} y^{k-m}\,w^m .$$
We now need to see how $\exp(c\,D_w)$ behaves on monomials $w^m.$ Expand the exponential and note how powers of $D_w$ behave:
$$ \big(D_w\big)^0 w^m = w^m $$
$$ \big(D_w\big)^1 w^m = \big(w\,\frac{d^2}{dw^2}+\frac{d}{dw}\big)w^m=m^2\, w^{m-1} $$
$$ \big(D_w\big)^2 w^m = \big(w\,\frac{d^2}{dw^2}+\frac{d}{dw}\big)\, m^2w^{m-1}=(m\,(m-1))^2\, w^{m-2} $$
$$ \big(D_w\big)^j w^m = \Big( j!\binom{m}{j} \Big)^2 \, w^{m-j} $$
Thus it is easy to see that
$$ \exp(c\,D_w) w^j = \sum_{j=0}^m c^j j! \binom{m}{j}^2 w^{m-j} =
m! \sum_{j=0}^m c^{m-j} \binom{m}{j} \frac{w^{j}}{j!} $$
where in the last step a series rearrangement has been performed. Use the definition of the Laguerre polynomials and we have
$$ \exp(c\,D_w) w^j = m!c^m \, L_m(-w/c) .$$ Thus we have a representation of a solution,
$$ \exp(c\,D_w) j_0(\sqrt{w+y}) = \sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)!} \sum_{m=0}^{k} \binom{k}{m} y^{k-m} c^m \, m!\,L_m(z) \, \, , (z=-w/c).$$
This representation can be transformed to a integral. Use the integral that can be found on Wolfram's website
$$ m!L_m(z) = e^z \int_0^\infty e^{-t} t^m J_0(2\sqrt{t\,z}) dt. $$
Pull the integral through the finite sum, then the infinite sum (justifiable by analytic continuation in a region that the integral at the end makes sense).
$$\exp(c\,D_w) j_0(\sqrt{w+y}) =e^z \int_0^\infty dt\,e^{-t} J_0(2\sqrt{t\,z}) \sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)!} \sum_{m=0}^{k} \binom{k}{m} y^{k-m} (ct)^m$$
$$=e^z \int_0^\infty dt\,e^{-t} J_0(2\sqrt{t\,z}) \sum_{k=0}^\infty \frac{(-1)^k}{(2k+1)!} \big(y + ct\big)^k
=e^z \int_0^\infty dt\,e^{-t} J_0(2\sqrt{t\,z}) \frac{\sin{(\sqrt{y+ct})}}{\sqrt{y+ct}} $$
In terms of the original variables this becomes
$$\exp{(i\lambda D_x)}\, j_0(x) = e^{\big((x^2-a^2)/(-4i\lambda)\big)}
\int_0^\infty dt\,e^{-t} J_0\Big(2\sqrt{t\,\frac{(x^2-a^2)}{-4i\lambda}}\Big) \frac{\sin{(\sqrt{a^2+4i\lambda t})}}{\sqrt{a^2+4i\lambda t}}.$$
Not knowing the domains of the parameters, I can't say if the integral is finite; it is probable that it is considering the nice functions in the integrand.
Best Answer
Quoting Wikipedia, when solving the Helmholtz equation in spherical coordinates by separation of variables, the radial equation has the form: $$ x^2 y'' + 2x y' + \left[x^2-n(n+1)\right] y = 0. \tag{1}$$ For $n=0$, it is straightforward to check that $y_0=\operatorname{sinc}(x)$ is a solution.
Suppose now that a solution of $(1)$ has the form $y=(-x)^n f(x)$. Then $f(x)$ satisfies the ODE:
$$ x\, f + 2(1+n)\,f' +x\, f'' = 0.\tag{2} $$ Can you prove now that a solution of $(2)$ is given by $\left(\frac{1}{x}\frac{d}{dx}\right)^n y_0$?