linear algebra - 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$$

Answer

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.$$