First order wave function in adiabatic approximation If the Hamiltonian is slowly varying in time and suppose the initial state is the n-th eigenstate of the initial Hamiltonian H(0), the adiabatic theorem says that the state will still on the n-th eigenstate of the Hamiltonian H(t) during the evolution, accumulating a dynamic phase and geometric phase.
The above theorem can be shown by considering only the zero-th order wave function, however, what does the first order wave function looks like? 
i.e. when we consider the the transition between the different energy levels, just like in the perturbation method in the time independent perturbation theorem.  
 A: To answer this question is helps to start with the adiabatic portion.
Consider a Hamiltonian $\hat{H}\left(\mathbf{P}\right)$, where $\mathbf{P}(t)$ is some time-varying parameter vector. This Hamiltonian has a set of time-dependent eigenstates $|n\left(t\right)\rangle$ with corresponding energies $E_n\left(t\right)$ which evolve via the Schrodinger's equation $i\hbar \partial_t|\Psi\left(t\right)\rangle = \hat{H}\left(t\right)|\Psi\left(t\right)\rangle$.
If $\mathbf{P}\left(t\right)$ is constant in time and the system starts out in one of the eigenstates $|n\left(0\right)\rangle$ at $t = 0$, the state accumulates a phase factor with time so that $|\Psi\left(t\right)\rangle = e^{-iE_nt/\hbar}|n\left(0\right)\rangle$.
For finite $\partial_t\mathbf{P}\left(t\right)$, the eigenstates and eigenvalues change in time. Still, let us set $|\Psi\left(t\right)\rangle = c_n(t)|n\left(t\right)\rangle$. Putting this into the Schrodinger's equation gives
\begin{equation}
 i\hbar\partial_t \left[ c_n(t)|n\left(t\right)\rangle\right] = E_n\left(t\right)c_n(t)|n\left(t\right)\rangle\,.
\end{equation}
Taking the time derivative of the left-hand side and then applying $\langle n\left(t\right)|$ to both sides, we get
\begin{equation}
 \dot{c}_n\left(t\right) = -\frac{i}{\hbar}c_n\left(t\right)E_n\left(t\right) -c_n\left(t\right)\langle n\left(t\right)|\dot{n}\left(t\right)\rangle  \,.
\end{equation}
This differential equation can be solved using the separation of variables:
\begin{equation}
 c_n\left(t\right) = \exp\left[-\frac{i}{\hbar}\int_0^t E_n\left(t'\right)dt' - \int_0^t\langle n\left(t'\right)|\dot{n}\left(t'\right)\rangle dt'\right]\,.
\end{equation}
It is possible to use the chain-rule to rewrite the second integrand as $\langle n\left(t'\right)|\dot{n}\left(t'\right)\rangle = \langle n\left(\mathbf{P}\right)|\nabla_\mathbf{P}n\left(\mathbf{P}\right)\rangle\cdot \dot{\mathbf{P}}\left(t'\right)$, where $\nabla_\mathbf{P}$ is the gradient in $\mathbf{P}$ space. Since $\dot{\mathbf{P}}\left(t'\right)dt' = d\mathbf{P}$, we have
\begin{equation}
 c_n\left(t\right) =e^{-\frac{i}{\hbar}\int_0^t E_n\left(t'\right)dt'}e^{i\int_{\mathbf{P}_i}^{\mathbf{P}_f}\mathcal{A}_n\cdot d\mathbf{P}}\,,
\end{equation}
where $\mathbf{P}_i = \mathbf{P}_i\left(0\right)$, $\mathbf{P}_f = \mathbf{P}\left(t\right)$, and $\mathcal{A}_n = i\langle n\left(\mathbf{P}\right)|\nabla_\mathbf{P}n\left(\mathbf{P}\right)\rangle$.
So, here you have the dynamic and the geometric (Berry) phase that you mentioned.
Up to now, we assumed that the system remains in the same eigenstate by asserting that the Hamiltonian changes very slowly and the system is always in the $n$th eigenstate. Let us now relax this assumption and write $|n\left(t\right)\rangle \rightarrow |n\left(t\right)\rangle + f\left(\dot{\mathbf{P}}\left(t\right)\right)|\delta n\left(t\right)\rangle$, where $f\left(\dot{\mathbf{P}}\left(t\right)\right)$ is some function. For $\dot{\mathbf{P}}\left(t\right)\rightarrow 0$, the second term vanishes as we restore the adiabatic regime. This means that, to the leading order, $f\left(\dot{\mathbf{P}}\left(t\right)\right)$ is linear in $\dot{\mathbf{P}}\left(t\right)$. By keeping only this linear term and absorbing it into the correction, we write the corrected state $|n\left(t\right)\rangle \rightarrow |n\left(t\right)\rangle + |\delta \left(t\right)\rangle$, where $|\delta \left(t\right)\rangle$ is linear in $\dot{\mathbf{P}}\left(t\right)$.
Above, we showed that the phase of the state does not depend on the rate of change of $\mathbf{P}$. So, we keep $c_n\left(t\right)$ for the state with the correction. Putting $c_n\left(t\right)\left[|n\left(t\right)\rangle + |\delta \left(t\right)\rangle\right]$ into the Schrodinger's equation and dividing by $c_n\left(t\right)$ yields
\begin{align}
   i\hbar\left\{\left( -\frac{i}{\hbar}E_n\left(t\right)+i\mathcal{A}_n\left(t\right) \cdot \dot{\mathbf{P}}\left(t\right) \right) \left[|n\left(t\right)\rangle + |\delta\left(t\right)\rangle\right] 
   + 
   \dot{\mathbf{P}}\left(t\right)\cdot \nabla_\mathbf{P} \left[|n\left(t\right)\rangle + |\delta\left(t\right)\rangle\right]\right\}
  &=
  \hat{H}\left(t\right)\left[|n\left(t\right)\rangle + |\delta\left(t\right)\rangle\right]\,,
\end{align}
where we have replaced the time derivative acting on $|n\left(t\right)\rangle + |\delta\left(t\right)\rangle$ by $\dot{\mathbf{P}}\left(t\right)\cdot \nabla_\mathbf{P}$. 
There are a few simplifications we make here. First, the $E_n\left(t\right)|n\left(t\right)\rangle$ term on the left gets cancelled by $\hat{H}\left(t\right)|n\left(t\right)\rangle$ on the right. Also, $\dot{\mathbf{P}}\left(t\right) |\delta\left(t\right)\rangle$ and $\dot{\mathbf{P}}\left(t\right) \cdot \nabla_\mathbf{P}|\delta\left(t\right)\rangle$ are quadratic in $\dot{\mathbf{P}}\left(t\right)$, so we eliminate them. The remaining terms are
\begin{align}
  i \hbar\dot{\mathbf{P}}\left(t\right)\cdot \left[
   i \mathcal{A}_n\left(t\right)  
    +
     \nabla_\mathbf{P} \right]
     |n\left(t\right)\rangle
  &=
  \left[\hat{H}\left(t\right) - E_n \left(t\right)\right]|\delta\left(t\right)\rangle \,.
 \end{align}
Finally, we apply $\langle m\left(t\right)|$ to both sides of this expression to get
\begin{equation}
 |\delta\left(t\right)\rangle = i\hbar 
 \sum_{m\neq n}|m\left(t\right)\rangle\frac{\langle m\left(t\right) |\partial_t n\left(t\right)\rangle }{E_m\left(t\right)-E_n\left(t\right)}\,,
 \label{eqn:delta_3}
\end{equation}
where we have recombined $\dot{\mathbf{P}}\left(t\right)\cdot \nabla_\mathbf{P}$ into $\partial_t$.
Using $\langle m\left(t\right)|\partial_t \left[\hat{H}\left(t\right)|n\left(t\right)\rangle\right] = \langle m\left(t\right)|\partial_t\hat{H}\left(t\right) |n\left(t\right)\rangle + \langle m\left(t\right)|E_m\left(t\right)|\partial_tn\left(t\right)\rangle = E_n\left(t\right)\langle m\left(t\right)|\partial_t n\left(t\right)\rangle$, we write the time-dependent state as
\begin{equation}
 |\Psi\left(t\right)\rangle =c_n\left(t\right)\left[|n\left(t\right)\rangle 
 -
 i\hbar 
 \sum_{m\neq n}|m\left(t\right)\rangle\frac{\langle m\left(t\right) |\partial_t \hat{H}\left(t\right)|n\left(t\right)\rangle }{\left[E_m\left(t\right)-E_n\left(t\right)\right]^2}\right]
 \,.
\end{equation}
