Let $f:P\in M_n\rightarrow P^2-P$ and $V$ be the variety $f^{-1}(0)$.

The derivative of $f$ in an idempotent matrix $P$ is:

$Df_P:H\in M_n\rightarrow PH+HP-H$, that is, $Df_P=P\otimes I_n+I_n\otimes P^T-I_n\otimes I_n$

when we stack the matrices row by row.
$Df_P$ is the sum of three matrice that pairwise commute; moreover, since any idempotent $P$ is diagonalizable, these matrices are diagonalizable and, consequently, are simultaneously diagonalizable. Let $spectrum(P)=(\lambda_i)_i$; the previous result implies that $spectrum(Df_P)=(\lambda_i+\lambda_j-1)_{i,j}$.

The dimension of $V$ is the max of the $dim(\ker(Df_P))$ when $P$ varies through the idempotent matrices, that is, the max of the number of zero eigenvalues of $\ker(Df_P)$ when $P$ varies.

Since $\lambda_i+\lambda_j-1=0$ iff we consider a couple $0,1$ of eigenvalues, we can convince ourselves that this max is obtained when the $0$'s and the $1$'s are about equal in number in $spectrum(Df_P)$.

Case 1. $n=2p$ is even. We consider some $P$ s.t. $spectrum(P)$ has $p\times 0$ and $p\times 1$. The number of couples $(0,1)$ is $p^2$ and the number of couples $(1,0)$ is also $p^2$. Then the maximal local dimension is $2p^2$.

Case 2. $n=2p+1$ is odd. We consider some $P$ s.t. $spectrum(P)$ has $p\times 0$ and $p+1\times 1$. The number of couples $(0,1)$ is $p^2+p$ and the number of couples $(1,0)$ is also $p^2+p$. Then the maximal local dimension is $2p^2+2p$.
In other words, $dim(V)$ is the integer part of $n^2/2$.

EDIT. Answer to @Jose Brox . We assume that the underlying field is $\mathbb{C}$. $V$ is not a pure variety -its local dimension is locally constant but not globally- but rather an algebraic set. Let $V_r$ be the set of idempotent matrices of rank $r$; then $V$ has $n+1$ connected components $V_0,\cdots,V_n$, each being a pure variety of dimension $2r(n-r)$ (cf. the above proof); in particular $V_0=\{0_n\},V_n=\{I_n\}$.

Another way to see that $dim(V_r)=2r(n-r)$ is as follows: to define an idempotent $P$ of rank $r$ is equivalent to choose two complementary vector spaces $W_0,W_1$ of dimension $n-r,r$: $W_0=\ker(P),W_1=im(P)$. Note that $W_0\in G_{n-r,n},W_1\in G_{r,n}$ where $G_{n-r,n},G_{r,n}$ are Grassmannian varieties of dimension$(n-r)r,r(n-r)$. Then $dim(V_r)=dim(G_{n-r,n}\times G_{r,n})=2r(n-r)$; indeed, if you are very unlucky, $W_0,W_1$ intersect; such a choice must satisfy a system of algebraic relations; then your choice is contained in a Zariski closed set of dimension $<2r(n-r)$.

Now, $dim(V_r)=2r(n-r)$ means that, for every $P_0\in V_r$, there are neighborhoods $R$ of $P_0$, $S$ of $0$ in $\mathbb{C}^{2r(n-r)}$
and a diffeomorphism $\Pi_{P_0}:R\rightarrow S$. Finally, $(\Pi_{P_0})^{-1}$ locally realizes a parametrization of $V_r$. That was the good news; the bad one, is that this result gives no explicit parameterization.

For example, consider the simple case when $n=2,r=1$; let $P_0=\begin{pmatrix}a_0&b_0\\c_0&1-a_0\end{pmatrix}\in V_1$ where $a_0(1-a_0)=b_0c_0$. If $b_0\not= 0$, then a local param. is $(a,b)\in \mathbb{C}\times \mathbb{C}^*\rightarrow \begin{pmatrix}a&b\\a(1-a)/b&1-a\end{pmatrix}$.
We similarly treat the case $c_0\not=0$; yet, what is an admissible param. when $P_0=diag(0,1)$ or $P_0=diag(1,0)$ ?

Of course, $Q\in GL_n\rightarrow Qdiag(0_{n-r},I_r)Q^{-1}\in V_r$ is not an admissible param.!!