Veo que ha etiquetado su pregunta con computación cuántica. Así que quizás te resulte más intuitivo si trabajamos en notación de Dirac. Es decir, dejamos que $\{| j \rangle\}$ sea una base ortonormal para el espacio vectorial subyacente (de modo que $\langle i|j\rangle =\delta_{i,j}$ ). Entonces lo que usted llama las matrices unitarias son equivalentes definidas como $$ E_{i,j} = |i\rangle\!\langle j|. $$ Es decir, tienen un 1 en la entrada en el $i$ 'th fila y $j$ y ceros en el resto. Ahora el conjunto de matrices $\{|E_{i,j}\rangle\!\rangle \}$ forman una base ortonormal para el espacio vectorial de matrices, donde el producto interior entre dos matrices $A$ y $B$ es $$\langle\!\langle A|B\rangle\!\rangle = \operatorname{Tr}(A^\dagger B). $$
Ahora podemos hacer el mismo truco. Definamos las "super" matrices unitarias como $$ F_{i,j;k,l} = |E_{i,j}\rangle\!\rangle\!\langle\!\langle E_{k,l}|. $$ Entonces, el $F$ cumplen las mismas propiedades que las matrices $E$ matrices.
Sin embargo, por tu última frase, no parece que te interesen estas matrices unitarias. Tu ecuación final sugiere la matriz que defines como $P$ es simplemente la matriz de 1's en cada entrada. Se puede escribir $$ P = \sum_{i,j} E_{i,j}. $$ A continuación, el cálculo directo conduce a $$ P^\dagger = \sum_{i,j} E_{i,j}^\dagger = \sum_{i,j} E_{j,i} = P, $$ y $$ 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. $$
