Let $P_n:\mathbb{R}\rightarrow\mathbb{R}$ be the $d$-dimensional Legendre polynomials, that is they are orthonormal polynomials w.r.t the probability measure $\mu_d$ on $[-1,1]$ given by $\mu_d= \frac{\Gamma\left(\frac{d}{2}\right)}{\sqrt{\pi}\Gamma\left(\frac{d-1}{2}\right)}\cdot (1-x^2)^{\frac{d-3}{2}}$. In other words, $\int_{-1}^1P_i(x)P_j(x)\mu_d(x) = \delta_{i,j}$.
Let $v_1,\dots,v_N\in \mathbb{S}^{d-1}$ be some vectors. I need to calculate the following product:
$$\int_{x\sim\mathcal{U}(\mathbb{S}^{d-1})} P_{i_1}(\langle v_1,x\rangle)\cdots P_{i_N}(\langle v_N,x\rangle) dx~,$$
where $\mathcal{U}(\mathbb{S}^{d-1})$ is the uniform distribution over the unit sphere.
I think that I can prove that for two Legendre polynomials, we have $\int_{x\sim\mathcal{U}(\mathbb{S}^{d-1})} P_i(\langle v_i,x\rangle)P_j(\langle v_j,x\rangle)dx = \delta_{i,j}P_i(\langle v_i,v_j\rangle)a_{i,d}$ where $a_{i,d}$ is some normalization term that depends on $i$ and $d$. However, even for three polynomials, I couldn't calculate this integral.
Also, is it possible to bound this integral by some constant independent of the degrees of the polynomials? For example, in the case of a product of two polynomials, I think this is bounded by $1$.
