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

\newcommand{\dd}[2]{\frac{\mathrm d #1}{\mathrm d #2}}
\newcommand{\ddh}[2]{\mathrm d #1/\mathrm d #2}
\newcommand{\dds}[2]{\frac{\mathrm d^2 #1}{{\mathrm d #2}^2}}
\newcommand{\ddhp}[3]{\mathrm d^{#3} #1/{\mathrm d #2}^{#3}}
\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{\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{\vecthree}[3]{\left[\begin{array}{c}#1\\#2\\#3\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{\B}{\mathbf B}
\newcommand{\F}{\mathbf F}
\newcommand{\x}{\mathbf x}
\newcommand{\p}{\mathbf p}
\newcommand{\vv}{\mathbf v}
\newcommand{\uv}{\mathbf u}

\begin{document}

\title{Off-Axis Magnetic Fields Extrapolated from On-Axis Multipoles}
\author{Stephen Brooks}
\maketitle

\section{Potential Term}
\newcommand{\er}{\mathbf e_r}
\newcommand{\et}{\mathbf e_\theta}
\newcommand{\ez}{\mathbf e_z}
In cylindrical polar coordinates $(r,\theta,z)$, consider the magnetic scalar potential
\[ \phi = \sin(n\theta+\psi) r^k f(z) \]
for some integers $n,k$, angle $\psi$ and function $f$.  The associated magnetic field is
\[ \B = \nabla\phi = \pdd\phi{r}\er + \frac1r\pdd\phi\theta\et + \pdd\phi{z}\ez = 
\vecthree{\sin(n\theta+\psi) kr^{k-1} f(z)}{n\cos(n\theta+\psi) r^{k-1} f(z)}{\sin(n\theta+\psi) r^k f'(z)}_{r\theta z},
\]
which automatically satisfies $\nabla\times\B = \nabla\times\nabla\phi = \mathbf 0$.  The only remaining condition on $\phi$ is
\[ 0 = \nabla\cdot\B = \nabla\cdot\nabla\phi = \pdds\phi{r} + \frac1r\pdd\phi{r} + \frac1{r^2}\pdds\phi{\theta} + \pdds\phi{z} \]
\[ = \sin(n\theta+\psi) \left(k(k-1)r^{k-2}f(z) + kr^{k-2}f(z) - n^2 r^{k-2} f(z) + r^k f''(z)\right) \]
\[ = \sin(n\theta+\psi) r^{k-2} \left( (k^2-n^2)f(z) + r^2f''(z) \right). \]

\subsection{Long Multipole}
If there is no $z$ behaviour ($f=1$, say), then for this to hold for all $r,\theta$ requires $k^2=n^2$.  Since $k\le0$ cases have a singularity at the origin, put $k=n\ge1$.  This gives the field for an infinitely long multipole:
\[ \phi = \sin(n\theta+\psi) r^n,
\qquad
\B = \vecthree{n\sin(n\theta+\psi) r^{n-1}}{n\cos(n\theta+\psi) r^{n-1}}{0}_{r\theta z}
 = nr^{n-1}\vecthree{\sin((n-1)\theta+\psi)}{\cos((n-1)\theta+\psi)}{0}. \]
$n=1$ corresponds to a dipole, $n=2$ to a quadrupole, etc.  Also $\psi=0$ gives these in their normal orientation and $\psi=\frac\pi2$ in their skew orientation.

\subsection{Normalisation}
Although the rest of the paper will evaluate fields for the potential containing $r^n$ above, this produces fields with magnitude $|\B|=nr^{n-1}$.  If multipole strengths are defined as the values of polynomial coefficients of the field function, then $|\B|=k_nr^n$ may be obtained from the potential $\phi = \sin(n\theta+\psi) \frac{k_n}{n+1}r^{n+1}$.  If strengths are defined as repeated derivatives of the field function, then $|\ddhp{\B}{x}{n}|=d_n$ (i.e. $|\B|=\frac{d_n}{n!}r^n$) may be obtained from $\phi = \sin(n\theta+\psi) \frac{d_n}{(n+1)!}r^{n+1}$.

\section{Series Solution}
If $f''$ is nonzero, the $\nabla\cdot\nabla\phi=0$ equation cannot be satisfied by a single term for all $r$.  However, consider a sum of such terms with the same $n$ but different $k$:
\[ \phi = \sin(n\theta+\psi) \sum_{k=n}^{\infty} r^k f_k(z). \]
Note the lowest term with $k=n$, which dominates the sum near the $r=0$ axis, corresponds to a long $2n$-pole modulated by $f_n(z)$.  The potential must satisfy
\[ 0 = \nabla\cdot\nabla\phi = \sin(n\theta+\psi) \sum_{k=n}^{\infty} r^{k-2} \left( (k^2-n^2)f_k(z) + r^2f''_k(z) \right) \]
for all $\theta$, so the RHS sum must be zero (for all $r,z$).  Equating coefficients of $r^k$ gives
\[ ((k+2)^2-n^2)f_{k+2}(z)+f''_k(z)=0 \]
for $k\ge n$ and the remaining coefficients of $r^{n-2}$ and $r^{n-1}$ are $0$ and $((n+1)^2-n^2)f_{n+1}(z)$ respectively, so $f_{n+1}(z)=0$.  The above relation gives $f_{k+2}$ as a scaled second derivative of $f_k$, so $f_{n+2j+1}(z)=0$ for all $j\ge0$ and 
\[ f_{n+2j}(z) = \left(\prod_{i=1}^j \frac{-1}{(n+2i)^2-n^2}\right) f_n^{(2j)}(z). \]
The coefficients will be written as
\[ C_{nj} = \prod_{i=1}^j \frac{-1}{(n+2i)^2-n^2} = \prod_{i=1}^j \frac{-1}{4(n+i)i} = \frac{(-\frac14)^j n!}{(n+j)!j!}. \]
Setting $f_n=f$, the full potential satisfying Maxwell's equations in free space is:
\[ \phi = \sin(n\theta+\psi) \sum_{j=0}^{\infty} r^{n+2j} f_{n+2j}(z) = \sin(n\theta+\psi) \sum_{j=0}^{\infty} C_{nj} r^{n+2j} f^{(2j)}(z). \]

\subsection{Magnetic Field}
Using the formula for the gradient of a single term ($\sin(n\theta+\psi) r^k f(z)$) given at the start of this note, the magnetic field associated with the series solution is
\[ \B = \nabla\phi = \sum_{j=0}^{\infty} C_{nj} 
\vecthree{\sin(n\theta+\psi) (n+2j)r^{n+2j-1} f^{(2j)}(z)}{n\cos(n\theta+\psi) r^{n+2j-1} f^{(2j)}(z)}{\sin(n\theta+\psi) r^{n+2j} f^{(2j+1)}(z)}_{r\theta z} \]
\[ = \sum_{j=0}^{\infty} C_{nj} r^{n+2j-1} 
\vecthree{(n+2j)\sin(n\theta+\psi) f^{(2j)}(z)}{n\cos(n\theta+\psi) f^{(2j)}(z)}{r \sin(n\theta+\psi) f^{(2j+1)}(z)}_{r\theta z} \]
\[ = \sum_{j=0}^{\infty} C_{nj} r^{n+2j-1} \left(
\vecthree{n\sin(n\theta+\psi) f^{(2j)}(z)}{n\cos(n\theta+\psi) f^{(2j)}(z)}{r \sin(n\theta+\psi) f^{(2j+1)}(z)}_{r\theta z} +
\vecthree{2j\sin(n\theta+\psi) f^{(2j)}(z)}{0}{0}_{r\theta z}
\right) \]
\[ = \sum_{j=0}^{\infty} C_{nj} r^{n+2j-1} \left(
\vecthree{n\sin((n-1)\theta+\psi) f^{(2j)}(z)}{n\cos((n-1)\theta+\psi) f^{(2j)}(z)}{r \sin(n\theta+\psi) f^{(2j+1)}(z)} +
\vecthree{2j\sin(n\theta+\psi)f^{(2j)}(z)\cos\theta}{2j\sin(n\theta+\psi) f^{(2j)}(z)\sin\theta }{0}
\right) \]
\[ = \sum_{j=0}^{\infty} C_{nj} r^{n+2j-1} 
\vecthree{(n\sin((n-1)\theta+\psi)+2j\sin(n\theta+\psi)\cos\theta) f^{(2j)}(z)}{(n\cos((n-1)\theta+\psi)+2j\sin(n\theta+\psi)\sin\theta) f^{(2j)}(z)}{r \sin(n\theta+\psi) f^{(2j+1)}(z)}. \]

\subsection{Computation Without Trigonometric Formulae}
\newcommand{\sfrac}[2]{{\textstyle\frac{#1}{#2}}}
\newcommand{\I}{\mathrm i}
The product formulae for $\sin$ and $\cos$ give 
\[ \sin(n\theta+\psi)\cos\theta = \sfrac12\sin((n-1)\theta+\psi)+\sfrac12\sin((n+1)\theta+\psi) \]
\[ \sin(n\theta+\psi)\sin\theta = \sfrac12\cos((n-1)\theta+\psi)-\sfrac12\cos((n+1)\theta+\psi), \]
which enables the Cartesian field formula to be written as 
\[ \B = \sum_{j=0}^{\infty} C_{nj} r^{n+2j-1} 
\vecthree{((n+j)\sin((n-1)\theta+\psi)+j\sin((n+1)\theta+\psi)) f^{(2j)}(z)}{((n+j)\cos((n-1)\theta+\psi)-j\cos((n+1)\theta+\psi)) f^{(2j)}(z)}{r \sin(n\theta+\psi) f^{(2j+1)}(z)}. \]
Note that every instance of $\cos$ or $\sin(n\theta+\psi)$ is multiplied by a large power of $r$ usually including $r^n$.  This can be used to convert fully to Cartesian coordinates using the complex formulae
\[ re^{\I\theta} = r(\cos\theta+\I\sin\theta) = x+\I y \]
\[ \Rightarrow\quad (re^{\I\theta})^n = r^ne^{\I n\theta} = r^n(\cos n\theta+\I\sin n\theta) = (x+\I y)^n, \]
whence setting $c_n=\mathrm{Re\,}(x+\I y)^n$ and $s_n=\mathrm{Im\,}(x+\I y)^n$ gives
\[ r^n\cos n\theta = c_n \qquad \mbox{and} \qquad r^n\sin n\theta = s_n. \]
These values can be calculated without trigonometry from the recurrence relation
\[ c_0=1; \qquad s_0=0; \qquad c_{n+1}=xc_n-ys_n; \qquad s_{n+1}=yc_n+xs_n, \]
which comes from the definition of complex multiplication $(x+\I y)^{n+1} = (x+\I y)(x+\I y)^n$ in terms of components.  The formula requires a slightly more general form
\[ \hat c_n = r^n\cos(n\theta+\psi) = \mathrm{Re\,} r^ne^{\I n\theta}e^{\I\psi}, \qquad \hat s_n = r^n\sin(n\theta+\psi) = \mathrm{Im\,} r^ne^{\I n\theta}e^{\I\psi}. \]
These satisfy a similar recurrence because the repeated multiplication by $re^{\I\theta} = x+\I y$ starts with $e^{\I\psi}$ rather than $1$:
\[ \hat c_0=\cos\psi; \qquad \hat s_0=\sin\psi; \qquad \hat c_{n+1}=x\hat c_n-y\hat s_n; \qquad \hat s_{n+1}=y\hat c_n+x\hat s_n. \]
The evaluations of $\cos\psi$ and $\sin\psi$ are not a problem because they can be precalculated for each multipole, unlike $\hat c_n$ and $\hat s_n$ that depend on $x,y$.  The field formula can now be rewritten as
\[ \B = \sum_{j=0}^{\infty} C_{nj} r^{2j} 
\vecthree{((n+j)\hat s_{n-1}+j\frac1{r^2}\hat s_{n+1}) f^{(2j)}(z)}{((n+j)\hat c_{n-1}-j\frac1{r^2}\hat c_{n+1}) f^{(2j)}(z)}{\hat s_n f^{(2j+1)}(z)} \]
\[ = \sum_{j=0}^{\infty} C_{nj}
\vecthree{((n+j)(r^2)^j\hat s_{n-1}+j(r^2)^{j-1}\hat s_{n+1}) f^{(2j)}(z)}{((n+j)(r^2)^j\hat c_{n-1}-j(r^2)^{j-1}\hat c_{n+1}) f^{(2j)}(z)}{(r^2)^j\hat s_n f^{(2j+1)}(z)}, \]
where powers of $r^2=x^2+y^2$ have been emphasised as they can be calculated without a square root.  Note that negative powers never arise from $(r^2)^{j-1}$ when $j=0$ because it is multiplied by $j$, so disappears.

\subsection{Computation Using a 2D Field Map}
The only parts of the above formula that depend on $\theta$ are the various $\hat c_n$ and $\hat s_n$, which are relatively quick (and unavoidable) to calculate for each evaluation point.  Calculating the repeated derivatives of $f$ and waiting for the infinite sum to converge are, on the other hand, quite slow but may be precalculated as three functions of just $r$ and $z$ in the following way:
\[ \B = \vecthree{\hat s_{n-1}F(r,z)+\hat s_{n+1}G(r,z)}{\hat c_{n-1}F(r,z)-\hat c_{n+1}G(r,z)}{\hat s_n H(r,z)}, \]
where
\[ F(r,z) = \sum_{j=0}^{\infty} C_{nj} (n+j)(r^2)^j f^{(2j)}(z); \qquad
G(r,z) = \sum_{j=1}^{\infty} C_{nj} j(r^2)^{j-1} f^{(2j)}(z); \]
\[ H(r,z) = \sum_{j=0}^{\infty} C_{nj} (r^2)^j f^{(2j+1)}(z). \]

\subsubsection{Universal Field Map}
Considering just one end of the magnet, if the fringe field is related to a universal fringe field function via $f(z) = \hat f(z/l)$, then $f^{(n)}(z) = l^{-n}\hat f^{(n)}(z/l)$ and
\[ F(r,z) = \sum_{j=0}^{\infty} C_{nj} (n+j)\left(\frac rl\right)^{2j} \hat f^{(2j)}\left(\frac zl\right) = \hat F_n\left(\frac rl,\frac zl\right); \]
\[ G(r,z) = \frac1{r^2} \sum_{j=1}^{\infty} C_{nj} j\left(\frac rl\right)^{2j} \hat f^{(2j)}\left(\frac zl\right) = \frac1{r^2}\hat G_n\left(\frac rl,\frac zl\right); \]
\[ H(r,z) = \frac1r \sum_{j=0}^{\infty} C_{nj} \left(\frac rl\right)^{2j+1} \hat f^{(2j+1)}\left(\frac zl\right) = \frac1r\hat H_n\left(\frac rl,\frac zl\right), \]
where the new functions $\hat F_n$, $\hat G_n$, $\hat H_n$ depend only on $n$ and the form of $\hat f$, not $l$.  In terms of these new functions,
\[ \B = \vecthree{\hat s_{n-1}\hat F_n(r/l,z/l)+\hat s_{n+1}\frac1{r^2}\hat G_n(r/l,z/l)}{\hat c_{n-1}\hat F_n(r/l,z/l)-\hat c_{n+1}\frac1{r^2}\hat G_n(r/l,z/l)}{\hat s_n\frac1r \hat H_n(r/l,z/l)}. \]
Unfortunately use of the universal function $\hat H_n$ incurs a square root to calculate $\frac1r$, which the non-universal forms do not have, assuming they are stored as functions of $r^2$ and $z$.

\subsubsection{Field Map for Magnet Ends Only}
In the interior of a long magnet, the field tends towards the expression given earlier:
\[ \B = nr^{n-1}\vecthree{\sin((n-1)\theta+\psi)}{\cos((n-1)\theta+\psi)}{0}
= \vecthree{n\hat s_{n-1}}{n\hat c_{n-1}}{0}. \]
This corresponds to $F(r,z)=n$ and $G(r,z)=H(r,z)=0$.  The lowest-order terms up to $r^2$ in these sums are:
\[ F(r,z) = C_{n0}nf(z) + C_{n1}(n+1)r^2f''(z) + \dots = nf(z) - \sfrac14 r^2f''(z) + \dots; \]
\[ G(r,z) = C_{n1}f''(z) + C_{n2}2r^2f''''(z) + \dots = \frac{-f''(z)}{4(n+1)} + \frac{r^2f''''(z)}{32(n+1)(n+2)} + \dots; \]
\[ H(r,z) = C_{n0}f'(z) + C_{n1}r^2f'''(z) + \dots = f'(z) - \frac{r^2f'''(z)}{4(n+1)} + \dots, \]
so when $f(z)=1$ and $f'(z)=f''(z)=0$ to some precision, the field calculation may be replaced by the simple multipole expression.  Similarly, when $f=f'=f''=0$ to a good approximation, the field may be given as zero.  This means $F$, $G$, $H$ only need to be calculated for the transition region where $0<f(z)<1$ for each magnet (or end of magnet).

\subsubsection{Magnets with Symmetrical Ends}
It is common to use the same fringe field for both ends of a magnet of length $L$, so that $f(z) = g(z)+g(L-z)$, where $g$ is a sigmoid function going from $-\frac12$ to $\frac12$, representing one end of the magnet (e.g. $g(z)=\frac12\tanh(z/l)$).  If $F$, $G$, $H$ are now calculated using $g$ in place of $f$, the magnetic field is given by
\[ \B = \vecthree{\hat s_{n-1}(F(r,z)+F(r,L-z))+\hat s_{n+1}(G(r,z)+G(r,L-z))}{\hat c_{n-1}(F(r,z)+F(r,L-z))-\hat c_{n+1}(G(r,z)+G(r,L-z))}{\hat s_n (H(r,z)-H(r,L-z))}, \]
noting the sign change for $H(r,L-z)$ because it contains odd derivatives of $g$.

For the regions where $g(z)$ is almost constant at $-\frac12$ or $\frac12$, $F(r,z)=-\frac12 n$ or $\frac12 n$ may be used respectively and $G(r,z)=H(r,z)=0$.  If both $g(z)$ and $g(L-z)$ are in a constant region and they have opposite sign (i.e. far outside the magnet), the $\B$ field will be zero and the calculation should be aborted before calculating $\hat s_{n-1}$ etc.

%\subsection{Series Convergence}

%\section{Solenoids and Other Fields}
% Edge fields?

\end{document}