Electrostatic ring for the pEDM madx, matrix and leapfrog tracking - - PowerPoint PPT Presentation

Electrostatic ring for the pEDM madx, matrix and leapfrog tracking

Mario Conte

University of Genova and INFN, Italy

Alfredo U. Luccio, Nicholas D’Imperio

Brookhaven National Laboratory, Upton, New York September 10, 2015 EUCard Workshop, Mainz, Sept./10-11 2015 1

Synopsis

We designed and tracked an electrostatic FODO ring for pEDM.

• 1. Basic design criteria for an EDM electrostatic ring;
• 2. Find the regions of stability and parameters for the ring

through a csh script that runs a modified version of CERN Madx with electric bends and quadrupoles;

• 3. Track orbits using (a)Madx matrices and then (b)to higher
• rder by differential equations kick integration with a sym-

plectic code based on Leapfrog algoritm;

• 4. Track spin dynamics by integration of the Thomas-BMT

equation with the Spink algorithm[9]. Address the issue of spin coherence; 1 EUCard Workshop, Mainz, Sept./10-11 2015

Basic design and tracking criteria

• 1. Design a ring lattice on the footprint of the AGS (800 m),

who seemed, for awhile to be a suitable site. The design can be scaled to any other site. A large footprint asks for a small value of the electric bend field;

• 2. Choose a basic FODO structure where the vertical tune is

much smaller than the horizontal, to optimize the measure- ment of the vertical spin component, proportional to the

• EDM. The FODO will not be symmetrical since cylindrical

electric bends produce an intrinsic horizontal focusing;

• 3. Design a lattice with a positive phase slip;
• 4. Adopt a convenient design package (say:

madx) to find tune islands of stability. Use a fast and symplectic tracking code (say: leapfrog);

• 5. Recur to parallel computing to calculate and optimize spin

coherence for a full beam of representative particles. 2 EUCard Workshop, Mainz, Sept./10-11 2015

Orbit/Spin tracking in an electric ring

Tracking of orbits and spin in an electrostatic storage ring for the EDM is important and should be done by more than one method to compare and benchmark. Tracking should be symplectic, stable in the long range and fast, because for EDM search ring turns will be in the billions. Keywords to keep in mind in tracking are accuracy and long range stability.. Codes proposed and used by various Authors for the tasks are based on 1-Runge-Kutta integration of differential equations for orbit (Lorentz) and spin (Thomas-BMT); 2-map description

• f machine elements (Madx) or the whole lattice (Cosy-infinity);

3-discrete kicks symplectic integration for propagation by kicks (Leapfrog, Teapot.) In this contribution we will describe a madx matrix and a Leapfrog

• rbit code of type 3, pllus the some features of a spin dynamics

code. 3 EUCard Workshop, Mainz, Sept./10-11 2015

Example of a ring designed using Madx

The code Madx distributed and maintained by CERN is a pack- age used by many accelerator designers, to optimize lattices. Madx was designed for magnetic bends. We modified it using matrices proposed by Mario Conte [1], to deal with electrostatic bends and quadruoiles. Once conceived the ring, we tested its stability by a Linux script to run Madx (by Nick D’Imperio [2]), to find islands of stability in tune space for different values of quadrupole strengths. An ex- ample table is shown. The values of βmax are growing at the edge

• f each island while the values of the betatron tune decreases.

4 EUCard Workshop, Mainz, Sept./10-11 2015

The magic condition

Spin dynamics is governed by the covariant Thomas-Bargman- Michel-Telegdi T-BMT equation ds dt = − q mγf × s, (1) where s is the real 3-dimensional spin vector of a 1/2-spin par- ticle, and f is a function of the position and the momentum of the particle and of the electric and magnetic field encountered by the particle along its trajectory. Spin is a passenger on the orbit, that has to be calculated first at each step. In a pure electrostatic ring, e.g. with no magnets

• r RF cavities, f reduces to

f = γ

• a −

1 γ2 − 1

E × v

c2 , (2) with a the spin anomaly. At the magic momentum pc = mc2/√a it is exactly f = 0 and the spin remains frozen in its direction at injection (longitudinal) respect to the orbit. 5 EUCard Workshop, Mainz, Sept./10-11 2015

Issues for electric accelerator lattices

An electrostatic lattice behaves differently than a classical mag- netic lattice. In an e.element the kinetic energy of a particle is modulated, while in a m.element it is not, since in the Lorentz equation of motion dp/dt = e(E + v × B) (3)

• nly the vector term parallel to the momentum appears while in

a m.lattice it is the vector product term present, with the force perpendicular to the momentum. In the present design we adopted simple cylindrical electrodes for the bends, that produce only a radial field far from edges. While conventional no-gradient magnetic bends don’t focus the beam, electric cylindrical bends produce a small horizontal focusing, so that to produce a FODO condition the focusing and defocusing quadrupoles should be slightly different. Vertical focusing will be obtained with electrostatic quadrupoles. An option would be to provide electric bends with also a ver- tical curvature of the electrodes. For the moment we are not considering this option because it is more hard and expensive to construct with the desired accuracy. 6 EUCard Workshop, Mainz, Sept./10-11 2015

We considered a ring of 800 m length (taylored on the BNL-AGS tunnel, that was for awhile proposed as a possible site for the pEDM ring.) 72 bends of 9 m length. 80 FODO quadrupoles of 2×0.5 m length, 4 drifts of 2×9 m. Values of basic parameters are listed below 7 EUCard Workshop, Mainz, Sept./10-11 2015

The ring, with 4×18 bends an 8 straights

February 2015

72 bends− 4 double drifts each 9 meter 80 FODO quadrupoles ’y2’ pEDM electrostatic ring total ring length = 800 m

D D F F D F D F F F D D D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D F D D D F D D F F D D F D F D F D F D F D F F F F

drift bend

Mario Conte Universit‘a di Genova, Italy Alfredo U Luccio Brookhaven National Laboratory

8 EUCard Workshop, Mainz, Sept./10-11 2015

Electrostatic bend

El.static bend matrix α =

• 2 − β2, a1 = cos(αθ), a2 = (ρ/α) sin(αθ), a3 = −(α/ρ) sin(αθ),

M =

a1

a2 a3 a1 1 Lb 0, 1

• , Lb = length of bend

. Other matrices, say: for drifts, are Madx’s 9 EUCard Workshop, Mainz, Sept./10-11 2015

Electrostatic quadrupole

gap = 2a = 10cm, Lq = 0.5m(quad length) El.gradient :GE = 2V0 a2 , k = e mc2 GE β2γ, V0 = a2 2 mc2 e (β2γ)k El.static quadrupole matrices MF =

• cos(

√ kLq) 1 √ k sin( √ kLq) − √ k sin( √ kLq) cos( √ kLq)

• MD =
• cosh(

√ kLq) 1 √ k sinh( √ kLq) √ k sinh( √ kLq) cosh( √ kLq)

• 10

EUCard Workshop, Mainz, Sept./10-11 2015

Value (MKSA) and dimension of all ring parameters are magic proton of β = 0.59837912., γ = 1.24810740, β2γ = 0.44689430, k = 0.043, V0 = 1.42245166.105V, Eqin the quads = 2.844903.106V/m, Ebin the bends = 2.54972867.106V/m Compared with the Stability Table, the above shows that the working γ of this particle is less than γT and the vertical betatron tune can be much smaller than the horizontal, as we want, for EDM measurements. About the field in the quads, note that for the optics the quantity

• f importance is

√ kLq with Lq the length of the quadrupole. The field in the quadrupole is proportional to k. Therefore increasing the length of the quadrupole, but at the same time decreasing k and keeping √ kLq constant, can effectively reduce the field. 11 EUCard Workshop, Mainz, Sept./10-11 2015

The matrix driven tracking program

We started with a first order madx matrices tracking

loop

MTRACK INIT DRIFT QUAD BEND DBDR BCELL DBDR ENERMOD MATRICES RAYTRANSFER (WRITE) FINIS GLOBAL ¨uber alles

MTRACK3

loop

flowchart 12 EUCard Workshop, Mainz, Sept./10-11 2015

Betatron oscillations by matrix tracking

13 EUCard Workshop, Mainz, Sept./10-11 2015

Orbit tracking by leapfrog

Going beyond first order orbit tracking, consider canonical in- tegration of the Lorentz diff. equation of motion, Eq.(6) by leapfrog or Verlet[4] kicks, method invented for astronomy by Delambre[5] in 1792, and adapted for accelerators by Ronald Ruth[6]. It is a kick integration method that interleaves drifts, where only space coordinates are advanced, with symplectic kick bends where the momentum components are advanced. Leapfrog is an algorithm accurate to 2.nd order in time step. Teapot, by R.Talman and L.Schachinger [6] it is similarly con-

• structed. See also [7].

Other integration algorithms, like Runge-Kutta are accurate to 4.th order in time. However they were written with mathemat- ical accuracy in mind, while the 2.nd order Leapfrog is exactly symplectic, i.e. was written with physical accuracy in mind. Symplectic Runge-Kutta has been discussed[8]. It makes a com- puter code slower to run, which defies our goal of short computer time for tracking an EDM ring for so many turns. 14 EUCard Workshop, Mainz, Sept./10-11 2015

Notes on the algorithm ’order’

For spin tracking, in particular to address the issue of spin co- herence (see later sections), orbit tracking algorithms would be

• f high order. Madx matrices are first order, leapfrog is of higher
• rder. How higher? Leapfrog works with calculation of the orbit

in the electric field of bends and quadrupoles, then we should use high order power expansions for the electric field, in (x, y, z) [11], and for the Hamiltonian, whose constant value should be used for a steady check of the symplecticity of the algorithm. The Hamiltonian contains the electric potential φ obtained as a power expansion solution of the Laplace equation ∇2φ = 0. (4) In the present simulation we expanded field and potential ex- pression to the third order. Also leapfrog integration should be done at high order, see R.Ruth [5]. Of course all this will slow down orbit calculation tracking and should be carefully considered. 15 EUCard Workshop, Mainz, Sept./10-11 2015

Leapfrog coordinates

While for madx matrix tracking we used Fr´ enet-Serret coordi- nates, for Leapfrog one uses Cartesian ”laboratory” coordinates (x, z, y),with ˆ y vertical axis, and time as the independent vari-

• able. Vertical electric field component is calculated by a power

expansion out of the ”horizontal” x, z plane of the ring The circular (plus straights) ring lattice shown is obtained by tracking a ”reference par- ticle” i.e at magic energy in- jected tangentially. Note that while a matrix method tracks trough a pre-designed lattice, leapfrog ”designs” the lattice by tracking a reference (magic) particle. 16 EUCard Workshop, Mainz, Sept./10-11 2015

Orbit Leapfrog formalism basics

M´ enagerie of quantities for the game ρ[m] = radius of curvature a = magnetic anomaly Uo[GeV ] = mc2, mass − energy ℘ ≡ pc[GeV ] = Uo/√a, momentum UT[GeV ] =

• ℘2 + U2
• , total energy

γ = UT/mc2, β =

• 1 − 1/γ2

Bρ[V · s/m] = 109℘/c, rigidity eE[eV/m] = (= ℘/r0)βc Electric bend field Leapfrog tracking conserves the value of the Hamiltonian, that is being continuously recalculated by series expansion from the Hamilton Equation ∇2φ = 0 during runs. H =

• (℘ − eA)2 + (mc2)2 + eφ.

(5) 17 EUCard Workshop, Mainz, Sept./10-11 2015

Will work on the EDM ring described above, Basic leapfrog cell for a bend is a sequence drift + momentum kick + drift Momentum kick are done by integration of the Lorentz equation, for an electric bend: dp dt = eE, with E = −∇φ (6) The potential, needed for the Hamiltonian, should obey the Laplace equation ∇2φ = 0. (7) Explicit expressions for φ and A are found by power expansion. The reference particle, around which the whole beam oscillates, is the magic particle whose spin would remain frozen in longitu- dinal position during the propagation. 18 EUCard Workshop, Mainz, Sept./10-11 2015

Leapfrog cell-Electric bends

Let us discuss what happens to a reference particle confined to the horizontal plane of the FODO. For simplicity we draw only 3 instead of 18 kick bends in a ring quadrant.

150714−1 θb ρ z x d d C A B a b a d

• u

b l e c e l l F E D G O

Fig.1 AOB → C→ BC,CD drift, kick-bend, LFdrifts C, center of curvature 19 EUCard Workshop, Mainz, Sept./10-11 2015

drift AOB

call ℘ = pc , in [GeV]. Start in A with initial coordinates (A) x = −xo, z = −zo, ℘x = 0, ℘z = ℘. where zo = −(d + ρ), see the figure. Eq’s for the drift, with a time step dt for the drift A→B: dx dt = ℘x Uoγc, dz dt = ℘z Uoγc, or x := x + ℘x/(Uoγ)c dt, z := z + ℘z/(Uoγ)c dt (8) using the identity ℘ = Uoβγ, we obtain at the kick bend C the new position (B) x = ro, z = βc dt, ℘x = 0, ℘z = ℘. 20 EUCard Workshop, Mainz, Sept./10-11 2015

kick in B

In B a kick is imparted to the momentum pc, using the Lorentz Equation, with a time step δt, different from the dt of the drift. ℘x := ℘x − eExc δt, ℘z := ℘z − eEzc δt (9) For cylindrical bend the field E is purely radial, with components eEx = −eE (ro/r) cos θ eEz = eE (ro/r) sin θ. (10) Now find the relation between dt and δt for leapfrog i.e:

• 1. Through the bend the value of the total momentum pc

must be conserved,

• 2. The trajectory in C should return tangent to the circle, as

in the figure. Namely: arccos (p(A) · p(C))/p2 = 2θ (11) If both conditions hold, the basic trajectory will be a polygon circumscribed to the circle. Other particles in the beam will dance around it in betatron oscillations. 21 EUCard Workshop, Mainz, Sept./10-11 2015

For condition (1): momentum conservation, combining the pre- ceding equations ℘x = −℘/r cos θ βc δt, ℘z = ℘ (1 − (1/r) sin θ βc δt) (12) then after kick (C): ℘2

x + ℘2 z = (pc)2

1 + ((βc/r)δt)2 − (2/r) sin θ βc δt . (13) Since: cos θ = z/r, sin θ = x/r, taking the value of x from Eq.(8), the term in [ ] in Eq.(13) above reduces to 1 when δt = 2 dt The time at the kick bend should be the same as for both straights across it ! For condition (2): new angle, it is calculated from the scalar product of the momentum before and after the kick

• (A) before kick: ℘x = 0, ℘z = ℘
• (C) after kick: ℘x = −(℘/r) cos θ βc δt, ℘z = ℘

1 − 2 sin2 θ angle = arccos ℘(A) · ℘(B) (pc)2 = arccos 1 − 2 sin2 θ = 2θ q.e.d. The constructed orbit remains tangent to the arc ! 22 EUCard Workshop, Mainz, Sept./10-11 2015

Reference Trajectory

Let us produce a reference trajectory on the horizontal plane by leapfrog tracking a magic particle along a polygonal pattern tangent to a structure made of straights (drifts) and circular arcs (bends). The leapfrog orbit is slightly longer than the reference

• rbit. The more kicks we put in a bend the lesser this difference

is. For our structure with 8 bends and 8 drifts of circa 270 m of total length, using 32 kicks in each bend of 103 m of radius, the difference in effective radius between the geometrical base line and the polygon is about 1 mm. The step is much larger than the required step of a solution by integration for similar accuracy, with a very large gain in computing speed. Tracking a reference particle will create a reference trajectory. An example is shown in the following picture. 23 EUCard Workshop, Mainz, Sept./10-11 2015

Reference Trajectory by tracking Ring with eight bends and eight straights

Fig.2 32 kicks per bend bend length=28.276 m drift length 2 × 2.83 m intra bend drift length = 0.44 m nominal curvature radius = 36 m Ecyl = −1.1647455107V/m 24 EUCard Workshop, Mainz, Sept./10-11 2015

Evaluation of the electric field

• d

D D d D d D d θ = π/4 θ = π/4 r

• r

Fig.3 In a general lattice the center

• f curvature for the calculation
• f the electric field

continuously changes and has to be re-evaluated every time (in the present case the center

• f curvature changes from quadrant

to quadrant) 25 EUCard Workshop, Mainz, Sept./10-11 2015

General tracking by leapfrog

The Leapfrog formalism extends to 3 dimensions and applies unchanged to particles that don’t have a magic energy or are injected in the lattice on a finite transverse emittance. Eqs.(8) and .(9) in 3 dimensions are

x := x + ℘x/(Uoγ)c dt

y := y + ℘y/(Uoγ)c dt z := z + ℘z/(Uoγ)c dt , ℘x := ℘x − eEx2c dt ℘y := ℘y − eEy2c dt ℘z := ℘z − eEz2c dt . (14) However, In a general case the leapfrog conditions (1) for mo- mentum and angle are not fully satisfied in a bend because, due to transverse oscillations, the particle sees a tangential compo- nent of the electric field that modulates the energy. During tracking the Hamiltonian is continuously calculated. It conserves its initial vaiue. 26 EUCard Workshop, Mainz, Sept./10-11 2015

Orbit coordinates

26 EUCard Workshop, Mainz, Sept./10-11 2015

Fig.4 x,y betatron oscillations vs. turn # 27 EUCard Workshop, Mainz, Sept./10-11 2015

Add a RF - Example of RF bucket

Fig.5 - Phase space of ∆ × pc for two particles, with dp/p = 1.e−4 and 2.e−4, with VRF = 1000V and h = 24. Number

• f turns for a complete oscillations is 335, corresponding

synchrotron frequency νs = 0.002985 oscillations per turn 28 EUCard Workshop, Mainz, Sept./10-11 2015

Code Spink for Spin Dynamics package

The package used for spin dynamics in our simulation is the code Spink [9], developed and used for many years for the AGS and RHIC at Brookhaven, that calculates the evolution of the particle spin. In Spink the spin in described as a 3-dimensional real vector (sx, sy, sz). The code uses a unitary matrix formalism, to in- tegrate by kicks the T-BMT covariant differential equation, Eq.(1), that we repeat here ds dt = − q mγf × s, (15) Elements of the function (matrix) f are expressed (non linearly) as a function of the instantaneous values of the position and momentum of the particle, so we say that spin dynamics rides

• n top of orbit dynamics.

Properties of the matrix, in particular the use of its eigenvectors and eigenvalues, as described in the following, are important to calculate spin coherence. 29 EUCard Workshop, Mainz, Sept./10-11 2015

Spin Dynamics: for EDM search

Spin kicks, applied at each Bend and Quad, follow the leapfrog pattern of the orbit. At the magic energy we have f = 0 and the spin remains frozen in its longitudinal direction from injection on. If the proton has an EDM, the spin is NOT completely frozen: by Relativistic transformation in the rest frame of the particle, the electric field appears as a magnetic field B’ ⊥ to E

B′ = γ

• B −
• β

c × E

• −

γ2 γ + 1( β · B) β. (16) Another small term is added to f in Eq.(15)

B′ = −γ

β × E. f := f + ηB′ × v. (17) The spin will make a precession around this magnetic field and a spin vertical component will appear, proportional to the electric dipole moment, that we want to measure. For a magic proton this is the only non vanishing additional spin component. 30 EUCard Workshop, Mainz, Sept./10-11 2015

Spin dynamics of a longitudinally frozen spin

Fig.8 - Longit. component of the frozen spin: red line in ac- celerator (Fr´ enet-Serret) coordinates, green line, in laboratory coordinates. The red line shows little wiggles because the re- sponsible proton is in this example on purpose not perfectly magic and there are also betatron oscillations. 31 EUCard Workshop, Mainz, Sept./10-11 2015

Spin coherence

For EDM measurements the spin of the beam should remain coherent (the beam remains polarized as an ensemble) for a measurement time of the order of milliseconds. Spin coherence is conveniently measured as the reciprocal of the width of a spin tune line, built up from the spin tune of a representative large ensemble of the particles in the beam (say, 256), tracked at the same time on a computer cluster. For this we wrote a Linux csh script to:

• 1. generate a distribution for a polarized particle beam of

given emittance, energy spread and polarization under dif- ferent distributions, depending on the beam injector;

• 2. run the tracking code in parallel for all the particles using

the Message Passing Interface (MPI) Library;

• 3. calculate as a running average of spin components the spin

tune line. Spin coherence time is the reciprocal of the spin tune linewidth. It is efficiently controlled using sextupoles in the ring lattice. We did a very successful simulation of spin coherence with the code UAL-Spink on polarized runs on COSY with a RF solenoid [10]. 32 EUCard Workshop, Mainz, Sept./10-11 2015

Spin tune

Spin tune, number of spin precession per machine turn, is cal- culated as an eigenvalue of a one-turn spin 3 × 3 matrix M, calculated in turn when the orbital 6 dimension phase space vector regains its starting value. Instead of waiting very many turns for this to happen, spin tune is equivalently calculated by running average of the spin tune in each orbit turn, over typically 15,000 turns. Calling T = M11 + M22 + M33 the trace of the spin matrix, the fractional part of spin tune Qspin or its complement to one is calculated as ˜ Qspin = µ 2π, µ = arccos

1 − T

2

• ,

(18) depending on the sign of µ. From other matrix elements, this sign, and the two Euler angles of the orientation of the spin precession axis, θ, latitude, and φ, longitude, are

    

φ = arctan(2) ((M12 − M21), (M23 − M32)) θ = arctan(2) ((M23 + M32) sin φ, (M13 + M31)) sin µ = M12 − M21 2 cos µ sign(µ) = − arctan(2)(sin µ, cosµ) (19) if(sign(µ) <= 1){ ˜ Qspin = 1 − ˜ Qspin} (20) 33 EUCard Workshop, Mainz, Sept./10-11 2015

References

• 1. M.Conte Electrostatic Storage Ring

University of Genova, PACS 29.20.db Oct.24, 2012

• 2. Nicholas D’Imperio (BNL) internal communication
• 3. L.Verlet Computer experiments in classical fluids

Phys.Rev. 159 98-103

• 4. J.B.Delambre Tables du soleil 1792
• 5. R.D.Ruth A Canonical Integration Technique

IEEE Trans. on Nuclear Sciences NS-3, No.4, August 1983

• 6. L.Schachinger and R.Talman

Teapot: A Thin-element Accelerator Program Particle Accelerators, 22 1987

• 7. J.M.Sanz-Serna

Symplectic integrators for Hamiltonian problems Acta Numerica (1991), 243-286 QA 297,A21 1992

• 8. W.H.Press et Al. Numerical Recipes Cambridge Uni.Press,1992
• 9. A.U.Luccio Brookhaven National Laboratory, July 30, 2007

Spink User’s Manual. Version v.2-21-beta

• 10. A.U.Luccio, E.Stephenson et Al

Proceedings of the the Bad Honnef Meeting on EDM, July 2012

• 11. Y.Orlov and Y.Semertzidis

private communication, internal notes 34 EUCard Workshop, Mainz, Sept./10-11 2015