What is the difference between $| E, t \rangle $ and $ | E, 0 \rangle $? It is written in Binney and Skinner's The Physics of Quantum Mechanics that states of well defined energy evolve in the following manner. $$| E_n, t \rangle = | E_n, 0 \rangle \mathrm{e}^{- \mathrm{i} E_n t / \hbar},$$ where $E_n$ is the $n$-th energy state and $t$ is time. Now, $$\lvert \mathrm{e}^{- \mathrm{i} E_n t / \hbar}\rvert = 1$$ and $| E_n, t \rangle$ points in the same direction as $| E_n, 0 \rangle$. Thus, $$| E_n, t \rangle = | E_n, 0 \rangle .$$ It is written in the book that the passage of time simply changes the phase of the ket at a rate $E_n / \hbar$. But since a vector is identified with its direction and modulus, both of which are the same for $| E_n, t \rangle$ and $| E_n, 0 \rangle$, both the states must be the same. So, what's the difference? 
Edit:
I think that both the states being the same makes sense physically as well. If we observe the system at some time and find that the energy is $E$ and again observe the system after some time, the energy $E$ is going to be the same; there is no change in the system at all! The phase in the equation is dependent on the time $t$, i.e. it depends on our choice when to make $t = 0$. And since our choice should not affect the physical properties of the state, they must be independent of time.
 A: First off,  $|\Psi(0)\rangle \neq |\Psi(t)\rangle$
While the magnitude of $e^{-iE_nt/\hbar}$ may be $1$, $|\Psi\rangle$ is a complex vector - you don't take the magnitude of $|\Psi\rangle$ until you are calculating the probability of certain states. As time evolves, the complex exponential revolves around the complex plane. We could expand the exponential using Euler's formula to see this:
$$|\Psi(t)\rangle = e^{-iE_nt/\hbar}|\Psi(0)\rangle = \cos(\omega t)|\Psi(0)\rangle - i\sin(\omega t)|\Psi(0)\rangle$$
From here we can see $|\Psi(t)\rangle$ has a real and an imaginary component to its amplitude. This is not the same as $|\Psi(0)\rangle$ which had either an entirely real amplitude or a constant complex amplitude.
Secondly, the example you gave is that of a stationary state which means that there is no measurable time evolution for that state because the time evolution of $\Psi$ does not affect the probability density. If it begins in that state it stays in that state. Physically, there is no detectable difference between $\Psi(0)$ and $\Psi (t)$ but there is a mathematical difference (as I showed) that makes a physical difference under certain circumstances.
The time factor does have a measurable effect when $|\Psi(0)\rangle$ is in a superposition of multiple states and we are seeking to measure an observable which itself does not share an eigenstate with $\Psi$.
Suppose $\Psi$ expresses a superposition of two energy states. Then initially,
$$|\Psi(0)\rangle = c_1|E_1\rangle + c_2|E_2\rangle$$
then
$$|\Psi(t)\rangle = c_1e^{-iE_1t/\hbar}|E_1\rangle + c_2e^{-iE_2t/\hbar}|E_2\rangle$$
$$|\Psi(t)\rangle= e^{-iE_1t/\hbar} \, \big[\, c_1|E_1\rangle + c_2e^{-i(E_2-E_1)t/\hbar}|E_2\rangle \, \big]$$
Now suppose we want to measure the probability of eigenvalue $a_1$ which is an eigenvalue  of measurable $A$. If in the energy basis $a_1$'s eigenstate is a superposition of energy states $|a_1\rangle = \alpha_1|E_1\rangle + \alpha_2|E_2\rangle$ then the probability of measuring $a_1$, $|\langle a_1|\Psi(t)\rangle|^2$, will have a time dependence with the angular frequency $\omega = \frac{E_2 - E_1}{\hbar}$ in it called the Bohr frequency. I'll leave it to you to work out the details of that time dependence.
A: To answer this question, we need a precise definition of what we mean when we say "state."  In many contexts, the word state is taken to mean an element of the Hilbert space $\mathcal H$ underlying the quantum mechanical system in question.  This works well enough, but it's not quite correct for several reasons.
The first reason is that given some $\psi\in \mathcal H$, it is obvious that $\lambda \psi \neq \psi$ for arbitrary $\lambda\in \mathbb C$, but $\lambda \psi$ and $\psi$ are identical in all measurable ways.  For instance, for some linear operator $A$, 
$$\langle A \rangle_{\lambda \psi} = \frac{\langle \lambda \psi, A(\lambda \psi)\rangle}{\Vert\lambda\psi\Vert^2} = \frac{|\lambda|^2\langle \psi,A\psi\rangle}{|\lambda|^2\Vert\psi\Vert^2} = \frac{\langle \psi, A\psi\rangle}{\Vert\psi\Vert^2}=\langle A \rangle_\psi$$
As a result, we should not consider $\lambda\psi$ and $\psi$ to be physically distinct states for the system to exist in, despite the fact that they are inescapably different elements of $\mathcal H$.
The second reason is that there are configurations of a quantum mechanical system which cannot be written as an element of $\mathcal H$ - the so-called mixed states, which are essentially statistical mixtures.  Given two elements $|0\rangle$ and $|1\rangle$ of $\mathcal H$ , there is a difference between the system definitely being in the superposition $|0\rangle+|1\rangle$, and the system being in state $|0\rangle$ or state $|1\rangle$ with 50% probability each.  The latter is a mixed state, and cannot be expressed as  an element of $\mathcal H$ in its own right.

Working physicists resolve these issues in a number of different ways.  


*

*Some simply don't worry about it.  If mixed states aren't an issue, then you can just normalize your states, focus on unitary transformations, and fix the arbitrary global phase.  This works well enough for quite a few people.

*Others are more careful, and note that while $\mathcal H$ underlies a quantum mechanical system, the physical states themselves live in $\mathcal P(\mathcal H) \equiv \mathcal H/\sim$, where $\sim$ is an equivalence relation on $\mathcal H$ whereby $\psi \sim \phi$ if $\psi = \lambda \phi$ for some non-zero complex number $\lambda$.  The elements of $\mathcal P(\mathcal H)$ (called the projective Hilbert space of $\mathcal H$) are equivalence classes of elements of $\mathcal H$ which differ from each other by some non-zero multiplicative factor.

*Others, in particular those who work with mixed states (e.g. quantum computing, quantum information), must be even more careful.  In these cases, one must define a state in terms of a density operator, which can represent both pure states and mixed ones.



To answer your question, in the spirit of (2) and (3) above, $|E,0\rangle$ and $|E,t\rangle$ as specified in the original question are different elements of $\mathcal H$ which correspond to precisely the same physical state, whether you consider the latter to be an element of $\mathcal P(\mathcal H)$ or a density operator.
A: As you already pointed out, these two vectors are not identical
$$ \tag{1}
| E_n, t \rangle = | E_n, 0 \rangle \mathrm{e}^{- \mathrm{i} E_n t / \hbar}
$$
they simply differ by a phase factor. Vectors that differ (globaly) by a phase factor are really physically indistinguishable. That's the reason books often say that system in eigenstate does not evolve (because the only factor that we get from the evolution is the phase factor that is not observable as we already said). However the phase in equation (1) can sometimes "give nontrivial contribution" if we have superposition of two (or more) states.
For example, consider the system which can be described by wavefunction $|\psi,t=0\rangle$ given by:
$$ \tag{2}
|\psi,t=0\rangle = |E_n,0\rangle + |E_m,0\rangle
$$
where $E_m \neq E_n$. It can easily be shown that in later time $t$, the system is in state:
$$ \tag{3}
|\psi,t \rangle = | E_n, 0 \rangle \mathrm{e}^{- \mathrm{i} E_n t / \hbar} +
| E_m, 0 \rangle \mathrm{e}^{- \mathrm{i} E_m t / \hbar}
$$
Now the phase factor $ \mathrm{e}^{- \mathrm{i} E_n t / \hbar}$ is very important and it gives nontrivial contribution to observable quantities because the phase $ \mathrm{e}^{- \mathrm{i} E_n t / \hbar}$ is not global in this case (we have different phases for the two terms in (3)). That is why it is important to keep in mind that we have the phase factor in (1) even though it gives no observable contribution in the case when the wavefunction is the eigenstate.
