I'm trying to understand a statement from a paper I'm reading. I'll first put forward the statement as I understand it, and then provide the context from the paper.
Statement
Let $G$ be a Lie group that is connected, simply connected, and finite dimensional. Furthermore, let $T$ be a unitary irreducible representation of $G$ acting on a Hilbert space $\mathcal{H}$. Let $|\psi_0\rangle$ be an element of $\mathcal{H}$ that we take as a fiducial state, and let $H\subset G$ be the subgroup of $G$ consisting of all elements $h$ that act trivially on $|\psi_0\rangle$, leaving it invariant up to a phase
\begin{align} T(h)|\psi_0\rangle&=e^{i\phi(h)}|\psi_0\rangle & h&\in H. \end{align}
Finally, let $X=G/H$ be the coset space which one can interpret as the orbit of $|\psi_0\rangle$ under the action of $G$ modulo phase. $X$ has an invariant measure $d\mu$, and the harmonic functions $Y_\nu$ (eigenfunctions of the Laplace-Beltrami operator on $X$) form a complete orthonormal basis for the Hilbert space $L^2(X,\mu)$.
The statement that I am having difficulty convincing myself of is that if I have a function $f:X\times X\to\mathbb{C}$ that is symmetric and positive, i.e. $f(\Omega,\Omega^\prime)=f(\Omega^\prime,\Omega)\geq0$, its expansion in terms of the harmonic functions must be of the form
\begin{align} f(\Omega,\Omega^\prime)&=\sum_\nu\tau_\nu Y_\nu^*(\Omega)Y_\nu(\Omega^\prime) \end{align}
where $0\leq\tau_\nu\in\mathbb{R}$.
This makes some intuitive sense to me, since a positive symmetric function seems like it ought to be diagonalizable in some basis. I'm curious to know why the harmonic functions are always the correct basis (why couldn't the eigenvectors be linear combinations of basis vectors for different irreducible representations?). I'm also curious how important the positivity of $f$ is for this statement. Does a similar expansion hold for arbitrary symmetric $f$?
I'm a physicist, and do not have much formal training in representation theory or harmonic analysis, so it could be that I'm missing something trivial, but I'd be much obliged to anyone who could point that trivial thing out to me!
Context
In case I've missed out on some essential detail for my question above, the context is the paper Phase-space formulation of quantum mechanics and quantum state reconstruction for physical systems with Lie-group symmetries, and I'm attempting to understand why they can write down (3.6), which specifically is
\begin{align} \big|\langle\Omega|\Omega^\prime\rangle\big|^2&=\sum_\nu\tau_\nu Y_\nu^*(\Omega)Y_\nu(\Omega^\prime) \\ &=\sum_\nu\tau_\nu Y_\nu^*(\Omega^\prime)Y_\nu(\Omega) \end{align}
The space $X$ is supposed to be a manifold of coherent states, which are the orbit of the vacuum state $|\psi_0\rangle$ under the irreducible unitary action of $G$. They go on to further constrain the coefficients $\tau_\nu$ using the invariance of the inner product under the group action, but I find myself confused by their starting point.
Work so far
I've tried to make sense of the situation by considering $\{Y_\nu^*(\Omega)Y_{\nu^\prime}(\Omega^\prime)\}_{\nu,\nu^\prime}$ as an orthonormal basis for $L^2(X\times X)$, and assuming $f(\Omega,\Omega^\prime)=f^*(\Omega^\prime,\Omega)$. This gives me
\begin{align} \sum_{\nu,\nu^\prime}\hat{f}_{\nu,\nu^\prime}Y_\nu^*(\Omega) Y_{\nu^\prime}(\Omega^\prime) &=f(\Omega,\Omega^\prime) \\ &=f^*(\Omega^\prime,\Omega) \\ &=\sum_{\nu,\nu^\prime}\hat{f}^*_{\nu,\nu^\prime} Y_\nu(\Omega^\prime)Y^*_{\nu^\prime}(\Omega) \\ &=\sum_{\nu,\nu^\prime}\hat{f}^*_{\nu^\prime,\nu} Y^*_\nu(\Omega)Y_{\nu^\prime}(\Omega^\prime) \end{align}
This means $\hat{f}_{\nu,\nu^\prime}=\hat{f}^*_{\nu^\prime,\nu}$, or if I define some matrix $[F]_{\nu,\nu^\prime}:=\hat{f}_{\nu,\nu^\prime}$, then $F=F^\dagger$.
If I play a bit fast and loose with the fact that $F$ is actually infinite-dimensional (I did mention that I'm a physicist), then I can diagonalize this matrix with a unitary $U$, so $F=U^\dagger DU$, $[D]_{\mu,\mu^\prime}=d_\mu\delta_{\mu,\mu^\prime}$. This lets me write $\hat{f}_{\nu,\nu^\prime}=u^*_{\mu,\nu}d_\mu u_{\mu,\nu^\prime}$ and
\begin{align} f(\Omega,\Omega^\prime) &=\sum_\mu d_\mu\left(\sum_\nu u^*_{\mu,\nu}Y^*_\nu(\Omega)\right) \left(\sum_{\nu^\prime} u_{\mu,\nu^\prime} Y_{\nu^\prime}(\Omega^\prime)\right) \\ &=\sum_\mu d_\mu V^*_\mu(\Omega)V_\mu(\Omega^\prime) \end{align}
where $\{V_\mu(\Omega)=\sum_\nu u_{\mu,\nu}Y_\nu(\Omega)\}_\mu$ is a new orthonormal basis for $L^2(X)$ (it is the eigenbasis for $f$). In this formulation it seems the important point above is not purely symmetry, but rather hermiticity (which seems rather obvious, in hindsight). I still am left with the question of why the $V_\mu$s don't mix $Y_\nu$s from different irreps (equivalently, why $U$ is block diagonal in the harmonic function labels). I have a pretty good feeling that I can show this must be the case when I invoke the covariance property of the inner product, but the paper left me with the impression that all that was needed was symmetry and positivity.