\section{Muon in a time varying magnetic field}
\label{sec:NMR}
We will begin by reviewing the dynamics of a muon at rest in a time varying field, and then adapt that well known 
formalism to a muon in motion in a fixed field\cite{griffiths}[p. 330] ({\it see Griffiths Quantum Mechanics 1st edition, problem 9.7 on page 305
and section 10.1.3 page 330)}.
The time dependence of the spin of a muon at rest in a magnetic field is given by
$$\der{\bfs}{t}=\bm\mu\times\bfB$$ where $\bm\mu=\frac{g e}{2mc}\bfs$. The system can equivalently be described in terms of an interaction energy ${\cal H}=-\bm\mu\cdot\bfB$.
The Hamiltonian is
$${\cal H} = -{\bm\mu}\cdot\bfB = -\frac{ge}{2mc}\frac{\hbar}{2}{\bm\sigma}\cdot \bfB =   -\frac{ge}{2mc}\frac{\hbar}{2} \begin{pmatrix} B_z & B_x - iB_y \cr B_x + i B_y & -B_z\end{pmatrix}$$   


%\subsection{At rest in field with time varying direction}
Suppose that the z-component of the magnetic field is fixed and that the magnetic field transverse to the z-axis rotates about that axis with frequency $\omega$ (see Figure~\ref{fig:rotating_field}), so that
$$\bfB = B_z\bfhz + B_{xy}(\cos(\omega t)\bfhx - \sin(\omega t)\bfhy)$$
\begin{figure}
\centering
\includegraphics{figures/rotating_field.pdf}
\caption{The z-component of the magnetic field is fixed. The transverse component of the field rotates about the z-axis with frequency $\omega$}
\label{fig:rotating_field}
\end{figure}

Then
\begin{equation}
{\cal H} =   -\frac{ge}{2mc}\frac{\hbar}{2} \begin{pmatrix} B_z & B_{xy}e^{i\omega t} \cr B_{xy}e^{-i\omega t} & -B_z\end{pmatrix}\equiv -\frac{\hbar}{2} \begin{pmatrix} \omega_0 & \Omega e^{i\omega t} \cr \Omega e^{-i\omega t} & -\omega_0\end{pmatrix}\label{eq:hamiltonian}
\end{equation}
where $\omega_0=\frac{ge}{2mc}B_z$ and $\Omega=\frac{ge}{2mc}B_{xy}$.
 The time evolution is governed by Schrodinger's equation. 
$$i\hbar\partder{}{t}\chi = {\cal H}\chi$$
For a spinor $\chi(t) = (a(t), b(t))$, we can write:
\begin{equation}
i\hbar \partder{}{t}\left(\begin{array}{c} a \cr b\end{array}\right)=-{\hbar \over 2}\left( \begin{array}{cc} 
    \omega_0& \Omega e^{i \omega t} \\
    \Omega e^{-i \omega t} & -\omega_0\end{array} \right)\left(\begin{array}{c} a \cr b\end{array}\right)
     \label{eq:schrodinger}
\end{equation}
\begin{equation}
 \left( \begin{array}{c} \dot{a}(t) \\ \dot{b}(t) \end{array} \right)
=  {i \over 2 }\left( \begin{array}{c} 
    a(t)\omega_0 + b(t)\Omega e^{i \omega t} \\
    a(t) \Omega e^{-i \omega t}  -b(t)\omega_0 \end{array} \right).\label{eq:coupleddiffequations}
\end{equation}
\input{rotatingframe}
(Alternatively see Appendix ~\ref{ss:labframesolution} for direct solution of the coupled differential equations~\ref{eq:coupleddiffequations})
\section{Observables and Polarization}
By whichever method we choose, we find the solution to ~\ref{eq:schrodinger}
\begin{equation}
\psi(t) = \begin{pmatrix} a \cr b\end{pmatrix} 
= \begin{pmatrix} \left(\cos({\omega^\prime t/2}) -{i(\omega-\omega_0)\over \omega^\prime}\sin({\omega^\prime t /2})\right)e^{i{\omega\over 2}t} & -\frac{i \Omega}{\omega^\prime}\sin(\omega^\prime t /2)e^{i{\omega\over 2}t} \cr
   -\frac{i \Omega}{\omega^\prime}\sin(\omega^\prime t /2)e^{-i{\omega\over 2}t}& \left(\cos({\omega^\prime t/2}) +{i(\omega-\omega_0)\over \omega^\prime}\sin({\omega^\prime t /2})\right)e^{-i{\omega\over 2}t}
 \end{pmatrix}
\begin{pmatrix} a_0 \cr b_0\end{pmatrix}\label{eq:umatrix}
\end{equation}
where 
\begin{equation}
\omega^\prime = [(\omega_0-\omega)^2 + \Omega^2]^{1/2}\label{eq:omegap}
\end{equation}
and $a_0$ and $b_0$ define the initial spin state.
The polarization (the observable),
\begin{equation}
\bfs = \onehalf\expect{\psi}{\bm\sigma}{\psi} = \onehalf\left((a^*b+b^*a)\bfhx + -i(a^*b-b^*a)\bfhy+(|a(t)|^2-|b(t)|^2)\bfhz\right)\label{eq:polarization}
\end{equation}
Since
\def\opt2{{\omega^\prime t/2}}
\def\op{{\omega^\prime}}
\def\opp{{\omega^{\prime\prime}}}
\begin{eqnarray*}
 a^* &=&
\left[\left(\cos(\opt2)+\frac{i(\omega-\omega_0)}{\op}\sin(\opt2)\right)a_0^* +\frac{i\Omega}{\op}\sin(\opt2)b_0^*\right]e^{-i\frac{\omega t}{2}} \\
b&=&\left[-\frac{i\Omega}{\op}\sin(\opt2)a_0 + \left(\cos(\opt2)+\frac{i(\omega-\omega_0)}{\op}\sin(\opt2)\right)b_0\right]e^{-i\frac{\omega t}{2}}
\end{eqnarray*}
we can write
\begin{eqnarray*}
a^* b&=&
\left(a_0^*b_0\cos^2(\omega^\prime t/2)-\frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2a_0^* b_0-\Omega^2b_0^* a_0-(\omega-\omega_0)\Omega(|a_0|^2-|b_0|^2)]\sin^2(\omega^\prime t/2)\right.\\
&+& \left.\frac{i}{\omega^\prime}\cos(\omega^\prime t/2)\sin(\omega^\prime t/2)\left(a_0^*((\omega-\omega_0)b_0-\Omega a_0)-b_0((\omega_0-\omega)a_0^*-\Omega b_0^*) \right)\right] e^{-i\omega t}
\end{eqnarray*}
If the intial state is $a_0=b_0=\frac{1}{\sqrt{2}}$ (an eigenket of $\sigma_x$), the intial state polarization is $s_x=1$.
Then
\begin{eqnarray*}
a^*b&=& \onehalf\left(\cos^2(\omega^\prime t/2)-\frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\omega^\prime t/2)\right.\\
&+& \left.\frac{i}{\omega^\prime}\cos(\omega^\prime t/2)\sin(\omega^\prime t/2)\left(((\omega-\omega_0)-\Omega )-((\omega_0-\omega)-\Omega ) \right)\right) e^{-i\omega t}\\
&=& \onehalf\left(\cos^2(\omega^\prime t/2)-\frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\omega^\prime t/2)+\frac{2i}{\omega^\prime}\cos(\omega^\prime t/2)\sin(\omega^\prime t/2)(\omega-\omega_0)\right) e^{-i\omega t}\\
\end{eqnarray*}
The polarization
\begin{eqnarray}
\vev{s_x} &=& \left(\cos^2(\omega^\prime t/2)-\frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\omega^\prime t/2)\right)\cos(\omega t)
+\left(\frac{2}{\omega^\prime}\cos(\omega^\prime t/2)\sin(\omega^\prime t/2)(\omega-\omega_0)\right) \sin(\omega t)\nonumber\\
\vev{s_y} &=& -\left(\cos^2(\omega^\prime t/2)-\frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\omega^\prime t/2)\right)\sin(\omega t)
+\left(\frac{2}{\omega^\prime}\cos(\omega^\prime t/2)\sin(\omega^\prime t/2)(\omega-\omega_0)\right) \cos(\omega t)\nonumber\\
\vev{s_z}&=& \onehalf(|a|^2-|b|^2) = \onehalf\left[\cos^2(\opt2)+\frac{(\omega-\omega_0+\Omega)^2}{\op^2}(\sin^2\opt2) -\left(\cos^2(\opt2)+\frac{(\omega-\omega_0-\Omega)^2}{\op^2}\sin^2(\opt2)\right)\right]\nonumber\\
&=& \frac{(\omega-\omega_0)\Omega}{\op^2}\sin^2(\opt2) \label{eq:spinvector}
\end{eqnarray}

\section{Effective Hamiltonian for Thomas equation}
\label{sec:bmt}
Of course we are not interested in the precession of muon at rest in a time varying magnetic field, but rather a muon in motion in a fixed field, as described by the Thomas-BMT\cite{jackson}[p. 559] equation.
\begin{equation}
\der{\bfs}{t}=\frac{e}{mc}\bfs\times\left[(a_\mu+\frac{1}{\gamma})\bfB -a_\mu\frac{\gamma}{\gamma+1}({\bm \beta}\cdot\bfB){\bm \beta}-\left(\frac{g}{2}-\frac{\gamma}{\gamma+1}\right){\bm\beta}\times\bfE\right]\label{eq:j170}
\end{equation}
or equivalently
$${\cal H} = -\tilde{\bm\mu}\cdot{\bf\cal B}$$ where $\tilde{\bm\mu} = \frac{e}{mc}\bfs$ and ${\bf\cal B}$ is an effective magnetic field,
\begin{equation}
{\bf\cal B} =  \left[(a_\mu+\frac{1}{\gamma})\bfB -a_\mu\frac{\gamma}{\gamma+1}({\bm \beta}\cdot\bfB){\bm \beta}-\left(\frac{g}{2}-\frac{\gamma}{\gamma+1}\right){\bm\beta}\times\bfE\right]\nonumber
\end{equation}



\begin{figure}
\centering
\begin{minipage}{.45\textwidth}
  \centering
  \includegraphics[width=1\linewidth]{figures/diagram.pdf}
  \captionof{figure}{The positive muon is injected at the top of the ring and rotates in the clockwise direction. The magnetic field is out of the page in the
positive z-direction.}
  \label{fig:diagram}
\end{minipage}%
\begin{minipage}{.45\textwidth}
  \centering
  \includegraphics[width=1\linewidth, trim =0 0 0 0, clip=.true.]{figures/spiral.pdf}
  \captionof{figure}{Muon trajectory in a uniform magnetic field.}
  \label{fig:spiral}
\end{minipage}
\end{figure}


\section{Hamiltonian for a Muon with a Vertical Component of Velocity} 
\label{sec:pitchanalysis}
\subsection{Fixed pitch angle (spiral)}
\label{ss:fixedpitchbetaz}
Consider a muon in a uniform magnetic field with a component of velocity parallel to the field. 
We work in a fixed cartesian coordinate system with $\bfB= B_z\bfhz$.
(Our coordinate system is shown in Figure~\ref{fig:diagram}.) According to the Lorentz force law\cite{jackson}[p. 559]
$$\der{\bm\beta}{t} = \frac{e}{\gamma m}[\bfE+\bm\beta\times\bfB-\bm\beta(\bm\beta\cdot\bfE)]$$
the particle trajectory is a clockwise spiral (see Figure~\ref{fig:spiral}), with velocity
$$\bm\beta = \beta_{xy}(\cos(\omega t)\bfhx - \sin(\omega t)\bfhy) +\beta_z\bfhz.$$ (For the time being we assume that $\bfE=0$.)

Then the effective magnetic field is
\begin{eqnarray*}
{\bf\cal B} &=&  \left[(a_\mu+\frac{1}{\gamma})B_z\bfhz -a_\mu\frac{\gamma}{\gamma+1}( \beta_z B_z)\left(\beta_z\bfhz+\beta_{xy}(\cos\omega t\bfhx -\sin\omega t\bfhy)\right)\right]\\
 &=&  B_z\left[\left((a_\mu+\frac{1}{\gamma}) -a_\mu\frac{\gamma}{\gamma+1} \beta_z^2\right)\bfhz
 -\left(a_\mu\frac{\gamma}{\gamma+1} \beta_z\beta_{xy}(\cos(\omega t)\bfhx -\sin(\omega t)\bfhy\right)\right]\\
\end{eqnarray*}
Define
\begin{eqnarray}
\omega_0 &\equiv&\frac{e}{mc}B_z\left[a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\beta_z^2\right]\label{eq:omegaodef}\\
\Omega&\equiv&-\frac{e}{mc}B_z\left[a_\mu\frac{\gamma}{\gamma+1}\beta_z\beta_{xy}\right]\label{eq:bomegadef}\\
\omega&=& \omega_c = \frac{eB_z}{mc\gamma}\label{eq:omegadef}
\end{eqnarray}
and our Hamiltonian assumes the form of Equation~\ref{eq:hamiltonian},
\begin{equation}
{\cal H} = -\tilde{\bm\mu}\cdot{\bf \cal B} = -\frac{e}{mc}\frac{\hbar}{2}\bm\sigma\cdot {\bf\cal B} = -\frac{\hbar}{2}\begin{pmatrix}\omega_0 & \Omega e^{i\omega t} \cr \Omega e^{-i\omega t} & -\omega_0\end{pmatrix}.\label{eq:pitchhamiltonian}
\end{equation}
%as does the solution.
In the rotating frame 
\begin{equation}
{\cal H}_{rot} = -\frac{\hbar}{2}\begin{pmatrix}\omega_0-\omega & \Omega  \cr \Omega & -\omega_0+\omega\end{pmatrix}\nonumber
\end{equation}
and the eigenvalue is the precession frequency in the rotating frame.
Now we simply substitute ~\ref{eq:omegaodef} and ~\ref{eq:bomegadef} into Equations~\ref{eq:spinvector}, to get the time dependence of the polarization.
\begin{figure}
\centering
\includegraphics{figures/polarization_analytic_0.pdf}
\includegraphics{figures/polarization_simulation_0.pdf}
\caption{Polarization of a muon following a circular trajectory, $\beta_z=0$. The analytic solution (Equations ~\ref{eq:spinvector}) is on the top, and simulation (integration of the equations of motion and
the BMT equation) at the bottom.}
\label{fig:polarizationnopitch}
\end{figure}
\begin{figure}
\centering
\includegraphics{figures/polarization_analytic_04.pdf}
\includegraphics{figures/polarization_simulation_04.pdf}
\caption{Polarization of a muon following a spiraling trajectory, $\beta_z=0.4$. The analytic solution (Equations ~\ref{eq:spinvector}) is on the top, and simulation (integration of the equations of motion and
the BMT equation) at bottom.}
\label{fig:polarizationwithpitch}
\end{figure}


The projection of the spin onto the direction of motion is given by 
\begin{eqnarray}
|\bm\beta|\bmbetahat\cdot s &=& \left(\cos^2(\opt2) -\frac{1}{\op^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\opt2)\right)\beta_{xy}(\cos^2(\omega t)+\sin^2(\omega t)) + \beta_z\frac{2(\omega-\omega_0)\Omega}{\op^2}\sin^2(\opt2)\nonumber\\
&=& \left(\cos^2(\opt2) -\frac{1}{\op^2}[(\omega-\omega_0)^2-\Omega^2]\sin^2(\opt2)\right)\beta_{xy} + \beta_z\frac{2(\omega-\omega_0)\Omega}{\op^2}\sin^2(\opt2)\nonumber\\
&=& \onehalf\left(1+\cos(\op t) -\frac{1}{\op^2}[(\omega-\omega_0)^2-\Omega^2](1-\cos\op t)\right)\beta_{xy} + \beta_z\frac{(\omega-\omega_0)\Omega}{\op^2}(1-\cos(\op t))\label{eq:sdotbpitch}
\end{eqnarray}
The thing we measure is the time rate of change of the projection of the polarization on the direction of motion
\begin{eqnarray*}
\der{}{t}(\bm\beta\cdot\bfs)&=& -\frac{\op}{2}\left[\left(1 +\frac{1}{\op^2}[(\omega-\omega_0)^2-\Omega^2]\right)\beta_{xy} - \beta_z\frac{(\omega-\omega_0)\Omega}{\op^2}\right]\sin(\op t)
\end{eqnarray*}
%The projection of the polarization on the direction of motion oscillates with frequency $\omega^\prim$. In the absence of a vertical
%component of velocity ($\beta_z=0$ and $\Omega=0$), $\omega^\prime = \omega_0-\omega$.
%For finite $\Omega$ (and recall that $\Omega$ is proportional to $\beta_z$), 
In general $\bmbetahat\cdot \bfs$ oscillates with frequency $\omega^\prime$, the eigenvalue of the Hamiltonian in the rotating frame. 

Finally, we use ~\ref{eq:omegaodef},\ref{eq:bomegadef},\ref{eq:omegadef} to evaluate $\omega^\prime$
\begin{eqnarray}
%\omega_0 &=&\frac{eB}{mc}\left(a_\mu+\frac{1}{\gamma} -a_\mu\frac{\gamma}{\gamma+1}\beta_z^2\right) = \frac{eB}{mc}\left(\frac{1}{\gamma} +a_\mu(\frac{1+\gamma(1-\beta_z^2)}{\gamma+1})\right)
%\label{eq:omegazero}\\
%\omega &=& \omega_c = \frac{eB}{mc\gamma}\nonumber\\
%\Omega &=& -\frac{e}{mc}\left(Ba_\mu\frac{\gamma}{\gamma+1}\beta_z(\beta^2-\beta_z^2)^{1/2}-\left(\frac{g}{2}-\frac{\gamma}{\gamma+1}\right)\bm\beta\times\bfE\right)\nonumber\\
%&\rightarrow& -\frac{e}{mc}\left(Ba_\mu\frac{\gamma}{\gamma+1}\beta_z-\left(a_\mu+\frac{1}{\gamma+1}\right)iE_z\right)(\beta^2-\beta_z^2)^{1/2}\bmbetahat\label{eq:bomega}\\
{\omega^\prime}^2 &=&(\omega_0-\omega)^2+\Omega^2= \left(\frac{eB}{mc}\right)^2\left[\left(\frac{1}{\gamma} +a_\mu(\frac{1+\gamma(1-\beta_z^2)}{\gamma+1}) -\frac{1}{\gamma}\right)^2 +
\left( a_\mu\frac{\gamma}{\gamma+1}\beta_z(\beta^2-\beta_z^2)^{1/2}\right)^2\right]\label{eq:exactfreq}\\
&=& \omega_a^2\left(\frac{1+\gamma(1-\beta_z^2))^2+ \gamma^2\beta_z^2(\beta^2-\beta_z^2)}{(\gamma+1)^2}\right)=\omega_a^2\left(\frac{(1+2\gamma(1-\beta_z^2)+\gamma^2(1-\beta_z^2)^2+ \gamma^2\beta_z^2(\beta^2-\beta_z^2)}{(\gamma+1)^2}\right)\nonumber\\
&=&\omega_a^2\left(\frac{1+2\gamma(1-\beta_z^2)+\gamma^2(1-2\beta_z^2)+ \gamma^2\beta_z^2\beta^2}{(\gamma+1)^2}\right)\nonumber\\
&=&\omega_a^2\left(\frac{(1+2\gamma(1-\beta_z^2)+\gamma^2(1-2\beta_z^2))+ \beta_z^2(\gamma^2-1)}{(\gamma+1)^2}\right)\nonumber\\
&=&\omega_a^2\left(\frac{(1+2\gamma(1-\beta_z^2)+\gamma^2(1-\beta_z^2))- \beta_z^2}{(\gamma+1)^2}\right)\nonumber\\
\omega^\prime&=& \omega_a\left(\frac{(1+\gamma)(1-\beta_z^2)^{1/2}}{\gamma+1}\right)\nonumber\\
\omega^\prime&=& \omega_a(1-\beta_z^2)^{1/2}\label{eq:fixedpitch}
\end{eqnarray}
$\op$ is the precession frequency measured in the rotating frame, the familiar pitch correction.

It turns out that the spiral trajectory corresponding to a fixed pitch can be transformed by a Lorentz boost ($\beta_z$)
into a circular trajectory with zero pitch. Frequencies measured in the two frames are related by the time dilation factor $\gamma=\sqrt{1-\beta_z^2}$.
The result is of course equivalent to~\ref{eq:fixedpitch}. 
See ~\ref{ss:lorentzboost} for details. 

%For $\Omega \ll \omega_0$, $$\omega^\prime&\sim&\frac{eB}{mc}a_\mu(1-\frac{1}{2}\beta_z^2)$$
\begin{eqnarray*}
%\omega^\prime &=& [(\omega_0-\omega)^2+\Omega^2]^{1/2}\sim (\omega_0-\omega)(1+\onehalf\frac{\Omega^2}{(\omega-\omega_0)^2}\\
%&\sim& (\omega_0-\omega)\left(1+ \onehalf\frac{\gamma^2\beta^2_z(\beta^2-\beta_z)}{(1+\gamma(1-\beta_z^2))^2}\right) =\lim_{\beta_z\ll\beta} \sim (\omega_0-\omega)(1+\frac{\gamma^2\beta_z^2\beta^2}{2(1+\gamma)^2})\\
%&\sim&(\omega_0-\omega)(1+\frac{(\gamma-1)\beta_z^2}{2(1+\gamma)})\\
%&\sim&\frac{eB}{mc}a_\mu\left(\frac{1+\gamma(1-\beta_z^2)}{\gamma+1}\right)(1+\frac{(\gamma-1)\beta_z^2}{2(1+\gamma)})\\
%&\sim&\frac{eB}{mc}a_\mu\left(\frac{1+\gamma(1-\beta_z^2)}{\gamma+1} + \frac{(\gamma-1)\beta_z^2+\gamma(\gamma-1)\beta_z^2(1-\beta_z^2)}{2(1+\gamma)^2}\right)\\
%&\sim &\frac{eB}{mc}a_\mu\left(\frac{1+\gamma(1-\beta_z^2)}{\gamma+1} + \frac{(\gamma-1)\beta_z^2}{2(1+\gamma)}\right)\\
%&\sim&\frac{eB}{mc}a_\mu\left(\frac{1+\gamma-\gamma\beta_z^2+\onehalf(\gamma-1)\beta_z^2}{\gamma+1}\right)\\
%&\sim&\frac{eB}{mc}a_\mu\left(\frac{1+\gamma-\frac{1}{2}\gamma\beta_z^2-\onehalf\beta_z^2}{\gamma+1}\right)=\frac{eB}{mc}a_\mu\left(\frac{(1+\gamma)(1-\frac{1}{2}\beta_z^2)}{\gamma+1}\right)\\
\end{eqnarray*}
For $\beta_z\ll \beta$, and $\beta\sim 1$, $\beta_z=\psi$ and $\op=\omega_a(1-\onehalf\psi^2)$ where $\psi$ is the pitch angle, in good agreement with 
Farley\cite{farley}, Silenko\cite{silenko}, Kim\cite{kim}, Miller\cite{miller_private} and others.

%sure enough consistent with the argument based on Lorentz dilation, Equation~\ref{eq:lorentz} in section~\ref{sec:lorentz}.

\begin{figure}
\centering
\includegraphics{figures/betadots_analytic.pdf}
\includegraphics{figures/betadots_simulation.pdf}
\caption{Projection of the polarization onto the direction of motion for a muon following circular ($\beta_z=0$ and  spiraling $\beta_z=0.4$ trajctories. 
The analytic solution (Equation ~\ref{eq:sdotbpitch}) is on the top, and simulation (integration of the equations of motion and
the BMT equation) at the bottom.}
\label{fig:betadots}
\end{figure}
%\section{QM with Efield}
\subsection{Equations of motion in a uniform magnetic field and vertically focusing electric field.}
\label{ss:lorentzforce}
The next step is to include the electric field and vertical oscillations. In order to construct the effective magnetic field in our Hamiltonian we need
to solve the equations of motion to determine the velocity. For motion in electromagnetic fields
\begin{equation}
\frac{d\bm\beta}{dt}=\frac{e}{\gamma m c}[\bfE + \bm\beta\times\bfB - \bm\beta(\bm\beta\cdot\bfE)]
\label{eq:eom}
\end{equation}
The electric and magnetic fields are both vertical (z-direction). With that stipulation~\ref{eq:eom} reduces to the coupled differential equations
\begin{eqnarray}
\frac{d\beta_z}{dt}&=&\frac{e}{\gamma m c}E_z(1 - \beta_z^2)\label{eq:eomez}\\
\frac{d\beta_x}{dt}&=&\frac{e}{\gamma m c}(\beta_y B_z - \beta_x\beta_z E_z)\label{eq:eomex}\\
\frac{d\beta_y}{dt}&=&-\frac{e}{\gamma m c}(\beta_x B_z + \beta_y\beta_z E_z)\label{eq:eomey}
\end{eqnarray}
The tranverse momentum (momentum in the x-y plane) is conserved so that $$\frac{p_{xy}}{m}=\gamma(t)\bm\beta_{xy}(t)=\gamma^0\beta^0_{xy}={\rm constant}$$
(Forces in the transverse plane are always perpendicular to the velocity.)
The total momentum depends on the z-component of velocity.
\begin{eqnarray}
\beta_{xy}&=&\beta_{xy}^0\frac{\gamma_0}{\gamma} = \beta_{xy}^0\left[\frac{1-(\beta_{xy}^2+\beta_z^2)}{1-{\beta^0_{xy}}^2}\right]^{1/2}\nonumber\\
\rightarrow \beta_{xy}^2(1-{\beta_{xy}^0}^2)&=&{\beta_{xy}^0}^2(1-(\beta_{xy}^2+\beta_z^2))\nonumber\\
\beta_{xy}^2&=&{\beta_{xy}^0}^2(1-\beta_z^2)\nonumber\\
\rightarrow \frac{\beta_{xy}}{\beta_{xy}^0}&=&(1-\beta_z^2)^{1/2} = \frac{\gamma_0}{\gamma}\label{eq:gamma0gammaz}
\end{eqnarray}
and in an obvious notation, $\gamma = \gamma_z\gamma_0$. With that handy relationship 
we rewrite the first of the three coupled equations as 
\begin{eqnarray}
\frac{d\beta_z}{dt}&=&\frac{e}{m c}E_z\frac{(1 - \beta_z^2)}{\gamma_0\gamma_z}\nonumber\\
&=&\frac{e}{m c}E_z\frac{(1 - \beta_z^2)^{3/2}}{\gamma_0}\label{eq:eomezred}
\end{eqnarray}
In the limit where $\beta_z\ll 1$, and assuming $E_z = -kz$ (vertically focusing linear restoring force and continuous quads)
\begin{eqnarray*}
\frac{d\beta_z}{dt}&\sim& \frac{e}{\gamma_0m c}E_z\\
\rightarrow\frac{1}{c}\der{^2z}{t^2}&\sim& -\frac{e}{\gamma_0 m c}kz\\
\rightarrow z&\sim& z_0\cos\omega_v t, \ \beta_z\sim-\frac{z_0}{c}\omega_v\sin\omega_v t\sim-\beta_z^0\sin\omega_v t
\end{eqnarray*}
where $$\omega_v = \sqrt{\frac{e k}{\gamma_0 m}}$$
We compute the next higher order term by substitution of the zeroth order solution into ~\ref{eq:eomezred} and integrating.
To order $\beta_z^2$,
\begin{eqnarray}
\beta_z &=& -\frac{ekz_0}{\gamma_0 m c}\int \cos\omega_v t\left(1-\left(\frac{z_0}{c}\omega_v\right)^2\sin^2\omega_v t\right)^{3/2} dt\nonumber\\
\beta_z &\sim& -\frac{ekz_0}{\gamma_0 m c}\int \cos\omega_v t\left(1-\frac{3}{2}\left(\frac{z_0}{c}\omega_v\right)^2\sin^2\omega_v t\right) dt\nonumber\\
 &\sim& -\frac{ekz_0}{\gamma_0 m c}\frac{1}{\omega_v} \left(\sin\omega_v t-\frac{1}{2}\left(\frac{z_0}{c}\omega_v\right)^2\sin^3\omega_v t\right) \nonumber\\
 &\sim& -\omega_v\frac{z_0}{ c}\sin\omega_v t[\cos\left(\left(\frac{z_0}{c}\omega_v\right)\sin\omega_v t\right)] \nonumber\\
 &\sim& -\beta_z^0\sin\omega_v t[\cos\left(\beta_z^0\sin\omega_v t\right)] \nonumber\\
 &\sim& -[\sin(\beta_z^0\sin\omega_v t)][\cos\left(\beta_z^0\sin\omega_v t\right)] \nonumber\\
 &\sim& -\onehalf\sin\left(2\beta_z^0\sin\omega_v t\right)] \sim -\sin(\beta_z^0\sin\omega_v t)\label{eq:betaz2o}
\end{eqnarray}
Using ~\ref{eq:gamma0gammaz} we write
\begin{eqnarray}
\beta_{xy}&=&(1-\beta_z^2)^{1/2}\beta^0_{xy} \sim\beta_{xy}^0\cos(\beta_z^0\sin\omega_v t)\label{eq:betax2o}
\end{eqnarray}
and to second order in the small parameter $\beta_z^0$,
\begin{equation}
\bm\beta = 
\beta^0_{xy}\cos(\beta_z^0\sin\omega_v t)(\cos(\omega t)\bfhx -\sin(\omega t)\bfhy) -\sin(\beta_z^0\sin\omega_v t)\bfhz
\end{equation}
which is confirmed by substitution into Equations~\ref{eq:eomez},\ref{eq:eomex} and~\ref{eq:eomey}. See section~\ref{ss:checkeom} for details.
\subsection{Hamiltonian for uniform magnetic field and vertically focusing electric field}
\label{ss:uniformbverticale}
The electric field appears in the BMT equation in the cross product with velocity.
\begin{eqnarray*}
\bm\beta\times \bf E&=&  \beta_{xy}(t)(-\sin(\omega t)\bfhx - \cos(\omega t))\bfhy)E_z(t)\\
&=& \beta^0_{xy}\gamma_z^{-1}(t) (\sin(\omega t)\bfhx + \cos(\omega t))\bfhy) k z_0(\cos\omega_v t)
\end{eqnarray*}
%Let's write $\Omega = \Omega_0 \sin\omega_v t$.
Including the electric field our effective magnetic field 
${\bf \cal B}$ becomes
\begin{eqnarray*}
{\bf\cal B} &=&  \left[(a_\mu+\frac{1}{\gamma})B_z\bfhz -a_\mu\frac{\gamma}{\gamma+1}( \beta_z B_z)\left(\beta_z\bfhz+\beta_{xy}(t)(\cos\omega t\bfhx -\sin\omega t\bfhy)\right)\right.\\
 & &\left. -\left(a_\mu +\frac{1}{\gamma+1}\right)\beta_{xy}(t) E_z(t) (-\sin(\omega t)\bfhx - \cos(\omega t))\bfhy)
\right]\\
 &=&  B_z\left[\left((a_\mu+\frac{1}{\gamma}) -a_\mu\frac{\gamma}{\gamma+1} \beta_z^2\right)\bfhz
 -\left(a_\mu\frac{\gamma}{\gamma+1} \beta_z\beta_{xy}(\cos(\omega t)\bfhx -\sin(\omega t)\bfhy\right)\right]\\
 & & -\left(a_\mu +\frac{1}{\gamma+1}\right)\beta_{xy}(t) E_z(t) (-\sin(\omega t)\bfhx - \cos(\omega t))\bfhy)
\end{eqnarray*}
Then 
\begin{eqnarray}
{\cal H}&=&-\tilde{\bm\mu}\cdot {\bf \cal B}\nonumber\\
&=& -\frac{\hbar}{2}\begin{pmatrix} \omega_0(t) & (\Omega(t) + i\Pi(t))e^{i\omega t}\cr (\Omega(t) -i\Pi(t))e^{-i\omega t} & -\omega_0(t)\end{pmatrix}\label{eq:ham}
\end{eqnarray} 
where
\begin{eqnarray}
\Pi(t)&=& -\frac{e}{mc}\left(a_\mu+\frac{1}{\gamma+1}\right)E_z(t)\beta_{xy}(t)\label{eq:pidef}
\end{eqnarray}
and generalizing  ~\ref{eq:omegaodef} and~\ref{eq:bomegadef}
\begin{eqnarray}
\omega_0(t) &\equiv&\frac{e}{mc}B_z\left[a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\beta_z(t)^2\right]\label{eq:bomegatimedef}\\
\Omega(t) &=&-\frac{e}{mc}\left( B_za_\mu\frac{\gamma}{\gamma+1}\beta_z(t)\beta_{xy}(t)\right)\label{eq:omegaotimedef}
\end{eqnarray}
%\section{Magnetic and Electric fields}
%Finally we restore the magnetic field. 
Substituting the explicit time dependence of $\beta_z(t)$ and $E_z(t)$ into ~\ref{eq:bomegatimedef} and~\ref{eq:pidef}
\begin{eqnarray}
\Omega(t)+i\Pi(t)&=&-z_0\frac{eB}{mc}\left(a_\mu\frac{\gamma}{\gamma+1}\frac{-\omega_v\sin\omega_v t}{c} + i\left(a_\mu+\frac{1}{\gamma+1}\right)\frac{\omega_v^2}{c\omega}\cos\omega_v t\right)\beta_{xy}(t)\label{eq:omegaandpi}
%\Omega(t)+i\Pi(t)&=&-\frac{z_0}{c}\frac{eB}{mc}\left(a_\mu\frac{\gamma}{\gamma+1}{\omega_v\sin\omega_v t} + i\left(a_\mu+\frac{1}{\gamma+1}\right)\frac{\omega_v^2}{\omega}\cos\omega_v t\right)
%\beta_{xy}^{0}(1-z_0^2\omega_v^2/c^2\sin^2\omega_v t)^{1/2}
\end{eqnarray}
where since $\omega_v^2 = \frac{ek}{\gamma m}$ and $\omega = \frac{eB}{\gamma mc}$, we can write $k = B\frac{\omega_v^2}{c\omega}$.
Also from ~\ref{eq:omegaotimedef}
\begin{eqnarray}
\omega_0(t) &=& \frac{eB}{mc}\left(a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\beta_z^2\right)\nonumber\\
&=& \frac{eB}{mc}\left(a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\left(\frac{z_0\omega_v\sin(\omega_v t)}{c}\right)^2\right)\nonumber\\
&=& \frac{eB}{mc}\left(a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\left(\beta_z^0\sin(\omega_v t)\right)^2\right)\label{eq:omegatime}
%&=& A+F\sin^2(\omega_v t)
\end{eqnarray}
%The cyclotron frequency will also assume a time dependence. The equation of motion,
%neglecting the last term in ~\ref{eq:eom}, noting that $\bfE$ is in the z-direction, 
%and with the approximation that $\bm\beta\perp\bfE$  we find that
%$$\omega\rightarrow \frac{e}{\gamma mc}\sqrt{B^2 + E^2/\beta^2}$$
%Then
%\begin{eqnarray}
%{\cal H^\prime}&=&\begin{pmatrix}A+F\sin^2\omega_v t -\omega & C\sin\omega_v t + iD\cos\omega_v t\cr
%                       C\sin\omega_v t -i D\cos\omega_v t & -A-F\sin^2 \omega_v t + \omega \end{pmatrix}\label{eq:hp}
%\end{eqnarray}

\begin{eqnarray*}
{\cal H} &=& -\bm\mu\cdot{\bf{\cal B}} = -\frac{e}{mc}\frac{\hbar}{2}\bm \sigma\cdot{\bf{\cal B}} = -\frac{\hbar}{2}
\begin{pmatrix}\omega_0(t) & (\Omega(t) + i\Pi(t))e^{i\omega t}\cr (\Omega(t) -i\Pi(t))e^{-i\omega t} & -\omega_0\end{pmatrix}\\
%&=& -\bm\mu\cdot{\bf{\cal B}} = -\frac{ge}{2mc}\frac{\hbar}{2}\bm \sigma\cdot{\bf{\cal B}} = -\frac{\hbar}{2}
%\begin{pmatrix}\omega_0(t) & (\Omega(t) + i\Pi(t))\cr (\Omega(t) -i\Pi(t)) & -\omega_0(t)\end{pmatrix}\\
\end{eqnarray*}
Having constructed the Hamiltonian, all that remains is to solve
 Schrodinger's equation
\begin{eqnarray}
i\hbar\partder{}{t}\psi &=& {\cal H}\psi\\
i\partder{}{t}\begin{pmatrix} a \cr b\end{pmatrix} &=&
  -\onehalf
\begin{pmatrix}\omega_0(t) & (\Omega(t) + i\Pi(t))e^{i\omega t}\cr (\Omega(t) -i\Pi(t))e^{-i\omega t} & -\omega_0(t)\end{pmatrix}
\begin{pmatrix} a\cr b\end{pmatrix}\label{eq:pitchingse}
\end{eqnarray}
\subsection{Time Dependent Perturbation Theory}
We solve Schrodinger's equation perturbatively to second order in the small parameter $\beta_z^0 = -z_0\omega_v/c$.
%\subsection{Computing a frequency shift with time dependent perturbation theory}
The reader is likely most familiar with using time dependent perturbation theory to compute transition rates between eigenstates of the unperturbed Hamiltonian.
Our purpose is to determine a frequency shift in the observable $\bm\beta\cdot\bfs$. To see how we will extract a frequency shift from the Dyson series, we will study a problem for which we know the exact solution,
namely, the Hamiltonian of Equation ~\ref{eq:pitchhamiltonian}. If we treat the off diagnonal element as a perturbation,
which is valid if the pitch angle is small, then we can compute the frequency shift due to the perturbation and compare to the exact result. And we do just that in Appendix~\ref{ss:tdp}.
\input{revisited}

%\input{TDPT}

\begin{comment}
We solve the Schrodinger equations with Hamiltonian~\ref{eq:hamiltonian} and $\Omega,\Pi,$ and $\omega_0$ from~\ref{eq:omegaandpi} and~\ref{eq:omegatime} with the help of perturbation theory.
We saw above that in the event that ${\cal H}^\prime$ is independent of time that the characteristic frequencies are the eigenvalues. In the limit where ${\cal H}^\prime(t)$ varies very slowly compared to 
the natural timescale, then we might use an adiabatic approximation. If on the other hand $\omega_v \sim \omega \sim\omega_0$, we use time dependent perturbation theory.
\subsection{Adiabatic approximation}
As long as $\Pi,\Omega,\omega_0$ vary slowly compared to the unperturbed frequencies, $\omega$ and $\omega_0$
might simply compute the time dependent version of the eigenvalues to determine the time dependence of the spin frequency. In that regime, where $\omega_v\ll\omega$ 
Equation~\ref{eq:umatrix} becomes
\begin{equation}
\begin{pmatrix} a \cr b\end{pmatrix} 
= \begin{pmatrix} \left(\cos({\omega^\prime t/2}) -{i(\omega-\omega_0)\over \omega^\prime}\sin({\omega^\prime t /2})\right)e^{i{\omega\over 2}t} & -\frac{i \Omega-\Pi}{\omega^\prime}\sin(\omega^\prime t /2)e^{i{\omega\over 2}t} \cr
   -\frac{i \Omega+\Pi}{\omega^\prime}\sin(\omega^\prime t /2)e^{-i{\omega\over 2}t}& \left(\cos({\omega^\prime t/2}) +{i(\omega-\omega_0)\over \omega^\prime}\sin({\omega^\prime t /2})\right)e^{-i{\omega\over 2}t}
 \end{pmatrix}
\begin{pmatrix} a_0 \cr b_0\end{pmatrix}\label{eq:umatrixem}
\end{equation}
and
\begin{eqnarray*}
{\omega^\prime}^2&=& (\omega_0(t) - \omega(t))^2 + \Omega(t)^2 + \Pi(t)^2\\
\end{eqnarray*}

\begin{eqnarray*}
{\omega^\prime}^2
&\sim& \left(\frac{e}{mc}\right)^2\left[(a_\mu B)^2(1-{\beta_z(t)}^2) + (a_\mu+\frac{1}{\gamma+1})^2 E^2\beta_{xy}^2\right]\\
&\sim& \left(\frac{e}{mc}\right)^2\left[(a_\mu B)^2(1-{\beta_z(t)}^2) + (a_\mu+\frac{1}{\gamma+1})^2 B^2\left(\frac{\omega_v^2}{c\omega}z_0\cos\omega_v t\right)^2\beta_{xy}^2\right]\\
&\sim& \left(\frac{e}{mc}\right)^2\left[(a_\mu B)^2(1-{\beta_z(t)}^2) + (a_\mu+\frac{1}{\gamma+1})^2 B^2\left(\frac{\beta_z^0\omega_v}{\omega}\cos\omega_v t\right)^2\beta_{xy}^2\right]\\
&\sim& \omega_a^2\left[(1-{\beta_z(t)}^2) + (1+\frac{1}{a_\mu(\gamma+1)})^2 \left(\frac{\beta_z^0\omega_v}{\omega}\cos\omega_v t\right)^2\beta_{xy}^2\right]
\end{eqnarray*}
At the magic momentum,
\begin{eqnarray*}
{\omega^\prime}^2&\sim& \omega_a^2\left(1-(\beta_z^0\sin\omega_v t)^2 + \gamma^2\left(\frac{\beta_z^0\omega_v}{\omega}\cos\omega_v t\right)^2\beta_{xy}^2\right)\\
&\sim& \omega_a^2\left(1-{\beta_z^0}^2\left[\sin^2\omega_v t - (\gamma^2-1)\left(\frac{\omega_v}{\omega}\cos\omega_v t\right)^2\right]\right)\\
&\sim& \omega_a^2\left(1-\frac{{\beta_z^0}^2}{2}\left[1-\cos 2\omega_v t - (\gamma^2-1)\left(\frac{\omega_v}{\omega}\right)^2(1+\cos 2\omega_v t)\right]\right)\\
&\sim& \omega_a^2\left(1-\frac{{\beta_z^0}^2}{2}\left[1-\left(\frac{\omega_v}{\omega}\right)^2(\gamma^2-1)-\cos 2\omega_v t(1+(\gamma^2-1)\left(\frac{\omega_v}{\omega}\right)^2)\right]\right)\\
\end{eqnarray*}
The polarization (see Equations~\ref{eq:polarizationefande} in~\ref{ss:checkeom}) is
\begin{eqnarray*}
\vev{s_x}&=& \left(\cos^2(\opt2) -\frac{(\omega_0-\omega)^2-\Omega^2+\Pi^2}{\op^2}\sin^2(\opt2)\right)\cos\omega t 
 +2\left( \frac{\Omega\Pi}{\op^2}\sin^2(\opt2)-\frac{\omega_0-\omega}{\op}\cos(\opt2)\sin(\opt2)\right)\sin\omega t\nonumber              \\
\vev{s_y}&=& \left(\cos^2(\opt2) -\frac{(\omega_0-\omega)^2-\Omega^2+\Pi^2}{\op^2}\sin^2(\opt2)\right)\sin\omega t 
 -2\left( \frac{\Omega\Pi}{\op^2}\sin^2(\opt2)-\frac{\omega_0-\omega}{\op}\cos(\opt2)\sin(\opt2)\right)\cos\omega t  \nonumber        \\
\vev{s_z}&=& \onehalf(|a|^2-|b|^2)= -\frac{(\omega_0-\omega)\Omega}{\op^2}\sin^2(\opt2) + \frac{\Pi}{\op}\cos(\opt2)\sin(\opt2)
\end{eqnarray*}


Let's expand to second order in $\omega$ and $\Pi$ and rearrange a bit

\begin{eqnarray*}
\vev{s_x}&\sim& \left(\cos^2(\opt2) -\frac{(\omega_0-\omega)^2-\Omega^2+\Pi^2}{(\omega_0-\omega)^2}(1-\frac{\Omega^2+\Pi^2}{(\omega_0-\omega)^2}\sin^2(\opt2)\right)\cos\omega t \\
& & +2\left( \frac{\Omega\Pi}{(\omega_0-\omega)^2}\sin^2(\opt2)-(1-\frac{\Omega^2+\Pi^2}{2(\omega_0-\omega)^2}\cos(\opt2)\sin(\opt2)\right)\sin\omega t\nonumber              \\
&\sim& \left(\cos^2(\opt2)-(1-\frac{2\Omega^2}{(\omega_0-\omega)^2})\sin^2(\opt2)\right)\cos\omega t \\
& & +
2\left( \frac{\Omega\Pi}{(\omega_0-\omega)^2}\sin^2(\opt2)-(1-\frac{\Omega^2+\Pi^2}{2(\omega_0-\omega)^2}\cos(\opt2)\sin(\opt2)\right)\sin\omega t\nonumber              \\
&\sim& \left(\cos(\op t)+\frac{2\Omega^2}{(\omega_0-\omega)^2})\sin^2(\opt2)\right)\cos\omega t \\
& & +
2\left( \frac{\Omega\Pi}{(\omega_0-\omega)^2}\sin^2(\opt2)-(1-\frac{\Omega^2+\Pi^2}{2(\omega_0-\omega)^2}\cos(\opt2)\sin(\opt2)\right)\sin\omega t\nonumber              \\
\vev{s_y}
&\sim& \left(\cos(\op t)+\frac{2\Omega^2}{(\omega_0-\omega)^2})\sin^2(\opt2)\right)\sin\omega t \\
& & 
-2\left( \frac{\Omega\Pi}{(\omega_0-\omega)^2}\sin^2(\opt2)-(1-\frac{\Omega^2+\Pi^2}{2(\omega_0-\omega)^2}\cos(\opt2)\sin(\opt2)\right)\cos\omega t\nonumber              \\
\vev{s_z}&\sim& \onehalf(|a|^2-|b|^2)= -\frac{\Omega}{\omega_0-\omega}\sin^2(\opt2) + \frac{\Pi}{\omega_0-\omega}\cos(\opt2)\sin(\opt2)
\end{eqnarray*}
We used $$\op = ((\omega_0-\omega)^2+\Omega^2+\Pi^2)^{1/2}\sim \sqrt{(\omega_0-\omega)^2 + \Omega^2} +\frac{\Pi^2}{2(\omega_0-\omega)}\equiv\opp + \frac{\Pi^2}{2(\omega_0-\omega)}$$
Next
\begin{eqnarray*}
\bm\beta\cdot\bfs&\sim& \left(\cos\op t +\frac{2\Omega^2}{(\omega_0-\omega)^2}\sin^2(\opt2)\right)\beta_{xy}^0\cos(\beta_z^0\sin\omega_v t)+
\frac{1}{(\omega_0-\omega)}\left(-\Omega\sin^2(\opt2) +\onehalf\Pi\sin(\op t)\right)\sin(\beta_z^0\sin(\omega_v t))\\
&\sim& \left(\cos\opp t - \sin\opp t\sin\frac{\Pi^2}{2(\omega_0-\omega)}t +\frac{2\Omega^2}{(\omega_0-\omega)^2}\sin^2(\opt2)\right)\beta_{xy}^0\cos(\beta_z^0\sin\omega_v t)\\+
& &\frac{1}{(\omega_0-\omega)}\left(-\Omega\sin^2(\opt2) +\onehalf\Pi\sin(\op t)\right)\sin(\beta_z^0\sin(\omega_v t))\\
\end{eqnarray*}

If $\omega_v = \omega$, which is not crazy, then
\begin{eqnarray*}
{\omega^\prime}^2&\sim& \omega_a^2(1-{\beta_z^0}^2+ \gamma^2\left(\beta_z^0\cos\omega_v t\right)^2)\\
{\omega^\prime}^2&\sim& \omega_a^2(1-{\beta_z^0}^2\left[2-\gamma^2-\gamma^2\cos2\omega_vt\right])
\end{eqnarray*}

%Or
%\begin{eqnarray*}
%{\omega^\prime}^2 &\sim&(a+b\sin^2\omega_v t-\omega)^2 +c^2\sin^2\omega_v t+d^2\cos^2\omega_v t\\
%\vev{{\omega^\prime}^2}&\sim&[(a+\onehalf b - \omega)^2 + \onehalf c^2 + \onehalf d^2]
%\end{eqnarray*}
\subsection{Time Dependent Perturbation Theory}
\end{comment}
\section{Radial Magnetic Field}
\label{sec:radial}
Suppose there is a uniform radial field (in addition to the uniform vertical magnetic field), and that the force on the muons due to the radial
field is balanced by an equal and opposite force due to a vertical electric field, presumably the quadrupole field.
\begin{eqnarray*}
\bfB &=& B_r(\sin\theta \bfhx - \cos\theta \bfhy) +B_z\bfhz\\
\bfE&=& E_z\bfhz
\end{eqnarray*}
The velocity is clockwise
$$\bm\beta = \beta(\cos\omega t\ \bfhx - \sin\omega t \ \bfhy)$$ 
so that
\begin{eqnarray*}
\bfB &=& B_r(\sin\omega t\bfhx + \cos\omega t \bfhy) +B_z\bfhz
\end{eqnarray*}
where $\theta=\omega t$, and
$$\bm\beta\times\bfE =E_z \beta(\sin\omega t\  \bfhx +\cos\omega t \ \bfhy)$$
Note that $\bm\beta\times\bfE$ is in the radial direction.
\begin{eqnarray*}
{\bf{\cal B}}&=& \left[(a_\mu+\frac{1}{\gamma}){\bfB} - (a_\mu + \frac{1}{\gamma+1})\bf\beta\times\bfE\right]\\
{\cal H} &=& -\tilde{\bm\mu}\cdot\bf{\cal B}\\
&=& -\frac{e\hbar}{mc}\frac{\bm\sigma}{2}\cdot\bf{\cal B}\\
&=& -\frac{\hbar}{2}\begin{pmatrix}\omega_0 & (\Omega_r-\Pi_r) e^{i\omega t}\cr (\Omega_r-\Pi_r) e^{-i\omega t} & -\omega_0\end{pmatrix}
\end{eqnarray*}
where
\begin{eqnarray*}
\omega_0 &=& \frac{e}{mc}B_z(a_\mu + \frac{1}{\gamma})\\
\omega &=&\frac{e}{mc\gamma}B_z\\
\Omega_r &=&-i\frac{e}{mc}B_r(a_\mu+\frac{1}{\gamma})\\
\Pi_r &=& -i\frac{e}{mc} E_z(a_\mu+ \frac{1}{\gamma+1})\beta
\end{eqnarray*}
By inspection
\begin{eqnarray*}
{\omega^\prime}^2 &=& (\omega_0-\omega)^2+|\Omega_r -\Pi_r|^2\\
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}(a_\mu+\frac{1}{\gamma}) - \frac{E_z\beta}{B_z}(a_\mu + \frac{1}{\gamma+1})\right)^2\right]
\end{eqnarray*}
 As the vertical forces sum to zero, $E_z=\beta B_r$ and
\begin{eqnarray*}
{\omega^\prime}^2
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}\right)^2\left(a_\mu+\frac{1}{\gamma} - \beta^2(a_\mu + \frac{1}{\gamma+1})\right)^2\right]\\
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}\right)^2\left(a_\mu+\frac{1}{\gamma} - \frac{\gamma^2-1}{\gamma^2}(a_\mu + \frac{1}{\gamma+1})\right)^2\right]\\
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}\right)^2\left(\frac{1}{\gamma} +\frac{a_\mu}{\gamma^2} - \frac{\gamma-1}{\gamma^2}\right)^2\right]\\
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}\right)^2\left(\frac{a_\mu}{\gamma^2} + \frac{1}{\gamma^2}\right)^2\right]\\
&=& \left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_r}{B_z}\right)^2\left(\frac{a_\mu+1}{\gamma^2}\right)^2\right]\\
&=& \left(\frac{eB_z}{mc}a_\mu\right)^2\left[1 + \left(\frac{B_r}{B_z}\right)^2\left(\frac{1+1/a_\mu}{\gamma^2}\right)^2\right]\\
\end{eqnarray*}
At the magic momentum, $a_\mu = 1/(\gamma^2-1)$ and
\begin{eqnarray*}
\omega^\prime&=& \left(\frac{eB_z}{mc}a_\mu\right)\left[1 + \left(\frac{B_r}{B_z}\right)^2\left(\frac{1+\gamma^2-1}{\gamma^2}\right)^2\right]^{1/2}\\
\omega^\prime&=& \omega_a\left[1 + \left(\frac{B_r}{B_z}\right)^2\right]^{1/2}
\end{eqnarray*}
as shown previously by Miller\cite{miller_long_rad} by direct solution of the Thomas-BMT equation.
\section{Longitudinal Magnetic Field}
\label{sec:longitudinalfield}
Consider the case where there is a longitudinal magnetic field, (in addition to the vertical field), $\bfB= B_l(\cos\theta \bfhx-\sin\theta\bfhy) + B_z\bfhz$. Since the velocity
is everywhere parallel to the longitudinal field the trajectory of the muon is unaffected. The velocity 
$$\bm\beta=\beta(\cos\omega t\bfhx-\sin\omega t\bfhy)$$
\begin{eqnarray*}
{\bf{\cal B}}&=&\left[(a_\mu+\frac{1}{\gamma})\bfB -a_\mu(\frac{\gamma}{\gamma+1})(\bm\beta\cdot\bfB)\bm\beta\right]\\
&=&\left[(a_\mu+\frac{1}{\gamma})(B_z\bfhz + B_l(\cos\theta \bfhx-\sin\theta\bfhy))
-a_\mu(\frac{\gamma}{\gamma+1})\beta^2 B_l(\cos\omega t \bfhx-\sin\omega t\bfhy)\right]\\
&=&\left[(a_\mu+\frac{1}{\gamma})(B_z\bfhz) + B_l(a_\mu+\frac{1}{\gamma} -a_\mu\frac{\gamma}{\gamma+1}\beta^2)(\cos\omega t \bfhx-\sin\omega t\bfhy)\right]\\
&=&\left[(a_\mu+\frac{1}{\gamma})(B_z\bfhz) + B_l(\frac{1}{\gamma} +a_\mu(1-\frac{\gamma}{\gamma+1}\beta^2))(\cos\omega t \bfhx-\sin\omega t\bfhy)\right]\\
\end{eqnarray*}
Then
\begin{eqnarray*}
{\cal H}&=&-\tilde{\bm\mu}\cdot{\bf{\cal B}}=-\frac{e}{mc}\frac{\hbar}{2}\begin{pmatrix}\omega_0 & \Omega_le^{i\omega t} \cr \Omega_le^{-i\omega t} & \omega_0\end{pmatrix}\\
\end{eqnarray*}
where
\begin{eqnarray*}
\omega_0 &=& \frac{e}{mc}B_z(a_\mu + \frac{1}{\gamma})\\
\omega &=& \frac{e}{mc\gamma}B_z\\
\Omega_l &=& \frac{e}{mc}B_l\left[\frac{1}{\gamma} +a_\mu(1-\frac{\gamma}{\gamma+1}\beta^2)\right]\\
\end{eqnarray*}
\subsection{Uniform Longitudinal Field}
\label{ss:uniformlongitudinal}
In the case that $B_l$ is independent of azimuth its affect on the precession frequency can be determined exactly. 
The eigenvalue in the rotating frame
\begin{eqnarray*}
{\omega^\prime}^2&=& (\omega_0-\omega)^2 + |\Omega_l|^2\\
&=&\left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_l}{B_z}\right)^2\left(\frac{1}{\gamma}+a_\mu(1-\frac{\gamma\beta^2}{\gamma+1})\right)^2\right]\\
&=&\left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_l}{B_z}\right)^2\left(\frac{1}{\gamma}+\frac{a_\mu}{\gamma}\right)^2\right]\\
\end{eqnarray*}
 
At the magic momentum
\begin{eqnarray}
{\omega^\prime}^2
&=&\left(\frac{eB_z}{mc}\right)^2\left[a_\mu^2 + \left(\frac{B_l}{B_z}\right)^2 (\gamma a_\mu)^2\right]\nonumber\\
&=&\omega_a^2\left[1+\left(\frac{B_l}{B_z}\right)^2\gamma^2\right]\nonumber\\
\rightarrow \frac{\Delta\omega^\prime}{\omega_a} & \sim & \onehalf \left(\frac{B_l}{B_z}\right)^2\gamma^2\label{eq:longitudinalexact}
\end{eqnarray}
as shown previously by Miller\cite{miller_long_rad} by direct solution of the Thomas-BMT equation.
Note that we define $\omega_a =\frac{eB_z}{mc}a_\mu$, that is with respect to the vertical component of the magnetic field
rather than the magnitude of the field $|B|=\sqrt{B_z^2+B_l^2}$. The dependence of precession frequency in the rotating frame
on longitudinal field is shown in Figure~\ref{fig:long_dw_dl}.
\begin{figure}
\centering
\includegraphics{figures/long_dw_bl.pdf}
\caption{Fractional shift in $\omega_a$ due to the ring average longitudinal field.}
\label{fig:long_dw_dl}
\end{figure}

The polarization is given by (see section~\ref{ss:linearalgebra} for details).
\begin{eqnarray*}
\vev{s_x}&=& \left(\cos^2(\opt2) -\frac{(\omega_0-\omega)^2-\Omega_l^2}{{\omega^\prime}^2}\sin^2(\opt2)\right)\cos\omega t 
 -2\left(\frac{\omega_0-\omega}{{\omega^\prime}}\cos(\opt2)\sin(\opt2)\right)\sin\omega t              \\
\vev{s_y}&=& -\left(\cos^2(\opt2) -\frac{(\omega_0-\omega)^2-\Omega_l^2}{{\omega^\prime}^2}\sin^2(\opt2)\right)\sin\omega t 
 -2\left(\frac{\omega_0-\omega}{{\omega^\prime}}\cos(\opt2)\sin(\opt2)\right)\cos\omega t              \\
\vev{s_z}&=& \onehalf(|a|^2-|b|^2)= \frac{(\omega-\omega_0)\Omega_l}{{\omega^\prime}^2}\sin^2(\omega^\prime t/2)
\end{eqnarray*}
The projection of the polarization on the direction of motion
\begin{eqnarray}
\bmbetahat\cdot s &=& \left(\cos^2(\opt2) - \frac{1}{{\omega^\prime}^2}[(\omega-\omega_0)^2-\Omega_l^2]\sin^2(\opt2)\right)
\label{eq:longitudinal}
\end{eqnarray}
is shown in Figure~\ref{fig:plot_longitudinal}.
\begin{figure}
\centering
\includegraphics{figures/plot_longitudinal1.pdf}
\includegraphics{figures/plot_longitudinal2.pdf}
\includegraphics{figures/plot_longitudinal3.pdf}
\caption{The projection of the polarization on the direction of motion ($\bmbetahat\cdot\bfs$) with longitudinal magnetic field 10\% of the vertical field and longitudinal field zero
is shown in each of the
three plots. The top plot is computed in simulation by integration of the equations of motion and the BMT equation. The middle plot is the analytic result (Equation~\ref{eq:longitudinal}).
The bottom plot is the simulation and analytic superimposed. The agreement is excellent.}
\label{fig:plot_longitudinal}
\end{figure}
\subsection{Longitudinal Harmonics}
In general the longitudinal field can be expanded in fourier harmonics. Then each harmonic is written as
 $$\bfB_n= b_n\cos n\theta (\cos \theta \bfhx-\sin\theta\bfhy) + B_z\bfhz$$
where $n=0$ is a uniform longitudinal field everywhere tangent to the trajectory and 
\begin{eqnarray*}
\bfB_n\cdot\bm\beta &=&b_n\cos n\omega t\\
{\bf\cal B}&=& (a_\mu+\frac{1}{\gamma})B_z\bfhz+b_n\cos n\omega t\left((\cos \omega t\bfhx-\sin\omega t\bfhy)(a_\mu+\frac{1}{\gamma})-a_\mu\frac{\gamma}{\gamma+1}\beta^2\cos(n\omega t)
(\cos\omega t\bfhx -\sin\omega t\bfhy\right)\\
\bm\sigma\cdot{\bf\cal B}&=&\begin{pmatrix}(a_\mu+\frac{1}{\gamma})B_z & 
b_n\left((a_\mu+\frac{1}{\gamma})\cos(n\omega t)e^{i\omega t}-a_\mu\frac{\gamma}{\gamma+1}\beta^2\cos(n\omega t) e^{i\omega t}\right)\cr
b_n\left((a_\mu+\frac{1}{\gamma})\cos(n\omega t)e^{-i \omega t}-a_\mu\frac{\gamma}{\gamma+1}\beta^2\cos(n\omega t) e^{-i\omega t}\right) &
-(a_\mu+\frac{1}{\gamma})B_z
\end{pmatrix}\\
\frac{e \hbar}{mc}\frac{\bm\sigma}{2}\cdot{\bf\cal B}&=&\frac{\hbar}{2}\begin{pmatrix}\omega_0 & \Omega e^{-i\omega t} \cr \Omega^* e^{i\omega t} & -\omega_0\end{pmatrix}\\
\end{eqnarray*} 
where
\begin{eqnarray*}
\omega_0 &=& B_z\frac{e}{mc}(a_\mu+\frac{1}{\gamma})\\
\Omega &=&\frac{e}{mc}b_n\left((a_\mu+\frac{1}{\gamma})\cos(n\omega t)-a_\mu\frac{\gamma}{\gamma+1}\beta^2\cos(n\omega t) \right)\\
&=&\onehalf X b_n\left( e^{-i n\omega t}+e^{in\omega t}\right)\\
\end{eqnarray*}
and $$X = \frac{e}{mc}(a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\beta^2).$$
Transforming to the rotating frame
\begin{equation}
{\cal H}=-\frac{\hbar}{2}\begin{pmatrix}\omega_0-\omega & \Omega \cr \Omega^*  & -\omega_0+\omega \end{pmatrix}
\end{equation}
In view of the time dependence of $\Omega$, we proceed using time dependent perturbation theory.
Solutions to the unperturbed Hamiltonian ${\cal H}_0 = \frac{\hbar}{2}\begin{pmatrix}\op & 0 \cr 0 & -\op\end{pmatrix}$, with $\op=\omega_0-\omega$
are
\begin{eqnarray*}
a(t)&=&a_0e^{-i\op t/2}\\
b(t)&=& b_0 e^{i\op t/2}
\end{eqnarray*}
The general solution to ${\cal H}$ is
\begin{eqnarray*}
\psi(t) &=& a(t)\psi_a + b(t)\psi_b \\
&=& c_a(t)e^{-i\op t/2}\psi_a + c_b(t)e^{i\op t/2}\psi_b 
\end{eqnarray*}
Suppose in the initial state that $c_a(0)=c_b(0)=\onroot2$. Then in zeroth order $c_a^{(0)}(t)=\frac{1}{\sqrt{2}},c_b^{(0)}(t)=\frac{1}{\sqrt{2}}$.
In first order
\begin{eqnarray*}
c_b^{(1)}(t) &=& -i  \int_0^t c_a^{(0)} \Omega^* e^{i\op \tp}d\tp = -i\frac{1}{\sqrt{2}}\int_0^t \onehalf X b_n\left(e^{-i(n\omega-\op) \tp} +  e^{i(n\omega +\op)\tp}\right) d\tp\\
&=&  -i\frac{1}{\sqrt{2}}\onehalf X b_n\left(\frac{e^{-i(n\omega-\op) t}-1}{-i(n\omega-\op)} +  \frac{e^{i(n\omega +\op)t} -1}{i(n\omega +\op)}\right)\\
c_a^{(1)}(t) &=& -i  \int_0^t c_b^{(0)} \Omega e^{-i\op \tp}d\tp\\
\end{eqnarray*}
%Note (for future reference) that the time average of $c_b^{(1)}$ is 
%$$ \vev{c_b^{(1)}}= -\frac{1}{\sqrt{2}}\onehalf X b_n\left(\frac{1}{-i(n\omega+\op)}+\frac{1}{i(n\omega-\op)}\right)$$
In second order
\begin{eqnarray*}
c_b^{(2)}(t) &=& -i \int_0^t c_a^{(1)}(\tp) \Omega^* e^{i\op \tp}d\tp\\
c_a^{(2)}(t) &=& -i \int_0^t c_b^{(1)}(\tp)  \Omega e^{-i\op \tp}d\tp\\
&=& -\frac{1}{4\sqrt{2}}\int_0^tX^2b_n^2\left(\frac{e^{-i (2n\omega -\op)\tp}-e^{-i n\omega \tp}}{-i(n\omega-\op)}
+ \frac{e^{i(2n\omega+\op) \tp}-e^{i n\omega \tp}}{i(n\omega+\op)}\right.\\
& &\left. +\frac{e^{i\op \tp} -e^{i n\omega \tp}}{-i(n\omega -\op)} +\frac{e^{i\op \tp}-e^{-i n\omega \tp}}{i(n\omega-\op)}\right)e^{-i\op \tp} d\tp\\
&=& -\frac{1}{4\sqrt{2}}\int_0^tX^2b_n^2\left(\frac{e^{-i (2n\omega )\tp}-e^{-i(n\omega +\op) \tp}}{-i(n\omega-\op)}
+ \frac{e^{i(2n\omega) \tp}-e^{i(n\omega -\op)\tp}}{i(n\omega+\op)}\right.\\
& &\left. +\frac{1 -e^{i (n\omega -\op)\tp}}{-i(n\omega -\op)} +\frac{1-e^{-i (n\omega+\op) \tp}}{i(n\omega+\op)}\right) d\tp\\
\end{eqnarray*}
\begin{eqnarray}
c_a^{(2)}(t)
&=& -\frac{1}{4\sqrt{2}}\int_0^tX^2b_n^2\left(\frac{e^{-i (2n\omega )\tp}-e^{-i(n\omega +\op) \tp}}{-i(n\omega-\op)}
+ \frac{e^{i(2n\omega) \tp}-e^{i(n\omega -\op)\tp}}{i(n\omega+\op)}
 +\frac{1 -e^{i (n\omega -\op)\tp}}{-i(n\omega -\op)} +\frac{1-e^{-i (n\omega+\op) \tp}}{i(n\omega+\op)}\right) d\tp\nonumber\\
&=& -\frac{1}{4\sqrt{2}}X^2b_n^2\left(\frac{1}{-i(n\omega-\op)}\left(\frac{e^{-i(2n\omega)t}-1}{-2in\omega}
-\frac{e^{-i(n\omega+\op)t}-1}{-i(n\omega+\op)}+ t -\frac{e^{i(n\omega-\op)t}-1}{i(n\omega-\op)}\right)\right.\nonumber\\
& &\left. +\frac{1}{i(n\omega+\op)}             
\left(\frac{e^{i(2n\omega)t}-1}{2in\omega}
-\frac{e^{i(n\omega-\op)t}-1}{i(n\omega-\op)}+ t -\frac{e^{-i(n\omega+\op)t}-1}{-i(n\omega+\op)}                
\right)\right)\\
&=& -\frac{1}{4\sqrt{2}}X^2b_n^2\left(\frac{1}{-i(n\omega-\op)}\left(\frac{e^{-in\omega t}\sin(n\omega)t}{n\omega}
-\frac{e^{-i(n\omega+\op)t}-1}{-i(n\omega+\op)}+ t -\frac{e^{i(n\omega-\op)t}-1}{i(n\omega-\op)}\right)\right.\nonumber\\
& &\left. +\frac{1}{i(n\omega+\op)}             
\left(\frac{e^{in\omega t}\sin(n\omega)t}{n\omega}
-\frac{e^{i(n\omega-\op)t}-1}{i(n\omega-\op)}+ t -\frac{e^{-i(n\omega+\op)t}-1}{-i(n\omega+\op)}                
\right)\right)\label{eq:nzero}\\
%\vev{c_a^{(2)}(t)}&=& -\frac{X^2b_n^2}{4\sqrt{2}}\left(\frac{2\op}{n\omega^2-\op^2}\vev{t}+ 0\right)
&=& \onroot2\left(-i\frac{\lambda}{2} t + C^{(2)}_{a-rem}\right)\nonumber
\end{eqnarray}
where $\lambda$ is the coefficient of imaginary time $it$ in Equation~\ref{eq:nzero}

\begin{enumerate}
\item[]
If $n=0$, $$\lambda(n=0)=2X^2b_0^2\frac{1 }{\eta},$$ as the terms proportional to $\frac{\sin(n\omega t)}{n\omega}$ reduce
to $t$ when $n=0$.
\item[]
If $n>0$, $$\lambda(n>0)=X^2b_n^2\frac{-\eta }{(n\omega)^2-\eta^2}$$
\end{enumerate}


%If $n=0$ (uniform longitudinal field)
%\begin{eqnarray*}
%c_a^{(2)}(t) &=& -\frac{1}{\sqrt{2}}\int_0^tX^2b_0^2\frac{1-e^{-i\op \tp}}{i\op} d\tp\\
% &=& -\frac{1}{\sqrt{2}}\frac{X^2b_0^2}{i\op}\left(t-\frac{e^{-i\op t}-1}{-i\op}\right) 
%\end{eqnarray*}
%The time average of $c_a^{(2)}$ 
%$$\vev{c_a^{(2)}}= 
%-\frac{1}{\sqrt{2}}\frac{X^2b_0^2}{\op}\left(i\vev{t}+\frac{1}{\op}\right)
%=-i\frac{\opp^2}{\op}\left(\vev{t}+\frac{i}{\op}\right)$$
%increases with time. 
Next recall
\begin{eqnarray*}
c_a(t)&\sim& \left(c_a^{(0)} + c_a^{(1)} + c_a^{(2)}+\ldots \right)\\
%&=&\onroot2\left(1 +\onehalf Xb_n\left(\frac{e^{-i(n\omega+\eta)t}-1}{-i(n\omega+\eta)}+\frac{e^{i(n\omega-\eta)t}-1}{i(n\omega-\eta)}\right)-i\frac{\lambda}{2}t +C^{(2)}_{a-rem}\right)\\
&=&\onroot2\left(1 +c_a^{(1)}-i\frac{\lambda}{2}t +C^{(2)}_{a-rem}\right)\\
&\sim& \onroot2\left(1+c_a^{(1)}+C^{(2)}_{a-rem}\right)e^{-i\lambda t/2}
\end{eqnarray*}
Note that there is no contribution to a frequency shift from $c_a^{(0)}$ or $c_a^{(1)}$ since they contain no term proportional to $t$, whereas $c_a^{(2)}$ does have such a term, 
with coefficient $\lambda/2$. $C^{(2)}_{a-rem}$ is all of the terms remaining in ~\ref{eq:nzero} 
that are not coefficients of $it$, which in turn is different for $n=0$ and $n>0$. Fortunately the value of that remaining term does not
contribute to the frequency shift.

The amplitudes to be in the states $\psi_a$ and $\psi_b$ are
\begin{eqnarray*}
a(t) &=&c_a(t)e^{-i\op t/2 } = \onroot2\left(1+c_a^{(1)}+C^{(2)}_{a-rem}\right)e^{-i(\op+\lambda)t/2 }\\
b(t) &=& c_b(t)e^{i\op t/2 } = \onroot2\left(1+c_b^{(1)}+C^{(2)}_{b-rem}\right)e^{i(\op+\lambda)t/2 }
\end{eqnarray*}
As shown in~\ref{eq:polarization}, the polarization $\expect{\psi}{\bm\sigma/2}{\psi}$ in the x and y directions is proportional to the real and imaginary parts of $a^*b$.
\begin{eqnarray*}
a^* b &=& \onehalf\left(1+c_a^{(1)}+C^{(2)}_{a-rem}\right)^*\left(1+c_b^{(1)}+C^{(2)}_{b-rem}\right)e^{i(\op+\lambda)t }
\end{eqnarray*}
For $n=0$, the precession frequency in the rotating frame at the magic momentum is
\begin{eqnarray}
\omega^\prime&=& \eta+\lambda(n) = \omega_0-\omega +\lambda(n)\nonumber\\
&=&\omega_0-\omega +2\frac{X^2b_n^2(\omega_0-\omega)}{((n\omega)^2-(\omega_0-\omega)^2)\nonumber}\\
&=& \omega_0-\omega + \onehalf\left(\frac{e}{mc}a_\mu\right)^2\left(1+\frac{\gamma^2-1}{\gamma}-\frac{\gamma-1}{\gamma}\right)^2\frac{1}{(\omega_0-\omega)}\nonumber\\
&=& \omega_a + \onehalf \omega_a^2\frac{b_0^2}{B_z^2}\gamma^2\frac{1}{\omega_a}\nonumber\\
&=& \omega_a(1 + \onehalf\frac{b_0^2}{B_z^2}\gamma^2)\nonumber\\
\frac{\Delta\omega^\prime}{\omega_a}&=&\onehalf\frac{b_0^2}{B_z^2}\gamma^2\nonumber
\end{eqnarray}
consistent (to second order) with the exact result in Equation~\ref{eq:longitudinalexact}. 
(Note that for clarity we here use the definition $\omega_a\equiv \frac{eB_z}{mc}a_\mu$ rather than $\omega_a\equiv \frac{e|\bfB|}{mc}a_\mu$.)
If $n>0$ the precession frequency in the rotating frame at the magic momentum is
\begin{eqnarray*}
\omega^\prime&=&  \omega_a - \omega_a^2\frac{b_n^2}{B_z^2}\frac{\gamma^2}{4}\frac{\omega_a}{((n\omega)^2-\omega_a^2)}\\
&=&  \omega_a\left( 1- \frac{b_n^2}{B_z^2}\frac{\gamma^2}{4}\frac{\omega_a^2}{(n\omega)^2-\omega_a^2}\right)\\
\frac{\Delta\omega^\prime}{\omega_a}&=& -\frac{1}{4}\frac{b_n^2}{B_z^2}\gamma^2\frac{\omega_a}{(n\omega)^2-\omega_a^2}
\end{eqnarray*}

\begin{comment}
%The only time dependence that does not average to zero on the time scale of $1/\op$ is that of $c_a^{(2)}(t)$ suggesting that we write,
%$$c_a(t) = e^{-i(\opp^2/\op) t}\left(c_a^{(0)} + c_a^{(1)}(t) + c_a^{(2)}(t)\right)$$
%which is correct to second order in $X b_0$.
%The amplitude $$\VEV{\psi_a}{\psi} = c_a(t)e^{-i\op t}\rightarrow e^{-i(\op^2+(2\opp^2/\op^2))^{1/2} t}\left(c_a^{(0)} + c_a^{(1)}(t) + c_a^{(2)}(t)\right)$$
%Evidently the frequency shift 
$$\Delta\op \sim \opp^2/\op = \frac{X^2b_0^2}{\op}$$
%(Of course the perturbative calculation is only valid as long as $tX^2b_0^2/\op \ll 2\pi$.)
%It is the term proportional to $t$ that determines the frequency shift.


\begin{eqnarray*}
\Delta\omega^\prime (n>0) &\sim& \frac{X^2b_n^2}{2\sqrt{2}}\left(\frac{\op}{(n\omega)^2-\op^2}\right)
\end{eqnarray*}
\end{comment}
It is amusing to see how the contribution to the precession frequency depends on the harmonic number.
To that end consider the ratio of the $n>1$ harmonics to the $n=0$ harmonic
$$\frac{\Delta\omega^\prime(n>0)}{\Delta\omega^\prime(n=0)} = \onehalf\frac{\omega_a^2}{(n\omega)^2-\omega_a^2}\frac{b_n^2}{b_0^2}.$$

\begin{comment}
Let's next evaluate
\begin{eqnarray*}
X b_0 &=& \frac{1}{\sqrt{2}}\frac{e b_0}{mc}(a_\mu+\frac{1}{\gamma} - a_\mu\frac{\gamma}{\gamma+1}\beta^2)\\
&=&\onroot2\omega_a \frac{b_0}{B_z}(1+\frac{1}{a_\mu\gamma} - \frac{\gamma}{\gamma+1}\frac{(\gamma^2-1)}{\gamma^2})\\
&=&\onroot2\omega_a\frac{b_0}{B_z}(1+\frac{1}{a_\mu\gamma} - \frac{(\gamma-1)}{\gamma})\\
&=&\onroot2\omega_a\frac{b_0}{B_z}(\frac{1}{a_\mu\gamma} + (1-\frac{(\gamma-1)}{\gamma})\\
&=&\onroot2\frac{\omega_a}{\gamma}\frac{b_0}{B_z}(1 + \frac{1}{a_\mu})\\
\end{eqnarray*}
At the magic momentum
\begin{eqnarray*}
X^2b_0^2 &=&\onehalf\left(\frac{\omega_a}{\gamma}\frac{b_0}{B_z}\right)^2(1 + \gamma^2-1)^2\\
 &=&\onehalf\left(\frac{\omega_a b_0}{\gamma B_z}\right)^2\left(\gamma^2\right)^2 =
\onehalf\left(\frac{\omega_a b_0}{B_z}\right)^2\gamma^2\\
\rightarrow \frac{\Delta\omega^\prime}{\op} &=&  \onehalf\frac{b_0^2}{B_z^2}\gamma^2
\end{eqnarray*}
in agreement with the expansion of the exact result (Equation~\ref{eq:longitudinalexact}) and shown as a function of $b_l/B_z$ in Figure~\ref{fig:long_dw_dl}.
\end{comment}

And if $b_n=b_0$
\begin{eqnarray*}
 \frac{\Delta\omega^\prime(n >0)}{\Delta\omega^\prime(n=0)}&=&\onehalf\frac{\omega_a^2}{(n\omega)^2-\omega_a^2}
\end{eqnarray*}
Using $\op = \frac{e}{mc}B_za_\mu$ and $\omega =\frac{e}{mc\gamma}B_z$ so that we write $\omega =\frac{ \op}{\gamma a_\mu}$
\begin{eqnarray}
 \Delta\omega^\prime(n>0)&=&\onehalf\frac{-(\gamma a_\mu\omega)^2}{(n\omega)^2-\omega^2(\gamma a_\mu)^2}\Delta\omega^\prime(n=0)\nonumber\\
&=&\onehalf\frac{-(\gamma a_\mu\omega)^2}{(n\omega)^2-\omega^2(\gamma a_\mu)^2}\Delta\omega^\prime(n=0)\nonumber\\
&=&\onehalf\frac{-(\gamma a_\mu)^2}{n^2-(\gamma a_\mu)^2}\Delta\omega^\prime(n=0)\nonumber\\
&=&\onehalf\frac{-\gamma^2}{n^2(\gamma^2-1)^2-\gamma^2}\Delta\omega^\prime(n=0)\label{eq:longitudinalratio}
\end{eqnarray}
(Note that the tuneshift for $n>0$ harmonics has the opposite sign of the $n=0$ term.)
The dependence of frequency shift on longitudinal harmonic is shown in Figure~\ref{fig:long_dw_dn_d0}. The points in 
the plot are from simulation (integration of the Thomas-BMT equation).
\begin{figure}
\centering
\includegraphics{figures/long_dw_dn_d0.pdf}
\caption{$\Delta \omega(n)$ is the frequency shift for the $n^{th}$ harmonic. $\Delta\omega(0)$ is the frequency shift for the $n=0$ harmonic.
The ratio is shown as a function of $n$. $b_n=b_0$ for all $n$. The points are from integration of the Thomas-BMT equation and the equations of motion. 
The open circles are computed for $b_l/B_z=0.01$ and the filled circles
for $b_l/B_z=0.005$. The line is Equation~\ref{eq:longitudinalratio} where $n$ is treated as a continuous variable.}
\label{fig:long_dw_dn_d0}
\end{figure}



