Approximation of Spherical Bessel function I am currently studying the CMB power spectrum from a numerical approach (easier than the analytical approach). In a Mathematica notebook that I am following, they work with spherical Bessel functions in order to free stream the multipole solution of the fluid equations in Fourier space. I understand the analytical implementation of the Bessel functions in the formula, but in the Mathematica code, they approximate these functions in a way which I have not been able to derive for myself or find online. The approximation of the Bessel function is
\begin{equation}
l^2 j_l^2(xl) = \frac{1}{2x\sqrt{x^2-1}}.
\end{equation}
I also have to use the derivative of the Bessel function which they approximate as
\begin{equation}
l^2 j_l'^2(xl) = \frac{1}{2x\sqrt{x^2-1}}\frac{x^2-1}{x^2}.
\end{equation}
This approximations are done in the limit where both $x$ and $xl$ are large.
I would very much appreciate anyone bringing some insights the derive these approximations!
 A: This is too long for a comment so I wrote this answer.
I looked in the obvious place,
G. N. Watson,
"Treatise on the Theory of Bessel Functions",
(Cambridge University Press,Cambridge,1980), second edition,
in section 8.12 he gives an expansion first derived by
Meissel for large order and $x$ times the order large.
Watson then discusses the stationary phase approximation in section 8.2.
Watson gives the Meissel series, where he says
this dominant term had been derived by L. Lorenz earlier,
\begin{equation}
J_\nu(x) \simeq \sqrt{\frac{2}{\pi\sqrt{x^2-\nu^2}}}
\cos\left (Q_\nu(x)-\frac{1}{4}\pi\right)
\end{equation}
\begin{equation}
Q_\nu(x) = \sqrt{x^2-\nu^2}-\frac{1}{2}\nu\pi+\nu\arcsin(\nu/x)
\end{equation}
If I substitute, $j_\ell (x) = \sqrt{\frac{\pi}{2x}} J_{\ell+1/2}(x)$,
\begin{equation}
j_\ell(\ell x) \simeq 
\sqrt{\frac{\pi}{2x\ell}}
\sqrt{\frac{2}{\pi\sqrt{\ell^2 x^2-(\ell+\frac{1}{2})^2}}}
\cos\left (Q_{\ell+1/2}(\ell x)-\frac{1}{4}\pi\right) \,.
\end{equation}
Simplifying, squaring, and multiplying by $\ell^2$,
\begin{equation}
\ell^2 j_\ell(\ell x) \simeq \frac{1}{x\sqrt{x^2-\left (\frac{2\ell+1}{2\ell}
\right)^2}}
\cos^2\left (Q_{\ell+1/2}(\ell x)-\frac{1}{4}\pi\right ) \,.
\end{equation}
If, as suggested by Emilio Pisanty in the first comment,
that the cosine squared is approximated by its average, $\frac{1}{2}$,
in the sums or integrals you are doing,
and you approximate
$\left (\frac{2\ell+1}{2\ell} \right)^2 \simeq 1$,
you get your result.
In the stationary phase approximation, it looks to me like the two
stationary phase points give integrals that give the
$\frac{1}{x\sqrt{x^2-\left (\frac{2\ell+1}{2\ell} \right)^2}}$ factor,
and their phases give the cosine term, but as I said, I'm too lazy to
spend the time to carefully calculate and check that term.
