Suppose $A$ is a square matrix whose eigenvalues are $\lambda_1, \lambda_2, \dots, \lambda_n$. The equation
$$\mbox{tr} \left( A^2 \right) = \mbox{tr} (A)^2$$
implies that $$\sum_{i\not=j}\lambda_i\lambda_j=0$$
What can we learn about the matrix $A$ from the equation above? What kind of matrix has this property?
This is really more a question about collections of numbers than about linear algebra.
As you know, $\operatorname{tr}A=\sum_i\lambda_i$ and $\operatorname{tr}A^2=\sum_i\lambda_i^2$. Let $\mu=\frac1n\sum_i\lambda_i$ and $\sigma^2=\left(\frac1n\sum_i\lambda_i^2\right)-\mu^2$ be the mean and variance of the eigenvalues. (The latter is actually the pseudo-variance if the eigenvalues are complex.) Then, $\operatorname{tr}A=n\mu$ and $\operatorname{tr}A^2=n(\mu^2+\sigma^2)$. So $(\operatorname{tr}A)^2=\operatorname{tr}A^2$ if and only if $n\mu^2=\mu^2+\sigma^2$, i.e. $$\mu=\pm\frac\sigma{\sqrt{n-1}}.$$
To me there does not seem to be anything very special about this property. For any matrix $A$ you can add a scalar multiple of the identity to obtain a matrix which satisfies it, i.e. $\operatorname{tr}^2(A+kI)=\operatorname{tr}(A+kI)^2$ for $k=\pm\sigma/\sqrt{n-1}-\mu$.