Evaluating sums with spherical Bessel functions

Question

I'm running Mathematica 13.0, student edition, on Windows. In physics, the total scattering cross-section of plane waves off a hard sphere of radius $1$ is

$$\sigma(k) = \frac{4 \pi}{k^2} \sum_{l=0}^\infty (2l+1) \left|\frac{j_l(k)}{h_l(k)} \right|^2$$

where $j_l$ is a spherical Bessel function and $h_l$ is a spherical Hankel function, both of the first kind.

I am hoping to evaluate this in Mathematica and to make a plot of $\sigma(k)$ against $k$ from $.1$ to $20$. However, I'm running into serious difficulties - Mathematica seems to hang, even if I truncate the upper limit of the sum from $\infty$ down to just $30$.

For example, for $k=1$,

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2 NSum[(2 l + 1) Abs[SBJ[l, 1]/SH1[l, 1]]^2, {l, 0, 30}]

hangs for a quite long time before spitting the error "Summand (or its derivative) [complicated expression for summand] is not numerical at point l=15".

Truncating the upper limit to just $10$ seems to run instantaneously, but I expect that to give a poor approximation at larger values of $k$; $k=20$, say.

I was under the impression that spherical Bessel and Hankel functions are actually rather simple - for example, $h_0(k) = -i \frac{e^{ik}}{k}$, $h_1(k) = (-\frac{i}{k^2}-\frac{1}{k})e^{ik}$, and so on. That is, the Hankel functions are just complex exponentials multiplying complex polynomials in $\frac{1}{k}$. I could see trouble brewing for small $k$, but above I was taking $k=1$.

---

What's causing the massive slowdown in calculating the sum, even when truncating the upper bound of the sum to a small value like $l_{max}=30$? How might I speed this up?

Answer 1

I changed Abs[...]^2 to ... Conjugate[...] and NSum to Sum[...]//N

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2   Sum[(2 l + 1)  # Conjugate[#] &[SBJ[l, 1]/SH1[l, 1]] , {l, 0, 30}] // N
(*10.6262 + 4.36409*10^-17 I*)

Answer 2

Abs is no problem, if you avoid interpolation of NSum and calculate NSum at the exact integer l values with the help of NSumTerms .

SBJ = SphericalBesselJ;
SH1 = SphericalHankelH1;
(4 Pi)/1^2 NSum[(2 l + 1) Abs[SBJ[l, 1]/SH1[l, 1]]^2, {l, 0, 30}, 
  NSumTerms -> 31]

(*   10.6262   *)