##### Solenoidal Fields for Ion Beam Transport and Focusing {#solenoidal-fields-for-ion-beam-transport-and-focusing .unnumbered}

Edward Lee and Matthaeus Leitner

Lawrence Berkeley National Laboratory

November 2007

**1. [Introduction and Contents]{.underline}**

In this report we calculate time-independent fields of solenoidal
magnets that are suitable for ion beam transport and focusing.

Electricity and Magnetism Textbooks present the formalism for magnetic
field calculations and apply it to simple geometries \[1-1\], but they
do not include enough relevant detail to be very useful for designing a
charged particle transport system. This requires accurate estimates of
fringe field aberrations, misaligned and tilted fields, peak fields in
coils and iron, external fields, and more. Specialized books on magnet
design, technology, and numerical computations \[1-2\] contain relevant
detail of this type, and some of that is presented here. The AIP
Conference Proceedings of the US Particle Accelerator Schools \[1-3\]
contain extensive discussions of design and technology of magnets for
ion beams -- except for solenoids. This lack may be due to the fact that
solenoids have been used primarily to transport and focus particles of
relatively low momentum, e.g. electrons of less than 50 MeV and protons
or H^-^ of less than 1.0 MeV, although this situation may be changing
with the commercial availability of superconducting solenoids with up to
20T bore field \[1-4\]. Internal reports from federal laboratories and
industry treat solenoid design in detail for specific applications. The
present report is intended to be a resource for the design of ion beam
drivers for Inertial Fusion Energy \[1-5\] and Warm Dense Matter
experiments \[1-6\], although it should also be useful for a wider range
of applications.

The field produced by specified currents and magnetized materials can
always be evaluated by solving Maxwell's equations numerically, but it
is also desirable to have reasonably accurate, simple formulas for
system design and fast-running beam dynamics codes, as well as for
general understanding. Most of this report is devoted to such formulas,
but an introduction to the Tosca code \[1-7\] and some numerical results
obtained with it are also presented. Details of design, fabrication,
installation, and operation of real magnet systems are not included.
Mathematical derivations are presented with a moderate amount of detail.
While there is no claim of originality, except for various numerical
approximations and an induction module design in section 20, many of the
results and discussions are hard to find elsewhere.

Our primary topic is axisymmetric solenoidal systems with no magnetic
materials. These simplifying features allow useful analytical
calculations, which occupy sections 2-13. Deviations from axisymmetry
are considered in sections 14, 15, 21, and 22, and the effects of
magnetic materials are treated in sections 16-20. Since magnetic
aberrations are mixed with geometric aberrations in computing ion
orbits, section 22 on the ion equations of motion is included.

**Contents**

1.  [Introduction and Contents]{.underline}

2.  [Solenoid Basics]{.underline}

    1.  Model assumptions and equations

    2.  Simple field estimate

    3.  Current layers, global picture, and moments

    4.  Vector potential and flux lines

    5.  Stress and energy

3.  [Near-Axis Field]{.underline}

    1.  On-axis field and near-axis expansion

    2.  Useful integrals

4.  [Design Formulas for the On-Axis Field]{.underline}

    1.  Simple cases

    2.  Thick wire layer

    3.  Uniform fields

5.  [Vector and Scalar Potentials]{.underline}

    1.  Vector potential

    2.  Scalar potential

    3.  Global scalar potential

6.  [Green Function Solutions]{.underline}

    1.  Green function formalism

    2.  Explicit field expressions

7.  [Approximation of Green Function Integrals]{.underline}

8.  [Thin Wire Layer]{.underline}

    1.  Semi-infinite layer

    2.  Wire end divergence

    3.  Finite length lens

    4.  Stored energy

9.  [Thick Wire Layer]{.underline}

    1.  Uniform thick layer

    2.  Approximation by a thin layer

    3.  General thick layer

10. [Periodic Thin Wire Layer]{.underline}

11. [Far Field and External Field]{.underline}

> 11.1 Far field expansion
>
> 11.2 External field expansion

12. [Field in the Wire Layer]{.underline}

13. [Current to Produce a Given Field]{.underline}

14. [Misaligned and Tilted Solenoid]{.underline}

15. [Multiple Channels -- Cross Talk]{.underline}

16. [Magnetic Materials Basics]{.underline}

> 16.1 Ferromagnetism
>
> 16.2 Magnetization formalism
>
> 16.3 B-H relations and jump conditions
>
> 16.4 Magnetization curve
>
> 16.5 Permanent magnet solenoid

17. [Wire Layer Embedded in a Highly Permeable Yoke]{.underline}

18. [Periodic Thin Wire Layer with a Highly Permeable Yoke]{.underline}

19. [Tosca Model]{.underline}

> 19.1 General solution strategy in 3D
>
> 19.2 Solution strategy in 2D
>
> 19.3 Example

20. [Application of Tosca to Induction Linac Design]{.underline}

21. [Biot-Savart Formula]{.underline}

> 21.1 General formula
>
> 21.2 Straight wire segment field

21.3 Helix field

22. [Single Particle Equations of Motion in an Arbitrary
    field]{.underline}

> 22.1 Time as the independent variable
>
> 22.2 Longitudinal distance as the independent variable
>
> 22.3 Paraxial approximation
>
> **References**

1.  Jackson, J.D., Classical Electrodynamics. 2^nd^ or 3^rd^ Edition,
    Wiley, New York, (1975 or 1999), chap.5.

> Feyman, R., Leighton, R., and Sands, M., The Feyman Lectures on
> Physics, Vol. II, Chs. 13, 14, 15, 36, 37, Addison -- Wesley, 1962.
>
> Panofsky, W.K.H., and Phillips, M., Classical Electricity and
> Magnetism, 2^nd^
>
> edition, Addison-Wesley, Reading, Mass. (1962), chaps. 7 & 8.
>
> Smythe, W.R., Static and Dynamical Electricity, 3^rd^ edition,
> McGraw-Hill, New
>
> York (1969), chap. 7.

2.  Wilson, M.N., Superconducting Magnets, Clarendon Press, Oxford
    (1983).

1-3 Magnet calculations are presented in many volumes of the AIP
Proceedings for

The US Particle Accelerator Schools (http://uspas.fnal.gov); one example
is:

> Willen, E., Dahl, P. and Herrera, J., Superconducting Magnets, AIP
> Conference Proceedings 153, pg 1228, American Institute of Physics,
> New York (1987).
>
> 1-4 For example: American Magnetics, Inc., 112 Flint Rd., Oakridge, TN
> 37831,

[www.americanmagnetics.com](http://www.americanmagnetics.com/)

5.  Lee, E., Solenoidal Transport for Heavy Ion Fusion, Nuclear
    Instruments and Methods in Physics Research, Volume 544, Issues 1-2,
    2005, 187.

6.  Barnard, J.J., et al., Accelerator and Ion Beam Tradeoffs for
    Studies of Warm Dense Matter, Proceedings of the 2005 Particle
    Accelerator Conference, Knoxville TN, 05/16-20/2005.

7.  Tosca is commercially available. Contact: Vector Fields Inc.,
    1700 N. Farnsworth Avenue, Aurora, IL 60505, USA,
    [http://www.vectorfields.com](http://www.vectorfields.com/).

> Another code widely used for magnet computations is Poisson Superfish,
> which is available from the Los Alamos Accelerator code group:
>
> Halback, K. and Holsinger, R.F. "SUPERFISH-A Computer Program for
> Evaluation of RF Cavities with Cylindrical Symmetry", Particle
> Accelerators, [7]{.underline} 213 (1976), and "Poisson Superfish
> Documentation" by James H. Billen and Lloyd M. Young, Los Alamos
> National Laboratory Report \# LA-UR-96-1834. The software can be
> downloaded at http://laacg1.lanl.gov.
>
> **Acknowledgements**
>
> The authors thank Ms. Lynn Heimbucher for her great effort and
> patience over several years in shepherding this report through its
> many revisions. Helpful discussions were held with Shlomo Caspi, Andy
> Faltens, Alex Friedman, Enrique Henestroza, Steven Lund, and William
> Waldron.

This work was supported by the Director, Office of Science, and Office
of Fusion Energy Sciences, of the U.S. Department of Energy under
Contract No. DE-AC02-05CH11231.

2.  **[Solenoid Basics]{.underline}**

> **2.1 Model assumptions and equations**

We wish to calculate the field of a system of solenoidal magnets
suitable for the transport and focusing of ion beams. Generally such a
field will be very strong, in the range 2-20 T in the magnet bore and
current-carrying wire. In this report we treat fields that are constant
in time and have strict solenoidal symmetry around the z-axis (except in
sections 14, 15, 21 and 22). Cylindrical coordinates (r,$\theta,$z) are
used for the general position vector $\overrightarrow{r}$ unless
indicated otherwise. The only component of wire current density is then
J${}_{\theta}$(r,z), which is the source of the two field components
$B_{r}$ (r,z) and $B_{z}$(r,z). This is true even when magnetic
materials are present since they have only a current density
$J_{\theta}^{\text{mag}}$(r,z) equal to the curl of the magnetization
density ($M_{r}$(r,z), $M_{z}$(r,z)). Axisymmetric induction cores, with
idealized current and field components
$\left( J_{r},J_{z},B_{q},M_{q} \right)$, are therefore not treated
here.

In general the magnetization density is not independent of wire current
except for the case of stiff permanent magnet material; this feature
(permeability) greatly increases the complexity of the analysis and
makes it necessary to use codes such as Poisson and Tosca \[1-7\] in all
but the simplest situations. Consideration of magnetic material
properties is delayed until section 16. In the absence of permeable
magnetic materials, solenoid fields may be treated analytically or with
elementary codes.

Wire leads and windings are assumed to be paired such that field
aberrations due to deviations from solenoidal symmetry are cancelled
down to a negligible level. This includes cancellation of the small z
component of current that is present in a single (helical) coil of wire
by returning the current in a second coil on top of it. Also, the
small-scale variations of field inside and close to individual wire or
cable sections are ignored, i.e. a smoothed model for $J_{\theta}(r,z)$
is adopted. In effect, the dominant, collective field from the current
density is calculated here, without any distinction between wire and
cable, and the small-scale variations are left to detailed magnet
computations. However, one must be careful when representing the actual
currents by a smoothed $J_{q}$(r,z). For example the effective length
and dipole moment of a magnet (see eqns. 2.11-2.14) should be
incorporated into the smoothed model.

SI units (Tesla, Ampere, meter) are used throughout this report. The
relevant Maxwell's equations in the absence of magnetic materials are

> $\overrightarrow{\nabla}x\overrightarrow{B} = \mu_{0}\overrightarrow{J},\begin{matrix}
>  & 
> \end{matrix}\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$,
> (2.1a,b)

where $m_{0} = 4px\text{10}^{- 7}T - m/A$. To convert any equation to
Gaussian (cgs) form, substitute $\mu_{0} \rightarrow 4\pi/c$. Unit
conversion factors are 1.0 Tesla=$\text{10}^{4}$Gauss, 1.0 Ampere =
2.99792...*x*$\text{10}^{9}$ statampere, and 1.0m=$\text{10}^{2}$cm.

Since the variable r denotes distance from the z-axis,
$\left| \overrightarrow{r} \right| = \sqrt{r^{2} + z^{2}}$ is the
distance from the coordinate origin. Assuming solenoidal symmetry, eqns.
(2.1 a,b) become

> $\frac{¶B_{r}}{¶z} - \frac{¶B_{z}}{¶r} = m_{0}J_{q}\begin{matrix}
> , & 
> \end{matrix}\frac{1}{r}\frac{¶}{¶r}\text{rB}_{r} + \frac{¶B_{z}}{¶z} = 0$.
> (2.2a,b)

From these coupled first-order equations one can derive the decoupled
second-order equations:

> $\frac{1}{r}\frac{¶}{¶r}r\frac{¶B_{z}}{¶r} + \frac{¶^{2}B_{z}}{¶z^{2}}\text{=-}m_{0}\frac{1}{r}\frac{¶}{¶r}\text{rJ}_{q}$,
> (2.3a)
>
> $$$\frac{\partial}{\partial r}\frac{1}{r}\frac{\partial}{\partial r}\text{rB}_{r} + \frac{\partial^{2}B_{r}}{\partial z^{2}} = \mu_{0}\frac{\partial J_{\theta}}{\partial z}$.
> (2.3b)

The operator acting on $B_{z}$ in eqns. (2.3a,b) is the scalar Laplacian
in cylindrical coordinates with $¶/¶q = 0$, while the operator acting on
$B_{r}$ is not quite of this standard form. In general the field
equations can be solved using the Green function for the Laplacian
operator. This may seem unnecessarily complicated since integrations are
over all three spatial variables (r,$\theta,$z) to obtain functions of
only the two variables (r,z). However for the case of greatest interest,
the thin wire layer, integrations can be easily reduced to the single
variable $\theta$, and the analysis is not much more complicated than
the textbook case of a single wire loop. It may seem simpler to use the
Biot-Savart formula (also a Green function solution), which yields
$\overrightarrow{B}$ from a single integration along the wire path - see
section 21, but this is somewhat misleading because the path may be very
long and must be specified by a more complicated formula than is used
for the smoothed $J_{\theta}$(r,z). On the other hand, the Biot-Savart
formula gives fine-scale detail of the field around the wire that is
missed by the Green function approach using smoothed sources, and is
therefore useful for detailed design. We also employ Fourier and Bessel
series and transformations, several types of magnetic potential, and for
really complicated cases a finite element code.

If cartesian components and coordinates are used then eqns. (2.3 a,b)
are replaced by the compact vector form

$Ñ^{2}\overrightarrow{B}\text{=-}m_{0}\overrightarrow{Ñ}x\overrightarrow{J}$
, (2.3c)

which can be immediately integrated using the Green function for the
Laplacian operator. It is then necessary to project $B_{r}$ from
$\overrightarrow{B}$; this is only a minor complication, and is
essentially equivalent to the direct solution of eqn. (2.3b) given in
section 6.

**2.2 Simple field estimate**

Generally we assume $J_{\theta} = 0$ at large
$\left| \overrightarrow{r} \right|$. However, to make a rough estimate
of field strength we examine the infinite, uniform, thin-current layer
carrying S Amperes/m:

$J_{q}(r,z) = Sd(r - R)$ all z, (2.4)

where $d(x)$ is the Dirac delta function. Eqns. (2.2 a,b) then give

$B_{z} = m_{0}\text{SH}(R - r),\begin{matrix}
 & 
\end{matrix}B_{r} = 0$. (2.5a,b)

Here H(x) denotes a positive unit step at x=0, with dH(x)/dx = $d(x)$:

$\begin{array}{r}
H(x) = \begin{array}{r}
\left\{ 1,\begin{matrix}
 & 
\end{matrix}x > 0, \middle| \right.\ 
\end{array} \\

\end{array}$

For a bore field of 1.0T we need

$S = \frac{B_{z}}{m_{0}} = \frac{1\text{.}0}{4px\text{10}^{- 7}} = \text{795}\text{.}8\frac{\text{kA}}{m}$
. (2.6)

In round numbers S=.80 MA/m for 1.0T (or 800 A/cm for 1.0 $k$G) for a
long magnet, but S must be somewhat higher for short magnets.
Superconducting wire or cable can carry very high current densities at
high fields, so a high field magnet can be made with a relatively thin
current layer. Examples of "engineering-averaged critical current
density" (see section 20) are:

> critical fields and current density at 4.2K

+-----------------------------------------------------------------------+
| superconductor wire field peak J                                      |
|                                                                       |
| $\text{Nb-Ti}$ 6.0T $\text{5x10}^{8}A/m^{2}$                          |
|                                                                       |
| $\text{Nb}_{3}\text{Sn}$ 10.0T $7\text{.}\text{5x10}^{8}A/m^{2}$      |
|                                                                       |
| $\text{Nb}_{3}\text{Sn}$ 15.0T $4\text{.}\text{8x10}^{8}A/m^{2}$      |
+-----------------------------------------------------------------------+

These current densities take into account space for Cu stabilizer,
insulation, and small gaps as well as actual superconductor. They are
lower than the critical current densities for the pure superconductor by
factors of 3 -10. For example, a bore field of 10.0T may require
S$$`<!-- -->`{=html}10MA/m, which can be made (in principle) using
$\text{Nb}_{3}\text{Sn}$ with current layer thickness
=$S/J = \text{10}^{7}/7\text{.}\text{5x10}^{8} = \text{.}\text{0133m}$.
When the radial distance from the magnet axis out to the superconducting
wire is greater than a few times the total layer thickness, we can model
the current density as an annulus of infinitesimal thickness for the
purpose of computing the field in the magnet bore. However, high field
solenoids, especially those using normal conductor, often have total
current layer thickness comparable to the bore radius or greater. This
more complicated geometry is analyzed in section 9.

**2.3 Current layers, global picture, and moments**

For a current layer of finite thickness it is useful to define the
cumulative current layer density

> $S(r,z)º\int_{}^{}{dr^{'}}J_{q}(r^{'},z)$, (2.7)

with net current layer density $S(z) = S(r = 0,z)$. If S(r,z) varys
slowly with z then eqns. (2.2a,b) suggest the approximate solution

$B_{z}»m_{0}S(r,z)$, (2.8a)

$B_{r}\text{»-}\frac{m_{0}}{r}\int_{}^{}{dr^{'}}r^{'}\frac{¶S(r^{'},z)}{¶z}$.
(2.8b)

However, most solenoids have abrupt wire layer ends, so eqns. (2.8a,b)
give only a rough approximation inside the magnet bore and wire, and
they are completely incorrect at large $r$. Sections 3-12 are devoted to
obtaining accurate formulas for all ranges of (r,z).

To avoid confusion over notation we mention here that the symbol $S$,
which denotes current layer density and is always a radial integration
of $J_{q}$, appears in four related forms depending on context: a simple
scale value $S$, the cumulative density $S(r,z)$, the net value $S(z)$,
and the cumulative density $S(r)$ which is present when $S(r,z)$ is
uniform in $z$ between specified magnet layer ends. Simple solenoids
composed of uniform thick or thin layers produce fields which are
suitable for beam transport. However, there are some applications that
require a more complicated layer function $S(r,z)$ within a single
magnet, for example to create a very uniform bore field. Such cases can
usually be constructed by superimposing the fields of several simple
layers. The symbol $I$ denotes either current in a wire or total
azimuthal current in a magnet, depending on context. In the latter case
$I = \int_{}^{}{\text{dz}\begin{matrix}

\end{matrix}S(z)}$.

Since $\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$ and
$\overrightarrow{B}$ vanishes at large
$\left| \overrightarrow{r} \right|$ (at least as fast as
$\left| \overrightarrow{r} \right|^{- 3}$), there must be a "return
flux". The total flux through any plane with fixed coordinate z is

> $\int_{}^{}{\text{dr}2pr}B_{z}(r,z) = 0$. (2.9)

All flux lines form closed loops at fixed values of $q$. For a very long
system the far field (i.e. the field at large
$\left| \overrightarrow{r} \right|$ greater than the system length) gets
mixed up with the Earth's field or the fields from other sources, so the
concept of return flux is actually somewhat ambiguous. However, the
external field produced close to the magnet system, usually with the
sign of $B_{z}$reversed from the adjacent bore field, is of great
interest -- see the figure below. For example, it may interfere with the
functioning of induction cores in a linac or with nearby
instrumentation. This topic is treated in several sections -- especially
section 11.

![](media/image1.pct){width="4.7362204724409445in"
height="3.416535433070866in"}

One rather expensive way to reduce or eliminate the far field is to
place an additional current layer with reverse polarity at larger radius
than the primary layer. For example a very long thin current layer of
density $S$ placed at radius $R$ can be surrounded by a second layer of
density $- S/4$ at $2R$, resulting in zero net flux in the region
$r < 2R$. This causes a reduction of the bore field to
$\text{.}\text{75}\begin{matrix}
m_{0}S
\end{matrix}$. A permeable annular yoke can also be used to confine the
return flux and reduce the external field, without any reduction of the
bore field -- see section 18.

For visualizing and calculating the field outside the magnet bore it is
useful to define "direct" and "residual" fields:

$\overrightarrow{B^{\text{direct}}}{ºm}_{0}S(r,z){\hat{e}}_{z}$, (2.10a)

> ${\overrightarrow{B}}^{\text{residual}}º\overrightarrow{B} - {\overrightarrow{B}}^{\text{direct}}$.
> (2.10b)

A "direct flux" is defined as

> $F\overset{\text{direct}}{(z)} = \int_{}^{}{\text{dr}2{\pi r\mu}_{0}}S(r,z) = \int_{}^{}{\text{dr}{\pi r}^{2}}\mu_{0}J_{\theta}(r,z)$,
> (2.11)

where we integrated by parts to obtain the last expression.

A single loop of wire of radius $R$ and carrying azimuthal current I has
magnetic dipole moment ${pR}^{2}I$. Magnetic moments add for a system
with solenoidal symmetry, so we have for the system's total magnetic
dipole moment $m_{1}$:

> $\begin{array}{r}
> m_{1} = \int_{}^{}\text{dz}\int_{}^{}{\text{dr}{pr}^{2}}J_{q} = \int_{}^{}\text{dz}\frac{F\overset{\text{direct}}{(z)}}{m_{0}} \\
> 2p\int_{\text{-¥}}^{\text{+¥}}\text{dz}\int_{}^{}{\text{drrS}(r,z)\begin{matrix}
> \text{.}
> \end{matrix}}
> \end{array}$ (2.12)

Assuming $m_{1}$ is non-zero, the field at very large
$\left| \overrightarrow{r} \right|$, i.e. the far field (see section
11), is

> $B_{r} = \frac{m_{0}m_{1}}{4p}\frac{3\text{rz}}{\sqrt{{r^{2} + z^{2}}^{5}}}$,
> (2.13a)
> $B_{z} = \frac{m_{0}m_{1}}{4p}\frac{2z^{2} - r^{2}}{\sqrt{{r^{2} + z^{2}}^{5}}}$
> . (2.13b)

Higher order moments, which generate field components that fall off more
rapidly with $\left| \overrightarrow{r} \right|$, are also present. The
coordinate origin can usually be positioned such that the magnetic
quadrupole moment
$\left( m_{2}\text{\textasciitilde}\int_{}^{}{\text{dz}\int_{}^{}{\text{drr}^{2}\text{zJ}_{q}}} \right)$
vanishes, but higher moments are always present.

So far we have made use of two radial moments of $J_{\theta}$:$$

> $S(z) = \int_{}^{}{\text{dr}\begin{matrix}
>
> \end{matrix}J_{q}(r,z)}$, (2.14a)
>
> $S(z)\overline{r^{2}}\overset{}{(z)} = \int_{}^{}{\text{dr}\begin{matrix}
>
> \end{matrix}r^{2}J_{q}(r,z) = 2\int_{}^{}{\text{drr}\begin{matrix}
>
> \end{matrix}S(r,z)}}$ . (2.14b)

Additional radial moments weighted by $J_{\theta}$ such as
$S\overline{r^{4}},S\overline{r^{6}}$, ... will be used in section (11)
to evaluate the external field. Note that
$\overline{r^{2}},\begin{matrix}
\overline{r^{4}}
\end{matrix}$, etc. can actually be negative if $J_{\theta}$ reverses
sign as a function of r, so they are not true mean values.

As mentioned, care is needed when representing a magnet by a smooth
current distribution. Consider a simple wire coil of $n$ turns at radius
$R$ with center-to-center longitudinal spacing $d$. If the wire current
is $I$ then the net azimuthal current is $\text{nI}$ and the net
magnetic moment is $m_{1} = {pR}^{2}\text{nI}$. However, the effective
coil length is $l = nd$, not $(n - 1)d$, since each loop would be
centered in a strip of length $d$ for an accurate finite element
calculation. The correct current layer density is
$S = \text{nI}/l = I/d$, which is consistent with the value
$m_{1} = {lpR}^{2}S$ obtained from eqn. (2.12). The physical length of
the coil alone, taking into account that it is actually a helix, is also
$l$ (plus one wire thickness). The longitudinal current of the coil is
$I$, which would be the source of an azimuthal field for $r > R$.
However, this is assumed to be cancelled by a second coil at slightly
larger radius.

4.  **Vector potential and flux lines**

It is often helpful to work with one of five different magnetic
potentials -- see sections 5 and 20. Here we only mention the vector
potential $\overrightarrow{A}$, with
$\overrightarrow{B} = \overrightarrow{Ñ}x\overrightarrow{A}$, which for
a solenoid may be reduced to the single component $A_{\theta}$(r,z):

$B_{r} = - \frac{\partial A_{\theta}}{\partial z},\begin{matrix}
B_{z} = \frac{1}{r}\frac{\partial}{\partial r}\text{rA}_{\theta}
\end{matrix}$. (2.15a,b)

A flux surface is given by

$\text{rA}_{q}(r,z) =$ constant. (2.16)

This is easily shown by considering any flux line that lies on a
particular flux surface. The line has constant $\theta$ and it is only
necessary to show that $\text{rA}_{\theta}$ is constant along it. Since
by definition the flux line element $d\overrightarrow{r}$ is parallel to
$\overrightarrow{B}$, the differential of $\text{rA}_{q}$ along the line
is

$d(\text{rA}_{\theta}) = (\frac{\partial}{\partial r}\text{rA}_{\theta})\text{dr} + (\frac{\partial\text{rA}_{\theta}}{\partial z})\text{dz} = r(B_{z}\text{dr} - B_{r}\text{dz}) = r{\hat{e}}_{\theta} \cdot \left( \overrightarrow{B}x\overrightarrow{d}r \right) = 0$.
(2.17)

If adjacent flux surfaces have constants $\text{rA}_{q} = C_{1}$ and
$C_{2}$, then the magnetic flux contained between the two surfaces is
$2p(C_{2} - C_{1})$. The mean field strength midway between the two flux
surfaces is approximately

$\left\langle B \right\rangle = \frac{C_{2} - C_{1}}{D \cdot r}$, (2.18)

where $D$ is the local normal distance between the surfaces and $r$ is
the radial distance to the midpoint between the surfaces. Due to the
factor of $r^{- 1}$ in eqn. (2.18), the magnitude of
$\overrightarrow{B}$ is not simply proportional to the inverse spacing
between flux surfaces $\left( D^{- 1} \right)$ as is sometimes assumed,
although this can be a useful visual aid in a region where $r$ is nearly
constant. The local density of flux lines through any plane normal to
$\overrightarrow{B}$ can be specified to be proportional to
$\overrightarrow{|B|}$, however this is not always done when making flux
line plots. A flux surface has the topology of a torus, and a system
made of more than one magnet may have several flux surfaces
corresponding to a given value of $\text{rA}_{\theta}$.

Another way to generate a flux line is to solve the equation

$\frac{dr}{\text{ds}} = \frac{\overrightarrow{B}(\overrightarrow{r})}{\left| \overrightarrow{B}(\overrightarrow{r}) \right|}$
, (2.19)

starting from some point on the line of interest. Here
$\overrightarrow{r}$(s) is simply the parameterization of the line, with
s denoting distance along the line. For example the starting point could
be on the surface of an acceleration gap where unwanted electrons are
emitted.

**2.5 Stress and energy**

Material stress is not the topic of this report, but it is of major
concern in the design of high field magnets, so a simple estimate is
presented here. The force density from the field in the current carrying
wire is

$\overrightarrow{J}x\overrightarrow{B} = {\hat{e}}_{r}J_{q}B_{z} - {\hat{e}}_{z}J_{q}B_{r}$.
(2.20)

Since the wire is in mechanical equilibrium, this force density must
equal the divergence of the material stress tensor, which itself is
proportional to components of the gradient of the microscopic mechanical
displacement field in the magnet. A detailed analysis is too complicated
to summarize here (see e.g. ref \[1-2\]), but it is seen from eqn.
(2.20) that the longitudinal magnetic force tends to compress a wire
layer while the radial magnetic force pushes it outwards against a
restraining collar of the magnet assembly. There can be no radial
pressure in the vacuum bore, so at the collar it is approximately

$P»\int_{}^{}\text{dr}J_{q}B_{z}»\frac{B_{w}^{2}}{2m_{0}}$, (2.21)

where $B_{w}$ is the longitudinal field at the inner wire edge of the
coil -- roughly the peak bore field. Suppose $B_{w} = 1\text{.}0T$, then
$P»\left( (1\text{.}0)^{2}/2x4px\text{10}^{- 7} \right)»4x\text{10}^{5}Pa»4\text{.}0\begin{matrix}
\text{Bar}
\end{matrix}$, which is quite modest. However a 15.0T bore field
produces roughly 900 Bar, which can crush some insulation.

The external field of a solenoid without a permeable return yoke extends
to a large distance compared to that of a quadrupole magnet. The
estimation of stresses, wire field, etc. of a system of solenoids must
therefore include the field contributions of nearby magnets as well as
that of the magnet of interest. This is particularly the case for a
transport lattice made of closely spaced solenoids.

Since solenoidal transport systems for energetic ions are expected to
have high fields in large-bore magnets, they also have large stored
field energy. In the absence of magnetic materials this is

$W = \int_{}^{}d^{3}r\frac{\left| \overrightarrow{B} \right|^{2}}{2m_{0}}$
. (2.22)

For example if a 1.0 $m^{3}$ volume is filled with a 15T field, then the
stored energy is $\text{15}/\left( 2x4px\text{10}^{- 7} \right)$ $»$
90MJ.

3\. [Near-Axis Field]{.underline}

For simulating beam dynamics it is often sufficient to accurately
approximate the magnetic field near the system axis, e.g. out to 50% of
the radius of the nearest wire layer. A general treatment is presented
here and some useful examples of the on-axis field are given in section
4.

3.1 On-axis field and near-axis expansion

Recall eqn. (2.3a) for $B_{z}$:

$Ñ^{2}B_{z}\text{=-}m_{0}\frac{1}{r}\frac{¶}{¶r}\text{rJ}_{q}$. (3.1)

This may be formally integrated using the Green function for the 3-d
scalar Laplacian operator (see section 6):

$B_{z} = \frac{m_{0}}{4p}\int_{}^{}d^{3}r^{'}$$\frac{\frac{1}{r^{'}}\frac{¶}{¶r^{'}}r^{'}J_{q}^{'}}{\left| r^{'} - \overrightarrow{r} \right|}$
, (3.2)

or explicitly in cylindrical coordinates,

$B_{z} = \frac{\mu_{0}}{4\pi}\begin{matrix}

\end{matrix}$ . (3.3)

Here and in similar integrations we may set $\theta = 0$inside the
integrand since by symmetry $B_{z}$ is a function only of r and z. After
an integration by parts in $r^{'}$ we get

$B_{z}(r,z) = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}J_{q}}\left( r^{'},z^{'} \right)\int_{}^{}{dq^{'}}\frac{\left( r^{'} - r\text{cos}\left( q^{'} \right) \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + r^{} + r^{2} - 2r^{'}r\begin{matrix}
\text{cos}\left( q^{'} \right)
\end{matrix}}^{3}}}$. (3.4)

This formula is completely general for a static solenoidal field, and it
gives the right answer even when the current layer extends to infinity
in $z$.

To get the near-axis field, inside the bore where $J_{\theta} = 0$, we
first define the on-axis field:

$B_{0}(z) \equiv B_{z}(r = 0,z) = \frac{\mu_{0}}{2}\int_{}^{}{dz^{'}\int_{}^{}{dr^{'}}}\frac{r^{}J_{\theta}\left( r^{'},z^{'} \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + r^{}}^{3}}}$
. (3.5)

If $B_{0}(z)$ is accurately determined from this expression or in some
other way, such as measurements on a real magnet assembly, then
$B_{0}(z)$ may be used to generate a power series solution of eqn.
(3.1), which is valid in the magnet bore. From

$\frac{1}{r}\frac{¶}{¶r}r\frac{¶B_{z}}{¶r}\text{=-}\frac{¶^{2}B_{z}}{¶z^{2}}$
,

we get by iteration:

$B_{z}(r,z) = B_{0} - \frac{B_{0}^{''}r^{2}}{4} + B_{0}^{''''}\frac{\begin{matrix}
r^{4}
\end{matrix}}{\text{64}} - \frac{B_{0^{\nu i}}r^{6}}{\text{2304}} + \text{.}\text{.}\text{.}$
. (3.6a)

The radial field component is easily obtained from
$\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$ by integrating in r:

$B_{r}(r,z)\text{=-}\frac{1}{r}\int_{}^{}{dr^{'}}r^{'}$$\frac{¶B_{z}\left( r^{'},z \right)}{¶z}\text{=-}\frac{B_{0}^{'}r}{2} + \frac{B_{0}^{'''}r^{3}}{\text{16}} - \frac{B_{0}^{\nu}r^{5}}{\text{384}} + \frac{B_{0}^{\nu\text{ii}}r^{7}}{\text{18432}} - \text{.}\text{.}\text{.}\text{.}\text{.}$
. (3.6b)

The leading term in each of the above expansions of $B_{z}$and $B_{r}$is
used when computing the linearized dynamics of beam ions. Higher order
(non-linear in $r$) terms are referred to as "fringe field aberrations"
in this context. A few terms from the expansions usually give an
accurate field representation inside the vacuum bore, especially near
the magnet center, but the expansion for $B_{z}$fails in the wire and
beyond, where the field component generally reverses sign. The expansion
for $B_{r}$ often gives a fair approximation in the wire except near the
magnet ends. When truncating the expansions after some power of r,
$B_{r}$should be represented by one power higher than $B_{z}$in order to
preserve $\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$.
Crossing a thin wire layer of total surface current layer density
$S(z)$, $B_{z}$ jumps by $- m_{0}S(z)$ but $B_{r}$ is continuous.
However $B_{r}$ has logarithmic singularities at the ends of the layer
-- see section 8. For large $|z|$, well beyond any magnet end, the power
series expansions are valid out to values of r comparable to the
longitudinal distance from the magnet end. But the expansions do not
correctly represent the reversal of $B_{z}$ at larger values of $r$.
Also, because of the high order derivatives, the expansions of $B_{r}$
and $B_{z}$ are difficult to apply using measured $B_{0}$.

By applying Stoke's theorem to
$\overrightarrow{Ñ}x\overrightarrow{B} = m_{0}\overrightarrow{J}$, we
obtain the useful general result:

> $\int_{}^{}{\text{dz}\begin{matrix}
>
> \end{matrix}B_{0}(z) = m_{0}\int_{}^{}{\text{dz}\begin{matrix}
> S(z)
> \end{matrix}}}$ . (3.7)

This may also be derived from eqn. (3.5) using the indefinite integral
(3.10a) given at the end of this section.

For the thin semi-infinite wire layer,

$J_{q} = \text{SH}( - z)d(r - R)$ , (3.8)

we get from eqn. (3.5) - see section 4:

$B_{0}(z) = \frac{\mu_{0}S}{2}\left( 1 - \frac{z}{\sqrt{z^{2} + R^{2}}} \right)$.
(3.9)

In this representative case fringe field aberrations are maximum at z=0
for $B_{r}$ and at $z \approx \pm \text{.}5R$ for $B_{z}$. Approximating
the fields by only two terms works well for r\<.5R, but for r=.8R (the
largest value one would consider for beam transport) the convergence is
slow, as shown in the following tables. Note that the approximated
$B_{z}$ oscillates as higher order terms are added.

+-----------------------------------------------------------------------+
| Number of expansion terms $\frac{B_{z}}{\mu_{0}}S$at r=.8R, z=.5R     |
| $\frac{B_{r}}{\mu_{0}}S$at r=.8R, z=0                                 |
|                                                                       |
| 1 .276393 .200000                                                     |
|                                                                       |
| 2 .207701 .248000                                                     |
|                                                                       |
| 3 .190116 .267200                                                     |
|                                                                       |
| 4 .189066 .276160                                                     |
|                                                                       |
| 5 .190881 .280676                                                     |
|                                                                       |
| exact. .191960 .286062                                                |
+-----------------------------------------------------------------------+

Number of expansion terms $\frac{B_{z}}{\mu_{0}}S$at r=.5R, z=.5R
$\frac{B_{r}}{\mu_{0}}S$at r=.5R, z=0

1 .276393 .125000

2 .249560 .136719

3 .246877 .138550

4 .246814 .138884

5 .246857 .138949

exact. .246867 .138967

**3.2 Useful integrals**

The following (easily verified) indefinite integrals are of use in
subsequent sections:

$\int_{}^{}\frac{\text{dx}}{\sqrt{{x^{2} + a^{2}}^{3}}} = \frac{x}{a^{2}\sqrt{x^{2} + a^{2}}}$,
(3.10a)

$\int_{}^{}\frac{\text{dx}}{\sqrt{x^{2} + a^{2}}} = \text{Log}\left( \frac{x + \sqrt{x^{2} + a^{2}}}{a} \right) = \text{sinh}^{- 1}(x/a)$
, (3.10b)

$\int_{}^{}{\frac{\text{dx}}{x^{2} - a^{2}}\text{=-}\frac{1}{a}}\text{tanh}^{- 1}\left( \frac{x}{a} \right)$
, (3.10c)

$\int_{}^{}\frac{\text{dx}}{\left( \text{ax}^{2} + b \right)\sqrt{\text{fx}^{2} + g}} = \frac{1}{\sqrt{b\left( \text{ag} - \text{bf} \right)}}\text{tan}^{- 1}\left( \frac{x\sqrt{\text{ag} - \text{bf}}}{\sqrt{b\left( \text{fx}^{2} + g \right)}} \right)\begin{matrix}
,
\end{matrix}$ (3.10d)

\[for $a,b,f,g > 0$ and $\text{ag} - \text{bf} > 0$\].

The surprising integral (3.10d) is used in section 9 (along with 3.10b
and 3.10c) to evaluate $B_{z}(r,z)$ for a uniform, thick current layer.

4\. [Design Formulas for the On-Axis Field]{.underline}

4.1 Simple cases

To calculate the on-axis field $B_{0}(z)$ we may replace a moderately
thick, uniform wire layer with a thin layer at its average radius (R)
with surprisingly little error -- see section 9. Making this
approximation results in a number of very simple formulas. The current
density has the simple form:

> $J_{q}(r,z) = S(z)d(r - R)$ , (4.1)

and eqn. (3.5) becomes

> $B_{0}(z) = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\frac{R^{2}S\left( z^{'} \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + R^{2}}^{3}}}$
> . (4.2)

[Case 1]{.underline} Constant $S(z) = S$ for all z:

Using the integral (3.10a) we recover eqn. (2.5a),

> $B_{0} = \frac{m_{0}S}{2}\frac{\left( z^{'} - z \right)}{\sqrt{\left( z^{'} - z \right)^{2} + R^{2}}}\underset{\text{-¥}}{\overset{\text{+¥}}{|}} = m_{0}S$
> . (4.3)

[Case 2]{.underline} Semi-infinite Layer:

> $S(z) = \begin{array}{r}
>  \\
>
> \end{array}$ (4.4)

A very useful formula is derived by changing the upper limit of
integration over $z^{'}$ to $z_{0}$ in eqn. (4.3):

> $B_{0}(z) = \frac{m_{0}S}{2}\frac{\left( z^{'} - z \right)}{\sqrt{\left( z^{'} - z \right)^{2} + R^{2}}}\underset{\text{-¥}}{\overset{z_{0}}{|}} = \frac{m_{0}S}{2}\left\lbrack 1 - \frac{\left( z - z_{0} \right)}{\sqrt{\left( z - z_{0} \right)^{2} + R^{2}}} \right\rbrack$
> . (4.5)

Note the symmetry around $z_{0}$, with
$B_{0}\left( z_{0} \right) = m_{0}S/2$. We may build up many other cases
of interest by superposition using this result.

[Case 3]{.underline} Lens of wire layer length $\ell$, centered at
$z = 0$:

We apply eqn. (4.5) with positive and negative current layers ending
respectively at $l/2$ and $- l/2$:

$B_{0}^{\text{lens}}(z) = \frac{m_{0}S}{2}\left\lbrack \frac{(z + l/2)}{\sqrt{(z + l/2)^{2} + R^{2}}} - \frac{(z - l/2)}{\sqrt{(z - l/2)^{2} + R^{2}}} \right\rbrack$
. (4.6)

![](media/image3.pct){width="4.381944444444445in"
height="1.2173611111111111in"}This formula is equivalent to the
well-known textbook result \[4-1\]:

> $B_{0} = \frac{m_{0}\text{NI}}{2}\left\lbrack \text{cos}\left( q_{1} \right) + \text{cos}\left( q_{2} \right) \right\rbrack$
> , (4.7)

where N is the number of wire loops per meter, each carrying current I,
and the angles are to points at the wire layer ends drawn from the point
of interest on the z axis (see figure). The field at the lens midpoint
$(r = 0,z = 0)$ is

> $B_{0}^{\text{lens}}(0) = \frac{m_{0}S}{\sqrt{1 + 4R{}^{2}/l^{2}}}$ ,
> (4.8)

so a lens of length $l = 2R$ has peak field on axis reduced
to$m_{0}S/\sqrt{2} = \text{.}\text{7071}m_{0}S$. For $|z|$ large, the
leading order expansion term of eqn. (4.6) is

> $B_{0} \rightarrow \frac{m_{0}\text{SR}^{2}l}{2|z|^{3}}$ , (4.9)

consistent with eqn. (2.13b). The on-axis field of a short lens

![](media/image4.pct){width="5.996527777777778in"
height="3.6381944444444443in"}$$is displayed below.

[Case 4]{.underline} Periodic system of lenses of wire layer length $l$,
gap length g, and period P=$l$+g:

> $B_{0}(z) = \sum_{n\text{=-¥}}^{¥}{}B_{0}^{\text{lens}}\left( z - \text{nP} \right)$
> (4.10)
>
> [Case 5]{.underline} Gap of length g, centered at z=0, between long
> solenoids:

We subtract a lens of length g from an infinitely long solenoid;

> $B_{0}^{\text{gap}}(z) = m_{0}S - B_{0}^{\text{lens}}(z)$ . (4.11)

At the center of the gap

$B_{0}^{\text{gap}}(0) = m_{0}S\left( 1 - \frac{1}{\sqrt{1 + 4R^{2}/g^{2}}} \right)$
. (4.12)

For g=2R, the field at the gap center is reduced to
$\left( 1 - 1/\sqrt{2} \right)m_{0}S = \text{.}\text{2929}m_{0}S$. To
prevent large variations of $B_{0}$ with z in a periodic system, the
gaps must be very short compared with 2R. However, this is not generally
necessary for the design of a good transport system if the ion betatron
wavelength is longer than 2P.

[Case 6]{.underline} Current ring of radius $R$:

> $J_{q} = Id(r - R)d(z)$ . (4.13)

From eqn. (3.5) or eqn. (4.2) we get

$B_{0}^{\text{ring}}(z) = \frac{m_{0}\text{IR}^{2}}{2\sqrt{{z^{2} + R^{2}}^{3}}}$
. (4.14)

4.2 Thick wire layer

The on-axis field of a uniform, thick wire layer can be computed from
any of the thin layer cases by a simple average over the radial position
$R$:

> $B_{0}^{\text{thick}} = \int_{}^{}\frac{\text{dR}}{R_{2} - R_{1}}\begin{matrix}
> B_{0}^{\text{thin}}(z;R)
> \end{matrix}$ , (4-15)

where $R_{1}$ and $R_{2}$ are the inner and outer radii of the layer.
For example eqn. (4.6) for a simple lens of length $l$ gives, using the
integral eqn. (3.10b):

$B_{0}(z) = \int_{}^{}\frac{\text{dR}}{R_{2} - R_{1}}\begin{matrix}
\frac{m_{0}S}{2}\left\{ \frac{(z + l/2)}{\sqrt{(z + l/2)^{2} + R^{2}}} - \frac{(z - l/2)}{\sqrt{(z - l/2)^{2} + R^{2}}} \right\}
\end{matrix}$

$= \frac{\mu_{0}S}{2\left( R_{2} - R_{1} \right)}\left\{ (z + \ell/2)\text{Log}\left\lbrack \frac{R_{2} + \sqrt{(z + \ell/2)^{2} + R_{2^{2}}}}{R_{1} + \sqrt{(z + \ell/2)^{2} + R_{1^{2}}}} \right\rbrack \right.\  -$

$- \left. \ (z - \ell/2)\text{Log}\left\lbrack \frac{R_{2} + \sqrt{(z - \ell/2)^{2} + R_{2^{2}}}}{R_{1} + \sqrt{(z - \ell/2)^{2} + R_{1^{2}}}} \right\rbrack \right\}\text{.}$
(4.16)

4.3 Uniform fields

[Case 1]{.underline} Helmholtz coil:

Two rings, each of radius $R$ and current $I$, which are separated by
distance $L$, have on-axis field, from eqn. (4.14):

$B_{0}(z) = \frac{m_{0}\text{IR}^{2}}{2}\left\lbrack \frac{1}{\sqrt{{(z - L2)^{2} + R^{2}}^{3}}} + \frac{1}{\sqrt{{(z + L2)^{2} + R^{2}}^{3}}} \right\rbrack$
. (4.17)

If $L = R$ then $B_{0^{¢¢}}(0) = 0$. This is a simple way to make a
uniform field close to the coil center, but it is not generally useful
for beam focusing since the field strength in the wire is much larger
than at the center.

[Case 2]{.underline} Nearly uniform field in a long solenoid:

A large volume of strong uniform field can be made by superposing simple
solenoids of various lengths, radii and current layer densities. A
mathematically simple (but unphysical) version consists of a pair of
thin-layer lenses \[see eqn. (4.6)\] placed symmetrically around $z = 0$
at the same radius $R$, but with unequal lengths
$\left( l_{1},l_{2} \right)$ and layer densities
$\left( S_{1},S_{2} \right)$. The derivatives $B_{0^{¢¢}}(0)$ and
$B_{0}^{iu}$ both vanish if we set

$\frac{l_{2}}{2R} = \sqrt{\frac{3\left( l_{1}/2R \right)^{2} + \text{10}}{4\left( l_{1}/2R \right)^{2} - 3}}$
, (4.18)

$S_{2}\text{=-}S_{1}\left( \frac{l_{1}}{l_{2}} \right)\sqrt{\frac{1 + \left( l_{2}/2R \right)^{2}}{1 + \left( l_{1}/2R \right)^{2}}^{5}}$
. (4.19)

Layer #2 is shorter than #1 if $l_{1}/2R\begin{matrix}
 > \sqrt{2\text{.}5}
\end{matrix}$, and has a reversed current. Both layers must have
$l/2R\begin{matrix}
 > \sqrt{\text{.}\text{75}}
\end{matrix}$. The combined layers produce a somewhat lower field
strength than that of layer #1 alone:

$B_{0}(0) = \frac{\mu_{0}S_{1}\ell_{1}/2}{\sqrt{\left( \ell_{1}/2 \right)^{2} + R^{2}}} + \frac{\mu_{0}S_{2}\ell_{2}/2}{\sqrt{\left( \ell_{2}/2 \right)^{2} + R^{2}}}$
. (4.20)

For example, taking $m_{0}S_{1} = \text{10}\text{.}0T$ and
$l_{1} = 4\text{.}0R$, we get
$m_{0}S_{2}\text{=-}3\text{.}\text{27097}T,\begin{matrix}

\end{matrix}$, producing $B_{0}(0) = 6\text{.}\text{35096}T$. The #1
solenoid alone would produce $B_{0}(0) = 8\text{.}\text{94427}T$.
On-axis fields for the double layer at representative points are

  --------------------- ------------------ ------------------ ------------------ ------------------ ------------------ ------------------
  $$\begin{matrix}      $$\begin{matrix}   $$\begin{matrix}   $$\begin{matrix}   $$\begin{matrix}   $$\begin{matrix}   $$\begin{matrix}
   &                     &                  \pm \text{.}2      \pm \text{.}4      \pm \text{.}6      \pm \text{.}8      \pm 1\text{.}0
  \end{matrix}z/R$$     \end{matrix}0$$    \end{matrix}$$     \end{matrix}$$     \end{matrix}$$     \end{matrix}$$     \end{matrix}$$

  $$\begin{matrix}      6.35096            6.35096            6.35083            6.34929            6.34063            6.30777
   &                                                                                                                   
  \end{matrix}B_{0}$$                                                                                                  
  --------------------- ------------------ ------------------ ------------------ ------------------ ------------------ ------------------

There is essentially no variation over 25% of the bore length in this
example and only .7% variation over 50% of the bore length. Transverse
variation is similarly small in these zones.

[Case 3]{.underline} Uniform field inside an ellipsoidal shell:

Although this is not particularly useful for beam transport, it is worth
mentioning that an ellipsoidal shell of constant current per meter (in
z) produces an exactly constant field inside the enclosed volume.
Denoting the cumulative current layer density

$S(r,z) = S \cdot H\left( 1 - \frac{r^{2}}{a^{2}} - \frac{z^{2}}{b^{2}} \right)$
, (4.21)

where a and b are ellipse radii, then

$J_{q}\text{=-}\frac{¶S(r,z)}{¶r} = \frac{2r}{a^{2}}S\begin{matrix}
d\left( 1 - \frac{r^{2}}{a^{2}} - \frac{z^{2}}{b^{2}} \right)
\end{matrix}$ . (4.22)

This current pattern can be created with closely spaced wire loops each
carrying current $I = SD$, where $D$ is the constant longitudinal
separation between loops.

Eqn. (3.5) yields

$B_{0}(z) = \frac{m_{0}\text{Sa}^{2}}{2}\int_{- b}^{b}{dz^{'}}\frac{\left( 1 - z^{'}2b^{2} \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + a^{2}\left( 1 - z^{'}2b^{2} \right)}^{3}}}$
, (4.23)

from which it may be numerically verified that $B_{0}(z)$ is actually
constant for $|z| < b$. $B_{z}$ is therefore also independent of $r$
inside the ellipsoid, and we may write

$B_{z}^{\text{inside}} = m_{0}\text{Sf}(a/b)$ , (4.24)

$f(a/b) = \frac{a^{2}}{b^{2}}\int_{}^{}\frac{\text{dy}\begin{matrix}
\left( 1 - y^{2} \right)
\end{matrix}}{\sqrt{{y^{2} + \frac{a^{2}}{b^{2}}\left( 1 - y^{2} \right)}^{3}}}$
. (4.25)

  ------------------- ------------------ --------------------------------- --------------------------------- -------------------------- -------------------------- -------------------------- ------------------
  $$\begin{matrix}    $$\begin{matrix}   $$\begin{matrix}                  $$\begin{matrix}                  $$\begin{matrix}           $$\begin{matrix}           $$\begin{matrix}           $$\begin{matrix}
   &                   &                  &                                 &                                 &                          &                          &                          & 
  \end{matrix}a/b$$   \end{matrix}0$$    \end{matrix}\text{.}\text{25}$$   \end{matrix}\text{.}\text{50}$$   \end{matrix}1\text{.}0$$   \end{matrix}2\text{.}0$$   \end{matrix}4\text{.}0$$   \end{matrix}¥$$

  $$\begin{matrix}    1.00000            .924593                           .826436                           .666667                    .472800                    .296359                    .000000
   &                                                                                                                                                                                          
  \end{matrix}f$$                                                                                                                                                                             
  ------------------- ------------------ --------------------------------- --------------------------------- -------------------------- -------------------------- -------------------------- ------------------

A proof that $B_{z}$ is constant is mathematically related to the
solution for the gravitational potential of an ellipsoid of constant
density, obtained in the Eighteenth Century \[4-2\]. It is found by
elementary integration that for $a/bºx < 1$

$f(x) = \frac{1}{1 - x^{2}}\left\lbrack 1 - \frac{x^{2}}{\sqrt{1 - x^{2}}}\text{Log}\left( \frac{\text{1+}\sqrt{\text{1-x}^{2}}}{x} \right) \right\rbrack$
, (4.26a)

and for $a/b = x > 1$ ,

$f(x) = \frac{1}{x^{2} - 1}\left\lbrack \frac{x^{2}}{\sqrt{x^{2} - 1}}\text{sin}^{- 1}\left( \frac{\sqrt{x^{2}\text{-1}}}{x} \right) - 1 \right\rbrack$
. (4.26b)

If several elliptical shells of the form (4.21 -- 4.22) are superposed
and all have the same ratio $a/b$, then the interior field is still
given by (4.24), but with $S\text{=å}\begin{matrix}
S_{i}
\end{matrix}$.

References

4-1 See e.g. Reference \[1-1\], J.D. Jackson, pg. 225.

4-2 Chandrasekhar, S., Ellipoidal Figures of Equilibrium, Chapter 3,
Dover Publications,

New York, 1969.

5\. [Vector and Scalar Potentials]{.underline}

The vector potential $A_{q}$ and four scalar potentials, denoted here by
$\varphi$, $j$, $F$ and $Y$, are of use in calculating and visualizing
solenoidal fields. Each potential has features that make it useful for
some applications and a poor choice for others. The scalar potentials
$f\begin{matrix}

\end{matrix}$ are defined below, while $F$ and$Y$ are special potentials
used in Tosca and described in section 19.

5.1 Vector potential ($A_{q}$)

As discussed in section 2, the vector potential is useful for defining
and plotting flux surfaces. Another property of $A_{\theta}$is that the
field components

> $B_{r}\text{=-}\frac{¶A_{q}}{¶z},\begin{matrix}
>  & 
> \end{matrix}B_{z} = \frac{1}{r}\frac{¶}{¶r}\text{rA}_{q}$ , (5.1 a,b)

automatically satisfy $\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$.
This is an advantage in beam dynamics applications, where the
conservation of canonical angular momentum

$P_{q} = r\left( g\text{Mv}_{q} + \text{qeA}_{q} \right)$ (5.2)

for an ion of $q -$ momentum ${gM\nu}_{q}$ and charge $\text{qe}$, might
otherwise be violated by an approximate representation of the field;
e.g. constant $P_{q}$ can be used as an equation of motion.

Unlike the scalar potential $f$ described below, the vector potential
provides a complete (or "global") characterization of the magnetic
field, which is valid everywhere. Its source is $J_{\theta}$, and if a
magnetization current density $J_{q}^{\text{mag}}$ is present, it simply
adds to the current density of the wire layer. The second order equation
for $A_{\theta}$ is

$\frac{¶}{¶r}\frac{1}{r}\frac{¶}{¶r}\text{rA}_{q} + \frac{¶^{2}A_{q}}{¶z^{2}} = \frac{¶B_{z}}{¶r} - \frac{¶B_{r}}{¶z}\text{=-}m_{0}J_{q}$
. (5.3)

As with eqn. (2.3b) for $B_{r}$, the operator acting on $A_{\theta}$ is
not quite the scalar Laplacian, but this does not cause a serious
problem in obtaining a Green function solution -- see section 6.

One drawback in using $A_{\theta}$(or the scalar potentials) to
calculate$\overrightarrow{B}$ is that it must be differentiated, which
is an extra operation that can produce unnecessarily complicated
expressions. Also, derivatives generate extra errors in discretized,
numerical models. Another drawback is that $A_{q}$is the azimuthal
component of a vector field and must therefore be handled carefully when
transforming from a tilted magnet.

5.2 Magnetic scalar potential ($f$)

The field in the bore of a magnetic dipole, quadrupole, or higher order
magnetic beam line element is often characterized by a magnetic scalar
potential $\varphi$. For a connected region of vacuum,
where$\overrightarrow{Ñ}x\overrightarrow{B} = 0$, we can set

$\overrightarrow{B} = m_{0}\overrightarrow{Ñ}f,\begin{matrix}
 & 
\end{matrix}Ñ^{2}f = \overrightarrow{Ñ} \cdot \overrightarrow{B}/m_{0} = 0$
, (5.4 a,b)

where the permeability of vacuum $m_{0}$ is inserted to be consistent
with standard convention. This is an attractive, compact formalism for
aberrations, which are represented as higher order multipoles, e.g.
$f\text{\textasciitilde}f_{n}(r,z)\text{cos}(nq)$, as well as fringe
field terms. A magnetic scalar potential may also be used to
characterize the bore field in a solenoidal system.

By itself $\varphi$ cannot be a global solution. It is necessary to
connect $\varphi$ to its current source through boundary conditions or
an already known component of $\overrightarrow{B}$. For example
$\varphi$ may be assumed constant on the surface of a highly permeable
material if the magnetic field is well below saturation, e.g.
$B \leq 1\text{.}0T$ for magnetically soft iron. In fact this formalism
is most useful when magnetic materials are present, with
$\overrightarrow{H} = \overrightarrow{Ñ}f$ everywhere outside the wire
and $\overrightarrow{B} = m_{0}\overrightarrow{H}$ in vacuum --see
sections 16 and 17. Another difficulty with $\varphi$ is that it is
usually multi-valued. This follows from
$\oint_{}^{}\overrightarrow{B} \cdot \overrightarrow{\text{dr}}¹0$
integrated around any loop that encloses a current layer. The idea is to
position the discontinuity (if any) somewhere such that a single branch
of $f$ can be applied in the region of interest.

Surfaces with $\varphi$ = constant have perpendicular intersections with
surfaces with $\text{rA}_{q}$ = constant. The constant $\varphi$
surfaces may be specified such that their separation is inversely
proportional to $\left| \overrightarrow{B} \right|$ throughout the
vacuum region where a scalar potential is valid.

5.3 Global scalar potential ($j$)

The limitations of the magnetic scalar potential may be circumvented for
a solenoid by defining a global scalar potential $j$. Let

$\overrightarrow{B} = m_{0}S(r,z){\hat{e}}_{z} + \overrightarrow{Ñ}j$,
(5.5)

where $S(r,z)$ is the cumulative current layer density defined in
section 2 and $\overrightarrow{J^{\text{mag}}} = 0$.

This formula for $\overrightarrow{B}$ satisfies
$\overrightarrow{\nabla}x\overrightarrow{B} = \mu_{0}\overrightarrow{J}$.
To make $\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$ we set

$Ñ^{2}j\text{=-}\frac{¶}{¶z}m_{0}S(r,z)$ . (5.6)

The term $m_{0}S(r,z){\hat{e}}_{z}$ in eqn. (5.5) is just the previously
defined direct field of the system (see section 2), and
$\overrightarrow{Ñ}j$ is the residual field (i.e. fringe field, external
field, etc.). It is clear from eqn. (5.5) that $j$ is single valued
since $\overrightarrow{B} - m_{0}S(r,z){\hat{e}}_{z}$ has vanishing
curl. Also, $j$ is a global solution since eqn. (5.6) is valid over all
space. In principle a global scalar potential can be defined for any
configuration of currents, but it is not useful unless
$\overrightarrow{Ñ}x\overrightarrow{B} = m_{0}\overrightarrow{J}$ has a
simple solution, as is the case for a solenoidal system.

Note that $\overrightarrow{Ñ}j$ is discontinuous at the ends of a magnet
bore and it generally subtracts from the direct field inside the bore.
This is a very awkward way to treat the field close to the axis, where
it is used for beam dynamics calculations, and continuity may be
important. However, $j$ gives an intuitively appealing picture of the
external field by an analogy with electrostatics, where the source on
the rhs of eqn. (5.6) acts like a charge density proportional to
$¶S(r,z)/¶z$. For a thin current layer at radius R, this source is a
pair of disks of opposite magnetic charge
$\left( \pm {pR}^{2}m_{0}S \right)$ at the solenoid ends. In section 11
the global scalar potential will be used to approximate the external
field, which will in turn be used to calculate the interference of
fields among off-axis solenoids in section 15.

**6. [Green Function Solutions]{.underline}**

**6.1 Green function formalism**

Any equation of the form (in 3-d)

> $Ñ^{2}f = g\left( \overrightarrow{r} \right)$ , (6.1)

with $g\left( \overrightarrow{r} \right) = 0$ at large
$\left| \overrightarrow{r} \right|$, may be formally integrated using
the Green function

> $G\left( \overrightarrow{r},\overrightarrow{r^{'}} \right) = \frac{- 1}{4p\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
> ; (6.2)
>
> $f\left( \overrightarrow{r} \right) = \int_{}^{}d^{3}r^{'}G\left( \overrightarrow{r},\overrightarrow{r^{'}} \right)g\left( \overrightarrow{r^{'}} \right)$
> . (6.3)

This works because

$Ñ^{2}G = Ñ^{}G = d\left( \overrightarrow{r^{'}} - \overrightarrow{r} \right)$
, (6.4)

where $d\left( \overrightarrow{x} \right)$ denotes the 3-d Dirac delta
function.

We have already used
$G\left( \overrightarrow{r},\overrightarrow{r^{'}} \right)$ to solve for
$B_{z}$ in section (3):

> $B_{z} = \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}}\left( \frac{1}{r^{'}}\frac{¶}{¶r^{'}}r^{'}J_{q}^{'} \right)\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|} = \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}}J_{q}^{'}\frac{r^{'} - r\text{cos}\left( q^{'} - q \right)}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|^{3}}$.
> (6.5)

Equation (5.6) for the gobal scalar potential $j$ is also of the proper
form for application of the Green function:

> $j = \frac{m_{0}}{4p}\int_{}^{}d^{3}r^{'}\left( \frac{¶S^{'}}{¶z^{'}} \right)\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
> . (6.6)

Equations (2.3b) for $B_{r}$ and (5.3) for $A_{q}$ do not have the
Laplacian operator, but they may be converted to the desired form by
multiplying both sides by $\text{cos}q$, i.e.

> $Ñ^{2}\left( B_{r}\text{cos}q \right) = \left( \frac{¶}{¶r}\frac{1}{r}\frac{¶}{¶r}\text{rB}_{r} + \frac{¶^{2}B_{r}}{¶z^{2}} \right)\text{cos}q = m_{0}\frac{J_{q}}{¶z}\text{cos}q$
> . (6.7)

This yields

> $B_{r}\text{=-}\frac{m_{0}}{4p}\int_{}^{}d^{3}r^{'}\left( \frac{¶J_{q}^{'}}{¶z^{'}} \right)\frac{\text{cos}\left( q^{'} - q \right)}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
> , (6.8)

where we have used the identity

$\text{cos}\left( q^{'} \right) = \text{cos}(q)\text{cos}\left( q^{'} - q \right) - \text{sin}(q)\text{sin}\left( q^{'} - q \right)$

inside the Green function integral. In similar fashion we get

> $A_{q} = \frac{m_{0}}{4p}\int_{}^{}d^{3}r^{'}J_{q^{¢}}\frac{\text{cos}\left( q^{'} - q \right)}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
> . (6.9)

The four Green function solutions given above all have denominators that
are odd powers of the function

> $\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right| =$. (6.10)

Several simplifications can be made. First, by symmetry, we may let
$\theta^{'} - \theta \rightarrow \theta^{'}$everywhere in an integrand.
Second, if a derivative with respect to $z$or $z^{'}$appears in the
integrand, we may use the identity

> $\frac{¶}{¶z}\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}\text{=-}\frac{¶}{¶z^{'}}\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
> , (6.11)

integrate by parts, and sometimes get a simplification. For example, a
unit step-function of $z^{'}$ can be turned into a delta function of
$z^{'}$.

**6.2 Explicit field expressions**

Writing everything explicitly, we have

> $B_{z} = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}}r^{'}J_{q}\left( r^{'},z^{'} \right)\int_{}^{}{dq^{'}}\frac{\left( r^{'} - r\text{cos}\left( q^{'} \right) \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + r^{'} + r^{2} - 2r^{'}r\text{cos}\left( q^{'} \right)}^{3}}}$
> . (6.12)

Denoting the denominator in eqn. (6.12)
by$\sqrt{{\text{.}\text{.}\text{.}\text{.}\text{.}}^{3}}$ we get for the
other fields

> $B_{r}\text{=-}\frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\frac{¶J_{q}\left( r^{'},z^{'} \right)}{¶z^{'}}}\int_{}^{}{dq^{'}}\frac{\text{cos}\left( q^{'} \right)}{\sqrt{\text{.}\text{.}\text{.}\text{.}\text{.}}}$
> , (6.13)
>
> $A_{q} = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}J_{q}\left( r^{'},z^{'} \right)}\int_{}^{}{dq^{'}}\frac{\text{cos}\left( q^{'} \right)}{\sqrt{\text{.}\text{.}\text{.}\text{.}\text{.}}}$
> , (6.14)
>
> $j = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\frac{¶S\left( r^{'},z^{'} \right)}{¶z^{'}}}\int_{}^{}{dq^{'}}\frac{1}{\sqrt{\text{.}\text{.}\text{.}\text{.}\text{.}}}$
> . (6.15)

It is clear that if we derive $B_{r}\text{=-¶}A_{q}/¶z$ from eqn. (6.14)
we will simply recover the result obtained by direct integration, eqn.
(6.13). We may also
derive$B_{z} = \frac{1}{r}\frac{\partial}{\partial r}\text{rA}_{\theta}$
using eqn. (6.14). This procedure is clearly correct, but it is not
obvious from inspection that the result is actually equivalent to the
direct integration, eqn. (6.12). We therefore have the alternative form:

> $B_{z} = \frac{1}{r}\frac{¶}{¶r}\text{rA}_{q} = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}J_{q}\left( r^{'},z^{'} \right)}$
> $\int_{}^{}{dq^{'}}\text{cos}\left( q^{'} \right)\left\lbrack \frac{1}{r\sqrt{\text{.}\text{.}\text{.}\text{.}\text{.}}} - \frac{r - r^{'}\text{cos}\left( q^{'} \right)}{\sqrt{{\text{.}\text{.}\text{.}\text{.}\text{.}}^{3}}} \right\rbrack$.
> (6.16)

It may be verified from numerical examples that the two formulas for
$B_{z}$ are the same. Mathematically, their equivalence can be reduced
to the identity

$\int_{}^{}{dq}\sqrt{1 - a\text{cos}(q)} = \left( 1 - a^{2} \right)\begin{matrix}
\int_{}^{}dq\frac{1}{\sqrt{{1 - a\text{cos}(q)}^{3}}}
\end{matrix}$ , (6.17)

which is also not obvious. However it can be shown that the two sides of
eqn. (6.17) have the same expansion coefficients multiplying all powers
of $a^{2}$. This identity will be used in section 7 and below in eqn.
(6.25d).

Not surprisingly, there is also a second integral expression for
$A_{\theta}$, which can be derived from eqns. (6.5) and (6.9):

$A_{\theta} = \text{rB}_{z} - r\frac{\partial A_{\theta}}{\partial r}$
(6.18)

> $= \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}\frac{J_{q}^{'}r^{'}r\text{sin}^{2}\left( q^{'} - q \right)}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|^{3}}}$
> (6.19)
>
> $= \frac{\mu_{0}}{4\pi}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}}r^{}\text{rJ}_{\theta}\left( r^{'},z^{'} \right)\int_{}^{}{d\theta^{'}}\frac{\text{sin}^{2}\left( \theta^{'} \right)}{\sqrt{{\text{.}\text{.}\text{.}\text{.}\text{.}}^{3}}}$
> . (6.20)

Similar to eqn. (6.17), the equivalence of this expression and eqn.
(6.14) implies a non-trival identity:

$\int_{}^{}{dq}\frac{\text{sin}^{2}(q)}{\sqrt{{1 - a\text{cos}(q)}^{3}}}\begin{matrix}

\end{matrix}$ . (6.21)

Numerical integration of the Green function solutions is greatly
simplified when the integrals over $\theta^{'}$ are worked in advance,
leaving only integrations over $r^{'}$ and $z^{'}$. This can be done
exactly using the complete elliptic integrals E and K, or in an accurate
approximation using analytical fits -- see section 7. First we isolate
the essential $\theta^{'}$ integrations; denote

> $a\left( r,z,r^{'},z^{'} \right)º\frac{2rr^{'}}{\left( z^{'} - z \right)^{2} + r^{2} + r^{}}$
> , (6.22)

which has range $0 \leq a \leq 1$, and define five functions:

$f_{0}(a) = \int_{}^{}\frac{dq}{2p}\sqrt{1 - a\text{cos}(q)}$ , (6.23a)

> $f_{1}(a) = \int_{}^{}\frac{dq}{2p}\frac{\text{cos}(q)}{\sqrt{1 - a\text{cos}(q)}}$
> , (6.23b)
>
> $f_{2}(a) = \int_{}^{}\frac{dq}{2p}\frac{1}{\sqrt{1 - a\text{cos}(q)}}$
> , (6.23c)

$f_{3}(a) = \int_{}^{}\frac{dq}{2p}\frac{\text{cos}(q)}{\sqrt{{1 - a\text{cos}(q)}^{3}}}$
, (6.23d)

$f_{4}(a) = \int_{}^{}\frac{dq}{2p}\frac{1}{\sqrt{{1 - a\text{cos}(q)}^{3}}}$
. (6.23e)

We express $f_{0}$ and $f_{2}$ in terms of E and K; the standard
definitions

\[6-1\] are

$E(x) = \int_{}^{}\text{dy}\sqrt{1 - x\begin{matrix}
\text{sin}^{2}
\end{matrix}(y)}$ , (6.24a)

$K(x) = \int_{}^{}\text{dy}\frac{1}{\sqrt{1 - x\begin{matrix}
\text{sin}^{2}
\end{matrix}(y)}}$ . (6.24b)

\[[Warning]{.underline}: these definitions of $E(x)$ and $K(x)$ are used
in Mathematica^®^ \[6-2\], but in some texts $x^{2}$ appears in the
integrands instead of $x$.\] Using these expressions and a change of
variable, we have

$f_{0}(a) = \frac{2\sqrt{1 + a}}{p}E\left( \frac{2a}{1 + a} \right)$ ,
(6.25a)

$f_{2}(a) = \frac{2}{p\sqrt{1 + a}}K\left( \frac{2a}{1 + a} \right)$ ,
(6.25b)

and the remaining$f_{i}(a)$ are found from $f_{0}$ and $f_{2}$ using the
relations

> $f_{1}(a) = \frac{1}{a}\left( f_{2}(a) - f_{0}(a) \right)$ , (6.25c)

$f_{4}(a) = \frac{f_{0}(a)}{1 - a^{2}}$ , (6.25d)

> $f_{3}(a) = \frac{1}{a}\left( \frac{f_{0}(a)}{1 - a^{2}} - f_{2}(a) \right)$
> . (6.25e)

The four Green function solutions, eqns. (6.12-6.15), are reduced to

$B_{z} = \frac{\mu_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{}J_{\theta}\left( r^{'},z^{'} \right)}\sqrt{\frac{a}{2rr^{'}}^{3}}\left\lbrack f_{4}(a) - \frac{r}{r^{'}}f_{3}(a) \right\rbrack$
, (6.26)

$B_{r} = \frac{- m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶J_{q}\left( r^{'},z^{'} \right)}{¶z^{'}}\sqrt{\frac{a}{2rr^{'}}}f_{1}(a)$
, (6.27)

$A_{q} = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}J_{q}\left( r^{'},z^{'} \right)}\sqrt{\frac{a}{2rr^{'}}}f_{1}(a)$
, (6.28)

$j = \frac{m_{0}}{2}$$\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\frac{¶S\left( r^{'},z^{'} \right)}{¶z^{'}}}\sqrt{\frac{a}{2rr^{'}}}f_{2}(a)$
. (6.29)

**References**

1.  Abramowitz, M. and Stegun, I.A., Handbook of Mathematical Functions,
    Dover, New York (1965).

> 6-2 Wolfram, S., The Mathematica^®^ Book, 3^rd^ edition, Cambridge
> University Press
>
> (1991).

**7. [Approximation of the Green Function Integrals ]{.underline}**

The complete elliptic integrals $E$and $K$, defined in section 6, may be
used to accurately evaluate the functions $f_{i}(a)$ that appear in the
general Green function solutions for $B_{z},B_{r},A_{q}$ and $j$.
However, the resulting formulas are complicated and not very revealing.
One might just as well evaluate the $f_{i}(a)$ numerically from their
defining integrals (see plots of $f_{i}(a)$ below). Except for $f_{0}$,
which is defined to help evaluate the other integrals, all of the
$f_{i}$ diverge as $a \rightarrow 1$, so their Taylor series expansions
are of limited value. However, since

$a \equiv \frac{2rr^{'}}{\left( z^{'} - z \right)^{2} + r^{} + r^{2}}$
(7.1)

is small when $r$ is either large or small, these expansions may be used
to evaluate both the near-axis field and far field. Expanding
$\sqrt{1 - a\text{cos}(q)}$ in the integrands and calculating the
resulting elementary integrals yields:

$f_{0}(a) = \int_{}^{}\frac{dq}{2p}\sqrt{1 - a\text{cos}(q)} =$ , (7.2a)

$f_{1}(a) = \int_{}^{}\frac{dq}{2p}\frac{\text{cos}(q)}{\sqrt{1 - a\text{cos}(q)}} = \frac{a}{4}$
, (7.2b)

$f_{2}(a) = \int_{}^{}\frac{dq}{2p}\frac{1}{\sqrt{1 - a\text{cos}(q)}} = \left( 1 + \frac{3}{\text{16}}a^{2} + \frac{\text{105}}{\text{1024}}a^{4} + \frac{\text{1155}}{\text{16384}}a^{6} + \cdot \right)$
, (7.2c)

$f_{3}(a) = \int_{}^{}\frac{dq}{2p}\frac{\text{cos}(q)}{\sqrt{{1 - a\text{cos}(q)}^{3}}} = \frac{3a}{4}\left( 1 + \frac{\text{35}}{\text{32}}a^{2} + \frac{\text{1155}}{\text{1024}}a^{4} + \cdot \right)$
, (7.2d)

$f_{4}(a) = \int_{}^{}\frac{dq}{2p}\frac{1}{\sqrt{{1 - a\text{cos}(q)}^{3}}} =$
. (7.2e)

These expansions suggest that $f_{1}$ and $f_{2}$diverge as Log
$\left( 1 - a^{2} \right)$, while $f_{3}$ and $f_{4}$ diverge as
$1/\left( 1 - a^{2} \right)$. This behavior is quantified below and
accurate approximate formulas are found that are valid through the
entire range $0 \leq a \leq 1$.

**8. [Thin Wire Layer -- Complete Solution]{.underline}**

**8.1 Semi-infinite layer**

In section 4 we derived the on-axis field for a thin semi-infinite wire
layer, which was then used to build up more complex layouts by shifting,
i.e. $z \rightarrow z - z_{0}$, and adding layers:

$J_{q}(r,z) = Sd(r - R)H( - z)$ , (8.1)

> $B_{0}(z) = \frac{m_{0}S}{2}\left( 1 - \frac{z}{\sqrt{z^{2} + R^{2}}} \right)$
> . (8.2)

Near-axis field components, including the lowest order fringe field
aberration, are

$\begin{array}{r}
B_{z} \approx B_{0} - \frac{B_{0^{\prime \prime}}r^{2}}{4} + \text{.}\text{.}\text{.} \\
\frac{\mu_{0}S}{2}\left\lbrack 1 - \frac{z}{\sqrt{z^{2} + R^{2}}} - \frac{3}{4}\frac{R^{2}\text{zr}^{2}}{\sqrt{{z^{2} + R^{2}}^{5}}} + \text{.}\text{.}\text{.} \right\rbrack\begin{matrix}

\end{matrix},
\end{array}$ (8.3a)

$\begin{array}{r}
B_{r}\text{=-}\frac{B_{0^{¢}}r}{2} + \frac{B_{0^{¢¢¢}}r^{3}}{\text{16}} - \text{.}\text{.}\text{.} \\
\frac{m_{0}S}{4}\left\lbrack \frac{R^{2}r}{\sqrt{{z^{2} + R^{2}}^{3}}} + \frac{3}{8}\begin{matrix}

\end{matrix}\frac{R^{2}\left( R^{2} - 4z^{2} \right)r^{3}}{\sqrt{{z^{2} + R^{2}}^{7}}} + \text{.}\text{.}\text{.} \right\rbrack\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (8.3b)

The near-axis results may be extended to any value of $r$ using the
Green function solutions. The resulting expressions involve only an
integration over $\theta$ because the integrations over $r^{'}$ and
$z^{'}$ are elementary. We have from eqns. (6.12) and (6.13):

$$\begin{array}{r}
B_{z} = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}}r^{'}\left\lbrack Sd\left( r^{'} - R \right)H\left( - z^{'} \right) \right\rbrack\begin{matrix}
\int_{}^{}{dq}r^{'} - r\frac{\begin{matrix}
\text{cos}(q)
\end{matrix}}{\sqrt{{\left( z^{'} - z \right)^{2} + r^{} + r^{2} - 2r^{'}r\begin{matrix}
\text{cos}(q)
\end{matrix}}^{3}}} & & \begin{matrix}
 & & & & 
\end{matrix} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p}\int_{}^{}{dz^{'}}\frac{R^{2} - \text{Rr}\begin{matrix}

\end{matrix}\text{cos}(q)}{\sqrt{{\left( z^{'} - z \right)^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}^{3}}} & 
\end{matrix} \\

\end{array}$$

$\begin{matrix}
 & & & 
\end{matrix} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p}\left( \frac{R^{2} - \text{Rr}\begin{matrix}

\end{matrix}\text{cos}q}{R^{2} + r^{2} - 2\text{Rr}\text{cos}q} \right)\left\lbrack 1 - \frac{z}{\sqrt{z^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}} \right\rbrack\begin{matrix}

\end{matrix},$ (8.4a)

$B_{r} = \frac{- m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\left\lbrack \frac{¶}{¶z^{'}}Sd\left( r^{'} - R \right)H\left( - z^{'} \right) \right\rbrack}\int_{}^{}{dq\frac{\text{cos}(q)}{\sqrt{\left( z^{'} - z \right)^{2} + r^{} + r^{2} - 2r^{'}r\text{cos}(q)}}}$

$\begin{matrix}
 & & & & 
\end{matrix} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p}\frac{R\begin{matrix}

\end{matrix}\text{cos}(q)}{\sqrt{z^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}}\begin{matrix}

\end{matrix}\text{.}$ (8.4b)

The first integral in eqn. (8.4a) for $B_{z}$ is simply

> $\int_{}^{}\frac{dq}{2p}\begin{matrix}
>
> \end{matrix}$ (8.5)

so we find the expected limits

$B_{z} \rightarrow \begin{array}{r}
\left\{ m_{0}\text{SH}(R - r) \middle| \right.\ 
\end{array}$ $\begin{array}{r}
z \rightarrow \text{-¥} \\
z \rightarrow \text{+¥}\begin{matrix}

\end{matrix},
\end{array}$ (8.6)

while at $z = 0$ there is the simple result

$B_{z} = \frac{m_{0}S}{2}H(R - r)$ . (8.7)

Equation (8.2) for $B_{0}(z)$ is recovered from eqn. (8.4a) by setting
$r = 0$. For general $(r,z)$ $B_{z}$ does not appear to be expressible
in terms of commonly used functions, but if we insert $r = R$, then we
obtain simply

$B_{z}(R,z) = \frac{m_{0}S}{4}\left\lbrack 1 - \frac{\text{zf}_{2}\left( 2R^{2}/\left( z^{2} + 2R^{2} \right) \right)}{\sqrt{z^{2} + 2R^{2}}} \right\rbrack$
. (8.8)

For $z < 0$ this is the longitudinal field along the
[center]{.underline} of the thin wire layer. The longitudinal field
immediately inside or outside the layer (denoted $r = R \pm$) is
therefore

$B_{z}(R \pm ,z < 0) = \frac{m_{0}S}{4}\left\lbrack 1 - \frac{\text{zf}_{2}\left( 2R^{2}/\left( z^{2} + R^{2} \right) \right)}{\sqrt{z^{2} + 2R^{2}}} \right\rbrack m\frac{m_{0}S}{2}$
, (8.9)

with $f_{2}$ defined in section (7):

$f_{2}(a) = \int_{}^{}\frac{dq}{2p}\frac{1}{\sqrt{1 - a\text{cos}(q)}}$
.

**8.2 Wire end divergence**

From equations (8.8) and (8.9) follows that very close to the wire layer
ends $(r = R \pm ,z = 0 \pm )$, $B_{z}$ is discontinuous as:

  ---------------------- ------------------------ ------------------------
  $r$                    $z$                      $B_{z}/m_{0}S$

  $R -$                  $0 -$                    $34$

  $R -$                  $0 +$                    $14$

  $R +$                  $0 -$                    $- 14$

  $R +$                  $0 +$                    $14$
  ---------------------- ------------------------ ------------------------

This seemingly peculiar behavior is easily seen by a numerical
evaluation of expression (8.4a) near the wire layer end. It can also be
deduced from the behavior of $B_{r}$ (see below) using the Maxwell
equation $\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$.

Expression (8.4b) for $B_{r}$ can be written in the compact form

$B_{r} = \frac{m_{0}S}{2}\begin{matrix}

\end{matrix}$ , (8.10)

with $bº2\text{Rr}/\left( z^{2} + r^{2} + R^{2} \right)$, and $f_{1}(b)$
defined in section 7 as:

$f_{1}(b) = \int_{}^{}\frac{dq}{2p}\frac{\text{cos}(q)}{\sqrt{1 - b\text{cos}(q)}}$
. (8.11)

Recall that as $b \rightarrow 1$,

$f_{1}(b) \rightarrow \frac{C}{4b}\left\lbrack D - \text{Log}\left( 1 - b^{2} \right) \right\rbrack$
, (8.12)

where $C = 2\frac{\sqrt{2}}{p}$ and
$D = \left( 6\text{Log}(2) - 4 \right)$. It is apparent that $B_{r}$ is
weakly divergent at the wire layer end $(r = R,z = 0)$. Explicitly,
using eqns. (8.10, 8.12),

$B_{r} \rightarrow \frac{m_{0}S}{4p}\left\{ D - \text{Log}\left\lbrack \frac{z^{2} + (r - R)^{2}}{R^{2}} \right\rbrack \right\}$
. (8.13)

Although $B_{r}$ becomes infinite for an infinitesimally thin current
layer, it is finite for a layer of finite thickness. In fact in a thin
layer $B_{r}$ exceeds the scale field strength $\left( m_{0}S \right)$
only when $\sqrt{z^{2} + (r - R)^{2}} \leq R/\text{400}$. For thickness
greater than $R/\text{20}$ the divergent behavior of the thin layer is
not representative. Interestingly, the logarithmic divergence does have
an effect in permanent magnets, which have essentially infinitesimal
surface current layer thickness, which causes a partial demagnetization
at magnet corners due to the non-ideal magnetization properties of the
material.

For $\overrightarrow{|r|}\text{>>}R$, eqns. (8.4a,b) may be expanded in
inverse powers of $\overrightarrow{|r|}$ to yield in leading order

$\overrightarrow{B}»m_{0}\text{SH}(R - r)H( - z){\hat{e}}_{z} + \frac{m_{0}\text{SR}^{2}\overrightarrow{r}}{4\left| \overrightarrow{r} \right|^{3}}$
. (8.14)

This is just a uniform direct field that terminates at $z = 0$ plus a
residual "magnetic monopole field" with net flux ${pR}^{2}m_{0}S$
centered at $\overrightarrow{r} = 0$. A generalization of this result is
presented in section 11.

**8.3 Finite length lense**

For a solenoid of finite length $\ell$ and wire radius $R$, centered at
$z = 0$, equations (8.4a,b) give by translation and superposition

$\begin{array}{r}
B_{z} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p} \\

\end{array}$ (8.15a)

$\begin{array}{r}
B_{r} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p} \\
 \cdot \begin{matrix}

\end{matrix}\left\lbrack \frac{R\begin{matrix}

\end{matrix}\text{cos}(q)}{\sqrt{(z - l/2)^{2} + r^{2} + R^{2} - 2\text{Rr}\begin{matrix}

\end{matrix}\text{cos}(q)}} - \frac{R\begin{matrix}

\end{matrix}\text{cos}(q)}{\sqrt{(z + l/2)^{2} + r^{2} + R^{2} - 2\text{Rr}\text{cos}(q)}} \right\rbrack\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (8.15b)

Plots of $B_{z}$ and $B_{r}$ for several values of $r$ are displayed
below for a lense of length $l = 4\text{.}\text{0m}$ and wire radius
$R = 1\text{.}\text{0m}$.

![](media/image5.pct){width="5.993055555555555in" height="6.30625in"}

The vector potential of a finite length solenoid can also be expressed
by a single

integral over $q$ using the Green function (6.14):

> $\begin{array}{r}
> A_{q} = \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}\frac{J_{q^{¢}}\text{cos}\left( q^{'} - q \right)}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}} \\
> \begin{matrix}
>  & & 
> \end{matrix} = \frac{m_{0}\text{SR}}{4p}\int_{}^{}{dq\begin{matrix}
>
> \end{matrix}\text{cos}(q)\int_{}^{}{dz^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2} + R^{2} - 2\text{Rr}\text{cos}(q)}}}
> \end{array}$

$\begin{matrix}
 & & 
\end{matrix} = \frac{m_{0}\text{SR}}{2}\int_{}^{}{\frac{dq}{2p}\begin{matrix}

\end{matrix}\text{cos}(q)\text{Log}\left\lbrack \frac{\sqrt{(z - l/2)^{2} + R^{2} + r^{2} - 2\text{Rr}\begin{matrix}

\end{matrix}\text{cos}(q)} - (z - l/2)}{\sqrt{(z + l/2)^{2} + R^{2} + r^{2} - 2\text{Rr}\begin{matrix}

\end{matrix}\text{cos}(q)} - (z + l/2)} \right\rbrack\begin{matrix}

\end{matrix}\text{.}}$ (8.16)

Equation (6.20) gives the alternative, equivalent form

$\begin{array}{r}
A_{q} = \frac{m_{0}\text{SR}^{2}r}{2}\int_{}^{}\frac{dq}{2p}\frac{\text{sin}^{2}(q)}{R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}\begin{matrix}

\end{matrix} \cdot \left\lbrack \frac{(z + l/2)}{\sqrt{(z + l/2)^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}} - \frac{(z - l/2)}{\sqrt{(z - l/2)^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}} \right\rbrack\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (8.17)

**8.4 Stored energy**

The magnetic field energy of a thin-layer lens is easily evaluated. In
the absence of

permeable materials we have in general,

$W = \int_{}^{}d^{3}r\frac{\left| \overrightarrow{B} \right|^{2}}{2m_{0}} = \int_{}^{}{d^{3}r\frac{\overrightarrow{B} \cdot \overrightarrow{Ñ}x\overrightarrow{A}}{2m_{0}}} = \int_{}^{}d^{3}r\frac{\overrightarrow{A} \cdot \overrightarrow{Ñ}x\overrightarrow{B}}{2m_{0}} = \int_{}^{}d^{3}r\frac{\overrightarrow{A} \cdot \overrightarrow{J}}{2}$
. (8.18)

For $J_{q} = S(z)d(r - R)$, eqn. (8.18) gives

$W = \pi R\int_{}^{}{\text{dz}\begin{matrix}

\end{matrix}S(z)A_{\theta}(R,z)}$. (8.19)

Inserting expression (8.17) for $A_{q}$ and working the elementary
integrations in $z$ gives for a lens of length $l$:

$W = \frac{{pR}^{2}l\left( m_{0}S \right)^{2}}{2m_{0}}f(x)$ , (8.20)

where $x = l^{2}/2R^{2}$ and

$f(x) = \frac{1}{\sqrt{x}}\int_{}^{}\frac{dq}{2p}\left( 1 + \text{cos}(q) \right)\left\lbrack \sqrt{x + 1 - \text{cos}(q)} - \sqrt{1 - \text{cos}(q)} \right\rbrack$.
(8.21)

This function can be expressed using the previously defined $f_{0}$ and
$f_{3}$ (section 7):

$f(x) = \sqrt{\frac{x + 1}{x}}\left\lbrack f_{0}\left( \frac{1}{1 + x} \right) - \frac{1}{3}\left( 1 - \frac{1}{(1 + x)^{2}} \right)f_{3}\left( \frac{1}{1 + x} \right) \right\rbrack - \frac{2}{3}\frac{C}{\sqrt{x}}$
, (8.22)

which slowly approaches unity for large $x$:

$f(x) \rightarrow 1 - \frac{2}{3}\frac{C}{\sqrt{x}} + \frac{1}{4x}$ .
(8.23)

For very large $x\text{>>}1$, i.e. $l\text{>>}\sqrt{2}R,\begin{matrix}
f(x) \rightarrow 1
\end{matrix}$ and expression (8.20) is just the energy one gets using
only the direct field of the magnet bore. For $x > \text{.}2$, which
covers all cases of practical interest,

$f(x)»\frac{1}{1 + \frac{\text{.}\text{63}}{\sqrt{x}}}$ , (8.24)

with error less that 1%, so we have in general for the thin layer lens

$W»\frac{{pR}^{2}l\left( m_{0}S \right)^{2}}{2m_{0}}\begin{matrix}
\frac{1}{1 + \text{.}\text{89}R/l}
\end{matrix}$ . (8.25)

It is somewhat surprising that the actual value of $W$for the
finite-length lens is smaller than the approximate value found using
only the direct field.

Suppose the current layer is made of $n$ turns of wire carrying current
$I$. Then $S = /\ell$ and the magnet's inductance is

$L = \frac{2W}{I^{2}} \approx \frac{\mu_{0}n^{2}\pi R^{2}/\ell}{1 + \text{.}\text{89}R/\ell}$
. (8.26)

![](media/image6.pct){width="5.986111111111111in"
height="8.097222222222221in"}

**9. [Thick Wire Layer ]{.underline}**

**9.1 Uniform thick layer**

While the thin wire layer formulas are sufficient for most beam dynamics
applications, thick layer results are needed for precise near-axis
fields, external fields, field in the wire, and mechanical stress. We
first treat the case of radially uniform current density

$J_{q} = \frac{S(z)}{R_{2} - R_{1}}\begin{matrix}
 & & & & 
\end{matrix}\left( R_{1} < r < R_{2} \right)$ . (9.1)

As in section 4, it is usually sufficient to derive the field for a
semi-infinite thick layer, $S(z) = \text{SH}( - z)$, and any system of
lenses may be built up from this by offsets
$\left( z \rightarrow z - z_{0} \right)$ and super-position.

Recall from section 8 the semi-infinite thin layer field
$\left\lbrack J_{q} = \text{SH}( - z)d(r - R) \right\rbrack$:

$B_{z}^{\text{thin}} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p}\frac{R^{2} - \text{Rr}\text{cos}(q)}{R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}\left\lbrack 1 - \frac{z}{\sqrt{z^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}} \right\rbrack$
, (9.2a)

$B_{r}^{\text{thin}} = \frac{m_{0}S}{2}\int_{}^{}\frac{dq}{2p}\frac{R\text{cos}(q)}{\sqrt{z^{2} + R^{2} + r^{2} - 2\text{Rr}\text{cos}(q)}}$.
(9.2b)

For the thick, uniform, semi-infinite layer we simply compute a radial
superposition of these thin layer fields:

$\left( B_{z},B_{r} \right)^{\text{thick}} = \int_{}^{}\frac{\text{dR}}{R_{2} - R_{1}}\left( B_{z},B_{r} \right)^{\text{thin}}$
. (9.3)

The on-axis field of the semi-infinite thick layer is

$\begin{array}{r}
B_{0}^{\text{thick}}(z) = \int_{}^{}\frac{\text{dR}}{R_{2} - R_{1}}\begin{matrix}
\frac{m_{0}S}{2}\left( 1 - \frac{z}{\sqrt{z^{2} + R^{2}}} \right)
\end{matrix}\begin{matrix}
 & & & & 
\end{matrix} = \frac{m_{0}S}{2}\left\lbrack 1 - \frac{z}{R_{2} - R_{1}}\text{Log}\left( \frac{R_{2} + \sqrt{z^{2} + R_{2}^{2}}}{R_{1} + \sqrt{z^{2} + R_{1}^{2}}} \right) \right\rbrack\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (9.4)

As mentioned in section 4, for a thick-layer lens of length $l$ centered
at $z = 0$, the on-axis field is

$\begin{array}{r}
B_{0}^{\text{Lens}}(z) = \frac{m_{0}S}{2\left( R_{2} - R_{1} \right)}\left\{ \left( \frac{l}{2} + z \right) \right.\ \text{Log}\left\lbrack \frac{R_{2} + \sqrt{\left( \frac{l}{2} + z \right)^{2} + R_{2}^{2}}}{R_{1} + \sqrt{\left( \frac{l}{2} + z \right)^{2} + R_{1}^{2}}} \right\rbrack \\
\begin{matrix}
 & & & & 
\end{matrix} + \left( \frac{l}{2} - z \right)\text{Log}\left\lbrack \frac{R_{2} + \sqrt{\left( \frac{l}{2} - z \right)^{2} + R_{2}^{2}}}{R_{1} + \sqrt{\left( \frac{l}{2} - z \right)^{2} + R_{1}^{2}}} \right\rbrack \\

\end{array}$ (9.5)

Near-axis values of $B_{z}$and $B_{r}$ can then be obtained using these
results and the radial expansions of section 3. For the field at the
magnet center $(z = 0)$, eqn. (9.5) reduces to

$B_{0}^{\text{Lens}}(0) = \frac{m_{0}Sl}{2\left( R_{2} - R_{1} \right)}\text{Log}\left( \frac{R_{2} + \sqrt{l^{2}/4 + R_{2}^{2}}}{R_{1} + \sqrt{l^{2}/4 + R_{1}^{2}}} \right)$
. (9.6)

The semi-infinite, thick-layer formulas for $B_{z}$and $B_{r}$at general
$(r,z)$ can also be written in terms of integrations over $\theta$. The
integrals over $R$ in eqn. (9.3) can actually be evaluated in terms of
elementary functions, although the results are somewhat messy. First we
define

$u(r,q) = R_{1} - r\begin{matrix}

\end{matrix}$ . (9.7a,b)

Then, using the integrals given at the end of section 3 and omitting
some intermediate steps, we get from eqns. (9.2, 9.3):

$B_{r}^{\text{thick}} = \frac{m_{0}S}{2\left( R_{2} - R_{1} \right)}\int_{}^{}\frac{dq}{2p}\text{cos}(q)\begin{array}{r}
\left\{ \middle| \right.\ \left\{ \sqrt{v^{2} + z^{2} + r^{2}\text{sin}^{2}(q)} \middle| \right.\ 
\end{array}$

> $$$\left. \  - \sqrt{u^{2} + z^{2} + r^{2}\text{sin}^{2}(q)} + r\text{cos}(q)\text{Log}\left\lbrack \frac{v + \sqrt{v^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}}{u + \sqrt{u^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}} \right\rbrack \right\}$
> , (9.8)

$$B_{z}^{\text{thick}} = \frac{m_{0}S(r)}{2} - \frac{m_{0}\text{Sz}}{2\left( R_{2} - R_{1} \right)}\int_{}^{}\frac{dq}{2p}\left\{ \text{Log}\left\lbrack \frac{v + \sqrt{v^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}}{u + \sqrt{u^{2} + z^{2} + r^{2} + \text{sin}^{2}(q)}} \right\rbrack \right.\ $$

$\begin{matrix}
 & 
\end{matrix} - \left( \frac{r\begin{matrix}

\end{matrix}\text{cos}(q)}{z} \right)\left\lbrack \text{tanh}^{- 1}\left( \frac{\sqrt{v^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}}{z} \right) - \text{tanh}^{- 1}\left( \frac{\sqrt{u^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}}{z} \right) \right\rbrack$

$\begin{matrix}
 & 
\end{matrix} - \left( \frac{r\begin{matrix}

\end{matrix}\text{sin}(q)}{z} \right)\left. \ \left\lbrack \text{tan}^{- 1}\left( \frac{z}{r\text{sin}(q)}\frac{v}{\sqrt{v^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}} \right) - \text{tan}^{- 1}\left( \frac{z}{r\text{sin}(q)}\frac{u}{\sqrt{u^{2} + z^{2} + r^{2}\text{sin}^{2}(q)}} \right) \right\rbrack \right\}\text{.}$
(9.9)

Here we have used

$S(r) = \begin{array}{r}
\left\{ S \middle| \right.\ \left\{ S \middle| \right.\ 
\end{array}$ (9.10)

**9.2 Approximation by a thin layer**

It was asserted in section 4 that the on-axis field of a thin layer
magnet, with small error could replace the on-axis field of a uniform
thick layer magnet. There are many ways to select the three parameters
(S, R, $l$) of the thin layer. Here we use the same S and $l$ as the
thick layer and set R to give the correct net magnetic dipole moment by
requiring

$\text{SR}^{2} = \int_{}^{}{\text{drJ}_{q}r^{2}} = \int_{}^{}\text{dr}\frac{S}{R_{2} - R_{1}}r^{2} = \frac{S}{3}\left( R_{2}^{2} + R_{2}R_{1} + R_{1}^{2} \right)$
. (9.11)

Defining a thickness parameter

$Dº\left( R_{2} - R_{1} \right)/R$ , (9.12)

we have from eqn. (9.11)

$R_{1} = R\left( \sqrt{1 - \frac{\Delta^{2}}{\text{12}}} - \frac{\Delta}{2} \right)$,
$R_{2} = R\left( \sqrt{1 - \frac{D^{2}}{\text{12}}} + \frac{D}{2} \right)$
. (9.13a,b)

We require that $D < \sqrt{3}$ to prevent $R_{1} < 0$. A somewhat
different formula for $R$ is obtained if we equate the focal lengths of
the magnets
$\left( f^{- 1}\text{\textasciitilde}\int_{}^{}\text{dzB}_{0}^{2}(z) \right)$.
The following table gives the maximum ratio of on-axis fields for the
thin layer eqn. (4.6), relative to the equivalent thick layer eqn.
(9.5), for various values of $D$ and $l/R$. Except for very small $l/R$
this maximum occurs in the fringe at
$z» \pm \left( l/2 + \text{.}8R \right)$.

  ---------- --------------------------------------------- ------------------- ------------------- --------------------
  $$l/R$$    $$\begin{array}{r}                            $$D = \text{.}2$$   $$D = \text{.}5$$   $$D = 1\text{.}0$$
             B_{0}^{\text{thin}}/B_{0}^{\text{thick}} \\                                           
             D = \text{.}1                                                                         
             \end{array}$$                                                                         

  1          1.00094                                       1.0038              1.024               1.101

  2          1.00088                                       1.0035              1.023               1.094

  4          1.00082                                       1.0033              1.021               1.087

  8          1.00079                                       1.0032              1.020               1.084

  16         1.00078                                       1.0031              1.020               1.083
  ---------- --------------------------------------------- ------------------- ------------------- --------------------

**9.3 General thick layer**

Finally we treat the thick wire layer with specified ends and general
current density profile $J_{q}(r)\text{=-}\text{dS}(r)/\text{dr}$
between the ends. The Green function solutions of section 6 can be used
for this, as can simple averages of eqns. (9.2). However, more useful
expressions can be derived using the Bessel transform. First, the
semi-infinite layer $J_{q} = J_{q}(r)H( - z)$ has
$B_{z}\left( r,z \rightarrow \text{-¥} \right) = m_{0}S(r)$, and by
symmetry

$B_{z}(r,z = 0) = \frac{m_{0}S(r)}{2}$ . (9.14)

For $z > 0,$

$B_{z}\text{\textasciitilde}J_{0}\left( \text{kr} \right)e^{- \text{kz}}$
, (9.15)

with $k > 0$ is a valid solution of $Ñ^{2}B_{z} = 0$, except for the
boundary condition (9.14). The correct solution for $z³0$ is given by
the Bessel transform:

$B_{z}(r,z³0) = \frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}
\text{kJ}_{0}\left( \text{kr} \right)e^{- \text{kz}}\overset{\sim}{S}
\end{matrix}(k)}$, (9.16)

$\overset{\sim}{S}(k) = \int_{}^{}{\text{dr}\begin{matrix}

\end{matrix}\text{rJ}_{0}\left( \text{kr} \right)S(r)}$ . (9.17)

We denote the n^th^ ordinary Bessel function by $J_{n}(x)$. A similar
calculation gives

$B_{z}(r,z \leq 0) = m_{0}S(r) - \frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}\text{kJ}_{0}\left( \text{kr} \right)e^{+ \text{kz}}\overset{\sim}{S}(k)}$
. (9.18)

Combining these results we have for the entire semi-infinite layer:

$B_{z} = \mu_{0}S(r)H( - z) + \frac{z}{|z|}\frac{\mu_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}\text{kJ}_{0}\left( \text{kr} \right)e^{- k|z|}\overset{\sim}{S}(k)}$
. (9.19)

$\overrightarrow{B}$ can also be derived from the modified scalar
potential (see section 5):

$B_{z}(r,z) = m_{0}S(r)H( - z) + ¶j/¶z$ . (9.20)

Equation (9.19) then gives for $j$ and $B_{r}$:

$j(r,z)\text{=-}\frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}J_{0}\left( \text{kr} \right)e^{- k|z|}\overset{\sim}{S}(k)}$
, (9.21)

$B_{r}(r,z) = \frac{\text{¶j}}{¶r} = \frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}\text{kJ}_{1}\left( \text{kr} \right)e^{- k|z|}\overset{\sim}{S}(k)}$
. (9.22)

Using superposition, a lens of length $l$, centered at $z = 0$ is seen
to have fields

$\begin{array}{r}
B_{z}(r,z) = m_{0}S(r)\left\lbrack H( - z + l/2) - H( - z - l/2) \right\rbrack \\
 + \frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}\text{kJ}_{0}\left( \text{kr} \right)\overset{\text{\textasciitilde}}{S}(k)\left\lbrack \frac{(z - l/2)}{|z - l/2|}e^{- k|z - l/2|}\begin{matrix}
 - \frac{(z + l/2)}{|z + l/2|}e^{- k|z + l/2|}
\end{matrix} \right\rbrack\begin{matrix}

\end{matrix},}
\end{array}$ (9.23)

$B_{r}(r,z) = \frac{m_{0}}{2}\int_{}^{}{\text{dk}\begin{matrix}

\end{matrix}\text{kJ}_{1}\left( \text{kr} \right)\overset{\text{\textasciitilde}}{S}(k)\left\lbrack e^{- k|z - l/2|} - e^{- k|z + l/2|} \right\rbrack}$
. (9.24)

It remains to evaluate $\overset{\sim}{S}(k)$. First we integrate by
parts in eqn. (9.17) using $J_{q} = - \text{dS}(r)/\text{dr}$:

$\overset{\sim}{S}(k) = \int_{}^{}{\text{drS}(r)}\frac{¶}{¶r}\frac{\text{rJ}_{1}\left( \text{kr} \right)}{k} = \int_{}^{}{\text{drJ}_{q}(r)\frac{\text{rJ}_{1}\left( \text{kr} \right)}{k}}$
. (9.25)

[Case 1]{.underline} (thin layer):

$J_{q}(r) = Sd(r - R)$ , (9.26)

$\overset{\sim}{S} = \frac{\text{RSJ}_{1}\left( \text{kR} \right)}{k}$ .
(9.27)

[Case 2]{.underline} ($J_{q}$ decreasing with $r$):

$J_{q}(r) = \frac{S}{r\begin{matrix}
\text{Log}\left( R_{2}/R_{1} \right)
\end{matrix}}\begin{matrix}
 & & 
\end{matrix}R_{1} < r < R_{2}$ , (9.28)

$\overset{\sim}{S}(k) = \frac{S}{k^{2}\begin{matrix}
\text{Log}\left( R_{2}/R_{1} \right)
\end{matrix}}\left\lbrack J_{0}\left( \text{kR}_{1} \right) - J_{0}\left( \text{kR}_{2} \right) \right\rbrack$
. (9.29)

[Case 3]{.underline} (uniform thick layer):

$J_{q} = \frac{S}{R_{2} - R_{1}}\begin{matrix}
 & & 
\end{matrix}R_{1} < r < R_{2}$ , (9.30)

$\overset{\sim}{S}(k) = \frac{S}{k^{3}\left( R_{2} - R_{1} \right)}\int_{\text{kR}_{1}}^{\text{kR}_{2}}\text{dy}\text{yJ}_{1}(y)\begin{matrix}

\end{matrix}\text{.}$ (9.31)

The integration over $y$, although not elementary, can be expressed as a
hypergeometric function:

$\int_{}^{}\text{dy}\text{yJ}_{1}(y) = \frac{x^{3}}{6}\begin{matrix}
F_{2}
\end{matrix}\left( \frac{3}{2};2,\frac{5}{2}; - x^{2}/4 \right)$ .
(9.32)

In Mathematica^®^ notation

$F_{2}(a;b,c;z) = \text{Hypergeometic PFQ}\left\lbrack \left\{ a \right\},\left\{ b,c \right\},z \right\rbrack$
(9.33)

$= 1 + \frac{a}{\text{bc}}\frac{z}{1!} + \frac{a(a + 1)}{b(b + 1)c(c + 1)}\frac{z^{2}}{2!} + \text{.}\text{.}\text{.}$
. (9.34)

Equation (9.31) becomes

$\overset{\sim}{S}(k) = \frac{S}{6\left( R_{2} - R_{1} \right)}$ .
(9.35)

**10. [Periodic Thin Wire Layer]{.underline}**

The field of an infinite, periodic wire layer is of interest for beam
transport. Its external field is very different from that of an isolated
magnet, particularly at large $r$, where it falls off exponentially
rather than as $r^{- 3}$. The return flux is diluted to zero intensity
in this case.

Let $\ell$denote an individual solenoid's wire layer length and $g$
denote gap lengths between adjacent wire layers. Assuming the solenoids
all have the same polarity, the lattice period length is$P = \ell + g$.
For one magnet of the lattice we set

$J_{q} = Sd(r - R)\begin{matrix}
 & & & 
\end{matrix}\text{for}\begin{matrix}

\end{matrix} - \frac{l}{2} < z < \frac{l}{2}$ , (10.1)

and extend this layer periodically. Then the Fourier expansion of $S(z)$
with period $P$ is

$S(z) = \frac{Sl}{P} + \sum_{n = 1}^{¥}S_{n}\text{cos}\left( \frac{2p\text{nz}}{P} \right)$
, (10.2)

with

$S_{n} = \frac{2S}{\pi n}\text{sin}\left( \frac{\pi n\ell}{P} \right)$ .
(10.3)

The longitudinal field can be similarly expanded:

$B_{z} = m_{0}S\frac{l}{P}H(R - r) + \sum_{n = 1}^{¥}{B_{n}(r)\text{cos}\left( \frac{2p\text{nz}}{P} \right)}$
, (10.4)

with

$B_{n}(r) = \left( m_{0}S \right)\text{sin}\left( \frac{pnl}{P} \right)\left( \frac{4R}{P} \right)K_{1}\left( \frac{2p\text{nR}}{P} \right)I_{0}\left( \frac{2p\text{nr}}{P} \right)\begin{matrix}
 & 
\end{matrix}0 < r < R\begin{matrix}

\end{matrix},$ (10.5a)

$B_{n}(r)\text{=-}\left( m_{0}S \right)\text{sin}\left( \frac{pnl}{P} \right)\left( \frac{4R}{P} \right)I_{1}\left( \frac{2p\text{nR}}{P} \right)K_{0}\left( \frac{2p\text{nr}}{P} \right)\begin{matrix}
 & 
\end{matrix}R < r\text{<¥}\begin{matrix}

\end{matrix},$ (10.5b)

where $I_{n}$ and $K_{n}$ are modified Bessel functions of the first and
second kind. The functions $B_{n}(r)$ satisfy

$\frac{1}{r}\frac{\partial}{\partial r}r\frac{\partial B_{n}}{\partial r} - \left( \frac{2\pi n}{P} \right)^{2}B_{n} = - \mu_{0}S_{n}\frac{1}{r}\frac{\partial}{\partial r}r\delta(r - R)$
. (10.6)

We have used the jump conditions for $B_{n}(r)$ at $R$, which may be
derived from eqn. (10.6):

$B_{n}(R + ) - B_{n}(R - ) = - \mu_{0}S_{n}$ , (10.7a)

$\frac{\partial B_{n}}{\partial r}(R + ) - \frac{\partial B_{n}}{\partial r}(R - ) = 0$,
(10.7b)

and the Wronskian identity,

$I_{0}(x)K_{1}(x) + I_{1}(x)K_{0}(x) = \frac{1}{x}$ .

In similar fashion we find

> $B_{r} = \sum_{n = 1}^{\infty}{C_{n}(r)}\text{sin}\left( \frac{2\pi\text{nz}}{P} \right)$
> , (10.8)

$C_{n}(r) = \left( m_{0}S \right)\text{sin}\left( \frac{pnl}{P} \right)\left( \frac{4R}{P} \right)K_{1}\left( \frac{2p\text{nR}}{P} \right)I_{1}\left( \frac{2p\text{nr}}{P} \right)\begin{matrix}
 & & 
\end{matrix}0 < r < R\begin{matrix}

\end{matrix},$ (10.9a)

$C_{n}(r) = \left( m_{0}S \right)\text{sin}\left( \frac{pnl}{P} \right)\left( \frac{4R}{P} \right)I_{1}\left( \frac{2p\text{nR}}{P} \right)K_{1}\left( \frac{2p\text{nr}}{P} \right)\begin{matrix}
 & & 
\end{matrix}R < r\text{<¥}\begin{matrix}

\end{matrix}\text{.}$ (10.9b)

This solves $\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$ with
$B_{z}$ given above. The components $C_{n}(r)$ are continuous at
$r = R$.

At large $r$the asymptotic expansion of $K_{0}(2pnrP)$ yields

$B_{n}\text{\textasciitilde}\frac{\text{exp}( - 2\pi nrP)}{\sqrt{r}}$ ,
(10.10)

showing the asserted exponential decrease of field strength outside the
wire. However, the calculated net flux in this example is
$\mu_{0}{S\pi R}^{2}\ell/P$, so the return flux must be regarded as
being diluted to zero magnitude. Plots of $B_{z}$and $B_{r}$ at various
$r$ are presented below for the parameters:

$m_{0}S\frac{l}{P} = \text{10}\text{.}0T,\begin{matrix}

\end{matrix}$

A total of 200 terms are used in the Fourier expansion. Note that while
$B_{z}$ reaches \~7.0T in the gaps, it drops to less than 2.0T in the
volume behind the wire where an induction core could be located. This
calculation is modified in section 18 with the addition of a permeable
yoke to reduce the external field.

> ![](media/image7.png){width="3.1798611111111112in"
> height="1.9638888888888888in"}

$B_{z}\begin{matrix}

\end{matrix}$

![](media/image8.png){width="3.0555555555555554in"
height="1.8756944444444446in"}

$B_{z}\begin{matrix}

\end{matrix}$

![](media/image9.png){width="3.1798611111111112in"
height="1.9590277777777778in"}

$B_{z}\begin{matrix}

\end{matrix}$

**11. [Far Field and External Field]{.underline}**

Because a solenoid system may have a large dipole moment, its field can
interfere with nearby experiments and the operation of accelerator
components. A short magnetic bend, which resembles a solenoid turned on
its side, also has a net dipole moment, but a permeable yoke is usually
used to confine its return flux. This technique is not as effective for
transport solenoids because the yoke must not block the vacuum bore.
Also, if annular yokes are used in an induction linac they may
considerably increase the diameter of adjacent induction cores and
thereby increase the overall cost. It is possible to greatly reduce the
far field by alternating the polarity of solenoid lenses, but this
reduces their effectiveness for beam transport and actually increases
the external field (just outside the wire). We are therefore motivated
to evaluate the far and external fields in the absence of return flux
confinement.

The far field and external field can be approximated using two different
series expansions of the global potential $j$. Both expansions may have
contributions from several or many individual magnets, so it may be
misleading to consider only one magnet at a time. Recall from section 5,

> $\overrightarrow{B} = m_{0}S(r,z){\hat{e}}_{z} + \overrightarrow{Ñ}j$,
> (11.1)

where S(r,z) is the cumulative surface current density:

$S(r,z) = \int_{}^{}{dr^{'}}J_{\theta}\left( r^{'},z \right)$ . (11.2)

The Green function solution for $j$ is, from section (6):

$j = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\int_{}^{}{dq^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{} + r^{2} - 2r^{'}r\text{cos}\left( q^{'} - q \right)}}\begin{matrix}
\text{.}
\end{matrix}$ (11.3)

**11.1 Far field expansion**

Consider first the far field expansion in inverse powers of
$\left| \overrightarrow{r} \right| = \left( r^{2} + z^{2} \right)^{12}$,
with r and z both considered to be large of order
$\left| \overrightarrow{r} \right|$, which is larger than
$\left| \overrightarrow{r^{'}} \right|$. We could laboriously expand
$\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|^{- 1}$ in
inverse powers of $\left| \overrightarrow{r} \right|$, work the
elementary integrals over $q^{'}$, and gather terms of the same order.
Fortunately this task has already been done for us in texts that treat
electrostatics using spherical harmonics. For
$\left| \overrightarrow{r} \right| > \left| \overrightarrow{r^{'}} \right|$,

$\int_{}^{}\frac{dq}{2p}\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|} = \sum_{n = 0}^{¥}{\frac{\left| \overrightarrow{r^{'}} \right|^{n}}{\begin{matrix}
\left| \overrightarrow{r} \right|^{n + 1}
\end{matrix}}\begin{matrix}
P_{n}\left( \text{cos}y \right)P_{n}\left( \text{cos \{}y \middle| \middle| ' \right)
\end{matrix}}$, (11.4)

where

$\text{cos}(y)º\frac{z}{\left| \overrightarrow{r} \right|},\begin{matrix}
 & 
\end{matrix}\text{cos}\left( y^{'} \right)º\frac{z^{'}}{\left| \overrightarrow{r^{'}} \right|}$
, (11.5)

and the $P_{n}$are the Legendre polynomials:

$\begin{array}{r}
P_{0}(x) = 1, \\
P_{1}(x) = x, \\
P_{2}(x) = \left( 3x^{2} - 1 \right)/2, \\
P_{3}(x) = \left( 5x^{3} - 3x \right)/2, \\
P_{4}(x) = \left( \text{35}x^{4} - \text{30}x^{2} + 3 \right)/8,\begin{matrix}

\end{matrix}\begin{matrix}
\text{etc}
\end{matrix}\text{.}
\end{array}$ (11.6)

Inserting expression (11.4) into eqn. (11.3) gives

$j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\frac{¶S^{'}}{¶z^{'}}}\left\{ \frac{1}{\left| \overrightarrow{r} \right|} + \frac{z^{'}z}{\left| \overrightarrow{r} \right|^{3}} + \frac{\left( 2z^{2} - r^{2} \right)\left( 2z^{} - r^{} \right)}{4\left| \overrightarrow{r} \right|^{5}} + \frac{\left( 2z^{3} - 3\text{zr}^{2} \right)\left( 2z^{} - 3z^{'}r^{} \right)}{4\left| \overrightarrow{r} \right|^{7}} + \text{.}\text{.}\text{.} \right\}\begin{matrix}
\text{.}
\end{matrix}$ (11.7)

The leading (monopole) term vanishes after integration by parts in
$z^{'}$ unless $S^{'}$ does not vanish as
$\left| z^{'} \right| \rightarrow ¥$:

$j\text{=-}\frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}S^{'}}\left\{ \frac{z}{\left| \overrightarrow{r} \right|^{3}} + \frac{\left( 2z^{2} - r^{2} \right)z^{'}}{\left| \overrightarrow{r} \right|^{5}} + \frac{\left( 2z^{3} - 3\text{zr}^{2} \right)\left( 6z^{} - 3r^{} \right)}{4\left| \overrightarrow{r} \right|^{7}} + \text{.}\text{.}\text{.} \right\}$
(11.8)

> =$j_{\text{dipole}} + j_{\text{quadrupole}} + j_{\text{sextupole}} + \text{.}\text{.}\text{.}$
> . (11.9)

Equation (11-8) may be written in the compact form

> $j\text{=-}\frac{m_{0}}{4p}\left\lbrack \frac{m_{1}P_{1}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{2}} + \frac{m_{2}P_{2}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{3}} + \frac{m_{3}P_{3}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{4}} + \text{.}\text{.}\text{.} \right\rbrack$
> , (11.10)

where the $m_{i}$are the multipole moments:

> $m_{1} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)$ , (11.11a)
>
> $m_{2} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)(2z)$ ,
> (11.11b)

$m_{3} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)\frac{3\left( 2z^{2} - r^{2} \right)}{2}$
, (11.11c)

etc. Here $m_{1}$ is the dipole moment with the conventional definition
given in section 2.

The quadrupole ($m_{2}$) and sextupole ($m_{3}$) moments have been
defined in a way that makes eqn. (11-10) simple in appearance. Their
general tensor forms are unnecessary with solenoidal symmetry.

For the far field expansion to be useful we must specify a coordinate
origin, which should be at the magnet system center in some weighted
sense. For a single lens the origin would usually be at the lens center.
More generally the quadrupole moment can usually be zeroed by the right
choice of origin. However, if $m_{1}$ vanishes, then $m_{2}$is
independent of the coordinate origin. In this case $m_{3}$ may be made
to vanish.

Using $\overrightarrow{B^{\text{residual}}} = \overrightarrow{Ñ}j$, the
dipole components of the far field are

$B_{r}^{\text{far}} = - \frac{\partial}{\partial r}\frac{\mu_{0}}{4\pi}\frac{m_{1}z}{\sqrt{{r^{2} + z^{2}}^{3}}} = \frac{\mu_{0}}{4\pi}\frac{3m_{1}\text{zr}}{\sqrt{{r^{2} + z^{2}}^{5}}}$,
(11.12a)

> $B_{z}^{\text{far}} = - \frac{\partial}{\partial z}\frac{\mu_{0}}{4\pi}\frac{m_{1}z}{\sqrt{{r^{2} + z^{2}}^{3}}}\text{=+}\frac{\mu_{0}}{4\pi}\frac{m_{1}(2z^{2} - r^{2})}{\sqrt{{r^{2} + z^{2}}^{5}}}$
> . (11.12b)
>
> **11.2 External field expansion**

To obtain an expansion for the external field we take $r^{'} < r$, but
make no assumption about $z$and $z^{'}$. Then
$r^{'}r < \left( z^{'} - z \right)^{2} + r^{2}$and we write equation
(11.3) as

> $j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\int_{}^{}\frac{dq^{'}}{2p}\frac{1}{\sqrt{1 - \left\lbrack \frac{{2r^{'}r\text{cos \{}q}^{'} - r^{}}{\left( z^{'} - z \right)^{2} + r^{2}} \middle| \right\rbrack}}$
> (11.13)
>
> $\begin{array}{r}
>  = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\int_{}^{}\frac{dq^{'}}{2p}\begin{matrix}
> \text{.}
> \end{matrix}\begin{matrix}
>  & & & & 
> \end{matrix}\left\{ 1 + \frac{1}{2}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack + \frac{3}{8}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{2} + \frac{5}{\text{16}}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{3} + \frac{\text{35}}{\text{128}}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{4} + \text{.}\text{.}\text{.} \right\},
> \end{array}$ (11.14)

where $\left\lbrack \text{.}\text{.}\text{.} \right\rbrack$ denotes the
bracketed expression in eqn. (11.13). After the integrals over $q^{'}$
are worked and terms of like order in $r^{'}$ are grouped we have

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}} \cdot \\
\left\{ \frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}} + r^{}\left\lbrack - \frac{1}{2}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}} + \frac{3}{4}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} \right\rbrack + r^{}\left\lbrack \frac{3}{8}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} - \frac{\text{15}}{8}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{7}}} + \frac{\text{105}}{\text{64}}\frac{r^{4}}{\sqrt{{\text{.}\text{.}\text{.}}^{9}}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}
\end{array}$ (11.15)

The radial moments defined in section 2 may now be used:

$S(z)\overline{r^{2}}(z) = 2\int_{}^{}{dr^{'}r^{'}S\left( r^{'},z \right)}$
, (11.16a)

$S(z)\overline{r^{4}}(z) = 4\int_{}^{}{dr^{'}r^{}S\left( r^{'},z \right)}$
, (11.16b)

$S(z)\overline{r^{6}}(z) = 6\int_{}^{}{dr^{'}r^{}S\left( r^{'},z \right)}$
; (11.16c)

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}\left\{ \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{2}}\left( z^{'} \right)}{2} \right\rbrack \cdot \frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}} + \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{4}}\left( z^{'} \right)}{4} \right\rbrack \cdot \right.\ } \\
\left. \ \left\lbrack - \frac{1}{2}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}} + \frac{3}{4}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} \right\rbrack + \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{6}}\left( z^{'} \right)}{6} \right\rbrack \cdot \left\lbrack \frac{3}{8}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} - \frac{\text{15}}{8}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{7}}} + \frac{\text{105}}{\text{64}}\frac{r^{4}}{\sqrt{{\text{.}\text{.}\text{.}}^{9}}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}
\end{array}$ (11.17)

This series appears to be complicated, but in fact we have again
generated Legendre polynomials. Let

$\text{cos}(g)º\frac{\left( z^{'} - z \right)}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}$
. (11.18)

Then eqn. (11.17) may be written

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\left\{ \frac{P_{0}\left( \text{cos}g \right)}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{2}\left( z^{'} \right)}}{2} \right\rbrack \right.\  \\
\left. \ \begin{matrix}
 & & & & 
\end{matrix} - \frac{1}{2}\frac{P_{2}\left( \text{cos}g \right)}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{4}\left( z^{'} \right)}}{4} \right\rbrack + \frac{3}{8}\frac{P_{4}\left( \text{cos}g \right)}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{6}\left( z^{'} \right)}}{6} \right\rbrack + \text{.}\text{.}\text{.} \right\}\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (11.19)

We treat the semi-infinite distribution

$S(r,z) = S(r)H( - z)$, (11.20)

with $S(r = 0)ºS\begin{matrix}

\end{matrix}$ from eqn. (11.17):

$j = \frac{- m_{0}S}{2}\left\{ \frac{\overline{r^{2}}}{2}\frac{1}{\left| \overrightarrow{r} \right|} + \frac{\overline{r^{4}}}{4}\left\lbrack \frac{- 1}{2\left| \overrightarrow{r} \right|^{3}} + \frac{3r^{2}}{4\left| \overrightarrow{r} \right|^{5}} \right\rbrack + \frac{\overline{r^{6}}}{6}\left\lbrack \frac{3}{8\left| \overrightarrow{r} \right|^{5}} - \frac{\text{15}r^{2}}{8\left| \overrightarrow{r} \right|^{7}} + \frac{\text{105}r^{4}}{\text{64}\left| \overrightarrow{r} \right|^{9}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}$
(11.21)

These are monopole, quadrupole, octopole, etc. terms. For a thin annular
layer at radius R,

$\overline{r^{2}} = R^{2},\begin{matrix}

\end{matrix}$ . (11.22)

For a uniform thick layer ($R_{1} \leq r \leq R_{2}$),

$\overline{r^{2}} = \frac{R_{2}^{3} - R_{1}^{3}}{3\left( R_{2} - R_{1} \right)},\begin{matrix}

\end{matrix}$ . (11.23)

The leading (monopole) term of eqn. (11.21) gives

$\overrightarrow{B^{\text{ex}}} = \overrightarrow{\nabla}\phi \approx \frac{\mu_{0}S\overline{r^{2}}}{4}\frac{\overrightarrow{r}}{\left| \overrightarrow{r} \right|^{3}} =$
, (11.24)

as expected.

**11. [Far Field and External Field]{.underline}**

Because a solenoid system may have a large dipole moment, its field can
interfere with nearby experiments and the operation of accelerator
components. A short magnetic bend, which resembles a solenoid turned on
its side, also has a net dipole moment, but a permeable yoke is usually
used to confine its return flux. This technique is not as effective for
transport solenoids because the yoke must not block the vacuum bore.
Also, if annular yokes are used in an induction linac they may
considerably increase the diameter of adjacent induction cores and
thereby increase the overall cost. It is possible to greatly reduce the
far field by alternating the polarity of solenoid lenses, but this
reduces their effectiveness for beam transport and actually increases
the external field (just outside the wire). We are therefore motivated
to evaluate the far and external fields in the absence of return flux
confinement.

The far field and external field can be approximated using two different
series expansions of the global potential $j$. Both expansions may have
contributions from several or many individual magnets, so it may be
misleading to consider only one magnet at a time. Recall from section 5,

> $\overrightarrow{B} = m_{0}S(r,z){\hat{e}}_{z} + \overrightarrow{Ñ}j$,
> (11.1)

where S(r,z) is the cumulative surface current density:

$S(r,z) = \int_{}^{}{dr^{'}}J_{\theta}\left( r^{'},z \right)$ . (11.2)

The Green function solution for $j$ is, from section (6):

$j = \frac{m_{0}}{4p}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\int_{}^{}{dq^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{} + r^{2} - 2r^{'}r\text{cos}\left( q^{'} - q \right)}}\begin{matrix}
\text{.}
\end{matrix}$ (11.3)

**11.1 Far field expansion**

Consider first the far field expansion in inverse powers of
$\left| \overrightarrow{r} \right| = \left( r^{2} + z^{2} \right)^{12}$,
with r and z both considered to be large of order
$\left| \overrightarrow{r} \right|$, which is larger than
$\left| \overrightarrow{r^{'}} \right|$. We could laboriously expand
$\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|^{- 1}$ in
inverse powers of $\left| \overrightarrow{r} \right|$, work the
elementary integrals over $q^{'}$, and gather terms of the same order.
Fortunately this task has already been done for us in texts that treat
electrostatics using spherical harmonics. For
$\left| \overrightarrow{r} \right| > \left| \overrightarrow{r^{'}} \right|$,

$\int_{}^{}\frac{dq}{2p}\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|} = \sum_{n = 0}^{¥}{\frac{\left| \overrightarrow{r^{'}} \right|^{n}}{\begin{matrix}
\left| \overrightarrow{r} \right|^{n + 1}
\end{matrix}}\begin{matrix}
P_{n}\left( \text{cos}y \right)P_{n}\left( \text{cos \{}y \middle| \middle| ' \right)
\end{matrix}}$, (11.4)

where

$\text{cos}(y)º\frac{z}{\left| \overrightarrow{r} \right|},\begin{matrix}
 & 
\end{matrix}\text{cos}\left( y^{'} \right)º\frac{z^{'}}{\left| \overrightarrow{r^{'}} \right|}$
, (11.5)

and the $P_{n}$are the Legendre polynomials:

$\begin{array}{r}
P_{0}(x) = 1, \\
P_{1}(x) = x, \\
P_{2}(x) = \left( 3x^{2} - 1 \right)/2, \\
P_{3}(x) = \left( 5x^{3} - 3x \right)/2, \\
P_{4}(x) = \left( \text{35}x^{4} - \text{30}x^{2} + 3 \right)/8,\begin{matrix}

\end{matrix}\begin{matrix}
\text{etc}
\end{matrix}\text{.}
\end{array}$ (11.6)

Inserting expression (11.4) into eqn. (11.3) gives

$j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}\frac{¶S^{'}}{¶z^{'}}}\left\{ \frac{1}{\left| \overrightarrow{r} \right|} + \frac{z^{'}z}{\left| \overrightarrow{r} \right|^{3}} + \frac{\left( 2z^{2} - r^{2} \right)\left( 2z^{} - r^{} \right)}{4\left| \overrightarrow{r} \right|^{5}} + \frac{\left( 2z^{3} - 3\text{zr}^{2} \right)\left( 2z^{} - 3z^{'}r^{} \right)}{4\left| \overrightarrow{r} \right|^{7}} + \text{.}\text{.}\text{.} \right\}\begin{matrix}
\text{.}
\end{matrix}$ (11.7)

The leading (monopole) term vanishes after integration by parts in
$z^{'}$ unless $S^{'}$ does not vanish as
$\left| z^{'} \right| \rightarrow ¥$:

$j\text{=-}\frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}S^{'}}\left\{ \frac{z}{\left| \overrightarrow{r} \right|^{3}} + \frac{\left( 2z^{2} - r^{2} \right)z^{'}}{\left| \overrightarrow{r} \right|^{5}} + \frac{\left( 2z^{3} - 3\text{zr}^{2} \right)\left( 6z^{} - 3r^{} \right)}{4\left| \overrightarrow{r} \right|^{7}} + \text{.}\text{.}\text{.} \right\}$
(11.8)

> =$j_{\text{dipole}} + j_{\text{quadrupole}} + j_{\text{sextupole}} + \text{.}\text{.}\text{.}$
> . (11.9)

Equation (11-8) may be written in the compact form

> $j\text{=-}\frac{m_{0}}{4p}\left\lbrack \frac{m_{1}P_{1}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{2}} + \frac{m_{2}P_{2}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{3}} + \frac{m_{3}P_{3}\left( \text{cos}y \right)}{\left| \overrightarrow{r} \right|^{4}} + \text{.}\text{.}\text{.} \right\rbrack$
> , (11.10)

where the $m_{i}$are the multipole moments:

> $m_{1} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)$ , (11.11a)
>
> $m_{2} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)(2z)$ ,
> (11.11b)

$m_{3} = 2\pi\int_{}^{}\text{dz}\int_{}^{}\text{drr}S(r,z)\frac{3\left( 2z^{2} - r^{2} \right)}{2}$
, (11.11c)

etc. Here $m_{1}$ is the dipole moment with the conventional definition
given in section 2.

The quadrupole ($m_{2}$) and sextupole ($m_{3}$) moments have been
defined in a way that makes eqn. (11-10) simple in appearance. Their
general tensor forms are unnecessary with solenoidal symmetry.

For the far field expansion to be useful we must specify a coordinate
origin, which should be at the magnet system center in some weighted
sense. For a single lens the origin would usually be at the lens center.
More generally the quadrupole moment can usually be zeroed by the right
choice of origin. However, if $m_{1}$ vanishes, then $m_{2}$is
independent of the coordinate origin. In this case $m_{3}$ may be made
to vanish.

Using $\overrightarrow{B^{\text{residual}}} = \overrightarrow{Ñ}j$, the
dipole components of the far field are

$B_{r}^{\text{far}} = - \frac{\partial}{\partial r}\frac{\mu_{0}}{4\pi}\frac{m_{1}z}{\sqrt{{r^{2} + z^{2}}^{3}}} = \frac{\mu_{0}}{4\pi}\frac{3m_{1}\text{zr}}{\sqrt{{r^{2} + z^{2}}^{5}}}$,
(11.12a)

> $B_{z}^{\text{far}} = - \frac{\partial}{\partial z}\frac{\mu_{0}}{4\pi}\frac{m_{1}z}{\sqrt{{r^{2} + z^{2}}^{3}}}\text{=+}\frac{\mu_{0}}{4\pi}\frac{m_{1}(2z^{2} - r^{2})}{\sqrt{{r^{2} + z^{2}}^{5}}}$
> . (11.12b)
>
> **11.2 External field expansion**

To obtain an expansion for the external field we take $r^{'} < r$, but
make no assumption about $z$and $z^{'}$. Then
$r^{'}r < \left( z^{'} - z \right)^{2} + r^{2}$and we write equation
(11.3) as

> $j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\int_{}^{}\frac{dq^{'}}{2p}\frac{1}{\sqrt{1 - \left\lbrack \frac{{2r^{'}r\text{cos \{}q}^{'} - r^{}}{\left( z^{'} - z \right)^{2} + r^{2}} \middle| \right\rbrack}}$
> (11.13)
>
> $\begin{array}{r}
>  = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}}\frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\int_{}^{}\frac{dq^{'}}{2p}\begin{matrix}
> \text{.}
> \end{matrix}\begin{matrix}
>  & & & & 
> \end{matrix}\left\{ 1 + \frac{1}{2}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack + \frac{3}{8}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{2} + \frac{5}{\text{16}}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{3} + \frac{\text{35}}{\text{128}}\left\lbrack \text{.}\text{.}\text{.} \right\rbrack^{4} + \text{.}\text{.}\text{.} \right\},
> \end{array}$ (11.14)

where $\left\lbrack \text{.}\text{.}\text{.} \right\rbrack$ denotes the
bracketed expression in eqn. (11.13). After the integrals over $q^{'}$
are worked and terms of like order in $r^{'}$ are grouped we have

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\int_{}^{}{dr^{'}r^{'}}\frac{¶S^{'}}{¶z^{'}} \cdot \\
\left\{ \frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}} + r^{}\left\lbrack - \frac{1}{2}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}} + \frac{3}{4}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} \right\rbrack + r^{}\left\lbrack \frac{3}{8}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} - \frac{\text{15}}{8}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{7}}} + \frac{\text{105}}{\text{64}}\frac{r^{4}}{\sqrt{{\text{.}\text{.}\text{.}}^{9}}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}
\end{array}$ (11.15)

The radial moments defined in section 2 may now be used:

$S(z)\overline{r^{2}}(z) = 2\int_{}^{}{dr^{'}r^{'}S\left( r^{'},z \right)}$
, (11.16a)

$S(z)\overline{r^{4}}(z) = 4\int_{}^{}{dr^{'}r^{}S\left( r^{'},z \right)}$
, (11.16b)

$S(z)\overline{r^{6}}(z) = 6\int_{}^{}{dr^{'}r^{}S\left( r^{'},z \right)}$
; (11.16c)

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}\left\{ \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{2}}\left( z^{'} \right)}{2} \right\rbrack \cdot \frac{1}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}} + \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{4}}\left( z^{'} \right)}{4} \right\rbrack \cdot \right.\ } \\
\left. \ \left\lbrack - \frac{1}{2}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}} + \frac{3}{4}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} \right\rbrack + \left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{6}}\left( z^{'} \right)}{6} \right\rbrack \cdot \left\lbrack \frac{3}{8}\frac{1}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}} - \frac{\text{15}}{8}\frac{r^{2}}{\sqrt{{\text{.}\text{.}\text{.}}^{7}}} + \frac{\text{105}}{\text{64}}\frac{r^{4}}{\sqrt{{\text{.}\text{.}\text{.}}^{9}}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}
\end{array}$ (11.17)

This series appears to be complicated, but in fact we have again
generated Legendre polynomials. Let

$\text{cos}(g)º\frac{\left( z^{'} - z \right)}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}$
. (11.18)

Then eqn. (11.17) may be written

$\begin{array}{r}
j = \frac{m_{0}}{2}\int_{}^{}{dz^{'}}\left\{ \frac{P_{0}\left( \text{cos}g \right)}{\sqrt{\left( z^{'} - z \right)^{2} + r^{2}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{2}\left( z^{'} \right)}}{2} \right\rbrack \right.\  \\
\left. \ \begin{matrix}
 & & & & 
\end{matrix} - \frac{1}{2}\frac{P_{2}\left( \text{cos}g \right)}{\sqrt{{\text{.}\text{.}\text{.}}^{3}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{4}\left( z^{'} \right)}}{4} \right\rbrack + \frac{3}{8}\frac{P_{4}\left( \text{cos}g \right)}{\sqrt{{\text{.}\text{.}\text{.}}^{5}}}\left\lbrack \frac{¶}{¶z^{'}}\frac{S\left( z^{'} \right)\overline{r^{6}\left( z^{'} \right)}}{6} \right\rbrack + \text{.}\text{.}\text{.} \right\}\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (11.19)

We treat the semi-infinite distribution

$S(r,z) = S(r)H( - z)$, (11.20)

with $S(r = 0)ºS\begin{matrix}

\end{matrix}$ from eqn. (11.17):

$j = \frac{- m_{0}S}{2}\left\{ \frac{\overline{r^{2}}}{2}\frac{1}{\left| \overrightarrow{r} \right|} + \frac{\overline{r^{4}}}{4}\left\lbrack \frac{- 1}{2\left| \overrightarrow{r} \right|^{3}} + \frac{3r^{2}}{4\left| \overrightarrow{r} \right|^{5}} \right\rbrack + \frac{\overline{r^{6}}}{6}\left\lbrack \frac{3}{8\left| \overrightarrow{r} \right|^{5}} - \frac{\text{15}r^{2}}{8\left| \overrightarrow{r} \right|^{7}} + \frac{\text{105}r^{4}}{\text{64}\left| \overrightarrow{r} \right|^{9}} \right\rbrack + \text{.}\text{.}\text{.} \right\}\text{.}$
(11.21)

These are monopole, quadrupole, octopole, etc. terms. For a thin annular
layer at radius R,

$\overline{r^{2}} = R^{2},\begin{matrix}

\end{matrix}$ . (11.22)

For a uniform thick layer ($R_{1} \leq r \leq R_{2}$),

$\overline{r^{2}} = \frac{R_{2}^{3} - R_{1}^{3}}{3\left( R_{2} - R_{1} \right)},\begin{matrix}

\end{matrix}$ . (11.23)

The leading (monopole) term of eqn. (11.21) gives

$\overrightarrow{B^{\text{ex}}} = \overrightarrow{\nabla}\phi \approx \frac{\mu_{0}S\overline{r^{2}}}{4}\frac{\overrightarrow{r}}{\left| \overrightarrow{r} \right|^{3}} =$
, (11.24)

as expected.

**12. [Field in the Wire Layer]{.underline}**

We now have several ways to calculate the field within the wire layer.
Near the layer ends an accurate calculation requires the general
formalism of section 6 or the thick layer solution (equation (9.3)),
plus the contribution from nearby magnets. Away from the wire ends the
external field expansion of section 11 plus the direct field
$m_{0}S(r,z){\hat{e}}_{z}$ is applicable.

The peak field in the wire layer is of great interest since it limits
the averaged critical current density (sections 2 and 20). [Usually this
peak field occurs at the inner edge of the wire layer at the middle of a
magnet]{.underline}. The logarithmic singularity of $B_{r}$ at the ends
of a thin layer is not representative of realistic cases. The peak field
can be accurately estimated using the first two terms of the near-axis
expansion if the magnet is not too short
$\left( R/l \leq \text{.}5 \right)$:

$B_{z}»B_{0}(z) - B_{0^{¢¢}}(z)\frac{r^{2}}{4}$ , (12.1a)

$B_{r}\text{»-}B_{0^{¢}}(z)\frac{r}{2} + B_{0^{¢¢¢}}(z)\frac{r^{3}}{\text{16}}$
. (12.1b)

We use the thick layer formula for a single magnet, eqn. (9.5), to
obtain $B_{0}$ and $B_{0^{¢¢}}$ at the magnet center (z=0), and by
symmetry $B_{r}$ vanishes there. At the inner wire radius $R_{1}$ this
gives for a magnet of length $l$:

$B_{\text{max}}@B_{0}(0) - \frac{B_{0^{¢¢}}(0)R_{1}^{2}}{4}\begin{matrix}

\end{matrix}$

$\left\{ \text{Log}\left\lbrack \frac{2R_{2}/\ell + \sqrt{1 + \left( 2R_{2}/\ell \right)^{2}}}{2R_{1}/\ell + \sqrt{1 + \left( 2R_{1}/\ell \right)^{2}}} \right\rbrack + \frac{1}{4}\left( \frac{2R_{1}}{\ell} \right)^{2}\left\lbrack \frac{\left( 2R_{2}/\ell \right)^{3}}{\sqrt{{1 + \left( 2R_{2}/\ell \right)^{2}}^{3}}} - \frac{\left( 2R_{1}/\ell \right)^{3}}{\sqrt{{1 + \left( 2R_{1}/\ell \right)^{2}}^{3}}} \right\rbrack \right\}$.
(12.2)

This formula may be written in dimensionless form using the definitions

$\overline{R} = \frac{R_{1} + R_{2}}{2},\begin{matrix}

\end{matrix}$ ; (12.3)

$\begin{array}{r}
B_{\text{max}}^{\text{approx}}@\left( m_{0}S \right)\frac{1}{\alpha D}\left\{ \text{Log}\left\lbrack \frac{\alpha\left( 1 + \frac{D}{2} \right) + \sqrt{1 + \alpha^{2}\left( 1 + \frac{D}{2} \right)^{2}}}{\alpha\left( 1 - \frac{D}{2} \right) + \sqrt{1 + \alpha^{2}\left( 1 - \frac{D}{2} \right)^{2}}} \right\rbrack \right.\  \\
\left. \  + \frac{1}{4}\alpha^{5}\left( 1 - \frac{D}{2} \right)^{2}\left\lbrack \frac{\left( 1 + \frac{D}{2} \right)^{3}}{\sqrt{{1 + \alpha^{2}\left( 1 + \frac{D}{2} \right)^{2}}^{3}}} - \frac{\left( 1 - \frac{D}{2} \right)^{3}}{\sqrt{{1 + \alpha^{2}\left( 1 - \frac{D}{2} \right)^{2}}^{3}}} \right\rbrack \right\}
\end{array}$

(12.4) (

> .

The degree of accuracy of eqn. (12.4) is indicated in the following
table, where the ratio $B_{\text{max}}/B_{\text{max}}^{\text{approx}}$
is given for several values of $D$ and $\alpha$.

  ----------------- ----------------- ----------------- -----------------
  $$D$$             1.0               .5                .25

  1                 .998320           .999746           .999990

  .5                .992932           .998618           .999948

  .25               .987348           .997435           .999904

  .10               .982344           .996471           .999868

  .04               .979903           .996029           .999851
  ----------------- ----------------- ----------------- -----------------

The exact values of $B_{\text{max}}/m_{0}S$ for several values of
$\alpha$ and $D$ are computed and tabulated below from

$\begin{array}{r}
B_{\text{max}} = \frac{m_{0}S}{R_{2} - R_{1}}\int_{}^{}{\text{dR}\int_{}^{}{\frac{dq}{2p}\left\lbrack \frac{R^{2} - \text{RR}_{1}\text{cos}(q)}{R^{2} + R_{1}^{2} - 2\text{RR}_{1}\text{cos}(q)} \right\rbrack}} \cdot \\
\begin{matrix}
 & & & & 
\end{matrix} \cdot \frac{l/2}{\sqrt{l^{2}/4 + R^{2} + R_{1}^{2} - 2\text{RR}_{1}\text{cos}(q)}}\begin{matrix}

\end{matrix}\text{.}
\end{array}$ (12.5)

  ----------------- ----------------- ----------------- -----------------
  $$D$$             1.0               .5                .25

  1.0               .742700           .897437           .968828

  .5                .776529           .907269           .971029

  .25               .798374           .912353           .971858

  .10               .812152           .915356           .972267

  .04               .817666           .916536           .972412
  ----------------- ----------------- ----------------- -----------------

**13. [Current to Produce a Given Field]{.underline}**

When designing a magnet system one usually asks for the current
distribution to produce a desired field. For solenoids it is clear that
any on-axis field profile with $B \leq \text{10}T$ can be made with
scale length variations down to a few centimeters. But for beam
transport and focusing this may also require an unacceptable small inner
radius for the wire layer. Increasing the wire layer radius increases
the scale length of field variations.

In the present calculation we assume the current is concentrated in a
thin layer of fixed radius $R$; $J_{q} = S(z)d(r - R)$. For given
$B_{0}(z)$ it is always possible to solve for $S(z)$, but the field at
the wire layer may be unacceptably large if the variations of $B_{0}(z)$
have scale length smaller than $R$. This problem may be evaded in some
applications by using small radius solenoids where the beam radius is
small. Recall the thin layer formula for the on-axis field (Section 3),

$B_{0}(z) = \frac{m_{0}R^{2}}{2}\int_{}^{}{dz^{'}}\frac{S\left( z^{'} \right)}{\sqrt{{\left( z^{'} - z \right)^{2} + R^{2}}^{3}}}$
. (13.1)

Defining the Fourier integral transformations

$\left\lbrack {\overset{\sim}{B}}_{0}(k),\overset{\sim}{S}(k) \right\rbrack = \int_{\text{-¥}}^{\text{+¥}}\text{dz}e^{- \text{ikz}}\left\lbrack B_{0}(z),S(z) \right\rbrack$
, (13.2)

we get from eqn. (13.1)

$\begin{array}{r}
{\overset{\sim}{B}}_{0}(k) = \frac{m_{0}R^{2}}{2}\int_{\text{-¥}}^{\text{+¥}}{dz^{'}S\left( z^{'} \right)}\int_{\text{-¥}}^{\text{+¥}}\text{dz}\frac{e^{- \text{ikz}}}{\sqrt{{\left( z^{'} - z \right)^{2} + R^{2}}^{3}}} \\
\frac{m_{0}}{2}\int_{}^{}{dz^{'}}e^{{- \text{ik \{}z}^{'}}
\end{array}S\left( z^{'} \right)\int_{}^{}\text{dx}\frac{e^{- \text{ikRx}}}{\sqrt{{x^{2} + 1}^{3}}},$
(13.3)

with $x = \left( z - z^{'} \right)/R$. Inserting the identity

$\int_{}^{}\text{dx}\frac{\text{cos}\left( \text{kRx} \right)}{\sqrt{{x^{2} + 1}^{3}}} = \left| \text{kR} \right|K_{1}\left( \left| \text{kR} \right| \right)$
(13.4)

gives the general relation

${\overset{\sim}{B}}_{0}(k) = m_{0}\overset{\sim}{S}(k)\left| \text{kR} \right|K_{1}\left( \left| \text{kR} \right| \right)$
, (13.5)

where $K_{1}$ is the modified Bessel function. Equation (13.5) gives
$\overset{\sim}{S}$ if ${\overset{\sim}{B}}_{0}$ is derived from a given
field.

For the specific example,

$B_{0}(z) = B_{\text{00}}e^{- z^{2}/L^{2}}$ , (13.6)

the transform of $B_{0}(z)$ and inversion for $S(z)$ are

${\overset{\sim}{B}}_{0}(k) = B_{\text{00}}\sqrt{p}\text{Le}^{- k^{2}L^{2}/4}$
, (13.7)

$m_{0}S(z) = \int_{\text{-¥}}^{\text{+¥}}\frac{\text{dk}}{2p}\begin{matrix}

\end{matrix}$ . (13.8)

Setting $y = \text{kL}$ gives the convenient form for computations:

$\frac{\mu_{0}S(z)}{B_{\text{00}}} = \frac{1}{\sqrt{\pi}}\frac{L}{R}\int_{}^{}\text{dy}\frac{e^{- y24}\text{cos}(yz/L)}{\text{yK}_{1}(yR/L)}$
. (13.9)

For $R/L \rightarrow 0$ and $\frac{|z|}{L}$ not too large, this formula
goes over to the expected form (using $\text{xK}_{1}(x) \rightarrow 1$),

$m_{0}S(z) \rightarrow \begin{matrix}
B_{\text{00}}e^{- z^{2}/L^{2}}
\end{matrix}$. (13.10)

However for finite $R/L$, $m_{0}S(z)$ exceeds $B_{\text{00}}$ at
$z = 0$, e.g. by a factor of 2.5 when $R/L$=1.0. Also, oscillations of
$S(z)$ appear; these provide guidance for the layout of wire that can
produce a desired on-axis field. The following figures plot
$m_{0}S(z)/B_{\text{00}}$ for various values of $R$, with
$L = 1\text{.}0$.

R/L=.01

![](media/image10.png){width="3.0in" height="1.84375in"}

R/L=.5

![](media/image11.png){width="3.0in" height="1.84375in"}

R/L=1

![](media/image12.png){width="3.0in" height="1.84375in"}

R/L=2

![](media/image13.png){width="3.0in" height="1.84375in"}

**14. [Misaligned and Tilted Solenoid]{.underline}**

So far we have only treated axisymmetric systems, with field components
$B_{r}(r,z)$ and$B_{z}(r,z)$. There are always small deviations from
this convenient symmetry because of imperfections in design and
manufacture, stress, and positioning errors. These imperfections can
have a significant effect on a transported ion beam, so it is necessary
to specify maximum allowed deviations and make precise field
measurements. However, the restricted problem of the effect of
positioning errors on the field can be treated analytically as presented
below.

The field of a single magnet is considered to be ideal and perfectly
characterized in its own coordinate frame. It greatly simplifies the
formalism to use Cartesian coordinates. For example if $B_{0}(z)$ is the
known on-axis field in the absence of a positioning error, then

$B_{z}^{I}(x,y,z) = B_{0}(z) - B_{0}^{''}(z)\frac{x^{2} + y^{2}}{4} + \text{.}\text{.}\text{.}$
, (14.1a)

$B_{x}^{I}(x,y,z) = - \frac{B_{0}^{'}(z)x}{2} + \frac{B_{0}^{'''}(z)x\left( x^{2} + y^{2} \right)}{\text{16}} - \text{.}\text{.}\text{.}$
, (14.1b)

$B_{y}^{I}(x,y,z) = - \frac{B_{0}^{'}(z)y}{2} + \frac{B_{0}^{'''}(z)y\left( x^{2} + y^{2} \right)}{\text{16}} - \text{.}\text{.}\text{.}$
(14.1c)

is the magnet's ideal near-axis field.

A misalignment is simply a translation of the magnet center by the
vector

$\overrightarrow{T} = \Delta_{x}{\hat{e}}_{x} + \Delta_{y}{\hat{e}}_{y} + \Delta_{z}{\hat{e}}_{z}$
, (14.2)

with displacements $D_{i}$ typically less than $\pm \text{10}^{- 3}m$.
The misaligned field is obtained from eqn. (14-1) simply by substituting
the displacement:

$\overrightarrow{B} = \overrightarrow{B^{I}}\left( x - D_{x},y - D_{y},z - D_{z} \right)$
. (14.3)

A second type of positioning error is a tilt, which holds a magnet's
center fixed and orients its axis at an angle $\alpha$ with respect to
the system axis, typically with $\alpha$ less than $\text{10}^{- 3}$
radians. The tilted magnet axis points in the direction of the unit
vector

$\hat{e} = \text{sin}(\alpha)\text{cos}(\beta){\hat{e}}_{x} + \text{sin}(\alpha)\text{sin}(\beta){\hat{e}}_{y} + \text{cos}(\alpha){\hat{e}}_{z}$
, (14.4)

where $\beta$ is the angle of the projection of $\hat{e}$ in the $x - y$
plane. The two angles $(\alpha,\beta)$ are the polar coordinates of
$\hat{e}$, the usual symbols $(q,j)$ having been used for other
quantities. A net azimuthal rotation of a magnet around its own axis
does not need to be considered since it does not change the field.

The coordinates of a fixed point in the system frame are denoted by
$(x,y,z)$ and in the translated and rotated magnet frame by $(X,Y,Z)$.
These coordinate sets are related by the inhomogeneous linear
transformation

$\begin{array}{r}
\left( X \middle| \right)\left( Y \middle| \right) \\

\end{array}$ , (14.5)

where $\left( R_{\text{ij}} \right)$is the orthogonal matrix

$\left( R_{\text{ij}} \right) = \begin{array}{r}
\left( \text{cos}(\alpha)\text{cos}^{2}(\beta) + \text{sin}^{2}(\beta),\left( \text{cos}(\alpha) - 1 \right)\text{sin}(\beta)\text{cos}(\beta), - \text{sin}(\alpha)\text{cos}(\beta) \middle| \right)\left( \left( \text{cos}(\alpha) - 1 \right)\text{sin}(\beta)\text{cos}(\beta)\text{,cos}(\alpha)\text{sin}^{2}(\beta) + \text{cos}^{2}(\beta), - \text{sin}(\alpha)\text{sin}(\beta) \middle| \right) \\

\end{array}$ . (14.6)

Eqn. (14.6) may be derived by multiplying the matrices for successive
rotations of the magnet's coordinate frame through angles
$\beta,\begin{matrix}

\end{matrix}$.

The field components of the rotated and translated magnet are then found
by applying the transpose of $\left( R_{\text{ij}} \right)$:

> $B_{x}(x,y,z) = R_{\text{11}}B_{x}^{I}(X,Y,Z) + R_{\text{21}}B_{y}^{I}(X,Y,Z) + R_{\text{31}}B_{z}^{I}(X,Y,Z)$
> , (14.7a)
>
> $B_{y}(x,y,z) = R_{\text{12}}B_{x}^{I}(X,Y,Z) + R_{\text{22}}B_{y}^{I}(X,Y,Z) + R_{\text{32}}B_{z}^{I}(X,Y,Z)$
> , (14.7b)
>
> $B_{z}(x,y,z) = R_{\text{13}}B_{x}^{I}(X,Y,Z) + R_{\text{23}}B_{y}^{I}(X,Y,Z) + R_{\text{33}}B_{z}^{I}(X,Y,Z)$
> . (14.7c)

The rotation angles $(\alpha,\beta)$ are general, but since $\alpha$ is
usually very small for tilt errors, the rotation matrix can be well
approximated by letting $\text{cos}(\alpha) \rightarrow 1$ and
$\text{sin}(\alpha) \rightarrow \alpha$. Defining the small tilt angles
in the $x$ and $y$ directions:

$q_{x} = \alpha\begin{matrix}
\text{cos}(\beta),\begin{matrix}
 & & & & 
\end{matrix}q_{y} = \alpha\begin{matrix}

\end{matrix}\text{sin}(\beta),
\end{matrix}$ (14.8a,b)

the approximate matrix is then

$\left( R_{\text{ij}} \right)»\begin{array}{r}
\left( 1\begin{matrix}
 & 
\end{matrix}0\begin{matrix}
 & 
\end{matrix} - q_{x} \middle| \right)\left( 0\begin{matrix}
 & 
\end{matrix}1\begin{matrix}
 & 
\end{matrix} - q_{y} \middle| \right) \\

\end{array}$ . (14.9)

This is not quite orthogonal and leads to small violations of
$\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$. Keeping only terms
which are linear in $$, we have the approximate near-axis field:

$B_{x}\text{»-}\frac{1}{2}B_{0}^{'}\left( z - D_{z} \right)\left\lbrack \left( x - D_{x} \right) - q_{x}\left( z - D_{z} \right) \right\rbrack + q_{x}B_{0}\left( z - D_{z} \right)$
, (14.10a)

$B_{y}\text{»-}\frac{1}{2}B_{0}^{'}\left( z - D_{z} \right)\left\lbrack \left( y - D_{x} \right) - q_{y}\left( z - D_{z} \right) \right\rbrack + q_{y}B_{0}\left( z - D_{z} \right)$
, (14.10b)

$B_{z}»B_{0}\left( z - D_{z} \right)$ . (14.10c)

**15. [Multiple Channels -- Cross Talk]{.underline}**

We have treated a system of solenoidal magnets that are centered on a
single straight axis, except for misalignments and tilts, which are
usually small. More complicated systems are of interest in which
solenoids are positioned along several axes, which may have differing
orientations. For example a multiple-beam linac may employ parallel
closely-packed solenoid channels, while a system of solenoids designed
to bring multiple beams to a single focal spot may have several
orientations. In such cases the problem of interference between channels
(cross talk) is severe, and it is important to make good estimates of
the unwanted field components. This does not mean that any significant
cross talk is unacceptable. Generally the effects of unwanted fields can
be greatly reduced by a symmetrical layout of channels and the addition
of weak bends and quadrupoles around each channel.

We need a simple formula for the residual field of a set of solenoids
with arbitrary positions and orientations. This is provided by a
generalization of equation (11.21) for the global potential of a
semi-infinite current layer:

$J_{\theta}(r,z) = J_{\theta}(r)H( - z)$, (15.1)

$S = \int_{}^{}\text{dr}J_{\theta}(r)$, (15.2)

$\phi \simeq - \frac{\mu_{0}S}{2}\left\{ \frac{\overline{r^{2}}}{2}\frac{1}{\left| \overrightarrow{r} \right|} + \frac{\overline{r^{4}}}{4}\left\lbrack \frac{- 1}{2\left| \overrightarrow{r} \right|^{3}} + \frac{3r^{2}}{4\left| \overrightarrow{r} \right|^{5}} \right\rbrack \right\}$.
(15.3)

Here we truncated after the lowest moments:

$$$$

$S\overline{r^{2}} = \int_{}^{}\text{drr}^{2}J_{q},\begin{matrix}
S\overline{r^{4}} = \int_{}^{}\text{drr}^{4}J_{q}(r)
\end{matrix}$. (15.4)

For a thin current layer with radius R, we have
$\overline{r^{2}} = R^{2}$ and $\overline{r^{4}} = R^{4}$. For a thick
uniform current layer between $R_{1}$ and $R_{2}$:

$\overline{r^{2}} = \frac{R_{2}^{3} - R_{1}^{3}}{3\left( R_{2} - R_{1} \right)},\begin{matrix}
\overline{r^{4}} = \frac{R_{2}^{5} - R_{1}^{5}}{5\left( R_{2} - R_{1} \right)}
\end{matrix}$. (15.5a,b)

Suppose a solenoid end, denoted by subscript i, is located at
$\overrightarrow{r_{i}}$ and has direction ${\hat{e}}_{i}$; then its
global potential is obtained by substituting into eqn. (15.3)
$\left| \overrightarrow{r} \right| \rightarrow \left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|$
and

> $r^{2} = \left| \overrightarrow{r} \right|^{2} - z^{2} = \left| \overrightarrow{r} \right|^{2} - \left( \overrightarrow{r} \cdot {\hat{e}}_{z} \right)^{2} \rightarrow \left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|^{2} - \left( \overrightarrow{r} \cdot {\hat{e}}_{i} - \overrightarrow{r_{i}} \cdot {\hat{e}}_{i} \right)^{2}$
> . (15.6)

The quantities $S,\overline{r^{2}},$ and $\overline{r^{4}}$are
unchanged; since $\phi$ is a scalar quantity we have

$j_{i}\text{@-}\frac{m_{0}S}{2}\left\{ \frac{\overline{r^{2}}}{2}\frac{1}{\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|} + \frac{\overline{r^{4}}}{4}\left\lbrack \frac{- 1}{2\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|^{3}} + \frac{3}{4}\frac{\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|^{2} - \left( \overrightarrow{r} \cdot {\hat{e}}_{i} - \overrightarrow{r_{i}} \cdot {\hat{e}}_{i} \right)^{2}}{\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|^{5}} \right\rbrack \right\}$.
(15.7)

Each magnet contributes two terms of this type. For example if the
downstream end is at $\overrightarrow{r_{2}}$ and the magnet wire layer
length is $\ell$, then the upstream end is at
$\overrightarrow{r_{1}} = \overrightarrow{r_{2}} - l\hat{e}$. The
overall minus sign in eqn. (15.7) is changed to a plus for the upstream
magnet end. The total residual field is simply

$\overrightarrow{B^{\text{residual}}} = \overrightarrow{\nabla}\sum_{}^{}\underset{i}{}\phi_{i}$
, (15.8)

with applicable values of
$\pm S,\overrightarrow{r_{i}},\overline{r^{2}},$ and $\overline{r^{4}}$.

We have included only the two lowest-order terms in the expansion
of$\phi$; this should be sufficient for nearly any estimate. In fact the
first (monopole) term,

$j_{i}\text{»-}\frac{m_{0}S_{i}\left( \overline{r^{2}} \right)_{i}}{4\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|}$
, (15.9)

should be a good rough guide. The field of the "monopole" term,

> $\overrightarrow{{B_{i}}^{\text{residual}}} = \overrightarrow{Ñ}j_{i}»\frac{m_{0}S_{i}\left( \overline{r^{2}} \right)_{i}\left( \overrightarrow{r} - \overrightarrow{r_{i}} \right)}{4\left| \overrightarrow{r} - \overrightarrow{r_{i}} \right|^{3}}$
> , (15.10)

is just a spherically symmetrical field with the same net flux as the
magnet's direct field.

**16.** **[Magnetic Materials Basics]{.underline}**

There are several reasons for placing highly permeable material around
solenoidal wire layers: confining return flux, shaping fields, and
reducing stored field energy. It is also possible to make a permanent
magnet solenoid. Another application of magnetic materials is in the
core of an induction module for acceleration, and also in limiting the
high frequency impedance of an acceleration gap. The latter two examples
involve time-dependent fields that do not have solenoidal symmetry and
are not considered here.

**16.1 Ferromagnetism**

We are primarily interested in highly-permeable ferromagnetic materials,
i.e. those where parallel alignment of atomic magnetic moments produces
strong fields in response to weak externally-applied fields.
Paramagnetic and diamagnetic materials are not of interest since their
induced fields are very small. While the subject of magnetic materials
is too vast (and beyond our competence) to even begin a summary here, a
few remarks are made for background and clarity. For details see
references \[16-1, 2, 3\].

At temperatures below their Curie points Fe, Co, and Ni display
spontaneous magnetization, in which the atomic magnetic moments line up
in parallel to minimize the local energy density, i.e. they have the
maximum value for projection in a particular direction. These atomic
moments are smaller by a factor of several in the ferromagnetic material
than in the isolated atoms. A typical magnetic domain, in which the
moments are all essentially parallel, is only microns in diameter, and a
bulk material sample is usually made up of many small crystals, each
containing many such domains. A macroscopic sample's total energy also
includes contributions from small-scale fields around the domains as
well as magnetostriction and domain wall structure. The small domain
size and random orientation minimizes the total energy in the absence of
an applied field. These features can be viewed as the result of a
competition between atomic and macroscopic (but small scale) forces.
However, the domains resist reorientation, so even a very soft (easily
magnetized) material can be in a persistent magnetized state of
non-minimum energy in the absence of any external excitation. The three
room-temperature ferromagnetic elements are often alloyed with each
other and/or other elements such as Cr, Si, Al, Mn, Cu, etc. to produce
desired characteristics while retaining ferromagnetic features. A small
amount of carbon is usually also present in iron as a separate phase
mixed with the magnetic crystals.

Below 912^o^ C the stable phase of pure iron is a body-centered-cubic
crystal with Curie point 770^o^ C. The magnetic domains are usually
aligned in the directions of the three cubic axes (say
${\hat{e}}_{x},\begin{matrix}

\end{matrix}$) This material is often referred to as "ferrite" -- not to
be confused with the ferrimagnetic iron oxide ceramics of the same name,
which are used in rapidly-pulsed transformer cores. Nearly pure iron
with about .1% carbon (low carbon steel) has the combined properties of
high permeability and high saturation field desired for application with
high field solenoids. An example is 1010 steel (.08-.13% $C$,
.3-.6%$M_{n}$), which is adopted for the linac module design presented
in section 20.

When an externally-generated field is applied to a permeable material
the domains tend to line up with it, again minimizing energy density. At
low-to-moderate fields this happens by domain wall movement of favorably
oriented domains at the expense of unfavorably oriented neighbors. For
high external fields domains rotate into alignment. The process of
domain growth by wall movement is resisted by internal fiction and
occurs in tiny jumps with dissipation of energy as heat. For an
excellent qualitative discussion of this process see reference \[16-3\].
Permanent magnets are manufactured with added elements that lock in an
induced field along a preferred crystal axis (the easy axis) as they
cool from a melt; examples are SmCo~5~, Sm~2~Co~17~, and NdBFe.
Amorphous ferromagnetic alloys such as metglas are produced as
uncrystalized tape by very rapid cooling of a melt, and as such they
have essentially zero anisotropy energy. Although expensive, the
insulated tape has ideal characteristics for transformers and pulsed
cores.

**16.2 Magnetization formalism**

If in a macroscopic region of a material the mean (vector) magnetic
dipole moment per atom is
$\overrightarrow{m}\left( \overrightarrow{r} \right)$, and
$n\left( \overrightarrow{r} \right)$ is the atomic number density, then

$\overrightarrow{M}\left( \overrightarrow{r} \right) = n\left( \overrightarrow{r} \right)\begin{matrix}
\overrightarrow{m}\left( \overrightarrow{r} \right)
\end{matrix}$ (16.1)

is called the magnetization density. For a highly permeable material
this is an average over many domains. We expect
$\left| \overrightarrow{m} \right|$ to be less than the maximum average
projected moment, which is on the order of a Bohr magneton (Fe:
2.218$m_{B}$, Co: 1.714$m_{B}$, Ni: .604$m_{B}$), but this limit is
approached at high applied fields. It was mentioned in section (2) that
for solenoidal symmetry $\overrightarrow{M}$ has the form

$\overrightarrow{M} = M_{r}(r,z){\hat{e}}_{z} + M_{z}(r,z){\hat{e}}_{z}$
, (16.2)

and that it provides a current density

$\overrightarrow{J^{\text{mag}}} = \overrightarrow{\nabla}x\overrightarrow{M} = {\hat{e}}_{\theta}\left( \frac{\partial M_{r}}{\partial z} - \frac{\partial M_{z}}{\partial r} \right)$.
(16.3)

The derivation of this relation is presented in many E&M textbooks and
is not repeated here. However, its plausibility is made clear by
considering a cylinder of length$\ell$and radius $R\text{<<}\ell$, with
uniform magnetization density $M{\hat{e}}_{z}$. Then eqn. (16.3) gives

$J_{\theta}^{\text{mag}} = M\begin{matrix}
\delta(r - R)
\end{matrix}$, (16.4a)

and the total circulating current is

$I^{\text{mag}} = \ell M$. (16.4b)

We verify that this gives the correct net magnetic moment:

moment = area x current = $\left( {pR}^{2} \right)(lM)$= volume x
$\begin{array}{r}
\left( \text{magnetization} \middle| \right) \\

\end{array}$$= \sum_{\text{atoms}}^{}m_{z}$ . (16.4c)

Eqn. (2.1a) is now generalized to

$\overrightarrow{Ñ}x\overrightarrow{B} = m_{0}(\overrightarrow{J} + \overrightarrow{J^{\text{mag}}}) = m_{0}(\overrightarrow{J} + \overrightarrow{Ñ}x\overrightarrow{M})$
, (16.5)

where $\overrightarrow{J}$ denotes the contribution from wire. It is
convenient and conventional to define a new field:

$\overrightarrow{H} \equiv \frac{\overrightarrow{B}}{\mu_{0}} - \overrightarrow{M}$
, (16.6)

so we write eqn. (16.5) as

$\overrightarrow{\nabla}x\overrightarrow{H} = \overrightarrow{J}$.
(16.7)

Eqn. (2.1b), $\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$, is
still true, so from eqn. (16.6) we have

$\overrightarrow{\nabla} \cdot \overrightarrow{H} = - \overrightarrow{\nabla} \cdot \overrightarrow{M}$.
(16.8)

Therefore $\overrightarrow{H}$ may be regarded as being generated by the
combination of wire current density $\overrightarrow{J}$ and the scalar
quantity $- \overrightarrow{\nabla} \cdot \overrightarrow{M}$ (sometimes
called the magnetic pole density). Applying a theorem of vector calculus
\[16-4\], $\overrightarrow{H}$may be derived from its curl and
divergence according to

$\overrightarrow{H} = \overrightarrow{\nabla}x\int_{}^{}\frac{d^{3}r^{'}}{4\pi}\frac{\overrightarrow{J^{'}}}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|} + \overrightarrow{\nabla}\int_{}^{}{\frac{d^{3}r^{'}}{4\pi}\frac{\overrightarrow{\nabla^{'}} \cdot \overrightarrow{M^{'}}}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}}$,
(16.9)

assuming the sources vanish at large $\left| \overrightarrow{r} \right|$
. It is useful to think of $\overrightarrow{H}$ as the part of
$\overrightarrow{B}$ that is produced by distributed sources
$\overrightarrow{J}\begin{matrix}

\end{matrix}$, while $\overrightarrow{M}$ is the local contribution to
$\overrightarrow{B}$.

We have followed SI convention in defining $\overrightarrow{H}$, so in
vacuum it differs from $\overrightarrow{B}$ by the factor
$\mu_{0}^{- 1}$. Unfortunately this causes $\overrightarrow{B}$and
$\overrightarrow{H}$ to have different units; these are Tesla or
Webers/m^2^ for $\overrightarrow{B}$ and Ampere/m for
$\overrightarrow{H}$ and $\overrightarrow{M}$. If
$\left| \overrightarrow{B} \right|$ is 1.0T in vacuum then
$\left| \overrightarrow{H} \right|$is 795775 Ampere/m. Sometimes
$\overrightarrow{H}$ is defined as
$\overrightarrow{H} = \overrightarrow{B} - m_{0}\overrightarrow{M}$
instead of eqn. (16.6) to avoid this inconvenience. If Gaussian units
are used then the factor $\mu_{0}$ in eqn. (16.6) is dropped and
$\overrightarrow{B}$, $\overrightarrow{H}$, and $\overrightarrow{M}$ all
have units of gauss (1.0 gauss = 10^-4^T), but for $\overrightarrow{H}$
this unit is called the oersted. Also, for Gaussian units factors of
4$\pi$and c appear as follows:

$\begin{array}{r}
\overrightarrow{J^{\text{mag}}} = c\overrightarrow{\nabla}x\overrightarrow{M},\begin{matrix}
\overrightarrow{H} = \overrightarrow{B} - 4\pi\overrightarrow{M},
\end{matrix}\overrightarrow{\nabla}x\overrightarrow{H} = \frac{4\pi}{c}\overrightarrow{J},\begin{matrix}

\end{matrix}\overrightarrow{\nabla} \cdot \overrightarrow{H} = - 4\pi\overrightarrow{\nabla} \cdot \overrightarrow{M}\text{.}
\end{array}$

There is considerable variety in the names given for
$\overrightarrow{B}$ and $\overrightarrow{H}$, although they tend to
reflect their respective properties or uses. Some of these names are

> $\overrightarrow{B}$: magnetic field, magnetic induction, magnetic
> flux density.

$\overrightarrow{H}$: magnetic field, magnetic intensity, magnetization
field.

**16.3 B-H relations and jump conditions**

So far we have not related $\overrightarrow{H}$ to $\overrightarrow{B}$
except by its definition, eqn. (16.6). In fact for a permeable
ferromagnetic material there is no completely fixed relation except

$\left| \overrightarrow{B}/\mu_{0} - \overrightarrow{H} \right| = \left| \overrightarrow{M} \right| < M_{0}$
, (16.10)

with $M_{0}$ the magnetization density within a single domain. For
example, pure iron (natural isotopic mix) at 0 Kelvin has

$\rho = 7\text{.}\text{87}x\text{10}^{3}\text{kg}/m^{3},$ atomic mass
$= 9\text{.}\text{273}x\text{10}^{- \text{26}}\text{kg},\begin{matrix}
m/m_{B} = 2\text{.}\text{218}\text{.}
\end{matrix}$

Inserting the natural constants

$m_{B} = 9\text{.}\text{274}x\text{10}^{- \text{24}}\text{Am}^{2},\begin{matrix}
m_{0} = 4px\text{10}^{- 7}T - m/A,
\end{matrix}$

We get

$M_{0} = \frac{\rho}{\text{mass}}\frac{m}{m_{B}}m_{B} = 1\text{.}\text{746}x\text{10}^{6}\frac{A}{m},$
$m_{0}M_{0} = 2\text{.}\text{194}T\text{.}$

Unalloyed iron with very small amounts of carbon can approach this level
of magnetization. An alloy of iron and cobalt (70% $\text{Fe}$, 30%
$\text{Co}$) can increase $M_{0}$ by 10%.

The value of $\overrightarrow{M}$in a macroscopic sample actually
depends on its history of exposure to externally generated fields as
well as the instantaneous value of $\overrightarrow{H}$. By contrast a
linear relation characterizes good permanent magnet blocks:

$\overrightarrow{B} = {\hat{e}}_{\text{II}}\left( B^{\text{rem}} + m_{0}m_{\text{II}}^{r}H_{\text{II}} \right) + {\hat{e}}_{\hat{}m_{0}m_{{\hat{}}^{r}H_{\hat{}}}}$,
(16.11)

where $\text{II}$and $$refer to directions along and normal to the easy
axis of magnetization. Here $B^{\text{rem}}$, called the remenant field,
is simply the "permanent" value of $m_{0}M$in a long sample, and the
relative permeabilities $\mu_{\text{II}}^{r}$ and $\mu_{}^{r}$ are
constants out to large values of $H$ and typically exceed unity by less
than 10%. For very large reverse $m_{0}H$ (depending on material but
ideally of order $B^{\text{rem}}$ or greater) a block's properties
become nonlinear and change irreversibly.

It is apparent that the commonly-assumed linear, isotropic relation

$\overrightarrow{B} = \mu\overrightarrow{H}$ (16.12)

is of limited validity except for paramagnetic and diamagnetic
substances at low fields. Nevertheless its use is sometimes a
justifiable approximation. For example cast iron excited from an initial
state of no remanant field exhibits such a linear relationship out to
$\left| \overrightarrow{B} \right| \approx \text{.}3T$, with
$m»\text{400}m_{0}$. Some types of steel have $\mu$greater than
$\text{10}^{4}\mu_{0}$. Confusing the issue, eqn. (16.12) is often
regarded as a definition of $\mu$for whatever $\overrightarrow{B}$ and
$\overrightarrow{H}$are present, as long as they are parallel. If we
assume the validity of eqn. (16.12) then the following definitions are
made:

$K^{m} \equiv \mu^{r} \equiv \mu/\mu_{0}$ (16.13)

is the relative permeability, and

$\chi^{m} = K^{m} - 1$ (16.14)

is the magnetic susceptibility. From the definition of
$\overrightarrow{H}$ we then have

$\overrightarrow{M} = \frac{\overrightarrow{B}}{\mu_{0}} - \overrightarrow{H} = K^{m}\overrightarrow{H} - \overrightarrow{H} = \chi^{m}\overrightarrow{H}$.
(16.15)

An improvement on (16.12), used in some computations, is

$\overrightarrow{B} = \overrightarrow{B_{0}} + m\left( \overrightarrow{H} - \overrightarrow{H_{0}} \right)$,
(16.16)

where $\overrightarrow{B_{0}}$ and $\overrightarrow{H_{0}}$are some
particular values and $\mu$ is a local coefficient.

For many computations it is simply assumed that
$\mu/\mu_{0} \rightarrow \infty$, i.e. $\overrightarrow{H}$ vanishes in
the highly permeable material. However, to compute $\overrightarrow{B}$
inside the material requires specification of details about the
permeability, even though it approaches infinity -- see section (18).

If $\overrightarrow{B}$and $\overrightarrow{H}$ are to be determined by
solution of the coupled underlying differential equations, then jump
conditions at vacuum-material interfaces are required. From eqns. (16.7)
and (2.1b) these are

$H_{t}$ is continuous, (16.17a)

$B_{n}$ is continuous, (16.17b)

where $t$ and $n$ refer to components tangential and normal to the
interface.

These jump conditions suggest the utility of a highly permeable material
for flux control or field shaping. Since $\mu_{0}\overrightarrow{H}$can
essentially vanish in the material, its tangential component on the
vacuum side of the boundary also vanishes, as does the vacuum tangential
component of $\overrightarrow{B}$. But the vacuum field can be derived
from a scalar potential,
$\overrightarrow{B} = m_{0}\overrightarrow{Ñ}f$, which is single valued
around the material if no wire is enclosed. Therefore $\varphi$ is
constant on the material surface. By analogy with electrostatics, the
material acts like a conducting surface, and strong normal fields may be
localized to particular locations where the surface is strongly curved.
Roughly speaking the field is channeled in and out of the material where
desired, and routed away from locations (such as an induction core)
where it could cause a problem. It is also apparent that stored energy
might be reduced in this way since energy density in the material is
roughly $\text{BH}/2 \approx 0$, although the concentration of flux in
particular locations outside the material tends to increase stored
energy density. The design of the poles and yoke of a cyclotron is an
application of these ideas.

Field calculations for highly permeable materials are complicated by
both non- linearity of $B - H$ curves and the phenomenon of hysteresis.
The latter is closely related to internal dissipation of energy during
excitation and de-excitation. Following the textbooks, we imagine that
in a macroscopic region of material $\overrightarrow{B}$ and
$\overrightarrow{H}$ are parallel (more about this assumption below) and
$H$ is increasing or decreasing slowly with time. Depending on the sign
of $\text{dH}/\text{dt}$ two different curves are generated. If $H$ is
cyclical, a closed loop in the $B - H$ plane is produced. If
sufficiently large values of $H$ are in the cycle, so that the material
is pushed well into saturation, then a maximal hysteresis loop is
generated, and all other loops lie somewhere inside of it. The maximal
hysteresis loop is the boundary of accessible points in the $B - H$plane
under the assumptions of parallelism and slow variation. This is
considerably more restrictive than the condition (16.10). The intercepts
of the maximal hysteresis curve with the $B$ and $H$ axes are called
respectively the remanence $\left( \pm B_{r} \right)$ and coercive force
$\left( \pm H_{c} \right)$. Some loops have zero bias current, i.e. no
mean $H$, and form a nested set around the origin $(B = H = 0)$. Other
(biased) loops are offset from the origin and may or may not encircle
it. In this picture we assume that the magnet is excited slowly enough
that eddy currents are not appreciable and that magnetic domains have
time to adjust their walls or orientations. All closed loops must cycle
in a counter-clockwise direction to avoid a violation of the second law
of thermodynamics. To see this, note that the differential work done by
an external circuit is $\text{HdB}$ per unit volume. A clockwise loop
would extract energy into the circuit and cool the material (eliminating
any need for controlled fusion). Through every point $(B,H)$ interior to
the maximal loop we may draw a pair of paths, which are followed
depending on whether $H$ is increasing or decreasing. A grid of crossing
paths covers the accessible $B - H$ space, and any point may be reached
by traversing several paths in their allowed directions. Complete
demagnetization can be done by cycling$H$with gradually decreasing
amplitude.

**16.4 Magnetization curve**

Since ferromagnetic materials have a nonlinear magnetization response
when driven over a large range, a generalization of eqn. (16.12) is
required for computations. Usually it is assumed that a material is
excited from an initially unmagnetized state and that $H$ only increases
in this process, so there is no ambiguity from hysteresis. Then
$B$follows a known "magnetization curve", i.e. $B = F(H)$ - see the
accompanying figure for examples (from Chao and Tigner, Handbook of
Accelerator Physics and Engineering, World Scientific, 1999.)
Computations may be made by setting

$\overrightarrow{B} = \overrightarrow{H}F\left( \overrightarrow{r},\left| \overrightarrow{H} \right| \right)/\left| \overrightarrow{H} \right|$
, (16.18)

which assumes$\overrightarrow{B}$and $\overrightarrow{H}$ remain
parallel everywhere. This assumption works if the response is linear,
i.e. $F = m\left( \overrightarrow{r} \right)H$. To demonstrate this we
write eqn. (16.9) in the form

$H_{i} = H_{i}^{J} + \int_{}^{}\frac{d^{3}r^{'}}{4p}M_{j}^{'}\frac{¶^{2}}{¶x_{i}¶x_{j}}\frac{1}{\left| \overrightarrow{r^{'}} - \overrightarrow{r} \right|}$
, (16.19)

with subscripts $(i,j)$ denoting Cartesian components with a sum over
$j$, and $\overrightarrow{H^{J}}$is the part of $\overrightarrow{H}$
produced by currents in wire. For the linear isotropic response we have

$\overrightarrow{M}\left( \overrightarrow{r^{'}} \right) = c^{m}\left( \overrightarrow{r^{'}} \right)\overrightarrow{H}\left( \overrightarrow{r^{'}} \right)$
, (16.20)

with $\chi^{m} = \mu/\mu_{o} - 1$ , and eqn. (16.19) can be reduced to a
set of coupled linear equations by dividing the entire space into small
volumes (finite element method). $\overrightarrow{H}$is then obtained
from $\overrightarrow{H^{J}}$by a single matrix inversion, and the
assumption of parallelism has lead to a unique solution. For a
non-linear magnetization curve, we can use

$\overrightarrow{M} = \frac{\overrightarrow{B}}{m_{0}} - \overrightarrow{H} = \left( \frac{F\left( \overrightarrow{r},\left| \overrightarrow{H} \right| \right)}{m_{0}\left| \overrightarrow{H} \right|} - 1 \right)\overrightarrow{H}$,
(16.21)

and it is necessary to iterate to solve eqn. (16.19). Unfortunately it
is not clear that this procedure leads to the physically correct
solution, although this is generally assumed. Consider a system that is
excited from the un-magnetized state by slowly increasing the wire
current. If the material response is linear then $\overrightarrow{B}$and
$\overrightarrow{H}$ simply scale up proportional to current without
changing direction. But the situation is changed when the material
saturates and
$\overrightarrow{B} \rightarrow \mu_{0}\overrightarrow{H}$; fields then
approximate a vacuum pattern determined from the wire current alone, and
in general they have a different direction than when the permeable
material is effective. During the current increase $\overrightarrow{B}$
and $\overrightarrow{H}$ have changed direction; but this is different
from the simple experimental conditions in which magnetization curves
are determined. It is conceivable that during a gradually changing
direction of $\overrightarrow{H}$, that $\overrightarrow{M}$ (and
therefore $\overrightarrow{B}$) stays parallel, but this is an
additional material property to be determined by experiment.

A mathematical model of hysteresis has been described in an interesting
series of reports - see \[16-5\]. Recall that every $B - H$point inside
the limiting cycle is crossed by two curves with slopes
$\text{dB}/\text{dH}$ depending on whether $H$ is increasing or
decreasing. These slopes are found to be well approximated by the
assumed formula

$\frac{\text{dB}}{\text{dH}} = \pm \alpha\left\lbrack f(H) - B \right\rbrack + g(H)$
, (16.22)

where $\alpha$is a constant and $(f,g)$ are material-dependent functions
with examples given in the cited report. The choice
$\pm \rightarrow + \begin{matrix}
\text{or -}
\end{matrix}$ depends respectively on $H$ increasing or decreasing. A
generalization that includes isotropic permanent magnets is written

> $\frac{\text{dH}}{\text{dB}} = \pm \alpha\left\lbrack \overset{\sim}{f}(B) - H \right\rbrack + \overset{\sim}{g}(B)$.
> (16.23)

The function $\overset{\sim}{g}$ can be modified to include dependence
on the rate of magnetization $\left( \text{dB}/\text{dt} \right)$. If
$\alpha$, $f$and $g$are known, then the magnetization curve $F(H)$ can
be found in numerical form by integrating eqn. (16.22) from the origin
$(B = H = 0)$with increasing $H$:

> $\frac{\text{dB}}{\text{dH}} + \alpha B = \alpha f(H) + g(H)$ ,
> (16.24)
>
> $F(H) = B = \int_{}^{}{dH^{'}}\left\lbrack \alpha f\left( H^{'} \right) + g\left( H^{'} \right) \right\rbrack e^{\alpha\left( H^{'} - H \right)}$
> . (16.25)

**16.5 Permanent magnet solenoid**

Finally we consider the simple permanent magnet solenoid consisting of
an annulus of constant magnetization $M{\hat{e}}_{z}$between radii
$R_{1}$and $R_{2}$ and with length $\ell$, centered at $z = 0$. The
magnetization currents are simply a pair of thin layers:

> $J_{q}^{\text{mag}}\text{=-}\frac{¶M_{z}}{¶r} = M\left\lbrack d\left( r - R_{2} \right) - d\left( r - R_{1} \right) \right\rbrack$
> . (16.26)

Applying the thin layer result, eqn. (4.5), we have the on-axis field

> $B_{0}(z) = \frac{\mu_{0}M}{2}\begin{array}{r}
> \left\lbrack \frac{(z + \ell/2)}{\sqrt{(z + \ell/2)^{2} + R_{2}^{2}}} - \frac{(z - \ell/2)}{\sqrt{(z - \ell/2)^{2} + R_{2}^{2}}} \middle| \right\rbrack \\
>
> \end{array}$ . (16.27)

This function is plotted below for $$. Note that $B_{0}$ is negative in
the solenoid bore and positive outside, and its maximum absolute value
is only $\text{.}\text{19}\mu_{0}M$. When such a solenoid is used as a
particle beam lens, its focal strength is proportional to the integral
of $B_{0}^{2}$, so it is not a very effective use of material, although
such lenses have been used to confine electron beams in Klystrons. It
can be easily verified that

$\int_{}^{}\text{dz}B_{0}(z) = 0$. (16.28)

In fact this relation is true for any solenoid made only of magnetic
materials since
$\overrightarrow{\nabla}x\overrightarrow{H} = \overrightarrow{J}$
requires that the line integral of $\overrightarrow{H}$ around any
closed loop must equal the enclosed wire current:

$\oint_{}^{}d\overrightarrow{r} \cdot \overrightarrow{H} = I = 0$ .
(16.29)

In the present case we may take the integration loop along the entire
axis and return at large radius, where there is a vanishing
contribution.

![](media/image14.pct){width="5.997222222222222in" height="7.48125in"}

mu0M = 1.0;

R2 = .5;

R1 = .25;

l = 1.0;

B0\[z\_\] = mu0M/

2\*((z + l/2)/((z + l/2)\^2 + R2\^2)\^.5 - (z -

l/2)/((z - l/2)\^2 + R2\^2)\^.5 - (z +

l/2)/((z + l/2)\^2 + R1\^2)\^.5 + (z -

l/2)/((z - l/2)\^2 + R1\^2)\^.5);

Plot\[B0\[z\], {z, -2, 2}\];

![](media/image15.png){width="4.0in" height="2.4583333333333335in"}

**References**

16-1 Bozorth, R.M., Ferromagnetism, IEEE Press, 1993.

16-2 Kittel, C., Introduction to Solid State Physics, 7^th^ ed., ch. 15,
John Wiley & Sons,

1996\.

16-3 Feynman, R., Leighton, R., and Sands, M., The Feynman Lectures on
Physics, Vol

II, ch 37, Addison -- Wesley, 1964.

16-4 Panofsky, W., and Phillips, M., Classical Electricity and
Magnetism, 2^nd^ ed., pp 2-5,

Addison -- Wesley, 1962.

16-5 Hodgdon, M., Applications of a Theory of Ferromagnetic Hysteresis,
IEEE

Transactions on Magnetics, Vol. 24, No 1, January 1988. 218-221.

**17.** **[Wire Layer Embedded in a Highly Permeable Yoke]{.underline}**

A simple geometry that displays flux confinement is a wire layer (of any
shape) that has net azimuthal current $I$ and is embedded in a long
permeable annulus with $m\text{=¥}$. The inner radius of the
annulus/yoke is at $r = b$, which is also the magnet's bore radius, and
a narrow slot connects the wire layer to the bore at $z = 0$ (see
figure).

![](media/image16.pct){width="5.666666666666667in"
height="2.4722222222222223in"}

Inside the bore, where

$\overrightarrow{H} = \frac{\overrightarrow{B}}{m_{0}} = \overrightarrow{Ñf}$,
(17.1a,b)

we have

$Ñ^{2}f = 0,\begin{matrix}
 & & 
\end{matrix}Ñ^{2}B_{z} = 0$ . (17.2a,b)

Along the inner surface of the annulus $f\begin{matrix}

\end{matrix}$is constant and $B_{z}$ vanishes, except at $z = 0$, where
$f\begin{matrix}

\end{matrix}$jumps and $B_{z}$ is infinite:

$f(b,z) = \begin{array}{r}
\left\{ - I/2\begin{matrix}
 & 
\end{matrix}z < 0 \middle| \right.\ 
\end{array}$ (17.3)

$B_{z}(b,z) = m_{0}Id(z)$ . (17.4)

Here the annulus is assumed to extend to $z = \pm ¥$ for computations.

To solve for $B_{z}$ in the magnet bore we apply a Fourier integral
transform to eqn. (17.2 b):

${\overset{\sim}{B}}_{z}(r,k) = \int_{}^{}\text{dz}\begin{matrix}
e^{- \text{ikz}}B_{z}(r,z)
\end{matrix}$ , (17.5)

$\frac{1}{r}\frac{¶}{¶r}r\frac{¶{\overset{\sim}{B}}_{z}}{¶r} - k^{2}{\overset{\sim}{B}}_{z} = 0$
. (17.6)

Equation (17.6) is satisfied by

${\overset{\sim}{B}}_{z} = f(k)I_{0}\left( |k|r \right)$ , (17.7)

where $I_{0}$ is the modified Bessel function and $f(k)$ is to be
determined from the boundry condition (17.4). Since we may write the
delta function as

$d(z) = \int_{}^{}\text{dk}\frac{e^{\text{ikz}}}{2p}$ , (17.8)

we have immediately from equations (17.4) and (17.7)

$B_{z} = \int_{}^{}\text{dk}\frac{e^{\text{ikz}}}{2p}{\overset{\sim}{B}}_{z} = m_{0}I\int_{}^{}\text{dk}\frac{e^{\text{ikz}}}{2p}\frac{I_{0}\left( |k|r \right)}{I_{0}\left( |k|b \right)}$
. (17.9)

A convenient form for computations, that displays scaling with $rb$ and
$zb$ is

$B_{z} = \frac{m_{0}I}{pb}\int_{}^{}{\text{dy}\frac{\text{cos}(yz/b)I_{0}\left( \text{yr}/b \right)}{I_{0}(y)}}$
, (17.10)

where $yº\text{kb}$. The radial field is obtained from eqn. (17.9) by
integrating $\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$:

$B_{r}\text{=-}m_{0}I\int_{}^{}\text{dk}\frac{k}{|k|}i\frac{e^{\text{ikz}}}{2p}\frac{I_{1}\left( |k|r \right)}{I_{0}\left( |k|b \right)}$
(17.11)

$\begin{matrix}
 & & 
\end{matrix} = \frac{m_{0}I}{pb}\int_{}^{}\text{dy}\begin{matrix}

\end{matrix}\text{sin}\left( \text{yz}/b \right)\frac{I_{1}\left( \text{yr}/b \right)}{I_{0}(y)}$
. (17.12)

The bore field is displayed for several values of $rb$ at the end of
this section. A strong (exponential) falloff with $|z|/b$ is apparent
and contrasts sharply with the $|z|^{- 3}$ fall off from a lens without
a yoke. This calculation is formally equivalent to that of an
electrostatic field in a conducting pipe that is split at $z = 0$ and
with a potential difference applied to the two sides.

The potential $f\begin{matrix}

\end{matrix}$can be derived from eqn. (17.9) by integration in $z$,
however at this point it is more instructive to use an expansion in
ordinary Bessel functions:

$f = \frac{I}{2} + \sum_{i = 1}^{¥}A_{i}\begin{matrix}
J_{0}\left( x_{i}rb \right)e^{- x_{i}z/b}
\end{matrix}$ , (17.13)

valid for $z > 0$. For $z < 0,$ $f\begin{matrix}

\end{matrix}$is determined by its antisymmetry:

$f = \left( r,\overset{}{-}z \right)\text{=-}f(r,z)$ , (17.14)

so the there is exponential decay in both directions. The$c_{i}$ are the
zeros of $J_{0}$, and the coefficients $A_{i}$ are determined from the
condition

$0 = f(r,0) = \frac{I}{2} + \sum_{i = 1}^{¥}A_{i}J_{0}\left( x_{i}r/b \right)$
; (17.15)

$A_{i}\text{=-}\frac{I}{x_{i}J_{1}\left( x_{i} \right)}$ . (17.16)

We have for $z > 0$,

$f = \frac{I}{2} - I\sum_{i = 1}^{¥}{}\frac{J_{0}\left( x_{i}rb \right)e^{- x_{i}zb}}{x_{i}J_{1}\left( x_{i} \right)}$
, (17.17)

$B_{z} = m_{0}\frac{\text{¶f}}{¶z} = m_{0}\frac{I}{b}\sum_{i - 1}^{¥}{}\frac{J_{0}\left( x_{i}rb \right)e^{- x_{i}zb}}{J_{1}\left( x_{i} \right)}$
. (17.18)

The first five values of $x_{i}$ and $J_{1}\left( x_{i} \right)$ are
tabulated to aid in computations:

$x_{i}$$J_{1}\left( x_{i} \right)$

$i$ $$

$$

The on-axis field can be calculated from either the Fourier integral or
Bessel series, although the latter is poorly convergent for $|z|/b»0$.
The leading term of the series is a good approximation for $|z|/b³1$:

$B_{0} \rightarrow \frac{m_{0}I}{b}\frac{e^{- x_{1}|z|/b}}{J_{1}\left( x_{1} \right)}$
, (17.19)

but is large by a factor of 2.9 at $z = 0$. An excellent approximation
for all $z$ is

$B_{0}»\frac{m_{0}I}{b}\frac{1}{J_{1}\left( x_{1} \right)}\frac{1}{\left\lbrack 2\text{cosh}\left( x_{1}|z|/\text{sb} \right) \right\rbrack^{S}}$
(17.20)

with $s =$`<!-- -->`{=html}1.538781. This has maximum relative error of
$\pm$`<!-- -->`{=html}2x10^-4^ at
$|z|/b = \left( 0, \pm \text{.}7b \right)$.

> ![](media/image17.png){width="4.0in" height="2.4583333333333335in"}

$B_{z}(r,z)\begin{matrix}

\end{matrix}$

**18.** **[Periodic Wire Layer with Highly a Permeable
Yoke]{.underline}**

Return flux can be channeled by a permeable annulus if the field there
is kept below saturation
$\left( B_{\text{sat}}»1\text{.}\text{0-2}\text{.}\text{0T for  iron} \right)$.
For a periodic system of solenoids all having the same polarity, the
field in the annulus or yoke is roughly

$B_{\text{yoke}} \approx - \frac{R^{2}B_{\text{bore}}\ell/P}{R_{2}^{2} - R_{1}^{2}}$
, (18.1)

where $R$ is the wire radius, $R_{1}$and $R_{2}$ are the yoke radii,
$\ell/P$ is the ratio of wire layer length to period length, and
$B_{\text{bore}}$ is the peak bore field. The volume of the yoke should
be at least $B_{\text{bore}}/B_{\text{sat}}$ times the total volume of
magnet bore to prevent stray flux from entering the outside world. In
practice there will be gaps in the yoke, so some flux will bulge out
anyway.

The field in the presence of a highly permeable yoke can be examined
with the simple periodic model:

$J_{q} = \begin{array}{r}
\left\{ Sd(r - R)\begin{matrix}
 & & 
\end{matrix} - l/2 < z < l/2, \middle| \right.\ 
\end{array}$ (18.2)

$m = \begin{array}{r}
\left\{ ¥\begin{matrix}
 & & & & 
\end{matrix}R_{1} < r < R_{2}, \middle| \right.\ 
\end{array}$ (18.3)

with $R < R_{1}$ and $J_{q}$ repeating with period P. The Fourier
expansion of the current density is

$J_{q} = d(r - R)\left\lbrack \frac{Sl}{P} + \sum_{n = 1}^{¥}d_{n}\text{cos}\left( k_{n}z \right) \right\rbrack ºd(r - R)S(z)$
, (18.4)

with

$d_{n} = \frac{2S}{pn}\text{sin}\left( \frac{pnl}{P} \right)$ ,
$k_{n} = \frac{2pn}{P}$ . (18.5a,b)

In the yoke $\overrightarrow{H} = \overrightarrow{B}/\mu \rightarrow 0$
by assumption, and since the tangential components of
$\overrightarrow{H}$ at $R_{1}$ and $R_{2}$are continuous, we conclude
the$B_{z} = \mu_{0}H_{z} = 0$ at $r = R_{1 -}$and$r = R_{2 +}$. Since
$\overrightarrow{B}$ must vanish as $r \rightarrow \infty$, it must also
vanish for all $r > R_{2}$. For $r < R_{1}$we must simply solve

$Ñ^{2}B_{z}\text{=-}m_{0}\frac{1}{r}\frac{¶}{¶r}\text{rJ}_{q}$ , (18.6)

with the boundary condition $B_{z}\left( R_{1 -} \right) = 0$.

The value of $\overrightarrow{B}$inside the yoke is indeterminate unless
more information is given about the permeability, which must be assumed
to be not quite infinite. For example, suppose we have a non-linear
magnetization curve

$\overrightarrow{B} = \frac{\overrightarrow{H}}{\left| \overrightarrow{H} \right|}F\left( \left| \overrightarrow{H} \right| \right)$
, (18.7a)

$\overrightarrow{H} = \frac{\overrightarrow{B}}{\left| \overrightarrow{B} \right|}F^{- 1}\left( \left| \overrightarrow{B} \right| \right)$
. (18.7b)

Then in the yoke we must solve
$\overrightarrow{Ñ} \cdot \overrightarrow{B} = 0$ along with

$\overrightarrow{Ñ}x\left\lbrack \frac{\overrightarrow{B}}{\left| \overrightarrow{B} \right|}F^{- 1}\left( \left| \overrightarrow{B} \right| \right) \right\rbrack = \overrightarrow{Ñ}x\overrightarrow{H} = 0$
. (18.8)

The boundary conditions for the yoke interior are that $B_{r}$be
continuously joined to its vacuum values at $R_{1}$and $R_{2}$ and that
the net flux of the system vanish. For the simple assumption that $\mu$
is very large but constant, the Cartesian components of the yoke field
satisfy

$Ñ^{2}\overrightarrow{B} =$$\overrightarrow{Ñ}\overrightarrow{Ñ} \cdot \overrightarrow{B}$$- \overrightarrow{Ñ}x\left( \overrightarrow{Ñ}x\overrightarrow{B} \right)\text{=-}\overrightarrow{Ñ}x\left( \overrightarrow{Ñ}xm\overrightarrow{H} \right) = 0$
. (18.9)

Returning to the analytical model and assuming $B/H$is infinite but
constant in the yoke, we have inside all four radial zones

$\nabla^{2}B_{z} = 0$. (18.10)

At zone boundaries $\left( R,R_{1},R_{2} \right)$ both $B_{r}$ and
$¶B_{z}/¶r$are continuous, while $B_{z}$ jumps across the wire:

$B_{z}\left( R_{+},z \right) - B_{z}\left( R_{-},z \right)\text{=-}m_{0}S(z)$
. (18.11)

At the yoke boundaries it is also found that $B_{z}$ jumps from zero in
vacuum to a finite value in the material that is found by solving the
model equations. We proceed by solving for the field produced by
individual terms of the expansion eqn. (18.4) for $J_{\theta}$. First,
the term $S\ell/P$produces

$B_{z} \rightarrow \begin{array}{r}
\left\{ \frac{m_{0}Sl}{P}\begin{matrix}
\begin{matrix}
\begin{matrix}
\begin{matrix}
\begin{matrix}
 & 
\end{matrix}
\end{matrix}
\end{matrix}
\end{matrix} & & & & 
\end{matrix}0 < r < R, \middle| \right.\ \left\{ \begin{matrix}
 & & 
\end{matrix}0\begin{matrix}
 & & & & 
\end{matrix}R < r < R_{1}, \middle| \right.\ \left\{ - \frac{R^{2}}{R_{2}^{2} - R_{1}^{2}}\frac{m_{0}Sl}{P}\begin{matrix}
 & & 
\end{matrix}R_{1} < r < R_{2}, \middle| \right.\ 
\end{array}$ (18.12)

This function has net zero flux, vanishes at $r = \infty$, and has no
on-axis contribution from the permeable material (as discussed in
section (16)). Next consider the field from any periodic term of eqn.
(18.4); this has the form

$S(z) \rightarrow d_{n}\text{cos}\left( k_{n}z \right)$ , (18.13a)

$B_{z} \rightarrow m_{0}f_{n}(r)d_{n}\text{cos}\left( k_{n}z \right)$,
(18.13b)

$B_{r} \rightarrow m_{0}g_{n}(r)d_{n}\text{sin}\left( k_{n}z \right)$.
(18.13c)

Inside each zone $f_{n}$ satisfies

$\frac{1}{r}\frac{d}{\text{dr}}r\frac{\text{df}_{n}}{\text{dr}} - k_{n^{2}}f_{n} = 0$,
(18.14)

and from $\overrightarrow{Ñ}´\overrightarrow{B} = 0$ , inside any zone

$$ $g_{n} = \frac{1}{k_{n}}\frac{¶f_{n}}{¶r}$. (18.15)

Equation (18.14) is a modified Bessel equation, so $f_{n}(r)$is a linear
combination of $I_{0}\left( k_{n}r \right)$ and
$K_{0}\left( k_{n}r \right)$ in every zone, and $g_{n}(r)$ is the
combination of $I_{1}\left( k_{n}r \right)$ and
$K_{1}\left( k_{n}r \right)$ determined from eqn. (18.15). Applying the
various boundary and jump conditions and finiteness at $r = 0$ we find,
suppressing the subscript $(n)$:

$$0 < r < R$$

$f(r) = \text{AI}_{0}\left( \text{kr} \right)$, (18.16a)

$g(r) = \text{AI}_{1}\left( \text{kr} \right)$ , (18.16b)

$R < r < R_{1}$

$f(r) = B\left\lbrack \frac{I_{0}\left( \text{kr} \right)}{I_{0}\left( \text{kR}_{1} \right)} - \frac{K_{0}\left( \text{kr} \right)}{K_{0}\left( \text{kR}_{1} \right)} \right\rbrack$,
(18.16c)

$g(r) = B\left\lbrack \frac{I_{1}\left( \text{kr} \right)}{I_{0}\left( \text{kR}_{1} \right)} + \frac{K_{1}\left( \text{kr} \right)}{K_{0}\left( \text{kR}_{1} \right)} \right\rbrack$,
(18.16d)

$R_{1} < r < R_{2}$

$f(r) = C\left\lbrack \frac{I_{0}\left( \text{kr} \right)}{I_{1}\left( \text{kR}_{2} \right)} + \frac{K_{0}\left( \text{kr} \right)}{K_{1}\left( \text{kR}_{2} \right)} \right\rbrack$,
(18.16e)

$g(r) = C\left\lbrack \frac{I_{1}\left( \text{kr} \right)}{I_{1}\left( \text{kR}_{2} \right)} - \frac{K_{1}\left( \text{kr} \right)}{K_{1}\left( \text{kR}_{2} \right)} \right\rbrack$,
(18.16f)

$$R_{2} < r < \infty$$

$f(r) = g(r) = 0$, (18.16g)

with coefficients

$A = \frac{\text{kR}}{I_{0}\left( \text{kR}_{1} \right)}\left\lbrack I_{1}\left( \text{kR} \right)K_{0}\left( \text{kR}_{1} \right) + K_{1}\left( \text{kR} \right)I_{0}\left( \text{kR}_{1} \right) \right\rbrack$,
(18.17a)

$B = \text{kRI}_{1}\left( \text{kR} \right)K_{0}\left( \text{kR}_{1} \right)$,
(18.17b)

$C = B\frac{\text{RI}_{1}\left( \text{kR} \right)}{R_{1}I_{0}\left( \text{kR}_{1} \right)}\left\lbrack \frac{I_{1}\left( \text{kR}_{1} \right)}{I_{1}\left( \text{kR}_{2} \right)} - \frac{K_{1}\left( \text{kR}_{1} \right)}{K_{1}\left( \text{kR}_{2} \right)} \right\rbrack^{- 1}$.
(18.17c)

We have used the Bessel function identity

> $I_{1}(x)K_{0}(x) + I_{0}(x)K_{1}(x) = \frac{1}{x}$ (18.18)

to simplify these expressions.

Plots.$B_{z}$ vs $z$ for various $r$ are given below for the case $$.
The suppression of the external field in the zone $R < r < R_{1}$ is
evident -- compare section (10).

![](media/image18.png){width="4.0in" height="2.4583333333333335in"}

![](media/image19.png){width="4.0in" height="2.4583333333333335in"}

![](media/image20.png){width="4.0in" height="2.4583333333333335in"}

### 19. [Tosca^©^ Model]{.underline}

Many computer programs have been written for the purpose of solving
Laplace\'s or Poisson\'s equation under given boundary conditions. A
number of different numerical solution methods exist for solving a
system of partial differential equations; one of the most commonly
applied techniques is the finite element discretization method. Finite
element analysis (FEA) requires special enhancements, which are
described briefly in this chapter, to make it applicable to
electromagnetic field calculations.

Computer programs most widely used in the accelerator engineering
community are the public domain code POISSON \[[^1]-7\], and the
commercial code Tosca^©^ \[1-7\]. POISSON can simulate only 2d plane or
axisymmetric geometries, whereas Tosca^©^ allows calculations in 2D and
3D. Since today\'s personal computers have become sufficiently powerful,
the fidelity of these codes is so superior that often simulations are
used to check measurements, and large magnetic devices are designed
without the need for building prototypes. For instance, a field accuracy
of the order of parts in 10^6^ is required for simulating shielded
superconducting NMR systems. Another application is to cyclotrons with
their stringent field accuracies (parts in 10^5^), which can nowadays be
entirely designed using Tosca^©^ alone.

FEA codes are especially well suited for modeling magnetic field
problems with non-linear materials present, e.g. magnet steel or
permanent magnets. Forces on conductors and pole pieces can be easily
calculated. In addition, harmonic or transient field problems with eddy
currents present can be solved. Eddy current heating as well as magnetic
forces can be coupled to thermal and structural FEA simulations for the
engineering design.

### 19.1 General solution strategy in 3D

In this section the notation for fields differs somewhat from that used
in the rest of the report in order to correspond to Tosca^©^
conventions. A detailed description of how to numerically solve
Maxwell\'s equations is beyond the scope of this report. Nevertheless,
since the computation of magnetic fields requires a few special numeric
treatments, its peculiarities are briefly mapped out here. We want to
discretize Maxwell\'s equations:

> $\overrightarrow{\nabla}x\overrightarrow{E} = - \frac{\partial\overrightarrow{B}}{\partial t}$
> , (19.1a)
>
> $\overrightarrow{\nabla}x\overrightarrow{H} = \frac{\partial\overrightarrow{D}}{\partial t} + \overrightarrow{j}$
> , (19.1b)
>
> $\overrightarrow{\nabla} \cdot \overrightarrow{D} = \rho$ , (19.1c)
>
> $\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$ , (19.1d)

with
$\overrightarrow{D} = e\overrightarrow{E},\overrightarrow{B} = m\overrightarrow{H},\text{ and  \{}\overrightarrow{j} = s\overrightarrow{E}$
(or is specified). In these equations $\overrightarrow{E}$and
$\overrightarrow{H}$ are the electric and magnetic fields and
$\overrightarrow{D}$and $\overrightarrow{B}$are the electric and
magnetic flux densities. The current density is denoted by
$\overrightarrow{j}$, and $\rho$ denotes the charge density. The three
material functions $e,m,\text{ and }s$(permittivity, permeability, and
conductivity) are in general functions of the spatial coordinates and
also nonlinear functions of the electromagnetic field strength, e.g.
saturation effects in metal. All three material properties may be
tensors, as in laminated transformer cores and permanent magnets.

![](media/image21.png){width="5.875196850393701in"
height="4.148425196850393in"} The usual procedure for solving Maxwell\'s
equations is to simplify them as far as possible, first setting some of
the quantities $e,\begin{matrix}

\end{matrix}$ equal to constants if possible, and then deducing a second
order differential equation. This second-order differential equation for
some intermediate function (typically a vector potential, scalar
potential, or both) can then be discretized and numerically solved for
realistically shaped structures by means of finite difference, finite
element, or integral methods. The finite element method is based on
division of the domain of these equations (volume of space in which the
equations are satisfied) into small volumes (the finite elements) as
shown in Figure 1. Within each finite element a simple polynomial is
used to approximate the solution.

Consider first a Poisson equation describing an electric potential
function $F$:

> $\overrightarrow{Ñ} \cdot \overrightarrow{E}\text{=- \{}\overrightarrow{Ñ} \cdot \overrightarrow{Ñ}F = \frac{\rho}{e_{0}}$
> . (19.2)

Such an equation,

> $- (\frac{¶^{2}}{¶x^{2}}F + \frac{¶^{2}}{¶y^{2}}F + \frac{¶^{2}}{¶z^{2}}F) = \frac{\rho}{e_{0}}$
> , (19.3)

can be easily discretized for computations. In order to define $F$,
boundary conditions are required; these may be either assigned values of
$F$ or its normal derivative $\partial\Phi/\partial x$ on a surface. In
all electrostatic field examples it is essential that the potential is
defined at least at one point in the domain, otherwise an infinite
number of solutions could be generated by adding an arbitrary constant
to a particular solution.

The situation is slightly different for magnetic problems. In the
magnetostatic limit, Maxwell\'s equations reduce to

> $\overrightarrow{\nabla}x(\frac{\overrightarrow{B}}{\mu}) = \overrightarrow{j}$
> , (19.4a)
>
> $\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$ . (19.4b)

Since $\overrightarrow{B}$ is not \"curl free\", it follows that the
magnetic flux density cannot be represented by the gradient of a scalar
potential as in the electrostatic case. It is convenient to split the
total field into two parts, a conductor source field and a gradient of a
scalar potential, in order to obtain a description of the field in terms
of a simple scalar potential:

> $\overrightarrow{H} = {\overrightarrow{H}}_{s} - \overrightarrow{\nabla}\Phi$
> with $\overrightarrow{\nabla} \cdot {\overrightarrow{H}}_{s} = 0$ .
> *Reduced scalar potential* (19.5)

The conductor source field intensity can be evaluated separately, for
example by the

Biot-Savart law

> $$ . (19.6)

Introducing the scalar permeability $m$ and combining
$\overrightarrow{\nabla} \cdot \overrightarrow{B} = 0$ with the above
two equations, we have the partial differential equation for $F$:

> $\overrightarrow{\nabla} \cdot \mu\overrightarrow{\nabla}\Phi - \overrightarrow{\nabla} \cdot \mu{\overrightarrow{H}}_{s} = 0$
> , (19.7)

with ${\overrightarrow{H}}_{s}$ calculated from eqn. (19.6).

The splitting of the magnetic field into a \"curly\" source part and a
reduced scalar potential part has finally reduced the problem to a
purely scalar one. Equation (19.7), like the Poisson equation for
electrostatic fields, can be solved using the finite element method.
Unfortunately, this method leads to large computational errors,
especially in volumes where ${\overrightarrow{H}}_{s}$and
$\overrightarrow{\nabla}\Phi$strongly cancel each other. This difficulty
can be avoided when currents are not flowing in the magnetic materials.
Exterior to the volumes where current flows the total field can be
represented using the \"total scalar potential\" $Y$:

> $\overrightarrow{H} = - \overrightarrow{\nabla}\Psi$ , *Total scalar
> potential* (19.8)

which satisfies:

> $\overrightarrow{\nabla} \cdot \mu\overrightarrow{\nabla}\Psi = 0$ .
> (19.9)

By combining the two representations (the total and the reduced scalar
potentials) cancellation difficulties can be completely avoided.
Therefore, the minimal combination for magnetic field simulations
consists of using the reduced potential plus $\overrightarrow{H_{s}}$
only inside volumes where currents flow and using the total scalar
potential everywhere else. On the interface between the total and
reduced potential spaces the two potentials are linked together
(internally in the computer code) by applying the conditions of normal B
and tangential H continuity. This procedure has important consequences
for setting up the simulation model. The reduced potential volume should
completely enclose any coils in such a way that it is not possible to
find a closed path in total potential which encloses a non-zero current.
In other words any closed contour integral of H through any total scalar
region must approach zero, since otherwise the solution would be
multi-valued.

### 19.2 Solution strategy in 2D

In 2D simulations, using either cylindrical or rectangular coordinates,
the solution strategy for Maxwell\'s equations can be dramatically
simplified compared to the general 3D strategy outline above. Since the
primary topic of this report is the calculation of axisymmetric solenoid
systems, we will describe the 2D simulation method in more detail, and
will give a few computation examples in the next chapter. Still, for a
solenoid-based particle accelerator design, 3D simulations would be
necessary in regions where several solenoid beamlines are transversely
merged into a single beamline, or where solenoid lead effects become
critical.

By using a vector potential,
$\overrightarrow{B} = \overrightarrow{\nabla}x\overrightarrow{A}$, we
can rewrite the relevant Maxwell\'s equation in the magnetostatic limit:

> $\overrightarrow{\nabla}x(\frac{1}{\mu}\overrightarrow{\nabla}x\overrightarrow{A}) = \overrightarrow{j}$
> . (19.10)

Imposing gauge condition
$\overrightarrow{\nabla} \cdot \overrightarrow{A} = 0$, and taking into
account the now two-dimensional symmetry of the physical geometry,
equation (19.10) can be rewritten:

> $\frac{\partial}{\partial z}\frac{1}{\mu}\frac{\partial A}{\partial z} + \frac{\partial}{\partial r}\frac{1}{\mu}\frac{1}{r}\frac{\partial(\text{rA})}{\partial r} = - j$
> , (19.11)

Where only the $\theta$ components (A, j) are present. In two dimensions
Maxwell\'s equations reduce to a scalar equation, which can be easily
solved by computational methods. Contrary to three-dimensional
simulations where the model geometry has to be divided into regions of
\"total\" and \"reduced\" scalar potentials to gain a numerical scalar
potential description, the equation to be solved in 2D is automatically
of an effectively scalar nature. No division into different regions has
to be performed in 2D. The magnetic field components are simply

> $B_{r} = - \frac{\partial A}{\partial z}$ , (19.12a)
>
> $B_{z} = \frac{1}{r}\frac{\partial(\text{rA})}{\partial r}$ . (19.12b)

This has an important consequence: Contrary to three dimensions, in a
planar two-dimensional simulation the equipotential lines of the vector
potential (which is a scalar quantity in 2D) are equal to the magnetic
field lines (for axisymmetric simulations a modified scalar potential
$r \cdot A$ has to be used). That fact makes two-dimensional simulations
much more intuitive and accessible from a practical point of view
compared to any three-dimensional simulation.

### 19.3 Example

Figure (2) shows a typical axisymmetric model setup of a solenoid with
an iron enclosure. Figure (3) shows the finite element mesh as generated
by Tosca^©^ -2D. The material properties have to be defined in the code.
For air and the copper conductor the relative permeability is set to 1,
for the return yoke a typical B-H magnetization curve for 1010 carbon
steel, as shown in figure 4, is used.

As shown in figure 5, correct boundary conditions have to be set to
solve the problem numerically. For axisymmetric simulations the center
axis requires a \"tangential field only\" (Dirichlet) boundary
condition. The air volume enclosing the region of interest has to be
chosen as large as possible to minimize the effect of the far-end
boundary conditions. It is always a good practice to check that the
far-field boundaries have been placed far enough by changing the normal
field to tangential field boundaries. The change in the magnetic field
in the region of interest should be minimal, otherwise the outer
boundaries have not been moved far enough out. In the example the
current density in the coil field region has been set to
$4\text{.}7x\text{10}^{6}$ A/m^2^. All that information is stored in a
binary database file which is processed by the analysis module of the
Tosca^©^ code.

![](media/image22.png){width="5.874409448818898in"
height="3.3751968503937007in"}

![](media/image24.png){width="5.874409448818898in"
height="2.999212598425197in"}

![](media/image25.png){width="5.874409448818898in"
height="2.999212598425197in"}

![](media/image26.png){width="5.874409448818898in" height="3.0in"}
Computation of the vector potential proceeds iteratively for each mesh
point. The computation proceeds in the iron region with the
determination of vector potential, magnetic induction, and permeability
at each mesh point. The vector potentials A are first computed at all
points by assuming a constant value of permeability, then the components
of the induction B are determined as partial derivatives of A, and the
permeability values as a function of the absolute values of B are read
from a B-H curve table. Then the whole process is repeated with adjusted
permeability values: cycling is continued until the changes in
permeability are all below a specified value.

Once the field computation has been performed the database is read into
a post-processor for further evaluation. Figure 6 shows the model
geometry with the magnetic flux lines displayed. Figure 7 shows contour
zones of the total magnitude of the magnetic flux density, and figure 8
shows a read out of the magnetic field on axis.

![](media/image27.png){width="5.874409448818898in" height="3.0in"}

![](media/image28.png){width="5.874409448818898in"
height="3.000787401574803in"}

[^1]:
    ### [20. Application of Tosca^©^ to Induction Linac Design]{.underline} {#application-of-tosca-to-induction-linac-design .unnumbered}

    ### 20.1 Introduction {#introduction .unnumbered}

    Recently there has been renewed interest in relatively low energy
    modular induction accelerators to drive close-coupled HIF targets
    (\~75 MV acceleration voltage of Ar^8+^, corresponding to 600 MeV
    kinetic energy). Similar, smaller-scale designs are suggested for a
    possible facility to perform high energy density physics
    experiments. These designs utilize solenoid beam transport because
    of its favorable scaling at the desired ion mass and kinetic energy.

    An example \"modular solenoid driver\" \[20-1\] would consist of 40
    beamlines (20 on each side of the final focus chamber). Each
    beamline would have a high current Ar^8+^ injector and a pre-bunch
    section to compress the beam for injection into a solenoid-based
    modular induction linac. The total linac length would be only 75
    meters. The accelerated beams, which have an approximately 15%
    energy tilt imposed, compress in a neutralized drift section (\~260
    meters long) before final radial compression inside a liquid-vortex
    flibe chamber with a solenoid focusing magnet set.

    Figure 1 shows the mechanical layout of the high-energy end of such
    a modular solenoid linac. The main dimensions of a linac cell are
    given in figure 2. The transport solenoids are superconducting
    Nb~3~Sn magnets producing 12 Tesla peak axial magnetic field.
    Finite-element based electromagnetic simulation tools are utilized
    in several design aspects of the linear accelerator:

    The solenoids are embedded inside induction cores of high
    permeability. To shield the solenoid return flux from the induction
    cores, the superconducting winding structure is enclosed by a thick,
    soft-magnetic steel shell. Finite element codes are reliable tools
    for calculating the effects of the permeable materials and allow
    optimization of the acceleration cell geometry.

    The ion beams exiting the 20 modular solenoid linacs must be matched
    into the focusing system leading to the fusion target.
    Three-dimensional finite element codes can reliably calculate the
    dipole and higher-order field contributions of the solenoid magnetic
    fields at the end of the linac structure, especially if permeable
    materials are present to channel the linac fields' return flux.

    Finite element codes are used to design the superconducting winding
    pack. They allow the exact coil dimensioning using critical current
    load lines.

    If one solenoid along a linac string quenches, the forces on the
    neighboring solenoids will be imbalanced. Finite element codes allow
    one to calculate these forces and to dimension a sufficient support
    system. The electromagnetic simulation models are directly linked
    with structural models to design the solenoid clamping structure.

    ![](media/image29.pct){width="5.99375in"
    height="8.213888888888889in"}

    ![](media/image30.pct){width="5.99375in"
    height="8.213888888888889in"}

    ![](media/image31.pct){width="6.0in" height="4.486805555555556in"}

    **20.2 Example**

    As an example of how electromagnetic modeling can be utilized to
    optimize the modular solenoid driver architecture, we look at the
    effect of adding iron outside the superconducting solenoid winding.
    We will also discus how to specify the superconducting winding
    structure.

    Figure 3 shows a simulation setup for an infinitely long solenoid
    string consisting of a (thin) superconducting wire pack surrounded
    by an iron cylinder. The rest of the model consists of air with a
    relative permeability equal to 1 the location of the induction cells
    is indicated, but we will not include their material properties in
    this example. The beam travels along the vertical axis, and the
    radial coordinate is the horizontal axis.

    Such a simulation should strongly take advantage of the model
    symmetry. First of all the simulation is 2D axisymmetric with the
    beam axis as symmetry axis. Second, because of the infinite nature
    of the problem the top and bottom model boundaries are symmetry
    planes through the centers of the solenoids. The correct numerical
    boundary conditions are shown in figure 4. The background air
    (vacuum) volume extends radially far outside to minimize the
    influence of the far field boundary condition, which was set to
    tangential $B_{r} = 0$.

    ![](media/image32.pct){width="6.0in" height="4.25in"}

    Figure 5 shows in more detail the region around one single solenoid
    and the finite element grid generated by the Tosca^©^ 2D simulation
    code. To enhance the accuracy of the simulation model, regions close
    to the beam axis - including the solenoid winding and the iron
    regions - are modeled using mapped hexagonal finite elements. The
    background air model uses tetragonal elements, which allow gradually
    increasing element sizes. In that way the model can have a very fine
    grid around the areas of interest and a coarser grid at the far
    field boundaries.

    After meshing the problem geometry the material and conductor
    properties have to be defined. Figure 6 shows the permeability curve
    for 1010 steel, a soft magnetic steel used for building
    electromagnets.

    Figure 7 shows a graph of critical current densities for different
    superconducting cable materials \[*2*\]. The blue curve with blue
    squares shows the critical current line for Nb~3~Sn at a temperature
    of 4.2 Kelvin, the temperature of the liquid helium coolant. The
    critical current density is the maximum current density a
    superconducting wire can transport without loosing its
    superconductivity, which is dependent on the applied magnetic field.
    The higher the magnetic field the lower the critical current
    density. Nb~3~Sn has similar critical current densities as NbTi but
    can be operated at much higher applied magnetic fields.

    It is important to realize that the critical current densities are
    shown for single stranded superconducting wire. The actual current
    densities used as input in the simulations to define the conductor
    regions are significantly lower. Figure 8 shows a cross sections of
    a few actual superconducting wires. It shows the superconducting
    wire strands embedded inside a copper matrix whose purpose is
    two-fold.

    ![](media/image33.pct){width="6.0in" height="7.166666666666667in"}

    First it helps to locally stabilize the superconducting wire. On a
    microscopic scale, the superconducting state is fluctuating, and a
    sufficiently sized copper matrix helps to stabilize the
    superconducting condition. Second, the wire has to be able to
    transport all the current in cases where the wire quenches and
    superconductivity is lost entirely. Enough copper must be available
    to transport the current without melting the wire.

    ![](media/image34.pct){width="5.434027777777778in"
    height="3.7243055555555555in"}

    Figure 6: BH curve for 1010 steel, a soft magnetic steel used for
    building electromagnets.

    In designing a magnet, what really matters is the overall
    \"engineering current density\". It is given by the superconducting
    strand critical current density times the superconducting material
    to copper ratio times the packing factor of the actually wound
    superconducting wire. Figure 8 demonstrates the packing of
    commercial round wires in a magnet winding.

    $J_{\text{overall}} = \text{packing factor} \cdot J_{\text{wire}} = \text{packing factor} \cdot \frac{1}{\text{Cu/SC ratio}}J_{c}$

    with

    > $J_{\text{overall}}$ \... \"engineering current density\"
    >
    > $\text{packing factor}$ \... takes account of space occupied by
    > insulation, eventual cooling channels, mechanical reinforcement,
    > epoxy, etc.
    >
    > $\text{Cu/Sc ratio}$ \... ratio of copper to superconducting
    > material\
    > Table 1 determines the maximum Cu/SC ratio by limiting the maximum
    > current density to be carried by the copper matrix to 1500 A/mm@
    > in case of a complete quench.
    >
    > $J_{c}$ \... critical current density for single superconducting
    > strand

    ![](media/image35.pct){width="5.998611111111111in"
    height="7.904166666666667in"}

    Table 1 summarizes the actual "engineering current densities used to
    specify the conductors in the simulations dependent on the applied
    magnetic field. An additional
    ![](media/image36.pct){width="5.998611111111111in"
    height="4.207638888888889in"}

    Figure 8: Cross section of various actual superconducting wires. The
    superconducting filaments, which are actually twisted, are embedded
    inside a copper matrix. Additional space around the wires is needed
    for insulation, structural epoxy, or cooling channels.

    30% safety margin is applied to the engineering current densities to
    account for winding strain and temperature fluctuations. Depending
    on the actual mechanical solenoid design an even larger temperature
    safety margin may be necessary at higher applied magnetic fields.

    Figure 9 displays the field lines of the resulting solenoid beam
    transport channel. Figure 10 shows the axial magnetic field strength
    along the beam center axis. The field ripple due to the fairly large
    magnet spacing, which is necessitated by the wide induction
    acceleration gap (1 MV/m gradient), is clearly visible.

    **B**

    **\[ T \]J~c~**

    **\[ A / mm^2^ \]Cu/SC Ratio ^(\*)^J~wire~**

    **\[ A / mm^2^ \]J~overall~**

    **\[ A / mm^2^ \]with 30% margin ^(\*\*)^**

    **\[ A / mm^2^
    \]**594546.301295906634677665.181257880616764314.291216851596853473.561171820574944462.9611227855501036892.4610667465231130482.0310057044931225001.679386564591320311.358636044231416311.097815473831512890.86693485340

    ^(\*)^ for a maximum current density of 1500 A/mm^2^ in the copper
    matrix

    ^(\*\*)^ margin for winding strain and temperature fluctuations

    Table 1: Engineering current densities (*J~overall~*) for Nb~3~Sn at
    different magnetic fields.

    ![](media/image37.pct){width="6.0in" height="4.209722222222222in"}

    ![](media/image38.pct){width="5.3909722222222225in"
    height="4.158333333333333in"}

    Figure 10: Magnetic field along the beam center axis.

    ![](media/image39.pct){width="5.157638888888889in"
    height="4.065277777777778in"}

    Figure 11: Comparison of \|B\| along the beam center axis and along
    the solenoid inner radius Indicated is the zone of highest magnetic
    field inside the solenoid winding.

    Figure 11 shows$B$ along the beam center axis and along the inner
    coil radius. The field non-uniformity must be included in ion optics
    simulations. Indicated in figure 11 is the location of the maximum
    magnetic field for a long solenoid. The solenoid current density
    must be below the \"engineering\" critical current density for the
    magnetic field encountered in that region.

    Finite element analysis plays an important role in reliably
    predicting magnetic fields inside permeable materials. In the above
    design the iron region around the magnet winding is fully saturated
    (magnetic field \> 2 Tesla), as can be seen in figure 12. Plotted in
    figure 13 is the magnetic field versus radius (at z=0.0, the middle
    of the solenoid, see figure 12) for different iron thicknesses. The
    magnetic field is constant within the solenoid bore, drops within
    the winding pack, and reverses its direction in the iron and the
    surrounding air. Since the iron is fully saturated it cannot
    completely shield the magnetic field return flux from the induction
    cells and the adjacent structures. With a 4cm return iron yoke the
    air magnetic field at the location of the induction cells would be
    around 0.2 Tesla (see figure 14).

    Increasing the iron thickness is the only option to further reduce
    the return field. That effect can be seen in figure 14. A
    significant thickness is required to shield a 12 Tesal solenoid
    field from the induction cells. Currently it is not clear if total
    shielding is

    required to maintain the full flux swing inside the induction cells
    since the induction cell magnetization is perpendicular to the
    solenoid magnetization direction. In addition, the induction cell
    magnetization is a transient effect compared to the DC magnetization
    caused by the solenoid.

    ![](media/image40.pct){width="6.0in" height="4.209722222222222in"}

    ![](media/image41.pct){width="4.666666666666667in"
    height="3.3333333333333335in"}

    Figure 13: Magnetic field (in solenoid center) versus radius for
    different iron thicknesses. See figure 14 for more details.

    ![](media/image42.pct){width="5.382638888888889in"
    height="3.859027777777778in"}

    Figure 14: Solenoid return field in the vicinity of the induction
    cells for different iron return yoke thicknesses.

    **References**

    20-1 Grant Logan: \"What product might a renewal of Heavy Ion Fusion
    development

    offer that competes with methane microbes and hydrogen HTGRs?"
    internal HIF-

    VNL report, May 1, 2006, LBNL-61905.

    20-2 Applied Superconductivity Center, University of
    Wisconsin-Madison,

    http://www.asc.wisc.edu/ or
    [http://magnet.fsu.edu/\~lee/plot/plot.htm\\](http://magnet.fsu.edu/~lee/plot/plot.htm\)

    **[Biot-Savart Formula]{.underline}**

    The Biot-Savart formula in its generalized form is an integration
    along the wire path, and is therefore an attractive tool when the
    path is mathematically characterized and there are no permeable
    materials. It may be used with a CAD program that lays out the wire
    path, which in general is a vector position $\overrightarrow{R}(s)$,
    where $s$ is path length measured along the wire. This is a good way
    to compute the fields of magnet leads and other deviations from
    ideal solenoidal symmetry. The Biot-Savart formula is often
    presented in texts as an experimental result for a straight wire and
    then generalized for an arbitrary wire path. Here we briefly derive
    it from the magnetostatic Maxwell's eqns. (2.1a,b). It is then used
    to obtain the fields of a straight wire segment of finite length and
    of a helix.

    > **21.1 General formula**

    In the absence of magnetic materials the Cartesian components of a
    static magnetic field satisfy

    > $Ñ^{2}\overrightarrow{B} =$
    > $\overrightarrow{Ñ}\overrightarrow{Ñ} \cdot \overrightarrow{B}$
    > $- \overrightarrow{Ñ}x\left( \overrightarrow{Ñ}x\overrightarrow{B} \right)\text{=-}m_{0}\overrightarrow{Ñ}x\overrightarrow{J}$
    > , (21.1)

    so the Green function solution is

    $\overrightarrow{B} = \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}}\frac{\overrightarrow{Ñ^{¢}}x\overrightarrow{J^{¢}}}{\left| \overrightarrow{r} - \overrightarrow{r^{'}} \right|}$
    . (21.2)

    An integration by parts gives

    $\overrightarrow{B} = \frac{m_{0}}{4p}\int_{}^{}{d^{3}r^{'}}\overrightarrow{J^{¢}}x\overrightarrow{Ñ^{¢}}\frac{1}{\left| \overrightarrow{r} - \overrightarrow{r^{'}} \right|}$

    > $= \frac{m_{0}}{4p}\int_{}^{}d^{3}r^{'}\begin{matrix}
    > \overrightarrow{J^{¢}}x\frac{\left( \overrightarrow{r} - \overrightarrow{r^{¢}} \right)}{\left| \overrightarrow{r} - \overrightarrow{r^{¢}} \right|^{3}}
    > \end{matrix}$ . (21.3)

    For a wire (of infinitesimal thickness) with path
    $\overrightarrow{R}(s)$ carrying current I, this may be written
    (Generalized Biot-Savart formula)

    $\overrightarrow{B} = \frac{m_{0}I}{4p}\int_{}^{}{\text{ds}\begin{matrix}
    \hat{e}
    \end{matrix}(s)x\frac{\left( \overrightarrow{r} - \overrightarrow{R}(s) \right)}{\left| \overrightarrow{r} - \overrightarrow{R}(s) \right|^{3}}}$,
    (21.4)

    where $\hat{e}(s)$ is the unit vector along the wire in the
    direction of the current:

    > $\hat{e}(s) = \frac{d\overrightarrow{R}(s)}{\text{ds}}$ . (21.5)

    That $\hat{e}$ is a unit vector is apparent from the definition of
    arc length,

    > $\left( \text{ds} \right)^{2} = \left| d\overrightarrow{R} \right|^{2} = \left( \text{dR}_{x} \right)^{2} + \left( \text{dR}_{y} \right)^{2} + \left( \text{dR}_{z} \right)^{2}$.
    > (21.6)

    **21.2 Straight wire field**

    First we recover the original result of Biot and Savart for an
    infinitely long wire along the z- axis. In this case we have simply

    $\hat{e}(s) = {\hat{e}}_{z},\begin{matrix}
    \overrightarrow{R}(s) = s{\hat{e}}_{z}
    \end{matrix}$. (21.7a,b)

    Eqn. (21.4) becomes

    $\overrightarrow{B} = \frac{\mu_{0}I}{4\pi}{\hat{e}}_{z}x\overrightarrow{r}\int_{- \infty}^{+ \infty}\text{ds}\frac{1}{\sqrt{{r^{2} + (z - s)^{2}}^{3}}}$
    , (21.8)

    where as usual $r^{2} = x^{2} + y^{2}$. The integration is
    elementary in this case:

    > $\overrightarrow{B} = \frac{m_{0}I}{4p}\frac{{\hat{e}}_{z}x\overrightarrow{r}}{r^{2}}\frac{(s - z)}{\sqrt{r^{2} + (s - z)^{2}}}\underset{\text{-¥}}{\overset{\text{+¥}}{|}}\begin{matrix}
    > \begin{matrix}
    >
    > \end{matrix}\frac{m_{0}I}{2p}\frac{{\hat{e}}_{z}x\overrightarrow{r}}{r^{2}} = \frac{m_{0}I}{2pr}{\hat{e}}_{q}
    > \end{matrix}$ . (21.9)

    For a semi-infinite wire running from $- \infty$ to 0 we have

    $\overrightarrow{B} = \frac{m_{0}I}{4p}\frac{{\hat{e}}_{q}}{r}\left( 1 - \frac{z}{\sqrt{r^{2} + z^{2}}} \right)$,
    (21.10)

    reminiscent of the on-axis field of the semi-infinite thin-layer
    solenoid treated in section 4. However in this case there must be an
    additional current path (unspecified here) leading away from the
    wire end at $z = 0$ to avoid a non-zero divergence of
    $\overrightarrow{J}$.

    A useful generalization of the above example is a straight wire
    segment between arbitrary points $\overrightarrow{r_{1}}$ and
    $\overrightarrow{r_{2}}$; we define

    > $\overrightarrow{C} = \frac{\overrightarrow{r_{1}} + \overrightarrow{r_{2}}}{2}$=
    > wire center, (21.11a)

    $L = \left| \overrightarrow{r_{2}} - \overrightarrow{r_{1}} \right|$=
    wire length, (21.11b)

    $\hat{e} = \frac{\overrightarrow{r_{2}} - \overrightarrow{r_{1}}}{L}$=
    wire direction. (21.11c)

    Then the wire path is

    > $\overrightarrow{R}(s) = \overrightarrow{C} + s\hat{e}$ , (21.11d)

    with $- L/2 < s < L/2$. Eqn. (21.4) becomes

    > $\overrightarrow{B} = \frac{m_{0}I}{4p}\int_{- L/2}^{L/2}\text{ds}\begin{matrix}
    > \hat{e}x\frac{\left\lbrack \overrightarrow{r} - \left( \overrightarrow{C} + s\hat{e} \right) \right\rbrack}{\left| \overrightarrow{r} - \left( \overrightarrow{C} + s\hat{e} \right) \right|^{3}}
    > \end{matrix}$ . $$ (21.12)

    The integration is again elementary, but it helps to employ the
    definitions

    > $\alpha = \hat{e} \cdot \left( \overrightarrow{r} - \overrightarrow{C} \right)\begin{matrix}
    >
    > \end{matrix}$ , (21.13)

    with
    $\left| \overrightarrow{r} - \overrightarrow{C}|^{2} = \alpha^{2} + |\overrightarrow{\beta} \right|^{2}$.
    Making the change of variable $t = s - \alpha$, we get from

    eqn. (21.12)

    $\begin{array}{r}
    \overrightarrow{B} = \frac{m_{0}I}{4p}\frac{\overrightarrow{\beta}}{\left| \overrightarrow{\beta} \right|^{2}}\frac{t}{\sqrt{\left| \overrightarrow{\beta} \right|^{2} + t^{2}}}\underset{- L/2 - \alpha}{\overset{L/2 - \alpha}{|}} \\
    \frac{m_{0}I}{4p}\frac{\overrightarrow{\beta}}{\left| \overrightarrow{\beta} \right|^{2}}\left\lbrack \frac{(L/2 - \alpha)}{\sqrt{\left| \overrightarrow{\beta} \right|^{2} + (L/2 - \alpha)^{2}}} + \frac{L/2 + \alpha}{\sqrt{\left| \overrightarrow{\beta} \right|^{2} + (L/2 + \alpha)^{2}}} \right\rbrack\begin{matrix}

    \end{matrix}\text{.}
    \end{array}$ (21.14)

    This formula may be used for the field from a straight magnet lead.
    Also, the field of a complete circuit with curved wired sections can
    be computed as the sum of fields from short straight sections, each
    with specified endpoints and defined quantities
    $\left( \overrightarrow{C},L,\hat{e},\alpha,\overrightarrow{\beta} \right)$.

    **21.3 Helix field**

    The path representation of a helix is more complicated; let $\rho$
    be its radius (wire distance from $z$-axis), $p$is the helix pitch
    angle restricted to $(0 < p < p2)$, and $t$ is a variable
    proportional to $s$ such that $0 \leq \tau \leq 2\pi$gives one loop.
    Then we may write

    > $R_{x} = \rho\text{cos}(t),\begin{matrix}
    >
    > \end{matrix}$. (21.15)

    We have
    $\left( \text{ds} \right)^{2} = \left| d\overrightarrow{R} \right|^{2} = \rho^{2}\left\lbrack 1 + \text{tan}^{2}(p) \right\rbrack(dt)^{2}$,
    so we set

    $t = \frac{s}{\rho\sqrt{1 + \text{tan}^{2}(p)}}\begin{matrix}

    \end{matrix}$ , (21.16)

    and $R_{z} = s\begin{matrix}
    \text{sin}(P)
    \end{matrix}$, as expected from the definition of $P$.

    The helix period length is $2\text{pρ}\text{tan}(p) = 1/N$, where
    $N$ is the number of turns per meter, so the mean $\theta$-component
    of surface current density is

    > $\left\langle S \right\rangle = \text{NI}$, (21.17)

    while the total z-current equals I.

    For an infinitely long helix, the on-axis $z$ component of the field
    is by symmetry

    > $B_{0}(z) = m_{0}\left\langle S \right\rangle = m_{0}\text{NI}$.
    > (21.18)

    For $r\text{>>}\rho$ we also have

    > $B_{\theta} \rightarrow \frac{\mu_{0}I}{2\pi r}$. (21.19)

    Close to the wire, which has infinitesimal thickness, the field
    diverges strongly. The on-axis transverse field does not vanish, but
    has a constant absolute value with direction rotating
    90­­${}^{\circ}$out of phase from the wire path. It is an interesting
    exercise to visualize the field lines.

    To proceed with the helix field calculation, first note that

    $\hat{e}(s) = \frac{d\overrightarrow{R}}{\text{ds}} = \frac{dt}{\text{ds}}\frac{d\overrightarrow{R}}{dt} = \text{cos}(p)\left\lbrack - {\hat{e}}_{x}\text{sin}(t) + {\hat{e}}_{y}\text{cos}(t) + {\hat{e}}_{z}\text{tan}(p) \right\rbrack$
    ,

    $\begin{array}{r}
    \hat{e}(s)x\left( \overrightarrow{r} - \overrightarrow{R}(s) \right) = \text{cos}(p)\left\{ {\hat{e}}_{x}\left\lbrack \left( z - R_{z} \right)\text{cos}(t) - \left( y - R_{y} \right)\text{tan}(p) \right\rbrack \right.\  \\
    \left. \  + {\hat{e}}_{y}\left\lbrack \left( x - R_{x} \right)\text{tan}(p) + \left( z - R_{z} \right)\text{sin}(t) \right\rbrack - {\hat{e}}_{z}\left\lbrack \left( y - R_{y} \right)\text{sin}(t) + \left( x - R_{x} \right)\text{cos}(t) \right\rbrack \right\} \\

    \end{array}$. (21.20)

    Each Cartesian field component may be projected from eqn. (21.4);
    for example

    > $B_{z} = \overrightarrow{B} \cdot {\hat{e}}_{z}\text{=-}\frac{m_{0}I}{4p}\int_{}^{}\text{ds}\frac{\text{cos}(p)\left\lbrack \left( y - R_{y} \right)\text{sin}(t) + \left( x - R_{x} \right)\text{cos}(t) \right\rbrack}{\sqrt{{\left( x - R_{x} \right)^{2} + \left( y - R_{y} \right)^{2} + \left( z - R_{z} \right)^{2}}^{3}}}$
    > . (21.21)

    This may be evaluated by a numerical integration package, using the
    given definitions of

    $\overrightarrow{R}(s)$. It helps to use $\tau = s\begin{matrix}
    \text{cos}\frac{(p)}{\rho}
    \end{matrix}$ as the integration variable. A helix might have for
    example m periods with $0 \leq \tau \leq 2\pi m$ and longitudinal
    dimension running from z=0 to z=$2\pi m\rho\text{tan}(p)$. The
    components $B_{x}$and $B_{y}$are found in similar fashion.

    For an infinite helix the on-axis longitudinal field is

    $B_{\text{zo}} = \frac{m_{0}I}{4p}\text{cos}(p)\int_{}^{}\text{ds}\frac{\left\lbrack R_{y}\text{sin}(t) + R_{x}\text{cos}(t) \right\rbrack}{\sqrt{{R_{x}^{2} + R_{y}^{2} + \left( z - R_{z} \right)^{2}}^{3}}}$
    (21.22)

    > $= \frac{m_{0}I}{4p}\text{cos}(p)\int_{}^{}\frac{\text{du}}{\text{sin}(p)}\frac{\rho}{\sqrt{{\rho^{2} + u^{2}}^{3}}}$
    > , (21.23)

    where $u = R_{z} - z = s\begin{matrix}
    \begin{matrix}
    \text{sin}(p) - z
    \end{matrix}
    \end{matrix}$. The integration is elementary:

    $B_{z_{o}} = \frac{m_{0}I\rho}{4p\begin{matrix}
    \text{tan}(p)
    \end{matrix}}\frac{1}{\rho^{2}}\frac{u}{\sqrt{\rho^{2} + u^{2}}}\begin{matrix}

    \end{matrix}$ , (21.24)

    as advertised in eqn. (21.18).

    The transverse fields of a finite or infinite helix are

    $B_{x} = \overrightarrow{B \cdot {\hat{e}}_{x} = \frac{m_{0}I}{4p}\int_{}^{}{\text{ds}\begin{matrix}

    \end{matrix}\frac{\text{cos}(p)\left\lbrack \left( z - R_{z} \right)\text{cos}(t) - \left( y - R_{y} \right)\text{tan}(p) \right\rbrack}{\sqrt{{\left( x - R_{x} \right)^{2} + \left( y - R_{y} \right)^{2} + \left( z - R_{z} \right)^{2}}^{3}}}}}$,
    (21.25)

    $B_{y} = \overrightarrow{B \cdot {\hat{e}}_{y} = \frac{m_{0}I}{4p}\int_{}^{}{\text{ds}\begin{matrix}

    \end{matrix}\frac{\text{cos}(p)\left\lbrack \left( x - R_{x} \right)\text{tan}(p) + \left( z - R_{z} \right)\text{sin}(t) \right\rbrack}{\sqrt{{\left( x - R_{x} \right)^{2} + \left( y - R_{y} \right)^{2} + \left( z - R_{z} \right)^{2}}^{3}}}}}$.
    (21.26)

    These components may also be evaluated numerically using $\tau$ as
    the variable of integration.

    The on-axis transverse fields for an infinite helix are

    > $B_{\text{xo}} = \frac{m_{0}I\begin{matrix}
    > \text{cos}(p)
    > \end{matrix}}{4p}\int_{\text{-¥}}^{\text{+¥}}\text{ds}\begin{matrix}
    > \frac{\left\lbrack \rho\text{sin}(t)\text{tan}(p) + \left( z - \text{tρ}\text{tan}(p) \right)\begin{matrix}
    > \text{cos}(t)
    > \end{matrix} \right\rbrack}{\sqrt{{\rho^{2} + \left( z - \text{tρ}\text{tan}(p) \right)^{2}}^{3}}}
    > \end{matrix}$, (21.27)
    >
    > $B_{\text{yo}} = \frac{m_{0}I\begin{matrix}
    > \text{cos}(p)
    > \end{matrix}}{4p}\int_{\text{-¥}}^{\text{+¥}}\text{ds}\begin{matrix}
    > \frac{\left\lbrack \left( z - \text{tρ}\text{tan}(p) \right)\text{sin}(t) - \rho\text{cos}(t)\begin{matrix}
    > \text{tan}(p)
    > \end{matrix} \right\rbrack}{\sqrt{{\rho^{2} + \left( z - \text{tρ}\text{tan}(p) \right)^{2}}^{3}}}
    > \end{matrix}$. (21.28)

    Making the change of variable $u = R_{z} - z = s\begin{matrix}
    \text{sin}(p) - z
    \end{matrix}$ we get

    $B_{\text{xo}} = \frac{m_{0}I}{4p\begin{matrix}
    \text{tan}(p)
    \end{matrix}}\int_{\text{-¥}}^{\text{+¥}}\text{du}$, (21.29)

    $B_{\text{yo}} = \frac{- m_{0}I}{4p\begin{matrix}
    \text{tan}(p)
    \end{matrix}}\int_{\text{-¥}}^{\text{+¥}}\text{du}\frac{\left\lbrack u\text{sin}\left( \frac{u + z}{\rho\text{tan}(p)} \right) + \rho\text{tan}(p)\text{cos}\left( \frac{u + z}{\rho\text{tan}(p)} \right) \right\rbrack}{\sqrt{{\rho^{2} + u^{2}}^{3}}}$.
    (21.30)

    It is seen that $B_{x_{0}}$ vanishes at $z = 0$,
    $\pm \text{πρ}\begin{matrix}
    \text{tan}(p)
    \end{matrix}$, $\pm 2\text{πρ}\begin{matrix}
    \text{tan}(p)
    \end{matrix}$, ... , consistent with the helical symmetry. At the
    coordinate origin we have, setting $u = \text{ρυ}$,

    $B_{y_{0}}(0)\text{=-}\frac{m_{0}I}{2\text{pρ}\begin{matrix}
    \text{tan}(p)
    \end{matrix}}\int_{}^{}{du}\left\lbrack \frac{u\text{sin}\left( \frac{u}{\text{tan}(p)} \right)}{\sqrt{{1 + u^{2}}^{3}}} + \frac{\text{tan}(p)\text{cos}\left( \frac{u}{\text{tan}(p)} \right)}{\sqrt{{1 + u^{2}}^{3}}} \right\rbrack$
    (21.31)

    > $= - \frac{\mu_{0}I}{2\text{πρ}\begin{matrix}
    > \text{tan}(p)
    > \end{matrix}}\left\lbrack \frac{Κ_{0}\left( 1/\text{tan}(p) \right)}{\text{tan}(p)} + Κ_{1}\left( 1/\text{tan}(p) \right) \right\rbrack$
    > , (21.32)

    where $Κ_{0}$ and $Κ_{1}$are the modified Bessel functions. Setting
    $1/\text{tan}(p) = 2\text{pρ}N$ we have

    > $B_{y_{0}}(0)\text{=-}m_{0}\text{IN}\left\lbrack \left( 2\text{pρ}N \right)K_{0}\left( 2\text{pρ}N \right) + K_{1}\left( 2\text{pρ}N \right) \right\rbrack$
    > . (21.33)

    Invoking periodicity yields

    > $\overrightarrow{B_{\hat{{}_{{}_{0}}} = B_{y_{0}}(0)}}$ . (21.34)

    Generally the parameter 2$\text{pρ}N$ is much larger that unity, and
    the Bessel functions may be expanded; the leading terms gives

    $B_{y_{0}}(0)\text{»-}\left( m_{0}\text{IN} \right)\sqrt{\frac{p}{2}}\left( \sqrt{2\text{pρ}N} + \frac{7}{8\sqrt{2\text{pρ}N}} \right)\begin{matrix}
    e^{- 2\text{pρ}N}
    \end{matrix}$ , (21.35)

    which is usually a very small on-axis effect. For an concrete
    example set $$, then $2\text{pρ}N = \text{12}\text{.}\text{566}$,
    and

    $\frac{B_{y_{0}}(0)}{- m_{0}\text{IN}}»\sqrt{\frac{p}{2}}\left( \sqrt{\text{12}\text{.}\text{566}} + \frac{7}{8\sqrt{\text{12}\text{.}\text{566}}} \right)e^{- \text{12}\text{.}\text{566}}@\text{.}\text{000017}$.
    (21.36)

    **22. [Single Particle Equations of Motion in an Arbitrary
    Electromagnetic Field]{.underline}**

    Although the topic of this report is the calculation of fields, we
    also include this brief section on the particle equations of motion.
    The goal is to have a convenient reference where these equations are
    written using longitudinal distance $z$ as the independent variable
    instead of time $t$. This form displays the geometric aberrations,
    which are mixed with field aberrations and are hidden in the more
    compact time-dependent formulation.

    **Time as the independent variable**

    With time as the independent variable, an ion's velocity
    $\overrightarrow{v}$ satisfies the Lorentz equation

    $\frac{d\left( g\overrightarrow{v} \right)}{\text{dt}} = \frac{\text{qe}}{M}\left( \overrightarrow{E} + \overrightarrow{v}x\overrightarrow{B} \right)$,
    (22.1)

    with relativistic factor

    $g = \left( 1 - v^{2}/c^{2} \right)^{- 12}$
    $= \left( 1 + g^{2}v^{2}/c^{2} \right)^{12}$ . (22.2)

    Here $q$ is the ion's charge state, $M$ is its mass,
    $v = \left| \overrightarrow{v} \right|$, and we use the natural
    constants

    $e = 1\text{.}\text{6022}\text{.}\text{.}\text{.}x\text{10}^{- \text{19}}C,\begin{matrix}
    c = 2\text{.}\text{9979}\text{.}\text{.}\text{.}x\text{10}^{8}m/s
    \end{matrix}$ .

    The electric and magnetic fields,
    $\overrightarrow{E}\left( \overrightarrow{r},t \right)$ and
    $\overrightarrow{B}\left( \overrightarrow{r},t \right)$, are assumed
    to be known or are computable.

    **22.2 Longitudinal distance as the independent variable**

    For time-dependent, multi-particle simulations it is convenient to
    use eqn. (22.1) as it stands, or with $g = 1$ for non-relativistic
    cases. However there are situations where it is preferable to use
    the longitudinal position $z$ as the independent variable. One
    example is the simulation of an ion source operating in steady
    state. In principle we may always use $z$ instead of $t$ if there is
    no turning back of ion orbits. In addition to the particle variables
    $\left\{ x(z),y(z),v_{x}(z),v_{y}(z),v_{z}(z),v(z),g(z) \right\}$ we
    also calculate $t(z)$, the time when an ion is at $z$. The fields
    the ion experiences are then

    $\overrightarrow{B} = \overrightarrow{B}\left( \overrightarrow{r}(z),t(z) \right),\begin{matrix}
     & 
    \end{matrix}\overrightarrow{E} = \overrightarrow{E}\left( \overrightarrow{r}(z),t(z) \right)\text{.}$
    (22.3a,b)

    The time derivative in eqn. (22.1) is converted to a $z$-derivative
    by the relation

    $\frac{d}{\text{dt}} = v_{z}\frac{d}{\text{dz}}$ . (22.4)

    Derivatives with respect to $z$ are denoted here by a "prime", i.e.
    $x^{'}(z) = \text{dx}/\text{dz}$, etc. Then we get

    $\begin{array}{r}
    v_{x} = \frac{\text{dx}}{\text{dt}} = v_{z}\frac{\text{dx}}{\text{dz}} = v_{z}x^{'}, \\
    v_{y} = \frac{\text{dy}}{\text{dt}} = v_{z}\frac{\text{dy}}{\text{dz}} = v_{z}y^{'}\text{.}
    \end{array}$ (22.5a,b)

    The absolute velocity is

    $v(z) = \sqrt{v_{x}^{2} + v_{y}^{2} + v_{z}^{2}},$ (22.6)

    so substituting from eqn. (22.5) (assuming $v_{z} > 0$):

    $v_{z} = \frac{v(z)}{\sqrt{1 + x^{} + y^{}}}$ . (22.7)

    The time derivative may therefore be written

    $\frac{d}{\text{dt}} = \frac{v(z)}{\sqrt{1 + x^{} + y^{}}}\frac{d}{\text{dz}}$
    . (22.8)

    The factor$\sqrt{1 + x^{} + y^{}}$ in eqns. (22.7, 22.8), which we
    will denote by $\sqrt{\text{.}\text{.}\text{.}}$, is the cause of
    the geometric aberrations for transport close to the $z$ axis.

    Eliminating time, eqn. (20.1) yields

    $\begin{array}{r}
    \frac{v}{\sqrt{\text{.}\text{.}\text{.}}}\frac{d}{\text{dz}}\frac{gv}{\sqrt{\text{.}\text{.}\text{.}}}\frac{\text{dx}}{\text{dz}} = \frac{\text{qe}}{M}\left( E_{x} + \frac{vy^{'}}{\sqrt{\text{.}\text{.}\text{.}}}B_{z} - \frac{v}{\sqrt{\text{.}\text{.}\text{.}}}B_{y} \right)\begin{matrix}
    ,
    \end{matrix}\frac{v}{\sqrt{\text{.}\text{.}\text{.}}}\frac{d}{\text{dz}}\frac{gv}{\sqrt{\text{.}\text{.}\text{.}}}\frac{\text{dy}}{\text{dz}} = \frac{\text{qe}}{M}\left( E_{y} + \frac{v}{\sqrt{\text{.}\text{.}\text{.}}}B_{x} - \frac{vx^{'}}{\sqrt{\text{.}\text{.}\text{.}}}B_{z} \right)\begin{matrix}

    \end{matrix}\text{.}
    \end{array}$ (22.9a,b)

    We also need an equation for $g$ (or equivalently $v$ or $gv$). From
    eqns. (22.1) and (22.2) we derive

    > $\frac{d}{\text{dt}}g = \frac{\text{qe}}{\text{Mc}^{2}}\overrightarrow{v} \cdot \overrightarrow{E,}$
    > (22.10)

    which becomes

    $\frac{dg}{\text{dz}} = \frac{\text{qe}}{\text{Mc}^{2}}\left( x^{'}E_{x} + y^{'}E_{y} + E_{z} \right)\begin{matrix}
    ,
    \end{matrix}$ (22.11a)

    or equivalently

    $\begin{array}{r}
    \frac{\text{dv}}{\text{dz}} = \frac{\text{qe}}{M}\frac{1}{g^{3}v}\left( x^{'}E_{x} + y^{'}E_{y} + E_{z} \right)\begin{matrix}
    ,
    \end{matrix}\frac{dgv}{\text{dz}} = \frac{\text{qe}}{M}\frac{1}{v}\left( x^{'}E_{x} + y^{'}E_{y} + E_{z} \right)\begin{matrix}

    \end{matrix}\text{.}
    \end{array}$ (22.11b,c)

    The equation for time is simply

    $\frac{\text{dt}}{\text{dz}} = \frac{1}{v_{z}} = \frac{\sqrt{\text{.}\text{.}\text{.}}}{v}$
    . (22.12)

    Eqns. (22.9a.b) are of an awkward form since, due to the
    differentiation of $\sqrt{\text{.}\text{.}\text{.}}$, the second
    derivatives $x^{''}$ and $y^{''}$ appear together. A simplification
    is realized when these equations are combined linearly to place
    these second derivatives in separate equations. It also helps to
    eliminate $d(gv)/\text{dz}$ using eqn. (22.11c). These steps finally
    yield after some algebra:

    $x^{''} = \frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left\{ \sqrt{{\text{.}\text{.}\text{.}}^{2}}\left( E_{x} - x^{'}E_{z} \right) + v\sqrt{\text{.}\text{.}\text{.}}\left\lbrack x^{'}y^{'}B_{x} - \left( 1 + x^{} \right)B_{y} + y^{'}B_{z} \right\rbrack \right\}$,
    (22.13a)

    $y^{''} = \frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left\{ \sqrt{{\text{.}\text{.}\text{.}}^{2}}\left( E_{y} - y^{'}E_{z} \right) + v\sqrt{\text{.}\text{.}\text{.}}\left\lbrack - x^{'}y^{'}B_{y} + \left( 1 + y^{} \right)B_{x} - x^{'}B_{z} \right\rbrack \right\}$.
    (22.13b)

    Eqns. (22.11), (22.12), and (22.13) are valid as long as
    $\text{dz}/\text{dt} > 0$; they contain all aberrations and may be
    used with tilted system axes and non-solenoidal fields.

    **Paraxial approximation**

    > Often a set of equations is used, in which terms such as $x^{}$
    > are known to be small

    and are dropped. This is the "paraxial approximation", with formal
    ordering $d/\text{dz}\text{\textasciitilde}e\text{<<}1$ and
    $(x,y)\text{\textasciitilde}1$:

    $\frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left( E_{z},\text{vB}_{z} \right)\text{\textasciitilde}e$
    , (22.14)

    > $\frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left( E_{x},E_{y},\text{vB}_{x},\text{vB}_{y} \right)\text{\textasciitilde}e^{2}$
    > . (22.15)

    Dropping terms which are explicitly small of order $e^{2}$, gives

    $c^{2}\frac{dg}{\text{dz}} = g^{3}v\frac{\text{dv}}{\text{dz}} = v\frac{dgv}{\text{dz}} = \frac{\text{qe}}{M}E_{z}$
    , (22.16)

    $\frac{\text{dt}}{\text{dz}} = \frac{1}{v}$ , (22.17)

    $x^{''} = \frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left\lbrack \left( E_{x} - x^{'}E_{z} \right) + v\left( y^{'}B_{z} - B_{y} \right) \right\rbrack$
    , (22.18a)

    $y^{''} = \frac{\text{qe}}{M}\frac{1}{{gv}^{2}}\left\lbrack \left( E_{y} - y^{'}E_{z} \right) - v\left( x^{'}B_{z} - B_{x} \right) \right\rbrack$
    . (22.18b)

    Field aberrations are still present in eqns. (22.16, 22-17, 22.18)
    and these include terms, which are small in powers of $e$, so the
    paraxial approximation by itself is not a consistent ordering.