A square-integrable function on a sphere $f(\hat n)$ can be expanded as a linear combination of spherical harmonics:
$$ f(\hat n) = \sum_{\ell, m} f_\ell^m Y_\ell^m(\hat n),$$
where $\hat n$ is a direction vector and it is usually parameterized by two angles $\theta$ and $\phi$. In the limit of small $\theta$, vector $\hat n$ starts to approach the $\hat z$ direction on the sphere and it is often schematically written that
$$ \sum_{\ell, m} f_\ell^m Y_\ell^m(\hat n) \to \int \frac{d^2\vec l}{(2\pi)^2}~ \tilde f(\vec l) e^{\vec l\cdot \vec\theta},$$
where the right-hand side it the 2D Fourier transform on a plane. I presume that the vector $\vec \theta$ lives now in this 2D plane and is given by the small amplitude $\theta$ and direction determined by $\phi$.
Can somebody show, or point me to, a robust derivation of this limit? I was able to find some intuitive explanations of how these two expressions above should correspond to each other in this limit, but this was not quite satisfactory. Is there a solid straightforward way to show this?