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}

I would very much appreciate anyone bringing some insights the derive these approximations!