Probabilistic lower and upper-bounds for a certain random quartic form involving gaussian random matrices Let $d,m \to \infty$ (integers) with $m/d \to \rho \in (0, \infty)$. Let $C$ be a $d \times d$ psd matrix with $trace(C)=\mathcal O(1)$, and let $w_1,\ldots,w_m$ be iid uniformly distributed on the unit-sphere in $\mathbb R^d$. Consider the quartic form
$$
F := \frac{1}{m}\sum_{j,\ell=1}^m (w_j^\top w_\ell)(w_j^\top C w_\ell).
$$

Question. What are good probabilistic lower and upper-bounds for $F$ only in terms of $\rho$ and the eigenvalues of $C$ ?

For example, the solution for the case where $C$ is diagonal will already be very helpful.
Isotropic example
Thanks to this post https://mathoverflow.net/a/334219/78539, we know that if $C = (1/d) I_d$, then $F = m^{-1}\|WW^\top\|_F^2 = m^{-1}\sum_{j}\lambda_j(W W^\top)^2\overset{a.s}{\to} \langle \lambda^2\rangle_{\text{MP}(1/\rho)}$ (if I haven't made some scaling errors), where $\text{MP}(\gamma)$ is the Marchenko-Pastur law with parameter $\gamma$.
 A: Assume iid $N(0,1)$ entries, assume $C$ diagonal, and focus first on the non-diagonal terms:
$G=\sum_j \sum_{l\ne j} w_j^Tw_l w_j^TCw_l
= \sum_{j\ne l, ik} w_{ji}w_{li} c_i w_{jk} w_{lk}$.
Write this quantity as
$$
\begin{split}
G=\sum_{j\ne l, i\ne k} w_{ji}w_{li} c_i w_{jk} w_{lk}
&+
\sum_{j\ne l, i=k} (w_{ij}^2-1)(w_{lj}^2 -1)c_i
\\
&+(w_{ij}^2-1)c_i + 
(w_{lj}^2 -1)c_i
+c_i
\end{split}
$$
This is a decomposition in uncorrelated polynomials (any two terms are uncorrelated), so that the second moment is
$$
E[(G-m(m-1)trace[C])^2]=\sum_{j\ne l, i \ne k} c_i^2
+
\sum_{j\ne l, i}(E[(Z^2-1)^2]^2 + 2 E[(Z^2-1)^2])c_i^2.
$$
$$= m(m-1)\|C\|_F^2((d-1)+E[(Z^2-1)^2]^2 + 2E[(Z^2-1)^2]).$$
The dominant term is of order $m^2d \|C\|_F^2$, while the mean is $m(m-1)trace[C]$. Hence $G/E[G]-1$ converges to 0 in probability (or in L2) provided that $E[G]^2 \gg Var[G]$, that is,
$$
m^2 trace[C] \gg \|C\|_F m \sqrt{d}.
$$
For the diagonal terms, we have
$\sum_j d w_j^TCw_j + \sum_{j} (\|w_j\|^2-d)w_j^TCw_j$.
The second term is negligible compared to the first one if you use $\chi^2$ concentration (e..g, Bernstein inequality) for $\|w_j\|^2-d$, while the first term has mean $md trace[C]$ and variance $2md^2\|C\|_F^2$.
Again, the mean dominates the standard deviation if and only if
$$
m d ~trace[C] \gg \sqrt m d \|C\|_F.$$ This is equivalent to the condition on the non-diagonal terms if $m\asymp d$.
Edit: since $\|C\|_F^2 \le trace[C]^2$ for $C$ psd, these conditions are always satisfied.
