Does $U^T U$ have real eigenvectors if $U$ is unitary? I was reading a paper by Kraus and Cirac on general two qubit gates. A crucial step is to decompose an arbitrary unitary operator acting on two qubit space into a special form:
$$ U_{AB} = U_A\otimes U_BU_d V_A\otimes V_B\tag{A1}$$
where $U_A,U_B,V_A,V_B$ are local unitary operators that only act on subspace $A$($B$). $U_d$ is a unitary operator of the form $$\exp[-i(\alpha \sigma_x \otimes \sigma_x + \beta \sigma_y \otimes \sigma_y + \gamma \sigma_z \otimes \sigma_z)]. $$
The proof of this decomposition is included in the appendix of the paper. And I got stuck with one step in the text between eq. (A7) and eq. (A8) claiming that for a general unitary operator $U$, $U^T U$ has  real eigenvectors. I couldn't figure out why this is the case. Is there any essential symmetry leading to this result? How do I prove that a matrix have all real eigenvectors? 
 A: Maybe I'm an idiot, but how about this: If $U$ is unitary then so is $U^T$. As the product of two unitaries is unitary, we have that    $U^TU$ is a complex-symmetric unitary matrix. Let's write $C\equiv U^TU = X+iY$ where $X$, $Y$ are real symmetric. Then $C^\dagger C=X^2+Y^2 +i[X,Y]$ is the identity and therefore real. This means that  $X$,$Y$ commute. A pair of real commuting symmetric matrices  can be simultaneously diagonalized by an orthogonal transformation, so $U^TU$ can be diagonalized by a real  orthogonal  transformation. I.e the eigenvectors (but not the eigenvalues) are real.
I'm sure Ill find my error in a few minutes!
A: To fix the Mathematica code in Emilio's answer I believe it should be CircularUnitaryMatrixDistribution instead of GaussianUnitaryMatrixDistribution.

GaussianUnitaryMatrixDistribution[σ,n], also referred to as a Gaussian unitary ensemble (GUE), represents a statistical distribution over the $n\times n$ complex Hermitian matrices, namely square complex matrices $A$ that satisfy $A=A^\dagger$, where $A^\dagger$ denotes the conjugate transpose of $A$.
CircularUnitaryMatrixDistribution[n], also referred to as the circular unitary ensemble (CUE), represents a statistical distribution over the $n\times n$ unitary complex matrices, namely complex square matrices $A$ satisfying $A.A^\dagger=I_n$.

And the eigenvectors are real for every test I tried.
A: It does look to me like there's something fishy going on. Quoting from the proof of Lemma 2:

Let us denote by $\{| \Psi_k\rangle\}$ the eigenstates of $U^TU$, where $U^T$ denotes the transpose of $U$ in the magic basis and $e^{2i\epsilon_k}$ are the corresponding eigenvalues. Note that the eigenvectors of the symmetric operator $U^TU$ are orthonormal and real, except for global phases.

I agree with you that it is pretty natural to read this text as implying that it is using the claim

the eigenvectors of any symmetric operator $A$ are orthonormal and real

as applied to $A=U^TU$. This is definitely false: this page provides the explicit example
$$
A 
= 
\begin{pmatrix} 2i & 1 \\ 1 & 0 \end{pmatrix}
=
\begin{pmatrix} i & 1 \\ 1 & 0 \end{pmatrix}
\begin{pmatrix} i & 1 \\ 0 & i \end{pmatrix}
\begin{pmatrix} i & 1 \\ 1 & 0 \end{pmatrix}^{-1}
$$
of a complex symmetric matrix $A = A^T$ that, being equivalent to a Jordan block, is not diagonalizable.
However, upon a closer inspection, the claim as stated is

the eigenvectors of the symmetric operator $U^TU$ are orthonormal and real except for global phases

and this is the result that actually gets used, i.e. relying on the additional structure of $A=U^TU$ with $U$ being unitary. To test this, here is some short Mathematica code to pull up a random unitary (from the Haar-invariant ensemble) and look at its eigenvectors:
Block[{U},
 U = RandomVariate[CircularUnitaryMatrixDistribution[4]];
 MatrixForm /@ (
   Chop[E^(-I Arg[#[[1]]]) #] & /@ Eigenvectors[
     Transpose[U].U
     ]
   )
 ]

The claim does appear to be true, at least on the level of the numerical evidence provided by the code above ─ it seems to invariably produce real-valued eigenvectors.
(A previous, incorrect version of this post claimed this was not the case. The mistake was using GaussianUnitaryMatrixDistribution and thinking that it would produce a gaussian distribution over the set of unitary matrices, without a careful enough use of the documentation or checks on the output to verify that this was the case. CircularUnitaryMatrixDistribution is the correct ensemble to use here. Apologies for the misinformation.)
A: I begin to think that this statement is actually valid. I will first illustrate a proof I found in this brief introduction to KAK decomposition (KAK decomposition is just a more formal way of saying the decomposition I referred to in the question). Then I will present some numerics I did with matlab. The results suggest that such a statement is valid, contrary to what Emilio found with Mathematica. 
This might not be the most concise way of proving it, but here it goes. 
The key idea is that for a general unitary matrix U, we can always found one singular value decomposition of it such that $U = Q_L e^{i\Theta} Q_R^T$, where $e^{i\Theta}$ is a diagonal unitary matrix, and $Q_R,Q_L$ are orthogonal matrices. If this is the case, our statement goes without saying. Proving this relies on one lemma and Eckart-Young theorem. 
Eckart-Young theorem 
Two real matrices A and B have simultaneous singular value decomposition (i.e. $A = Q_LD_1 Q_R^T, B = Q_LD_2 Q_R^T,D_1$ and $D_2$ being real,$Q_L$ and $Q_R$ being orthogonal) iff $AB^\dagger$ and $A^\dagger B$ are symmetric matrices .
Lemma
Define the real part and imaginary part of a unitary matrix U as $U_R = \frac{U+U^*}{2}$ and $U_I = \frac{U-U^*}{2i}$ ($U^*$ being the complex conjugate of U), we have $U_I U_R^T$ and $U_I^T U_R$ real and symmetric. 
The proof can be found in the reference. The theorem and the lemma imply that we can do singular value decomposition to the real part and imaginary part of U simultaneously, thus U has the singular value decomposition as $U = Q_L e^{i\Theta} Q_R^T$. One can then derive that both $U^T U$ and $U U^T$ has real eigenvectors. 
The numerics are presented below: 
I use matlab code to generate random unitary matrices, the function is adopted from this quora answer. 
function [U] = Unitary(n)
% generate a random complex matrix
X = complex(rand(n), rand(n))/sqrt(2);
$ factorize the matrix
[Q,R] = qr(X);
R = diag(diag(R)./abs(diag(R)));
% unitary matrix
U = Q*R;
end

The test I did then is
randomU = Unitary(4);
U1 = transpose(randomU)*random(U);
U2 = randomU*transpose(randomU);
[eigvec1,~] = eig(U1);
[eigvec2,~] = eig(U2);
disp(eigvec1);
disp(eigvec2);

All the eigenvectors have zero imaginary part to the precision limit, and some of the results are:

