Unit Matrices for $M_n(M_n(\mathbb{R}))$ In the case of $M_n(\mathbb{R})$ we can define unit matrices $E_{ij}$ in the following way,
$$E_{ij}f_l = \delta_{lj}f_i$$
where the $\{f_i\}$ is a basis for our space.  One consequence of this definition is that the $(i,j)$ element of $E_{ij}$ is $1$ and everything else is $0$; additionally
$$E_{ij}E_{lk} = \delta_{jl} E_{ik}$$
.... okay, anyway I'm not sure how unit matrices look in $M_n(M_n):=M_n(M_n(\mathbb{R}))$.
And how do the above properties carry over?
Also, I don't understand why the following line in a proof is true (as per the above):
Consider $P=(E_{ij}) \in M_n(M_n)$.  Then $P=P^\dagger$, and 
$$P^2 =  (\sum_{k=1}^n E_{ik}E_{kj})(nE_{ij}) = nP$$
 A: I see that you have tagged your question with quantum computation. So perhaps you may find it more intuitive if we work in Dirac notation. That is we let $\{| j \rangle\}$ be an orthonormal basis for the underlying vector space (so that $\langle i|j\rangle =\delta_{i,j}$).  Then what you call the unit matrices are equivalent defined as
$$
E_{i,j} = |i\rangle\!\langle j|.
$$
That is, they have a 1 in the entry at the $i$'th row and $j$'th column and zeros everywhere else.  Now the set of matrices $\{|E_{i,j}\rangle\!\rangle \}$ form an orthonormal basis for the vector space of matrices, where the inner product between two matrices $A$ and $B$ is $$\langle\!\langle A|B\rangle\!\rangle = \operatorname{Tr}(A^\dagger B).
$$
We can now play the same trick.  Define the "super" unit matrices as 
$$
F_{i,j;k,l} = |E_{i,j}\rangle\!\rangle\!\langle\!\langle E_{k,l}|.
$$
Then, the $F$ matrices satisfy the same properties as the $E$ matrices. 
From your last sentence, however, it does not seem like you are interested in these unit matrices.  Your final equation suggests the matrix you define as $P$ is simply the matrix of 1's in every entry.  This can be written
$$
P = \sum_{i,j} E_{i,j}.
$$
Then direct calculation leads to
$$
P^\dagger = \sum_{i,j} E_{i,j}^\dagger = \sum_{i,j} E_{j,i} = P,
$$
and
$$
P^2 = \sum_{i,j,k,l} E_{i,j}E_{k,l} = \sum_{i,j,k,l} E_{i,l}\delta_{j,k} = n \sum_{i,l} E_{i,l} = n P.
$$
