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 on the mathematica code, they approximate these functions in a way which I have not been able to derive 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 the approximate in the following way,
\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!