I would like to find a check whether a given matrix is unitarily equivalent to a permutation matrix. I split this question into two parts below. I expect the solution will use the second part, but if not that's fine.
${\bf(1)}$ Let $M$ be a unitary matrix. How can I check if there exists a unitary $U$ such that $U M U^\dagger = P$ is a permutation matrix?
I know that a permutation matrix has eigenvalues that form complete sets of roots of unity. The question then boils down to.
${\bf(2)}$ Given a set of complex numbers on the unit circle $S = \{z_1,z_2,z_3,\ldots\}$, is there a simple test to check if this contains complete sets of roots of unity. In other words, check that $S= S_1 \cup S_2 \cup \cdots$, where $S_i = \{e^{2\pi i \frac{n}{q}} | n \in \{0,\ldots,q-1\} \}$, for $q\in \mathbb{N}$.
Examples
$\{1,1,1,-1\} = \{1\} \cup \{1\} \cup \{1,-1\}$ and so this does have complete sets of roots of unity.
$\{1,i,-i,1\}$ does not contain complete sets of roots of unity.
Main ideas
i.) Unitary matrices (including permutation matrices as a special case) are normal, so they are unitarily similar to a diagonal matrix $D$, which is itself unitary. Thus to test whether we can have $UMU^* = P$, with given $M$ and some $\text{n x n}$ Permutation matrix $P$ it suffices to check that there is some (unitary) diagonal matrix $\Lambda$ that they could both be similar to. Since diagonalizability is guaranteed, it suffices to check some algebraic invariant like the characteristic polynomial.
ii.) With $M = QDQ^*$ and $P= V\Lambda V^*$ The algebraic invariant of interest is in fact the trace. I.e. examine $\text{trace}\big(M^k\big)= \text{trace}\big(D^k\big)$ for $k\in\{1,2,...,n\}$ and checking whether it could be equal to $\text{trace}\big(P^k\big) = \text{trace}\big(\Lambda^k\big)$ where $\Lambda$ is diagonal and has the constraint that it only has complete sets of $r_i$th roots of unity, with each $r_i \in \{1,2,...,n\}$. The trace is of interest since it counts closed walks on graphs and any permutation matrix $P$ is an adjacency matrix for a graph. Of course $P$, if it exists, may be written in a particularly simple form-- counting the number of eigenvalues of 1 tells us how many communicating classes it has (up to graph isomorphism this should be thought of as disjoint blocked cyclic permutation submatrices); this result follows e.g. from Perron-Frobenius Theory. We can even do better and construct a simple $P$ here, if it exists, by paying attention to the specific rth roots of unity and associated 'multiplicities'. The algorithm below actually gives a way to construct a 'nice' Permutation matrix $P$ that is similar to $M$ (and fails if some similarity is not possible).
Part (2) of the OP is a good question in that it maps question (1) down to a simpler problem, but implies the result for (1). I give an algorithmic solution to (2), below. Once that is understood, one can re-run the below algorithm verbatim, except compute $\text{trace}\big(M^k \big)$ and $\text{trace}\big(M^k - \Lambda^k\big)$ instead of $\text{trace}\big(D^k\big)$ and $\text{trace}\big(D^k - \Lambda^k\big)$.
For avoidance of doubt for (2), instead of a "set" I interpret this as a "multi-set" $S$ and $\big \vert S\big \vert = n$. With part 1 in mind this means we are interested solely in finding complete sets of $r_i$th roots of unity for $r_i \in \{1,2,...,n\}$.
The below algorithm assumes we are working with exact precision (or symbolically). Issues surrounding round-off for arbitrarily large matrices... are way outside the scope.
algorithm for (2) which implies (1)
Suppose that your multiset of eigenvalues / numbers on the unit circle correspond to a Permutation matrix. Place them in diagonal matrix $D$ and focus on taking the trace for powers of $D$. Also construct n x n diagonal matrix $\Lambda$ where $\lambda_i = 0$ for $i\in\{1,2,..,n\}$ at initialization (i.e. it is initialized to be the zero matrix).
Now compute $\text{trace}\big(D\big) = \alpha$. If $\alpha =0$ proceed to next step. If $\alpha$ is not a natural number you can reject the hypothesis that $D$ is similar to a permutation matrix, and the algorithm terminates. If $\alpha\gt 0$ (i.e. $\alpha$ is a positive integer), then set the first available $\alpha$ components on the diagonal of $\Lambda$ equal to one, i.e. $\lambda_i :=1$ for $i\in\{1,2,...,\alpha\}$, and now proceed to the next step.
Compute
$\text{trace}\big(D^k - \Lambda^k\big) = \text{trace}\big(D^k\big) - \text{trace}\big(\Lambda^k\big)$
for $k=2,3,...,n$.
It should be zero or a positive integer throughout -- if you see anything else, the algorithm terminates and we reject the eigenvalues as belonging to any $\text{n x n}$ permutation matrix. (This is essentially an assert line in a for loop.)
As soon as you see some positive integer, you want to 'nullify it'. (The idea to 'nullify' comes from the idea of reducibility of a graph in Perron Frobenius theory.) Put differently, we have $\text{trace}\big(D^1 - \Lambda^1\big)=0$ and we want to replicate this behavior for higher powers.
Suppose e.g. $\text{trace}\big(D^2-\Lambda^2\big)=5$. This implies $5$ 2nd roots of unity that did not 'show up' in $\text{trace}\big(D^1\big)$ -- this is a problem... because if they are complete sets we need $\text{5%2} =0$. So if $\text{trace}\big(D^2- \Lambda^2\big)=\alpha$ and $\alpha \text{% 2} \neq 0$ the algorithm terminates and we reject the eigenvalues as being associated with any $\text{n x n}$ permutation matrix $P$. Otherwise if $\alpha \text{% 2} = 0$ and $\alpha$ is a positive integer, proceed to 'nullify' these roots of unity. Say e.g. $\alpha = 6$. Then place 3 complete sets of 2nd roots of unity, in 3 contiguous 'diagonal blocks' in the the next available slots of $\Lambda$ (i.e. first remaining zeros on the diagonal) .
Now proceed to compute $\text{trace}\big(D^3-\Lambda^3\big)$ and repeat the above process, then $\text{trace}\big(D^4-\Lambda^4\big)$ and repeat above process, then $\text{trace}\big(D^5-\Lambda^5\big)$, and so on.
After this has been done all the way up to and including $k=n$ without raising an error, then we know
$\text{trace}\big(D^k-\Lambda^k\big) = 0 \longrightarrow \text{trace}\big(D^k\big) = \text{trace}\big(\Lambda^k\big)$
for $k \in\{1,2,...,n\}$
By Newton's Identities, $D$ has the same characteristic polynomial (and hence eigenvalues) as $\Lambda$. I.e. up to graph isomorphism $D = \Lambda$.
The underlying idea here is that sums of sets of complete $r_i$th roots of unity ($r_i \geq 2)$ telescope for powers less than $r_i$ so each time we append new $r_i$th roots of unity to $\Lambda$ we don't impact the $\text{trace}\big(\Lambda^k\big)$ for $k\in \{1,2,...,r_i-1\}$, and this allows us to exploit a kind of optimal substructure.
closing ideas:
Newton's Identities implies $\big \vert\det\big(\Lambda \big )\big \vert = \big \vert\det\big(D\big )\big \vert = 1$ -- so if the algorithm made it through $k=n$ without raising an error, then $\Lambda$ has no zeros on its diagonal. Furthermore, we constructed all the non-zero diagonal components of $\Lambda$ to come in diagonal 'blocks' $r_i$ at a time corresponding to complete sets of $r_i$th roots of unity. Thus $\Lambda$ can be separated into $m$ distinct complete sets of $r_i$ roots of unity for $r_i \in \{1,2...,n\}$ for $i\in\{1,2,....,m\}$ where $m$ is the number of 1's in $\Lambda$ -- and e.g. Perron Frobenius Theory tells us the associated permutation matrix $P$ has $m$ distinct communicating classes, each of size $r_i$.
This motivates a different approach to the algorithm -- rather than constructing $\Lambda$ as we go, we could construct a particularly nice $P$ as the algorithm proceeds. e.g. if at step one we had $\alpha = 5$, then the next step had $\alpha = 4$, and then $\alpha = 3$, then after 3 iterations we'd have a candidate $P= \begin{bmatrix} I_{5} & \mathbf 0 & \mathbf 0& \mathbf 0 & \mathbf 0 & \cdots & \mathbf 0 \\ \mathbf 0& C^{(2)} & \mathbf 0& \mathbf 0 & \mathbf 0 & \cdots & \mathbf 0 \\ \mathbf 0& \mathbf 0& C^{(2)}&\mathbf 0 & \mathbf 0 & \cdots & \mathbf 0 \\ \mathbf 0& \mathbf 0& \mathbf 0& C^{(3)} & \mathbf 0 &\cdots & \mathbf 0 \\ \mathbf 0& \mathbf 0& \mathbf 0& \mathbf 0& \mathbf 0 & \cdots & \mathbf 0 \\ \vdots& \vdots& \vdots& \vdots& \vdots & \ddots & \vdots \\ \mathbf 0& \mathbf 0 & \mathbf 0& \mathbf 0& \mathbf 0 & \cdots& \mathbf 0 \end{bmatrix}$
where $C^{(j)}$ is the $\text{j x j}$ 'cyclic permutation matrix' which is also the Companion matrix associated with the polynomial $x^j -1$. Of course, re-visiting this entire process, it becomes clear that we never needed the diagonal matrix $D$ that our unitary matrix $M$ is similar to -- we could have simply worked with traces of powers of $M$.