\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins}

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\pdd}[2]{\frac{\partial #1}{\partial #2}}
\newcommand{\pdds}[2]{\frac{\partial^2 #1}{\partial{#2}^2}}
\newcommand{\pddd}[3]{\frac{\partial^2 #1}{\partial #2\partial #3}}

\newcommand{\de}{\,\mathrm d}

\newcommand{\rowfour}[4]{\left[\begin{array}{cccc}#1&#2&#3&#4\end{array}\right]}
\newcommand{\vectwo}[2]{\left[\begin{array}{c}#1\\#2\end{array}\right]}
\newcommand{\mxtwo}[4]{\left[\begin{array}{cc}#1&#2\\#3&#4\end{array}\right]}

\newcommand{\sq}[1]{\left(#1\right)^2}
\newcommand{\floor}[1]{\left\lfloor#1\right\rfloor}

\newcommand{\rem}[1]{}

\newcommand{\bb}{\mathbf b}
\newcommand{\B}{\mathbf B}
\newcommand{\s}{\mathbf s}
\newcommand{\cc}{\mathbf c}

\newcommand{\e}{\mathrm e}

\begin{document}

\title{Fringe Fields for VFFAG Magnets with Edge Angles}
\author{Stephen Brooks}
\maketitle

\section{Conventions, Assumptions and Definitions}
The beam is travelling in the $+z$ direction and the machine closed orbit shifts in the $+y$ direction (`upwards') with increasing energy.  The standard hard-edged scaling VFFAG field is $B_y=B_0 \e^{ky}$ within the magnet body on the $x=0$ mid-plane \cite{HB2010}, where $k$ has units of inverse length.

For a skew magnet, the edge angle is parametrised by $\tau = \tan\theta$, where $\theta$ is the angle between the magnet edge and the $y$ axis.  Defining $\zeta=z-\tau y$, the function $f(\zeta)$ determines the field fall off at both ends of the magnet, approaching $f=1$ within the magnet and $f=0$ outside.

\section{Fringe Field Issues}
A standard Maxwellian extrapolation \cite{PAC09} away from the $x=0$ plane is insufficient because it only satisfies the three Maxwell equations involving $\partial_x$ and there is also the curl condition $(\nabla\times\B)_x=\partial_yB_z-\partial_zB_y=0$ to satisfy within the plane itself.  This is not a problem for conventional dipoles where only $B_x$ (the perpendicular field component) is nonzero on the midplane.

\subsection{A General Solution for Zero Edge Angle}
Given a desired $B_y$, the curl condition can be enforced explicitly via
\[ \partial_yB_z-\partial_zB_y=0 \qquad\Rightarrow\qquad \partial_yB_z=\partial_zB_y \qquad\Rightarrow\qquad B_z=\int^y \partial_zB_y \de y, \]
where the lower $y$-limit of the integral can be chosen as desired for each value of $z$.  If the goal is to have $B_y=f(z)$, the field components
\[ B_y=f(z), \qquad B_z=yf'(z) \]
satsify $\partial_yB_z=\partial_zB_y=f'(z)$.  For the more general situation $B_y=g(y,z)$, the field
\[ B_y=g(y,z), \qquad B_z=\int^y \partial_zg(y,z) \de y \]
satisfies $\partial_yB_z=\partial_zB_y=\partial_zg(y,z)$.

\subsection{Extension to Angled Fringe Fields}
Although the solution in the previous section can be used with $g(y,z)$ equal to the desired slanted field $f(\zeta)$, this produces fringe fields that propagate vertically rather than parallel to the magnet edge.  However, consider the field
\[ B_y=f(\zeta)-\tau yf'(\zeta), \qquad B_z=yf'(\zeta). \]
Since $\partial_y f(\zeta)=-\tau f'(\zeta)$ and $\partial_z f(\zeta)=f'(\zeta)$, the curl condition is satisfied with $\partial_yB_z=\partial_zB_y=f'(\zeta)-\tau yf''(\zeta)$.  Can this be generalised in an analogous way to the zero edge angle case, by replacing a multiplication by $y$ with an integral in $y$?  For a function $g(y,\zeta)$ this would suggest the field
\[ B_y=g(y,\zeta)-\tau\int^y \partial_\zeta g(y,\zeta)\de y, \qquad B_z=\int^y \partial_\zeta g(y,\zeta)\de y. \]
Checking the derivatives gives 
\[ \partial_yB_z=\partial_zB_y = \partial_\zeta g(y,\zeta) - \tau \int^y \partial^2_\zeta g(y,\zeta)\de y, \]
noting firstly that $\zeta$ depends on $y$, so $g(y,\zeta)$ depends on $y$ twice and secondly that $\partial_z=\partial_\zeta$ since $y$ is held constant in both cases.  The lower limit of the $B_z$ integral is still free to be chosen, this time as a function of $\zeta$ and its evaluation can be reused in $B_y=g(y,\zeta)-\tau B_z$.

\subsection{Scalar Potential Formulation}
The fringe field issue is caused by having to satisfy $\nabla\times\B=\mathbf0$ in free space.  However if $\B=\nabla\phi$, the curl of a gradient is zero so this equation always holds.  The four fields given in this section have the scalar potentials below (valid only on the $x=0$ plane).
\begin{eqnarray*}
\phi_{f0} &=& yf(z) \\
\phi_{g0} &=& \int^y g(y,z) \de y \\
\phi_{f\theta} &=& yf(\zeta) \\
\phi_{g\theta} &=& \int^y g(y,\zeta)\de y
\end{eqnarray*}

\section{The Scaling VFFAG Mid-Plane Field}
Using the angled fringe field formula with $g(y,\zeta) = B_0 \e^{ky} f(\zeta)$ yields
\[  B_z=\int^y B_0 \e^{ky} f'(\zeta)\de y = B_0 f'(\zeta)\ \frac{\e^{ky}-\e^{kY(\zeta)}}{k}, \]
where $Y(\zeta)$ is the lower limit of the integral.  Because the integral is bounded in the negative $y$ direction (assuming $k>0$), the most elegant choice here is $Y(\zeta)=-\infty$, so that
\[ B_z=B_0 \e^{ky} \frac{f'(\zeta)}{k}, \qquad B_y = B_0 \e^{ky} f(\zeta) - \tau B_z = B_0 \e^{ky}\left(f(\zeta) - \frac{\tau f'(\zeta)}{k}\right). \]

\section{Off-plane Field Expansion}
This section adapts the magnetic field expansion used in \cite{PAC09} for the vertical orbit-excursion case.  Maxwell's equations in free space, $\nabla\cdot\B=\nabla\times\B=0$, can be rearranged to give
%\setlength{\arraycolsep}{3pt}
\[\partial_x\B = \left[\begin{array}{ccc}0&-\partial_y&-\partial_z\\\partial_y&0&0\\\partial_z&0&0\end{array}\right]\B
= \left[\begin{array}{cc}0&-\nabla_{y,z}^T\\\nabla_{y,z}&0\end{array}\right]\B
, \]
from which a Taylor expansion in $x$ yields
\[ B_x(x,y,z)=\sum_{n=0}^\infty \frac{x^{2n+1}}{(2n+1)!}(-\nabla_{y,z}^2)^n(-\nabla_{y,z}\cdot\B_{y,z}(0,y,z)) \]
\[ \B_{y,z}(x,y,z)=\B_{y,z}(0,y,z) + \sum_{n=1}^\infty \frac{x^{2n}}{(2n)!}\nabla_{y,z}(-\nabla_{y,z}^2)^{n-1}(-\nabla_{y,z}\cdot\B_{y,z}(0,y,z)). \]
Here it has been assumed that $B_x(0,y,z)=0$.  Since our field has a scalar potential on the midplane, $\B_{y,z}(0,y,z)=\nabla_{y,z}\phi(y,z)$, which simplifies the formulae further: 
\[ B_x(x,y,z)=\sum_{n=0}^\infty \frac{x^{2n+1}}{(2n+1)!}(-\nabla_{y,z}^2)^{n+1}\phi(y,z) \]
\[ \B_{y,z}(x,y,z)=\sum_{n=0}^\infty \frac{x^{2n}}{(2n)!}\nabla_{y,z}(-\nabla_{y,z}^2)^n\phi(y,z). \]
For the VFFAG mid-plane field, $\phi(y,z) =  B_0 \e^{ky} \frac{f(\zeta)}k$.

\subsection{The Iterated 2D Laplacian}
The bulk of the field calculation is now in evaluating $(-\nabla_{y,z}^2)^n\phi(y,z)$ from the summands in the previous section.  Assuming $f$ and its derivatives are evaluated at $\zeta$ unless otherwise stated, it is useful to evaluate
\begin{eqnarray*}
 -\nabla_{y,z}^2 ( \e^{ky} f^{(n)} ) &=& -\partial_y\partial_y( \e^{ky} f^{(n)} ) - \partial_z\partial_z( \e^{ky} f^{(n)} ) \\
&=& -\partial_y( k\e^{ky} f^{(n)} - \tau\e^{ky} f^{(n+1)} ) - \partial_z( \e^{ky} f^{(n+1)} ) \\
&=& -( k^2\e^{ky} f^{(n)} - 2k\tau\e^{ky} f^{(n+1)} + \tau^2\e^{ky} f^{(n+2)}) - \e^{ky} f^{(n+2)} \\
&=& \e^{ky}(-k^2 f^{(n)} + 2k\tau f^{(n+1)} - (1+\tau^2) f^{(n+2)}).
\end{eqnarray*}
Therefore the general form is
\[ (-\nabla_{y,z}^2)^n\phi(y,z) = B_0 \e^{ky} \sum_{j=0}^{2n} a_{nj} f^{(j)} \]
and the recurrence relation is
\[ a_{n+1,j} = -k^2 a_{nj} + 2k\tau a_{n,j-1} - (1+\tau^2) a_{n,j-2}, \]
where the definition is made that $a_{n,-1}=a_{n,-2}=0$.  Initially, $a_{00}=\frac1k$ and $a_{0j}=0$ elsewhere.

\subsection{Field Evaluation Formulae for General $f(\zeta)$}
Noting that
\[ \nabla_{y,z}(-\nabla_{y,z}^2)^n\phi(y,z) = B_0 \e^{ky} \sum_{j=0}^{2n} a_{nj} \vectwo{kf^{(j)}-\tau f^{(j+1)}}{f^{(j+1)}}, \]
the field formulae become
\[ B_x(x,y,z)=B_0 \e^{ky}\sum_{n=0}^\infty \frac{x^{2n+1}}{(2n+1)!}\sum_{j=0}^{2n+2} a_{n+1,j} f^{(j)}, \]
\[ \B_{y,z}(x,y,z)= B_0 \e^{ky}\sum_{n=0}^\infty \frac{x^{2n}}{(2n)!}\sum_{j=0}^{2n} a_{nj} \vectwo{kf^{(j)}-\tau f^{(j+1)}}{f^{(j+1)}}. \]
Here, $a_{nj}$ can be precalculated for each magnet (specific $k$, $\tau$ values) but the derivatives $f^{(j)}$ need to be found for each point ($\zeta$ value).  After the scaling factor $B_0 \e^{ky}$, the rest of the formula is a function of $x$ and $\zeta$ only, so it suffices to produce a 2D field map and use the scaling law
\[ \B(x,y,z) = \e^{ky}\ \B(x,0,z-\tau y). \]


\section{$\tanh$ Fringe Field Derivatives}
Given a fringe field extent $l$, the field fall off function can have the form
\[ f(\zeta) = \sigma(\zeta/l)-\sigma((\zeta-L)/l), \]
where $\sigma(x)$ is a unit sigmoid function going from $0$ to $1$.  The $n^\mathrm{th}$ derivative is
\[ f^{(n)}(\zeta) = l^{-n}\left(\sigma^{(n)}(\zeta/l)-\sigma^{(n)}((\zeta-L)/l)\right). \]
A common sigmoid function with exponential fall off characteristics (e.g.\ for iron fringe field shielding) is $\sigma(x)=\frac12+\frac12\tanh x$.  \rem{Since $\tanh x = (\e^x-\e^{-x})/(\e^x+\e^{-x})$, this makes $\sigma(x)=\e^x/(\e^x+\e^{-x})=\e^{2x}/(\e^{2x}+1)$.  }
With this choice of $\sigma(x)$, the constant offset of $\frac12$ cancels in the expression for $f$:
\[ f^{(n)}(\zeta) = \frac{l^{-n}}{2}\left(\tanh^{(n)}(\zeta/l)-\tanh^{(n)}((\zeta-L)/l)\right). \]

\subsection{Computational Simplification}
Usefully, $\tanh$ satisfies the relation
\[ \dd{\tanh x}{x} = 1 - \tanh^2 x \quad \Rightarrow \quad \dd{\tanh^n x}{x} = n \tanh^{n-1} x (1 - \tanh^2 x) = n \tanh^{n-1} x - n \tanh^{n+1} x, \]
which can be used to generate a recurrence for $\tanh^{(n)}x$ in terms of powers of $\tanh x$:
\[ \tanh^{(n)} x = \sum_{j=0}^{n+1} t_{nj} \tanh^j x \]
\[ t_{01}=1, \quad t_{0j}=0 \quad (j\ne1), \qquad t_{n+1,j}=(j+1)t_{n,j+1}-(j-1)t_{n,j-1}. \]
Substituting this into the formula for $f^{(n)}$ gives
\[ f^{(n)}(\zeta) = \frac{l^{-n}}{2}\sum_{j=0}^{n+1} t_{nj}\left(\tanh^j(\zeta/l)-\tanh^j((\zeta-L)/l)\right) = \frac{l^{-n}}{2}\sum_{j=0}^{n+1} t_{nj}(T_1^j-T_2^j), \]
where $T_1=\tanh(\zeta/l)$ and $T_2=\tanh((\zeta-L)/l)$ can be precalculated for each $\zeta$.

\subsection{Calculating the Field Map}
Precalculate $t_{nj}$ to the desired order and $a_{nj}$ for this specific magnet.  For each value of $\zeta$ calculate $T_1$ and $T_2$, then the vector of $T_1^n-T_2^n$ values that gives the vector of $f^{(n)}$ after a linear transformation using $t_{nj}$.  Calculate two further linear transformations using $a_{nj}$:
\[ g_n=\sum_{j=0}^{2n} a_{nj} f^{(j)} \qquad \mbox{and} \qquad h_n=\sum_{j=0}^{2n} a_{nj} f^{(j+1)}. \]
Assuming the field map is being made at $y=0$, for each $x$ value, the field components are:
\[ B_x(x,0,z)=B_0 \sum_{n=0}^\infty \frac{x^{2n+1}}{(2n+1)!}\ g_{n+1} \]
\[ \B_{y,z}(x,0,z)= B_0 \sum_{n=0}^\infty \frac{x^{2n}}{(2n)!} \vectwo{kg_n-\tau h_n}{h_n}. \]


\begin{thebibliography}{9}   % Use for  1-9  references
%\begin{thebibliography}{99} % Use for 10-99 references
\newcommand{\link}[1]{\texttt{#1}}
\bibitem{HB2010} \emph{Vertical Orbit Excursion FFAGs}, S.J. Brooks, Proc. HB2010.
\bibitem{PAC09} \emph{Extending the Energy Range of 50\,Hz Proton FFAGs}, S.J.~Brooks, Proc. PAC 2009.
\end{thebibliography}

\end{document}