When it comes to decoherency it is often said, that the off-diagonal elements of the density matrix are responsible for coherency. They vanish iff there is no coherency. I spent a lot of time trying to understand this, and I'm still not shure if my explanation is correct. Let's see:
Imagine we'd be able to clone an ensemble discribed by $\rho$ before making measurements on it. This magic technique would allow us to obtain probabilities for measurement results without destroying the original ensemble. In the first run, we would determine the probability for $|\psi_0\rangle$ and $|\psi_1\rangle$ which should be $\langle \psi_0 | \rho | \psi_0 \rangle$ and $\langle\psi_1|\rho|\psi_1\rangle$ respectively. For the next experiments we change the measuring device to measure another observable. Lets assume that the set of orthonormal eigenvectors of this second observable contains an element of the form $$|\sigma\rangle := \alpha |\psi_0\rangle + \beta|\psi_1\rangle.$$ If there is no coherency in the ensemble, we would expect the probability for $|\sigma\rangle$ beeing
$$
|\alpha|^2 \langle\psi_0|\rho|\psi_0\rangle + |\beta|^2 \langle \psi_1 | \rho| \psi_1 \rangle.
$$
But if our experiments show a different probability, one would likely say: "bazinga, there is is some quantum mystery (aka. coherency) going on". Quantum theory says, that the probability for $|\sigma\rangle$ is
$$
\langle \sigma |\rho |\sigma \rangle = |\alpha|^2 \langle\psi_0|\rho|\psi_0\rangle + |\beta|^2 \langle \psi_1 | \rho| \psi_1 \rangle + 2\Re\big( \alpha^*\beta \langle \psi_0|\rho|\psi_1\rangle\big),
$$
which differs from the above iff the off-diagonal element $\langle\psi_0|\rho|\psi_1\rangle$ does not vanish. Summing up we see decoherency if the density matrix has non vanishing off-diagonal elements.
There is another point, I'd like to mention in this context. Lets say our density matrix has all vanishing off-diagonal elements. Does this mean, that there is no coherency in the ensemble? I think the right answer to this question is "no, but...". Lets take a closer look. For a density matrix
$$
\rho = \sum_k p_k |\varphi_k\rangle\langle\varphi_k|,
$$
where $p_k$ is the probability for a system beeing in the state $|\varphi_k\rangle$, vanishing off-diagonal elements only mean that
$$
0 = \sum_k p_k \langle \psi_0 | \varphi_k \rangle \langle \varphi_k | \psi_1 \rangle,
$$
i.e. the elements $p_k \langle \psi_0 | \varphi_k \rangle \langle \varphi_k | \psi_1 \rangle$ sum up to zero. To surely say that there is no coherency, we need all elements $p_k \langle \psi_0 | \varphi_k \rangle \langle \varphi_k | \psi_1 \rangle$ being zero, but this is not what follows from vanishing off-diagonal elements. The key point is, that we have no chance to distinguish between different ensembles which are described by the same density matrix. So for our experiments, an ensemble where all $p_k \langle \psi_0 | \varphi_k \rangle \langle \varphi_k | \psi_1 \rangle$ sum up to zero behaves exactly like an ensembly where all $p_k \langle \psi_0 | \varphi_k \rangle \langle \varphi_k | \psi_1 \rangle$ vanish, i.e. an example without coherent states. Summarizing we can say, that an ensemble described by a density matrix with vanishing off-diagonal elements behaves like if there would be no coherent states. There may still be coherent states in our ensemble, but we have no chance to detect them.
Indirectly.
To start with, set the pesky and useless constants $\hbar/m\omega \mapsto 1 $, to work in natural units. You then recognize the matrix, as Eric emphasizes by identifying to his (5), as written in the number basis of oscillator discrete energy eigenstates $|n\rangle$. Now for a celebrated basis change. Infinity should not bother you.
Consider the state
$$
|q\rangle\equiv \sum_n \psi^*_n(x) |n\rangle ,
$$
where $\psi_n(x)= \langle x|n\rangle $, are (real) Hermite functions, eigenfunctions of the oscillator Schroedinger operator in the x representation-- dot with $\langle m|$: get it?
Dotting with $\langle y|$ and exploiting the completeness of said Hermite functions,
$$\langle y|q\rangle =\sum_n \psi^*_n(x) \psi_n(y)=\delta(x-y),
$$
so you identify this state with $|q\rangle=|x\rangle$.
The pros normally use simpler operators instead of complete sums; that is, they sum the sum to
$$
|x\rangle =\frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2}|0\rangle ,
$$
so that , as you should check,
$$
\hat{x} |x\rangle= \frac{(a+a^\dagger)}{\sqrt 2} ~ \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle= x~\left ( \frac{e^{x^2/2}}{\pi^{1/4}} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle \right ) =x|x\rangle ~.
$$
So the answer to your question is the eigenvalues are all continuous values x, for eigenvectors $e^{x^2/2} e^{-(a^\dagger-\sqrt{2} x)^2/2} |0\rangle /\pi^{1/4} $.
To check your algebra, try the evident unnormalized eigenvector corresponding to x=0, $(1,0,-1/\sqrt{2},0,\sqrt{3/8},0,-\sqrt{5}/4,0,...)$, associated to Hermite numbers.
For a quick crib-sheet on the number basis ops cf Messiah, Quantum Mechanics vI Ch XII ยง5.
Best Answer
The exponential of an operator is generally defined via Taylor series as: $$e^{\hat A} \equiv \sum_{n=0}^\infty \frac{\hat A^n}{n!}$$ Which is completely independent of the matrix representation of $\hat A$ in any basis. Now if $\hat A$ is diagonalizable, as in the case of the hamiltonian, one can show that the above Taylor series definition boils down to exponentiating the diagonal elements of the operator in its eigenbasis (in which it is diagonal), which is what you have written. For the general matrix elements of $e^{\hat A}$ in an arbitrary basis, one can first write the matrix elements in its own eigenbasis, and then use a basis transformation.Explicitly, the spectral decomposition of $e^{\hat A}$ is: $$e^{\hat A} = \sum_k e^{A_k} |A_k \rangle\langle A_k|$$ Now imagine you want the matrix elements in the basis $\{|e_i \rangle \}_i$, this can easily be calculated with the above relation: $$\langle e_i|e^{\hat A}|e_j \rangle = \sum_k e^{A_k} \langle e_i|A_k \rangle\langle A_k|e_j \rangle$$ What the above sum actually is depends on the specific operator in question. However, note the important fact that in general, $\langle e_i|e^{\hat A}|e_j \rangle \neq e^{\langle e_i|\hat A|e_j \rangle}= e^{A_{ij}}$. Matrix elements generally can not be exponentiated for non-diagonal matrices.