I conjecture that in the complex vector space $\mathbb{C}^{2^N \times 2^N}$, where $N$ is a positive integer, there is a basis $\mathcal{B}$ whose elements are Hermitian and unitary. However, I don't know how to prove (or counterprove) it.
Given a number $N$, it is always possible to find $2N$ distinct Hermitian and unitary matrices other than the identity $\{\gamma_1. \cdots. \gamma_{2N}\}$ such that each pair of them anti-commutes: $$\gamma_i \gamma_j + \gamma_j \gamma_i = 0 \ \text{if} \ i \neq j.$$ Physicists usually call those $\gamma$ Majorana operators, and the reason why it is always possible to find $2N$ such matrices is based on the relation between fermionic operators and Majorana operators. For example, if $N = 1$, we define a so-called fermionic annihilation operator as $$d = \begin{bmatrix} 0 & 1\\ 0 & 0 \end{bmatrix}. $$ Then the matrix $$\gamma_1 = d + d^\dagger = \begin{bmatrix} 0 & 1\\ 1 & 0 \end{bmatrix} $$ and the matrix $$\gamma_2 = i(d^\dagger - d) = \begin{bmatrix} 0 & -i\\ i & 0 \end{bmatrix} $$ are Majorana operators.
The set of $2N$ Majorana operators can be extended to a set of $4^N$ Hermitian and unitary matrices: $$\mathcal{B}=\left\{i^{\frac{n(n-1)}{2}} (\gamma_{a_1} \cdots \gamma_{a_n}): \ n, a_i \in \{1,\cdots,2N\}, \, a_i < a_{i+1}\right\} \bigcup \{1\}.$$ For instance, when $N = 1$, $$\mathcal{B}=\left\{1, \gamma_1,\gamma_2, i \gamma_1 \gamma_2\right\}.$$
A direct calculation shows $$|\mathcal{B}| = 1 + \sum_{n=1}^{2N} \begin{pmatrix} 2N\\n \end{pmatrix} = \sum_{n=0}^{2N} \begin{pmatrix} 2N\\n \end{pmatrix} = 2^{2N} = 4^N = \dim\left(\mathbb{C}^{2^N \times 2^N}\right).$$ I feel this set $\mathcal{B}$ is a basis for the complex vector space. Indeed, it is true at least when $N = 1$: $\mathcal{B}$ is the set of Pauli matrices. Nevertheless, I don't know how to prove it in higher dimensional spaces.
Thanks in advance for your help. :)
I worked it out myself. This is the link to my proof.