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). 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 symmetric unitary matrices $A_1,\dots,A_n$ 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+1}{2}$ orthogonal symmetric $n\times n$ unitary matrices for any odd $n>3$.
