An interesting Markov chain with uniform marginals Consider the Markov chain $(\theta_n, \phi_n)$ on $S^1 \times S^1$ constructed in the following way. For $\xi_n$ a sequence of i.i.d. normal random variables and $\kappa > 0$ a fixed number, we set
$$
\theta_{n+1} = \theta_n + \kappa \xi_n\;,\qquad \phi_{n+1} = \arg((3+2\sqrt 2)\cos \phi_n,\sin \phi_n) + \theta_n\;,
$$
where $arg(v)$ denotes the angle of the vector $v$ with $(1,0)$. For every $\kappa > 0$, it has a unique invariant measure $\mu_\kappa$ and it is obvious that its first marginal $\pi_\theta^* \mu_\kappa$ is just Lebesgue measure $\lambda$ for every $\kappa$.
Of course, one has $\lim_{\kappa \to \infty} \mu_\kappa = \lambda \otimes \lambda$, so that $\lim_{\kappa \to \infty} \pi_\phi^*\mu_\kappa =\lambda$. More surprisingly, it follows from an explicit but lengthy calculation that one also has $\lim_{\kappa \to 0} \pi_\phi^*\mu_\kappa =\lambda$. This begs the question: is it true that $\pi_\phi^*\mu_\kappa =\lambda$ for every $\kappa > 0$? I've checked it numerically to rather high accuracy for several values of $\kappa$ and it seems to be the case. Given the simplicity of the statement, it feels like there should be a rather simple argument if it is true, but it eludes me at the moment.
According to computer simulations, the statement seems independent of the specific value $3+2\sqrt 2$ appearing in the definition above although i) it appears to fail when the value is replaced by a negative number and ii) the specific value $3+2\sqrt 2$ simplifies the calculation of the limit $\kappa \to 0$ somewhat.
 A: Since I was requested to elaborate, here goes. First, let's look at the automorphism of the unit circle induced by this mapping (written in the least revealing way). With $z=e^{it}$, as usual, we have $2\cos t=z+z^{-1}, 2i\sin t=z-z^{-1}$, so for positive $3+\sqrt 2$ (I absolutely loved this red herring) the direction is that of $z+\delta z^{-1}$ with $|\delta|<1$. Normalizing, we get 
$$
\Psi(z)=\frac{z+\delta z^{-1}}{|z+\delta z^{-1}|}
$$
and, most importantly,
$$
\Psi(z)^2=\frac{(z+\delta z^{-1})^2}{|z+\delta z^{-1}|^2}=
\frac{z+\delta z^{-1}}{z^{-1}+\delta z}=\frac{z^2+\delta}{1+\delta z^{2}}\,.
$$
Now the full transformation for random variables $U=e^{i\theta}$, $W=e^{i\varphi}$ is $U,W\mapsto ZU, U\Psi(W)$ where $Z=e^{i\xi_n}$ is something smeared a bit and independent of $U,W$. Since it doesn't matter for the convergence to the limiting distribution with what joint distribution to start, we'll start with independent $U,W$ uniformly distributed over the circle.
Claim 1: All four pairs $(\pm U,\pm W)$ have the same distribution at every iteration step.
Proof: It is clearly true in the beginning and, since $\Psi$ is odd, this property is preserved by the mapping.
This alone immediately kills all moments $\mathcal E [U^kW^\ell]$ with $k$ or $\ell$ odd.
Claim 2: If $k,\ell\ge 0$ with $k+\ell>0$, then $\mathcal E [U^{2k}W^{2\ell}]=0$.
Proof: It is clearly true in the beginning. Now just recall that $\Psi^2$ is even analytic, so $\Psi(z)^{2\ell}=\sum_{m\ge 0} a_{\ell,m}z^{2m}$ and, thereby,
$$
\mathcal E [U_{\text{new}}^{2k}W_{\text{new}}^{2\ell}]=
\mathcal E [(ZU)^{2k}(U\Psi(W))^{2\ell}]=
\mathcal E [Z^{2k}]\sum_{m\ge 0}a_{\ell,m}\mathcal E [U^{2k+2\ell}W^{2m}]=0\,.
$$
In particular $\mathcal E[W^{2\ell}]=0$ for $\ell>0$ (and, by conjugation, for $\ell<0$ too) and we are done. 
The above recursion can also be used to find the Fourier coefficients of the limiting distribution with $k<0,\ell>0$ row by row.
