For a positive integer $n$, let $S_n$ denote the set of $n\times n$ symmetric matrices over $\mathbb{C}$. As a complex vector space, this set has dimension $\mathrm{dim}(S_n)=\binom{n+1}{2}$.  The standard Hilbert-Schmidt inner product
$$
\langle A,B\rangle = \mathrm{Tr}(AB^*) 
$$
for $A,B\in S_n$ turns this space into an inner product space.

**Question.** Does there always exist an orthogonal basis $\{U_1,U_2,\dots,U_{\frac{n(n+1)}{2}}\} \subseteq S_n$ of symmetric matrices such that each $U_i$ is unitary? 

-------------
**Case for $n=2$**

The case $n=2$ is simple, since we may take the matrices
$$
U_1 = \begin{pmatrix} 1&0\\0&1\end{pmatrix},\quad
U_2 = \begin{pmatrix} 1&0\\0&-1\end{pmatrix},\quad\text{and}\quad
U_3 = \begin{pmatrix} 0&1\\1&0\end{pmatrix}.
$$
These matrices are all symmetric and unitary, and they satisfy $\mathrm{Tr}(U_iU_j^*) = 0$ whenenver $i\neq j$. Moreover, these matrices span the space of $2\times 2$ symmetric matrices.

---------
**Case for $n$ even**

I've come up with a way to construct such a basis for any even $n$. 

Consider the matrices $\{E_{i,j}\,:\, i,j\in\{1,\dots,n\}\}$, where $E_{i,j}$ is the $n\times n$ matrix that has a 1 in the $i$th row and $j$th column with zeros elsewhere. One orthogonal basis for the set of $n\times n$ symmetric matrices is the set 
$$
\{H_{i,j}: i,j\in\{1,\dots,n\},\, i\leq j\}
$$
where we define
$$
H_{i,j} = \left\{\begin{array}{ll}E_{i,i} & \text{if }i=j\\
\frac{1}{\sqrt{2}}(E_{i,j}+E_{j,i}) & \text{if }i\neq j\end{array}\right.
$$
Consider the complete graph $K_n$ with $n$ vertices. We may identify the edges of this graph with the set
$$
\mathcal{E} = \{H_{i,j} : i,j\in\{1,\dots,n\},\, i<j\}
$$
 where the edge connecting vertices $i$ and $j$ (with $i<j$) is denoted $H_{i,j}$. Note that $\mathcal{E}$ has $\binom{n}{2}=\frac{n(n-1)}{2}$ elements.

If $n$ is even, there is a 1-factorization of this graph ([see here][1]). A 1-factorization corresponds to a partitioning the edges $\mathcal{E}$ into $n-1$ subsets, $\mathcal{F}_1,\dots,\mathcal{F_{n-1}}$, each with $n/2$ edges, such that, for each $k$, no two edges in $\mathcal{F_k}$ are adjacent and
$$
\mathcal{F_1}\cup\cdots\cup\mathcal{F_{n-1}} = \mathcal{E}.
$$
For each $k\in\{1,\dots,n-1\}$, we can label the elements of $\mathcal{F_k}$ as
$$
\mathcal{F_k} = \{F_{k,1},\dots,F_{k,n/2}\}.
$$
Let $\omega = e^{i2\pi/n}$ denote the $n$th root of unity, and for each $k\in\{1,\dots,n-1\}$ and $\ell\in\{1,\dots,\frac{n}{2}\}$ we define the matrix
$$
G_{k,\ell} = 2\sum_{a=1}^{n/2} \omega^{2\ell a} F_{k,a}.
$$
It can be verified that $G_{k,\ell}$ is symmetric and unitary. Finally, define the symmetric unitary matrices $A_1,\dots,A_n$ defined by
$$
A_j = \sum_{a=1}^n \omega^{ja} E_{a,a}
$$
for each $j\in\{1,\dots n\}$. 
It can be shown that the set
$$
\{A_1,\dots,A_n\}\cup\{G_{k,\ell}\, :\, k\in\{1,\dots,n-1\},\, \ell\in\{1,\dots,n/2\}\}
$$
is an orthogonal basis of $S_n$. 

---------

**Case when $n=3$**

The case for $n$ odd seems to be a bit trickier. I've at least been able to construct a basis of the desired form when $n=3$ as follows:
\begin{align*}
U_1 &= \frac{\sqrt{3}}{2}\begin{pmatrix} \frac{2}{\sqrt{3}} & 0 & 0\\ 0 &\frac{1}{\sqrt{3}} & 1\\ 0 & 1 & -\frac{1}{\sqrt{3}}\end{pmatrix} 
& U_2&=\frac{\sqrt{3}}{2}\begin{pmatrix} \frac{2}{\sqrt{3}} & 0 & 0\\ 0 &\frac{1}{\sqrt{3}} & -1\\ 0 & -1 & -\frac{1}{\sqrt{3}}\end{pmatrix}\\
U_3 &= \frac{\sqrt{3}}{2}\begin{pmatrix} \frac{-\overline{\beta}}{\sqrt{3}} & 0 & 1\\ 0 &\frac{\alpha}{\sqrt{3}} & 0\\ 1 & 0 & \frac{\beta}{\sqrt{3}}\end{pmatrix} 
& U_4&=\frac{\sqrt{3}}{2}\begin{pmatrix} \frac{-\overline{\beta}}{\sqrt{3}} & 0 & -1\\ 0 &\frac{\alpha}{\sqrt{3}} & 0\\ -1 & 0 & \frac{\beta}{\sqrt{3}}\end{pmatrix} \\
U_5 &= \frac{\sqrt{3}}{2}\begin{pmatrix} \frac{i\overline{\beta}}{\sqrt{3}} & 1 & 0\\ 1 &\frac{i\beta}{\sqrt{3}} & 0\\ 0 & 0 & \frac{i\alpha}{\sqrt{3}}\end{pmatrix}
& U_6&=\frac{\sqrt{3}}{2}\begin{pmatrix} \frac{i\overline{\beta}}{\sqrt{3}} & -1 & 0\\ -1 &\frac{i\beta}{\sqrt{3}} & 0\\ 0 & 0 & \frac{i\alpha}{\sqrt{3}}\end{pmatrix}
\end{align*}
where $\alpha = \sqrt{\frac{27}{8}} - i\sqrt{\frac{5}{8}}$ and $\beta=\sqrt{\frac{3}{8}} + i\sqrt{\frac{5}{8}}$. It can be verified that the matrices $\{U_1,\dots,U_6\}$ are symmetric, unitary, and pairwise orthogonal.

---------

**Existence for odd $n$ greater than $3$?**

Numerically, I've been able to find an orthogonal basis of symmetric unitary matrices for odd dimensions up to $n=11$. But I found no discernible pattern that allowed me to construct an exact solution.

I'd like to know if there exists a collection of $\binom{n}{2}$ orthogonal symmetric $n\times n$ unitary matrices for any odd $n>3$.


  [1]: https://en.wikipedia.org/wiki/Graph_factorization#Complete_graphs