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 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!