\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{\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{\floor}[1]{\left\lfloor#1\right\rfloor}
\newcommand{\mod}[1]{\left|#1\right|}

\newcommand{\rem}[1]{}

\newcommand{\B}{\mathbf B}
\newcommand{\J}{\mathbf J}
\newcommand{\x}{\mathbf x}

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

\begin{document}

\title{Magnetic Field of a Winding Sheet defined by a Contour Function}
\author{Stephen Brooks}
\maketitle

\section{Problem}
The source of the magnetic field is a flat winding sheet in the $z=0$ plane with current vector $\J$ parallel to the plane:
\[ \J(x,y,z) = \delta(z)\vecthree{j_x(x,y)}{j_y(x,y)}{0}. \]
The magnetic field $\B$ is governed by Maxwell's equations with a current source but no other materials ($\mu=\mu_0$):
\[ \nabla\cdot\B=0, \qquad \nabla\times\B=\mu_0\J. \]

\section{Current Contour Function}
Conservation of current requires that $\nabla\cdot\J=0$.  One way of enforcing this for this problem is to make the (2D) current flow around the contour lines of a function $\phi(x,y)$ with density equal to the density of the contours.  In other words, the current is $\nabla_{x,y}\phi$ rotated by 90 degrees:
\[ j_x=-\partial_y\phi, \qquad j_y=\partial_x\phi. \]

\subsection{Units}
$\J$ has units of A/m$^2$, $j_x$ and $j_y$ have units of A/m and $\phi$ has units of Amps.

\section{Field Evolution in $z$}
Rearranging Maxwell's equations to bring the $z$ derivatives all on one side gives:
\[ \partial_z\B = \vecthree{\partial_xB_z+\mu_0J_y}{\partial_yB_z-\mu_0J_x}{-\partial_xB_x-\partial_yB_y}, \]
with the fourth equation becoming a consistency condtion within the plane $\partial_xB_y-\partial_yB_x=\mu_0J_z$, which equals zero in this problem.

\subsection{Source Discontinuity}
The delta function means the 3D current density $\J$ becomes `infinite' around $z=0$.  The equations can be made proper by integrating over a small region around zero:
\[ \int_{0^-}^{0^+} \partial_z\B \de z = \int_{0^-}^{0^+} \vecthree{\partial_xB_z+\mu_0J_y}{\partial_yB_z-\mu_0J_x}{-\partial_xB_x-\partial_yB_y} \de z. \]
Here, the left side is $\B(x,y,0^+)-\B(x,y,0^-)$, the difference in magnetic field across the discontinuity, which shall be written $\Delta\B(x,y)$ from now on.  On the right side, any finite-valued functions such as $\partial_xB_z$ within the integral will integrate to zero as the region is small; the terms including $\J$ will give finite results as $\int_{0^-}^{0^+}\delta(z)\de z=1$.  Thus,
\[ \Delta\B(x,y) = \mu_0\vecthree{j_y}{-j_x}{0} = \mu_0\vecthree{\partial_x\phi}{\partial_y\phi}{0}. \]

\subsection{Solution for $z\ne0$}
Consider a single Fourier mode of the contour function $\phi(x,y)=\sin(ax)\sin(by)$.  Any other function can be built from these modes by linear superposition (using different $a$ and $b$) and translation in $x$ and $y$.  The discontinuity is
\[ \Delta\B(x,y) = \mu_0\vecthree{a\cos(ax)\sin(by)}{b\sin(ax)\cos(by)}{0}. \]
Without current sources, the free-space evolution in $z$ is a differential operator:
\[ \partial_z^\mathrm{FS}\B = \vecthree{\partial_xB_z}{\partial_yB_z}{-\partial_xB_x-\partial_yB_y}. \]
Applying this repeatedly to the two-dimensional form of our $\B$ field gives terms like:
\[ \partial_z^\mathrm{FS}\vecthree{a\cos(ax)\sin(by)}{b\sin(ax)\cos(by)}{0} = \vecthree{0}{0}{(a^2+b^2)\sin(ax)\sin(by)} \]
and
\[ \partial_z^\mathrm{FS}\vecthree{0}{0}{\sin(ax)\sin(by)} = \vecthree{a\cos(ax)\sin(by)}{b\sin(ax)\cos(by)}{0}. \]
Consider a 3D magnetic field with the following form:
\[ \B(x,y,z) = f(z)\vecthree{a\cos(ax)\sin(by)}{b\sin(ax)\cos(by)}{0} + g(z)\vecthree{0}{0}{\sin(ax)\sin(by)}. \]
Maxwell's equations in free space will be satsified if $\partial_z\B=\partial_z^\mathrm{FS}\B$,
\[ \partial_z^\mathrm{FS}\B = f(z)(a^2+b^2)\vecthree{0}{0}{\sin(ax)\sin(by)} + g(z)\vecthree{a\cos(ax)\sin(by)}{b\sin(ax)\cos(by)}{0}, \]
so equating coefficients gives: 
\[ f'(z)=g(z), \qquad g'(z)=(a^2+b^2)f(z). \]
The general solution is then $f(z)=Ae^{\sqrt{a^2+b^2}z}+Be^{-\sqrt{a^2+b^2}z}$.

\subsection{Whole Solution}
\newcommand{\sgn}{\mathrm{sgn}}
To avoid fields that increase to infinity, parts of the general solution can be chosen either side of the $z=0$ plane as follows:
\[ f(z<0)=Ae^{\sqrt{a^2+b^2}z}, \qquad f(z>0)=Be^{-\sqrt{a^2+b^2}z}. \]
At $z=0$, the discontinuity also fits the same general form of $\B$ provided that:
\[ f(0^+)-f(0^-)=\mu_0, \qquad g(0^+)=g(0^-). \]
The first equation implies that $B-A=\mu_0$.  The second implies that $g=f'$ is continuous, which happens if $A=-B$.  Therefore
\[ f(z<0)=-\frac{\mu_0}{2}e^{\sqrt{a^2+b^2}z}, \qquad f(z>0)=\frac{\mu_0}{2}e^{-\sqrt{a^2+b^2}z}; \]
\[ g(z)=-\frac{\mu_0}{2}\sqrt{a^2+b^2}e^{-\sqrt{a^2+b^2}|z|}. \]
Putting this all together, the solution can be written
\[ \B(x,y,z) = \frac{\mu_0}{2} \vecthree{\sgn(z)a\cos(ax)\sin(by)}{\sgn(z)b\sin(ax)\cos(by)}{-\sqrt{a^2+b^2}\sin(ax)\sin(by)} e^{-\sqrt{a^2+b^2}|z|}, \]
where $\sgn(z)=-1$ for $z<0$ and $1$ for $z>0$.

\section{Application: Double Layer}
The field of a single winding sheet may not have many applications but interesting controlled fields can be formed between a pair of sheets.  If two sheets with identical parallel currents are stacked above each other in $z$, the midplane field has only a $B_z$ component, as the $\sgn(z)$ terms in $B_x$ and $B_y$ cancel:
\[ B_{z,\mathrm{parallel}}(x,y,0) = -\mu_0\sqrt{a^2+b^2}\sin(ax)\sin(by)e^{-\sqrt{a^2+b^2}Z}, \]
where the sheets at $z=\pm Z$ have $\phi_{\pm Z}(x,y)=\sin(ax)\sin(by)$.

If the otherwise-identical sheets are given opposing current directions, the midplane field has $B_z=0$ and
\[ \B_\mathrm{opposing}(x,y,0) = \mu_0 \vecthree{-a\cos(ax)\sin(by)}{-b\sin(ax)\cos(by)}{0} e^{-\sqrt{a^2+b^2}Z}, \]
where the sheet at $z=Z$ has $\phi_Z(x,y)=\sin(ax)\sin(by)$ and $\phi_{-Z}(x,y)=-\sin(ax)\sin(by)$.

\subsection{Small Separation Approximation}
If the layer separation $Z$ is much smaller than the wavelengths $2\pi/a$, $2\pi/b$, then the exponent $\sqrt{a^2+b^2}Z\simeq0$.  Going down to two dimensions ($x,y$), the midplane fields can be written
\[ B_{z,\mathrm{parallel}}(x,y) \simeq -\mu_0\sqrt{a^2+b^2}\phi_Z, \]
\[ \B_\mathrm{opposing}(x,y) \simeq \mu_0 \vectwo{-a\cos(ax)\sin(by)}{-b\sin(ax)\cos(by)} = -\mu_0\nabla_{x,y}\phi_Z. \]
The $\B_\mathrm{opposing}$ field is approximately the (2D) gradient of $\phi$ (this is now true for any general function $\phi$ that is a superposition of Fourier modes).  Given a goal $\B$ field, such a potential $\phi$ can always be found because the in-plane consistency condition $\partial_xB_y-\partial_yB_x=0$ from Maxwell's equations guarantees the field has no curl in this plane.

\subsection{Iteration for Opposing Current Sheets}
Suppose $\B_\mathrm{goal}(x,y)$ is known and the sheets have opposing currents.  Using the small separation approximation in the previous section, make a first guess at the required winding potential
\[ \phi_{Z,1}(x,y) = \frac{-1}{\mu_0}\int_{-\infty}^{y}B_{\mathrm{goal},y}(x,\hat{y}) \de\hat{y}. \]
Work out the true field $\B_1(x,y)$ from this potential without using any approximations.  Then define new winding potentials via
\[ \phi_{Z,n+1}(x,y) = \phi_{Z,n}(x,y)-\frac{1}{\mu_0}\int_{-\infty}^{y}B_{\mathrm{goal},y}(x,\hat{y})-B_{n,y}(x,\hat{y}) \de\hat{y}. \]
This sort of iteration may be able to correct for other small differences to the idealised case, such as curvature in a sector magnet.

\subsection{Direct Evaluation via Fourier Transform}
The field formula for sheets with parallel currents can be inverted for each Fourier mode:
\[  B_{z,\mathrm{parallel}} = -\mu_0\sqrt{a^2+b^2}e^{-\sqrt{a^2+b^2}Z}\phi_{\pm Z} 
\qquad \Rightarrow \qquad
\phi_{\pm Z} = \frac{-1}{\mu_0}\frac{e^{\sqrt{a^2+b^2}Z}}{\sqrt{a^2+b^2}}B_{z,\mathrm{parallel}}. \]
If a given $B_{z,\mathrm{goal}}$ can be expressed as a sum of Fourier modes, the required winding potential can be obtained as follows:
\begin{eqnarray*}
B_{z,\mathrm{goal}} &=& \sum_{a,b,\alpha,\beta} c_{ab\alpha\beta}\sin(ax+\alpha)\sin(by+\beta) \\
\qquad \Rightarrow \qquad
\phi_{\pm Z} &=& \frac{-1}{\mu_0} \sum_{a,b,\alpha,\beta} c_{ab\alpha\beta} \frac{e^{\sqrt{a^2+b^2}Z}}{\sqrt{a^2+b^2}}\sin(ax+\alpha)\sin(by+\beta).
\end{eqnarray*}

For sheets with opposite currents, if the $\B_\mathrm{goal}$ field parallel to the midplane is achievable in free space, there must be some ideal potential $\Phi$ (for the midplane field rather than the windings) such that $\B_\mathrm{goal} = -\mu_0\nabla_{x,y}\Phi$.  Going back to a single Fourier mode $\phi_Z(x,y)=\sin(ax)\sin(by),$ 
\[ \B_\mathrm{opposing} = \mu_0 \vectwo{-a\cos(ax)\sin(by)}{-b\sin(ax)\cos(by)} e^{-\sqrt{a^2+b^2}Z}
 = -\mu_0\nabla_{x,y}\phi_Z e^{-\sqrt{a^2+b^2}Z} , \]
so that if $\phi_Z=\Phi e^{\sqrt{a^2+b^2}Z}$, then $\B_\mathrm{opposing} = \B_\mathrm{goal}$ for the amplitude of this particular mode.  This defines a `filter' in frequency space where the winding potential $\phi_Z$ can be obtained by multiplying each frequency component of the ideal midplane potential $\Phi$ by $e^{\sqrt{a^2+b^2}Z}\ge1$.

In theory, this will give the perfect field in one step, but higher frequency components or large separations $Z$ will cause the amplitude multiplication factor to become exponentially large, leading to an impractical winding configuration.  In particular, discontinuities in $\B_\mathrm{goal}$ or its derivatives will produce a spectrum with high frequency tails, which may not even converge when amplified by an exponential of frequency.  One way to ensure good convergence is to apply a Gaussian blur to the required field, which will multiply frequency tails by a function like $e^{-kf^2}$.  Another way is to explicitly construct the goal field as a Fourier series only containing lower frequencies.

\subsection{Opposing Current Sheets between Infinite Parallel Iron Plates}
The frequency filter that takes $\phi_Z$ to $\Phi$ (windings to midplane field) has an impulse response function proportional to $Z/(Z^2+d^2)^{3/2}$ where $d=\sqrt{x^2+y^2}$ is transverse distance.  This means a local change to the winding potential has an effect with a long $d^{-3}$ tail in the midplane potential, which is a $d^{-4}$ tail in the field.  Magnets designed using the Fourier method above often have many weak windings at large distances from the main magnet to cancel this tail, which is normally not a part of the desired fringe field functions.

If the current sources are placed between two parallel (perfect $\mu=\infty$) iron plates, the natural fringe decay rate becomes exponential, which is a better fit for most designs.  If the plates are at $z=\pm2Z$, the image currents will form a regular pattern, with parallel-current images at $\pm3Z$, images of the opposing-current side at $\pm5Z$, then second images of that at $\pm7Z$ and continuing to infinity in the pattern $++--++--$ including the $\pm Z$ original winding at the start.  This means the formula for getting the midplane potential from a Fourier mode of a single winding is:

\[ \Phi_\mathrm{plates} = \phi_Z \sum_{n=0}^\infty (-1)^n\left(e^{-\sqrt{a^2+b^2}(4n+1)Z}+e^{-\sqrt{a^2+b^2}(4n+3)Z}\right). \]

Letting $u=e^{-\sqrt{a^2+b^2}Z}$, the sum is $\sum_{n=0}^\infty(-1)^n(u^{4n+1}+u^{4n+3}) = (u+u^3)\sum_{n=0}^\infty (-u^4)^n = \frac{u+u^3}{1+u^4} = \frac{u^{-1}+u}{u^{-2}+u^2}$.  Thus,

\[ \Phi_\mathrm{plates} = \phi_Z \frac{e^{\sqrt{a^2+b^2}Z}+e^{-\sqrt{a^2+b^2}Z}}{e^{2\sqrt{a^2+b^2}Z}+e^{-2\sqrt{a^2+b^2}Z}}
 = \phi_Z \frac{\cosh(\sqrt{a^2+b^2}Z)}{\cosh(2\sqrt{a^2+b^2}Z)}. \]

%An interesting piece of trivia says that sech is its own Fourier transform, so the impulse response for winding changes now looks like sech as well and has exponential tails.

\subsection{Asymmetric Sheet Pairs for General Fields}
If $\B_\mathrm{goal}$ has all three components, a superposition of parallel and opposing current sheets can generate the field.  If parallel current sheets with potentials $\phi_Z^\mathrm{p}=\phi_{-Z}^\mathrm{p}$ generate only the $B_{z,\mathrm{goal}}$ component and opposing current sheets with potentials $\phi_Z^\mathrm{o}$ and $\phi_{-Z}^\mathrm{o}=-\phi_Z^\mathrm{o}$ generate the $\B_{x,y,\mathrm{goal}}$ components parallel to the midplane, then the superimposed winding potentials $\phi_Z=\phi_Z^\mathrm{p}+\phi_Z^\mathrm{o}$ and  $\phi_{-Z}=\phi_Z^\mathrm{p}-\phi_Z^\mathrm{o}$ will generate the complete field.  In general, these two windings will not be symmetrical.

%\subsection{Correction of General Field Errors}

\end{document}