Let $O$ be an orthogonal matrix, $O^T O = I$, thus its eigenvalues lie on the unit circle, $\lambda(O)=e^{i\theta}$. Furthermore, assume the form 
$O = X Y$, where both matrices satisfy $X^2 = I$ and $X=X^T$, and also $Y^2=I$ and $Y=Y^T$, thus $\lambda(X),\lambda(Y) \in \{-1,1\}$.
(These matrices do not commute $[X,Y] \ne 0$, $[O,X]\ne 0$.)


What can be said about the eigenvalues of the matrix $P = O + \delta X = X(Y + \delta I)$? Numerically it is possible to verify that up to some value of $\delta$, the *complex eigenvalues* of $P$ still falls on a circle whose radius now depends on $\delta$. Is it possible to prove this?


Any ideas or techniques on how to approach this problem?



About the Eigenvalues of Orthogonal Matrix plus Perturbation
The reason is very simple, actually. Consider the following tautological identity, for some $0\leq z\in \mathbb{R}$,
$$
P^2-2\frac{P+z P^{-1}}{2}P+z I=0.
$$
Denote $\frac{P+z P^{-1}}{2}=B$, then find
$$
P^2-2BP+z I=0.
$$
Note that $[B,P]=0$ and thus if we can diagonalize $P$ and $B$, we can do so simulatneosly, and then for respective eigenvalues we will find
$$
\lambda\_P=\lambda\_B\pm i \sqrt{z-\lambda\_B^2}.
$$
Assume that $\lambda\_B^2\in\mathbb R$ and $\lambda\_B^2<z$, then obviously
$$
|\lambda\_P|=\sqrt z\quad\text{and}\quad \mathrm{Re}\,\lambda\_P=\lambda\_B
$$
If $\lambda\_B^2\in\mathbb R$ and $\lambda\_B^2\geq z$, then 
$$\lambda\_P\in \mathbb{R}.$$


Now the idea is that if we pick $z=1-\delta^2\in \mathbb{R}$, then $B$ is real symmetric and thus diagonalizable with real eigenvalues. Indeed, 
$P^{-1}=\frac{Y-\delta I}{1-\delta^2}X$, and so
$$
2B=P+(1-\delta^2)P^{-1}=XY+YX=O+O^T.
$$


Even more, $B$ does not depend on $\delta$, and its spectrum is just the real part of the spectrum of $O$, so one can write down an explicit form for the spectrum of $P$ in terms of that of $O$. For example, the OP observation that the real parts of the complex eigenvalues are those of $O$ is immediately explained.


Edit: corrected minor errors as pointed out by the OP


