I think there's a typo here. I think it's supposed to say that without loss of generality $S$ is a rotation by angle $\theta$ (and hence your proposed parametrisation of $S$ is correct). The most direct evidence I have for this is the eigenvalues of $U$: if $U$ were a rotation by an angle $\theta$, its eigenvalues would be of the form $e^{\pm i\theta}$.
Given this change, how does the calculation proceed? First, verify that the two eigenvectors you're given really are eigenvectors with eigenvalue 1. There are some tricks that help with this. Personally, I rewrote the first eigenvector (for example) as
$$
(I\otimes H)(|01\rangle+|10\rangle)/\sqrt{2}.
$$
Then I applied $U$ to this, and noted that:
- if I commute Hadamard through controlled-not, it changes into controlled-phase
- if I commute Hadamard through $S(\theta)$, I have $S(\theta)H=HS(-\theta)$ because $HYH=-Y$.
- $S(\theta)\otimes S(-\theta)(|01\rangle+|10\rangle)=(|01\rangle+|10\rangle)$.
So, let's assume we know that two of the eigenvalues of $U$ are 1. What are the others? We can easily calculate that $\text{det}(U)=1$ since $\text{det}(S)=1$. This means that our other two eigenvalues are of the form $e^{\pm i\alpha}$ since this is the only way they can multiply together to give 1. Then we calculate the trace, which is equal to the sum of the eigenvalues:
$$
2+2\cos\alpha=4\cos^4\theta
$$
(I just brute-forced this in Mathematica.)
This is just the same as
$$
4\cos^2\frac{\alpha}{2}=4\cos^4\theta.
$$
Since we're looking for $\pm\alpha$ solutions. Taking the square root leaves us where you need to be.
The second eigenvector: If you look at the second stated eigenvector, it only appears to be a 1-qubit state. So there must be another typo. I believe it should be
$$
\sin(\theta)|0\rangle|-\rangle-|1\rangle(\sin(\theta)|+\rangle+2\cos(\theta)|-\rangle)
$$
Here are the calculations I did in Mathematica. First, define the matrix and verify the first eigenvector (should give all-zeros vector):
rootU = KroneckerProduct[{{Cos[q], -Sin[q]}, {Sin[q],
Cos[q]}}, {{Cos[q], -Sin[q]}, {Sin[q], Cos[q]}}].{{1, 0, 0,
0}, {0, 1, 0, 0}, {0, 0, 0, 1}, {0, 0, 1, 0}};
U = rootU.rootU;
(U - IdentityMatrix[4]).{1, -1, 1, 1} // FullSimplify
The problem is if you just ask it to solve for the two eigenvectors of +1 eigenvalue, it doesn't know how to pick them within the degenerate subspace. So, we help. Remove the one eigenvector we already know about, so that there is only one eigenvalue of 1:
V = U - {{1, -1, 1, 1}, -{1, -1, 1, 1}, {1, -1, 1, 1}, {1, -1, 1, 1}}/
4;
-(Sin[q] - 2 Cos[q]) Eigenvectors[V, 1] // FullSimplify
The form of the other two eigenvalues does not seem to drop out nicely, which is why I went the direction I did to reproduce the results of the paper.