Phase factors of eigenstates for a time-dependent Hamiltonian

For a time-dependent Hamiltonian, the Schrödinger equation is given by

$$i\hbar\frac{\partial}{\partial t}|\alpha;t\rangle=H(t)|\alpha;t\rangle,$$

where the physical time-dependent state $|\alpha;t\rangle$ is given by

$$|\alpha;t\rangle = \sum\limits_{n}c_{n}(t)e^{i\theta_{n}(t)}|n;t\rangle$$

and

$$\theta_{n}(t)\equiv -\frac{1}{\hbar}\int_{0}^{t}E_{n}(t')dt'.$$

Here $e^{i\theta_{n}(t)}$ is the phase factor that has been pulled out from the eigenstate-expansion coefficients of $|\alpha;t\rangle$.

Why is $\theta_{n}(t)$ given by

$$\theta_{n}(t)\equiv -\frac{1}{\hbar}\int_{0}^{t}E_{n}(t')dt'?$$

The choice is motivated by trying to generalize the phase factor $e^{-iE_n t/\hbar}$ from the time-independent theory, and more specifically it is the solution to $i\hbar \partial_te^{i\theta_n(t)} = E_n(t)$. This means, in turn, that it is chosen so that the combination $e^{i\theta_n(t)}|n;t⟩$ obeys the Schrödinger equation - so long as you forget that $i\hbar \partial_t |n;t⟩$ might also be nonzero, or in other words that it obeys the Schrödinger equation if you ignore nonadiabatic effects.
However, it is crucial to realize that the choice carries no physical meaning, since that phase is only one of the factors of the actual coefficient: $$\tilde c_n(t) = c_n(t) e^{i\theta_n(t)}.$$ It is only the full $\tilde c_n(t)$ that gives you the probability amplitude contribution of $|n;t⟩$ to $|\alpha;t⟩$, and those are the physical dynamics. The factorization into a known factor $e^{i\theta_n(t)}$ and a variable $c_n(t)$ is for convenience: it allows us to put in explicitly the dynamics which we have already solved, and therefore it allows us to use the $c_n(t)$ to solve for the remaining unknown dynamics.
It is also important to note that the $c_n(t)$ are not trivial quantities, either: they represent non-adiabatic transitions as well as non-adiabatic phases that act in addition to the $e^{i\theta_n(t)}$. However, defining the $\theta_n(t)$ that way allows us to use the $c_n(t)$ to focus explicitly on these (complicated) nonadiabatic effects, and put the (solved) adiabatic dynamics to bed early.
HINT: You can find $\phi_n(t)$ by substituting $|\alpha; t\rangle$ into the first equation (equation of motion). The fields $c_n(t)$ obeye the relation $$i \hbar \frac{d}{dt}c_n(t) + E_n(t)c_n(t) = H(t)c_n(t) \tag{\ast}$$ for all $n \in \mathbb{N}$ (by comparison of coefficients). In the case where $c_n$ is independent on time $(\ast)$ is the well-known eigenvalue equation for stationary quantum mechanics.