The function $f$ of two 3-dimensional unit vectors $\hat r_1$ and $\hat r_2$ only depends on the angle from one vector to another, so we can write $f(\hat r_1, \hat r_2)=g(\theta)$, where $\theta$ is the angle from $\hat r_1$ to $\hat r_2$. We can further decompose the function $g(\theta)$ into angular harmonics as \begin{equation} g(\theta)=\sum_mg_me^{im\theta} \end{equation} Notice we do not assume $g(\theta)=g(-\theta)$. Here the $-\theta$ can be interpreted as the angle from $\hat r_2$ to $\hat r_1$.
Now consider the integral \begin{equation} I(\hat r_1, \hat r_2)=\int d\hat r'f(\hat r_1, \hat r')f(\hat r', \hat r_2) \end{equation} where the integral is over all 3-dimensional unit vectors, $\hat r'$. I would expect this integral to still only depend on the angle from $\hat r_1$ to $\hat r_2$, so it also allows a decomposition into angular harmonics: \begin{equation} I(\hat r_1, \hat r_2)=\sum_mI_me^{im\theta} \end{equation} where $\theta$ is still the angle from $\hat r_1$ to $\hat r_2$.
The question is how to write $I_m$ in terms of $g_m$.