\documentclass[a4paper,11pt]{article}
\parskip=1.5ex
\usepackage{verbatim,a4margins,amsmath,amssymb,bm,graphicx,doi}

\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{\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{\abs}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\A}{\mathbf A}
\newcommand{\E}{\mathbf E}
\newcommand{\B}{\mathbf B}
\newcommand{\F}{\mathbf F}
\newcommand{\J}{\mathbf J}

\newcommand{\x}{\mathbf x}
\newcommand{\y}{\mathbf y}
\newcommand{\z}{\mathbf z}

\newcommand{\p}{\mathbf p}
\newcommand{\q}{\mathbf q}
\newcommand{\s}{\mathbf s}
\newcommand{\rv}{\mathbf r}

\newcommand{\uv}{\mathbf u}
\newcommand{\vv}{\mathbf v}

\newcommand{\e}{\mathbf e}
\newcommand{\f}{\mathbf f}
\newcommand{\g}{\mathbf g}

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

\newcommand{\st}{$^\mathrm{st}$ }
\newcommand{\nd}{$^\mathrm{nd}$ }
\newcommand{\rd}{$^\mathrm{rd}$ }
\newcommand{\thh}{$^\mathrm{th}$ }

\newcommand{\divv}{\nabla\cdot}
\newcommand{\curl}{\nabla\times}

\begin{document}

\title{Off-Axis RF Fields Extrapolated from On-Axis Values}
\author{Stephen Brooks}
\maketitle

\section{Maxwell's Equations for an RF Field}
Start with Maxwell's equations
\begin{eqnarray*}
\divv\E &=& \rho/\epsilon_0 \\
\curl\E &=& -\pdd{\B}{t} \\
\divv\B &=& 0 \\
\curl\B &=& \mu_0\left(\J + \epsilon_0\pdd{\E}{t}\right),
\end{eqnarray*}
set the sources $\rho,\J$ to zero and look for a periodic solution with $\E=\E_0\cos\omega t$, $\B=\B_0\sin\omega t$ where $\E_0,\B_0$ do not depend on time.  This gives
\begin{eqnarray*}
\divv\E_0 &=& 0 \\
\curl\E_0\cos\omega t = -\pdd{\B}{t} = -\omega\B_0\cos\omega t  \quad\Rightarrow\quad \curl\E_0 &=& -\omega\B_0 \\
\divv\B_0 &=& 0 \\
\curl\B_0\sin\omega t = \mu_0\epsilon_0\pdd{\E}{t} = \frac1{c^2}(-\omega\E_0\sin\omega t) \quad\Rightarrow\quad \curl\B_0 &=& -\frac{\omega}{c^2}\E_0.
\end{eqnarray*}
Taking the curl of the second equation, notice that 
\[ \curl(\curl\E_0) = -\omega\curl\B_0 = \frac{\omega^2}{c^2}\E_0. \]
It turns out that finding an $\E_0$ that solves this equation is sufficient to satisfy all four, because $\B_0$ can be defined as $-\frac{1}{\omega}\curl\E_0$ and then both $\E_0,\B_0$ are expressed as the curl of another vector field, making their divergences zero, satisfying the other two equations.

\section{Case with No Angular Dependency}
The aim is to build a solution for $\E_0$ out of a sum of simpler functions.  Defining $r=\sqrt{x^2+y^2}$, consider a longitudinal ($z$ directed) vector field $\mathbf a$ with $r^k$ scaling in the radial direction:
\[ a_z=r^kf(z), \quad a_x=a_y=0. \]
The curl of this has two nonzero components:
\begin{eqnarray*}
(\curl\mathbf a)_x &=& \partial_ya_z = kr^{k-1}\pdd{r}{y}f(z) = kr^{k-2}yf(z), \\
(\curl\mathbf a)_y &=& -\partial_xa_z = -kr^{k-1}\pdd{r}{x}f(z) = -kr^{k-2}xf(z). \\
\end{eqnarray*}
Now consider a similar term $\mathbf b$ that points azimuthally around the axis:
\[ b_x=r^{k-1}yf(z), \quad b_y=-r^{k-1}xf(z), \quad b_z=0. \]
Its curl components are:
\begin{eqnarray*}
(\curl\mathbf b)_x &=& -\partial_zb_y = r^{k-1}xf'(z), \\
(\curl\mathbf b)_y &=& \partial_zb_x = r^{k-1}yf'(z), \\
(\curl\mathbf b)_z &=& \partial_xb_y - \partial_yb_x \\
&=& -r^{k-1}f(z) -(k-1)r^{k-3}x^2f(z) - r^{k-1}f(z) - (k-1)r^{k-3}y^2f(z) \\
&=& -2r^{k-1}f(z) -(k-1)r^{k-1}f(z) \\ &=& -(k+1)r^{k-1}f(z).
\end{eqnarray*}
Finally, define a radial term $\mathbf c$:
\[ c_x=r^{k-1}xf(z), \quad c_y=r^{k-1}yf(z), \quad c_z=0, \]
that has curl components
\begin{eqnarray*}
(\curl\mathbf c)_x &=& -\partial_zc_y = -r^{k-1}yf'(z), \\
(\curl\mathbf c)_y &=& \partial_zc_x =r^{k-1}xf'(z), \\
(\curl\mathbf c)_z &=& \partial_xc_y - \partial_yc_x \\
&=& (k-1)r^{k-3}xyf(z) - (k-1)r^{k-3}xyf(z) \\
&=& 0.
\end{eqnarray*}
As a shorthand, subscripts will indicate substitutions of $k$ or $f$ in these terms, omitting any that are unchanged.  In this notation,
\begin{eqnarray*}
\curl\mathbf a &=& k\mathbf b_{k-1}, \\
\curl\mathbf b &=& -(k+1)\mathbf a_{k-1} + \mathbf c_{f'}, \\
\curl\mathbf c &=& -\mathbf b_{f'}.
\end{eqnarray*}
Now the double curls may be evaluated:
\begin{eqnarray*}
\curl(\curl\mathbf a) &=& k\left( -k\mathbf a_{k-2} + \mathbf c_{k-1,f'} \right) \\
&=& -k^2\mathbf a_{k-2} + k\mathbf c_{k-1,f'}, \\
\curl(\curl\mathbf b) &=& -(k+1)\left( (k-1)\mathbf b_{k-2} \right) + \left( -\mathbf b_{f''} \right) \\
&=& -(k^2-1)\mathbf b_{k-2} - \mathbf b_{f''}, \\
\curl(\curl\mathbf c) &=& -\left( -(k+1)\mathbf a_{k-1,f'} + \mathbf c_{f''} \right) \\
 &=& (k+1)\mathbf a_{k-1,f'} - \mathbf c_{f''}.
\end{eqnarray*}
Analogously with the magnetic field solutions in \cite{BrooksOffAxis}, seek a solution of the form
\[ \E_0 = \sum_{k=0}^\infty \mathbf a_{k,f_k} + \mathbf b_{k,g_k} + \mathbf c_{k,h_k} \]
for some functions $f_k,g_k,h_k$.  The double curl is
\begin{eqnarray*}
\curl(\curl\E_0) &=& \sum_{k=0}^\infty 
-k^2\mathbf a_{k-2,f_k} + k\mathbf c_{k-1,f'_k} \\
&-&(k^2-1)\mathbf b_{k-2,g_k} - \mathbf b_{g''_k} \\
&+&(k+1)\mathbf a_{k-1,h'_k} - \mathbf c_{h''_k}.
\end{eqnarray*}
Equating this to $\frac{\omega^2}{c^2}\E_0$ and only taking the order $k$ terms gives
\begin{eqnarray*}
 \frac{\omega^2}{c^2}( \mathbf a_{f_k} + \mathbf b_{g_k} + \mathbf c_{h_k} ) 
&=& -(k+2)^2\mathbf a_{f_{k+2}} + (k+1)\mathbf c_{f'_{k+1}} \\
&-&((k+2)^2-1)\mathbf b_{g_{k+2}} - \mathbf b_{g''_k} \\
&+&(k+2)\mathbf a_{h'_{k+1}} - \mathbf c_{h''_k}.
\end{eqnarray*}
Equating the longitudinal, azimuthal and radial vector coefficients gives
\begin{eqnarray*}
\frac{\omega^2}{c^2}f_k &=&  -(k+2)^2f_{k+2} + (k+2)h'_{k+1}, \\
\frac{\omega^2}{c^2}g_k &=& -((k+2)^2-1)g_{k+2} - g''_k, \\ 
\frac{\omega^2}{c^2}h_k &=& (k+1)f'_{k+1} - h''_k.
\end{eqnarray*}
At this point, the first two equations can be rearranged to give $f_{k+2}$ and $g_{k+2}$ in terms of lower-numbered functions, as required to calculate them all via a recurrence.  However, the highest function in the third equation is $f'_{k+1}$ and not $h$.  Fortunately, some careful cancellation can get rid of both the $f'_{k+1}$ and $h''_k$ terms.

Taking the first equation for the $k-1$ case and taking the $z$ derivative gives
\[ \frac{\omega^2}{c^2}f'_{k-1} =  -(k+1)^2f'_{k+1} + (k+1)h''_k. \]
Multiplying the third equation by $-(k+1)$ gives
\[ -(k+1)\frac{\omega^2}{c^2}h_k = -(k+1)^2f'_{k+1} + (k+1)h''_k. \]
The right hand sides are now identical, implying that 
\[ f'_{k-1} =  -(k+1)h_k, \]
which allows $h$ to be eliminated entirely.  In particular,
\[ f''_k =  -(k+2)h'_{k+1}, \]
which allows writing the first two equations as:
\begin{eqnarray*}
\frac{\omega^2}{c^2}f_k &=&  -(k+2)^2f_{k+2} - f''_k, \\
\frac{\omega^2}{c^2}g_k &=& -((k+2)^2-1)g_{k+2} - g''_k. 
\end{eqnarray*}
Now everything can be written as simple recurrence formulae
\begin{eqnarray*}
f_{k+2} &=& \frac{-1}{(k+2)^2}\left(\frac{\omega^2}{c^2}f_k + f''_k\right), \\
g_{k+2} &=& \frac{-1}{(k+2)^2-1}\left(\frac{\omega^2}{c^2}g_k + g''_k\right), \\
h_{k+1} &=& \frac{-1}{k+2}f'_k. 
\end{eqnarray*}
To avoid discontinuities, the lowest allowed nonzero functions are $f_0$ for the on-axis electric field and $g_1$ and $h_1$.

\subsection{3D Field from On-Axis Electric Field}
In the case where just the on-axis electric field $f_0(z)=f(z)$ is specified, azimuthal components $g_k$ can be set to zero as their recurrence relation is independent of the rest.  $f_k$ and $h_k$ are all linear combinations of derivatives of $f$.  Specifically, for all $k$, $f_{2k+1}=0$, $h_{2k}=0$ and 
\begin{eqnarray*}
f_{2k} &=& \left(\prod_{j=1}^k\frac{-1}{(2j)^2}\right)\left(\frac{\omega^2}{c^2}+\pdds{}{z}\right)^kf \\
 &=& \frac{(-1)^k}{4^k(k!)^2}\sum_{j=0}^k\binom{k}{j}\left(\frac{\omega^2}{c^2}\right)^{k-j}f^{(2j)}, \\
h_{2k+1} &=& \frac{-1}{2k+2}f'_{2k} \\
 &=& \frac{(-1)^{k+1}}{2^{2k+1}k!(k+1)!}\sum_{j=0}^k\binom{k}{j}\left(\frac{\omega^2}{c^2}\right)^{k-j}f^{(2j+1)}.
\end{eqnarray*}
The field shape explicitly is:
\begin{eqnarray*}
 \E_0 &=& \sum_{k=0}^\infty \mathbf a_{k,f_k} + \mathbf c_{k,h_k} \\
  &=& \sum_{k=0}^\infty [r^{k-1}xh_k(z),r^{k-1}yh_k(z),r^kf_k(z)].
\end{eqnarray*}
The magnetic component is:
\begin{eqnarray*}
\B_0 &=& -\frac{1}{\omega}\curl\E_0 \\
&=&  -\frac{1}{\omega}\sum_{k=0}^\infty \curl\mathbf a_{k,f_k} + \curl\mathbf c_{k,h_k} \\
&=&  -\frac{1}{\omega}\sum_{k=0}^\infty k\mathbf b_{k-1,f_k} - \mathbf b_{k,h'_k} \\
&=&  -\frac{1}{\omega}\sum_{k=0}^\infty (k+1)\mathbf b_{k,f_{k+1}} - \mathbf b_{k,h'_k} \\
&=&  -\frac{1}{\omega}\sum_{k=0}^\infty ((k+1)f_{k+1}(z)-h'_k(z)) [r^{k-1}y,-r^{k-1}x,0].
\end{eqnarray*}
The last expression can be simplified slightly using the recurrence relations for $f$ and $h$:
\begin{eqnarray*}
(k+1)f_{k+1} &=&  \frac{-1}{k+1}\left(\frac{\omega^2}{c^2}f_{k-1} + f''_{k-1}\right) \\
h'_k &=& \frac{-1}{k+1}f''_{k-1} \\
\Rightarrow \qquad (k+1)f_{k+1} - h'_k &=& \frac{-1}{k+1}\frac{\omega^2}{c^2}f_{k-1}.
\end{eqnarray*}
Thus,
\[ \B_0 = \frac{\omega}{c^2}\sum_{k=0}^\infty \frac{f_{k-1}(z)}{k+1} [r^{k-1}y,-r^{k-1}x,0]. \]

The full time-dependent field can be recovered with $\E=\E_0\cos\omega t$ and $\B=\B_0\sin\omega t$.

\section{Case with Angular Dependency}
Defining cylindrical coordinates with $x=r\cos\theta$ and $y=r\sin\theta$, consider a longitudinal ($z$ directed) vector field $\mathbf a$ with $n$-fold rotational symmetry in $\theta$ and $r^k$ scaling in the radial direction:
\[ a_z=e^{in\theta}r^kf(z), \quad a_r=a_\theta=0. \]
Here, the field components are complex so the cosine or sine part can be extracted by taking the real or imaginary component respectively.  These components still satisfy all the differential equations.  The curl in cylindrical polar coordinates $(r,\theta,z)$ has two nonzero components:
\begin{eqnarray*}
(\curl\mathbf a)_r &=& \tfrac1r\partial_\theta a_z = ine^{in\theta}r^{k-1}f(z), \\
(\curl\mathbf a)_\theta &=& -\partial_ra_z = -e^{in\theta}kr^{k-1}f(z).
\end{eqnarray*}
It is sufficient to consider two other single-component terms $\mathbf b$ and $\mathbf c$ where
\[ b_\theta=e^{in\theta}r^kf(z), \quad b_r=b_z=0, \]
\[ c_r=e^{in\theta}r^kf(z), \quad c_\theta=c_z=0. \]
The curl of $\mathbf b$ has nonzero components
\begin{eqnarray*}
(\curl\mathbf b)_r &=& -\partial_z b_\theta = -e^{in\theta}r^kf'(z), \\
(\curl\mathbf b)_z &=& \tfrac1r\partial_r(rb_\theta) = e^{in\theta}(k+1)r^{k-1}f(z).
\end{eqnarray*}
The curl of $\mathbf c$ has nonzero components
\begin{eqnarray*}
(\curl\mathbf c)_\theta &=& \partial_z c_r = e^{in\theta}r^kf'(z), \\
(\curl\mathbf c)_z &=& -\tfrac1r\partial_\theta c_r = -ine^{in\theta}r^{k-1}f(z).
\end{eqnarray*}
As a shorthand, subscripts will indicate substitutions of $k$ or $f$ in these fields, omitting any that are unchanged.  In this notation,
\begin{eqnarray*}
\curl\mathbf a &=& in\mathbf c_{k-1} - k\mathbf b_{k-1}, \\
\curl\mathbf b &=& -\mathbf c_{f'} + (k+1)\mathbf a_{k-1}, \\
\curl\mathbf c &=& \mathbf b_{f'} - in\mathbf a_{k-1}.
\end{eqnarray*}
Now the double curls may be evaluated:
\begin{eqnarray*}
\curl(\curl\mathbf a) &=& in(\mathbf b_{k-1,f'} - in\mathbf a_{k-2}) - k( -\mathbf c_{k-1,f'} + k\mathbf a_{k-2}) \\
&=& (n^2-k^2)\mathbf a_{k-2} + in\mathbf b_{k-1,f'} + k\mathbf c_{k-1,f'}, \\
\curl(\curl\mathbf b) &=& -(\mathbf b_{f''} - in\mathbf a_{k-1,f'}) + (k+1)(in\mathbf c_{k-2} - (k-1)\mathbf b_{k-2}) \\
&=& in\mathbf a_{k-1,f'} - \mathbf b_{f''} - (k^2-1)\mathbf b_{k-2} + in(k+1)\mathbf c_{k-2}, \\
\curl(\curl\mathbf c) &=& (-\mathbf c_{f''} + (k+1)\mathbf a_{k-1,f'}) - in(in\mathbf c_{k-2} - (k-1)\mathbf b_{k-2}) \\
&=& (k+1)\mathbf a_{k-1,f'} + in(k-1)\mathbf b_{k-2} - \mathbf c_{f''} + n^2\mathbf c_{k-2}.
\end{eqnarray*}
Seek a solution of the form
\[ \E_0 = \sum_{k=0}^\infty \mathbf a_{k,f_k} + \mathbf b_{k,g_k} + \mathbf c_{k,h_k} \]
for some functions $f_k,g_k,h_k$.  The double curl is
\begin{eqnarray*}
\curl(\curl\E_0) &=& \sum_{k=0}^\infty 
(n^2-k^2)\mathbf a_{k-2,f_k} + in\mathbf b_{k-1,f'_k} + k\mathbf c_{k-1,f'_k} \\
&+& in\mathbf a_{k-1,g'_k} - \mathbf b_{g''_k} - (k^2-1)\mathbf b_{k-2,g_k} + in(k+1)\mathbf c_{k-2,g_k} \\
&+& (k+1)\mathbf a_{k-1,h'_k} + in(k-1)\mathbf b_{k-2,h_k} - \mathbf c_{h''_k} + n^2\mathbf c_{k-2,h_k}.
\end{eqnarray*}
Equating this to $\frac{\omega^2}{c^2}\E_0$ and only taking the $r^k$ terms gives
\begin{eqnarray*}
 \frac{\omega^2}{c^2}( \mathbf a_{f_k} + \mathbf b_{g_k} + \mathbf c_{h_k} ) 
&=&  (n^2-(k+2)^2)\mathbf a_{f_{k+2}} + in\mathbf b_{f'_{k+1}} + (k+1)\mathbf c_{f'_{k+1}} \\
&+& in\mathbf a_{g'_{k+1}} - \mathbf b_{g''_k} - ((k+2)^2-1)\mathbf b_{g_{k+2}} + in(k+3)\mathbf c_{g_{k+2}} \\
&+& (k+2)\mathbf a_{h'_{k+1}} + in(k+1)\mathbf b_{h_{k+2}} - \mathbf c_{h''_k} + n^2\mathbf c_{h_{k+2}}.
\end{eqnarray*}
Equating the longitudinal, azimuthal and radial vector coefficients gives
\begin{eqnarray*}
\frac{\omega^2}{c^2}f_k &=& (n^2-(k+2)^2)f_{k+2} + ing'_{k+1} + (k+2)h'_{k+1}, \\
\frac{\omega^2}{c^2}g_k &=& inf'_{k+1} - g''_k - ((k+2)^2-1)g_{k+2} + in(k+1)h_{k+2}, \\ 
\frac{\omega^2}{c^2}h_k &=& (k+1)f'_{k+1} + in(k+3)g_{k+2} - h''_k + n^2h_{k+2}.
\end{eqnarray*}
The highest function in the first equation is $f_{k+2}$, which is good for setting up a recurrence.  However, this time the $g$ and $h$ equations have a combination of $g_{k+2}$ and $h_{k+2}$ in each of them and the two combinations are linearly dependent (note $(k+2)^2-1=(k+1)(k+3)$), so there is not enough information to solve for $g_{k+2}$ and $h_{k+2}$ separately.  On the other hand, the linear dependence allows one equation to be obtained where these terms are removed entirely.

To simplify manipulation, define the mixtures
\[ m_k=ing_k + (k+1)h_k \qquad \mbox{and} \qquad p_k=(k+1)g_k - inh_k, \]
where $p$ is `perpendicular' to $m$ in $g$--$h$ space.  Now the three equations can be written
\begin{eqnarray*}
\frac{\omega^2}{c^2}f_k &=& (n^2-(k+2)^2)f_{k+2} + m'_{k+1}, \\
\frac{\omega^2}{c^2}g_k &=& inf'_{k+1} - g''_k - (k+1)p_{k+2}, \\ 
\frac{\omega^2}{c^2}h_k &=& (k+1)f'_{k+1} - h''_k + inp_{k+2}.
\end{eqnarray*}
Further mixing the second and third equations yields
\begin{eqnarray*}
\frac{\omega^2}{c^2}m_k &=& ((k+1)^2-n^2)f'_{k+1} - m''_k, \\ 
\frac{\omega^2}{c^2}p_k &=& - p''_k + (n^2-(k+1)^2)p_{k+2}.
\end{eqnarray*}
The $p$ equation is now completely decoupled.  Taking the $z$ derivative of the $f$ equation for the $k-1$ case (a trick also used in the no angular dependency case) gives
\[ \frac{\omega^2}{c^2}f'_{k-1} = (n^2-(k+1)^2)f'_{k+1} + m''_k, \]
which equals $-1$ times the $m$ equation, so
\[ m_k=-f'_{k-1} \qquad \Rightarrow \qquad m'_{k+1}=-f''_k. \]
This allows the $f$ equation to be decoupled:
\[ \frac{\omega^2}{c^2}f_k = (n^2-(k+2)^2)f_{k+2} - f''_k. \]
A similar decoupling is possible for the $m$ equation but then the relation $m_k=-f'_{k-1}$ would be lost.  Now everything can be written as simple recurrence formulae
\begin{eqnarray*}
f_{k+2} &=& \frac{1}{n^2-(k+2)^2}\left(\frac{\omega^2}{c^2}f_k + f''_k\right), \\
p_{k+2} &=& \frac{1}{n^2-(k+1)^2}\left(\frac{\omega^2}{c^2}p_k + p''_k\right), \\
m_{k+1} &=& -f'_k.
\end{eqnarray*}
$g$ and $h$ may be recovered with
\[ g_k=\frac{inm_k + (k+1)p_k}{(k+1)^2-n^2} \qquad \mbox{and} \qquad  h_k=\frac{(k+1)m_k - inp_k}{(k+1)^2-n^2}. \]
To keep the polynomial nature of the functions near the axis, $k\ge n$ in $r^ke^{in\theta}$, so the lowest allowed nonzero functions are $f_n$ for the on-axis electric field and $g_{n+1}$ and $h_{n+1}$.

\section{Other Approaches}
It is also possible to express $\curl(\curl\E_0) = \nabla(\divv\E_0) - \nabla^2\E_0$ where $\nabla^2\E_0$ is the vector Laplacian that must be chosen appropriately for the coordinate system being used.  This has a simpler form than the double curl, so assuming that $\divv\E_0=0$ by one of Maxwell's equations, it gets to the second-order recurrence relations more quickly.  The expression $\nabla^2\E_0=-\frac{\omega^2}{c^2}\E_0$ is known as the Helmholtz equation.  The downside of this approach is that the first-order recurrence relation is not recovered, because it comes from the assumption that $\divv\E_0=0$, so it must be added manually by substituting the series for $\E_0$ into that equation.

\section{Implementation Notes}
It is possible to set $\omega=0$ in the above formulae to get an entirely electrostatic solution.

These formulae work well when the repeated derivatives of the on-axis function(s) such as $f^{(n)}$ are easy to calculate.  One example is when the entrance and exit ramps of a field area are constructed from $\tanh(z/l)$ for some fringe field length $l$.  As $\dd{}{z}\tanh z=1-\tanh^2 z$, the $n^\mathrm{th}$ derivative is an $(n+1)^\mathrm{th}$ order polynomial in $\tanh z$.  For a field of length $L$, the exit ramp looks like $\tanh((L-z)/l)$, which has similar repeated derivatives except the odd order ones have changed sign.  This corresponds to changing the sign of the $h$ terms in the no angular dependency case.

Once an array of repeated derivatives $[f,f',f''\cdots f^{(N)}]$ has been calculated, the core operation is applying the operator $\frac{\omega^2}{c^2}+\pdds{}{z}$ to it repeatedly.  This is easy, since scalar multiplication works as usual and differentiation shifts the array to the left.  Additionally, the second entry remains the derivative of the first, which is useful since both $f_{2k}$ and $f'_{2k}$ are required after each application of the second-order operator.


\begin{thebibliography}{9}

	\bibitem{BrooksOffAxis}
	S.~J.~Brooks,
	\textquotedblleft{Off-axis magnetic fields extrapolated from on-axis multipoles}\textquotedblright,
	BNL, Upton, NY, USA, Rep.~BNL-223622-2022-TECH, 2013.
	\doi{10.2172/1895080}

\end{thebibliography}

\end{document}