\def\op{{\omega^\prime}}
\def\opt{{\omega^\prime t}}
\def\opt2{{\omega^\prime t/2}}
%\section{Longitudinal Magnetic Field}
%\label{sec:longitudinalfield}
We consider the case where in addition to the nominal vertical magnetic field there is a longitudinal component. The total field $\bfB= B_l(\cos\theta \bfhx-\sin\theta\bfhy) + B_z\bfhz$. 
Since the velocity of a muon circulating at the magic radius
is everywhere parallel to the longitudinal field, the particle trajectory is unaffected. The velocity 
$$\bm\beta=\beta(\cos\omega t\bfhx-\sin\omega t\bfhy)$$ is determined by the vertical component.
(The circular orbit is in the x-y plane, $\theta=\omega t$.)
\section*{Uniform longitudinal field}
If $B_l$ is uniform, that is independent of azimuth,
 the Thomas-BMT equation can be solved exactly (see \href{https://gm2-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=25607}{GM2-doc-25607}).
The precession frequency measured 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*}
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*}
 
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}
Note that here 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 (solution to Thomas-BMT) is given by 
\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}
\section*{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, discussed above.
As shown in \href{https://gm2-docdb.fnal.gov/cgi-bin/sso/ShowDocument?docid=25607}{GM2-doc-25607}, 
the Thomas-BMT equation can be derived from an effective Hamiltonian ${\cal {H}}$ representing the interaction energy of magnetic moment and magnetic fields.
The solution to Schrodinger's equation 
$$\i\hbar\partder{}{t}\psi = {\cal H}\psi$$
is equivalent to the solution of the Thomas-BMT equation. 
If the longitudinal field is uniform, ${\cal H}$ is independent of time (in the rotating frame) and the Schrodinger equation can be solved exactly.
If time dependent (nonuniform longitudinal field), it
can be solved using time dependent perturbation theory.

The frequency shift, to second order in the perturbation $b_n/B_z$ is (there is no first order contribution to the precession frequency)
\begin{enumerate}
\item[]
If $n=0$, $$\lambda(n=0)=2X^2b_0^2\frac{1 }{\eta},$$ 
\item[]
If $n>0$, $$\lambda(n>0)=X^2b_n^2\frac{-\eta }{(n\omega)^2-\eta^2}$$
\end{enumerate}
where $$X=\frac{e}{mc}(a_\mu+\frac{1}{\gamma}-a_\mu\frac{\gamma}{\gamma+1}\beta^2).$$
and $\eta =\omega_0-\omega= \omega_a$

%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. 
For $n=0$, the precession frequency in the rotating frame at the magic momentum is
\begin{eqnarray}
\omega^\prime&=& \omega_a+\lambda(n) = \omega_0-\omega +\lambda(n)\nonumber\\
&=&\omega_a +2\frac{X^2b_n^2(\omega_0-\omega)}{((n\omega)^2-(\omega_0-\omega)^2)\nonumber}\\
&=& \omega_a + \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}

\section*{Summary}
For a magic momentum muon circulating on its closed orbit, the shift in frequency due to a uniform ($n=0$) longitudinal field is
$$\frac{\Delta\op}{\omega_a}\sim\onehalf\frac{b_0^2}{B_z^2}\gamma^2$$
and the frequency shift for $n>0$ 
\begin{eqnarray*}
\frac{\Delta\op}{\omega_a}&\sim&-\frac{1}{4}\frac{b_n^2}{B_z^2}\gamma^2\frac{\omega_a^2}{(n\omega)^2-\omega_a^2}\\
&\sim&-\frac{1}{4}\frac{b_n^2}{B_z^2}\gamma^2\frac{\gamma^2}{(n^2(\gamma^2-1)^2-\gamma^2}\\
&\sim&-\frac{1}{4}\frac{b_n^2}{B_z^2}\frac{1}{n^2}
\end{eqnarray*}
where
 $$\bfB= \sum_n b_n\cos n\theta (\cos \theta \bfhx-\sin\theta\bfhy) + B_z\bfhz$$


