Sum of Square of the Eigenvalues of Wishart Matrix Let $A\in\mathbb{R}^{m\times d}$ matrix with iid standard normal entries, and $m\geqslant d$, and define $S=A^T A$. 
I want to have a tight upper bound  for $\sum_{k=1}^d \lambda_k^2$, where $\lambda_1,\dots,\lambda_d$ are the eigenvalues of $S$. 
What I tried:


*

*We know that (see e.g. Corollary 5.35 in Vershynin's notes), for $A\in\mathbb{R}^{m\times d}$, for any $t\geqslant 0$, with probability at least $1-2\exp(-\Omega(t^2))$, it holds: 
$$
\sqrt{m}-\sqrt{d}-t \leqslant \sigma_{min}(A)\leqslant \sigma_{max}(A)\leqslant \sqrt{m}+\sqrt{d}+t.
$$
Simply ignoring $\sqrt{d},t$ terms (say I am in the regime $m\gg d,t$), this yields $\lambda_i(A)<m^2$, and thus, the sum above is upper bounded by $m^2d$. 

*We also have the following:
$$
\sum_{k=1}^d (\lambda_k - m) = \sum_{i =1}^m \sum_{j=1}^d (A_{ij}^2-1),
$$
which is sum of sub-exponential random variables, and thus, by a Bernstein-type bound, $\sum_{k=1}^d \lambda_k \leqslant md+\omega(\sqrt{md})$, for some function $\omega(\sqrt{md})$ growing faster than $\sqrt{md}$. 

*The sum above is simply the trace of $S^2=A^TAA^TA$.


I'm new to random matrix business, so any help is greatly appreciated.
 A: I will assume $m=\alpha_d d$ with $\alpha_d\to \alpha \in [1,\infty)$ independent of $d$. The case $\alpha\to\infty$ is actually easier.
Define $Z=d^{-1} m^{-2} \sum_{i=1}^d \lambda_i^2$. Then $Z$ converges a.s. to $\int x^2 d\mu_\alpha(x)$ where $\mu_\alpha$ is the Pastur-Marchenko distribution of parameter
$\lambda=1/\alpha$, see  https://en.wikipedia.org/wiki/Marchenko%E2%80%93Pastur_distribution
A: The problem can be solved using only elementary arguments (i.e without RMT).

Claim. In the limit $d,m \to \infty$ such that $m/d \to \rho \in (0,\infty)$, it holds that
$$
m^{-1}d^{-2}\sum_{i}\lambda_i(S)^2 \overset{a.s}{\to} 1+\rho.
$$

Indeed, one may write
$$
\sum_{i}\lambda_i(S)^2 = \sum_{i,j=1}^m s_{i,j}^2 = \sum_{i,j} (a_i^\top a_j)^2 = \sum_{i=1}^m(\|a_i\|^4 + \sum_{j\ne i}^m (a_i^\top a_j)^2).
\tag{1}
$$
Now, by the Law of Large numbers, it's clear that $d^{-2}\|a_i\|^4 = (d^{-1}\|a_i\|^2)^2 \overset{a.s}{\to} 1^2 = 1$ for all $i \in [m]$. On the other hand, if $i \ne j$, then $d^{-2}(a_i^\top a_j)^2$ has a beta distribution with parameters $\alpha=1$ and $\beta=d-2$ (see this post https://mathoverflow.net/a/227156/78539), and expected value $\alpha / (\alpha + \beta) = 1/(d-1)$. Combining with (1) via another application of LLNs (to handle the second term) then completes the proof after noting that $(m-1) / (d-1) \to \rho$.

Question. Comparing with the accepted answer, is true that $\langle \lambda^2\rangle_{MP(1/\rho)} = 1 + \rho$ ?

A: The joint probability distribution of the eigenvalues of $S$ is proportional to
$$ \rho(S)=e^{-{\rm Tr}(S)}\det(S)^{a/2}|\Delta(S)|,$$
where $a=m-d-1$ and $\Delta(S)$ is the Vandermonde. The average value of any symmetric function $f$ of the eigenvalues can be computed exactly by writing $f$ as a linear combination of Zonal polynomials, $f(S)=\sum_\lambda c_\lambda Z_\lambda(S)$, and then using the explicit result of the Selberg-like integral
$$ \int_0^\infty Z_\lambda(S)\rho(S)dS,$$
which can be found e.g. in The Importance of the Selberg Integral, by P. Forrester and O. Warnaar.
