Magnetic Multipole Fields

John Jowett

Introduction

This notebook is largely of a tutorial nature and shows how some basic accelerator physics related to magnetic multipoles can be implemented in Mathematica.  It includes generation of algebraic expressions and plots for the potentials, a method for identification of the multipole content of an arbitrary harmonic potential function and manipulations of the Hamiltonian for betatron motion to derive the tune-dependences on close orbit, momentum deviation (chromaticity) and action variables.  This is done for arbitrary combinations of multipoles including feed-down effects of closed orbit.

History

Originally to make pictures for Summer Students in 1996.

Added multipole-coefficient finder after query from Jean-Pierre Koutchouk, 10/2/2000.

Adding calculations of detuning with amplitude, etc. following DAΦNE MAP3, 6/2/2001.

Future: probably worth encapsulating the functions from the notebook in a package.  Implement the canonical transformation from first principles.  The detuning function is not as elegant as it should be though it does the job.

Updated, replacing Notation[ ] cells 6/8/2012.

Setup and definitions

Load packages

Needs["Accelerator`ConstantsUnits`"]

Needs["Graphics`ColorListPlot`"]

Some notations

(Mathematica'a Notation package is loaded by our ConstantsUnits package.)

The Notation function seems to have changed since I wrote the ConstantsUnits palette.   It must now be entered using the notation templates.

Once you have the palette on your screen, it’s a good idea to install it from the Palettes menu so that you can easily find it again.

Notation[(NotationTemplateTag[(Subscript[J, x])]) ⟺ (NotationTemplateTag[Jx])];
IntroduceSymbol[Jx,"is the horizontal betatron action.",Meter];

Notation[(NotationTemplateTag[(Subscript[J, y])]) ⟺ (NotationTemplateTag[Jy])];
IntroduceSymbol[Jy,"is the vertical betatron action.",Meter];

Notation[(NotationTemplateTag[(Subscript[ϕ, x])]) ⟺ (NotationTemplateTag[phix])];
IntroduceSymbol[phix,"is the horizontal betatron phase.",1];

Notation[(NotationTemplateTag[(Subscript[ϕ, y])]) ⟺ (NotationTemplateTag[phiy])];
IntroduceSymbol[phiy,"is the horizontal betatron phase.",1];

Notation[(Subscript[Q, x]) ⟺ Qx];
IntroduceSymbol[Qx,"is the horizontal tune.",1];

Notation[(Subscript[Q, y]) ⟺ Qy];
IntroduceSymbol[Qy,"is the horizontal tune.",1];

Notation[(NotationTemplateTag[(Subscript[β, x])]) ⟺ (NotationTemplateTag[betx])];
IntroduceSymbol[betx,"is the horizontal optical functio.",Meter];

Notation[(NotationTemplateTag[(Subscript[β, y])]) ⟺ (NotationTemplateTag[bety])];
IntroduceSymbol[bety,"is the vertical optical functio.",Meter];

Notation[(NotationTemplateTag[(Subscript[D, x])]) ⟺ (NotationTemplateTag[Dx])]
IntroduceSymbol[Dx,"is the horizontal dispersion function.",Meter];

Notation[(NotationTemplateTag[(Subscript[D, y])]) ⟺ (NotationTemplateTag[Dy])]
IntroduceSymbol[Dy,"is the vertical dispersion function.",Meter];

Functions for multipole fields

normalPotential[n_Integer]:=ComplexExpand[((Re[(((x+I y))^n)])/n)]

skewPotential[n_Integer]:=ComplexExpand[Im[((((x+I y))^n)/n)]]

P[x_,y_]=normalPotential[3]

F=-∇P ; (*   *)

D[P[x,y],x]
D[P[x,y],y]

Formulas for Some Multipoles

TableForm[Table[normalPotential[n],{n,1,10}]]

TableForm[Table[skewPotential[n],{n,1,10}]]

Potential  and Field Pictures for Multipole Magnets

Colouring functions

The use of a quadratic colour function sometimes helps to emphasise the middle range of values in contour plots (otherwise they tend to have large areas of green).  Choose one of the following

SetOptions[ContourPlot,ColorFunction->(Hue[.7#1]&)];

SetOptions[ContourPlot,ColorFunction->(Hue[.7 ((#1)^2)]&)];

SetOptions[ContourPlot,ColorFunction->(Hue[.7(1-((#)^2))]&)];

SetOptions[ContourPlot,ColorFunction->(Hue[.7(Cos[π #/2])]&)];

Dipole

ContourPlot[normalPotential[1],{x,-1,1},{y,-1,1}]

ContourPlot[skewPotential[1],{x,-1,1},{y,-1,1}]

VectorPlot[Evaluate[(∂_({{x,y}}))normalPotential[1]],{x,-2,2},{y,-2,2}]

Quadrupole

ContourPlot[normalPotential[2],{x,-1,1},{y,-1,1}]

ContourPlot[skewPotential[2],{x,-1,1},{y,-1,1}]

VectorPlot[Evaluate[(∂_({{x,y}}))normalPotential[2]],{x,-2,2},{y,-2,2}]

Sextupole

ContourPlot[normalPotential[3],{x,-1,1},{y,-1,1}]

VectorPlot[Evaluate[(∂_({{x,y}}))normalPotential[3]],{x,-2,2},{y,-2,2}]

Octupole

Timing[ContourPlot[normalPotential[4],{x,-1,1},{y,-1,1},PlotPoints->40]]

Combined function example

SetOptions[ContourPlot,Contours->30,ContourLines->False,PlotPoints->50];

Increase the number of contours.  This will take longer!

Timing[ContourPlot[normalPotential[2]+2skewPotential[3],{x,-1,1},{y,-1,1}]]

Timing[ContourPlot[normalPotential[1]+2normalPotential[2],{x,-1,1},{y,-1,1}]]

Although we did not bother in the earlier examples, we can speed up these plots dramatically using Evaluate so that the form of the potential function will only be evaluated once.  It makes a difference now that we are using many more contours and plot points.

Timing[ContourPlot[Evaluate[normalPotential[2]+.4 skewPotential[6]],{x,-1,1},{y,-1,1}]]

Now we can really increase the number of contours and plot points to make a smooth plot.

Timing[ContourPlot[Evaluate[normalPotential[2]+.4 skewPotential[6]],{x,-1,1},{y,-1,1},PlotPoints->50,Contours->100]]

Timing[ContourPlot[Evaluate[normalPotential[3]+4normalPotential[10]+.4 skewPotential[6]],{x,-1,1},{y,-1,1},PlotPoints->50,Contours->100]]

Identifying  the Multipoles in an Arbitrary Magnetic Field

An example to work on

Generate a random multipole field to develop the method on

ϕ=(∑_(k=2)^7)(RandomReal[{-10,10}] normalPotential[k]+RandomReal[{-10,10}] skewPotential[k])

We can still see the answer of course!  Expand out all the brackets to merge the monomials and forget.

ϕ=Expand[ϕ]

We can plot this of course

ContourPlot[ϕ,{x,-1,1},{y,-1,1},PlotPoints->60]

Development of the method

Using the Apply function to replace the head of the expression with List  separates the monomials.

List@@ϕ

Try to isolate the order 4 terms.  The trick of multiplying by (x^2)(y^2) makes the pattern-matching much simpler (we don't need special instances for 0 and 1 powers).

ord4=Plus@@((Cases[(x^2) (y^2) List@@ϕ,a_Real (x^(m_)) (y^(n_))/;m+n==4+4] )/((x^2)(y^2)))

Here is a way to extract the (Subscript[a, 4]) and (Subscript[b, 4]).

Solve[{
(ord4==a4 normalPotential[4]+b4 skewPotential[4]/.{x->1,y->2}),(ord4==a4 normalPotential[4]+b4 skewPotential[4]/.{x->2,y->1})
}
,{a4,b4}]

Laplacian

Now we know what to do.  But first pause to define a 2D Laplacian since we must check the function is harmonic

laplacian2D[ϕ_,{x_,y_}]:=Chop[(∂_({x,2}))ϕ+(∂_({y,2}))ϕ]

laplacian2D[ϕ,{x,y}]

Encapsulate the method as a function for a given order

Assemble the above method into a function that will extract the multipole coefficients for a given order.

multipoles[ϕ_,{x_,y_},n_Integer]:=If[PolynomialQ[ϕ,{x,y}]&&laplacian2D[ϕ,{x,y}]==0,
Module[{ordn,an,bn},
ordn=Plus@@((Cases[(x^2) (y^2) List@@ϕ,a_ (x^(m_)) (y^(k_))/;m+k==n+4] )/((x^2)(y^2)));
{an,bn}/.First[Solve[{
(ordn==bn normalPotential[n]+an skewPotential[n]/.{x->1,y->2}),(ordn==bn normalPotential[n]+an skewPotential[n]/.{x->2,y->1})
}
,{an,bn}]]
],Print["multipoles: not a harmonic function"];Null]

Test the function

multipoles[ϕ,{x,y},3]

multipoles[ϕ,{x,y},4]

multipoles[ϕ,{x,y},5]

multipoles[ϕ,{x,y},6]

multipoles[(1/x)+(x^2)-(y^2),{x,y},2]

Should really use the Message[] thing.

Overload the function to treat all orders at once

Overload the function to get all the multipoles.  First find the maximum order we need.

Max[Dimensions[CoefficientList[ϕ,{x,y}]]]

multipoles[ϕ_,{x_,y_}]:=Table[Flatten[{p,multipoles[ϕ,{x,y},p]}],{p,1,Max[Dimensions[CoefficientList[ϕ,{x,y}]]]}]

This version returns a list of triples of the form {n,(Subscript[a, n]),(Subscript[b, n])}.

Tests of the extended function

multipoles[ϕ,{x,y}]

TableForm[multipoles[ϕ,{x,y}],TableHeadings->{{},{"n","(Subscript[a, n])","(Subscript[b, n])"}}]

Check that it also works with symbolic coefficients

ψ=Expand[(f/((u-1))) normalPotential[3]+(((g-u))^3) skewPotential[4]]

multipoles[ψ,{x,y}]

Performance analysis

Actually this is not very efficient since the testing is done too often but never mind.  Have a quick look at how the timing goes as we go to potentials of higher and higher order.  Here we plot the time taken in Second as a function of the maximum multipole order in a randomly generated potential function.

ListPlot[Table[ϕ=Expand[(∑_(k=2)^p)(RandomReal[{-10,10}] normalPotential[k]+RandomReal[{-10,10}] skewPotential[k])];{p,First[Timing[multipoles[ϕ,{x,y}]]]},{p,2,50}],Prolog->PointSize[0.2`]]

So the performance is quite good enough unless we really have to do a lot of these.  So I won't bother to optimise.

Tune-shifts with Closed-orbits and Amplitude, Chromaticity

Vector Potential

Use the same vector notation for multipoles as in MAD9.

As[KN_List,KS_List]:=Plus@@MapIndexed[#1 normalPotential[First[#2]]&,KN]+Plus@@MapIndexed[#1 skewPotential[First[#2]]&,KS]

Basic quadrupole

As[{K1,K2}]

Add skew-octupole

As[{K1,K2},{0,0,Ks3}]

As[{b1 (Lb1/Circ),b2 (Lb2/Circ)},{0,a2 (La2/Circ),0,0,a5 La5/Circ}]

Hamiltonian for betatron motion

Hamiltonian[x_,px_,y_,py_,KN_,KS_]:=((px^2)/2)+((py^2)/2)+(1/2) ((G^2)+K1) (x^2)-((K1 (y^2))/2)+As[KN,KS]

Cheat a wee bit since we know how it goes with the linear motion.  Come back and implement the canonical transformation from first principles sometime.

Take[Flatten[Solve[{
((Qx Jx)/R)==Hamiltonian[x,px,0,0,{0},{0}],
((Qy Jy)/R)==Hamiltonian[0,0,y,py,{0},{0}]},{px,py}]],-2]

actionAngle={px->(Sqrt[((2 Jx Qx)/R)-(G^2) (x^2)-K1 (x^2)]),py->(Sqrt[((2 Jy Qy)/R)+K1 (y^2)]),x->xc+(Subscript[D, x]) δ+(Sqrt[2 (Subscript[β, x]) Jx]) Cos[(Subscript[ϕ, x])],y->yc+(Subscript[D, y]) δ+(Sqrt[2 (Subscript[β, y]) Yj]) Cos[(Subscript[ϕ, y])]}

Simplify[Hamiltonian[x,px,y,py,{0,0,b3},{}]//.actionAngle]

Hamiltonian[x,px,y,py,{0,b2,b3},{}]//.actionAngle//Expand

Tune-shifts with amplitude and closed orbit shifts: implementation

The procedure for transforming the Hamiltonian to the action angle variables of the linear motion, performing an average over the phases and differentiating with respect to the action variables to get the detuning with amplitude.  Use dummy variables I1 and I2 in arguments so that the final derivative is computed properly and to avoid small problem with notation.

detuning::"usage"="detuning[xc,yc,δ,Jx,Jy,{b1,b2,b3,...},{a1,a2,a3,...}] returns a expression for the dependence of betatron tunes on transverse closed orbit changes xc,yc, momentum deviation δ and betatron action variables Jx, Jy in a ring with normalised normal and skew multipole components {b1,b2,b3,...} and {a1,a2,a3,...} in addition to the dipole and quadrupole fields determining the tune on the closed orbit. The result is to first order in these field perturbations.  The first component of the multipole field vectors corresponds to a dipole field.  The result contains a factor Circ which can be formally interpreted as an integral of the other factors over the circumference.";

This is coded in an ugly way for efficiency.

detuning[xc_,yc_,δ_:0,J1_:Jx,J2_:Jy,KN_List:{},KS_List:{}]:=
Module[{H,actionAngle},
actionAngle={px->(Sqrt[((2 Jx Qx)/R)-(G^2) (x^2)-K1 (x^2)]),py->(Sqrt[((2 Jy Qy)/R)+K1 (y^2)]),
x->xc+(Subscript[D, x]) δ+(Sqrt[2 (Subscript[β, x]) Jx]) Cos[(Subscript[ϕ, x])],y->yc+(Subscript[D, y]) δ+(Sqrt[2 (Subscript[β, y]) Jy]) Cos[(Subscript[ϕ, y])]};
H=Hamiltonian[x,px,y,py,KN,KS];
H=Expand[TrigToExp[H//.actionAngle]];
H=ExpToTrig[H]/.{Sin[(Subscript[ϕ, x])]->0,Cos[(Subscript[ϕ, x])]->0,Sin[(Subscript[ϕ, y])]->0,Cos[(Subscript[ϕ, y])]->0,
Sin[(Subscript[ϕ, x])+(Subscript[ϕ, y])]->0,Cos[(Subscript[ϕ, x])+(Subscript[ϕ, y])]->0,Sin[(Subscript[ϕ, x])-(Subscript[ϕ, y])]->0,Cos[(Subscript[ϕ, x])-(Subscript[ϕ, y])]->0,
Sin[_Integer (Subscript[ϕ, x])+(Subscript[ϕ, y])]->0,Cos[_Integer (Subscript[ϕ, x])+(Subscript[ϕ, y])]->0,
Sin[(Subscript[ϕ, x])+_Integer (Subscript[ϕ, y])]->0,Cos[(Subscript[ϕ, x])+_Integer (Subscript[ϕ, y])]->0,Sin[_Integer (Subscript[ϕ, x])]->0,
Cos[_Integer (Subscript[ϕ, x])]->0,Sin[_Integer (Subscript[ϕ, y])]->0,
Cos[_Integer (Subscript[ϕ, y])]->0,Sin[_Integer (Subscript[ϕ, x])+_Integer (Subscript[ϕ, y])]->0,
Cos[_Integer (Subscript[ϕ, x])+_Integer (Subscript[ϕ, y])]->0};
Expand[R {(∂_Jx)H,(∂_Jy)H}/.{Jx->J1,Jy->J2,R->(Circ/(2 π))}]]

Lots of standard results can be derived with just this function.  See next section.

Some standard results for tune-shifts, chromaticities and detuning with amplitude

?detuning

Zero-amplitude tunes

detuning[0,0,0,0,0]

Tune-shift from a quadrupole perturbation, (Subscript[K, 1]) being replaced by (Subscript[K, 1])+(Subscript[k, 1])

detuning[xc,yc,δ,(Subscript[I, x]),(Subscript[I, y]),{0,(Subscript[k, 1])}]

Chromaticity due to a quadrupole perturbation.  Note that we should include the chromatic dependence of the multipole fields.

detuning[xc,yc,δ,Jx,Jy,(({0,k1})/(1+δ))]

Expand to first order in δ

Series[detuning[0,0,δ,0,0,(({0,k1})/(1+δ))],{δ,0,1}]

To derive the natural chromaticity of a ring, we can look only at the chromatic perturbation from the normal quadrupoles

Series[detuning[0,0,0,0,0,(({0,K1})/(1+δ))-{0,K1}],{δ,0,1}]

Chromaticity from a skew quadrupole perturbation

detuning[xc,yc,δ,(Subscript[I, x]),(Subscript[I, y]),{},(({0,Ks1})/(1+δ))]

Chromaticity from sextupoles to first order in δ

Series[detuning[0,0,δ,(Subscript[I, x]),(Subscript[I, y]),(({0,0,K2})/((1+δ))),{}],{δ,0,1}]

Tune-shift from closed orbit in sextupoles (so-called "feed-down")

detuning[xc,yc,0,(Subscript[I, x]),(Subscript[I, y]),{0,0,K2},{}]

Detuning with amplitude from an octupole

detuning[0,0,0,(Subscript[I, x]),(Subscript[I, y]),{0,0,0,K3 },{}]

Second-order chromaticity from an octupole

Series[detuning[0,0,δ,0,0,(({0,0,0,K3 })/(1+δ)),{}],{δ,0,2}]

Chromaticity from orbit feed-down in an octupole

Series[detuning[xc,yc,δ,0,0,(({0,0,0,0 K3 })/(1+δ)),{0,0,0,KS3}],{δ,0,2}]

Extract the chromaticity

Coefficient[#,δ  ]&/@%

Tune-shifts with closed orbit: the DAΦNE wigglers

In the DAΦNE wigglers there is a non-zero horizontal closed orbit in the wigglers. The experiments reported in February 2001 varied the closed orbit around this value by an amount (Subscript[x, c1]) and measured the resulting tunes.    (Similar measurements at RHIC later in the year.)

Notation[(NotationTemplateTag[(Subscript[x, c0])]) ⟺ (NotationTemplateTag[xc0])]

Notation[(NotationTemplateTag[(Subscript[x, c1])]) ⟺ (NotationTemplateTag[xc1])]

With the hypothesis of a decapole field in the wigglers, the tune change expected should be

detuning[xc0+xc1,0,0,0,0,{0,0,0,0,b5}]

In these formulas the factor Circ, appearing in any term, should be interpreted as an integral of the other factors in that term around the circumference.

The coefficients of the quadratic terms are

Coefficient[#,(xc1^2)]&/@%

Translate into the notation used by Miro Preger in his transparencies to show that we agree

%/.{b5->b4,xc0->X}

A skew decapole would have no effect

detuning[xc0+xc1,0,0,0,0,{},{0,0,0,0,a5}]

What else could give us a similar effect?  Consider all possible normal and skew multipoles up to  7th order (hide the result)

{dqxc,dqyc}=detuning[xc0+xc1,0,0,0,0,{b1,b2,b3,b4,b5,b6,b7},{a1,a2,a3,a4,a5,a6,a7}];

and show the results in tidy form:

Collect[dqxc,{Circ,(Subscript[x, c1])}]

Collect[dqyc,{Circ,(Subscript[x, c1])}]

Extracting the coefficient of the second-order term in the closed-orbit shift, we find

Coefficient[dqxc,(xc1^2)]

So the effect could come from a pure octupole, the decapole with first-order feed-down from the initial closed orbit or higher order feed-down from higher normal multipoles.  Skew multipoles cannot give such a term.  No surprises! The effect on the vertical tune is the same with opposite sign and exchange of β-functions.

Coefficient[dqyc,(xc1^2)]

Coefficient[dqxc,(xc1^2)]==(-Coefficient[dqyc,(xc1^2)]/.{(Subscript[β, x])->(Subscript[β, y]),(Subscript[β, y])->(Subscript[β, x])})