Eigenvalue pattern We consider a matrix
$$M_{\mu} = \begin{pmatrix} 1 & \mu & 1 & 0 \\ -\mu & 1 & 0 & 1 \\ -1 & 0 & 0 & 0 \\ 0 &-1 & 0 & 0 \end{pmatrix}$$
One easily checks that $\operatorname{det}(M_{\mu})=1$.
I however noticed something peculiar:
Consider a sequence of real numbers $\mu_i$ then the four eigenvalues $\lambda_1,..,\lambda_4$ of
$$ A=\prod_{i=1}^n M_{\mu_i}$$
have the property that they can be chosen to satisfy $\lambda_1 = \overline{\lambda_2}$ and $\lambda_3 = 1/\lambda_1$ and consequently by the determinant $\lambda_4 = 1/\lambda_2.$
I tried to explicitly study the product and see if I can see some structure explaining all this, but so far only with limited luck.
Question: Show that the eigenvalues of $A$ for every $n$ and $\mu_i \in \mathbb R$ satisfy $$\lambda_1 = \overline{\lambda_2 } = 1/\lambda_3 = 1/\overline{\lambda_4}.$$
Update: What we know so far is that $\lambda_1 = \frac{1}{\lambda_3}$ and $\lambda_2 = \frac{1}{\lambda_4}$ and if any of the eigenvalues is non-real then we have it. However, we cannot exclude $\lambda_1>\lambda_2$ both real.
 A: EDIT: just a partial answer that does not settle the question completely.
This is a variant of symplectic matrices. Your matrices $M_\mu$ are orthogonal wrt the indefinite scalar product induced by
$$
\Omega = \begin{bmatrix}
0 & 0 & 0 & 1\\
0 & 0 & 1 & 0\\
0 & -1 & 0 & 0\\
-1 & 0 & 0 & 0
\end{bmatrix},
$$
i.e., they satisfy the property
$$
M^{*}\Omega M = \Omega.
$$
One can see that if $(\lambda, v)$ is an eigenpair for $M$, then $(\frac{1}{\lambda}, \Omega v)$ is an eigenpair for $M^*$, and hence $\frac{1}{\overline{\lambda}}$ is in the spectrum of $M$:
$$ \frac{1}{\lambda}\Omega v = \frac{1}{\lambda}M^*\Omega Mv = \frac{1}{\lambda}M^*\Omega \lambda v = M^*\Omega v.$$
This property, together with the fact that real matrices have eigenvalues that come in complex conjugate pairs, puts strong constraints on  which eigenvalues you can get.
There is also the possibility that the eigenvalues come in the form $\lambda_1,\lambda_2, \frac1{\lambda_1}, \frac1{\lambda_2}$, for $\lambda_1,\lambda_2\in \mathbb{R}$. This case actually happens: for instance, $M_2M_{-2}$ has four real eigenvalues, although in this case I get $\lambda_1 = \lambda_2$.
Four purely imaginary eigenvalues can also appear: for instance $M_{-2}M_{1/2}$ has eigenvalues $\{\pm 2i, \pm \frac12 i\}$.
EDIT: on further thought, I still cannot exclude that the eigenvalues are in the form suggested by Sascha in a comment; this is not a counterexample as the real eigenvalues have multiplicity 2. Sorry :(
A: Your matrix $M_\mu$ is symplectic: $M_\mu^T\Omega M_\mu=\Omega$ where
$$\Omega=\begin{pmatrix} 0_2 & Y \\ -Y & 0_2 \end{pmatrix},\qquad Y=\begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}.$$
Then every product $A$ is still symplectic. This implies $A^T\Omega=\Omega A^{-1}$. because the spectrum of $A^T$, equal to that of $A$, is the inverse of that of $A^{-1}$, this shows that the spectrum is stable under $z\mapsto 1/z$. Because the matrix is real, the spectrum is invariant under $z\mapsto \bar z$ as well.
There remains to eliminate the possibility of a spectrum $\lambda_1,\lambda_2,\lambda_1^{-1},\lambda_2^{-1}$ with both $\lambda_j$ real... TODO
A: The explanation is pretty simple with a suitable change of basis.
Letting
$$B = 
\begin{pmatrix} 
1 & 0 & 1 & 0 \\
i & 0 & -i & 0 \\
0 & 1 & 0 & 1 \\ 
0 & i & 0 & -i 
\end{pmatrix}$$
we have
$$B^{-1}M_{\mu}B = 
\begin{pmatrix} 
1+i\mu & 1 & 0 & 0 \\
-1 & 0 & 0 & 0 \\
0 & 0 & 1-i\mu & 1 \\ 
0 & 0 & -1 & 0 
\end{pmatrix}$$
Letting $N_\mu = \begin{pmatrix} 
1+i\mu & 1 \\
-1 & 0 \\
\end{pmatrix}$, we thus have
$$B^{-1}AB = \begin{pmatrix} 
\prod N_{\mu_i} & 0 \\
0 & \overline{\prod N_{\mu_i}}
\end{pmatrix}$$
where the bar denotes entry-wise complex conjugation. Thus the eigenvalues of $A$ are those of $\prod N_{\mu_i}$ plus those of $\overline{\prod N_{\mu_i}}$, which are their complex conjugates.
Moreover, since $N_\mu$ has determinant $1$, so does $\prod N_{\mu_i}$, so its two eigenvalues are inverses of each other.
