This answer is incomplete, but it should provide an answer for almost all symmetric states (*i.e.* it suffices for all but a set of symmetric states having measure zero).

The symmetric subspace is spanned by product states. We may then consider different ways in which a particular symmetric state decomposes into symmetric products; in particular, if any choice of decomposition naturally gives rise to an invariant.

A greedy way to go about decomposing a symmetric $\def\ket#1{|#1\rangle}\ket\psi$ state into symmetric products would be to simply look for the symmetric product $\ket\phi\ket\phi\cdots\ket\phi$ with which $\ket\psi$ has the overlap of the largest magnitude. Let $\ket{\phi_0}$ be the single-spin state satisfying this, and $$\def\bra#1{\langle#1|}\alpha_0 = \Bigl[\bra{\phi_0}\cdots\bra{\phi_0}\Bigr]\ket\psi$$
which without loss of generality is positive. With probability 1, the state $\ket{\phi_0}$ is unique (in that the set of symmetric states for which it is not unique has measure zero). Let $\ket{\psi_1}$ be the projection of $\ket{\psi}$ into the orthocomplement of $\ket{\phi_0}^{\otimes n}$: this is another symmetric state. So, we define $\ket{\phi_1}$ to be the (again with probability 1, unique) single spin state such that $\ket{\phi_1}^{\otimes n}$ has maximal overlap with $\ket{\psi_1}$; we let $\alpha_1$ be that overlap; and we define $\ket{\psi_2}$ to be the projection of $\ket{\psi}$ onto the orthocomplement of $\mathrm{span}\{\ket{\phi_0}^{\otimes n}\!\!,\;\; \ket{\phi_1}^{\otimes n}\}$. And so on.

By continually projecting $\ket{\psi}$ onto the orthocomplement of spans of larger and larger sets of symmetric products, we ensure that the resulting projections $\ket{\psi_j}$ will not have maximal overlap with any product state that came before, or more generally which can be spanned by the preceding symmetric products. So at each iteration we obtain a single-spin state $\ket{\phi_j}$ such that $\mathrm{span}\{\ket{\phi_0}^{\otimes n}\!\!,\;\ldots\,,\; \ket{\phi_j}^{\otimes n}\}$ has dimension one larger than in the previous iteration. In the end, we will obtain a collection of symmetric products which, if they don't span the symmetric subspace, at least contain $\ket{\psi}$ in their span. So we obtain a decomposition
$$ \ket{\psi} \;=\; \sum_{j = 0}^\ell \alpha_j \ket{\phi_j}^{\otimes n} $$
where the sequence $\alpha_0, \alpha_1, \ldots$ is strictly decreasing. Let us call this the symmetric product decomposition of $\ket{\psi}$. (It wouldn't take too much to generalize this representation to one where the states $\ket{\phi_0}$, $\ket{\phi_1}$, etc. are non-unique; but for what comes next uniqueness will be important.) Given the symmetric product decomposition of $\ket{\psi}$, it is trivial to describe the corresponding representation for $U^{\otimes n} \ket{\psi}$: just multiply each of the $\ket{\phi_j}$ by $U$. And in fact, if you computed $U\ket{\psi}$ and then determined its symmetric product decomposition, the decomposition
$$ U^{\otimes n}\ket{\psi} = \sum_{j=0}^\ell \alpha_j \Bigl[ U\ket{\phi_j} \Bigr]^{\otimes n} $$
is exactly what you would find: $\ket{\phi_0}^{\otimes n}$ has maximal overlap with $\ket{\psi}$ if and only if $[ U \ket{\phi_0} ]^{\otimes n}$ has maximal overlap with $U^{\otimes n} \ket\psi$, and so on. So to show that two symmetric states are LU-equivalent, it suffices to show that the sequence of amplitudes $\alpha_j$ are the same, and that the sequence of single-spin states $\ket{\phi_j}$ are related by a common unitary.

The last part can be most easily done by finding a normal form for the states which are equal, for sequences of single-spin states which are related by a common single-spin unitary. We can do this by finding a unitary $T$ which

- maps $\ket{\phi_0}$ to $\ket{0}$,
- maps $\ket{\phi_1}$ to a state $\ket{\beta_1}$ in the span of $\ket{0}$ and $\ket{1}$, with $\bra{1}\beta_1\rangle \geqslant 0$,
- and for each subsequent $j > 1$, maps $\ket{\phi_j}$ to some state $\ket{\beta_j}$ which is in the span of standard basis states $\ket{0}, \ldots, \ket{b_j}$ for $b_j$ as small as possible, with $\bra{b_j}\beta_j\rangle \geqslant 0$ if possible. (For any state $\ket{\phi_j}$ which is in the span of preceding states $\ket{\phi_h}$, the state $\ket{\beta_j}$ will similarly be determined by the states $\ket{\beta_h}$ for $h < j$, in which case we do not have a choice of the value of $\bra{b_j}\beta_j\rangle$.)

We then have a unitary such that $T \ket{\phi_j} = \ket{\beta_j}$; and for any two sequences of states $\ket{\phi_j}$ and $U\ket{\phi_j}$, we should obtain the same sequence of states $\ket{\beta_j}$. You can then determine that two symmetric states are equivalent if they give rise to the same sequence of amplitudes $\alpha_j$ and the same "normal form" single-spin states $\ket{\beta_j}$.

In the case that there is not a unique state $\ket{\phi_j}^{\otimes n}$ which has maximal overlap with $\ket{\psi_j}$ in the construction of the symmetric product decomposition, the problem is then in defining the normal form states $\ket{\beta_j}$. However, so long as the states $\ket{\phi_j}$ are unique, which happens with measure 1, you should have a polynomial-size invariant (up to precision limitations) for determining if two symmetric states are the same.

This post has been migrated from (A51.SE)