This is a partial answer but too long for a comment.

Let us first prove the property stated i.e.
$$
\mathbf{A} \mbox{ unitary }\Longrightarrow \mbox{ all roots of }p(z) \mbox{ are of modulus } 1
$$
Firsly $|p(0)|=|det(-\mathbf{A})|=1$,
so $p$ has not zero as a root.

Second, if $z\not=0$, one can write
$$
det(\mathbf{D}(z)-\mathbf{A})
=det(\mathbf{D}(z))det(I-\mathbf{D}(z)^{-1}\mathbf{A})
$$
so $p(z)=0$ is equivalent to $1\in sp(\mathbf{D}(z)^{-1}\mathbf{A})$ and to the existence of $\mathbf{v}\not=0$ such that
$$
\mathbf{D}(z)^{-1}\mathbf{A}\mathbf{v}=\mathbf{v}
$$
for such $\mathbf{v}$ one has
$
\mathbf{A}\mathbf{v}=\mathbf{D}(z)\mathbf{v}
$
and, writing $z=\rho e^{it}$ one gets
$
\mathbf{A}\mathbf{v}=\mathbf{D}(\rho)\mathbf{D}(e^{it})\mathbf{v}\ .
$
Now, $\mathbf{D}(e^{it})$ being unitary, one gets finally
$$
||\mathbf{v}||_2=||\mathbf{D}(e^{-it})\mathbf{A}\mathbf{v}||_2=||\mathbf{D}(\rho)\mathbf{v}||_2\qquad \mbox{(*).}
$$
As it is easy to check that for all $\mathbf{v}$ and $\rho>0$,
$$
||\mathbf{D}(\rho)\mathbf{v}||_2\geq \rho ||\mathbf{v}||_2 \mbox{ if } \rho>1\ ;\ ||\mathbf{D}(\rho)\mathbf{v}||_2\leq \rho ||\mathbf{v}||_2 \mbox{ if } \rho<1\ ,
$$
the result follows from (*).

For $\alpha=(m_1,\cdots ,m_N)\in (\mathbb{N}_{\geq 1})^N$, let
$$
\mathbf{D}(\alpha,z)=diag(z^{m_1}, \dots, z^{m_N})
$$
for $z$ fixed, one has
$$
\mathbf{D}(\alpha,z)\mathbf{D}(\beta,z)=\mathbf{D}(\alpha+\beta,z)
$$
these matrices form a semigroup.

For $\alpha$ fixed

$$
\mathbf{D}(\alpha,z_1)\mathbf{D}(\alpha,z_2)=\mathbf{D}(\alpha,z_1z_2)
$$
for $z\not=0$ these matrices form a group. The set of all these matrices is normalised by the monomial matrices, i.e. the semi-direct product of the (Weyl) group of permutation matrices and of diagonal matrices as, if $W=W(\sigma)$ is a permutation matrix and if $D$ is a diagonal (regular) matrix, one has
$$
WD\mathbf{D}(\alpha,z)D^{-1}W^{-1}=\mathbf{D}(\alpha_\sigma,z)
$$

As was remarked by Victor, matrices such as unitary (see above) and upper (or lower) triangular with unitary diagonal possess the property.

Their conjugates through the **monomial group** possess also the property.

**How to test (algorithmically) that a matrix is conjugated of a unitary matrix through the monomial group ?**

Firstly, as was remarked as the permutation matrices and as the diagonal ones with unitary spectrum are unitary, to be such is equivalent of being conjugated of a unitary matrix through the diagonal group of matrices with strictly positive eigenvalues i.e. for a matrix $B$ test whether it exists a unitary matrix and $R=diag(r_1,\cdots ,r_N)$ such that
$$
B=R^{-1}AR
$$
(one can even restrict to the special group of them, but we will not use this)

**(analysis)** suppose it were the case, then $B^*R^2B=R^2$

**(synthesis)**

- find all the diagonal matrices $D$ which fulfil $B^*DB=D$ (it is a linear system with $N$ variables, i.e. diagonal eigenvectors for the eigenvalue $1$ of the linear transformation $D\rightarrow B^*DB$).
- among them select, if possible, a $D$ with strictly positive spectrum and set $R=\sqrt{D}$ then $A=R^{-1}AR$ is unitary.

**Remark** It can happen that the transformation $D\rightarrow B^*DB$ admit diagonal eigenvectors for the eigenvalue $1$ none of which is strictly positive. As the procedure provides a necessary and sufficient condition, the corresponding matrix is **not** diagonally conjugate to a unitary matrix.

**(Counter)-example** Set
$$
B=
\begin{pmatrix}
\sqrt{2} & 1\cr
1 & \sqrt{2}
\end{pmatrix}
$$
then, solving
$$
B^*\begin{pmatrix}
x & 0\cr
0 & y
\end{pmatrix}
B
=
\begin{pmatrix}
x & 0\cr
0 & y
\end{pmatrix}
$$
yields $x=-y$ so there are eigenvectors as
$$
\begin{pmatrix}
1 & 0\cr
0 & -1
\end{pmatrix}
$$
but none of them is strictly positive.