\subsection{Determining the frequency shift with time dependent perturbation theory}
\label{ss:tdp}
Our goal is to solve the Schrodinger equation, defined by the time dependent Hamiltonian (Equation~\ref{eq:ham}) using standard time
dependent perturbation theory.  But before we address that problem let's try something simpler, namely the
case with fixed $\beta_z$ as in ~\ref{eq:pitchhamiltonian}. 
We already know the exact solution and we can expand the exact frequency shift to second order in the small parameter $\Omega$.
(Recall that $\Omega$ is proportional to the velocity in the direction of the magnetic field $\beta_z$.)
We can alternatively solve the Schrodinger equation for ~\ref{eq:pitchhamiltonian} to second order in $\Omega$ using perturbation theory\cite{griffiths}[p. 302]
and check for consistency.
({\it Following Griffiths QM, 1st edition, section 9.1.2, page 302, problem 9.4}).
%For the dipole in a magnetic field, $E_a= -\hbar\frac{\omega_0}{2}$ ,$E_b= +\hbar\frac{\omega_0}{2}$ and $V=-\hbar\frac{\Omega}{2}e^{i\omega t}$, where
%$B_z = k\omega_0$ and $B_x +i B_y = k\Omega$. Then
\def\calH{{\cal H}}
\begin{eqnarray*}
{\cal H} &=& -\frac{\hbar}{2}\begin{pmatrix}\omega_0 &\Omega e^{i\omega t} \cr \Omega e^{-i\omega t} & -\omega_0\end{pmatrix}\\
&=&{\cal H}_0+{\cal H}^\prime\\
&=&-\frac{\hbar}{2}\begin{pmatrix} \omega_0 &0 \cr 0 & -\omega_0\end{pmatrix}
-\frac{\hbar}{2}\begin{pmatrix} 0 & \Omega e^{i\omega t} \cr \Omega e^{-i\omega t} & 0\end{pmatrix}
\end{eqnarray*}
The Schrodinger equation is
\begin{eqnarray*}
i\hbar\partder{}{t}\psi(t)&=&\calH\psi(t)\\
\end{eqnarray*}

$\psi_a$ and $\psi_b$ are eigenkets of the unperturbed Hamiltonian with eigenvalues $E_a=-\frac{\hbar}{2}\omega_0=-E_b$
\begin{eqnarray*}
\calH_0\psi_a &=& E_a\psi_a,\ \ \ \calH_0\psi_b = E_b\psi_b
\end{eqnarray*}
and since $\psi_a$ and $\psi_b$ form a complete set, the general solution for the exact Hamiltonian is a linear combination of the two
\begin{eqnarray}
\psi(t)&=& c_a(t)\psi_a e^{-iE_z t/\hbar}+c_b(t)\psi_b e^{-iE_b t/\hbar}\nonumber\\
&=& c_a(t)\psi_a e^{i\omega_0 t/2}+c_b(t)\psi_b e^{-i\omega_0 t/2}\label{eq:generalsolution}
%\expect{\psi}{H}{\psi}&=&\begin{pmatrix}\expect{\psi_a}{\calH}{\psi_a} & \expect{\psi_a}{\calH}{\psi_b} \cr \expect{\psi_b}{\calH}{\psi_a} & 
%\expect{\psi_b}{\calH}{\psi_b} 
%\end{pmatrix}=\begin{pmatrix}E_a & V \cr V^* & E_b\end{pmatrix}
\end{eqnarray} 
\def\op{{\eta}}
\def\opt2{{\eta t/2}}
Substitution of ~\ref{eq:generalsolution} into the Schrodinger equation yields, after some manipulation 
\begin{eqnarray*}
\dot c_a&=&-\frac{i}{\hbar}c_b\calH^\prime_{ab}e^{-i\omega_0 t}\\
&=&-ic_b\frac{\Omega}{2} e^{i\omega t}e^{-i\omega_0 t}
\end{eqnarray*}
or
\begin{eqnarray}
i\dot c_a&=& c_b \frac{\Omega}{2} e^{i\op t}\nonumber\\
i\dot c_b&=& c_a \frac{\Omega}{2} e^{-i\op t}\label{eq:tdpde}
\end{eqnarray}
where $\eta=\omega-\omega_0$.
If the muon is initially (at $t=0$) in an eigenstate of $\sigma_x$ the zeroth order solution to the Schrodinger equation is,
\begin{eqnarray*}
c_a(0)&=& \frac{1}{\sqrt{2}}\\
c_b(0)&=& \frac{1}{\sqrt{2}}
\end{eqnarray*}
In order to compute the first order correction we substitute the zeroth order solution into the right hand side of ~\ref{eq:tdpde}
and integrate
\begin{eqnarray*}
c^1_a(t)&=&\frac{-i}{\sqrt{2}}\int_0^t Ve^{i\op t}= \frac{-i}{\sqrt{2}}\frac{V}{i\op}(e^{i\op t}-1)=-\frac{1}{\sqrt{2}}\frac{V}{\op}(e^{i\op t}-1)\\
c^1_b(t)&=&\frac{-i}{\sqrt{2}}\int_0^t Ve^{-i\op t}= \frac{-i}{\sqrt{2}}\frac{V}{-i\op}(e^{-i\op t}-1)= \frac{1}{\sqrt{2}}\frac{V}{\op}(e^{-i\op t}-1)
\end{eqnarray*}
where $V=\frac{\Omega}{2}$.
Now we can write $c_a(t), c_b(t)$ to first order
\begin{eqnarray*}
c_a(t)&\sim&c_a^0(t)+c^1_a(t)=\frac{1}{\sqrt{2}}- \oneonr2\frac{V}{\op}(e^{i\op t}-1)\\
c_b(t)&\sim&c_b^0(t)+c^1_b(t)=\oneonr2+ \oneonr2\frac{V}{\op}(e^{-i\op t}-1)
\end{eqnarray*}
To compute the second order term
substitue $c^1_a(t)$, $c^1_b(t)$ into \ref{eq:tdpde} and integrate again
\begin{eqnarray}
i\dot c^2_a&=& c^1_b Ve^{-i\op t} =\frac{1}{\sqrt{2}}\frac{V}{\op}(e^{-i\op t}-1) Ve^{i\op t}\nonumber\\
 c_a^2 &=& -\frac{i}{\sqrt{2}}\frac{V}{\op}\int_0^t(e^{-i\op \tp}-1) Ve^{i\op \tp}d\tp\nonumber\\
 c_a^2 &=& -\frac{i}{\sqrt{2}}\frac{V}{\op}\int_0^t(1-e^{i\op \tp}) Vd\tp\nonumber\\
 c_a^2 &=& -\frac{i}{\sqrt{2}}\frac{V}{\op}(t -\frac{e^{i\op t}-1}{i\op}) V\nonumber\\
&=& \frac{1}{\sqrt{2}}\frac{V^2}{\op}\left(\frac{e^{i\op t}-1}{\op}-it\right)\label{eq:itterm}\\
c_b^2&=&\frac{1}{\sqrt{2}}\frac{V^2}{\op}\left(\frac{e^{-i\op t}-1}{\op}+it\right)\nonumber
\end{eqnarray}
Following Griffith's notation where $a(t) = c_a(t)e^{-i\omega_0 t/2}$ and $b(t) = c_b(t)e^{i\omega_0 t/2}$ we can write
\begin{eqnarray}
a(t)&\sim& \frac{1}{\sqrt{2}} e^{i\omega_0 t/2 }\left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1-i\op t\right)+\ldots\right)\label{eq:itfreqshift}
\end{eqnarray}
Rearranging 
\begin{eqnarray*}
a(t)&\sim& \frac{1}{\sqrt{2}} e^{i\omega_0 t/2 }\left([1 -i\frac{V^2}{\eta}t]-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\\
\end{eqnarray*}
Then to order $V^2$
\begin{eqnarray}
a(t)&\sim&\frac{1}{\sqrt{2}} e^{i\omega_0 t/2 }[1 -i\frac{V^2}{\eta}t]\left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\nonumber\\
&\sim&\frac{1}{\sqrt{2}} e^{i\omega_0 t/2 }e^{ -i\frac{V^2}{\eta}t}\left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\nonumber\\
&\sim&\frac{1}{\sqrt{2}} e^{i(\omega_0 -\omega)t/2 }e^{i\omega t/2}e^{ -i\frac{V^2}{\eta}t}\left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\nonumber\\
&\sim&\frac{1}{\sqrt{2}} e^{-i\lambda t/2 }e^{i\omega t/2}\left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\label{eq:tdptest}
\end{eqnarray}
%This is a good approximation as long as $\frac{V^2}{\op}t\ll 1$.
where $\lambda/2 = \frac{\op}{2}+\frac{V^2}{\op}$.
Finally with the substitutios  $\eta=\omega-\omega_0$ and $V=\Omega/2$
\begin{eqnarray}
a(t)&\sim&\frac{1}{\sqrt{2}} e^{-i(\omega-\omega_0+\frac{\Omega^2}{2(\omega-\omega_0)}) t/2 }e^{i\omega t/2}\left(1-\frac{\Omega}{4(\omega-\omega_0)}(e^{i(\omega-\omega_0) t} -1)+ \frac{\Omega^2}{4(\omega-\omega_0)^2}\left(e^{i(\omega-\omega_0) t}-1\right)+\ldots\right)\nonumber\\
&\sim&\frac{1}{\sqrt{2}} e^{-i(\omega-\omega_0+\frac{\Omega^2}{2(\omega-\omega_0)}) t/2 }e^{i\omega t/2}\left(1-\left(\frac{\Omega}{(\omega-\omega_0)}-\frac{\Omega^2}{2(\omega-\omega_0)^2}\right)i e^{i(\omega-\omega_0) t/2}\sin((\omega-\omega_0)t/2)+\ldots\right)\nonumber\\
&\sim&\frac{1}{\sqrt{2}} e^{-i(\frac{\Omega^2}{2(\omega-\omega_0)}) t/2 }e^{i\omega t/2}\left(e^{-i(\omega-\omega_0)t/2}-\left(\frac{\Omega}{(\omega-\omega_0)}-\frac{\Omega^2}{2(\omega-\omega_0)^2}\right)i\sin((\omega-\omega_0)t/2)+\ldots\right)\label{eq:tdpt3}\\
&\sim&\frac{1}{\sqrt{2}}e^{i\omega t/2}\left(e^{-i(\omega-\omega_0+\frac{\Omega^2}{2(\omega-\omega_0)})t/2}-\left(\frac{\Omega}{(\omega-\omega_0)}-\frac{\Omega^2}{2(\omega-\omega_0)^2}\right)i\sin((\omega-\omega_0)t/2)+\ldots\right)\label{eq:tdpt4}
\end{eqnarray}
To get from~\ref{eq:tdpt3} to~\ref{eq:tdpt4} we drop terms of order $>\Omega^2$.

Let's compare~\ref{eq:tdpt4} to the exact solution ~\ref{eq:umatrix} with initial conditions $a_0=b_0=\onroot2$, expanded to second order in $\Omega$.
\begin{eqnarray*}
a&=&\onroot2\left(\cos(\omega^\prime t/2)-\frac{i(\omega-\omega_0)}{\omega^\prime}\sin(\omega^\prime t/2) -\frac{i\Omega}{\omega^\prime}\sin(\omega^\prime t/2)\right)e^{i\omega t/2}
\end{eqnarray*}
Note $\omega^\prime=\sqrt{(\omega-\omega_0)^2+\Omega^2}\sim \omega-\omega_0 + \frac{\Omega^2}{2(\omega-\omega_0)}$. With that approximation
\begin{eqnarray}
a&\sim&\onroot2\left(\cos(\omega^\prime t/2)-i(1-\frac{\Omega^2}{2(\omega-\omega_0)^2})\sin(\omega^\prime t/2) -\frac{i\Omega}{\omega-\omega_0}\sin(\omega^\prime t/2)\right)e^{i\omega t/2}\nonumber\\
&\sim&\onroot2\left(e^{-i(\omega-\omega_0 +\frac{\Omega^2}{2(\omega-\omega_0)}) t/2}+i(\frac{\Omega^2}{2(\omega-\omega_0)^2})
\sin((\omega-\omega_0) t/2) -\frac{i\Omega}{\omega-\omega_0}\sin((\omega-\omega_0) t/2)\right)e^{i\omega t/2}\label{eq:exact1}
\end{eqnarray}
Sure enough the exact solution expanded to second order in $\Omega$ (Equation~\ref{eq:exact1}) is consistent with the perturbative solution (Equation~\ref{eq:tdpt4}).
The important take away is that the second order frequency shift is the coefficient of the term linear in imaginary time ($it$). (See Equation~\ref{eq:itterm} -~\ref{eq:tdptest})

\def\opp{{\lambda}}

\begin{comment}
\begin{eqnarray*}
a(t)&\rightarrow &\frac{1}{\sqrt{2}} e^{-i(\lambda/2)t}e^{-i\omega t/2} \left(1-\frac{V}{\op}(e^{i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{i\op t}-1\right)+\ldots\right)\\
b(t)&\rightarrow& \frac{1}{\sqrt{2}}e^{i(\lambda/2)t }e^{i\omega t/2}\left(1+\frac{V}{\op}(e^{-i\op t} -1)+ \frac{V^2}{\op^2}\left(e^{-i\op t}-1 \right)+\ldots\right)
\end{eqnarray*}
To determine the polarization, as given in Equation~\ref{eq:polarization} we evaluate the real and imaginary parts of $a^* b$ and $|a|^2-|b|^2$.
\begin{eqnarray*}
a^* b&\sim& e^{i\lambda t}\left(1-\frac{V^2}{\omega_0^2}\left(-4e^{-i\op t}\sin^2\opt2\right) - \frac{V^2}{\omega_0^2}4i\sin\omega t/2 e^{-i\opt2})\right)
\end{eqnarray*}

\begin{eqnarray*}
\vev{s_x}=\Re(a^*b)&=& \cos(\omega) t(1+{\cal O}\frac{V^2}{\op^2}\\
\vev{s_y}=\Im(a^*b)&=& \sin(\omega) t(1+{\cal O}\frac{V^2}{\op^2}\\
\vev{s_z}=|a|^2-|b|^2&=&2\frac{V}{\op}(\cos\op t -1)
\end{eqnarray*}
Finally 
\begin{eqnarray*}
\bmbetahat\cdot\bfs&=& \beta_{xy}(\cos\omega t \vev{s_x} + \sin\omega t\vev{s_y} ) + \beta_z \vev{s_z}\\
&=& \beta_{xy}\cos\opp t +2\beta_z\frac{V}{\op}(\cos\op t -1)\\
&\rightarrow&\beta_{xy}\cos\opp t
\end{eqnarray*}
and in the limit where $\beta_z V/\op\ll 1$, the polarization projected along the direction of motion oscillates with frequency $\opp$,
where 
\end{comment}
The frequency derived from the perturbation theory
\begin{equation}
\lambda =  \op+ 2\frac{V^2}{\op} = (\omega_0-\omega)+\onehalf\frac{\Omega^2}{(\omega_0-\omega)}
\end{equation}
and the exact frequency (Equation~\ref{eq:exactfreq})
\begin{equation}
\omega^\prime = \left((\omega_0-\omega)^2+\Omega^2\right)^{1/2}\sim (\omega_0-\omega) +\onehalf\frac{\Omega^2}{(\omega_0-\omega)^2}
\end{equation}
%We see that it is the term proportional to $it$ that appears in the second order correction (Equation~\ref{eq:itterm}) that corresponds to the frequency shift
%that is defined by Equation~\ref{eq:itfreqshift}.

\begin{comment}  %Below is the evaluation of higher order terms
Substitue $c^2_a(t)$, $c^2_b(t)$ into \ref{eq:tdpde} and integrate again
\begin{eqnarray*}
i\dot c^3_a&=& c^2_b Ve^{-i\omega_0 t} = \frac{1}{\sqrt{2}}\frac{V^2}{\omega_0}\left(\frac{e^{i\omega_0 t}-1}{\omega_0}-it\right)                   Ve^{-i\omega_0 t}\\
 c_a^3 &=& -\frac{i}{\sqrt{2}}\frac{V^3}{\omega_0}\int_0^t\left(\frac{e^{i\omega_0 \tp}-1}{\omega_0}-i\tp\right) e^{-i\omega_0 \tp} d\tp\\
 &=& -\frac{i}{\sqrt{2}}\frac{V^3}{\omega_0}\left(\frac{e^{-i\omega_0 t}-1}{i\omega^2_0}+\frac{t}{\omega_0}-i\int \tp e^{-i\omega_0 \tp} d\tp\right)\\
 &=& \frac{i}{\sqrt{2}}\frac{V^3}{\omega_0}\left(\frac{e^{-i\omega_0  t}-1}{i\omega^2_0}+\frac{t}{\omega_0}-i  \left(\frac{e^{-i\omega_0 t}-1}{\omega^2_0}+i\frac{t e^{-i\omega_0 t}}{\omega_0}\right)        \right)\\
 &=& \frac{i}{\sqrt{2}}\frac{V^3}{\omega_0}\left(-2i\frac{e^{-i\omega_0  t}-1}{\omega^2_0}+\frac{t}{\omega_0}(1 +  e^{-i\omega_0 t})\right)      \\
 &=& \frac{1}{\sqrt{2}}\frac{V^3}{\omega_0}\left(2\frac{e^{-i\omega_0  t}-1}{\omega^2_0}+i\frac{t}{\omega_0}(1 +  e^{-i\omega_0 t})\right)      \\
c_b^3 &=& -\frac{1}{\sqrt{2}}\frac{V^3}{\omega_0}\left(2\frac{e^{i\omega_0  t}-1}{\omega^2_0}-i\frac{t}{\omega_0}(1 +  e^{i\omega_0 t})\right)      \\
\end{eqnarray*}
where
\begin{eqnarray*}
\int t e^{i\omega t} dt&=& -i\frac{d}{d\omega}\int e^{i\omega t} dt\\
&=&  -i\frac{d}{d\omega}\frac{e^{i\omega t}}{i\omega} \\
&=&  -\frac{d}{d\omega}\frac{e^{i\omega t}}{\omega} \\
&=&  -\left(-\frac{e^{i\omega t}}{\omega^2_0}+i\frac{t e^{i\omega t}}{\omega}\right)\\
&=&  \left(\frac{e^{i\omega t}}{\omega^2_0}-i\frac{t e^{i\omega t}}{\omega}\right) 
\end{eqnarray*}
In general
\begin{eqnarray*}
\int t^n e^{i\omega t} dt&=& \frac{1}{i^n}\frac{d^n}{d\omega^n}\int e^{i\omega t} dt\\
\int t^2 e^{i\omega t} dt&=& -\frac{d^2}{d\omega^2}\int e^{i\omega t} dt\\
&=& -\frac{d^2}{d\omega^2}\left(\frac{ e^{i\omega t}}{i\omega}\right)\\
&=& -\frac{d}{d\omega}\left(\frac{ it e^{i\omega t}}{i\omega}-\frac{e^{i\omega t}}{i\omega^2}\right)\\
&=& -\left(\frac{ -t^2 e^{i\omega t}}{i\omega} -\frac{te^{i\omega t}}{\omega^2}-\left(\frac{ite^{i\omega t}}{i\omega^2} -2\frac{e^{i\omega t}}{i\omega^3} \right)     \right)\\
&=& \left(\frac{ -it^2 e^{i\omega t}}{\omega}+2\frac{te^{i\omega t}}{\omega^2} +2i\frac{e^{i\omega t}}{\omega^3}    \right)\\
\end{eqnarray*}

Substitue $c^3_a(t)$, $c^3_b(t)$ into \ref{eq:tdpde} and integrate again

\begin{eqnarray*}
i\dot c^4_a&=& c^3_b Ve^{-i\omega_0 t} = \left(     -\frac{1}{\sqrt{2}}\frac{V^3}{\omega_0}\left(2\frac{e^{i\omega_0  t}-1}{\omega^2_0}-i\frac{t}{\omega_0}(1 +  e^{i\omega_0 t})\right)\right)   Ve^{-i\omega_0 t}\\
c^4_a&=& -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\int_0^t\left(2\frac{e^{i\omega_0  \tp}-1}{\omega_0}-i\tp(1 +  e^{i\omega_0 \tp})\right)e^{-i\omega_0 \tp}d\tp\\
&=& -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\int_0^t\left(2\frac{1-e^{-i\omega_0  \tp}}{\omega_0}-i\tp(e^{-i\omega_0 \tp})-i\tp\right) d\tp\\
&=& -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\left(2\frac{\tp}{\omega_0}-2\frac{e^{-i\omega_0  \tp}}{-i\omega^2_0}-i\left(\frac{e^{-i\omega_0 \tp}}{\omega^2_0}+i\frac{\tp e^{-i\omega_0 \tp}}{\omega_0}\right) -i\frac{\tp^2}{2}\right)|_0^t\\
&=& -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\left(2\frac{t}{\omega_0}-2\frac{e^{-i\omega_0  t}-1}{-i\omega^2_0}-i\left(\frac{e^{-i\omega_0 t}-1}{\omega^2_0}+i\frac{t e^{-i\omega_0 t}}{\omega_0}\right) -i\frac{t^2}{2}\right)\\
&=& -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\left(\frac{t}{\omega_0}(2+e^{-i\omega_0 t})-3i\frac{e^{-i\omega_0  t}-1}{\omega^2_0} +\frac{t^2}{2}\right)\\
c^4_b&=&-\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\left(-\frac{t}{\omega_0}(2+e^{i\omega_0 t})-3i\frac{e^{i\omega_0  t}-1}{\omega^2_0} +\frac{t^2}{2}\right)\\
\end{eqnarray*}
Substitue $c^4_a(t)$, $c^4_b(t)$ into \ref{eq:tdpde} and integrate again

\begin{eqnarray*}
i\dot c^5_a&=& c^4_b Ve^{-i\omega_0 t} = -\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\left(-\frac{t}{\omega_0}(2+e^{i\omega_0 t})-3i\frac{e^{i\omega_0  t}-1}{\omega^2_0} +\frac{t^2}{2}\right)   Ve^{-i\omega_0 t}\\
 c^5_a&=& i\frac{1}{\sqrt{2}}\frac{V^4}{\omega^2_0}\int_0^t\left(-\frac{\tp}{\omega_0}(2+e^{i\omega_0 \tp})-3i\frac{e^{i\omega_0  \tp}-1}{\omega^2_0} +\frac{\tp^2}{2}\right)   Ve^{-i\omega_0 \tp} d\tp\\
&=& i\frac{1}{\sqrt{2}}\frac{V^5}{\omega^2_0}\int_0^t\left(-\frac{\tp}{\omega_0}(2e^{-i\omega_0 \tp}+1)-3i\frac{1-e^{-i\omega_0  \tp}}{\omega^2_0} +\frac{\tp^2e^{-i\omega_0 \tp}}{2}\right) d\tp\\
&=& i\frac{1}{\sqrt{2}}\frac{V^5}{\omega^2_0}\int\left(-2\frac{\tp}{\omega_0}e^{-i\omega_0 \tp}-\frac{\tp}{\omega_0}-3i\left(\frac{1}{\omega_0^2} -\frac{e^{-i\omega_0  \tp}}{\omega^2_0}\right) +\frac{\tp^2e^{-i\omega_0 \tp}}{2}\right) d\tp\\
&=& i\frac{1}{\sqrt{2}}\frac{V^5}{\omega^2_0}\left(-\frac{2}{\omega_0}\left(\frac{e^{-i\omega_0 \tp}}{\omega_0^2}+i\frac{\tp e^{-i\omega_0 \tp}}{\omega_0}\right)
-\frac{\tp^2}{2\omega_0}-3i\left(\frac{\tp}{\omega_0^2} -\frac{e^{-i\omega_0  \tp}}{-i\omega^3_0}\right) 
+\onehalf  
\left(\frac{ i\tp^2 e^{-i\omega_0 \tp}}{\omega_0}+2\frac{\tp e^{-i\omega_0 \tp}}{\omega_0^2} -2i\frac{e^{-i\omega_0 \tp}}{\omega_0^3}    \right)
\right)|_0^t\\
&=& i\frac{1}{\sqrt{2}}\frac{V^5}{\omega^2_0}\left(-\frac{5}{\omega_0}\left(\frac{e^{-i\omega_0 \tp}}{\omega_0^2}\right)+\left(-2i\frac{\tp e^{-i\omega_0 \tp}}{\omega_0^2}\right)
-\frac{\tp^2}{2\omega_0}-3i\left(\frac{\tp}{\omega_0^2}\right) 
+\frac{ i\tp^2 e^{-i\omega_0 \tp}}{2\omega_0}+\frac{\tp e^{-i\omega_0 \tp}}{\omega_0^2} -i\frac{e^{-i\omega_0 \tp}}{\omega_0^3} 
\right)|_0^t\\
&=& i\frac{1}{\sqrt{2}}\frac{V^5}{\omega^5_0}\left[\left(-5 -2i\omega_0\tp
+\frac{ i\tp^2\omega_0^2}{2}+\omega_0\tp -i\right)e^{-i\omega_0 \tp} 
-\frac{\omega_0^2\tp^2}{2}-3i\omega_0\tp\right] 
|_0^t\\
\end{eqnarray*}
\end{comment}
