# Integrating a function with the spherical Bessel function of the first kind

I want to integrate a function with the spherical Bessel function of the first kind from $0$ to $\infty$. I tried in many different ways, but didn't get any appropriate answer; if anyone can help me I am very thankful.

psi[n_?NumericQ, x_?NumericQ] := x SphericalBesselJ[n, x];

Result = NIntegrate[psi[n, x], {x, 0, ∞}, MaxRecursion -> 20,
PrecisionGoal -> 10, AccuracyGoal -> 15]

• Note that $x j_n(x) \sim\cos \left(\frac{\pi}{2}(n+1)-x\right)$ as $x\rightarrow\infty$, so these integrals don't converge for any real values of $n$. If you plot the integrand, you'll see why. – Peter Richter Jan 1 '18 at 4:15
• If you multiply the integrand by Exp[-a x], integrate, and take the limit as a -> 0, you can obtain Binomial[n - 1/2, n - 1]/2 for odd n and 1/Binomial[n - 3/2, n - 1] for even n (assuming n is a nonnegative integer). Caveat: I did this with Integrate, which gave me an answer in terms of Hypergeometric2F1Regularized, but I haven't checked for correctness/reasonableness. – Michael E2 Jan 1 '18 at 7:29

Let me formalize Michael's proposal in the comments. What he is suggesting is that you take the Laplace transform of your integrand, and then take the limit of that transform at $0$ (i.e., evaluating the divergent oscillatory integral in the Abel sense).

Thus:

LaplaceTransform[x SphericalBesselJ[n, x] // FunctionExpand, x, h] // FullSimplify
2^(-1 - n) h^(-2 - n) Sqrt[π] Gamma[2 + n]
Hypergeometric2F1Regularized[(2 + n)/2, (3 + n)/2, 3/2 + n, -1/h^2]


Unfortunately, attempting to take the limit of this as $h\to0$ does not yield anything useful, so we take a circuitous route.

The idea is to convert the expression to a Meijer $G$-function and hope the result is a bit more amenable.

Using MeijerGReduce[] on the hypergeometric function, we have

MeijerGReduce[Hypergeometric2F1Regularized[(2 + n)/2, (3 + n)/2, 3/2 + n, -z], z]
MeijerG[{{-n/2, -1/2 - n/2}, {}}, {{0}, {-1/2 - n}}, z]/
(Gamma[(2 + n)/2] Gamma[(3 + n)/2])


(I have suppressed the inactivated form for readability in the previous and in all subsequent outputs.)

Combining the prefactors of this result with the prefactors of the bigger expression, we have

2^(-1 - n) h^(-2 - n) Sqrt[π] Gamma[2 + n]/(Gamma[(2 + n)/2] Gamma[(3 + n)/2])
// FullSimplify
h^(-2 - n)


so our function under consideration is now

h^(-2 - n) MeijerG[{{-n/2, -1/2 - n/2}, {}}, {{0}, {-1/2 - n}}, 1/h^2]


We now employ a few identities satisfied by the $G$ function. First, we use the "absorption identity", which yields

Inactive[MeijerG][{{-n/2, -1/2 - n/2}, {}} + (n + 2)/2,
{{0}, {-1/2 - n}} + (n + 2)/2, 1/h^2] // Simplify
MeijerG[{{1, 1/2}, {}}, {{1 + n/2}, {1/2 - n/2}}, 1/h^2]


Finally, employing the "flip identity" gives

Inactive[MeijerG][1 - {{1 + n/2}, {1/2 - n/2}}, 1 - {{1, 1/2}, {}}, h^2] // Simplify
MeijerG[{{-n/2}, {(1 + n)/2}}, {{0, 1/2}, {}}, h^2]


And now, we let h -> 0:

 Activate[% /. h -> 0]
(Sqrt[π] Gamma[1 + n/2])/Gamma[(1 + n)/2]


or alternatively, n/2 Beta[1/2, n/2].

You can do this to see what is going on. Mathematica can get the indefinite integral.

intp = Integrate[x SphericalBesselJ[n, x], x]
(*Sqrt[Pi] 2^(-n - 2) x^(n + 2)Gamma[n/2 + 1]
HypergeometricPFQRegularized[{n/2 + 1}, {n + 3/2, n/2 + 2}, -(x^2/4)]*)


I know you don't specify n to be an integer, but we can look at the first 5 integers for n, just to get a look.

tab = Table[intp, {n, 5}];


At x = Infinity

Limit[tab, x -> Infinity] // N
(*{Interval[{0.570796,2.5708}],Interval[{1.,3.}],Interval[{1.35619,3.35619}],Interval[{1.66667,3.66667}],Interval[{1.94524,3.94524}]}*)


The oscillation of the trig functions give a value range at Infinity.

At x = 0, Limit[tab, x -> 0] yields 0 for the 5 cases.

In many cases Mathematica fails to produce a result with NIntegrate, but can produce results with Integrate, so it is wise to check out the other possibility when one fails.

Integral regularization in simpler way:

f=Integrate[ x^-e*SphericalBesselJ[n, x], {x, 0, ∞},
Assumptions -> {e >= 0, n > 0, n ∈ Integers}][[1]]

(*(2^(-1 - e) Sqrt[π] Gamma[1/2 (1 - e + n)])/Gamma[1/2 (2 + e + n)]*)

Limit[f, e -> -1]

(* (Sqrt[π] Gamma[(2 + n)/2])/Gamma[(1 + n)/2] *)


$$\frac{\sqrt{\pi } \Gamma \left(\frac{n+2}{2}\right)}{\Gamma \left(\frac{n+1}{2}\right)}$$

• Oh, you took the Borel route. :D – J. M.'s torpor Sep 30 '18 at 14:46