Integration on the unit sphere I have an integral on the unit sphere as follows. 
$$I(\mathbf{s}_1, \mathbf{s}_2) = \int_{\mathbb{S}^2} f(\mathbf{x} \cdot \mathbf{s}_1)f(\mathbf{x}\cdot\mathbf{s}_2)d\mathbf{x} $$
where the integral is on the whole unit sphere $\mathbb{S}^2$. $\mathbf{s}_1$ and $\mathbf{s}_2$ are also unit-length vectors on the unit sphere. For example, $\mathbf{s}_1$ is the vector from the origin to a point on the sphere $\vec{O s_1}$.
Intuitively, since the integral is over the whole sphere, $I(\mathbf{s}_1, \mathbf{s}_2)$ should only depend on the relative angle between $\mathbf{s}_1$ and $\mathbf{s}_2$. This is actually done here. 
Now I want to go further and do a Taylor expansion in terms of the relative spherical angle between $\mathbf{s}_1$ and $\mathbf{s}_2$, say $\theta$, so that $\theta = \angle (\mathbf{s}_1, \mathbf{s}_2)$. When $\mathbf{s}_1$ is close to $\mathbf{s}_2$, $\theta$ will be small and I want to do an expansion in terms of $\theta$. The problem is that I don't know how to explicit write $I(\mathbf{s}_1, \mathbf{s}_2)$ as a function of $\theta$, so I am stuck. 
I need help on write this $I(\mathbf{s}_1, \mathbf{s}_2)$ as an Taylor expansion in terms of $\theta$, the relative angle between vectors $\mathbf{s}_1$ and $\mathbf{s}_2$. Any comments will be greatly appreciated. 
 A: Here is a way how to explicitly write $I$ in terms of the inner product.
For such integrals it is usually most convenient to use the addition theorem of the spherical harmonics
$$P_l( \mathbf{x}\cdot\mathbf{y} ) = \frac{4\pi}{2l+1}\sum_{m=-l}^l Y_{l m}^*(\theta_\mathbf{x},\varphi_\mathbf{x}) \, Y_{l m}(\theta_\mathbf{y},\varphi_\mathbf{y})$$
valid for two vectors $\mathbf{x}, \mathbf{y}$ on the unit sphere; here $(\theta_\mathbf{x},\phi_\mathbf{x})$  are the spherical coordinates of $\mathbf{x}$, respectively.
Due to the completeness of the Legendre polynomials, we can introduce
$$f_l =\frac{2}{2l+1}\int_{-1}^1 f(x) P_l(x) dx$$
such that $$f(x) = \sum_{l=0}^\infty f_l P_l(x).$$
Due to the integral over the unit sphere, we can choose $\mathbf{s}_1$ to be along the $z$ axis and $\mathbf{s}_2$ in the $xz$ plane. Then
$$\begin{align}I(\mathbf{s}_1, \mathbf{s}_2) &= \int f(\cos \theta_\mathbf{x})f(\mathbf{x}\cdot\mathbf{s}_2)d\Omega_{\mathbf{x}}=
\sum_{l l'} f_l f_{l'} \int P_l(\cos \theta_\mathbf{x}) P_{l'}(\mathbf{x}\cdot\mathbf{s}_2)d\Omega_{\mathbf{x}}\\
&=\sum_{l l'} f_l f_{l'} \sqrt{\frac{4\pi}{2l'+1}}\sum_{m'=-l'}^{l'}Y_{l'm'}(\theta_{\mathbf{s}_2},0)\underbrace{\int Y_{lm}(\theta_\mathbf{x},0) Y^*_{l'm'}(\theta_\mathbf{x},\varphi_\mathbf{x})d\Omega_{\mathbf{x}}}_{\delta_{ll'}\delta_{m'0}}\\
&=\sum_l f_l^2 P_l(\cos \theta_{\mathbf{s}_2})\\
&=\sum_l f_l^2 P_l(\mathbf{s}_1\cdot\mathbf{s}_2),\end{align}$$
where in the last step we have reintroduced a coordinate independent way of writing the expression.
