Solution (source code)

= Solution

The <multiplicative order> $r$ makes the states $|a^k\bmod N\rangle$, $0\leq k<r$, distinct and cyclic under $U_a$. Therefore
$$
\begin{aligned}
U_a|\psi_s\rangle
&=\frac1{\sqrt r}\sum_{k=0}^{r-1}
e^{-2\pi isk/r}|a^{k+1}\bmod N\rangle\\
&=e^{2\pi is/r}\frac1{\sqrt r}\sum_{j=0}^{r-1}
e^{-2\pi isj/r}|a^j\bmod N\rangle.
\end{aligned}
$$
Thus each $|\psi_s\rangle$ is an <eigenvector> with
$$
\boxed{U_a|\psi_s\rangle=e^{2\pi is/r}|\psi_s\rangle}.
$$
These are the Fourier eigenvectors of the cyclic modular-multiplication orbit.