TL;DR
I'm trying to prove the following conjecture:
For every positive integer $N$, there exists a unitary matrix $M=[m_{ij}]_{N\times N}$ (i.e., $M^* M = MM^* = I$) such that every element of $M$ has absolute value equal to $1/\sqrt{N}$ (i.e., $|m_{ij}|=1/\sqrt{N}$).
I don't know if that statement is true or not, but I have reasons to believe it is. I'm particularly interested in odd values of $N$.
Any kind of insight or reference (Wikipedia article, etc.) is much appreciated.
Examples
I was able to find examples of such matrices for odd $N$ up to $9$:
- $N=1$ (trivial):
$$M_1 = \begin{bmatrix}1\end{bmatrix}$$
- $N=3$:
$$a = e^{2\pi i/3}$$
$$M_3 = \frac{1}{\sqrt{3}}\begin{bmatrix}1&a&\bar{a}\\ 1&\bar{a}&a\\ 1&1&1\end{bmatrix}$$
- $N=5$:
$$a = e^{2\pi i\cdot1/5}$$ $$b = e^{2\pi i\cdot2/5}$$
$$M_5 = \frac{1}{\sqrt{5}}\begin{bmatrix}1&a&\bar{a}&b&\bar{b}\\ 1&\bar{a}&a&\bar{b}&b\\ 1&b&\bar{b}&\bar{a}&a\\ 1&\bar{b}&b&a&\bar{a}\\ 1&1&1&1&1\end{bmatrix}$$
- $N=7$:
$$a = e^{2\pi i\cdot1/7}$$ $$b = e^{2\pi i\cdot2/7}$$ $$c = e^{2\pi i\cdot3/7}$$
$$M_7 = \frac{1}{\sqrt{7}}\begin{bmatrix}1&a &\bar{a}&c &\bar{c}&b &\bar{b}\\ 1&\bar{a}&a &\bar{c}&c &\bar{b}&b \\ 1&b &\bar{b}&\bar{a}&a &\bar{c}&c \\ 1&\bar{b}&b &a &\bar{a}&c &\bar{c}\\ 1&c &\bar{c}&b &\bar{b}&\bar{a}&a \\ 1&\bar{c}&c &\bar{b}&b &a &\bar{a}\\ 1&1 &1 &1 &1 &1 &1\end{bmatrix}$$
- $N=9$:
$$a = e^{2\pi i\cdot1/9}$$ $$b = e^{2\pi i\cdot2/9}$$ $$c = e^{2\pi i\cdot3/9}$$ $$d = e^{2\pi i\cdot4/9}$$
$$M_9 = \frac{1}{\sqrt{9}}\begin{bmatrix}1&a &\bar{a}&d &\bar{d}&\bar{b}&b &c &\bar{c}\\ 1&\bar{a}&a &\bar{d}&d &b &\bar{b}&\bar{c}&c \\ 1&b &\bar{b}&\bar{a}&a &\bar{d}&d &\bar{c}&c \\ 1&\bar{b}&b &a &\bar{a}&d &\bar{d}&c &\bar{c}\\ 1&c &\bar{c}&c &\bar{c}&c &\bar{c}&1 &1 \\ 1&\bar{c}&c &\bar{c}&c &\bar{c}&c &1 &1 \\ 1&d &\bar{d}&\bar{b}&b &a &\bar{a}&c &\bar{c}\\ 1&\bar{d}&d &b &\bar{b}&\bar{a}&a &\bar{c}&c \\ 1&1 &1 &1 &1 &1 &1 &1 &1\end{bmatrix}$$
It took me a long time to find those matrices and I can't see any pattern that would make it easy to generalize the construction for arbitrary odd values of $N$.
Context
I was studying the solution to the Schrödinger equation for the hydrogen atom and I noticed that the usual expressions for the probability distribution for the orbitals higher than s didn't have spherical symmetry:
This is perfectly reasonable if you measure the angular momentum of an electron and choose a reference frame such that one of the axis is aligned with the angular momentum vector. But before measuring anything, the system has spherical symmetry, so the probability distribution also has to be spherically symmetric.
For a single electron, it's easy to see how this works: the electron is in a superposition of all the orbitals of the same type. For example, for the orbital $2p$:
$$\Psi_{2p} = \frac{1}{\sqrt{3}}|2p_x\rangle + \frac{1}{\sqrt{3}}|2p_y\rangle + \frac{1}{\sqrt{3}}|2p_z\rangle$$
where the wave functions for each state $|2p\rangle$ are given by:
$$\Psi_{2p_x}(r,\theta,\phi) = R(r)\sqrt{\frac{3}{8\pi}}\sin\theta \;e^{-i\phi}$$
$$\Psi_{2p_y}(r,\theta,\phi) = R(r)\sqrt{\frac{3}{8\pi}}\sin\theta \;e^{i\phi}$$
$$\Psi_{2p_z}(r,\theta,\phi) = R(r)\sqrt{\frac{3}{4\pi}}\cos\theta$$
where
$$R(r) = \frac{1}{\sqrt{24a_0^3}}\cdot\frac{r}{a_0}\cdot e^{-r/2a_0}$$
The probability distribution is given by:
$$\begin{array}{} P_{2p}(r,\theta,\phi) &= \langle\Psi_{2p}|\Psi_{2p}\rangle \\ &= \frac{1}{3}\langle\Psi_{2p_x}|\Psi_{2p_x}\rangle + \frac{1}{3}\langle\Psi_{2p_y}|\Psi_{2p_y}\rangle + \frac{1}{3}\langle\Psi_{2p_z}|\Psi_{2p_z}\rangle \\ & = R^2(r)\end{array}$$
So the electron is indeed spherically distributed around the nucleus.
For more electrons, things get much more complicated. Because of the Pauli exclusion principle, the wave functions need to be orthogonal to each other (ignoring spin):
$$\langle\Psi_{i}|\Psi_{j}\rangle = \delta_{ij}$$
For the $2p$ orbital with 3 electrons (ignoring spin), one solution is:
$$\Psi_{2p_1} = \frac{1}{\sqrt{3}}|2p_x\rangle + \frac{e^{2\pi i/3}}{\sqrt{3}}|2p_y\rangle + \frac{e^{-2\pi i/3}}{\sqrt{3}}|2p_z\rangle$$
$$\Psi_{2p_2} = \frac{1}{\sqrt{3}}|2p_x\rangle + \frac{e^{-2\pi i/3}}{\sqrt{3}}|2p_y\rangle + \frac{e^{2\pi i/3}}{\sqrt{3}}|2p_z\rangle$$
$$\Psi_{2p_3} = \frac{1}{\sqrt{3}}|2p_x\rangle + \frac{1}{\sqrt{3}}|2p_y\rangle + \frac{1}{\sqrt{3}}|2p_z\rangle$$
You can see that the problem of finding the coefficients for each eigenfunction is the same as finding the matrices $M$ described above. In the case of the $2p$ orbital, for example, the coefficients are the elements of $M_3$:
$$\begin{bmatrix} \Psi_{2p_1}\\\Psi_{2p_2}\\\Psi_{2p_3}\end{bmatrix} = M_3 \begin{bmatrix} |2p_x\rangle\\|2p_y\rangle\\|2p_z\rangle \end{bmatrix}$$
Related Questions
Counting symmetric unitary matrices with elements of equal magnitude
Well found on your own! It is the Discrete Fourier Transform The orthonormal colomn vectors are: $$u_k = \frac{1}{\sqrt{N}} \left(\exp(\frac{2\pi i }{N} kn)\right)_{0\leq n<N}$$ for $k=0,...,N-1$