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?
The particular case $n=2$ has a nice specific proof.
As the characteristic equation of $A$ can be written under the form
$$\lambda^2- \operatorname{tr}(A)\lambda +\det(A)=0,$$
Cayley–Hamilton theorem gives
$$A^2-\operatorname{tr}(A) A +\det(A)I=0.$$
Taking the trace (linear operator) :
$$\operatorname{tr}(A^2)-\operatorname{tr}(A)^2=-2\det(A)$$
Therefore the issue is simply to look for matrices with zero determinant, i.e., matrices such that $\lambda_1\lambda_2=0$.
See this excellent document by Darij Grinberg.
See as well : (https://math.stackexchange.com/q/1267700.)
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$.
$\textbf{Part 1.}$ Let $Z=\{X\in M_n(\mathbb{C});\;(*)\;tr(X^2)=(tr(X))^2\}$; since the studied relation $(*)$ is homogeneous, $Z$ is a complex algebraic cone of dimension $n^2-1$ (that is, $n^2-1$ independent complex parameters).
$\textbf{Part 2.}$ The relation $(*)$ can be written $tr(X^2-tr(X)X)=0$, that is equivalent to
there are $U,V\in M_n(\mathbb{C})$ s.t. $(**)$ $X^2-tr(X)X=UV-VU$.
Conversely, if you want to construct such matrices $X$ without using its eigenvalues, then you can do as follows -when $n\geq 3$ (the case $n=2$ has been solved by Jean Marie)-
i) Randomly choose $U,V$
ii) Solve the equation $(**)$. In general (generic choice of $U,V$) , this equation has $2,8,22,52$ solutions (pairwise opposite) when $n=3,4,5,6$.
Essentially, we have the following equation in $\mathrm x \in \mathbb C^n$
$$\mathrm x^\top \left( \mathrm I_n - 1_n 1_n^\top \right) \mathrm x = 0$$
Note that matrix $\mathrm I_n - 1_n 1_n^\top$ has eigenvalue $1$ with multiplicity $n-1$ and eigenvalue $1-n$ with multiplicity $1$. Hence, the equation above can be rewritten in the eigenvector basis as follows
$$y_1^2 + y_2^2 + \cdots + y_{n-1}^2 - (n-1) y_n^2 = 0$$
Note that $\mathrm y = 1_n$ is a solution.
OP's condition $${\rm tr}(A^2)~=~{\rm tr}(A)^2\tag{1}$$ always holds for $n=1$.
Let us assume $n\geq 2$ from now on. OP's condition (1) is then equivalent to the vanishing of the coefficient in front of the third-highest power $\lambda^{n-2}$ in the characteristic polynomial $$\begin{align}p_{A}(\lambda)~:=&~\det(A-\lambda \mathbb{1})\cr ~=&~(-\lambda)^n~+~ {\rm tr}(A)(-\lambda)^{n-1}~+~\frac{{\rm tr}(A)^2-{\rm tr}(A^2)}{2} (-\lambda)^{n-2}~+~\ldots .\end{align}\tag{2}$$ For the case $n=2$ this agrees with Jean Marie's answer.
