Let $\mathcal P$ the set of all probability distributions on the sphere $\Bbb S^n\subset \Bbb R^{n+1}$ that have smooth densities with respect to the Lebesgue surface measure. For all $P\in\mathcal P$ let $$ \rho(P)=\max_{\lambda \in \Bbb S^n}\lambda^T Cov(P) \lambda $$ the largest eigenvalue of the covariance matrix $$ Cov(P)=\left(\int_{\Bbb S^n} x_ix_j P(dx)-\int_{\Bbb S^n} x_i P(dx)\int_{\Bbb S^n} x_j P(dx)\right)_{i,j=1,\ldots,n+1}. $$ It seems intuitively clear to me that if $\mathcal U$ denotes the uniform distribution on $\Bbb S^n$, then we should have $$ \sup_{P\in\mathcal P}\rho(P)=\rho(\mathcal U). $$ I am looking for a reference for this statement.
I thought this might be found in Directional Statistics by Kanti Mardia and Peter Jupp, but so far I've had no luck finding it there.
PS: I am aware that $Cov(\mathcal U)=(n+1)^{-1}\cdot 1_{(n+1)\times(n+1)}$ and hence $\rho(\mathcal U)=(n+1)^{-1}$. That is not part of my question.
The claim $$ \sup_{P\in\mathcal P}\rho(P)=\rho(\mathcal U) $$ is false.
Let $p=(1,0,\ldots,0)\in\Bbb S^n$ and $$ P=\frac 12 (\delta_p+\delta_{-p}). $$ Then $$ \int_{\Bbb S^n} x_ix_j P(dx)-\int_{\Bbb S^n} x_i P(dx)\int_{\Bbb S^n} x_j P(dx)= \begin{cases} 1, & i=j=1, \\ 0, & \text{else,}\end{cases} $$ and hence $$ \rho(P)=1. $$ If $(P_m)_{m\in\Bbb N}$ is a sequence of smooth approximations of $P$, i.e. $P_m\to P$ for $m\to\infty$ in the sense of weak convergence of probability measures, then (by compactness of $\Bbb S^n$) we have $$ Cov(P_m)\to Cov(P), $$ and therefore $$ \rho(P_m)\to\rho(P)=1. $$ Hence $$ \sup_{P\in\mathcal P}\rho(P)\ge 1>\frac{1}{n+1}=\rho(\mathcal U). $$ (Since $1$ is trivially also an upper bound, we have $\sup_{P\in\mathcal P}\rho(P)=1$.)