The set of size-$n$ unitary matrices span $\Bbb C^{n \times n}$ (this can be proven nicely using polar decomposition). If we select a maximal linear subset of unitary matrices, then we have a basis of $\Bbb C^{n \times n}$ consisting of $n^2$ unitary matrices. My question is whether we can find a basis that satisfies the additional constraint of *orthogonality*. That is:

Does there exist a basis $\mathcal B$ of $\Bbb C^{n \times n}$ such that every $P \in \mathcal B$ is unitary (that is, $P^*P = I$) and for all distinct $P,Q \in \mathcal B$, we have $\langle P,Q \rangle = 0$?

Here, $\langle \cdot , \cdot \rangle$ refers to the Frobenius (AKA Hilbert-Schmidt) inner product, namely $\langle P, Q \rangle = \operatorname{trace}(PQ^*)$.

When $n = 2$, the Pauli matrices provide a convenient solution. That is, we can take
$$
\mathcal B = \{I,\sigma_1,\sigma_2,\sigma_3\} \subset \Bbb C^{2 \times 2}.
$$
We can use this to produce a solution whenever $n = 2^k$. In particular, if we define $\sigma_0 = I$ for convenience, we can take
$$
\mathcal B = \{\sigma_{m_1} \otimes \cdots \otimes \sigma_{m_k} : 0 \leq m_j \leq 3\} \subset \Bbb C^{2^k \times 2^k}.
$$
Could we come up with a basis for any other $n$? Could we do so for *every* $n$?

Some observations so far:

- Without loss of generality, we can assume that $\mathcal B$ contains the $n \times n$ identity matrix $I$. If $I$ is an element of the basis, it follows that the remaining matrices form a basis for the subspace of all trace-$0$ matrices.
- A commuting set of matrices spans at most an $n$-dimensional subset, so there must be elements of $\mathcal B$ that fail to commute