Why do degenerate normal modes appear as complex conjugate pairs? Can anyone prove this? 
THE DETAILS
Suppose we have a system with n components, (i.e. $|~{\psi}(t)\rangle=\sum \langle x_i|~\psi(t)\rangle ~|~x_i\rangle$) where our equation of motion is described by 
 $$\frac{d^2}{dt^2}|~\psi(t)\rangle=\hat{M}^{-1}\hat{K}~|~\psi(t)\rangle$$
And say we are looking for solutions $|~\omega\rangle$ such that $\frac{d}{dt}|~\omega\rangle=\omega|~\omega\rangle$, and thus our equation simplifies to 
 $$\hat{M}^{-1}\hat{K}~|~\omega\rangle=\omega^2|~\omega\rangle\tag{1}$$
(Note that $\hat{M}^{-1}\hat{K}$ is Hermitian). Is there any way to prove that if $|~\omega_1\rangle,|~\omega_2\rangle$ are linearly independent and both satisfy (1), then $|~\omega_1\rangle=|~\omega_2\rangle^*$, i.e. the components of one vector are the complex conjugates of the components of the other vector?
I ask this because this keeps on popping up in our class study of coupled oscillators and I want to know why this is happening
Thanks!
 A: Just for some closure, I'll summarize Qmechanic's answer:
Since $T$ is self-adjoint, note that if
$$ Tv=\lambda v$$
Then, the corresponding equation in the dual space holds:
$$T^*v^*=\lambda^*v^*=\lambda v^*\tag{$\lambda\in \mathbb{R}$ because $T$ is self-adjoint}$$
But $T^*=T$, so $Tv^*=\lambda v^*$. So for each eigenvector of $T$, its complex conjugate is also an eigenvector of $T$ with the same eigenvalue.
A: I'd like to close off on David's Answer and Joshphysic's comment: 

You should still be careful. You have only shown that given an eigenvector, it's complex conjugate is also an eigenvector with the same eigenvalue. You have not shown that degenerate, linearly independent eigenvectors must be complex conjugates of one another. In fact, I'm skeptical that's even true; consider, for example, the eigenvectors $v+v^∗$ and $v−v^∗$.

because together they answer the OP's question fully. Namely, the OP's result is NOT generally true. What is true, as David has proven, is that degenerate eigenvectors can be arranged to be complex conjugate pairs by convention. So if $v$ is an eigenvector corresponding to eigenvalue $\lambda$, then so is $v^*$. The space $\ker(T - \lambda \,\text{id})$  of eigenvectors corresponding to eigenvalue $\lambda$ is always a linear space, and you can use whatever linearly independent basis to span it that takes your fancy, so in general there will be no special relationship between two linearly independent members of this space. Actually David's proof is often used in optical waveguide theory: as joshphysics proved you can take any member $v\in\ker(T - \lambda \,\text{id})$  and form $v+v^∗$ and $(v−v^∗)/i$ and these eigenvectors are both real. So David has proven that $\ker(T - \lambda \,\text{id})$ can always be spanned by purely real eigenvectors. Thus bound fibre modes (or radiation modes from the fibre's continuous spectrum, for that matter) are always taken, by convention, to be real (if the fibre's material is lossless they are eigenfunctions of an Hermitian operator). Otherwise put: the wavefronts for eigenmodes are transverse planes advancing steadily along the fibre: the phase of a conventionally chosen mode is constant over a transverse plane and advances in that direction through the phase factor $e^{i\,\beta\,z}$ where $\beta$ is the eigenvalue (called the propagation constant) and $z$ is the translationally-invariant direction. Quantum mechanics is simply four dimensional waveguide theory, with $t$ being the translationally-invariant direction, hence the phase of the conventionally chosen conventionally chosen eigenfuncgtion is constant over a "transverse plane" (a section of constant time) and advances in that direction through the phase factor $e^{-i\,E\,t/\hbar}$, where $E$ is the energy eigenvalue. 
Another convention about eigenvectors in QM (and optical waveguide theory) is that the basis of $\ker(T - \lambda \,\text{id})$  is chosen, by the Gramm-Schmidt orthogonalisation procedure, so that it is orthogonal. Thus, by introducing a convention, we can say "eigenvectors are always orthogonal", but we keep in mind that for many-dimensional linear spaces of degenerate eigenvectors we have brought a particular choice to bear to "make" the statement true. Any two randomly chosen linearly dependent members of $\ker(T - \lambda \,\text{id})$ will not in general be orthogonal.
A: There are a few problems in your problem statement and your own answer. For instance, the matrix $M^{-1}K$ is not Hermitian if $M$ and $K$ are to be interpreted as the usual mass matrix and spring stiffness matrix. Also, it is not clear to me what OP means. Taking OP's question literally, then the answer is no. Degenerate normal modes not always come in complex conjugation pairs in classical oscillators. A counter-example (essentially the example given by @V.Moretti) is a set of two independent oscillators with equal frequency.
Still, let me try to sum up what have been said and make them more precise.
To avoid mathematical non-sense, we should stick to the simplest case: a finite collection of mechanical oscillators coupled to each other through springs.
The equation of motion is given by
$$
M\ddot{q} = Kq.
$$
Here $q$ is the (generalized) coordinates, which is a real, $n$ dimensional column vector. $M$ is the $n\times n$ mass matrix, and $K$ is the standard, $n\times n$ spring stiffness matrix. Both $M$ and $K$ are real, symmetric, and positive-definite. The positive-definiteness is required by physics as both the kinetic and potential energy must be bounded from below. (Strictly speaking, semi-positive-definiteness is allowed. Here I assume the positive-definiteness to make my life a bit easier.)
As a side note, we see that the matrix $M^{-1}K$ is not symmetric since $M$ and $K$ don't commute in general.
Now we make a change of coordinates by defining the following:
$$
x = M^{1/2}q.
$$
The meaning of $M^{1/2}$ is the following: since $M$ is a real, symmetric, positive-definite matrix, it can be brought to the diagonal form: $M=V\Lambda_M V^t$.$V$ is a $n\times n$ real orthogonal matrix, and $\Lambda_M = diag(m_1, m_2, m_3,\cdots m_n)$ is a $n\times n$ diagonal matrix. In particular, $m_i>0$ for all $1\le i\le n$. Therefore, we can define a new matrix, $M^{1/2}$, as
$$
M^{1/2} = V\left(\begin{array}{cccc}
\sqrt{m_1} & 0 & \cdots & 0\\
0 & \sqrt{m_2} & \cdots & 0\\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & \sqrt{m_n}
\end{array}\right)V^t.
$$
It is easy to check that $M^{1/2}M^{1/2}=M$, and it is real, symmetric, and positive-definite by construction. Likewise, $M^{-1/2}$ can be defined, and it is a real, symmetric, and positive-definite matrix, and $M^{-1/2}M^{1/2}=1_n$.
The equation of motion can then be written in terms of the new coordinates:
$$
\ddot{x} = M^{-1/2}KM^{-1/2}x \equiv Hx.
$$
Here $H = M^{-1/2}KM^{-1/2}$. We see that $H$ is real and symmetric. More importantly, it is positive-definite.
The rest is then clear: since $H$ is a real, symmetric, and positive-definite matrix, it can be diagonalized by a real, orthogonal transformation. The eigenvalues of $H$ correspond to the natural frequencies squared of the system. Again, it is important to have a positive-definite $H$; otherwise, the natural frequencies could be imaginary, which is unphysical in our current model. Most importantly, we see clearly from the analysis that, there is generically no degeneracy in the system. And, if there is any degeneracy, we shouldn't expect them to be related by complex-conjugation.
Still, the classical oscillators system possesses time-reversal symmetry. Thus, it is expected that is $x(t)$ is a solution and then $x(-t)$ is a solution. This is related to the fact that, for each normal mode, there are two solutions, $e^{i\omega t}$ and $e^{-i\omega t}$ (or, if you like, $\cos(\omega t)$ and $\sin(\omega t)$). This can be seen for a single oscillator. But this is not to be confused with the degeneracy of normal modes.
