NICH AN S E UNIVERSITY LIBRARIES IlllllllWilli!llm1‘l‘lll ”I I 3 1293 01028 9217 LIBRARY Michigan State University I This is to certify that the dissertation entitled RIGOROUS BOUNDS 0N SURVIVAL TIMES IN CIRCULAR ACCELERATORS AND EFFICIENT COMPUTATION OF FRINGE-FIELD TRANSFER MAPS presented by Georg Heinz Hoffstétter has been accepted towards fulfillment of the requirements for Ph.D. Physics degree in H owls“ 1% Major professor 9-27-1994 Date MS U is an Affirmative Action/Equal Opportunity Institution 0-12771 PLACE N RETURN BOXtomnovothi-chockwttmn yourncord. TO AVOID FINES return on or baton date duo. DATE DUE DATE DUE DATE DUE MSU ioAn Affinndivo Adlai/Equal Opportunity Instituion RIGOROUS BOUNDS ON. SURVIVAL TIMES IN CIRCULAR ACCELERATORS AND EFFICIENT COMPUTATION OF FRINGE—FIELD TRANSFER MAPS By Georg Heinz Hoffstatter A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Physics and Astronomy 1994 ABSTRACT RIGOROUS BOUNDS ON SURVIVAL TIMES IN CIRCULAR ACCELERATORS AND EFFICIENT COMPUTATION OF FRINGE-FIELD TRANSFER MAPS By Georg Heinz Hoffstatter Analyzing stability of particle motion in storage rings contributes to the general field of stability analysis in weakly nonlinear motion. A method which we call pseudo invariant estimation (PIE) is used to compute lower bounds on the survival time in circular accelerators. The pseudo invariants needed for this approach are computed via nonlinear perturbative normal form theory and the required global maxima of the highly complicated multivariate functions could only be rigorously bound with an extension of interval arithmetic. The bounds on the survival times are large enough to be relevant; the same is true for the lower bounds on dynamical apertures, which can be computed. The PIE method can lead to novel design criteria with the objective of maximizing the survival time. A major effort in the direction of rigorous predictions only makes sense if accurate models of accelerators are available. Fringe fields often have a significant influence on optical properties, but the computation of fringe—field maps by DA based integration is slower by several orders of magnitude than DA evaluation of the propagator for main—field maps. A novel computation of fringe—field effects called symplectic scaling (SYSCA) is introduced. It exploits the advantages of Lie transformations, generating functions, and scaling properties and is extremely accurate. The computation of fringe—field maps is typically made nearly two orders of magnitude faster. This thesis is dedicated to Johann Sebastian Bach and to his inspiring dedications. iii ACKNOWLEDGEMENTS All too often acknowledgements seem to be misused to emphasize the quality of the presented work by mentioning invitations, journeys, and support; mostly before thankfulness towards the most important supporters of private life is expressed. In order to clarify priorities and to separate my gratitude for selfless work of others from acknowledgements of professional support, which tends to shed light on my own work, I will distinguish between private and professional acknowledgements. Private: Deep reverence goes to my wife Anke, who, with happiness, endured my education and the inconveniences it brought with it. Without doubt, the liveliness she and our three children, Samuel, Lydia, and Tabea bring to our home contribute much to my work. The time and energy spent by my parents Heinz and Elisabeth for fostering my general curiosity and my interest in history is deeply appreciated. Although his name has the touch of professionalism, I want to express my gratitude towards Prof. Dr. Martin Berz here for being a personal friend and advisor and for shaping my scientific personality in many ways. Much of the idealistic and financial support mentioned below was initiated by him. It is due to his support that my family could settle down in Michigan quickly, and the academic group gathered around him contributed much to my stimulating experience at the National Superconducting Cyclotron. The members of this group, Meng Zhao, Khodr Shamseddine, Weishi Wan, Dr. Ralf Degenhardt, and Kyoko Fuchi all share a part of the credit of the presented work due to their discussions and friendship. I am also deeply indebted to Prof. Dr. Harald Rose for the energy and happiness with which he introduced me to physics research and guided me towards the field of particle optics. iv Professional: The financial support of The German National Scholarship Founda- tion (Studienstiftung des Deutschen Volkes) has made my education, starting from undergraduate times, much smoother than it would have been otherwise. And it is also from their generous support that I could visit Prof. Dr. Martin Berz’s group, which finally lead to the completion of this dissertation. This support was followed up by generous funds due to a National Superconducting Cyclotron Fellowship and a Scholarship of the College of Natural Science at Michigan State University. I am also thankful for the Quadrille Ball Fellowship of the Germanistic Society of Amer- ica. The financial support for participating in the conference on Numerical Analysis with Automatic Result Verification (Lafayette, LA 1993) and in the Fourth Interna- tional Conference on Charged Particle Optics (Tsukuba, Japan 1994) granted by the conference committees is appreciated. The visit to the 1992 United States Particle Accelerator Summer School at Stan- ford University and the received instruction and encouragement by Prof. Dr. Karl Brown and Prof. Dr. Roger Servranckx became important to my work. The real fun of a scientific career starts when others become interested in your work and I am therefore very thankful to the conference committee of the 1994 Crystal City APS meeting to invite me for a talk on the computation of fringe—field effects. Many thanks go to Dr. Desmond Barber, who invited me to DESY to introduce him into the code COSY INFINITY and to give a talk about my work. For the same reason thanks go to Prof. Dr. Harald Rose, who invited me to a seminar talk on the fringe—field aspects of the work which will be presented. Many aspects of interval arithmetic became clear to me in discussions with Dr. Dietmar Ratz and Olaf Kniippel from Karlsruhe and Hamburg respectively, who I thank for the time they took to introduce me into the tricks of interval optimization. Contents LIST OF TABLES LIST OF FIGURES Foreword 1 Normal Form Theory 1.1 Normal Form Theory for Order n Symplectic Maps .......... 1.2 First Order Transformations ....................... 1.3 Nonlinear Transformations ........................ 1.4 Invariants ................................. 1.5 Real Normal Forms ............................ 1.6 Action—Angle Variables .......................... 1.7 Pseudo Hamiltonians ........................... 1.8 Parameter Dependent Normal Forms .................. 1.9 General Normal Form Theory ...................... 2 Long Term Estimates for Weakly Nonlinear Motion 2.1 Introduction to the PIE Method ..................... vi ix xi 11 1'2 25 27 35 42 44 46 50 51 2.1.1 History ............................... 2.2 Pseudo Invariant Estimation (PIE) ................... 2.3 Normal Form Transformations and Pseudo Invariants ......... 2.4 Maps to Test the Method ........................ 2.4.1 The Physical Pendulum . . . .' ................. 2.4.2 Henon Map for One Degree of Freedom ............. 2.4.3 Coupled Pendulum ........................ 2.4.4 The IUCF Ring .......................... 2.4.5 The PSR II Ring ......................... 2.4.6 The Demo Ring .......................... 2.5 The Pseudo Invariant and How to Parameterize Regions ....... 2.6 Analysis of Pseudo Invariants ...................... 2.7 Influence of Resonances .......................... 2.8 Symplectic Representations ....................... Optimization by Scanning 3.1 Refinements and Examples ........................ 3.1.1 Dividing Phase Space ....................... 3.1.2 Time Evaluation ......................... 3.1.3 Parameter Dependence ...................... Rigorous PIE with Interval Arithmetic 4.1 Interval Arithmetic ............................ vii 52 54 58 62 63 65 66 66 96 97 101 102 4.2 Interval Chains and Optimization .................... 106 4.3 Parameterizing Regions for the IC—PIE Method ............ 112 4.4 Comparison Between Intervals and Interval Chains .......... 113 4.5 Refinement of the Rigorous Estimates .................. 116 4.5.1 Analysis of Multi—Turn Maps .' ................. 117 4.5.2 Random Parameters ....................... 123 4.6 Using Taylor Maps with Remainder Bound ............... 124 5 Symplectic Scaling (SYSCA) 126 5.1 Integration in DA ............................. 129 5.2 Evaluating Propagators in DA ...................... 130 5.3 Desirable Properties of an Arbitrary Order‘ Fringe—Field Approximation 135 5.3.1 Order n Symplecticity ...................... 136 5.3.2 Accuracy for a Wide Range of Apertures ............ 138 5.3.3 Usability for Arbitrary Orders .................. 139 - 5.4 The Principles and Usefulness of SYSCA ................ 139 5.5 Examples ................................. 148 5.6 The Lie Exponent ............................. 154 5.7 The Generating Function for the Linear Map .............. 156 5.8 Transformation between Cartesian and Canonical Coordinates . . . . 159 LIST OF REFERENCES 162 viii List of Tables 2.1 2.2 2.3 2.4 2.5 2.6 2.7 3.1 3.2 3.3 4.1 4.2 4.3 4.4 4.5 Smallest resonance denominators for the physical pendulum and the Henon map with a tune of 0.379. .................... Resonance denominators for the map of the coupled pendulum. Resonance denominators for the map of the IUCF ring ......... Resonance denominators for the map of the PSR 11 ring. ....... Resonance denominators for the map of the Demo ring ......... Maximum of the deviation function after symplectifying the map by adding higher orders. ........................... The maximum of the deviation function as a function of the turns when the map is represented by a generating function ............. Minimum number of turns required to move from the initial region to the forbidden region. The initial emittance was chosen to be half the acceptance. The maps were evaluated in order eight, whereas the pseudo invariants were computed to the indicated order. Scanning was performed with 20" points for k relevant dimensions. ......... Number of map applications used for improving the estimation. Lower bounds on the turns of particles for initial emittance of one half the acceptance and lower bounds on the stable emittances for 108 turns obtained by various variations of the PIE method. .......... Predictions of the number of stable turns as a function of the order of the polynomials describing the normal form transformation for the physical pendulum with a maximum elongation of 1/10 rad. Because of energy conservation, the map is known to be permanently stable. . Predictions of the number of stable turns for the Henon map at tune 0.13, strength parameter 1.1, and starting position (:r,a) = (0.01, 0) as a function of the order of the normal form transformation. ...... Predictions of the number of stable turns as a function of the order of the approximate invariant for the Los Alamos PSR II storage ring for the motion in a phase space of 100 mm mrad. ............. Decreasing width of the remainder interval with evaluation order. Remainder intervals for multi—turn maps ................. ix 74 76 77 78 79 89 90 94 97 100 114 115 115 121 122 4.6 4.7 4.8 5.1 Lower bounds on the turns of particles for an initial emittance of one half the acceptance, obtained by the IC—PIE method. ........ 123 Maximum error of the Taylor map in the phase space region of interest. 125 Results of the IC—PIE bounds on the survival time of particle motion for rigorous description of the systems by a Taylor map with remainder intervals ................................... 125 Tilt angle and opening aberration for various fringe—field models. . . 151 List of Figures 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 a) The initial region 0 and the allowed region A of phase space 7’ with O C A C 'P. b) The gap Af that has to be bridged ........... Phase Diagram for 2000 turns in an accelerator for four initial condi- tions. The left picture shows the motion displayed in standard particle optical coordinates x and a, and the right picture shows the same mo- tion in normal form coordinates ...................... 9 particles tracked for 500 turns through the Henon map. The starting conditions were a: = 0.1 ~j, j E {1, . . . ,9} and the kick strength was 1.1. The tune in figure a) is 0.255 and in figure b) 0.29 ......... 9 particles tracked for 500 turns through the coupled pendulum. The starting conditions were angles :c/lo = 0.01 -j forj 6 {1, . . . , 9}. This corresponds to a maximum emittance of 4050mm mrad. The horizontal and vertical tunes in figure a) are 0.36 and 0.78, and in figure b) 0.37 and 0.78 ................................... a) Layout of the IUCF ring, only dipoles and quadrupoles are shown. b) :r—pl. phase space positions for 500 turns. The initial conditions are (a: = 5j,y = 5) ' 10“ withj E {1,...,10}. ............... a) Layout of the PSR 11 ring, only dipoles and quadrupoles are shown. b) x—p, phase space positions for 500 turns. The initial conditions are (:1: = 5-j,y =5-j)-10‘3 withj E {1,...,8}. ............. a) Layout of the Demo ring, only dipoles and quadrupoles are shown. b) z—px phase space positions for 500 turns. The initial conditions are (a: = 5-j,y =5-j)-10'3 withj E {1,...,8}. ............. The motion on nonlinear invariants in the phase space section x-px in figure a) and y—py in figure b). The allowed and the forbidden region and the definition of Hill is depicted in figure c) ............. Fluctuation of the pseudo invariant on the surface of the allowed region. Elliptic beam shapes (3 = o) and rectangular beam shapes (3 = CI) were used for the left and the right picture respectively. The four cases and the emittances used are the coupled pendulum (£1 = 0.66, 62 = 0.46), the IUCF ring (61 = 0.36, 62 = 0.76), the PSR II (61 = 0.056, 62 = 0.956), and the Demo ring (61 = 0.46, 62 = 0.66).' ............ Invariant defect for the physical pendulum on an invariant ellipse. The coordinates are 0 to 6 from left to right and 0 to 27r from front to back. The range of the depicted function is [—2.36 - 10"”,563- 10‘12]. . . . xi 57 60 63 64 66 67 68 71 72 74 2.11 2.12 2.13 2.14 2.15 2.16 2.17 2.18 2.19 2.20 2.21 3.1 3.2 Invariant defect for the Henon map on an invariant ellipse. The coor- dinates are 0 to 6 from left to right and 0 to 21r from front to back. The range of the depicted function is [—5.59- 10’8,5.56 - 10‘8]. . . . . 75 Invariant defect for the coupled pendulum on a pseudo invariant torus. The coordinates are in [0,21r] x [0,27r]. The range of the depicted function is [-2.55 - 10‘”, 7.36 - 10‘13] ................... 77 Invariant defect for the map of the IUCF ring on a pseudo invariant torus. The coordinates are in [0, 27r] x [0, 27r]. The range of the depicted function is [—6.84 - 10‘9,6.00 - 10‘9]. .- ................. 78 Invariant defect for the map of the PSR II ring on a pseudo invariant torus. The coordinates are in [0, 27r] x [0, 27r]. The range of the depicted function is [—1.76 - 10'9, 1.65 - 10‘9]. .................. 79 Invariant defect for the map of the Demo ring on a pseudo invariant torus. The coordinates are in [0, 27r] x [0, 27r]. The range of the depicted function is [—7.22 - 10'9, 6.95 - 10‘9]. .................. 80 Regions of resonances for two degrees of freedom up to order a) 6, b) 12. 82 Variation of the maximum 6 of the deviation function d f with tune for the Henon map of a) order 8, b) order 4. The scale is logarithmic and inverted. Staying out of resonance holes yields long survival time predictions with the PIE .......................... 84 Variation of the maximum 6 of the deviation function (1; with tune for the 8th order map of the physical pendulum. The scale is loga- rithmic and inverted. Due to the existence of an invariant, resonance denominators do not cause a problem during normal form computation. 85 Variation of the maximum 6 of the deviation function df with tune for the 8th order map of the PSR II. The scale is logarithmic and inverted. The whole horizontal and vertical tune space is covered: V; x 11,, 6 [0,1] x [0,1]. The range of d; is from 2.2 - 10'10 to 1.9- N”. 86 Variation of the maximum 6 of the deviation function 61; with tune for the 8") order map of the Demo ring. The scale is logarithmic and inverted. The range of d; is from 1.3 - 10‘9 to 1.9 - 10". ....... 88 Quality of normal form invariants when the map is represented by a generating function for the four examples with two degrees of free- dom: a) the coupled pendulum, the range of the displayed function is [-2.15, 9.45]-10‘13, b) IUCF: [—1.16,1.11]-10“8, c) PSR II: [—2.63, 2.63]- 10‘9, d) Demo: [—6.36,8.49] - 10"). For the order 8 symplectic Taylor maps, the corresponding figures were displayed previously. ...... 91 Figure a) depicts F (:r) and figure b) D(a:) in the allowed region. The variation of the pseudo invariant relative to Hill is shown in figure c). 93 Dividing phase space approximates the area under the displayed curve by several rectangles to obtain the bound on the survival time. Without this improvement, the PIE method would only give the area of the drawn rectangle. ............................. 96 3.3 3.4 4.1 5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11 a) slow increase of the pseudo invariant over many turns, b) periodic change of the pseudo invariant after relatively few turns. ....... The figures display the variation of the deviation function in the al- lowed region for the six examples. From top left to bottom right these are a) Pendulum, b) Henon map, 6) Coupled pendulum, d) IUCF, e) PSR II, f) Demo. ............................. Top left: The function of which an upper bound should be found. Others: Covering the region of interest_with smaller intervals reduces the overestimation of the maximum .................... Definition of fringe—field and main—field maps. ............. Time consumption of integration in DA relative to main—field map computation with the propagator for different expansion orders. From top left to bottom right: dipole, quadrupole, hexapole, octupole. . . . Accuracy of SCOFF for different orders, a random choice of Taylor coefficients would be close to one. From top left to bottom right: dipole, quadrupole, hexapole, octupole .................. Tracking through a 15° dipole, displayed in conventional coordinates on the left and in normal form coordinates on the right. The figures on the top display accurate tracking data, those in the middle were created by the 5th order Taylor map, and the bottom pictures were produced by 5th order symplectic generating function tracking ..... The coordinates of the particle trajectories in two elements scale with the factor a if the elements scale with the factor a and the fields scale with the factor l/a ............................. Schematic outline of the procedures leading to the reference represen- tation and from the reference representation to the general transfer map. .................................... left: (zlxza) for a quadrupole as a function of the field at the pole tip. right: Error A of the approximation of (:rlxma) with different expansion orders for the reference representation at B=2T ....... Factor of time advantage of SYSCA to numerical integration with ac- curacy of 10'”9 as a function of the expansion order. From top left to bottom right: dipole, quadrupole, hexapole, octupole .......... Relative deviation of predicted field settings with SCOF F and SYSCA from the correct settings for five quadrupoles. The standard fringe— field approximation of TRANSPORT is given as a reference; the devi- ation is mainly due to the neglect of quadrupole fringe fields. Beam spots with SYSCA (left) and SCOFF (right) approximation. The plot produced with the exact fringe fields can not be distinguished from the plot produced with SYSCA. Note the difference in scale. . 500 turn tracking without (left) and with (right) fringe fields trough an 1Gev proton storage ring ........................ xiii 97 98 107 127 132 134 137 142 147 148 150 151 152 152 5.12 5000 turn tracking with fringe fields obtained by numerical integration (left), SYSCA (middle), and a non—symplectic fringe—field approxima- tion (right). The initial position of the particle is (x, y) = (3cm,3cm) with no initial inclinations z’ and y’. .................. xiv Foreword Predictions of accelerator performance require an accurate model of particle optical devices and accurate computation schemes to evaluate the motion of particles in each element of the accelerator. Usually it is considered an accurate description of an accelerator when the motion of a single particle can be computed with a small error. However, having such an accurate model available only allows one to evaluate the accelerator’s influence on a limited number of particles. Thus, any approach computing single particle motion can only evaluate the action of the accelerator in a subset of phase space of measure zero. Completely rigorous predictions of accelerator performance can therefore only be achieved when the accelerator’s action can be evaluated on the entire relevant phase space region. The presented work deals with both problems, rigorous predictions and accurate computation of charged particle motion in optical systems. The predictions concern the survival time of particles in storage rings. This dissertation is naturally divided into three interrelated parts. These are a detailed chapter on normal form theory; the description, application, and rigorous evaluation of what we call the pseudo invariants estimation (PIE) method for survival time bounds; and finally a chapter on the DA based tool of symplectic scaling. Each of these parts has a detailed introduction including references and comments about the history and purpose of the approaches which will be developed. Please refer to the pages (5), (50), and (126) for these specific introductions. 1 All parts of this work use transfer maps NI of dynamical systems that map initial phase space coordinates 2',- into coordinates 2'; which are considered final in some sense. This transfer map can also depend on 72,, parameters 6;, i E {1,...,np} of the physical system involved, leaving :3"; = M(§},6.). In the applications which will be presented, the dynamical system is a particle moving in an accelerator or storage ring. Following the belief that the quality of a Ph.D. thesis increases with a decrease of the presentation of well known material, the so—called DA method of manipulating and computing all derivatives of a transfer map up to order n, will not be intro- duced in detail but will only be described shortly. It is based on the possibility of introducing “=,,” as an equivalence relation, equating functions which have the same partial derivatives at the origin up to order n. Taking derivatives at the origin is not a restriction, since the origin can be shifted to any reference point by transformation. Traditionally, equivalence in this sense is expressed by considering functions as equal up to order 12, since the partial derivatives specify the Taylor expansion up to order n. The theory of constructing a differential algebra on the set of equivalence classes ” is introduced with rigor in [Ber92a]. When it is relevant to the created by “=n presented work, it will be used and restated that agreement up to order n can be expressed by an equivalence relation. Often we will refer to Taylor maps, which rep- resent an equivalence class. The terminology of Taylor maps rather than equivalence classes is used, since many readers will be acquainted with manipulating Taylor maps up to order n. It should, however, be realized that no approximation is made when Taylor maps are used. To the contrary, all partial derivatives which are equated with “2,,” are exact. Therefore, computations within the differential algebra of derivative classes leads to all partial derivatives up to order n without any approximation. In normal form theory, which is treated in chapter (1), a map [If will be trans- formed by means of a. nonlinear transformation A into a normal form map N =n A0117 0 A“. The Taylor expansion of IV has as many vanishing Taylor coefficients as possible. As usual, the symbol “0” describes the composition of maps. Although A is a Taylor map and all arguments are performed with respect to “2.3,”, strict state- ments about the partial derivatives, which are the Taylor coefficients, can be drawn. Normal form theory will serve to obtain functions f which are invariants of the map A? up to order n + l and thus f 0 AZ =n+1 f. For weakly nonlinear dynamics, these functions are approximate invariants, or pseudo invariants, of the transfer map. This is, for instance, the case for particles moving close to the central periodic orbit in a storage ring, which is represented by the origin of phase space. For the PIE method, a pseudo invariant f will be found which is suitably chosen to let f (5') describe the distance between a particle with coordinates 2’ and the central orbit. At this orbit, f(5') = 0. A = {E'If(5) S 6} is a volume in phase space which contains the origin. In an application of the map [Ff on 5, the so described distance from the origin changes by df(§') = (f 0 MM; — f (2') Where evaluating a function f at 2' is written as f '3, which is equivalent to f (2') If the maximum 6 of d; on the phase space volume A can be found, it can rigorously be said how far a particle can move away from the closed orbit in one turn. It can also be stated rigorously that particles in the phase space volume 0 with 0 = {5:1 f (5) S 6 -— N6} will not leave the volume A for N applications of the map. Analyzing the prospects of this approach and bounding the maximum of d; rigorously by an extension of interval arithmetic is the main subject of the chapters (2), (3), and (4). First results of this method were reported in [B H94b]. This approach is only rigorous for single particle motion, since the one turn map is computed for a single particle. Particle loss due to residual gas scattering, beam- strahlung in the collision region, space charge effects, and other collective effects are not considered. A substantial effort for making predictions of accelerator performance rigorous, like the introduction of interval chains in section (4.2) or the development of differen- tial algebra with remainder intervals (RDA) in [BH94a], is only useful if the transfer maps of accelerators can be described accurately. In chapter (5) it will be shown that often the influences of fringe fields are very relevant, especially for higher order deriva- tives of the transfer map. The well established DA method of obtaining Taylor maps of particle optical elements by evaluating Lie derivatives can only compute transfer maps of main field regions, and numerical integration in DA is very time consuming. For the method of symplectic scaling (SYSCA), it will be analyzed how known scaling properties of maps in geometric coordinates can be used to scale symplectic repre- sentations of the transfer map of a Hamiltonian system. The representations used are various generating functions and exponential Lie operators. To apply SYSCA, a symplectic representation is computed by integration in DA as a function of the magnetic field of one specific particle optical element and a specific particle. The Taylor expansion of this function is stored to a file. Scaling this symplectic repre- sentation of the fringe-field map for a different particle or for a different size of the optical element is much faster than the conventional methods which leads to accurate fringe—field maps and it is also much more accurate than a number of conventional approximations which will be mentioned. Additionally, this method leads to maps which are symplectic up to order 72. Several examples will demonstrate the usefulness of this method. Chapter 1 Normal Form Theory The idea of normal form transformations for simplifying differential equations is very old (according to [Nay93] it goes back to Euler’s days). After Poincaré’s work on maps in nonlinear dynamics [Poi99], normal form transformations were used for maps by Birkhoff [Bir27] and many others after him. Reference [Tur91] gives a detailed account of the history of normal forms. For symplectic maps, the possibility of normal form transformations was mentioned in [DF76] and computationally performed for Taylor maps in high orders in [FBI89]. In these approaches the symplecticity was guaranteed by the use of exponential Poisson—bracket operators. In [Ber93b] a method of computing the normal form of maps was demonstrated which can be applied. to non-symplectic maps as well. This approach can be evaluated in a straightforward way with DA—based programs and was implemented in the code COSY INFINITY [Ber92b]. This chapter will explain this method, and it will be shown that it can be used to obtain pseudo invariants for symplectic maps. We will use the symbol “=n” to indicate that the left and right hand side of an equation have the same partial derivatives at the origin up to order n; in other words, the Taylor expansions of the left and the right hand side agree up to order n. It can be shown that “=n” is an equivalence relation and therefore establishes equivalence classes. A differential algebra of these equivalence classes can be constructed [Ber92a] 5 which can be used to formulate relationships for the partial derivatives of functions. This theory establishes in a rigorous way that the well known method of evaluating a function of a Taylor expansion leads to the correct Taylor expansion of the resulting function and therefore to correct partial derivatives. All arguments will be performed with Taylor maps, but it is worthwhile to keep in mind that an nth order Taylor map, i.e. a map whose vector components are polynomials of order n, is one representative of an equivalence class established by “ =n”. What is found for Taylor maps is also true for partial derivatives. 1.1 Normal Form Theory for Order 72. Symplectic Maps Conventions: Some notations will be used often, therefore we want to define them now. J24 is the 2d x 261 matrix with 1 fort oddandj=i+1 J2d,ij= —1 forieven andi=j+1 . (1.1) 0 else The identity function will be symbolized by 5', i.e. 21(5) = x; for all if. Often an incorrect notation is used by writing expressions like (f-Tfiyi', while 2' is understood to be an element of R“ and f is a function R24 —r 1%“. This must be wrong, since 2' E R“ can not be differentiated. Sometimes this problem is corrected by using the identity f to write ”‘75)”... Since this notation can become very peculiar and unconventional, a slightly different but very useful notation is advocated. The clarity of the incorrect notation can be conserved by simply choosing 2' to be the identity map. This convention will be used consequently. Elements of B“ will usually be denoted by 5. At this point attention should also be drawn to another incorrect use of notation. For a function g : B“ -+ B, 9(f) is incorrectly written, since 9 acts on elements of It“, not on functions like f“. The correct notation for a concatenation of functions is g o f. This clear notation can become important to avoid misunderstandings. An example is the ambiguity of g(f+ 5), which could either mean 9 o (f-F 5') or g - (f+ 5). Strictly speaking it is also incorrect to write 9 o :i’ for f E It“, since “0” concatenates functions, and f is not a function. In this chapter the notation 9(3) is usually used; whenever this would cause confusion, 9],? is written for the evaluation of the function g at 53'. Now enough about these trivialities, which are only mentioned to advocate correct use of notations. The partial derivative of f with respect to its i‘h variable is written as 0,- f. If f only depends on one variable, a f will denote the derivative. Multidimensional maps will be denoted by capital vectors, e.g. M, and their 2d X 2d Jacobi matrix by the capital letter, e.g. M; so Mg,- = a,M.-. The Gauss bracket or integer part “[...]” will often be used. All the maps we will mention will be assumed to be sufficiently differentiable, such that all partial derivatives at the origin equated by “=n” will be defined. All the maps we will mention will be origin preserving. Since this restriction will be implied for every map, we will not mention it further. Concatenation of functions is denoted by ”0”. Before we state the theorem about normal forms, we will define some terminology and review some properties of symplectic matrices. Definition 1.1.1 (Order n Symplecticity) A 261 dimensional map M, with a Jacobi matrix M that satisfies the symplectic con- dition up to order n — 1, MJMMT =n-l sz , (1-2) is called order n symplectic. The set of all order n symplectic maps which are 0'“ —* 8 0'“ is denoted by Spffld'). Correspondingly, Spfflfl) symbolizes the set of all order n symplectic maps from B“ into B“. A Jacobian satisfying IVIJMAIT =n_1 Jgd satisfies det(M)2 =n_1 1. Therefore, an inverse of the matrix M (0) exists. All operations which lead to a matrix inverse via Kramer’s rule can be performed with the Taylor expansions of the elements of 1M. Neglecting all orders higher than n — 1 yields a matrix M; with the property AIM] =n_1 [HIM =n_1 Z, which symbolizes the unit matrix. The unit matrix is denoted by Z, since it is the Jacobian of the identity map 2'. MTJng =,,_1 —J2dM,(MJ2dMT)J2dM =,,_1 .12., (1.3) Therefore, just as there are two conditions for symplectic matrices, there are two equivalent conditions for order n symplectic maps, {MIMTJMM .—_,,_1 .12... M} = {MIMJMMT =.._, MM} , (1.4) where the maps are either B“ —+ IR“ or 02" —> 0“. Theorem 1.1.2 For F being either It ord’, the concatenation of/f E SP3:(F) and B E S’Pfflf’) is order min(m,n) symplectic. Proof: Let us first establish the Jacobian C of the concatenation C = A o B, Cij = (9,A,-(B) = ZBsz{(81A.-)o B} or C = A(B)B. (1.5) 121 Now the theorem follows immediately, because CJgdCT = A(B)BJ2dBTAT(B) =,, (ATJgdA) o [3’ =,,,,,,(.,,,,., .12.1 . (1.6) Definition 1.1.3 (Non—Degenerate Maps) The eigenvalues of the Jacobian at the origin M(0) ofa map [Ff are called the linear eigenvalues of M. We call a map non—degenerate if all its linear eigenvalues have multiplicity one. Theorem 1.1.4 IfA is an eigenvalue ofa symplectic matrix.~ N, then also A” is an eigenvalue of N and they both have the same multiplicity. Proof: The symplectic condition N szN T = J24 yields N"1 = szNTszf, which shows that N ‘1 and N T are similar matrices, since J27} = —J2d. Therefore, N and N ‘1 have the same characteristic polynomial; this proves theorem (1.1.4). This also shows that non—degenerate symplectic matrices cannot have 1 as an eigenvalue. We will always arrange the eigenvalues of symplectic matrices to get A2,_1/\2,- = 1. Definition 1.1.5 (Resonance of Order m) The 2d eigenvalues /\( of a non-degenerate symplectic matrix, which are ordered to give A2;_1)\2.- = 1, are said to be in resonance of order m if there is a set ofd integers k.- with 23;. lkil = m +1 and Iii-L. A’s: = 1. Definition 1.1.6 (Order n Inverse) A map ALL” with A‘l'" o A = A0 A‘l’” =n 2' is called an order n inverse of A. -0 Theorem 1.1.7 Let order n symplectic map M be origin preserving, i.e. 112(0) 2 0. Then M has an order n symplectic order n inverse. Proof: First it will be established that every map with a non—singular Jacobian at the origin has an order n inverse. We use the linear map M1 =1 M, which has an inverse, and we write AI = M 0 Nil-1 - iwhich has (5N1)I5 = 0 for all its components. Now, M-lr z, [{1,-1 0 (5+ Zenith) , (1.7) i=1 10 where raising a map to power i means concatenating it i times with itself. Inserting shows that M ‘1'” is indeed an order n inverse of NI: M-lvn o M =,, til-1 0 (2+ (—1)n1\7"+1)M, =,, 2' (1.8) n MOM-1v“ =,, (3+N)o(3+(—1)"21V‘)=,,Z. (1.9) i=1 Since the Jacobian of NI at the origin is a symplectic map, its determinant is plus one [MH91], and there is an order n inverse of NI. Since the unit matrix Z is the Jacobian of the identity 2', it can be written Z =n..1 (M‘l'n o ADM and therefore {(M“"‘ o M)MJ2dMT(M‘1'"T o in} o M-l'" =,,._, M-lqu-MT =,,_l J... , (1.10) which finally proves order n symplecticity of the order n inverse. It is worthwhile to remark that the equivalence classes of origin preserving order 7, n symplectic maps form a group, if “2,, is introduced as an equivalence relation. With these definitions and with the aforementioned conventions, we will set out to prove the main theorem about symplectic normal forms. This theorem deals with general order n symplectic maps from B2" into R“ and leads to a complex normal form map; therefore it is not very intuitive. For the special case of maps Ad with linear eigenvalues of modulus one, theorem (1.5.4) on real normal forms will be shown later, which has a graphic interpretation. In this case the real normal form map produces rotations in d two dimensional subspaces up to order n. The amplitudes of these rotations will be seen to be invariants of M up to order n+1. However, such invariants can also be found for the generalicase. The general theorem will be formulated first and specific cases will be considered later. Theorem 1.1.8 (Symplectic Normal Forms) A non—degenerate map [If E SpidMZ) with linear eigenvalues which are not in res- onance to any order m S n can be transformed by means ofa map A and an order 11 n inverse A‘l'" E S'Pffld') to N =n A0 NI 0 ALL" 6 5133:1(0') such that there are 2d polynomials f1, 1 E {1, . . . ,2d}, of order [(n — 1)/2] which satisfy the conditions f2i—1f2i =[(n-1)/2] 1 and f2i8if'2j =[(n—3)/2] f2jajf2i VA]. E {la--'ad} (L11) and describe the 2d components of IV by IV] :71 21(1) 0 El.) With a; = 2211.422; Vi E {1,“ .,d} . (1.12) Maps N and A with the given properties are called a normal form map of [VI and a normal form transformation for NI respectively. In special cases, further restrictions are imposed on f1, which will be the subject of another theorem. We will prove that in first order such a transformation exists and then proceed by induction to higher orders. This proof is chosen because it illustrates a method with which a normal form map and a normal form transformation can be computed. 1.2 First Order Transformations The map NI, which has no constant part MRI), will now be transformed to have a diagonal linear part. Let M (0) denote the Jacobi matrix at the origin. This matrix can be diagonalized, since it is non-degenerate by the assumption of the theorem. The linear transformation that diagonalizes M (0) will be denoted by A1. In order to illustrate that A] can always be chosen symplectic, we will review some properties of symplectic matrices. A discussion of such properties is, for example, given in [MH91]. Theorem 1.2.1 The eigenvectors if; and 27;, corresponding to eigenvalues A; and M, of a symplectic matrix with AM}, at 1 are .12.; orthogonal, i.e. finds. = 0. Proof: (AMI: — llfiTszfi'k = 176(NTJ24N - Juli-1'1: = 0 (1-13) 12 Theorem 1.2.2 If a non-degenerate symplectic matrix B has only eigenvalues of multiplicity one, then there is a symplectic transformation C, possibly complex, to diag()\1, . . . , A“) which satisfies A2,_1)\2, = 1. Proof: Since the matrix B is non—degenerate, the eigenvectors 6', are linearly in- dependent. Because det(J2d) 31$ 0, the vectors 112,27, are also linearly independent. Therefore, with theorem (1.2.1), 27,"szka 76 0 when AM), 2 1. Let us order the eigen- values to have A2,_1A2, = 1 and scale the eigenvectors to get 13’,1}_,szer, = 1; then the matrix C"1 with column vectors if, is symplectic. According to the proof of theorem (1.1.4) its inverse C is also symplectic. Furthermore CBC"1 = diag(/\1, . . . , AM) with Azi—ibzi = 1- Let us choose [II-1 = 2122117121 with the eigenvectors if, of 113(0) according to theorem (1.2.2). After the transformation A1, the Jacobi matrix of the nth order Taylor map N1 with N1 =,, A, 0 NI 0 Afl has the property N,,,,.(6) = my. Wk 6 {1, . . . .24} and A,.-_, 4,. = 1 . (1.14) With this choice, A1, Afl, N1 6 S’Pffld'); N has the structure of the theorem about symplectic normal forms to first order, i.e. N, =1 21(f, o it) with f; = A]. I 1.3 Nonlinear Transformations We will break down the rest of the proof, which addresses the nonlinear structure, in three parts. 1. Prove that there are transformations A and A‘l'" which transform [If to normal form structure N1 = 21(f1 o ('1') for some polynomial f1. This part was already proven in the literature [Ber93b]. For completeness and since it does not take much space, it will be proven here as well. 13 2. Show that requiring order n symplecticity for N is equivalent to the conditions for f, which are imposed by the theorem on symplectic normal forms. 3. Prove that A and ALL" can be chosen order n symplectic. After having performed the first order transformation, we proceed by applying maps A2, . . . , An which successively transform NI from second order up to order n, Nm_1 =,, Am_1o...oAloA/IoA1-lo...oA,',',l_"f. (1.15) Every map A, is chosen so as not to have an influence on any order below i. Since we are only interested in partial derivatives up to order n, we choose all maps Nm, Am, and A?” to be polynomial maps of order n. We assume that the components of Nm_1 already have the structure Nm_1,) =m_1 31(f10c'i) to order m— 1. The map Nm_1 is then said to be in normal form up to order m — 1. We now look for a transformation Am which only affects the orders greater or equal to m and transforms the map into normal form to order m. The linear part of the map All, which also is the linear part of the maps N.- with i E {1, . . . ,m — l}, is denoted by Z. According to section (1.2), the Jacobian L is a diagonal matrix. 2", Nm_1+TmoZ—L0Tm. (1.16) Tm is a polynomial map which only has m‘h order contributions. This map Tm has to be chosen suitably in order to eliminate as many Taylor terms of NW4 as possible. To find the proper choice for Tm, we write it as a Taylor expansion. For simplification 14 of notation, we use the norm ”fill = XE, lkll and denote the l” component of Tm by “Ml:m TM: 2 (Tm,,|i€)zfl .....zggd. (1.17) d I: With this notation, it then follows that the l” component of the commutator of I: and Tm has the Taylor expansion Llon—TmJoI: ukn=m = E (Tm,,|k){1,zf1-... 2;;- 0.121),“-...-()\2d22d)k“} I? = ”in T,,,,|k)( A,— 1’," .3302 {=1-...-z§;d. (1.18) 1; For simplicity we will employ the notation Al‘: H,“ A,’ . Most Taylor coefficients in .9 Nm_1 are cancelled when (TWA/.510, — Xi) = (Nm_,,,|i£). (1.19) 01(6) = (A) —— F) is called the resonance denominator. The terms that can be eliminated and the proper choice of Tm depends on this denominator. We will show that with the choice (TmJll-C.) :{ (Nm—l,ll:)/Dl(k) :: 3% (1.20) for all orders up to n, N has normal form structure. Later we will prove that this definition of Tm allows Am and 5;,1'" to be chosen 6 877?. It is worthwhile to note that the choice of 0 in the case of D106) = 0 is arbitrary. Other choices do not influence the m‘h order due to cancellation in the commutator. However, these terms do influence the order m + 2 and can further be chosen to cancel higher orders of the map. Since this process, called minimal normal form [MW92, MW93], does not conserve the order n symplecticity of A, it will not be analyzed here. 15 All terms in IVm-1 for which the resonance denominator does not vanish can be eliminated. Two cases can be distinguished in which the denominator vanishes. 1. D,,_,(E) = o if 12H — 12, = 5., and 0,.(1'5) = 011k,,_, — 1,,- = —5,,, because A2,-1A2, : 1. The corresponding Taylor coefficient cannot be eliminated, except by deliberately choosing maps [II with coefficients (Nm_1,le-i) that vanish for these Ii. 2. Even if the first case does not apply, 01(6) can vanish. For this to happen, 21:, A;}_, must have nontrivial solutions for integers k,- with 221:, [16,-] S m + 1. Such a case is called resonant and is excluded by the conditions in the symplectic normal form theorem. All coefficients in the polynomial Nm_1,2,-_1, except those with k2,-1 — k2, = 6.5, will be eliminated by the transformation Am. With an equivalent consideration for Nm_1,2,-, this means IIkIISm Nm.2i—1 =m Z (Nm.2i—llk)22i—l(zlz2)kl---(z2d—lz2d)k2d_l 3 (1-21) I? IIkIISm _, Nm,2i =m Z (Nm.2i lk)z2i (2122),"---(Z2d—122d)k2d‘1- (1-22) ii To order m, the map N", has normal form structure Nmy =m z,(f1 0 ii), which is the structure we strive to obtain to order n. Performing these steps up to order n will lead to a map N = N” with normal form structure up to the evaluation order n. Now we come to the second task of establishing the conditions on the functions f1. Theorem 1.3.1 Let the 2d components ofN have the form 1V1 =1, 2((f105)with a,- = 22;..122.’ V2 6 {1,. . . ,d} . (1.23) 16 Then the conditions f21—1f2.‘ =[(n—1)/2] 1 and fZiaifZJ' =[(n-3)/2] f2jajf2i Vid 5 {1,...,d} (1.24) are necessary and sufiicient for IV to be order n symplectic. Proof: We use the functions g) = f, 0 ii and write the Jacobian N of N as N11 Nln N: g g (1.25) N,, N..., with the 2 X 2 matrices 621-1(Z2i—lg2i—1) 321(22,_1g2,~_1) ) N,- = . . I ( 821-1(22, 921' ) 321(22, 92g ) (1 26) With this notation, we can write the symplectic condition as NJngT =.,__l .12, 4:» (1.27) iNngNJ-T,‘ =n_1 J26,j Vi,j E {1,...,n} . (1.28) (=1 For every combination i, j, this 2 x 2 matrix equation yields component wise four equations. Altogether there are 4n2 equations. We will denote the relation in the 611‘" row and the )6“ column by < 0,6 >. Case <1,1>: ELI 621—1(22£-192i—1)aQIIZZj—IQZj—l) - 621(221-1921-1)321-1(22j—192j-1) = 1:1 221—122j—1(321—1921—1321921'4 - 321921—1321—1921—1) ' ZZj-lg2i-162ig2j-l — 22i-192j—132j92i-1 = i=1 22i—122j-1221—122l{(alf2i-1alf2j—l — alei—lalf2j—l) 0 d} + Z2i—lz2j—1{(f2i-laif2j—l — f2j-16jf21—1) 0 CI} =n-1 0 t: f2i—laif2j—l =[(n—3)/2] f2j-lajf2i—l (1-29) 17 Case < 1,2 >: [i=1 ail—1(221—1921—1)321(22192)) — 021(221—1921—1)321—1I2’21921) = 1:1 221—1221(821—192i—162192j — 32192140214923”) 221921—1821921 + 221—192102j-1921-1 + 5721—1921511 = i=1 221—1221221—122zf(31f21—131f2)‘ - alei—lalej) 0 5} 221-1221{(f21—131f21 + f2jajf2i—l) 0 5} + 9214921611 =n—1 <51) 4:) 8,(z,f2,-_1f2,-) =[(,,_1)/2] 1 and (1.30) fZi-laiij + f2jajf'2i—l =[(n-3)/2] 0 for i 76 .7 (1-31) ¢=> f2i-1f2i =[(n-1)/2] 1 a f21—131f21 + f‘Zjajf2i-l =[(n-3)/2] 0 (13?) Equation (1.31) can be derived from f2;-1f2, =[(,,_1)/2] 1 and < 1,1 >, f2i-laif2j =[(n—1)/2] f2i—1f2j(f2j-laif2j) =[(n—3)/2] f21—1f21(-f2j31f21-1) :[(n—3)/2] _f22j(f2i—laif2j—l) :IIn-BVZ] —f22j(f2j—lajf2i—1) =[(n—3)/21 -f2,-6,-f2.-_1 - (133) Case < 2,1 >: This equation of the 2 x 2 matrix can be obtained from the < l, 2 > case by exchanging i and j and multiplying by —1. Since i,j E {1,...,n} are arbitrary, this leads to conditions equivalent to equation (1.32). Case < 2,2 >:. The condition of the < 2,2 > position can be obtained from the case < 1,1 > by exchanging 2i — 1 by 2i and 2j — 1 by 2j, f2iaif2j :[(n—3)/2] f2jajf2i . (1-34) 18 The condition of < 1,1 > can be obtained by f2,_1f2,- =[(,,_1)/2] 1 and this line, f21—13if2j—1 =[(n—l)/2] fir—ifuf'zj—1 (f2131f21-1) =[(n—3)/2] fit—1f21f21—1I—f2j—1aifzj) :[(n—3)/2l ‘fii-ifij—ifziaifzj :[(n-3)/2l -f22i-1f22j—1f2jajf2i =[(n—3)/2] f2i-1f2j-l("f2i-lajf2i) =[(n—3)/2] f2j-lajf2i—l - (1.35) This proves that the conditions (1.34) and f2,-_1f2,- =[(n_1)/2] 1 are necessary and sufficient for the order n symplecticity of N. Finally to the third task of establishing the order n symplecticity of the normal form transformation. To show that the map N is order n symplectic and therefore, with theorem (1.3.1), satisfies all the claims of the theorem of symplectic normal forms, we are left with showing that the transformations Am can be chosen order n symplectic with the choice of Tm presented in equation (1.20). Definition 1.3.2 A map T for which the Jacobian of Jng is symmetric is called Hamiltonian. Let a map F be 0'“ X ”23’ 4C“ and let Ft. :0!“ 40'“ with F,-(i:’) = 13(5, t") V5 6 C“ be Hamiltonian for all t" 6 R3”. Then a differential equation for g: ”23' AC“ of the form as: Fe (37,1) (1.36) with the identity t : If; —> R3“ is called Hamiltonian. If the general solution of a Hamiltonian differential equation NI :C“ x R3 40'” with 82.14.1114 = FO(M, t) and M(§:’,0) = 5? V56 0“ exists, it is called a Hamiltonian flow. Questions of existence and uniqueness are, for instance, expressed by the theorem of Picard—Lindelbf and will not be stressed here. 19 Theorem 1.3.3 If A? :0'“ X R"; —1 (PM is a Hamiltonian flow, then the map Nit. : 4 0'“ —> C“, which is defined by MAE) = 1W(:E,t“) Vii E C“, is symplectic for all t“ E El- Proof: We introduce the 2d X 2d matrix AI :(l'zd X R; —+ C“? with Aid 2 81111,. For a given t" 6 B3", M,- :0“ 40'“? defined by M,.(a':') = M(:i:',t"') Vi" EC“ is the Jacobian of Nip. Similarly we introduce the 2d X 2d matrix F with ng = 81:17,. Then, 2d 32d+1Mij = 3j32d+1M1= ajFi°(M1t) = ZajM1{(alFi) CIA/1115)} a l=1 82d+1M = {F o (M,t)}M. (1.37) We define another matrix P = lI/IJngT; then 82.1111) 2: {F o (1171, t)}MJ2dMT.+ MszMT{FT 0(11-/l,t)} (1.38) GP+PGT with G: “(11,11). (1.39) Inserting any if E 0}“ into the first 2d components yields a linear ordinary differential equation for the elements of W : E21“ -—+ C“? with W(t") = P(.i:',t') Vt“ 6 ”2+. W(0) = sz, because [Flo = 2'. Such an ODE has an unique solution. Since F,- is Hamiltonian for all t‘ 6 E23”, szG is symmetric and therefore all components of Gsz + $2.10“ vanish. This establishes that W = J2d is the unique solution for all if E C24 and therefore P = Jgd. This finally implies the symplecticity of My, M,T.J2,1M,. = .72, Vt‘ 6 r2; . (1.40) Theorem 1.3.4 For a C°° function p : CM 60' and with the Poisson—bracket oper- ator : p: acting on a difierentiable function h via : p : h = ngJgddh, exp(: p z)? is a symplectic map. 20 Proof: If y : ESL -> 0'“ satisfies the autonomous Hamiltonian differential equation ar=—(J.15p)or (:P=5)037, 9(0) =92. (1.41) then for any differentiable function q :0“ —+ C, we have 8(qo y) = mow} =(51‘p1215q10 i=(zpzq1o 17. (1.42) An example for such a function q o 37 would be the right hand side of equation (1.41). The differential equation (1.41) can be solved by power expansion with respect to re 62;, :2 .270") = 17(0) +t'(= p : at. + %{= p = (= p: 2111.7. +... (1.43) = {exp(t‘:119=)2"}|.70 . (1.44) in the range of convergence of this sum. y : RT -> 0'“ defined by the sum can be inserted into the differential equation (1.42) to show that it is a solution. This gives y'(t‘) for all initial conditions 370 and for all t“. Therefore, exp(t : p z)? is a Hamiltonian flow and, with theorem (1.3.3), exp(: p z)? is symplectic. The following theorem will only be needed later, but it is suitable to prove it here, since the necessary tools have just been assembled. Theorem 1.3.5 Let two functions p,g :03“ —>d' be such that the Hamiltonian ODE 837 = (: p : :7?) 037 has a solution for all starting conditions 37(0) 2 370 and such that the Taylor expansion ofg o 17: ”2+ —+ 0' around 0 converges at 1 to (g o 37)]1. Then exp(: 1) =)g = 9 0(exp(= 10:)5') . (1-45) Proof: The exponential operator is defined by its power series, #2 {exp(t‘ zp :)9}|ro = 9(370) + t"(= p = gm. + 1t7{=p:(=10=g)}lgo +-~ - (1-46) 21 H With equation (1.42), we can reformulate the right hand side to *2 9(a) + t'{3(9037)}|o+ %‘ia2(9°§)}lo +... (1.47) = (9 0 Mllwom) = g 0 {exp(t'u 3? =)5}l170 (1-48) for t" = 1. Since this holds for all 370, it proves the theorem. Theorem 1.3.6 lfT is a homogeneous Taylor map of order m and 5+ T is order m symplectic, then there is an exponential Poisson—bracket operator with exp(: tm :)E=m 5+ T. (1.49) Proof: With equation (1.4), we can write the symplectic condition as (Z + T)TJ2.1(Z + T) =m-1 J21 7:» T7172.) = —J21T (1.50) This is equivalent to the statement that Jng is symmetric or that the Jacobian of .1ng is symmetric. Therefore, the potential problem 5t", 2 Jng has a solution tm which is a homogeneous polynomial of order m + 1 and exp(: t... :)2’=m 5— Jgdgtm =,., 2+ T'. (1.51) In order to guarantee that Am is in S’Pffld’), we only have to show that I? + Tm is order m symplectic and choose Am =,, exp(: tm :)2'. A‘l'" =n exp(— : tm :) will be 4 an order n Inverse of Am. We assume that all the ILA?” E S’Pffld') for i E {1,...,m — 1} and, with theorem (1.1.2), Nm_1 E 3793:1(0'). Before we use the symplectic condition to check the symplecticity of 5+ Tm, we analyze restrictions on the mth order of Nm_1, which is closely related to Tm by the equation (1.20). We separate the map Nm_1 into two parts. The first part F contains all contributions which already have normal form 22 structure; therefore, it contains all m — 1 orders and the normal form part of order m, F) = 21(f1 0 ii). Since it is order m — 1 symplectic, the functions f; satisfy the condition (1.24) according to theorem (1.3.1) for n = m - 1. The second contribution H contains the rest of order m, 1\7,,,_, =m F + R. (1.52) We write the Jacobi matrices of H and H as F and R respectively and employ the 2 X 2 matrices notation which was used in equation (1.27). The linear part of the map Nm_1 has the diagonal Jacobian L with 2M 0 ) (I...) and we write 823.4122” 32j132;_1 a2i—1F2i-1 aZjF2i—l R J ( 021-4122.- 02,122, ) an J ( 021—11531 621' F 21' ( ) The symplectic condition will again be analyzed for the four equations in the 2 X 2 matrix separately, 71 Zara-1+ 12.1)J2(Ffi+ Rfl) =,,._.l J26,,- v1.) 6 {1, . . . ,n} (1.55) l=1 ti? 2511111736 + [111172332- + Rij-IzLii =m—1 J2 - (1°56) l=l Parts of the required manipulations are equivalent to steps in the proof of theorem (1.3.1); we will therefore skip some lines and refer to the equations (1.29) and (1.32). Case <1,1>: 221—122j—1{(f2£—131f2j—1 — ij—lajfZi—I) 0 5} +A2i—la2iR2j-l - A2j—132jR21—1 m-l 0 (1.57) 4: f21—131f21—1 - f2j-lajf2i—1 =[(m—3)/2] 0 (1-58) and A2i—182iR2j-1 - A2,_162,~R2.-_1 =m—l 0 - (1-59) 23 This separation can be performed, because the polynomial on the left hand side of equation (1.59) has no terms of the structure 22,-_122,--1(h 0 ii) for any monomial h of order [(m - 1) / 2]. Such terms cannot occur, since H does not contain terms of normal form structure. However, the rest of equation (1.57) has only terms of this structure. Case <1,2 >: z2i-122j{(f2i—laif2j+ f2131f21-1) 0 5} +6ij(f2i—1f2i) 0 5 '1' 3211-1321sz + Away-13214 =m—1 511' (1-60) 4:) f2i—1f2i =[(m—1)/2] 1 1 f2i—1aif2j + f2jajf2i—1 =[(m_3)/2] 0 (1.61) and A21-1621sz + A2.70214111321—1 =m—1 0 - (1-62) This separation can be performed with reasoning corresponding to that of the < 1, 1 > case. The polynomial of the left hand side in equation (1.62) has no terms of the structure zg,_122j(h 0 ii) for any monomial h of order [(m — 1) / 2], the rest in equation (1.60), however, only have terms of this structure. Case < 2,1 >: This equation of the 2 X 2 matrix can be obtained from the < 1,2 > position by exchanging i by j and vice versa and multiplying by —1. Since i, j E {1, . . . ,n} are arbitrary, this leads to conditions identical to < 1, 2 >. Case < 2, 2 >: The condition of the < 2, 2 > position can be obtained from the position < 1, 1 > by exchanging 2i — 1 by 2i and 2j — 1 by 2j and vice versa. The conditions (1.58) and (1.61) establish that the map H, which simply consists of all terms in Alm_1 with normal form structure, is symplectic up to order m. This is a nontrivial statement, since symplecticity of F was only assumed to order m — 1. The equations concerning the rest of the map H are essential for establishing order m symplecticity for the map 2' + Tm. We use the symplectic condition, which we write 24 with the Jacobi matrix T of Tm. T is a 2n X 2n matrix and, as before, the T,, are 2 X 2 matrices, T11 T1,, We can write the symplectic condition as TJngT 2 J2.) 4:) iTflJQTJ’If 2 J26,,- Vi,j E {1,...,n}. (1.64) 1:1 Since we just performed a similar computation for Nm_1 , we can obtain the symplectic conditions by setting F = 2' and H = Tm in the equations (1.56). Case < 1,1 >: 82; T2j_1 — 32,- T2,-1 =m_1 0 Case < 1,2 >: 82, ng + 021-1T2,_1 =m_! 0 Case < 2,1 >: —82,-_1T2j_1 — 82,- T2, =m_1 0 Case < 2,2 >: 82,-4ng -— 02,-_1T2, =m_1 0 Since the third and fourth condition follow from the first and second by changing indices as pointed out earlier, and since the same symmetries with respect to indices apply to the corresponding equations for H, it is sufficient to show that the first two equations always hold by means of the above equations for H. Comparing the Taylor coefficients in the equations for T yields Case < 1,1 >: Pi = l91+ 51,21 1 (I: = k: + 61.2j , (1.65) 0 = (T2j-llmp2i ’ (T'Zi-llq-‘quj 1:* (1-66) 0 = (R21-1|17)p2.-(A2.~-1 - xii/\zjl - (RZi-II‘TlQ2j(/\2j—l — ALE/Mi) :* (1-67) 0 = (1 — X£A2i/\2j){/\2i-l(R2j-llfilp2i — A2j—1IRzi—1l‘fl‘12j} - (1-68) The last line is guaranteed by equation (1.59). 25 Case < 1,2 >: P1 = (61+ 51,21 , Q: = 161+ 51,21—1 . (1-69) 0 = (szlmpzi - (Tzi—il‘flqv-l <==> (1-70) 0 = (RumpuUm—i - Fin—1)+(Rzi_1|€D<121‘—1()‘21— XII/Vt) <=> (1.71) 0 = (1 - XfA21A21—1){A21—1(sz—1llflpzi— A2j-1(R2i—ll(flq2j} . (1-72) The last line is guaranteed by equation (1.62). This finally shows that Am can be chosen order n symplectic. Am is not unique, since one could add any polynomial of orders higher than m + l to tm; this would change the terms in Am from order m + 1 up to order n. This implies the possibility that A might not be unique. However, the question if there are normal form transfor- mations other than the one which is constructed here and how these would be related to A will not be analyzed here. Performing this procedure order by order to the truncation order 71 gives the normal form map N E 5733,“!(0'), the normal form transformation A = An 0 . . . 0 A1 6 S’Pifid’), and its order n inverse ff” 6 S’Pffld') with -o NznAoMoA"l’". (1.73) 1.4 Invariants Definition 1.4.1 (Order n Invariants) A function I is called an order n invariant ofa map M if] o M =,, I. We want to show that normal form transformations can be used to find d order n + 1 invariants, which are not related to each other by functional dependences. To formulate this clearly, we introduce the concept of independent functions. 26 Definition 1.4.2 (1 functions I,- :d'd —+d’ are called independent if there is no func- tionfzd'd’l ->d' with I]:f0(11,”.,I;_1,11+1,...,Id) for anylE {1,...,d}. It can easily be observed that the d components of the identity function 5' are independent, since 21 can have different values for the same values of 51,, k 74 I when these functions are applied to appropriate :7." E C .d. Theorem 1.4.3 Ifa map fzd'd -—+d’d is locally invertible in some volume, then the components off are independent functions. Proof: Assume there is a function f :d'd‘l —> d’ which satisfies [1=f0([1,...,I(_1,I(+1,...,Id) . (1.74) Applying the local inverse 111 to the right and the left hand side yields Z:=f0(Zia---,Zl—1,Zz+1,---,Zd) (1-75) for the domain of the inverse, which, as just observed, cannot be true for any f. Theorem 1.4.4 (Invariants of Order Symplectic Maps) For a non—degenerate map M E SPZdUR) with linear eigenvalues which are not in resonance to any order m _<_ n, d independent order n + 1 invariants of M can be found. Proof: With theorem (1.1.8), a normal form transformation can be found to order n and the d polynomials of order n + 1 given by I.- =n+1 Ark--1142.- are order n + 1 invariants, because I.- o M =n+1 (A2.-. o Nix/42.- o M) (1.76) =.+1 (Nu-1° ff)(N2.- o 21') (1.77) = A2,_1A2.-({( 13,-.1 05002,. 0 5)} o A‘) =,,..l I.- . (1.78) 27 The agreement up to order n+1 in line (1.77) follows from the fact that A2,_1 o M =n 1V2,-_1 o {If and that N2, 0 161. has no constant part, N2,- 0 if =0 0. Similarly in line (1.78), {(f2,-_1or'i)(f2,-oc'i)}ox1o =n_1 1 and A2,-1A2; =1 0. Since the polynomials 1.- =1 2,, the Jacobian of the map f at the origin has determinant one and this map is therefore invertible in a volume that contains the origin, due to the implicit function theorem. With theorem (1.4.3), one concludes the independence of the functions 1,. 1.5 Real Normal Forms For the case that the linear eigenvalues of an order n symplectic map are either real or have modulus one, it is possible to find additional properties of the functions ff and 1?. Due to these properties, it will be possible to formulate a real normal form transformation to a real normal form map. This special case is of particular importance, since all maps with 2 x 2 block diagonal structure in the J acobian at the origin fall into this category. Conventions: We will order the linear eigenvalues to obtain MA = 1 for l in {1,... ,2r} and A). E R for k in {2r + 1,.. . ,2d}. An important step towards real normal forms will be the next theorem. To express this theorem, we have to define a certain subset of 0'“, C = {@1376 072d X ”22(d-r), ij-l = iygjvj E {1: ° ° - ,27‘}} ' (179) It is useful to define the set of maps which take C into C , s = {MlM‘m e c vge e}. (1.80) Theorem 1.5.1 If each linear eigenvalue of M in the theorem for symplectic normal forms is either real or has unit modulus, then the normal form transformation 1‘1. and 28 the normal form map A7 have the property 11.:1de—+C and IVES. (1.81) Proof: The eigenvalues of modulus one satisfy A2111 = A3,» In order to make fill symplectic, all eigenvectors are scaled to 0bS€fV€.1-)g;_1.]2d772j = 1. With M (6)1};j_1 = Mpg-4 for the first 2r eigenvalues, one gets 17214 = (127%. The factor a has to be purely imaginary, since -1 a :0-1 -oT ." . _ 4T. -0 . _ 4T 7” _ -1. Scaling the vectors 272,- to give @‘Jgd-figj = i requires a = -i. Since for every eigenvector with real eigenvalue its real part is also an eigenvector, the eigenvec- tors 6,, with k > r are chosen to be real. We can write Afl = 2,23,532, and 141.1“ 0 {41.1 = XE, 171/11,; = 5'. Multiplying with 27.34sz from the left reveals Am, = rig-4.12.13:—iv’£‘J2d5=iA;,2J-_l VjE{1,...,r} and (1.83) A13]; = figk_1J2dEE B Vk E {7‘ +1,...,d} . (1.84) This is equivalent to the statement 11(5) 6 C for all real vectors :3. The inverse Iffl must take elements of C into IR“. From these properties and the fact that M takes real vectors into reals vectors, we can conclude that All 6 S. The next step is to show that all film, which are obtained from [ill by the steps of the normal form transformation, are in 8. Assume that NW4 takes elements of C into elements of C. Writing the components of Nm_1 as polynomials, this means E< -* k... r 1: Eli “‘n (Nm-1,2,-_1|Ic)y{‘1 yé" . - - 1123-3 .1133 ---y231“ (1-85) . i; < -o k I: k k _ k . 2r ’6 =1 22“” (Mani-loving!...y2221y23'marl—z) These lines entail relations for the Taylor coefficients, which can be expressed most 29 conveniently with the 2d x 2d matrix S, 1 ifl=k+1andk£2rodd 1 ifk=l+1andk_<_2reven 5’“ ‘ 1 iszl andk>2r (1'86) 0 else Equation (1.85) reveals (N..-1,2.--.IE) = «Mart-1513‘ - (4)241“ . (1.87) Since this condition does not relate coefficients of different orders, an n‘h order Taylor map which is in S has this property for any truncation order smaller than n. According to equation ( 1.20), Tm has this property also, since for the resonance denominator a similar equation holds, namely Din--10?) = (A2.-.—X‘)=(A2.—XS‘)'=02.(SE)'. (1.88) To finally show that Am =n exp(: tm :)5 has the same property, we prove a helpful lemma for which we need another 2d x 2d matrix ifl=lc+1andlc$2rodd iflc=l+1andlc_<_2reven ifk=l andlc>2r else Ski = (1.89) OHN.N. and the corresponding map 3' with §(g‘) = S37. Lemma 1.5.2 m“ e s and F’ e s, then :FTé‘F’ e s. Proof: It will be useful to remember that, for a map 112 and its Jacobian M, 5{goM} =MT{(5g)oA/T}. (1.90) 30 The function which takes all 37 E C into 15(3))“ is denoted by 13" ‘. Due to membership in S, T and 1:". are functions on C with f 2 ST" and 13" = Sfi‘. Accordingly, 13" o 3 takes 37" into 1.7. (37)" for all j E C. From this information, the subsequent relations follow: (5% = (15F: o silt-r = («516' o 3}o my)" (1.91) = ({3T5F‘Hy)‘ = ({STé'S-lfingr . (1.92) Now it can be concluded that -o TT(.27)(5 ”)Ig = T'T*(v)S‘T{(ST53-‘F)Ig}' (1.93) = 3-“{(T'T§T*§T5F)Ig}'=é{(“T5“)Ig}*, (1.94) which proves the lemma. To analyze the exponential Poisson—bracket operator, we observe :tm :1? = 87th2d51-7" = T1574. (1.95) All terms in the expansion of the exponential map have this structure, exp(: 1,.. )5: 2+ T1,. + game... +31—!(r,:5)2rm + imam + . .. (1.96) and, due to lemma (1.5.2), with Tm, also Am is an element of 8. Also the order n inverse 4:” -—--1, exp(— : tm z)? is in 5. Since the concatenation of two elements of S is again in 8, also 1V", =n Am 0 lilm_1 o A?“ is in S. In this last step it is again crucial that truncating Taylor maps preserves membership in 5. Performing the normal form transformation to the evaluation order n yields therefore [\l E S, which proves the first part of theorem (1.5.1). The performed transformation can be written as -o -o -o AznAno...oAgoA1. (1.97) 31 A1 takes real vectors into C and the concatenation of all nonlinear maps if", is in S. This proves the second part of the theorem (1.5.1), namely that A takes elements of B“ into C. Theorem (1.5.1) has strong implications for the functions f, of the normal form map in theorem (1.1.8). Theorem 1.5.3 Let A? be a map satisfying the conditions of theorem (1.1.8) on symplectic normal forms. If each of the r complex linear eigenvalue offld has unit modulus, then there are r functions (bi : H’ x ”2“" ——> B with f2j-1 :[(n—1)/2] exp(i¢j) and (1-98) fgj =[(n_1)/2] exp(—i¢j) for]. E {1,. . .,7‘} , (1.99) fk : I’de-r—ilt forkE {2r+1,...,2d}. (1.100) The phases 451' satisfy the condition 81¢,- =[(n_3)/2] 8,45]; Proof: From theorem (1.5.1) and with N1 = 21(fgoc‘i) where (lj = 2211122” it follows that {22j-1(f23'-1° Ellis? =i({223'(f21'° alllili V376 C- (1-101) Realizing that 6(37) 6 H' X lid” for 37 E C, one concludes 1314(5) = {M5}: vt’e r x IR“ . (1.102) In order to take advantage of the requirement f2j_1f2j =[(n_1)/2] 1, we write f2j..1 = exp(itbj) + Pj and f2,- = exp(-i¢j) + Pf, (exp(icbj) + Pj)(exp(-i¢j) + Pf) =[(n—1)/21 1 - 0103) The constant part of this equation shows that Pj has no constant part. For comparing higher orders, we write Pj exp(—i¢j) + P; exp(iqu) + 131in =[(n_1)/2] 0 (1.104) 32 and assume that R, has no contributions up to orders m — 1. Since the constant part of f1 equals )1), we get up to order m that Pjh; + QJ'AJ' =m 0 . (1.105) Therefore, the mth order of Pi must have the structure Milt,“ with Rm.) : H' x Rd" —+ B. Writing the non—constant part of (15]: as $5") shows that f,,_, zm A,{exp(z‘¢§"’) + mm} =m A,- exp(igbg-n) + 112...,» . (1.106) The last equation holds, since Rm) has no contributions to orders lower than m and therefore exp(Rm,,-) =m 1+Rm,j. ng can be chosen to be 0 by taking its contribution into of”. Performing this argument up to the evaluation order n yields that up to order n, f2j_1 has the structure of exp(idj) for some function ¢j : H' x Rd" -+ ”2. According to the theorem on symplectic normal forms (1.1.8), f2z'aif2j :[(n-3)/2] f2j6jf2i ViJ E {1,-Md} - (1-107) Writing the f, as given by equation (1.99) yields 6,625]- =[(n_3)/2] 81¢“ which proves the theorem. This finally puts us into the position of formulating a theorem on real normal form transformations. Theorem 1.5.4 (Real Symplectic Normal Forms) A non—degenerate map M E S’Pfflft) with linear eigenvalues which are not in reso- nance to any order m S n and are either real or have modulus one can be transformed by a transformation 1? E S’Pidwt) to If 2,, Bo [Flo é'l’" E S’Pfflfl) where all these maps are I?“ —+ B“ and R2j_1 =n 22j_1COS{’(/Jj ob} + zgj sin{wj o b} , (1.108) 33 jo =n —22,-_1 sin{t,l’j o b} + 22]“ cos{wj o b} , (1.109) R2k—1 =n Z2k-lgk 0 5, (1-110) R2,. =,, 22k/ 9,. o 13', In.) = 1, A2,. 6 H. (1.111) The function bis given by by = 231-4 +231 for IAQJ-l = 1 and bk = 221,422), for A2), 6 B. The functions 11),- and gk are Rd —+ H2. Proof: We will again choose n‘h order polynomial maps for If, If, and 3"”. A linear transformation 6 and its inverse are needed to obtain these maps from the normal form theorem (1.1.8), 1 —i7r 1 in 021-1 = 7512 ”(2214+ 221‘) , 023' = 728 /2(—sz—1+ 221‘), - 1 i'l' ' — 1 in - 0231.1 = 7.26 ,/2 (Z2j—1 + 222]”) a 2f = $6 /2( ZZj—l - 2221') ) C'k=z;c and C',:l=z;c for j€{1,...,r}, k>r. (1.112) The easiest way to see that C" "1 is the inverse of C: is to multiply the Jacobi matrices. Since these matrices have block diagonal structure, the result is evident. All 2 x 2 submatrices on the diagonal have determinant one and therefore are symplectic, which _9 .9 _. establishes the symplecticity of C". Let us now prove the claim it = C o N o C“. First evaluate E = ('1' 0 0‘1, bj and Ck = 22111-1221: = bk . (1.113) 2 i Ci = “(Zij-1 + Zij) = i 2 Let us write a scaling transformation 5' by 31- = 52,- forj S r and 3;, = 2k forj > r (1.114) and introduce the new phases «[2,- = (bk 0 5'. The functions b, o C: o A are order n + 1 invariants of the map M, since they are related to the I; = a; o A of theorem (1.4.4) 34 by 5‘ o b = [i o (‘1. Performing the transformation C.7 on if leaves 122,-.l = 021—1 oN'oé-l (1.115) 1 . _. .. ~- = 7.26_”/2{32)-1(f2j—10 a) + 221'”sz all 0 C 1 (L116) e"”/2{(z..-. + 2..) sow. o a) | 3 3'... [\D + i(zgj_1 —22,-)sin((15j06)}06—l (1.117) .9 .9 2,, 22,--1 cos(wj o b) — 22,- sin(wj o b) . (1.118) Similarly we obtain R2, = 02,0N'oé-1 (1.119) 1 . =n Ee‘”/2{(—22j-1 + 221') cos(¢j o ('1') — l(22j_1+ 22]) sin(¢j o 5)} o 5‘1 (1.120) =n 22j_1 sin('z,l)j o b) + 22, cos(1,bj o b) . (1.121) For real eigenvalues A2). one gets Rzk—i =n Z2k-1(gk 0 3) , 32}: =7; sz/(gk 0 5) With 9k = f2k—l 0 57- (L122) Finally we have to establish that fit?) = (éoA)|g is a real vector. If takes 53' E B“ into C and Cl takes elements of C into R“, which establishes that B is R“ -+ R“. The phases 11),- are real, since 5(5) 6 I' x Rd" and is taken into R by <15 j, according to theorem (1.5.3); similarly gk(.i:') = f2k_1 0 5(3) is real. Interpretation: With the help of a 2d x 2d matrix which depends only on b, we can formulate the real normal form map as R2j—l __ cosh/2 o b) — sin(1/) o b) sz-l ( Rzi ) —" ( sin(1,bob) cos(wob) ) ( zgj ) v (1.123) 35 for complex /\2j. For real eigenvalues /\2k, a diagonal matrix can be used, R2k-l 9k 03 0 ZQk—l . :n ‘* . .1./ (11.. l (0 1/(gkob))(z2k l (1 2“) The motion generated by E is described by rotations on a circle for every [)1 J-l = 1 and by motion on a hyperbola for every A). E B. The speed of this motion only depends on the b1, the order n + 1 invariants of the mapff. These invariants are the square amplitudes z§i_l + 2%,- of the rotations and half the square minimum distances from the origin 22;...1z2k of the hyperbolas. 1.6 Action—Angle Variables Definition 1.6.1 (Order n Action—Angle Variables) Given a Taylor map 1?; we call maps faction variables and 6; angle variable oflil if 1. the map (To?) is order n symplectic, 2. the functions J) are order n + 1 invariants of N, 3. up to order n — 1, the functions 01' o M — a,- can be expressed as functions of the j only. The rotations in real symplectic normal form space remind one strongly of action— angle variables for periodic motion in classical mechanics. We will prove that they are indeed order n action—angle variables. Theorem 1.6.2 Action-angle variables of normal form maps For a normal form map Al as given in theorem (1.1.8), there are ordern action—angle variables j and 62. 36 Following, we introduce variables and show that they satisfy the properties of order n action-angle variables for the normal form map If. Let 1 J1: 321—1221 7 01' = 5108(221/221—1) , (1.125) where we restrict the domain of the action—angle variables to {my 6 C, y; 74 0}. The Jacobian of (£62) is denoted by K and .has 2 x 2 block diagonal form. Like before, we represent these blocks by K11 with 19,1 6 {1, . . . ,d}, i.e. K11: ( _ZL 2:1) . (1.126) 2221-1 2221: Since the determinant of each 2 X 2 block is one, the matrix K is symplectic. The J,- are order n + 1 invariants of the map 117. Let us analyze the change of the functions 01' under action of N .. 1 01° N - 09' = §{log(Nn,2,-_1/z 71.2)) - 10g(221—1/221)} (1-127) = $110.11..-.)—log(f..>}oa=.-.«w (1.128) It has to be noted that the domain of the functions (11 o [\7 — oi does not have to be restricted. This shows that all three conditions to call (j, 62) order n action—angle variables of IV are satisfied. In spite of this theorem, it becomes clear that the usual concept of action angle variables is not suited for discussions in the context of Taylor expansions. This was seen in some intermediate steps when the origin had to be excluded and rises from the fact that the polar angle cannot be defined as a differentiable function at the origin. It is therefore much more suitable to stay in the real normal form space and to consider rotations and hyperbolas as the most basic concept of motion. Definition 1.6.3 (Normal Modes) [fa map R2"! —+ It“ can be transformed by symplectic transformations to rotations 37 and hyperbolic motion in d decoupled spaces, then these decoupled motions are called the normal modes of the map. The square of the radius of rotation and half the square of the distance of the hyperbolas from the origin are called normal invariants of the map. It is worth mentioning that every map which— can be transformed to action—angle variables also has normal modes. The transformation between these two notions of elementary motion can be taken from the proof of theorem (1.6.2). By this definition, a symplectic map in its normal modes has the structure of If with 321-1 _ cos(0,-) —sin(0,~) 221.4 ( sz ) _ ( sin(0j) cos(0,-) ) ( 22]. ) v (1.129) for the r rotations and with a diagonal matrix for the d — r cases of motion on a R2 —1 h 0 22 —1 (R2: )2” ( 0k 1/h,, ) ( 2,: ) - (1.130) The functions 0,- and h;c are It“ -) R, and the d components of bwith b, = z§j_1 +sz hyperbola and bk 2 221,422,. are invariants of R. Given a symplectic map M that can be Taylor expanded at the origin; its Taylor map Mn up to order n is order n symplectic. Our previous discussion was based on the condition that the linear eigenvalues are not in resonance up to order n. In the case that the map A? has d normal modes, it is a reasonable question to ask if the resonance condition could avoid the possibility of computing the normal form map of Mn. A statement about this problem can be proved most easily with the help of the following two theorems. Theorem 1.6.4 Symplectic maps preserve the area in all of the d 22j_1 x 22,- sub- spaces. 38 Proof: First we need to establish the conservation of the Poisson brackets : f : g under symplectic maps. We can write the derivative with the Jacobian M as 50011?) = MT{(5f)oM}. :fo (kl : (g o M) {(5Tf)011/7}1WJ2d1VIT{(5g)01VT} (1.131) = (5TfJ2dgg)o11-l = (: f : g) 0 NT. (1.132) When new coordinates are defined by the map, Q = MM) with 11‘, Q E B“, then the change of an area element is given by the change of the differential 2-form dQ2j-1dQ2j dQ2j—1 dQ2j szj—1 /\ szj = dqzj“ A dq2j( (1921-1 (1921‘ — (1421' (1421—1) (1133) = 4421-1 A @2211 11426-1 1 M-zjllq‘ (L134) = dqz.-. A cum-{c <12.-. z (121') o MM.- (1.135) = dQ2j_1 /\ dqgj . (1.136) This is true for every of the d subspaces and proves the theorem. This proof and also another more intuitive proof can be found in the second chapter [SSC94]. Theorem 1.6.5 Let a symplectic map in its basic modes be written as equation {1.129) and {1.130). Then there are functions #5 : W —+ B and gk : B“ —+ B with 0]- : -wjoband hj = gj ob. Proof: First we will give a graphic proof using theorem (1.6.4), after that an al- ternative more formalistic proof will be given. Let us analyze what conditions the theorem about area preservation in every 221;; x 22,- coordinate plane imposes on normal modes. Consider first the case of rotations. An area located between two similar radii and two initial angles 451,,- and 432,.- is transformed into an area between the same radii and two angles (151,, and (152,, by one application of the map M. If area preservation is to hold, the difference between the angles is invariant. Therefore, the 39 phase advance can only depend on the invariants of the map. A similar argument is possible if the motion is confined to hyperbolas. If the factors 9 in the 2 x 2 matrix describing the motion, 221-1 _ g 0 sz—1 . (221 )2*(0 1/9)(221 )1’ (HM) depend only on invariants of motion, the map is area preserving. If the factors would depend on the position on the normal invariants, this condition would be violated, since not all area elements located between two hyperbolas would be conserved under application of M. As a conclusion it can be stated that the functions 0]- and h. in equation (1.129) and (1.130) only depend on the invariants of the map, which are the components of b. For the alternative more formalistic proof, we use C" and 6“ as defined in equation (1.112) to obtain the map [)7 = Cl“ 0 If 0 6 which is an element of 5. Elements of S are maps which take elements of C into elements of C as defined in equation (1.79). This means that iN2j_1(f/') = N2,(37)*, Vj E {1, . . . ,r}. With the Jacobian of 1V written as 32 “-1N2i—1 321N2i-1 . N,~ = J , 1.138 J ( 82j—1N21 3210/21 ( ) the symplectic condition will again be written in 2 x 2 matrix notation as : N2i-1 1N2j_1 I N2,’_1 I lng .._ T .._ n . I: J26U " (NJ2dN )1] - g‘NtIszvjl ( :N2i :N2j-l :N2i :sz ) a (1139) where the notation : f : g = 8-7 f .12459 for the Poisson bracket between f and g was used, which was already defined in theorem (1.3.4). The conditions of the < 1, 2 > po- sition and of the < 2, 1 > position lead to equivalent conditions due to antisymmetry of the Poisson bracket. We are therefore left with the conditions 2N2g_12N2j_1 = 0, 2N2i I sz = 0, 1N2,‘_1 I Ngj =1 6,") . (1.140) 40 From this it follows that (21V2,_1222j_122j) = (I N2,_1 2 1V2j_11'V21') (1.141) = N2,_1(:N2,_1:1V2,*)+(:N2,-_1 31V2j_1)1V-2j (1.142) and therefore 22j_102j_1lV2,'_1 — ZQJ'BleVQi—1 = (5ng23‘_1 . (1.144) If we write N1 = 2m; and insert into the second line, we get z2j—162j—ln2i—1 — 223'321‘7121—1 = 0 , (1-145) from which we conclude that there is a function g : C ——» B with h = g 0 ii where {i is again defined by a,- = 22,422,. Similar to the evaluation in the equations (1.113) to (1.122), we conclude that If = C." o [VG-L1 has the structure RZj—l _ cos(i,b,- 0 ii) —- sin(7,/),- o 1;) 221-1 ( Raj ) — ( sin(wj o b) cosh/2, 03) ) ( 22,- ) a (1.146) for rotations and R2k-l 9k 05 0 22k—1 ) (R21: ) (0 1/(gkob))(22k ( ) for motion on hyperbolas. Theorem 1.6.6 If a symplectic map 114 can be transformed to normal modes by a symplectic transformation E, and M as well as 8’ can be Taylor expanded up to order n, then the Taylor map Mn can be transformed to normal form by means ofthe method presented in the proof of theorem (1.1.8). Proof: The map M can be written with the existing transformation if and If as _. _. A zg—loRog. (1.148) 41 One can compute 11' = 5‘1 0 T1 with C“ from equation (1.112). Since with the given assumptions, all the involved maps have Taylor expansions to order n, this relation holds with respect to =m for all m E {1, . . . , n}. Therefore, --0 it =,. X-‘oNo/i' (1.149) where, due to theorem (1.6.5), [\7 has the structure of a normal form map given in theorem (1.1.8). Assigning the linear transformation [1.1 :1 f1 gives the linear transformation for the theorem (1.1.8) on order n symplectic normal forms. This transformation can be used to compute All =,, .41 o M o [fl-1. The second order expansion of A. o 41-1 yields the map 5." + T2, which eliminates all second order terms of All which do not have normal form structure. It is order 2 symplectic and can therefore be used to compute A} = exp(: t2 :)5. Similarly all subsequent orders are treated with HT... =,,, [to "';lo...oA;,‘_"{. (1.150) Continuing this strategy to order n shows that all the T... exist and that there is no problem with vanishing D100), even when there are resonances of the linear eigenval- ues. This shows that under the assumed conditions all the (Nmylli) in equation (1.20) vanish whenever 01(5) vanishes due to resonances described in part 2 after equation (1.20). Theorem (1.6.6) gives rise to a very useful corollary which can be used as the basis for computer proofs to establish that a system does not have d invariants of motion or to establish that a Hamiltonian system is not integrable. Corollary 1.6.7 (Excluding Invariants) Given a Taylor—expandable symplectic map M : R“ —-> R“. If at any order during the normal form transformation a division by a vanishing resonance denominator becomes necessary and the numerator at that point does not vanish, then the map M does not have d invariants of motion. 42 Corollary 1.6.8 (Excluding Integrability) [fa Hamiltonian equation of motion has an unique Taylor erpandable time step one map M and the normal form transfor- mation procedure would at some order require a division ofa non-vanishing coefficient by a vanishing resonance denominator, then the Hamiltonian system is not integrable. Constructing computer proofs, of course, requires completely rigorous arithmetic. Such computations can be achieved with interval arithmetic techniques [Moo88]. 1 .7 Pseudo Hamiltonians Theorem 1.7.1 For a non-degenerate map All 6 Spfflflt) with linear eigenvalues which are not in resonance to any order m S n and all have modulus one, a function H : B“ -—> R can be found with M =n exp(: H :)2'. (1.151) Proof: According to theorem (1.5.3), there are d polynomials d),- of order [(n - 3) / 2] with the property that the J acobian of (if is symmetric. This implies that the potential problem 5(1) = if has a solution ‘13 which is a polynomial of order [(n - 1) / 2]. 4 In order to prove N =n exp(—i : o ('1' :)2', the following observation will be instrumental; d :0 o Zi:{22,-_1(h o 5)} = 221-1 Z{0,,.-1( o (7)824}; o 5') — 02m 0 5)a,,._,(h o 5)} - az,(<1>k:5)(h o 5') (1.152) = 2d: 22,_122k_122k{(8k<1>01h)o 5 —(a,.<1>a.h)o 5'} - 2.1—10619) o 5}(h o 5) (1.153) =n 22j_1{(—(bjh) 05‘} , (1.154) 43 and, due to antisymmetry of the Poisson bracket, :Qoa': {z2,(ho&')} =22,{(¢,h)oa'}. (1.155) With these formulas the evaluation of the exponential operator becomes obvious; exp{—i : (I) o {i :}22,-_1 =,, exp{ iqbo {i}22,_1 , (1.156) exp{—i : (I) o if 2):»:2, =,, exp{—id> o 5}22,~ . (1.157) The original map M can therefore be represented by .4 A"1'"o(exp(-i:(1>oa':)é')0/1.=n [VT . (1.158) In the next step, the conservation of Poisson—bracket operators under order n symplectic maps is needed. Theorem 1.7.2 For a map M E SP:d(F) and two functions g,f E CS; :foM:(goM)=n_1(:f:g)oh7. (1.159) The field F is either It ofd' and the functions are either real or complex respectively. Proof: A corresponding proof for symplectic maps was given in equation (1.132). Now we merely have to check what changes by considering order n symplectic maps. With the derivative of concatenated maps written with the help of the Jacobian, we prove 5(f o M) = MT{(5'f) o M}. :foM:(goM) = {(5Tf)0M}1WJ2dWIT{(5g)01\T} (1.160) =n_1 (nygdgg) o 11.4 =n_1 (: f : g) 0171 . (1.161) With this theorem and with theorem (1.3.5), we can reformulate equation (1.158) to a single exponential Poisson—bracket operator. The notation H = —i o ("i 0 f1 yields M =n exp(: H :)2'. (1.162) 44 H is called the pseudo Hamiltonian of the map 191. It is an order n + 1 invariant of If}, since it is a function of the invariants I,- = a, 0 4 from theorem (1.4.4). This fact also becomes clear by virtue of theorem (1.3.5): Ho AT =,, H 0 (exp(: H :)E') =,, exp(: H :)H =n H. (1.163) Finally we want to show that H is R“ -) Ht. Since the functions (9', are H' x ”2"" —+ B and 8,0 = 05,-, we get (I) : H" x Rd" ——» I. Theorem (1.5.1) specifies that lf. is It“ —> C, 61' in turn takes elements of C into 11" x Rd". Taking all these facts together we find that H = —-i o [i o A. is a real valued function of R“. 1.8 Parameter Dependent Normal Forms The presented theorems also apply to maps 11:1 : RM” —+ B“ that depend on p parameters 5. For such maps, every Taylor coefficient can be viewed as a function of fixed but arbitrary parameters. All proofs are performed as before, and finally these functions of parameters in [if and A are substituted by their Taylor expansions to the evaluation order to obtain the parameter dependent Taylor maps. For computations, one starts with a Taylor map with respect to the coordinates as well as the parameters. Then one performs all computations according to the described normal form method while taking the partial derivatives with respect to the parameters into account when equating with “=n”. It ought to be mentioned that in a given order m there will be contributions of coordinates to lower orders i < m, when parameters occupy m - i orders in the corresponding term of the Taylor polynomial. It can be shown by reviewing the previous sections that in all the given proofs it was never required that a given order only contains the same order in the coordinates. Due to “2””, we were only concerned about the order of partial derivatives, which can be taken with respect of coordinates as well as parameters. 45 One exception applies. In the parameter-free case the map was assumed to be origin preserving 114(0) = 0. In the parameter dependent case 171010.) = 0 will be required. However, C((f‘) = 114(0,5‘) with a 6" 6 1R” does not have to be 0; C" is a parameter dependent constant part, which has to be eliminated before performing the normal form transformation. The procedure, described subsequently, can be found in [Ber93b]. This simplification is achieved by transforming the map to a parameter dependent fixed point. To do this, we have to extend the involved maps to be R2“? —+ R24”, _0 e.g. (If/1,6). Following, (5,6) will describe the identity function. After the required transformation 141.0, the nth order polynomial map lilo =,, A}, o (M, ) 0 (1f; 1,6) will have the property -0 ...—o -o No(0,6‘) =,. (0', *) V6" 6 12”, (1.164) whereas at the parameter dependent fixed point if“, : Ht” —* R“ of the original map, 17f has the property 23,,(6‘) =.. 12(5,,,(5’*),6’*) V6" e 112? . (1.165) The fixed point is a function of the parameters and can be expressed by the following lines, -. _. -. (6, 6*) :73 (l _ 5’6)|(2I'1(8’.)’g‘) 1 (1.166) (1 —z,6')-1vn|(6,~.) =,, (23.46“), 1.) vs‘eznp. (1.167) The required inverse exists, since the linear part of the map [Fl does not have 1 as an eigenvalue. To completely understand this notation, one has to keep in mind that 6‘ is a vector of It”, whereas (2', 6) is the identity map in RM”. The transformations 1‘10 and [161 are therefore 14.0:5—5fgrOé, K51=5+E;306, (1.168) 46 which are symplectic maps. 1.9 General Normal Form Theory If a map is not order n symplectic or if its linear eigenvalues are degenerate, one can still try to transform it in such a way that its Taylor map to a given order has a simpler form, meaning that it has less non—vanishing coefficients. Such a transformation can be obtained from the general normal form procedure which will be lined out. Definition 1.9.1 Let a map NT :(l'd —> (l'd be given and let its Jacobian at the origin be diagonalizable. The successive performance of the following n steps produces the general n‘h order normal form transformation If of M 1. Find a linear map A) such that the Jacobian of .41 o lid 0 1‘11“] is diagonal. 2. From NW4 2 Aim“ 0 . . . o 1‘11 0 171 o 14.1" o . . . o 1422?, find 1‘1", with the help of Tm, which is given analogously to equation (1.20) by (N..-1..Ik’)/D.(E) if D.( ii) at 0 0 if 0,05) = 0 ’ (”69) (ME) = { where the previously used notation for Taylor coefficients is applied. 141.", is given by 44—. Tim =m exp(TZgfi , 4;,“ =m exp(—T,,T,8)5 , (1.170) where the exponential Lie operator is defined by the power expansion of the exp function. To order m this yields -o Am =,,, 5+ 71, A"’" =,,, 2— Tm. (1.171) As shown in section (1.3), equations (1.170) and (1.171) guarantee that 1‘1", changes [\7 _1 in a suitable fashion to eliminate the highest possible number of 47 mth order Taylor coefficients in KmNm_14;1'“. Use this iteration from m = 2 tomzn. In general this method leads to an nth order polynomial map if 22,, 4,, o . . . o A} which simplifies the Taylor expansion of 11:1 by 1V 2,, 14.0 M 0 ALL". For two special cases, we will specify further properties. Definition 1.9.2 (General Resonance of Order m) The d eigenvalues of a di- agonalizable d x d matrix are said to be in general resonance of order m if there is a set ofd integers with 221:, Ila-l = m + 1 and [Ii-1:1 A!“ = 1. Following we will call a map diagonalizable when its Jacobian at the origin is diago- nalizable. Theorem 1.9.3 (Nonlinear Diagonalization) If a diagonalizable map 1131 :(l'd —> 0'“ has linear eigenvalues which are not in general resonance to any order m S n, then there is a map and an order n inverse If and ALL" such that the polynomial map of order n with If" =n 141.0 M 0 14‘1”” is given by N) = 112). Proof: Due to the condition, there are no general resonances up to order n, and therefore all resonance denominators 01(3) do not vanish and the normal form process creates a map with no Taylor coefficients between order 2 and order n. Theorem 1.9.4 For a non—degenerate map M E 5793:1(112) with linear eigenvalues which are not in resonance to any order m S n, the formalism demonstrated in the proof of theorem (1.1.8) is equivalent to the general normal form transformation, if the linear matrix in step 1. of definition (1.9.1) is chosen to be symplectic and the eigenvalues are ordered such that 12,1112,- = 1. 48 Proof: The diagonalization with a symplectic matrix is identical to the diagonal- ization performed in section (1.2). To show that also step 2, the iteration to higher orders, is identical to the process in section (1.3), we have to show that _._. exp(Tnfa)5= exp(: tm :)5 (1.172) with 5t", 2 Jngm. This, however, is obvious, since the Poisson bracket is defined as :tm : f = Hrthgdgf = T38). One can also restrict the formulated method to avoid small denominators when perturbations of order n symplectic maps are analyzed. Definition 1.9.5 Let a diagonalizable map Ad : R“ —+ B“ have linear eigenvalues which are not in general resonance to any order m S n. The successive performance of the following 71 steps produces the restricted nth order normal form transformation 1‘1 of NT. 1. Find a linear map 41 such that the Jacobian of 41 011? o 41" is diagonal and /\2,-_1 = A3, for complex eigenvalues. 2. From Nm_1 = fink] o . . . o 14.1 o M o If? 0 . . . 0 14:3, find Aim with the help of Tm, which is given analogously to equation {1.20) by If k2j_1 — kgj = 621'...” for Odd l 0 (nglk) "—' 0 If kzj — k2j_1 = 62L] for even I (1.173) (Nm_1,glk)/D((k) else and choose the maps [1,, and 14;,” as in equation (1.170). This choice of Tm corresponds to the choice for order n symplectic maps given in equation (1.20), since here (Tmlli) is chosen 0 whenever D)(l-i) is 0 for an order n symplectic map. When the transformation procedure is performed successively up to 49 order n, the only Taylor coefficients which are not eliminated are such that N) = 21(hlo it) with a,- = 22,422,. However, there are no special restrictions on the polynomials h,- of order [(n — 1) / 2], since M was not assumed to be order n symplectic. It is worthwhile to introduce a real restricted normal form. To do this, we first observe that [V E S with the definition of S and C given in equation (1.80). The formal proof is lengthy and was already performed for order n symplectic maps when theorem 1.5.1 was proved. Therefore, it will only be mentioned that every step in that proof can be performed for the non—symplectic case, when the linear map 14.1 is chosen to be R“ —> C. With the reasoning used for the order n symplectic case it . follows also that H = C 0 A7 0 6'1 is a real map and we can write H with the structure R2j—l _ -' cosh/13' o b) — sin(z,bj o 1;) 2214 ( R2,“ ) "n (r O b) ( sin('ll’j o b) 608((b, o g) ) ( 22], ) (1.174) for complex eigenvalues 12,11 and 12,-, and R2k—l " 9k 0 f; 0 22k-1 :n b .. 1.175 (12... 1 ’""° (0 1/(g,ob))(z.. ) ( l for real eigenvalues. b again is given by b, = 231-4 + 2%- and bk 2 22),-122k. Let M 6 57’3“! satisfy the conditions of theorem 1.5.4 on real symplectic normal forms and a perturbation of this map 11.4e be given. The restricted normal form procedure will obtain a normal form map 11". of AT. which is a perturbation of the order n symplectic normal form map Fl. Then the T'j are close to one. In the case of complex eigenvalues, which have modulus one for M, the motion described by N: will not be on a circle but will slowly spiral away from the a circle. Chapter 2 Long Term Estimates for Weakly Nonlinear Motion Estimating the time of stable motion for planetary systems has first started the in- terest in the stability of weakly nonlinear mechanical systems. In accelerator physics this question became important with the introduction of storage rings. In large stor- age rings particles often have to be kept in the accelerator for up to a billion turns or more. The presented work contributes to the complex subject of stability analysis by providing a method which allows one to compute rigorous lower bounds on the time of stability in weakly nonlinear motion. In the past, the question of long term stability in storage rings has been analyzed by various methods including kick tracking [Tal91], element by element tracking and one—turn map tracking [Ber88b, Yan91, KSZ92], symplectic long term generating function tracking [Ber88a, Ber91b, Yan93, Gja93], approximately symplectic track- ing [KSYZQI], evaluation of Lyapunov exponents and tune shift analysis [Sch91], as well as Nekhoroshev estimates [Tur90]. The principle underlying the proof of the Nekhoroshev estimate [Nek77] was evaluated numerically to obtain lower bounds for the survival time. We will call this the pseudo invariant estimation (PIE) method [WR92]. Although some of these methods are useful analysis tools, they all fail to 50 51 give mathematically rigorous lower bounds on the time particles stay inside the ac- celerator when the motion is described by a nonlinear map. We will introduce the PIE method by using pseudo invariants of normal form theory. This method comes close to giving guaranteed lower bounds on the survival time. To make the bounds completely rigorous, an extension of the PIE method is introduced which can ob- tain such a rigorous bound on the turn number. The pseudo invariant needed for this method is computed via nonlinear normal form theory, which was described in detail in chapter (1). The bounds are made completely rigorous by performing the required optimizations with interval methods. The use of interval arithmetic seems imperative for any rigorous treatment of the stability problem, since any tracking method only tests a small part of phase space of measure zero. The functions that have to be optimized are far more complex than typical applications of interval op- timization. Through the introduction of arithmetic on the new structure of interval chains (IC), one can exploit the special properties of the problem and enable interval optimization. Computations in this structure are performed by introducing a new data type into the FOXY language, which is the input language of COSY INFINITY [Ber90b, Ber92b, BZWH91, Ber93a, Ber94]. 2.1 Introduction to the PIE Method After a short review of stability analysis in weakly nonlinear systems, the PIE method will be introduced in detail. A section about normal form theory describes our choice of getting pseudo invariants. These nearly invariant functions of the one-turn map are essential for defining the beam region and for bounding the survival time. The importance of resonances and their influence on estimates of the survival time is discussed. Two refinements will be introduced which increase the obtainable bounds on long term stability. Up to that point, the PIE method assumes that the one—turn 52 map of the storage ring in question is well known. Since this is rarely the case, the theory has been extended to maps which depend on an unknown parameter. The PIE method hinges on efficient global optimization; when the required opti- mizations are performed by scanning the relevant volumes in phase space, the method cannot be completely rigorous. An introduction into interval analysis and rigorous global optimization will be given, followed by the definition of interval chains. The concept of interval chains takes advantage of the special structure of the problem and is faster than conventional interval arithmetic optimization by many orders of magnitude for the special functions that have to be optimized. Most of the demonstrated examples are described by polynomial maps. To ana- lyze general weakly nonlinear motion, it is necessary to find a bound on the Taylor remainder of a weakly nonlinear system. Wewill estimate a bound on such remain- ders by comparing tracking and one—turn maps. Such an estimate, however, is not completely rigorous. Interval chains are a subcase of a method called differential algebra with remainder or RDA [BH94a]. Like DA, RDA enables computation of the Taylor expansion of a function, but in addition automatically gives a bound on the remainder of the Taylor expansion in a given interval of the function’s domain. Toward the end, outlining possible improvements, it will be described how RDA can be used to also make the predictions for general maps completely rigorous. 2.1.1 History The investigation of the stability of planetary motion has been an important ques- tion for over a century. After early attempts by Laplace and Lagrange to under— stand the stability of the solar system, Poincaré [Poi99], Birkhoff [Bir27], and Siegel [Sie52, Sie56], among others, investigated the problem in detail. Usually the problem of planetary motion was analyzed by considering it as a perturbation of a known and 53 solvable Hamiltonian system. Innovative investigations of this problem were achieved by Kolmogorov [K0154], Arnol’d [Arn63], and Moser [M0362]. Nekhoroshev formu- lated a theory which estimates the time of stability of a system with a perturbation strength proportional to e by an exponential estimate. In reference [Nek77] the fol- lowing theorem is proven: (citation from p. 4) “Suppose that Ho satisfies certain steepness conditions, Then there are positive constants a, b, and co with the following property. Let 0 < e < 60. Then for every solution I (t), ¢(t) of the system with the Hamiltonian Ho(I) + eH1(I,¢), |I(t) — I(O)I < e” for all t E [0,T], where = fexp(;‘;).” This theorem is proven by performing canonical transformations ( I , 4)) —* (J, 21)) in order to minimize the dependence of the Hamiltonian on (l) as much as possible, thus bringing the new coordinates J as close to invariants of motion as possible, which is called creating “almost integrals” in [Nek77] on page 21. The exponential estimate is established by a detailed analysis of this canonical transformation, which is performed in a perturbative way in respect to e, and by finding the optimum order to which the transformation should be performed. For certain problems concerned with general Hamiltonians, celestial mechanics, and also single particle motion in accelerators, the Nekhoroshev method of exponential estimates has been used by finding values a and b for the specific problem. Examples for accelerator physics can be found in [Tur90]. The idea of the proof of the Nekhoroshev estimate has prompted an analysis of Stability of the nonlinear motion in particle accelerators by analyzing it in normal form space, a space in which the Hamiltonian has as little dependence on i/J as pos- sible. In this space the change of the “almost integrals” or pseudo invariants is not estimated by bounding the series of canonical transformations but by performing the canonical transformations on the computer and then evaluating their effect on 54 the pseudo invariants. This possibility was first mentioned and programed by R. L. Warnock for maps obtained by interpolation of individual tracking points. Later it was realized that Hamiltonians are not needed when the one—turn map or Poincaré map of a storage ring is known [WRGE89, WR89]. One only has to find the maxi- mum change 6 of the nearly invariant function during one application of the transfer map of the accelerator. This maximum change over the relevant regions of phase space bounds the change of the pseudo invariant for the entire particle motion in this region. Several other improvements on the method of pseudo invariants were made; they include using maps which describe many turns in the accelerator and differ- ent means of finding canonical transformations to the pseudo invariant coordinates [WR91, War91]. Most of this work was done by R. L. Warnock and R. D. Ruth in the Stanford Linear Accelerator Center, the fullest account is found in reference [WR92]. Following, our approach to the PIE method will be described in detail and applied to several examples, thereafter it will be shown how the predictions can be made completely rigorous and examples will demonstrate the applicability of the method. Our approach was described in [HBQ2b, HB93a, HB93b, BH94b]. 2.2 Pseudo Invariant Estimation (PIE) Storage rings are designed to hold particles for a long time. For example, the Large Hadron Collider (LHC) at CERN will have to allow particles to circle the 27 km long tunnel for one day at one millionth of a percent less than the speed of light in order to make effective high energy physics experiments; this corresponds to 109 orbits around the ring. It was pointed out in [WR92] that, if one considers the effect of every magnet as a perturbation comparable to perturbations encountered during one year of planetary motion, then this stability requirement corresponds to 1012 years of 55 stability in the solar system, far more than the estimated age of the universe. To ensure that the machine design is capable of holding the required numbers of particles for such a long time, it is important to develop methods which find out how long a particle with a given initial condition will remain inside the ring. This could be done by tracking the orbits of particles through 108 turns, which is far too time consuming with today’s computing power to be performed accurately. Another disadvantage of this approach is that the stability of motion can only be checked for a limited number of particles. Furthermore, the computational inaccuracies can build up to an untolerable amount. There are, however, some programs available which follow this approach using kick approximations for the optical elements to speed up the computation [Ta191, Sch91]. Other approaches look at the one turn transfer map that relates initial phase— space coordinates 2', to final coordinates after one turn 2'; = liq/1(5). This one—turn map contains all information about particle motion after many turns, since many turns are described by successive action of the one—turn map. The transfer map can be approximated in different ways. Recently an appropriate choice of B-spline functions and Fourier series has been applied [BWRF 93]. More commonly the Taylor expansion of the function is used. This has been done for light optics by Hamilton [Pra33] and was used for charged particle optics since the 19308. TRANSPORT and many other more recent computer codes follow this concept. In [Ber87] it was recognized that Taylor arithmetic in a DA framework allows to do this to ar- bitrary order. For the theory behind the DA technique, please refer to the references [Ber92a, Ber91a, Ber90a, Ber89, Ber87, Ral81], and for different DA programs, refer to the references [Ber94, Mic94, Yan94b, vZ94]. Time considerations often restrict calculations to about order 12 [Yan94a]. The accuracy of the Taylor map approach increases with proximity to the closed orbit. Once the one turn map is obtained, 56 particles can be tracked through the map to find out how long they stay inside the accelerator. Applying high order maps the required number of turns can still be very time consuming and, as in the case of element by element tracking, the stability can only be checked for a very limited number of particles. The PIE method analyzes the one—turn map directly without tracking through it several times, which in particular avoids computational inaccuracies. Furthermore this method does not only test single particles but provides information about all particles in a given region of phase space. We assume that there is a closed orbit in the ring. Particles with phase space coordinates near the closed orbit will not be lost, particles which are too far away from the closed orbit will be lost during their motion around the ring. We therefore divide the phase space 7’ into the allowed region A and the forbidden region ’P\A. The question we want to answer is: How many turns does a particle which origi- nates in a given region of phase space (9 circle the ring without leaving the accelerator. We therefore look for the number Nu,“ = max{n|A-f"((’)) Q A} , (2.1) where hf”(0) = {1VT"(5)|5 E O}, and 114W?) stands for n applications of M. The different regions are shown in figure (2.1a). With the following method we will find a strict lower bound N for me. If we find a real valued test function f that does not have common values in (9 and in 'P\A, then successive action of the map must bridge a gap A f as shown in figure (2.1b) in order to map a 2' E 0 into ’P\A. Particles start to bridge this gap by entering the phase—space region S.- = M(O)\O. The gap is bridged when a particle has reached the region 8f = M(A)\A. If S,- or S; are empty, particles in 0 will never leave A. If they are not empty, the gap goes from f.- to ff with f,- = max{f(z?')l§' E 5,} 57 0 A ‘6 Figure 2.1: a) The initial region 0 and the allowed region A of phase space P with O C A C ’P. b) The gap Af that has to be bridged. and ff = min{f(5)l§'€ 5;}. The function at. = 1(1)?) - f (2.2) describes how much f deviates from being an invariant of the map M; d; is called deviation function. When a phase space point 2' is mapped through 64 once, the gap A f is diminished by d {(27) If we assume ff 2 f,, the step from f,- towards f f is always smaller or equal to 6 = max{df(Z)IE°E (A\O)} . (2.3) A particle that starts in 0 therefore survives at least N turns, where N = [If—6:51 s M... . (2.4) We are thus left with four problems: 1. finding a suitable test function f such that N becomes favorable, 2. finding f,-, the maximum of f on Si, 3. finding ff, the minimum of f on Sf, 58 4. finding 6, the maximum of the deviation function in the appropriate region. To make the desired estimate as large as possible, we should find a function f which increases between the allowed and the forbidden region and should, at the same time, be close to an invariant of the one—turn map to make 6 as small as possible. The remaining three problems are concerned with finding maxima. These maxima can be found in a mathematically rigorous way using interval arithmetic, which will be explained in chapter (4). Doing this would, together with equation (2.4), give a mathematically rigorous estimate for the survival time of particles propagated by the given one turn Taylor map; this will be described in a later section. First we have to choose a suitable pseudo invariant f, and we are faced with the problem of describing the phase—space regions 0 and A in a sensible way. We will use nonlinear normal form theory on Taylor maps to solve both problems. 2.3 Normal Form Transformations and Pseudo Invariants Normal form theory for general order n symplectic maps is treated in detail in chapter (1). In the context of the PIE method we are only interested in the special case of stable linear motion, which means that all linear eigenvalues have modulus one. The phases of the linear eigenvalues are the tunes of the system, with 12,11 = exp(i27ruj), j E {1, . . . , d} for d degrees of freedom. For this special case, the definition (1.1.5) of resonances reduces to the well known resonance conditions for the tunes. Definition 2.3.1 (Resonance of Order m for linearly stable systems) The d tunes V: of a non-degenerate stable symplectic 2d x 2d matrix, are said to be in resonance of order m if there is a set ofd integers k,- with 221:, Ik,| = m + 1 and 1:11am: 0 (mod 1). 59 Theorem 2.3.2 (Normal Forms for Linearly Stable Maps) A non—degenerate linearly stable map M E 573:“(12) with tunes which a." not in resonance to any order S n can be transformed by a transformation H E S’Pfflfli) to If =n H 0 NT 0 H‘l’" E SP3,“(R) and ( 221—1 ) =7. ( C05f¢j(§ll “Sigh/JAE” ) (22] 221’— 1 ) _ (2.5) 22' sink/119)} COS{¢.(5)} The map 6: it" —+ E?" is given by b,- = 23,4 +23, These functions bj are order n + 1 invariants of the map If and the polynomial of order 72 +1 given by I,- =n+1 b,(H) are order n + 1 invariants of the map Ni. Please refer to the section on symplectic normal forms for the proof of this theorem, which is just a restricted form of theorem (1.5.4). It can immediately be seen that 1).-(ii) =.+. b.- . (2.6) which means that the functions bj are order n + 1 invariants of H. In a similar vain we write Ij(M)-=n+1 bjOBOM=n+1bjOROB=n+1bjOB-—n+1Ij. (2.7) Up to order n, the motion in normal form coordinates is described by rotations. The amplitudes of these rotations are invariant up to order n + 1. The proof given in chapter (1) demonstrates a method which allows direct computation of the normal form maps. During the calculations, Taylor terms of M) have to be divided by the resonance; denominators D1(lli) with D (ii): xp(i27rz/)) — exp( (i21r 2 10.11,) . (2.8) i=1 Since the resonance denominators vanish at resonances below order n+1, the theorem on symplectic normal forms excludes tune resonances. However, it can be shown that 60 there is no problem with resonances of the linear eigenvalues when the map A7! is the Taylor map of a system which has d exact invariants of motion. To illustrate the normal form transformation, the motion in phase space for 2000 turns in a typical accelerator is shown in the left part of figure (2.2). For each turn, the horizontal position a: as well as its canonical conjugate momentum a is displayed. The finite width and the irregular structure of the band is a result of nonlinear effects and of coupling to the other degree of freedom, the motion in vertical direction. From the picture one can see that the particle positions are bounded for the number of turns shown. However, it is very difficult to estimate if particles are on average moving away from the _origin or not and what would happen if the number of turns were increased. The right picture in figure (2.2) shows the same motion after transformation by the normal form map. The motion now has nearly circular shape, which will make it much easier to estimate the long term stability. “a. Figure 2.2: Phase Diagram for 2000 turns in an accelerator for four initial conditions. The left picture shows the motion displayed in standard particle optical coordinates x and a, and the right picture shows the same motion in normal form coordinates. 61 2.4 Maps to Test the Method We will use six different maps to test the normal form invariants and the PIE method with and without interval optimization. Two of these example maps have one degree of freedom and four have two degrees of freedom. One map on two and one on four dimensional phase space is known to have completely stable motion; these maps are suitable to check what PIE yields for systems where infinite survival times can be guaranteed by analytic calculations. 2.4.1 The Physical Pendulum A pendulum of length l and mass m in an uniform acceleration g is used as a model for stable motion in two dimensional phase space. The chosen position coordinate :r is the arc length of the pendulum’s elongation. The Hamiltonian is given by NIH p2 H = -— — mglcos( ) . (2.9) 2m In order to compare this motion to motion in an accelerator, we transform from the time t to a length 3 as independent variable, pmcm is the maximum momentum for bound motion, 3 = tpm” = 2t\/;l. (2.10) m Furthermore, we introduce the dimensionless canonical momentum a = p/pmax. With these coordinates and the independent variable 3, the motion is canonical for the new Hamiltonian 2 1 G: ZxHhmpmx): %_ZCOS(T)' (2.11) The motion is governed by Hamilton’s equations d d _x = 6.0, —a = —6.G , (2.12) ds ds 62 and the linear map becomes 1:; _ cos(21ru) 2l sin(27r1/) 1:,- _ _s_ ( a; ) — ( —%sin(27ru) cos(21rz/) a,- ’ V _ 47rl ' (213) The phase u corresponds to the tune in an accelerator and the invariant ellipse of linear motion is given by 2 Il(.7:,a) =2 % + 21a2 2 e. (2.14) when the emittance is 67F. This invariant of the linear map agrees up to order 2 with the order n + 1 invariant 11 of the nth order map. In the examples to follow, we chose an emittance of e = 10,000 mm mrad and a tune of u = 0.379 if not stated specifically. The nonlinear map is computed by evaluating the exponential Poisson—bracket operator acting on the identity up to order n: M 2,, exp(— :47r1/lG :)E'. (2.15) As introduced in theorem (1.3.4), : f : g denotes the Poisson bracket between f and g. If not specifically stated, we will use order 8 in the examples. 2.4.2 Henon Map for One Degree of Freedom The Henon map [LL83] is a standard test case for the analysis of nonlinear motion, because it exhibits many phenomena encountered in Hamiltonian nonlinear dynamics. These include stable and unstable regions, chaotic motion, and elliptic fixed points. The Henon map can even serve as a very simplistic model of an accelerator under the presence of sextupoles for chromaticity correction. The figures (2.3a,b) show typical tracking pictures for the Henon map $f _ x o cos(27ru);r,- + sin(27ru)ag (9 16) a; _ a + kmz — sin(27ru).v.- + cos(27ru)a,- ’ -. which is a composition of a kick map and a rotation. 63 Figure 2.3: 9 particles tracked for 500 turns through the Henon map. The starting conditions were a: = 0.1 ~j, j E {1, . . . , 9} and the kick strength was 1.1. The tune in figure a) is 0.255 and in figure b) 0.29. Since kicks and rotations are symplectic, the Henon map is symplectic. One can find Hamiltonians which have the Henon map as time step one map, for example, H = 1r1/{:1:2 + a2 — 2(tlc)a:c2 + (tk)2:c4} — glen? . (2.17) However, it is not possible to find a time independent Hamiltonian, since in general the Henon map does not have an invariant of motion. To compare with accelerators, the tune corresponds to V and the linear invariant ellipse for emittance err is 11(3, a) =2 x2 + a2 = e. In the examples the emittance will be chosen to be 10, 000mm mrad and the tune will be V = 0.379 if nothing different is stated specifically. 2.4.3 Coupled Pendulum Coupling two pendulums is not very suitable for comparison with an accelerator, since in linear approximation the motion does not decouple for a natural choice of coordi- nates, whereas in accelerators the linear motion in horizontal and vertical direction typically decouples. We therefore choose a physical pendulum with an elastic string, 64 a bungee jumping point mass, so to say, which behaves corresponding to mid—plane symmetry. The figure (2.4a) shows typical tracking pictures for this pendulum. The motion must be bound for small enough amplitudes, since the energy is a Lyapunov function [Lya92]. o . . on" ‘ .Iu : g o. 0.. : o . a. Q. . o '.‘ o. ‘. “5."? .l-VNW‘“ . gm, . 'Q’. ~ {l}: a - ' ":----— I - ~:::: \-" U ‘WHM‘,W “WW "‘41. 31mm...“ "aw-m... ..... g"; . W “l... . ‘.' 'f V" ‘ Q .ar ‘. r...‘.l ..”"~ 32"...” ‘I. ‘1‘.“'.';H‘.M; 5.0 kiwi-if” Figure 2.4: 9 particles tracked for 500 turns through the coupled pendulum. The starting conditions were angles x/lo = 0.01 - j for j E {1, . . . , 9}. This corresponds to a maximum emittance of 4050mm mrad. The horizontal and vertical tunes in figure a) are 0.36 and 0.78, and in figure b) 0.37 and 0.78. The pendulum has a relaxation length lo and is elongated at rest by A to the rest length lo = lo + A. The canonical positions are :r, the angle times lo, and y, the elongation from the rest length. Lagrange’s formalism leads to the conjugated momenta p: = msi:(1 + y / lo)2 and p1, = mg]. The Hamiltonian is given by H = —(p§ + ——l’£—) — mg(lo + 31) cos(% o"‘)+ 2 —(y + A)2 . (2.18) (1 + 341012 A The frequencies of linear motion are u), = ‘fg/Io and my = ‘fg/A. Like for the physical pendulum in subsection (2.4.1), normalized coordinates will be introduced. The maximum momentum for bound motion is pm” = 2m\/glo and szte-M a=—?-”— b=—l’l’— G: mH. (2.19) 3 t 7 m pmax pmaz: pmaz 65 The new Hamiltonian is 1 1 G = — 2———— 2(a (1 + 31/10)? (31 + AV _ EFL/’2 C ______._ 4 2? (,2 + ) 10 os( )+ (2.20) As in equation (2.15), the Taylor map is computed by M =n exp(— : 80 :)E'. In linear approximation, the map is given by it cos(21er) 210 sin(27er) x,- ( a: ) = ( —§%;sin(27rV,,) cos(27er) ) ( at. ) , (2.21) 333' g, cos(21rVy) 2 10A sin(27rVy) a,- = 1 . (2.22) b, —m sm(27rVy) cos(21rVy) y,- bi The tunes depend on the independent parameters and their ratio depends on the 3 lo V;- — 47m) , Vy — V,\/;. (2.23) The invariant ellipses corresponding to the emittances earn and eyrr are parameter of the pendulum, $2 2 _. I 2 = 1:, _y 2 = . . 210 + 2 0a 6 2m + 2V loAb Cy (2 24) If not mentioned differently, the pendulum will have length lo = 1m, tunes V1- : 0.17, Vy = 0.91, emittances of 30001rmm mrad, and its map will be evaluated up to order 8. 2.4.4 The IUCF Ring As an example storage ring, we used the ring of the Indiana University Cyclotron Facility (IUCF). The electron cooling device and the RF system were not used, since we want to compute the stability of motion when no such devices are present. The device is usually used for emittances of 0.31r mm mrad, which are made so small by electron cooling. Since we want to analyze operation without cooling, we assumed 66 ex = 3.71rmm mrad and «51’ = 2.21rmm mrad if not otherwise stated. The uncertainty of the first magnet’s strength is assumed to be 0.01% and the linear tunes are chosen to be V;- = 0.7727 and V,, = 0.6650 in the examples shown. «(f-ll II -illth? .3. ‘ | 3 /§\ \ \\ H (7’)] \l \ it x t '/ 79‘ f? / Vx.’ \ \\S/\ s/ 9'4 I. . xx” \Q'.‘——ju:\/ Figure 2.5: a) Layout of the IUCF ring, only dipoles and quadrupoles are shown. b) :c— Px phase Space positions for 500 turns. The initial conditions are (x = 5j, y = 5) . 10“ withj E {1,...,10}. 2.4.5 The PSR II Ring A second example storage ring is the PSR II, which was designed as a possible upgrade of PSR at the Los Alamos Meson Physics Facility (LAMPF). This device was analyzed for emittances of 401rmm mrad and for a uncertainty in the strength of the first magnet of 0.01%. The linear tunes are chosen to be VI = 0.2313 and V3, = 0.2705. 2.4.6 The Demo Ring The third storage ring which will be analyzed is a simple example ring. It consists of three identical cells. Every cell is mirror symmetric in the elements which influence the first order. The dispersion also is symmetric, which makes each cell an achromat. 67 Figure 2.6: a) Layout of the PSR 11 ring, only dipoles and quadrupoles are shown. b) :r—px phase space positions for 500 turns. The initial conditions are (a: = 5 - j, y = 5-j)-10‘3 withj E {1,...,8}. The chromaticities, which describe the linear dependence of the tunes on energy, are canceled by two hexapoles and the effect of nonlinear resonances is reduced by two further hexapoles. These hexapoles are adjusted by reducing the nonlinear resonance strength which can be computed by the normal form method. The 12 hexapoles are not shown in figure (2.7). This device was analyzed for emittances of ex = 57r mm mrad ey = 71r mm mrad and for an uncertainty in the first quadrupole of 0.01%. The linear tunes are V; = 0.37 and V3, = 0.67. If not otherwise stated, the transfer map of the Demo ring will be evaluated to 8th order. 2.5 The Pseudo Invariant and How to Parame- terize Regions Describing the initial region and the allowed region is essential to finding 5.- and Sf and therefore to finding a function that changes substantially between those two regions. Accelerators, usually have mid-plane symmetry, which implies that the linear map 68 ‘ Figure 2.7: a) Layout of the Demo ring, only dipoles and quadrupoles are shown. b) :r—p, phase space positions for 500 turns. The initial conditions are (x = 5 - j,y = 5-j)-10-3 with j e {1,...,8}. does not couple the x~px and y—py component of motion. The projection of the linear motion of beam particles in the :r—px subspace lies on invariant ellipses. The area of these ellipses is called the a: and y emittances of the beam. Since the product of two ellipses is topologically a torus, the linear motion is said to lie on an invariant torus. The allowed region for a beam in a storage ring is typically given by the acceptance in the x—px phase space and the y-—py phase space. The a: and y acceptances are defined as the largest invariant ellipses that can be transported through the accelerator in the horizontal and vertical plane respectively. The reason for giving the acceptance as an invariant ellipse is described as a peel— off effect. After every turn, a particle is on the same invariant ellipse in the .2: and y section of phase space. Turn by turn the particle’s position rotates around these ellipses. The angles of these rotations, which are 2n times the tunes, are chosen not to be a fraction of 21r to avoid resonances. If some obstacle in the beamline touches an ellipse, sooner or later all particles on this ellipse will hit this obstacle and get lost. 69 Any obstacle therefore peels off particles lying outside an invariant ellipse. The rotations in nonlinear normal form theory correspond to the concept of ro- tations on the invariant ellipses. It is therefore natural and follows from the same argument that in a nonlinear theory the allowed region or nonlinear acceptance .4 should be given by a nonlinear invariant of the map, which can be computed by normal form theory. According to theorem (1.4.4), normal form theory gives, one invariant circle for each of the (1 degrees of freedom. The pictures a) and b) in figure (2.8) describe the nonlinear invariants which specify the boundary of the allowed region. Since the linear contribution in the map dominates, they are close to invariant ellipses of linear motion. There are many different invariant surfaces which can be described by the d non- linear invariants I,-. The two most suitable ways will be discussed here: d .n==mimmmsn, (ma i=1 An = {£11,(5')Se,,ViE{1,...,d}}. (2.26) The corresponding beam shape is expressed in figure (2.8c) by drawing the biggest allowed a: and y coordinates. Elliptic beam shapes correspond to A0 and rectangular beam shapes correspond to An. To keep the notation simple, we describe the ini- tial region 0 in a similar way with nonlinear emittances 0637?, where a < 1. The boundaries are most easily described when a norm is introduced which measures the distance from the closed orbit according to the invariant torus on which the parti- cle moves in nth order approximation. If we want to represent the beam by a round shape, we choose ”5] lo, while if we want to describe it by a rectangle, we choose Hi! In, where ‘1. Mt=zfl m; “mammaLam €{LMJH. on) 70 In real normal form space, to which the transformation 3 in theorem (2.3.2) leads, the beam region can be suitably represented by a set of invariant tori. These tori can be transformed into phase space by the order n inverse 13““. We parameterize such a curve by the emittances e, which describe the torus in normal form space by writing 7(6) ={é-1vn(2)|( 22f“ ) = Jail 3%) ) ,¢. 6 [0.2”],vi e {1,...,d}}, (2.28) With this notation, we can introduce regions in phase space which can describe the initial region and the allowed region. Let 7:460 d . {515€7(77),ESZ%SC}, (2.29) i=1 ' $0030 = {515ET(77),€S%SC,W€{1,-mall}, (2-30) then the regions of interest are given by 0, = f,(0,a) , A, = f,(0,1) , s E {0,0} . (2.31) Since 3’” is an order n inverse of E, this definition is approximately the same as Hill, S C or llfll, S C with C E {0,1}; therefore, the pseudo invariant llEll, fluctuates very little on the surfaces of 0, and $1,. Phase space points in the regions of interest are easily parameterized by the choice of the (b,- and 77,. It is worthwhile to note that the first order inverse 13"” is the exact inverse of the first order of H. 5"” transforms circles into invariant ellipses of linear motion. Therefore, when n = 1 is chosen in equation (2.28), then the conventional definition of the acceptance is obtained. Since the polynomial map H’l'“ is continuous, it maps closed regions of normal form space into closed regions of phase space. Therefore, this definition of the acceptance is meaningful; particles can never leave the region of 173(0, C) without crossing the surface fs(C2 C)- 71 To make the desired estimate as large as possible, we should find a function f which tends to increase when the norm increases and should at the same time be close to an invariant to make 6 in equation (2.3) as small as possible. The appropriate choice is fo(5) = ”Ella and fu(5)=||51|u- There are three reasons which suggest the use of f0. First, for d degrees of freedom, evaluating fa takes d times longer than evaluating f0, since d pseudo invariants have to be evaluated, whereas for f,, the polynomials Ig/ég can be summed before evaluation. Second, since the beamline is generally circular, it seems more appropriate to choose this notation. A third reason comes from a somewhat heuristic argument. The figures (2.9) shows that the average invariant defects [All-f) — I,- are of the same order of magnitude on the surface of A0, whereas they fluctuate substantially on the surface of Ag. If one accepts the heuristic view that bad pseudo invariants in the normal form picture are an indication of a low live time of beam particles, then it is more appropriate to specify the beam shape by llE’llo. Following, the quantities with the circular subscript will be used and the subscript will be dropped. A py A J]:- AP J5 HZHU=1 \\ P— \ / :3 \\§’ \ ) llzllo=1\ \ / ‘ «5‘ «I: 4> Figure 2.8: The motion on nonlinear invariants in the phase space section :c—p, in figure a) and y—-py in figure b). The allowed and the forbidden region and the definition of Hill is depicted in figure c). 72 0 7r/2 0 7r/2 {> D Figure 2.9: Fluctuation of the pseudo invariant on the surface of the allowed region. Elliptic beam shapes (3 = o) and rectangular beam shapes (3 = Cl) were used for the left and the right picture respectively. The four cases and the emittances used are the coupled pendulum (£1 = 0.66, 62 = 0.46), the IUCF ring (61 = 0.36, 62 = 0.76), the PSR II (61 = 0.056, 62 = 0.956), and the Demo ring (61 = 0.46, 62 = 0.66). 73 2.6 Analysis of Pseudo Invariants The PIE method relies on the choice of the function f, which should be as close as possible to invariant under application of the map M. In order to get a sense how well the order n + 1 invariants of the normal form transformation are invariant under Ad, typical cases were tested. These cases are the six systems described in section (2.4). In figure (2.10) the invariant defect is shown for the physical pendulum. The deviation function d f = f (M ) — f with the pseudo invariant f was plotted over the region in phase space parameterized with equation (2.29). The coordinate axes in the picture are polar coordinates (fl and (b. In high orders the invariants are extremely accurate. This accuracy is only possible because the physical pendulum has an exact invariant, the new Hamiltonian G of equation (2.11). In table (2.1) the smallest resonance denominators in every order up to order 10 are displayed for the Physical pendulum and the Henon map. The denominators are equivalent for these two examples, since we chose V = 0.379 for both. All reso- nance denominators have absolute values which are sufficiently big to avoid divisions by dangerously small numbers. However, the normal form transformation requires successive divisions by these denominators and after multiple divisions, very big and inaccurate coefficients can occur in the normal form map. This fact however will not have an influence on the reliability of predictions given by the PIE method, since no exact invariants are required. Also for depicting the invariant defect for the Henon map, polar coordinates (fe- and (b were used. The accuracy of the pseudo invariants is less than for the pendulum, since no exact invariants exist. 74 Table 2.1: Smallest resonance denominators for the physical pendulum and the Henon map with a tune of 0.379. ///”. ’ ”I“. l .0333” It ll .. airmail “Mitt" ll l ’ ninth“ . «3,. t «Nil 14;", , .....i n \\ \ am I \\ 0,134.?“ 4". 0 n. O Figure 2.10: Invariant defect for the physical pendulum on an invariant ellipse. The coordinates are 0 to 6 from left to right and 0 to 27r from front to back. The range of the depicted function is [—2.36 - 10‘12,5.63' 10‘12]. 75 ..... ."o"“ | \ \ —-::.—.;:; «ravi‘ax‘tma‘a, 7"“ " o‘.‘c‘ .“‘o°“‘\‘ ‘ - .0 .0 .o‘ .o‘ - ;. ::.:::.-.;.;-::.‘m‘\‘\\\\\\ . Z 0. .... o... o:::0:‘::0:‘ . . ::o:’:o:’:.:.... “““ , ......;..;.;;;.:.;::: ‘ :-:-:'-:-‘:::“‘I‘I‘I“.\\‘“\\\\‘““l » a}... . - -e._-' «.53.: ‘:.‘.§9.‘3s\‘“‘\\\6 .\“‘I| <.‘ ' - . - - . _, "o"- -‘ .“. . ..‘. 9‘. — -'- ."O ‘0 . .-,. Figure 2.11: Invariant defect for the Henon map on an invariant ellipse. The coordi- nates are 0 to 6 from left to right and 0 to 27r from front to back. The range of the depicted function is l—5.59 ' 10‘8, 5.56 - 10‘s]. 76 The other four examples are systems for four dimensional phase space. We will therefore depict the invariant defect on a torus 7(6). The coordinates in the figures are the angles 451 and d); which parameterize this surface. In figure (2.12) the invariant defect is shown for the coupled pendulum. In high orders (here eight) the invariants are extremely accurate. When comparing the accuracy of the coupled pendulum to the other systems, it should be noted, that the emittance of the pendulum is about 1000 times bigger. This accuracy is only possible, because the physical pendulum has an exact invariant, the new Hamiltonian G of equation (2.20). In table (2.2) the resonance denominators up to order 10 are displayed. In figure (2.13), the invariant defect is shown for the IUCF ring. The invariants cannot reach the accuracy of the coupled pendulum, since no exact invariants of motion exist. In table (2.3) the resonance denominators up to order 10 are displayed. The tables (2.4) and (2.5) show the smallest resonance denominators for every order for the PSR II and the Demo ring; the figures (2.14) and (2.15) show the invariant defects for these systems. Order 1 1‘71 k2 emu, _ ei21r(k1 U1+l€2 V2) 2 l 2 0 0.42830104E-01 3 1 3 0 011284027 4 1 3 —1 0.64523782E-01 5 1 4 -1 0.11918375E-01 6 1 5 —1 0.12615773E-01 7 1 6 -1 0.83713329E-O2 8 1 1 —7 0.10957355E-01 9 1 2 —7 0.17465246E-01 10 1 4 6 0.41817583E-02 Table 2.2: Resonance denominators for the map of the coupled pendulum. 2.7 Influence of Resonances In the chapter about normal form theory it was mentioned that in general a nor- mal form transformation is only possible if no resonance condition up to order n is 77 u {l . \ WWII“ .’ \ e~ I m "ii‘i'l‘ llllvi“ “am“‘h \\‘u no...» (illlilhltttill‘lllllllllllli ‘. /‘ 99‘9 ' 3;.“ / It) «I: ittvlliiiiil‘x‘tt‘t‘ ('2. :l‘ I it“) = \“““:.3:42.1'346‘ \‘ “\\\\)\‘_.‘:‘ {3:7, ‘I‘ ‘I‘ ‘|““‘\““ ‘\\“ ““““ “II “II \ \\) one/’0, ‘I‘III IIIIIII “I9“ " ‘\\\‘\M \\ “I“ “I ““\\\\l\\‘“‘i..“‘o n it’xl‘,“ \\\\“‘.:2“‘2§i \lllil““2‘““c« \\\\ ‘\‘::’ :;;':g| \\\\\\\\\\\\\\\l “‘9. at“ ”twink ::«“‘ \“\\“ \‘|‘“““\\\\ “an“ Will 1‘ ‘9‘.II““\\\ \“ “:u‘“ ““\\\‘\\ M‘ “|“‘,‘ “9“,. (a; 9“ “I\ ““9 “ mu" :a‘u‘I “N Ill Figure 2.12: Invariant defect for the coupled pendulum on a pseudo invariant torus The coordinates are in [0, 27r] x [0, 2r]. The range of the depicted function is [—2.55- 10-13, 7.36- 10-13]. Table 2.3: Resonance denominators for the map of the IUCF ring. 78 /’ uh .I’III‘I‘IW‘ //’t ‘\‘ N)“ ‘9‘“ W liti’t‘ill) h) Iiil‘llll‘ll «Kg/MW“ Om‘)‘ llli ““‘,'/,,I‘I\I\\\)\\\ . lift ' lli‘) (III, .: I)“ “)“)\‘“ III“ \\\\\\\\\)~) ’l”. ‘l‘(:\‘ “\ W “‘“\‘\\\\\0 "I/“I‘)“)“‘ ‘illjvil’fl’i'lfk‘Il‘I‘l‘l‘ll‘lll “\\\ ,’:“‘\ \‘\‘\\‘\)G‘ " / \‘ “ III, ’IMI‘ W“ \\\\)\\\\\I“o"o'IIHtI )‘ 353.06". I, ’.9‘\\\\ i) I” “II N“ ,ff I,” / ’IIIII”“:\:‘\\. ‘I‘Il’ I’éQI‘I\““\ l\‘\)\\\\\\\\\“\‘\‘\“ I. III/09,601” .|\\\‘\( (\\\\\\‘\,I/ 33‘“ \\‘\\\\ III". ”I’ll, 333‘) ‘il‘lllztlfi ‘\‘\\) ('l'”IIIuII“‘\\‘ \\ If”. ‘9“““\\\) ‘““ l‘“ “‘ N“ t\ H‘l\\ )\‘\\\\\\\h \l\ ‘\ 30’” ,,”I’ I‘\\‘\“‘)‘ \‘\‘“’ ”9‘ \\‘\l\‘ ('5‘)‘ “\\I‘\‘,-;~ , 0’" \II‘IINI;I,'I/,/,£I’I:::\‘\\‘\‘l‘ “‘19:...“ ‘3‘)“ \“)“ mm“ 3?, VIII,” “IlIIj’l/I “\‘\‘\\\‘\\‘ MON.” ‘7’- \\ “‘3’,th III II a“: \‘. .4 l \\ \Q’. ”A Figure 2.13: Invariant defect for the map of the IUCF ring on a pseudo invariant torus. The coordinates are in [0,21rl x [0,27r]. The range of the depicted function is [—6.84 - 10‘9, 6.00 - 10‘9]. Table 2.4: Resonance denominators for the map of the PSR 11 ring. 79 .— ... "h, _ "lll/I/[I II] [I’ll/{35W VIII, ’I’II’II‘” ‘ ‘ . 9 ‘H ,....,,,3 9,, m“ :‘t’i ‘.::‘:,> ‘3: 0 00 0‘9‘o“ 5 :E:.: “9‘.0.‘9“ 5.... . ’:.: t. I. C ::4, 'v “9 \‘:t§!h_.l i W‘I‘I‘I‘. 00 I. IIIIIII ,,'0'I9 h.“ IIIW ’4’?“ ,0,00 “""'Ii’/’ 3,3,, ‘53,), ,3,” 09,3I,90,I,‘,‘5959‘I‘ m’; ‘9,,.),‘ m9;9,’3’9,,,,,3 53:53,. 5‘I:,II:,III‘ :,:II, 3303:), “I: II “xx/3 ,O,I, ”3 9,9,9“. 353'533335I55III,I‘5 99 9 0’ “) ‘IIIII MI“ )9)‘ I939 ”"9" “9‘99“0‘9“ ““‘5‘9 9,99" I‘ ‘“‘I“9“.“..5 909,9,‘9‘,“‘i‘~‘“ 0‘9‘9‘35‘5‘.) 0’9 ‘IIIIII, II‘9‘III ”““‘9‘“‘)‘I‘9 “II‘” 9,9, )5. 55,999. 55 ,3,,,,5555999“9‘9“1II9I9,0999 III, \\\‘ w“°“ ‘I‘ ::: ::. {II ..... 009 33,313,, 5,)“, ”)H II.) ‘II... 3,, 9999995., )9 ‘9‘ ‘0“9I‘5“ ‘9),. ‘W‘ 5 ‘5‘9‘9‘9“ I N 9959“ ‘I“““‘I ““ ‘9‘9‘99999‘ ,‘9“‘5“‘9 9‘“ 3 H’I‘I‘N )9) 999 ““ ,3399999‘9‘5g3359 9.9 ,5 9955555,: 0,30,39,I3.3 ‘55 5,), ...... “9“9‘519‘ \‘\\)9 “I A \efit: 0'. , . ...”! ... Figure 2.14: Invariant defect for the map of the PSR II ring on a pseudo invariant torus. The coordinates are in [0,21r] x [0,27r]. The range of the depicted function is [—1.76-10‘9,1.65-10‘9]. Table 2.5: Resonance denominators for the map of the Demo ring. 80 ”/0”; “\ /’ ’ ‘ O‘m //// I ’I ”I' ‘-“‘\\ / ’I” Iceman ,..;.(////, ’///,'I,,,,,. ”I ’I “‘\\\ \\\\\ II’ “Vt" // //////,’o'o , .. ‘“\‘\‘ ,////’I’fi'o:.‘\‘\\‘\\\;\g§\\u\ / /’ ix...‘\\f:,/////{f{g,fo:¢:. ' III ”'0’ III’IJ’I’IMeé’/” g \ I , ;., a»; ‘\ \‘,‘/ 3:. O / triti‘fiitvie‘éfii‘s I;I/I.I.‘v.‘eé'Iqil/éifzstt'o'o'o'tt';I.6.I.W.«\\ \ ,////,,,I,I,;;., "I", “.““\\\\‘\““~“ ’I ' ‘ \ ’60. .Q 9” ’ ‘4'; OOOgfl“ ““ ‘|“.o‘ a, //// ll, ’11" w \ \ “w, I/ l I \\~ Mm mu ’IIIIc¢\I¢g.\\ .«w , //// ”II/I . “-3333“? I, ’ ’ 'O. ‘ ‘ {5‘ O. Q‘ .. ”//I ’0. ‘ ‘\‘\:‘.~‘. \ \ \“ ‘:’:’;/.v1/// Ill/”#19" W O I um \m‘h‘m't // 0’6}. - NW ‘3: WW?» I 1%“. t ..N H to.“ l!‘.\“\\ .\|“ “\‘,o‘.o::///// //// bog-J“ “ ‘v, y 9". 9\\ 1,, 5,99... . OOO' ’I’II".:“‘\ ‘07“..‘““ \\ ““‘“.:’.ot';' ////// l/Il’;0.0;.v. ‘\“‘“ ‘1.””’ O l.“\ ,;:.0:.. .‘ ‘. ‘. NH 5',’Iog\\‘\‘\u;:.\\\\\ \\ . 5.3%,.”- ,‘ //// ”We, ""0“?!” ”i," ‘H§\\\‘«:.‘.«“.t ‘| “ ”Wu-"I \H‘x‘x‘uee. ;‘ \w\\‘«“"'”’///// ////////,”’I/o'o“o‘m"fl O" NH “Wk“ « 6 N6: /9‘\‘ \‘o‘ ‘\“\\\\“\‘ ““ Wm // .. ”MM‘N’MH N N t“\‘\\o‘\‘m\‘«‘o m" \s‘s. \\\\ “v.3 $5,049,» ////I 14,599... M"! OW Mb” w‘mm M. ,1 ’ ~\\\\\\“v‘&' ‘7’ 'm's‘\"" 4'90"“ WNW N 0'! 0.‘~.~\\\\\‘\m\“‘ -.-. -“‘I"/’l"N’ v '-:‘::«\“N mm H i I... ‘ a: ‘9' ' "I ’ ' . “\\“‘\“““ ‘ .. ” ’l’ 'I’ \\ ‘9'6/ 9"¢I0 MHHONI.II m, WW iii'fl‘i‘iiifi:333v35‘3o‘o66f93'3"."" \: . Woiti'ti'o.‘twtxx\“.‘xé‘s'xxef"-" "~" 0 m m" ’ “ ““ ‘:\‘“‘::::E ””16”,." ‘33::3/ Figure 2.15: Invariant defect for the map of the Demo ring on a pseudo invariant torus. The coordinates are in [0, 21r] x [0, 21r]. The range of the depicted function is [-7.22 - 10'9, 6.95 - 10‘9]. 81 satisfied. It was also mentioned that this transformation can always be performed, even when resonances are present, if the motion is integrable, meaning there are as many invariants of motion as degrees of freedom. In general this is not the case, since probably most real systems are non—integrable. We therefore have to face the fact that resonances up to evaluation order have to be avoided in order to perform the normal form transformations needed for the PIE method. The figures (2.16) show the tune space for two degrees of freedom. In an accelerator these freedoms are the motion in vertical and in horizontal direction. The resonance conditions create lines in tune space on which the tunes are not allowed to lie. These lines are shown in figure (2.16a) for resonances up to order 6 and in (2.16b) for resonances up to order 12. Similar graphs are shown in most elementary accelerator physics text books and elementary papers [C858, E893]. In these graphs the tunes vary between 0 and 1 / 2. This quadrant of the complete tune space is depicted, since the set of resonance lines is mirror symmetric in respect to the lines V1 = 1/2 and 11,, = 1/2. If the tunes lie close to a line, small denominators 01(5) occur in the normal form calculation, the transformation becomes inaccurate and the quality of the pseudo invariants decrease. The reason for this fact is not only a computational problem but there are physical reasons, since at resonances, good pseudo invariants do not exist. To analyze this fact, we computed the effect of the tune on the pseudo invariants and thus on the long term estimates of the PIE method. For the Henon map, we scanned the tune in 200 steps from 0 to l. The variation of the deviation function is recorded in figure (2.17a), where a logarithmic scale was used. The rapidly increased deviation from invariance and thus the rapidly decreased survival times close to tune resonances can be clearly seen for every single resonance up to order 8, which was the evaluation order of the normal form transformation. This figure suggests the well known fact that it is advisable to keep the tune of storage rings out of “resonance 82 Figure 2.16: Regions of resonances for two degrees of freedom up to order a) 6, b) 12. holes” . In figure 2.18 the tune of the physical pendulum was scanned from 0 to 1 in 200 steps and no resonance problem can be observed. Normal form theory is able to ap- proximate the exact invariant given by the Hamiltonian in equation (2.11) even when the resonance denominators vanish. This also indicates that in the case of the Henon map, the problems at resonances are not only due to the resonance denominator problem. In order to see if this effect is only of a computational nature or if the physical reason of non—existing pseudo invariants of good quality is important, We used the Henon map and performed only a fourth order normal form transformation. The re- sult is depicted in figure (2.17b). The invariance defects tend to become larger, since now the invariant surfaces are only approximated to much lower order. However, resonances of orders higher than the evaluation order cannot be observed. If the de- creased invariance of f were only due to physical reasons, such resonances should be 83 seen. It can also be observed that with higher orders the invariance decreases at res- onances. This is not too surprising, since nonlinear normal form transformations try to approximate a non—existing invariant. We can therefore conclude that the normal form invariants are the Taylor expansion of true invariants if these exist, however we cannot conclude that they are the best possible pseudo invariants, especially close to resonances. The figures (2.17a) and (2.17b) also show that there is an optimum order for normal form transformations, which can be quite low if the system is close to low order resonances. For the example system we analyze, the tunes were sufficiently far away from resonances, such that this maximum order was not yet reached at order 8. This can be seen in the tables (3.1a—d). Figure (2.19) shows a corresponding picture for the deviation function of the PSR II ring. Calculations were performed in order 8. The tunes xxx and 11,, were scanned from 0 to 1 in 100 x 100 steps. In order to obtain comparable systems, the map of the PSR II was computed as discussed in section (2.4) and then composed with a linear map to obtain a symplectic map with the desired tunes. Again the figure is shown with a logarithmic scale. The boundary of the tune space has zero tune and is therefore excluded. Resonance lines up to order 3 can be clearly seen. The lines which can be observed are one second order resonance line, Vx+uy = 1, and four third order lines, ux+2uy = 2, ux -- 21!” = —-1, 11,, + 2V1, = 1, and V1, — 2Vy = 0. Influences of other resonances cannot be observed. This indicates two things. First, not only the tune is responsible for the quality of pseudo invariants, and second, the increased deviation from invariance is not only an artifact of the normal form method but reflects properties of the system. Figure (2.20) shows many more resonance lines. Resonances up to order 5 can be seen. Again it is apparent that not all resonances are observed. It can be imagined that survival time versus tune plots of this or similar kind can be helpful for machine 84 A 4.33.10—13— 1.00 1.94-107 All! 11]-]; l nJLL A1 II A D 1111 12 l§2§ l iééz §§ $§§Z 8765 473857 2 758374 5678 434-10—9— 0.14 1.00 -1IL LA [.11 1 nl-‘ ll ILA- D 1.1.1; 1.2 l§2§ l fiéél §§ 1§§Z 8765 473857 2 758374 5678 Figure 2.17: Variation of the maximum 6 of the deviation fllIlCthIl‘df with tune for the Henon map of a) order 8, b) order 4. The scale is logarithmic and inverted. Staying out of resonance holes yields long survival time predictions with the PIE. 85 A I6.94 - 10"l7 1.88 . 10"14 l> con—l «ll—‘1- Ozlh‘ " cult-H- Ah-o #- NICO!- wiv— r- WIODI culm- le" NIH Il- Nib I- O‘lw I- Glen» (JIM- film '- biw— CMA -- Oblon- NIIO) '- com» Figure 2.18: Variation of the maximum 6 of the deviation function d; with tune for the 8“ order map of the physical pendulum. The scale is logarithmic and inverted. Due to the existence of an invariant, resonance denominators do not cause a problem during normal form computation. 86 I I’d-l I ”I I'I'II‘O' ll ff!!! 7‘ III! III I .I all! '5‘ I! ‘00! 7‘ II!!! Variation of the maximum 6 of the deviation function d f with tune for the 8th order map of the PSR II. The scale is logarithmic and inverted. The whole Figure 2.19 : 11,, x 11,, 6 [0,1] x [0,1]. The range of d} horizontal and vertical tune space is covered is from 2.2 - 10—10 to 1.9 - 10‘7. 87 analysis and optimization. However, more experience with this technique is needed to make detailed predictions. 2.8 Symplectic Representations The fact that the normal form method yields invariants hinges critically on order n symplecticity of the normal form map. This fact and also the notion that symplecticity guarantees area conservation in every 225-1 x 2;,- plane in phase space suggests to look for completely symplectic transfer maps to describe an accelerator. This can be done by computing an order n generating function for the transfer map. The procedure for this computation was mentioned in the references [Ber88a, Ber91b, Ber90a]. The Taylor expansions of generating functions of symplectic map can be computed with DA based programs, generating functions give implicit defi- nitions of symplectic maps, whereas the nth order Taylor expansion of a symplectic map is only order n symplectic. The implicit definition of the transfer map by a gen- erating function is usually inverted numerically with a Newton method. This method was implemented in COSY INFINITY. As mentioned in [Ber92b], six different modes of symplectic tracking are possible. The map can be represented by the four differ- ent generating functions F1, F2, F3, and F4, or the linear part of the map can be factored out of the map and the mixed variables generating functions F2 or F3 can be used to approximate the nonlinear part of the map by a completely symplectic transformation. We used this procedure to check if it is reasonable to compute pseudo invariants by nonlinear normal form theory and to represent the map by a completely symplectic representation. The analysis revealed that the normal form invariants I ,- of chapter (1) are not better pseudo invariants for a symplectic map generated by either of the six 88 r «.9 'lf‘lh“ .l’" u. “‘13. r .lI,‘ We. ' I'M l ,2 h 5: / I.“ 3': ill“, ‘ ' ' .. ‘ .5 l. i If; Iii-1”“: l Figure 2.20: Variation of the maximum 6 of the deviation function d f with tune for the 8“ order map of the Demo ring. The scale is Iogarithmic and inverted. The range of df is from 1.3 .10’9 to 1.9-10“. 89 generating function. In fact, usually the normal form invariants are better invariants for the Taylor maps than for any of the symplectic representations. We believe that the reason for this behaViour is the order n + l invariance of the functions I,- under transformation of the Taylor map. The relation 11(112) :n+l I} (2.32) only holds for the Taylor map, not for a symplectic map defined by a generating func- tion, since the partial inversion by a Newton method changes all orders of the map, not only orders higher than n. The Taylor map M is the Taylor expansion of a sym- plectic map. Therefore, increasing the order of the expansion makes 11? become closer to a symplectic map. In table (2.6) the invariance defect of the coupled pendulum for the 9‘” order polynomial f is shown for symplectification of the map by evaluating it to different orders. The quality of the invariants does not change substantially, since the first 9 orders cancel with and without this kind of symplectification. Order Coup. Pend. (10‘13) IUCF (10‘9) PSR II (10‘9) Demo (10‘9) 8 5.858924 5.746786 1.597229 8.808577 9 9.721390 5.717355 1.599737 9.759362 10 9.677259 5.721978 1.599741 9.735593 11 9.744150 5.722880 1.599745 9.137200 12 9.745815 5.722909 1.599745 9.867690 Table 2.6: Maximum of the deviation function after symplectifying the map by adding higher orders. In figure (2.21) the deviation function is given for the three example accelerators. For the emittances used, the different generating functions yielded graphs which dif- fered less than printer resolution. Therefore, the generating function F1 was used for all figures. The domain in phase space on which the deviation functions were evaluated is the same as in the figures (215,213,214, and 2.15). In all cases the 90 function f is a better pseudo invariant for the Taylor map than for the symplectic map represented by the generating function. Even for many turns, where satisfying the symplectic condition could be superior, our calculations show that the normal form invariants are better invariants for the Taylor map, as can be seen in table (2.7). C. Pend. (10'14) IUCF (10'11) PSR II (10‘10) Demo (10“3) Turns F1 117 F1 117 F1 A7 F1 A7 10 293 276 881 823 126 115 122581 129261 20 534 508 969 988 232 223 97341 87990 30 723 705 1191 1496 321 310 77279 67182 40 866 840 1091 1396 379 373 21238 36807 50 1002 949 1497 1219 462 449 43855 61792 60 1232 1214 1178 1334 571 561 128 125 70 1476 1444 1427 1028 674 666 162 153 80 1656 1637 1187 934 761 754 172 177 90 1803 1778 859 1134 820 814 194 186 100 1941 1890 1007 721‘ 856 853 209 210 Table 2.7: The maximum of the deviation function as a function of the turns when the map is represented by a generating function. 91 -\ t ,u, a“ fillllllu ““94“!“ . \ .I“ I,‘n\\\l;.xl,\ul‘-}f¢. ’I‘lllllli‘l‘, Aim, " ‘ I“ “m‘ I““l\\l In“ 111103103111: \\. 111m ‘ fills" t‘lu‘l'? m Juli w l ‘ l f, . “I ‘ 111,511," ‘ ‘f‘llllllm' ' . ‘Il‘\\\\\\;‘.I.u“,\‘c:.~, l“ ‘nux. 2‘0. tiny.“ “\I 7“ “1‘“. , t, l\\\\. w". \l '31“ “1‘.'¢|‘||‘\‘|‘\‘,: ’7‘ NH“ ' Will i-; . .‘ ,m 'u ‘\\\h, \u uu‘ ‘n I . ‘ \ . u . ‘l H‘s"... \ . . I I “0:59 “‘1“, m| u \ \ ‘1 \u‘ I‘ t \I‘l , "t 10““! ll, 7‘1““ " ’r'n‘n‘m, . ‘1 1““ i“ . t \ |\\ . \m uh “0'“ .. ’1';:~‘:~ . ~'. ’:.‘.:~ . —_‘--‘ ‘ | N... I u' ‘| 0.. ‘I Figure 2.21: Quality of normal form invariants when the map is represented by a generating function for the four examples with two degrees of freedom: 3.) the cou- pled pendulum, the range of the displayed function is [—2.15,9.45] - 10'”, b) IUCF: [—1.16,1.11]-10‘8, c) PSR II: [—2.63,2.63]-10‘9, d) Demo: [—6.36,8.49]-10‘9. For the order 8 symplectic Taylor maps, the corresponding figures were displayed previously. Chapter 3 optimization by Scanning On page (57) four problems were mentioned. The first problem, finding a suitable pseudo invariant f, has been discussed. The remaining three problems are connected to finding the minimum of f on 8], the maximum of f on 8;, and the maximum of df on A\0. The regions 5,- = M(0)\0 and Sf = M(.A)\A cannot be represented as clearly as the regions 0 = {ila 2 Hill} and A = {ill 2 Hill}. This does not lead to a problem when phase space regions are used which contain 5,- and 8f. When the maximum 50 of the deviation function on O and the maximum 6 on A\0 is known, it is sufficient to choose 5,-=.7:'(a,a+6o), Sf=f(1,l+6) (3.1) with the phase space domain .7(£, C) defined in equation (2.29). The functions f and d; have some properties which allow to make sensible sim- plifications. Those properties will be demonstrated for the proposed PSR II. The evaluation order is 6 and the acceptances are c, = 100mm mrad. As shown in figure (3.1a), the function F(x) = max{f(i)lf(x,x)} is typically growing monotonously with a: so that the maximum of f(i) on 5; occurs on f(a + 6o,a + 60), which is approximately described by Hill = a; and the minimum of f(i) on 8; occurs on f(1,1), where Hill is approximately 1. The function 0(1) 2 max{df(i)|f'(x,a:)} is 92 93 3 AH __1+ 5.0- 10-13 LmullllllllHH _1 — 4.3-10‘l3 a 1.1: $ Figure 3.1: Figure a) depicts F (2:) and figure b) D(:c) in the allowed region. The variation of the pseudo invariant relative to Hill is shown in figure c). also typically growing monotonously as shown in figure (3.1b). Therefore, the maxi- mum 6 occurs also at the border of the allowed region, where Hill is approximately 1. Furthermore, figure (3.1c) shows that the variation of f (i) on S,- and 8f is much smaller than A f = f f — f,- which therefore is close to 1 — a. We obtain the estimate 1—a N = maxidfmlfm» ’ (3.2) which involves finding only one maximum on a subspace with nearly constant Hill. The figure in section (2.6) show the range of the functions d f on the border of the allowed region A for all the example systems. The function df does not have sharp maxima so that sampling with 20 steps in each direction gives a good approximation of the maximum value. Table (3.1) displays N for different systems and for different evaluation orders for a = 1 / 2. Due to energy conservation, the pendulum and the coupled pendulums are stable for all times. The quality of our estimate is shown by the big numbers N which we obtain for those cases, in spite of the fact that the emittances are chosen to be extremely large and nonlinearities quite important. The evaluation of the functions f and df, using interval arithmetic to find the max- imum values and therefore establishing a mathematically strict lower bound for the turn number Nmm, will be analyzed in a later section. First we want to demonstrate 94 Order Pendulum Henon Map Coupled Pendulums 2 434 6 43 3 434 41 915 4 1,039,578 1,109 85,907 5 1,039,578 7,149 2,577,221 6 455,537,706 27,556 61,418,923 7 455,537,706 176,827 1,535,527,685 8 92,114,163,553 1,474,124 29,750,319,370 9 92,114,163,553 9,133,037 357,584,630,384 Order IUCF PSR II Demo 2 9 806 6 3 321 831 129 4 1,288 252,893 1,220 5 19,995 235,650 25,657 6 370,294 6,977,545 84,087 7 3,265,268 8,255,710 1,320,751 8 11,277,884 65,472,668 4,554,994 9 65,734,218 76,092,850 55,548,695 Table 3.1: Minimum number of turns required to move from the initial region to the forbidden region. The initial emittance was chosen to be half the acceptance. The maps were evaluated in order eight, whereas the pseudo invariants were com- puted to the indicated order. Scanning was performed with 20" points for 1: relevant dimensions. 95 some possibilities of improving the obtained estimates. 3.1 Refinements and Examples Two methods will be introduced that increase the quality of the bounds on long term stability. One is connected to separating phase space in appropriate regions, the other involves multi—turn maps. Furthermore, unknown parameters of the system will be included in the estimates. The results for the six example maps are accumulated in table (3.3). In this table the guaranteed number of turns for the assumed acceptance c1r, which is given in section (2.4) for each example, and an initial beam with emittance cvr / 2 is given. For the coupled pendulum and for the IUCF ring, the emittance for which a beam can be guaranteed to survive 108 turns is also shown. To limit the computation time, we evaluated 10“ points for k relevant phase space dimensions. The parameter interval was also covered by 10 points. 3.1.1 Dividing Phase Space Because of the rapid increase of 6 with increasing Hill, it is appropriate to separate the regions between Hill = a and Hill = l by surfaces f(a,,a,-), i E {0,...,k}, where do = a and a}, = 1. With u,- = max{f(i)li E f(ag,a,-)} and l,- = min{f(i)li E, .77 ((1,, 05)}, the lower bound on the number of turns becomes 1‘ . _ . N = 2 L351 with 6,- : max{df(i)lf(a,-,a,-)} . (3.3) i=1 i The turn numbers obtained from this technique and the transportable emittances are given in table (3.3) in the third line. In our experience, this separation of phase space can improve the estimate by up to a factor of 10. The potential of this approach can be seen when using D(:r) = max{d(i)lf(i) = 9:} . (3.4) 96 We estimate the change that occurs in x by k map applications as 1-(n + k) — x(n) = D(:::(n + k)) - k and approximate for big turn numbers N as M.) = 0. 5% . (3.5) In figure (3.2) n(:1:) corresponds to the area under the curve. The PIE method without dividing phase space, however, only gives the area in the rectangle as guaranteed survival time. Am 2.26 - 109 1.18-108 Figure 3.2: Dividing phase space approximates the area under the displayed curve by several rectangles to obtain the bound on the survival time. Without this improve- ment, the PIE method would only give the area of the drawn rectangle. 3.1 .2 Time Evaluation Figure (3.3a) shows how close the pseudo invariants are to invariants of motion. Over many turns the quantity of f changes but it always oscillates around a mean which grows very slowly. The normal form method, however, bounds the growth of the pseudo invariant by considering the largest growth that can happen in one application of the map. As demonstrated in figure (3.3b), the biggest change that can 97 Turns 9500 A 5-10-9 WW” —2 - 10-9 Turns 0 Figure 3.3: a) slow increase of the pseudo invariant over many turns, b) periodic change of the pseudo invariant after relatively few turns. be generated by N applications of the map is usually much smaller than N times the largest growth that can happen during a single turn. Therefore, it is advantageous to consider the maximum growth that can occur when the map is applied several times. In order to find the optimum number of map applications, the maximum and minimum of d(i) over the allowed region A was plotted for many turns. The fig- ures (3.4a—f) display the variation of d(i) for our examples. For the improvements illustrated in table (3.3), the number of map applications used are displayed in table (3.2). [l Pend. Henon Coup. Pend. IUCF PSR III DEMO |[ 14 37 6 89 247 100 Table 3.2: Number of map applications used for improving the estimation. 3.1.3 Parameter Dependence So far, the normal form method assumes that the one—turn map of the storage ring in question is well known. Since this is rarely the case, the theory has been extended 98 Ale-12 Ale-8 WW 0 WW” "18.47 0 100 > Ale-11 A 10-8 2.29 o —2.13 0 200 {> o 300 "> 0 150 {> Figure 3.4: The figures display the variation of the deviation function in the allowed region for the six examples. From top left to bottom right these are a) Pendulum, b) Henon map, c) Coupled pendulum, d) IUCF, e) PSR 11, f) Demo. 99 to maps which depend on an unknown parameter. Neither particle energy nor the magnet parameters are known exactly and have to be treated as parameters which can not be accurately specified. We computed the map as a function of a parameter of interest and then scanned not only through the interesting region of phase space, but also through the interval in which the parameter could lie in reality. The unknown parameter for each example is mentioned in seétion (2.4), where the examples are introduced. If the map NI is a function of a parameter, then the nonlinear normal form trans- formation 8 also depends on a parameter. Using DA programs allows, as mentioned in chapter (1), to compute the Taylor expansion of this parameter dependent map 3: The pseudo invariant f then also depends on the parameter. Changing the parameter changes the transfer map and the pseudo invariant simultaneously, such that f stays a good pseudo invariant for a wide range of the parameter. This is apparent when one considers that f0 M — f :n+1 0 (3.6) also holds for parameter dependent normal form transformations, where partial deriva- tives in respect to the parameter and in respect to the coordinates are considered when equating with “=,,”. We belive that the described method usually gives a very reliable lower bound on the number of stable turns. However, to make completely rigorous statements about guaranteed bounds, several steps have to be performed in a more rigorous way. In particular, it will not be sufficient to approximate the maximum of 6 by scanning. Interval arithmetic methods which guarantee a global maximum have to be used. Utilizing this guaranteed optimization requires a short introduction into interval arithmetic. 100 Lower bound for Pendulum Henon map Simplest application 92,114,163,553 9,133,037 Length and Strength uncertainty 1% 23,556,300,993 8,167,533 Divided phase space 452,868,876,965 44,999,781 Multi—turn maps 388,804,862,856 1,456,171,297 Both 1,895,348,117,634 4,779,711,057 Lower bound for PSR II Demo Simplest application 76,092,850 55,548,695 Quad field uncertainty 0.01% 47,166,060 51,963,620 Divided phase space 373,642,327 284,008,517 Multi—turn maps 21,172,838,624 1,042,575,616 Both 18,731,455,785 5,121,716,506 Coup. Pendulum predicted turns stable emittance Simplest application 357,584,630,384 5.00 x 5.001rm mrad Length uncertainty 1% 144,173,434,143 3.40 X 3.407rm mrad Divided phase space 1,765,031 ,547,898 11.3 x 11.37rm mrad Multi—turn maps 1,029,815,934,687 6.00 X 6.007rm mrad Both 5,087,629,041,331 14.5 x 14.51rm mrad IUCF ring predicted turns stable emittance Simplest application 65,734,218 3.3 X 2.07rmm mrad Quad field uncertainty 0.01% 18,535,102 2.9 x 1.77rmm mrad Divided phase space 335,420,083 4.8 x 2.97rmm mrad Multi—turn maps 5,804,832,818 8.1 x 4.87rmm mrad Both 27,745,480,680 13 x 7.71rmm mrad Table 3.3: Lower bounds on the turns of particles for initial emittance of one half the acceptance and lower bounds on the stable emittances for 108 turns obtained by various variations of the PIE method. Chapter 4 Rigorous PIE with Interval Arithmetic The introduced method of pseudo invariant optimization has the potential to be com- pletely rigorous. Equation (2.4) is a strict inequality, if all the maxima and minima can be strictly bounded. In the field of numerical analysis with automatic result ver- ification, a number of methods have been developed which lead to rigorous bounds on the extrema of functions. Below, the basics of interval arithmetic will be intro- duced and it will be explained how interval arithmetic can lead to rigorous bounds on extrema. In preparation of the presented work, several methods of standard inter- val optimization have been used. However, for the complicated function of; involved, which is a multidimensional polynomial of order n(n + 1) when n is the evaluation order of the normal form transformation, conventional techniques turned out to be fair to slow to be useful at current computer speed. Therefore, these standard meth- ods will only be shortly mentioned and the reader is referred to the references. In the weakly nonlinear systems of interest, the inefficiency of the conventional methods is mainly due to the low order contributions in the polynomial df. In order to take advantage of the special structure exhibited by the deviation function, namely 4 df=f(A/)—f=n+10, (4.1) 101 102 we introduced an innovative concept which we shall call interval chains. Interval chains allow one to bound the contribution to a polynomial of order higher than n + 1. Using interval chains therefore completely eliminates the inefficiencies caused by low order contributions. After this concept is introduced, several examples will be presented. The example maps used are the same as those of the previous chapter and the refinements used correspond to those in section (3.1). 4.1 Interval Arithmetic Interval arithmetic was developed because of the fact that the set of numbers which can be expressed on a computer is finite. Therefore, computers cannot store accurate information about a real number, except that it is located within an interval of width equal to or greater than the computer accuracy. Thus, if completely accurate arith- metic should be performed on a computer, an arithmetic for intervals has to be de- veloped. Such an arithmetic is also necessary to formulate accurate statements using results of measurements, because measurements usually do not yield real numbers but intervals, often denoted by error bars [AH83, Ku189, May89, M0066, M0079, M0088]. Following, intervals will be symbolized by capitol letters and real numbers will be denoted by small letters. The upper and the lower bound of an interval will be denoted by underscore and overscore respectively, X = [LY] = {xIX S x S Y} - (4-2) If f : B H B then f (X ) is defined as the set of possible results for all elements of the interval X, f(X) = {ala = f(x),:1: E X}. (4.3) It will be our goal to describe a method to find an interval which contains all possible results. Such an interval that encloses f (X ) will be denoted by F (X ) 2 f (X ) 103 Elementary Operations If we are concerned with functions which can be constructed by a finite number of elementary operations, it is sufficient to ensure that the result of every elementary operation on intervals is an interval containing all possible results. If this is true for any elementary operation, it will infallibly be true for the result of the complete function evaluation. This thought leads us to the need of establishing the result of elementary operations between intervals. The required property of an elementary operation between n numbers, described by the function e : R" H B, is given by E(X1,...,X,,) 2e(X1,...,X,,). (4.4) In order to achieve results F (X ) which approximate f (X ) as closely as possible, one tries to define the elementary operations such that the equality holds in equation (4.4). At first we choose addition, subtraction, multiplication, and division as elementary operations: X+Y = [x+i:,7+7'], (4.5) X — Y = [Kim + [‘73 ‘31] i (46) X - Y = [min{x7, 317,771 ,77}, max{Xlfi “117,11 , Y}] , (4.7) X/Y = [_X,mo[1/17,l/leor0¢l’. (4.8) The addition and the multiplication are commutative and associative. It is important to note that intervals with addition do not form a group, since there is no inverse of addition, which follows immediately from equation (4.10). In particular, —X :- [—X, —X] is not the additive inverse of X, because X — X = [K - X,X - X] 94 0. There is also no distributivity law, but subdistributivity X-(Y+Z)§X~Y+X-Z. (4.9) 104 The diameter of an interval is defined as d(X) = X — 2g, and the absolute value is defined as IX l = max{l_)£ I, le}. It is helpful to know the diameter and absolute value of results of elementary operations d(X:l:Y) = d(X)+d(Y), (4.10) d(X-Y) g d(X)-|Y|+d(Y)-|X|, (4.11) |X+Y| g |X|+|Y|, (4.12) lX-Yl = IXI'IYI- (4-13) If more complicated functions have to be evaluated, which cannot be established by a finite number of those four operations, then further elementary functions have to be defined. Examples are sine and cosine functions of intervals. Interval Blow—up The function f: x H x-x yields f(X) = 0 and F(X) = X —X = [X—X,X—X_]. This readily establishes that interval computations, albeit yielding all possible results, can overestimate the set of possible results substantially. In this example the set of possible results {0} is contained in F (X ), which is however a great overestimation. The blow—up occurs because the constructed arithmetic treats the two intervals X to the right and to the left of the operator - as two unrelated intervals. Another example is f : a: H :c*a:. The results for an interval X = [-—a, a] can be in the interval f (X ) = [01, 02], whereas the arithmetic yields a bigger interval F (X ) = [-a2, a2]. There are three major ways to minimize blow—up: i) Reduction of the number of elementary operations. This can be achieved by restructuring the evaluation of the function. For example, f (3:) = a: — a: can be restructured to f(x) = 0. Because of the subdistributive property of interval 105 arithmetic (4.9), a reduction of the number of multiplications by restructuring parenthesis will always be useful. 2 ii) Introduction of new elementary operations. For example, 8(1)) = a: can be introduced as new elementary operation which is evaluated such that E( X) = e(X) by E(X) = [max{0, (X - X)}, (X - X)]. iii) Decreasing the absolute value of intervals involved in a multiplication decreases the resulting diameter according to equation (4.11). Therefore, it is sometimes helpful to use a so—called centered form, which can be obtained from the inter- mediate value theorem. The intermediate value theorem states for differentiable functions on X: Vac E X 3 5 between a: and c such that f(x) = f(c) + f’(§) - (a: - c) . (4.14) Any value f(x) with a: E X is in the interval f(c) + F’(X) - (:1: — c). Therefore, we can write the centered form F C(X ) 2 f (X) as ‘FC(X)=f(c)+F'(X)-(X—c) (4.15) and have achieved a multiplication with an interval X —c, which has the smallest absolute value when c is chosen to be the center of X. The center of an interval X is defined as m(X) = (X+ _X_)/2. The possibility of computing an interval F(X) 3 f(X) (4.16) in a straight forward way for complicated functions f implies the possibility of finding upper bounds on f in the interval X. F(X) is such an upper bound. If F(X) overestimates the maximum of f on X substantially, the interval X can be subdivided 106 in several subintervals X,- with UgX,‘ = X and a better bound can be found by f(x) S fm = max(F(X,-)) . (4.17) When the width of X,- is made small, the bound will become tight, due to decreased blow—up. This interval rastering method automatically gives a measure for the over- estimation fm — f S max(F(X,-)) — max(E_(X,-)) . (4.18) The series of graphs in figure (4.1) shows how the blow—up reduces when the region of interest is covered with smaller intervals. For further refinement of global optimization with intervals see [Jan91, JK92, Jan92a, Jan92b, Han79, Han80, RR88, WH885, Rat92, Cse91, Eri91, Han88, IF79, MR88, Moh90]. One important refinement is the successive decrease of the width of the subintervals X, after exclusion of subintervals X k which certainly do not contain the maximum due to F(Xk) < max(£(X,-)) . (4.19) Conventional local optimization can find intervals X,- with big E(X,-), which will effectively exclude many intervals X k by criterion (4.19). 4.2 Interval Chains and Optimization The deviation function 4 .7 _. df=f R can be written as p : p(i) = Z_0 p,-(x ). For i greater than m, p,- is chosen to be zero. Definition 4.2.2 Call an interval chain P(Ak) = (P0,...,P,,+1) an interval chain (IC) evaluation ofa polynomial p of order m on. the interval boa: A’c = A1 x x A), ifP, 2 {p,(iHi E Ak}, 0 S i S n and Pn+12{2:;n+1p,(i)lie Ak}. Thus Pu.” bounds all contributions to the polynomial of order higher than n. Theorem 4.2.3 If F(Ak) is an [0 evaluation of a polynomial f on A" and C(A") is an [0 evaluation ofa polynomial g on A", then F(Ak) + C(Ak), F(A") - C(A"), and rF(A"), r E R are [C evaluations off + g, f - g, and rf on Al“, respectively. Proof: Given the polynomials f : f(x )= 210 (i) and g :g(:i:') = 26:0 g,(i’), then by definition Fi2ffi(f)|5€A"} . Gi2{9i(5)|564k} OSiSn. Fn+1 2i 2 fiff )ll‘ 6 Ak} aGn+1 2f 2 gifl‘ 5H3? 6 Ak} (4-23) i=n+1 i=n+1 from this we infer (F(A") + G(A")).- = F.- + G.- 2 {(f +9).(&:‘)l5:' e Ak}, 0 s 2' s n +1. (rF(Ak)),- = rF, 3{(rf) (1:)l €,Ak} 0SiSn+1, (F(A*)-G(A")).- = ZFG.-.;>{ZI. x)g.-_.-( (”)lxeAk} -- {(fg).-(:2')I:i’ e Ak}. 0 s 2' s n. n+1 n+1 (FfAkl’G(Ak))n+1 = 203‘ Z Gj) i=0 j=n+1-i 110 = :[m 2 G,+Gn+1)]+Fn+1(zn:Gj +Gn+1) i=0 j=n+1-i J=0 2 [21]“ Z gj+ Z ngl+ 2 fi) (:)gj‘l' :9 gill )llAk) i=0 j=n+1-—i j=n+l i: n+1 j: j=n+1 a 3 = [20.- 2 g.)1(A*) i=0 j=max(n+l-i,0) 0+5 0+5 = [2(Zf.-.g.)1(A")={Z()fg-(a'i)la:EA’°} (4.24) i=n+l J=0 i=n+l The expression in the fourth line from the bottom contains the expression in the next line, since it is an interval evaluation of the function in the third line from the bottom. Some notations will become useful: Let I be an interval, then [I ; k] is an interval chain with the k“ order being I and all other orders being 0. We will only use [I;n+1] =(0,...,0,I) and [1; 1] = (0,1,0...,0). (4.25) Given an IC by [c 2 (IO, II, . . . , In“), then the interval [10])‘ = I), is the kth order of 1,. Theorem 4.2.4 Let now A be any algorithm based on addition, scalar multiplication, and multiplication that evaluates a polynomial p at i = (x1, . . . , wk). Then performing A on 1'; = [(0,A,,0, . . . ,0),. (.0 A.,0, . 0)]: M"; 1], (4.26) a vector of interval chains, yields an IC evaluation ofp on 11". Proof: For all I E {1, . . . , k}, the interval chain (0, 141,0, . . . ,0) is an IC evaluation of the identity polynomial i1 : i;(i) = x1 and P(A") is the result of finitely many elementary operations which evaluate p and therefore, using theorem (4.2.3), we infer that P(A") is an IC evaluation of p on A". 111 Bearing this in mind, it is easy to see that for an IC evaluation D(Ak) of the deviation function d; on A" it holds that {d;(i)|i E A"} Q Du“, since df is known to have no contributions with orders lower than n + 1 if the map and the normal form transformation are evaluated to order n — 1. In this approach cancellations up to order it do not contribute and blow—up caused by such cancellations is completely avoided. In fact, subtraction of f is not necessary at all any more, since D... = [4112)]... = [f o 47(4))... —U(f.)1..4 = [f 0 Main... , (4.27) with the polynomial f of order n. But even with these simplifications, the resulting objective functions have a ten- dency to exhibit interval blow-up because of complexity, while the bounds of the function have to be determined rather tightly in order to guarantee many turns. F ur- thermore, the functions have a large number of local maxima. All these effects make the exclusion of intervals rather difficult and not practically possible unless in the or- der of 10“ intervals per dimension are used. These large numbers make bookkeeping of intervals for later exclusion rather cumbersome if not impossible. Because of this situation, the conventional methods discussed in subsection (4.1), which are based on disposing of intervals that can be excluded and halving of the remaining ones combined with occasional local optimization in real arithmetic, are not directly applicable. For this reason, we restricted ourselves to a mere rastering of the objective function with a large number of intervals of equal size. The results in the next section were obtained by choosing 630 intervals for the examples with one degree of freedom and 1,000,188 for the example with two degrees of freedom. Because of limitations in computation time, it was still assumed that the maximum of the deviation function occurs on the surface of the allowed region. Without the concept of interval chains a realization of the described method was virtually impossible. For 112 examples with two degrees of freedom approximately 1012 times more intervals would have been needed for similar results. 4.3 Parameterizing Regions for the IC—PIE Method To avoid interval blow—up, it is advisable to perform as few operations as possible. 7 This fact becomes obvious when subdistributivity in equation (4.9) is considered. We therefore try to minimize the computations required to represent the initial and the allowed region. Previously the invariant tori of order n were represented by equation (2.28). The conventional emittance and acceptance of linear motion is given by the invariant tori of first order motion, which are given by n = 1 in equation (2.28), and by construction of the sets f({, C) in equation (2.29). We now write the tori in linear normal form space as 7mg, 2 {£1 ( 2..-. ) = fi( $58)) ),q), e [0,24],v.'e {1,...,d}} (4.28) ~ . ‘2: and the regions in linear normal form space as rims C) = {212' <—: mom 5 )3 —" s C} . (429) Then the relevant phase space regions are given by HM) = {B'r‘wnze mam} . (4.30) where B," is the linear part of the inverse normal form transformation B“. To find the maximum deviation function 6 on the surface of the acceptance, we have to find 6 = max{(fo)v7—f)oB',-‘|f~p(1,1)} (4.31) = max{foBfloB1oMoBfl—foBfIIFNFUJH. 113 The map El 0 NI 0 Bfl is the transfer map in the first order normal form space and is to first order a rotation in every 22,_1 x 2;,- plane in phase space. Call this linear rotation map If and let N = Bl 0 NI 0 Bfl o R"; N is the identity to first order. Write g = f o B," to get 6: max{goNoR—glfnp(1,1)} . (4.32) g is a polynomial of order n if the evaluation order for [V7 is n - 1. g therefore does not contribute to the order n + 1 of the interval chain which is used to bound d; and the second appearance of g in equation (4.32) can be omitted. The rotation If leaves the tori invariant and can therefore be avoided; we finally get the simple interval evaluation of 6 = max{g o NIpr(1,1)} . (4.33) Now the intervals [0, 2n] in the definition of TNF in equation (4.28) have to be covered by many small intervals. The maximum upper bound of the intervals [g o N I Icln+1 has to be found for all interval chains 16 that cover fNF(1,1). This yields a rigorous upper bound for 6. All blow—up due to linear transformations is avoided in this way. 4.4 Comparison Between Intervals and Interval Chains Several nonlinear systems were studied using the interval chain rastering methods to provide upper bounds for the invariant defects. In order to get a sense for the quality of these upper bounds, the numbers were compared with approximations for the max- imal invariant defects obtained by a rather tight rastering in real arithmetic. Because of the large number of local maxima, this method proved to be the most robust non— interval approach to estimate the absolute maxima of the functions involved. Lower 114 Order of Interval Interval Conventional Invariant Bounding Chains Rastering (guaranteed) (guaranteed) (optimistic) 3 11,252 743,667 849,195 4 11,252 743,667 849,195 5 11,306 876,059,284 982,129,435 6 11,306 876,059,284 982,129,435 7 11,306 432,158,877,713 636,501,641,854 8 11,306 432,158,877,713 636,501,641 ,854 Table 4.1: Predictions of the number of stable turns as a function of the order of the polynomials describing the normal form transformation for the physical pendulum with a maximum elongation of 1/ 10 rad. Because of energy conservation, the map is known to be permanently stable. bounds on the number of stable turns obtained by conventional intervals are given in the tables (4.1), (4.2), and (4.3) in order to illustrate the usefulness of interval chains. When conventional intervals were used, the deviation function was simplified as much as possible by accounting for cancellations up to second order analytically. The number of conventional intervals and the number of interval chains used in the bounding are equivalent. As the first example to check the method, we used a one—dimensional physical pendulum. This is a good test case, since energy conservation requires the nonlinear motion to be stable. Table (4.1) shows the results of the stability analysis for this case. As is to be expected, the number of stable turns predicted increases with the order and hence accuracy of the approximate invariants. While the approximate scanning method can take full advantage of this increased accuracy, the interval bounding method shows a saturation at 11, 306 turns. This asymptotic behavior is connected to the size of the intervals because of the unavoidable blow—up of intervals in the process of cancellation of large terms. The blow-up in third order dominates the calculation, causing the higher order improvements to not materialize. The method of interval 115 chains takes care of all the cancellations up to evaluation order and consequently the estimate is much better. Order of Interval Interval Conventional Invariant Bounding Chains Rastering (guaranteed) (guaranteed) (optimistic) 2 895 891 1,086 3 1,736 9,926 11,450 4 1,668 54,016 65,667 5 1,674 678,725 809,612 6 1,670 3,389,641 4,351,679 7 1,671 42,640,927 52,474,387 8 1,671 192,650,961 263,904,035 Table 4.2: Predictions of the number of stable turns for the Henon map at tune 0.13, strength parameter 1.1, and starting position (:r,a) = (0.01,0) as a function of the order of the normal form transformation. As another example, we chose the Henon map described in section (2.4). The results of these calculations are shown in table (4.2). Similar to the previous case, the number of predicted turns increases with order. In the case of interval bounding, the number of guaranteed turns shows asymptotic behavior limited by blow-up. Again the superiority of strict bounding with interval chains is obvious. Order of Interval Interval Conventional Invariant Bounding Chains Rastering (guaranteed) (guaranteed) (optimistic) 3 179 16,137 38,385 4 179 18,197 38,857 5 173 309,356 560,309 6 173 347,312 613,135 7 171 925,531 2,184,998 8 171 1,004,387 2,248,621 Table 4.3: Predictions of the number of stable turns as a function of the order of the approximate invariant for the Los Alamos PSR II storage ring for the motion in a phase space of 100 mm mrad. In the final example, we study a realistic accelerator, the Los Alamos PSR II 116 already described above. The same data are shown as for the two previous, more aca- demic examples. To limit the calculation time, the intervals used for the optimization were 5 times as wide as the intervals used for the previous two tables. The results are shown in table (4.3). 4.5 Refinement of the Rigorous Estimates In section (3.1) the long term estimation by pseudo invariants was improved by three methods: 1. analysis of multi-turn maps, 2. separation of phase space, 3. consideration of random parameters. The improvements obtained by separating phase space in smaller regions can increase the bound on the survival time by about a factor of ten. With interval arithmetic this method becomes impracticable. As observed, using arithmetic with interval chains, one can bound the maximum of the deviation function d f rigorously and rather tightly. This is only possible, since all low order contributions to the deviation function vanish. When separate phase space regions are used, one has to find several maxima and minima u,- = max{f(i)|iE f(a,,a,-)} and l,- = min{f(i)|iE f(a,,a,-)} as mentioned in equation (3.3). The polynomial f does not have the property that low orders vanish and therefore the advantages of ICs do not materialize. At current computer speed, these maxima and minima can therefore not be bounded rigorously with the required tightness. In order to minimize interval blow—up, it is desirable to minimize the number of operations needed. For the rigorous computation with intervals, we chose therefore 117 the first order acceptances to describe the allowed and the initial region, and we performed all first order transformations analytically without blow-up, as already described in section (4.3). Therefore n = 1 was used in equation (2.28) when representing .7: (C , C ) by equation (2.30). For optimization by scanning, it was necessary to use the nonlinear acceptance to make I,- — u,_1 as big as possible when separation of phase space was used. Because of the described problems with finding extrema of f, it is not possible at current computer speed to find the required maximum on S, and minimum on S I rigorously. In the examples given below we performed these optimizations by scanning. The more critical deviation function d I however was bounded rigorously. As mentioned before, the introduction of interval chains improved the computation time at least by a factor of 1012. However, even after this reduction of CPU time, at current computer speed the needed computations take very long. In the examples a restriction to the surface of the allowed region was necessary when df was optimized with interval chains. One can try to use multi—turn maps and random parameters to refine the IC- PIE method. For the simple application of PIE, the surface of the allowed region was completely covered in order to rigorously bound the maximum of the deviation function. In the case of two dimensional phase space, the surface was covered by 800 intervals; for four dimensional phase space, 8,000,000 intervals were used. For the parameter dependent case, these numbers were 160,000 and 50,000,000, respectively. 4.5.1 Analysis of Multi—Turn Maps It was discussed in section (3.1) and displayed in figure (3.3) that the quality of the pseudo invariants fluctuates depending on how often the transfer map is applied. One 118 can try to use multi—turn maps for rigorous estimation with the use of ICs. Usually the map has to be applied 16 times for every phase space interval, in order to evaluate the k—turn map. Computing the nth order Taylor approximation of the k— turn map would only take approximately [log2(k)] Taylor map compositions and would save a lot of CPU time, however, using this approach would not be rigorous, since the orders higher than the evaluation order n, which are created by the compositions, would be neglected. Interval chains allow one to take advantage of the speed increase achieved by Taylor map compositions. This is due to the fact that we can bound the high order contributions of map compositions by applying ICs. In the subsequent section (4.6), Taylor maps with remainder intervals will be used to describe non—polynomial maps rigorously. For this we denote the map by 114 and its nth order Taylor map by {M}n. The remainder only has contributions of orders higher than n. An interval bound on the remainder of {ll-4}" over a phase space region A will be denoted by the vector of intervals R1(A), and then _9 q 114(2) 6 {M},,(E') + 61(4) Via A. (4.34) With interval chains, a remainder of the two—turn map 472 can be computed if a region Bl of phase space can be found with 81 2 A7 (A) One way to find such a region is the interval evaluation of 37,, by -0 31 = {M}n(A) + 3104) 2 MM) - (435) Let R1(Bl) be the remainder of the Taylor map in the region B]. Then the Taylor expansion of the two—turn map is bounded by m M212) {M}.({M}.(2) + 6.04)) + 6.03.) {lenIE'l + [{MInIIMInIl‘li 1l) + [Rid/1); n '1' ll)ln+1 + RN31) (1) (WM?) + 62(4), (4.36) 119 and one is left with an interval remainder for the nth order expansion of the two— turn map. [B1(A);n + 1] could only be written because the remainder only has contributions of order higher than n. The necessity to know a bound on the region NRA) in order to compute B204) can make computations very involved. This becomes apparent when a remainder interval for the 2k—turn map is computed from a given 2k‘1-turn map with remainder interval. Let equation (4.34) be given and assume that from this information a bound on the 2""1 turn map was computed with ICs to get M2*"(5) e {M2*“}.,(2) + R2444) vz‘e A . (4.37) The demonstrated method computes the remainder of the 2k—turn map by RNA) = [UV-I?"-1 }n({M2k-l}n(lA; ll) + [Eh—1(4); n +lllln+1+ Riv-4334) (4-38) with the region 324—1 2 1172“] (A), which can again be found by interval arithmetic evaluation of {ll-42k-l },,(A) + Rye—l (A). However, finding the remainder B21.-1(B21_1) requires to repeat all the procedure of obtaining bounds of compositions from M up to 1172“], now for 821—1 instead of A. This general problem can be overcome easily in case of the PIE method. The reason for the simplifications is that the simple application of the PIE method for the one—turn map gives bounds on all the needed regions 8,. If an interval remainder of the 2'—turn map is desired, then let B = f(0,7) and let the maximum 6 of the deviation function in this region be given. Choose 7 such that A = f(0,7 - 2l - 6) is the phase space region for which remainder intervals are to be found. B therefore contains all regions M‘(A) for i E {1, . . . ,2’}. We can compute all remainders over that region B to get bounds for map applications in the region A. The step from the 2"”1—turn map to the 2k—turn map simply reads 2172*(2) e {W}. (4.39) 120 + 1102*" 1.0474“ }.([A; 1)) + [6.144183 n +11)).+.+ 14,.-.(3) . There is now no need to go back to a lower turn number for obtaining bounds on a different region than B. All the required interval chain evaluations can be performed on a set of smaller intervals covering B. As before, the union of all remainders found on the subintervals is then used as an interval remainder. The size of remainder intervals decreases drastically with the evaluation order used. This is shown in table (4.4), where the union of the remainder intervals for the different phase space coordinates is displayed. To allow comparison to the previously demonstrated calculations, we continue to use order 8. In table (4.5) the bounds on the remainder of the multi-turn maps for various turn numbers is displayed for four examples. For the Pendulum and the Henon map, 700 intervals were used to cover the region of interest. For the coupled pendulum, the number of intervals was 1,000. For the PSR II ring, 1,000,000 intervals were used to reduce blow—up. The deviation function d; is bounded by d;(z“) e [f({M2"}.([A;1])+[R’;n+11)]..4.+{f(I)72")}.(z*) (4.40) = [f({Mzk}.([A;1])+[R';n+10].... (4.41) The term { f (ll-4' 2" )}n vanishes, since f is an order n invariant of M as well as of 11772“. The higher order Taylor coefficients of the composed map 1172’“ increase with k, and with the increased nonlinearities also the interval blow—up increases. This blow— up counteracts the advantages of applying multi—turn maps, which do not materialize with the current evaluations of intervals and at the evaluation order used. Increasing the evaluation order substantially, however, is not possible at current computer speed. Table 4.4: 121 Order Pendulum Henon 6 902880.288 - 10--5 90.450,0.483j - 10-4 7 904550455 -10-'7 90.997,0.724j -10-5 8 90.456, 0.456 - 10-7 90.189,0.199j -10-5 9 901120.112 - 10-8 90.600,0.562j .10-6 10 901120113 .10-8 901820.179 -10-6 11 90.347,0.347j -10-lo 90.450,0.500j -10-7 12 90.350,0.350j .10-10 90.145,0.140j -10-7 13 90.121,0.121j .10-ll 90.433.0462 -10-8 14 901220.122 .10-11 901420.150 -10-8 15 904220.422 .10-13 90.481,0.461j .10-9 16 904250.425) .10-13 901520.153 -10-9 Order Coup. Pend. PSR II ring 6 —0. 401, 0. 401j 10-5 —0. 451, 0. 451j 10-5 7 90 499, 0. 499 - 10-6 90. 370, 0. 370 - 10-6 8 901120113 -10-6 903520.358 - 10-6 9 901520152 - 10-7 90.289,0.289j -10-7 10 903320.332 - 10-8 902780.278 -10--7 11 904820.482 - 10--9 90.217, 0.217j .10-8 12 90.982,0.982j -10-1° 90.207,0.207j -10-8 Decreasing width of the remainder interval with evaluation order. 122 Turns Pend. Henon 2 :—0.659,0.659: - 10‘11 0.000, 0.000: 4 :—0.354,0.354: . 10’10 :—0.128, 0.114: - 10'8 8 :—0.261,0.261: - 10‘9 :—0.732,0.870: - 10'7 16 :—0.353,0.353: - 10‘8 :—0.296,0.339: - 10‘6 32 :—0.457,0.457: - 10‘7 :—0.190,0.199: - 10‘5 64 :—0.791,0.791: - 10‘6 :-0.221,0.239: - 10'4 128 :—0.124,0.124: - 10“ :—0.214,0.216: - 10‘3 256 :—0.218,0.218: - 10‘3 :—0.361,0.320: - 10‘2 Turns Coup. Pend. PSR II ring 2 :-0.247,0.247: - 10‘11 :-0.367,0.367: - 10‘11 4 90.108,0.108j -10-1° 90.144,0.144j -10-'9 8 :—0.509,0.517: - 10"10 :—0.212,0.212: - 10‘8 16 903720.374} . 10-9 902420240 -10-7 32 :—0.449,0.449: - 10’8 :—-0.44l,0.441: - 10"6 64 [—0.586,0.592: - 10"7 :-0.789,0.789: - 10‘5 128 [-—0.705,0.705: - 10‘6 :-0.932,0.931: -10“ 256 :—0.117,0.115: - 10“ :—0.880,0.880: - 10'2 Table 4.5: Remainder intervals for multi—turn maps. 123 4.5.2 Random Parameters We want to find a lower bound on the survival time for a given accelerator including uncertainty of a parameter in a given tolerance. The same settings as in table (4.6) were used. Instead of the length or field strength, also other parameters could have been used, like the uncertainty of the particle’s energy or uncertainty in the length of an element. The pseudo invariants are computed by parameter dependent nonlinear normal form theory, which is described at the end of chapter (1). While a different parameter produced a different transfer map, it also produces a different pseudo invariant which changes very little during one application of that map. This is due to the fact that f o M - f =n+1 0 (4.42) holds for parameter dependent map, when “=,,” equates partial derivatives in respect to parameters as well as in respect to coordinates, and f (5,6) as well as 1136,45) are functions of the parameter 6. Lower bound for Pendulum Henon map Simplest application 52,366,096,777 5,830,904 Length and Strength uncertainty 1% 14,051,907,204 3,123,964 Lower bound for - Coup. Pendulum IUCF Simplest application 56,917,938,176 19,500,358 Length/ Field uncertainty 1% / 0.01% 27,080,626,416 10,768,020 Lower bound for PSR II Demo Simplest application 58,680,622 10,106,151 Length/ Field uncertainty 1% / 0.01% 45,819,009 4,831,120 Table 4.6: Lower bounds on the turns of particles for an initial emittance of one half the acceptance, obtained by the IC—PIE method. 124 4.6 Using Taylor Maps with Remainder Bound So far we were only concerned with nonlinear motion which is described by Taylor maps. The number of interest was the survival time of particles in an accelerator. This time was formulated as the number of map applications for which no phase space point of the initial beam distribution is mapped- into a forbidden region. A method was presented with which rigorous lower bounds on this number can be obtained. In the phase space regions which we analyzed, the Taylor maps usually describe the accelerator well and the limits obtained are valuable for storage rings, however, strictly speaking they are not rigorous, since the transfer map of an accelerator is not a polynomial map. This problem could be overcome if a bound on the remainder of the Taylor map would be known. The method of RDA, described in reference [BH94a], potentially provides a way of rigorously bounding the Taylor remainder of functions which can be evaluated with DA based programs. As on page (108) for ICs, also for this method elementary operations are introduced such that an operation between two Taylor maps with remainder intervals yields again a Taylor map with remainder interval. Furthermore, function evaluations on Taylor maps with remainders are introduced, as well as a derivative. With these operations, all functions which can be evaluated with DA programs to obtain their Taylor map, can also be evaluated with RDA programs to obtain their Taylor map and a bound on the remainder. The implementation of the rigorous RDA method is a bigger project and will be future work. Here we will estimate remainders by comparing the 8th order map with the 19“ order map for the two dimensional case and with the 12‘“ order map for the cases of four dimensional phase space. The results for the examples are shown in table (4.7). The long term analysis will be performed with the Taylor map and the remainder 125 Pendulum Henon Coup. Pend. [—0.5l7,0.517] - 10‘12 [ 0000,0000] [—0.268,0.268] ' 10‘9 IUCF ring PSR II Demo ring [—0.248.0.248] - 10“l3 [—0.282,0.282] - 10“ll [—0.727,0.727] - 10'8 Table 4.7: Maximum error of the Taylor map in the phase space region of interest. term by evaluating 41(3) e {M31}. +1f(M(1A;11)+1R‘;n +1111... (443) = 0070411) + [Rm + 19].“ (4.44) with the remainder R. Calculations equivalent to those of section (4.5) were performed with this complete map and the results are shown in table (4.8). Guaranteed Turns With Unknown Parameter Pendulum 21,517,254,240 10,064,641,158 Henon 5,831,178 3,123,964 Coup. Pend. 25,021,201 24,506,730 IUCF 19,441,816 10,749,711 PSR 11 49,241,881 39,571,014 Table 4.8: Results of the IC—PIE bounds on the survival time of particle motion for rigorous description of the systems by a Taylor map with remainder intervals. ...A‘Fu—T‘? Chapter 5 Symplectic Scaling (SYSCA) The map which relates coordinates 2',- in the initial plane of an optical system to coordinates 2?} in the final plane contains all information about optical properties. This map depends on the parameters of the system and is expressed by M : R22“? —) B“ for 61 degrees of freedom and p system parameters 6,, such that -o 1 =M(2’.,6). (5.1) “1 It occurred first to Hamilton that optical systems could be analyzed by computing the Taylor expansion of their transfer map [Pra33]. He utilized what he called the characteristic function to compute this expansion to third order for rotationally sym- metric optical systems. This concept was also used in charged particle optics. The motion of the particles is described by a set of ordinary differential equations 8) . (5.2) where s is the coordinate along the central trajectory and 2' describes the usual particle optical coordinates. The general solution of this set of differential equations yields the transfer map of the system. Starting from the 19303, many references can be given which carried the computation of Taylor maps to higher orders. Some important contributions can be found in [8834, BBB64, Wol65, PR71, R0887]. The Taylor 126 127 coefficients are typically described as multiple integrals over powers of solutions of the linearized equation of motion. The kernels of these integrals can be computed with the program MOPS [Pre92] to an arbitrary order; above order seven the formulas tend to become very lengthy. If the fields involved do not change with s and consequently the equation of motion (5.2) is autonomous, then the linearized equation of motion can be solved in basic functions and also the integrations can be performed analytically. In this case the Taylor coefficients can be computed with the program COSY 5.0 up to fifth order [BHW87]. / 3- 30 8+ m 2. % fringe field >9 w Figure 5.1: Definition of fringe—field and main—field maps. Because of the simplifications involved in the regions where the field does not depend on s, traditionally the field of a particle optical device is separated into the main—field region and the fringe—field region as shown in figure (5.1). This figure also describes the definition of the fringe—field and the main—field map. To compute the main—field map Mmf, one assumes that the field ELI, y, s) between the effective field boundaries does not change along the central trajectory and has the structure of the 128 field at the middle 5 = 3m of the element. Outside the effective field boundary the field is assumed to vanish. _, { 63(x,y,sm) for 3 1n the main field (53) Bm a. a = ' ' f(x y 8) for s outsrde the main field The abrupt change of the field violates the Laplace equation and can therefore not represent a physical system. To describe a realistic system, the main-field map has to be composed with fringe-field maps, which describe the connection of the main-field to the field—free region outside the element. The fringe—field map M” can be decomposed into three maps. To illustrate this, let so denote the effective field boundary at the entrance of the element, 3_ be a position so far before the optical device that the field can be neglected, and let 3+ be a position so far inside the element that 3(3, y, 3) changes very little with s. The map Mmf,,o_..,+ describes the particle motion through the main field given in equation (5.3) from the effective field boundary to 3+. The fringe—field map is constructed in such a way that a drift 13,32“, from s- to the effective field boundary composed with M” and then with Mmf,,o_.,+ yields the transfer map M,__.,+ from s- to 3+ by 1W3__.,+ = mf,,0_.,+ 0 Mff O D,__.,,o . (5.4) This leaves the fringe—field map as M” = I" o M,__,,+ o 1.5—1 . (5.5) mf,so->s+ s_—vso Computing fringe-field maps requires the computation of the map M,__.,+ of a sys- tem where 3(3, y, 3) depends on s, and therefore the right hand side in the equation of motion (5.2) depends on s. For this, methods to find general solutions of non— autonomous differential equations are needed. The method of numerical integration in DA will be described, which yields Taylor maps to non—autonomous ordinary dif— ferential equations for an arbitrary expansion order. 129 5.1 Integration in DA Since the theory of differential algebra (DA) and its applications is sufficiently dis- cussed in the literature [Ber92a, Ber91a, Ber90a, Ber89, Ber87, Ra181], only the un- derlying idea will be lined out here. As mentioned previously, one can introduce the notion that functions agree up to order n as an equivalence relation. Equivalence is denoted by “=,,”, which equates functions which have the same partial derivatives at the origin up to order n. In other words, “=n” equates maps which have the same Taylor map up to order n. A differential algebra of the corresponding equivalence classes can be established. Since a Taylor expansion of any map in an equivalence class represents the whole class, this concept leads to a rigorous description of the manipulation of Taylor maps. Functions which can be represented by a computer consist of finitely many elementary operations and elementary function evaluations. The set of these functions is too big to be handled on a computer. The success which the DA concept experienced in the last 8 years is due to the fact that there is a homeomorphism from this set of functions into the differential algebra. If f is a function of this set, then we symbolize the corresponding element in the differential algebra as [f],,. Let A; be an algorithm evaluating a function f : Rd -) [R then Afr? = f(i) for all :1" 6 ”2“. When Eis the identity function in Ed, then Afiz f. Due to the homeomorphism, we can evaluate the algorithm in the differential algebra via «412': f ==> 41121. =.. [f]. . (5.6) “ 7’ The equivalence class to 2,, of the Taylor expansion of f is denoted by [f]... To obtain all partial derivatives of f up to order n, it is therefore sufficient to know the identity and to perform the elementary operations and the elementary function evaluations which contribute to A; in the differential algebra. In the next section 130 we will deal with an algorithm which also contains derivatives. Under conditions explained in that section, also derivatives can be performed in the differential algebra. Generally the equation of motion is solved by numerical integration. Starting with an initial coordinate :2",- and a choice of parameters 6,- with z' E {1,...,np}, the integration produces a final coordinate 2}, which means that the combination of elementary numerical steps performed by the. integrator define an algorithm AM which evaluates the transfer map M(Z,6). According to equation (5.6), evaluating the numerical integrator in the differential algebra yields the Taylor map to any order n. In our applications, the Taylor expansion is computed with respect to 5 as well as to 6. With the Taylor maps, all partial derivatives up to order n are computed. Because of this fact, the described function evaluation in DA is often called automatic differentiation of algorithms. In connection with this method, it is often referred to [Ral81], the basic idea, however, is much older. 5.2 Evaluating Propagators in DA Let L f denote the operator 5Tf+ 0, f. Then the equation of motion dE/ds = f(é', 6, 3) allows one to compute all derivatives of 2' with respect to s by [{2-. -o-o -o -o _ds: = fraf+af=Lf , (5.7) d"" .. _, 23‘: = (Lf~)"“f=(Lf~)"z. (5.8) Let the central trajectory 2’ = 0 be a solution of the equation of motion for the system with parameters 6 = 0, then f(0, 0, s) = 0. The Taylor expansion of f has no constant part and is known to order 77., therefore fr 5 f can be computed up to order n, in spite of the fact that 57 can only be computed up to order n — 1. However, 8,f can only be known to order n — 1. Therefore only maps to autonomous equations of motion, which have 3,}? = 0, can be evaluated in DA by the propagation operator, for which 131 we write 1L 2'1 N1 00 00 M =n exp(L ~) )5: 2L =;L( (5.9) In writing this equation, we assume that all the involved sums converge. To use the propagator for non—autonomous differential equations, an extension of DA would have to be introduced, in which the coordinate s can be evaluated to orders higher than n. In the main—field region the equation of motion is autonomous, and using the propagator in equation (5.9) to compute the map is very efficient. Figure (5.2) shows the time advantage as a function of the evaluation order for dipoles, quadrupoles, hexapoles, and octupoles. Also the following table makes the time consumption of integration in DA apparent and shows the computation speed which is accomplished by the symplectic scaling approximation (SYSCA). It displays the time needed for fitting the multipole fields of the S800 beamline and the S800 Spectrograph, which are under construction at the NSCL, to satisfy 14 conditions on first order Taylor coefficients and 6 conditions on second order. Only main fields with propagator 51 sec. Fringe fields with DA integration 7 hours 10 min. 50 sec. Fringe fields with SYSCA 6 min. 38 sec. In order to analyze if the transfer map can be represented well by the main-field map, monomials are symbolized by if = ?_1 21““; 6f““ with natural numbers kj, j E {1, . . . ,2d + np}. The order of such a monomial is given by III-r.” = 23:?” hi. Now the 7'“ component of the transfer map M can be written as 4 71 ll ll=m _, _, M,- = Z Z (M,|2*) 2* (5.10) i; m=l A I: “"1; ‘ 132 2130 . _ 1840f _fi 1700. A ___— 14701 1280 . —— 1100 . 853. —I 735. 427. 368. 1 . 1 . 1 a . . 1 . . 1 . . r 0 1 2 3 4 5 6 7 0 2 3 4 5 6 1790. - _, ' 605. ._ 1440 . _, 4841 +‘__ 1080. 363. _. 718. 243. _l w_, 360. —_l.—L , 122. 1....l.f.l.l 1...... 3 4 5 6 7 Figure 5.2: Time consumption of integration in DA relative to main—field map com- putation with the propagator for different expansion orders. From top left to bottom right: dipole, quadrupole, hexapole, octupole. 133 with symbols (rylgli'f) for the Taylor coefficients. Figure (5.3) shows the normalized average deviation Am between the Taylor coefficients of the transfer map and those of the main—field map in order m for four typical charged particle optics devices. The wedge dipole in the examples of this section has radius 2m, an angle of 30°, and an aperture of 2.54cm. All other multipoles have: length 41.9cm, a pole—tip field of 2T, and an aperture of 2.54cm. The chosen ion is 1.603" with an energy of 25MeV per nucleon. The fringe fields used are those of the Enge model [Ber92b, KE87, BSBI]. Am is computed by 2 ”1%:le I(Milzdf)‘ (Mmfilzffl 111)i=1 1? I(Mgli")l-l~l( anli'kfl (5.11) If the coefficients (M,|z"f) and (Mmf,,-|5f) are randomly chosen, then Am is on average one. This can be seen by the following integration. If the coefficients are allowed to vary in an arbitrary interval between —N and N, one can write N ... l l _ 16(2) = [:1] '“ b'dadb=1+2/O/O '“ b'dadb (5.12) N|a|+|b|4N2 a+b 4 — 1+ / /:+ba;2bdadb— / [Ma 2b)dadb (513) _ 2 . = %(1+ / {1-2b— 2616(14 +b)}db) (5.14) 1 , 1+b) = _ _ __ .1 2(1 [b1n(—— )16— / 1+bulb) (5 5) 1 1 = §{1+ln(2)—/0(1-1—+—b)db}= 1n2(). (5-16) Obviously, some trivial computational steps have been skipped. Figure (5.3) shows that the so called “sharp cut off fringe field” (SCOF F ) approxi- mation, which replaces the transfer map by the main-field map, is very inaccurate for nonlinear coefficients. Some averages are even larger than one, and thus on average 134 1.090 0.873 0.655 0.436 0.218 0.970 0.776 0.582 0.388 0.194 0.000 . . 4 4 4 0.000 1 . . 2 3 5 6 7 3 5 6 1.130 1.... 0.926 —— r— 0.904 _ 0.741 0.678 0.556 _— 0.452 0.370 0.226 0.185 0.000 A . 4 4 . 0,000 1 J . 2 3 5 6 7 3 5 6 Figure 5.3: Accuracy of SCOFF for different orders, a random choice of Taylor co- efficients would be close to one. From top left to bottom right: dipole, quadrupole, hexapole, octupole. 135 the approximation is worse than a random choice of coefficients. More examples of the insufficient accuracy of the SCOFF approximation will be demonstrated in section (5.5) where several examples are discussed. These results show that the advantages of very quick DA evaluation of propagators can only be taken advantage of if an efficient fringe—field approximation can be found. 5.3 Desirable Properties of an Arbitrary Order Fringe—Field Approximation After stating three properties which we consider desirable for a fringefield approx- imation, the importance of each property will be explained in separate subsections. The approximation should: 1. lead to order n symplectic maps, 2. represent the fringe effect well for a wide range of apertures, 3. be usable for arbitrary orders. The simplest approximation, as already described, is SCOFF, where fringe fields are simply ignored. As illustrated in section (5.2), this method strongly violates point 2. The impulse approximation [Hel63] used in the code TRANSPORT [BRCI77] violates the points 2 and 3 and the method of fringe—field integrals [HBW90, HIW93, Wol65] used in the computer code GIOS [WHB88, W0192] violates the points 1 and 3. 136 5.3.1 Order 12 Symplecticity The necessity for order n symplectic maps is especially apparent when long term behaviour in storage rings is to be analyzed. If the Taylor map does not approximate a symplectic map very well, then, according to theorem (1.6.4), the phase space area is not conserved in every 22,-_1 x 22,- phase space section. This is the case for order n symplectic maps, when the evaluation order n is small. In [Yan91], for instance, it is illustrated that in case of the SSC lattice, expansion to 12” order is needed to have a sufficient degree of symplecticity. If the fringe—field map is not computed to be order n symplectic, even very high evaluation orders can not lead to trustworthy predictions of long term behaviour. In figure (5.4) the tracking picture of a 15° dipole for the horizontal plane is shown. On the left side the coordinates used are the conventional coordinates x—px, whereas on the right side the phase space points are displayed in normal form coordinates. Normal form coordinates are used since violation of the symplectic symmetry can most easily be seen in these coordinates, where the motion lies on circles if the map is symplectic. The top two figures illustrate the action of the transfer map of the dipole. Since 24 applications of the map yield the identity, the motion is stable for all times. Approx- imating this map by a fifth order Taylor map does not approximate the symplectic transfer map sufficiently and the tracking points move further and further away from the origin, as if the motion were unstable. This is shown in the third and fourth figure in conventional and in normal form coordinates respectively. To make the map completely symplectic, one can compute the Taylor expansion of a generating function from the Taylor map. The generating function can then be used for tracking particle coordinates through a completely symplectic approximation of 137 . . . .. (n) ................ O ...... ...._ /’:-_'.'\ //""- “k //.". T. t ?\ ’I o ' - ”\Q / o . . \ I, ' o . . v: r \ ”' . o . .s \\ o . o . . ‘\ I], O \ \ I": . .‘ .‘\\ “... ...\\ \\..:,.:.:,.()..:3.:. ...... ..,.I\\ \‘,..........O.:.:......1\ \‘°. \" "' \\\ . .. 0" ’/ \\‘.. 0” \\~‘~'_- -' .5 >//’/ \45. . . v.,'// Figure 5.4: Tracking through a 15° dipole, displayed in conventional coordinates on the left and in normal form coordinates on the right. The figures on the top display accurate tracking data, those in the middle were created by the 5th order Taylor map, and the bottom pictures were produced by 5“ order symplectic generating function tracking. 138 the transfer map. Strategies for symplectic tracking are discussed in [Ber88a, Ber91b, Yan93, Gja93]. This method relies on the order n symplecticity of the underlying map and therefore only a fringe—field approximation which leads to order n symplectic maps can be used for symplectic tracking. The two graphs on the bottom of figures (5.4) show tracking with the fifth order approximations of a generating function. Again the left involves conventional generating function tracking, whereas the right displays symplectic normal form tracking. It is apparent that the transfer map which creates the top two figures is only approximated by the generating function which creates the bottom two figures, but the permanent stability of motion is correctly represented due to the fact that the generating function represents a symplectic map. Still another important reason for striving to find order n symplectic fringe—field maps is the interrelations between different Taylor coefficients which are imposed by symplecticity. Some systems rely on these interrelations to correct aberrations. Examples are the high order achromates described in the references [WGBQ3, WGBQ4] and the opening aberration correction for electron microscopes by means of hexapoles described in [Ros90, Hof91], to mention only a small sample of the wealth of systems which were designed with the symplectic condition in mind. All these systems can not be analyzed well if fringe—field maps are used which are not order n symplectic. 5.3.2 Accuracy for a Wide Range of Apertures Designing spectrographs, simulating electron microscopes, and analyzing beamlines and high energy storage rings areall topics whiCh rely on the same principle of calcu- lating transfer maps from the equation of motion. It therefore seems highly desirable to formulate calculation methods which work equally well for all these subfields of particle optics. The effect of the fringe fields in these areas is of different importance. In high energy accelerators, fringe fields typically have less influence, however, also 139 in this field substantial influences can be observed. Most fringe-field approximations assume that the region of the fringe field is much shorter than the region of the main field. For large aperture spectrographs and for electron microscopes, this is typi- cally not the case. Therefore, to compute fringe—field maps of general usefulness, this assumption has to be avoided. 5.3.3 Usability for Arbitrary Orders With the described DA methods it became easily possible to compute Taylor maps to very high order. The limit is only given by computer memory and computation time. To allow the potentials of the DA approach to be used with fringe fields, the com- putation of fringe—field maps should work for arbitrary evaluation orders. Especially since an efficient way of computing fringe—field maps was sought for the arbitrary order code COSY INFINITY, this requirement was important for the development of symplectic scaling (SYSCA). 5.4 The Principles and Usefulness of SYSCA At first the basic ideas of symplectic scaling will be outlined and it will be explained why one can expect that our approach will yield accurate fringe—field maps. This method was first outlined in [HBQL HBQ3C]. The accuracy and speed of the approx- imation will be illustrated on several examples including linear design, high order effects, and long term tracking. These examples are also analyzed in [HBQ2a, H394]. All the steps needed for symplectic scaling will be described in detail; since they can not only be used for fringe—field maps, the scaling principles will be formulated for general maps. Following, we will use geometric coordinates for the map which means the variables 140 of motion are the cartesian coordinates and slopes x, :1", y, y', the path length I, and the momentum p. Those six quantities form the vector 59 = (:r, :r’ ,y, y’ , l, p), which is transformed by a transfer map which depends on 71,, parameters Dj, *1 6,; = 417(6613) . (5.17) N; t Here and in the following sections, parentheses sometimes symbolize functional depen- dences and sometimes they symbolize function evaluations. For example, [PI (59, D) is a map ”22“"? —+ R“, whereas 0269.213) = NILE}, D)l(é‘g,.-J3) in equation (5.17) is an element of R“. Whether a map or a real vector is described will become clear from the context. The parameters we want to consider here are the mass m and the charge q of the particle, the size A of the optical element, and the pole—tip field B. When the size A is changed, the length and the aperture of an element change simultaneously. Call 112(59, m, q, A, B) : R“ x R4 —> If“ the general parameter dependent transfer map of an optical element with the vector of parameters D = (m, q, A, B). From now on we will restrict ourselves to magnetostatic elements, although parts of the procedure are applicable to electrostatic elements, too. Well known scaling properties allow to compute the general parameter dependent map M(Eg,m,q, A, B) from the specific transfer map Mmi’qi'Ai’B'Qfl : R“ ——» B“ for one particle type with mass m." and charge q“ and one parameter setting for the optical element of size A" and field B‘. This is due to the fact that maps in geometric coordinates observe two scaling properties. The bending radius of the path of a particle with momentum p at field B is P = — 5.18) (131 ( where B .L denotes the field component perpendicular to the momentum of the particle. All maps that describe particles with equivalent bending radii along their path are 141 identical. If the specific map Afmi’q'AhBXEg) is known for geometric coordinates 59 = (0:, cr’,y,y’, l,p), then it is known for all kinds of particle momenta p. If the specific map for any other parameter vector (m, q, A‘, B) is desired, one only has to choose p appropriately to create a map which describes the trajectory that corresponds to this parameter setting by ‘ I qB M77436.) = me‘s‘fi‘w, y, y', I, p- 80 ) . (5.19) This is just another way of saying that the map depends only on the ratio of field to magnetic rigidity. Therefore the first scaling property is called scaling with magnetic rigidity. Strictly speaking this property only holds if the field structure B(:r, y, s) / B does not change with the pole-tip field, which means that saturation effects are not important. However, even when saturation effects are important, SYSCA is very accurate over a substantial range of the field strength. One can also compute the field dependent map Mmi’qi’AXEg, B) : I?“ X R -> B“ from the specific map Mmi'qi'Ai'BXEg) by scaling with magnetic rigidity. This map will become important for geometric scaling. From the field dependent map Mm."7"‘4'(é'g, B), we can compute the specific map for any different size A of the element. For that, let us consider two similar magneto- static elements that differ only by a scaling factor a. If the bending radii also differ by a factor of a, the maps are similar. Equation (5.18) shows that this is the case whenever the increase in size by a factor a is accompanied by a decrease in the field strength by the same factor. After scaling the coordinates x, y, I, we obtain 7 41:14:96.1) 644201434) 1 E .q 1413(5'9) M: 'q ’A (2.913%) Mm.’q.’A 8(59) : 7447Mm. (I. ”(591375) (5.20) M7"“""B(2g) ”"F’" 'A (56.8%) \ Mzm.’q.’A’B(gy) A (rm-1.1.10) \ ffflfnt’q. A (2},Bfi) ) ($"§"i"§‘p) ' 142 / / I - l/ '3 ——- ~19 - Aal/ l/ Jl B/a A 5 Figure 5.5: The coordinates of the particle trajectories in two elements scale with the factor a if the elements scale with the factor a and the fields scale with the factor 1/ a. Geometric scaling (5.20) and scaling with magnetic rigidity (5.19) together allow one to compute any specific map Mm'q'A'B(2'g). Since this holds for any parameter vector D = (m,q, A, B), the general map A269, m, q, A, B) can be computed from a «‘4'48751. specific map M m. To conclude, we state that the knowledge of the transfer map Mmi'qi'AXZg, B) as a function of the field strength at the pole tip for particles with a specific magnetic rigidity and for an element with a specific aperture is sufficient to know all transfer maps of similar elements for all energies, masses, and charges. In fact the map does not even have to be known as a function of the field because the dependence on the field can be obtained by equation (5.19) from the dependence of the map on the momentum. -o In accelerator physics, the phase space coordinates are often chosen to be 2 = (22,:6’, y,y’,6;, 6p). 61 is the difference between the length of the particle’s trajectory 143 and the reference particle’s trajectory. 6p is the relative deviation from the momen- tum p of the reference particle. It should be noted that in the literature the reference momentum is often denoted by p0; this would lead to confusing notations in the for- mulas which are to follow and therefore p is used to describe the reference momentum. The particles momentum, which was called p in the geometric notation, is given by p(1+ 6,0). The coordinates 2}, are introduced to ensure that particles which move close to the reference particle have small phase space coordinates. Then a Taylor expansion of the map is usable to describe the weakly nonlinear motion. Often this notation is called TRANSPORT notation [BRCI77], it can be transformed to the geometric notation and the scaling properties can be applied. In general it is useful and customary to work with the canonical coordinates [Ber90a] 55 = (x, a, y, b, 7', 613), which are described in detail in subsection (5.8); 65 is the relative deviation from the reference energy E. These coordinates are also small when the particle is close to the reference particle and Taylor expansions can be used. Maps in canonical notation have to be transformed to TRANSPORT notation before the two scaling laws can be applied. If D is a parameter and the Taylor expansion around the reference parameter D" is needed, we use 60— = (D - D‘)/D"' as parameter. Given a reference pole—tip field B’, then 63. can be used to expand the field dependence of a map around B'. COSY INFINITY readily computes the Taylor expansion of a field dependent map Mmf’qi'A'GE” 63.) at reference values indicated by “at”. This Taylor map has to be stored to a file in order to compute maps of similar fields and all kind of beams by the previously mentioned scaling method. The accuracy with which the Taylor expansion approximates the field dependent map depends on 144 the chosen order of the expansion and on the difference between the reference field B“ and the field B, to which one has to scale. Unfortunately, the symplectic structure of the map would not be conserved in this process. While Mm"q.'A° (EE- , 6 3°) is symplectic for all 63. and the corresponding Tay- lor map is order n symplectic, this is not the case for the Taylor map Mmi’qi'AXE'E. , 63.) 5B, which is B“ -> R“ and is obtained when 63 E R is inserted in the field depen- dence. Order 71 symplecticity, however, is an intrinsic symmetry of canonical motion that arises from the special structure of Hamilton’s equations and should not be violated as mentioned in subsection (5.3.1). This drawback can be eliminated by transforming the field dependent Taylor map Mm',q',A'(EE., 63.) to a symplectic representation, either in the form of a generating function, or in the form of a Lie exponent Mm.’q.’A.(Z-..EO, 63‘) = L(68.)6:P(EEo,6Bo)zf (521) with the identity map L In higher orders the representation via generating functions is slow, because a map inversion is required. The Lie representation has the disadvantage that the matrix L(6B.) 58 is not exactly symplectic. A combination of both methods is most efficient. We represent the nonlinear part by P(EE.,63.) and the linear part by the generating function F (6113-) that is most accurate for the given matrix L(6B~). The left of figure (5.6) describes the procedure of computing the symplectic ref- erence representation, which is stored in a reference file. First the Taylor map Mmi'qi'Ai’B'GEo) for a certain size of the optical element and for a certain parti- cle is computed in canonical variables. In our implementation this is done with the code COSY INFINITY. Scaling with magnetic rigidity yields the field dependent map Mm"°"A‘(EEo,6B-). From this map the symplectic representations F (63-) and P(EE-,63.) are computed and stored in a file [see subsections (5.6) and (5.7)]. 145 The Taylor expansion of the general parameter dependent map NIGE, m, q, A, B) around (m,q,A, B) E R4 is written as 117(55,6m,6q,6,4,63). The goal is to compute this Taylor map from the reference representation. The right side of figure (5.6) il- lustrates the procedure which will achieve this goal. After reading the reference file, the first step is to insert a value 63. into the field dependence of the symplectic rep- resentation. 63. is chosen in such a way that the symplectic representation describes particle motion which differs only by a scale in geometrical size from the trajecto- ries of the desired map Mm’q'A’B(E'E). If A = A‘, then equation (5.19) yields that B, = 8&3. Then scaling with magnetic rigidity will create the map of a particle with p and q for the field B. If the desired element is bigger by a factor f7, then the field B, has to be chosen bigger by the same factor, due to equation (5.20). Then geometric scaling to size A will bring the field back to 351%. Therefore 63. = (B p M — B’)/B" (5.22) m A has to be inserted. From the so scaled symplectic representation one computes the order n symplectic Taylor map [see subsection (5.6) and (5.7)] which is then transformed to TRAN S- PORT coordinates [see subsection (5.8)]. The transformation map T(E,m) from canonical to transport coordinates depends on the properties of the particle which’s motion the map describes: ... ... Map, )= T(E,m) o MGE, )o T—1(E,m). (5.23) The map in TRANSPORT coordinates 5,, can then be scaled to the desired field by scaling with magnetic rigidity and to the correct size of the element by geometric scaling. As mentioned before, we can obtain the general parameter dependent map from the specific map by scaling. With DA this procedure automatically leads to the Taylor expansion with respect to all the parameters. 146 Finally one obtains the required general parameter dependent order n symplectic Taylor map 113(53, 6m,6q, 6 .4, 63) by transforming back to canonical coordinates [see subsection (5.8)]. Symplectic scaling can be applied to any map, but it is especially useful for fringe- field maps. The following examples show the accuracy and speed of obtaining fringe- field maps by symplectic scaling of a stored symplectic reference representation. 147 -o M(ZE76m36Q16A763) H (5.23) at E, m 19.1121), 6m26Q9 6A, 6B) 1 Mmc’qo,Ao’Bo _. . -o , 7 (ZE ) Mm'q’A’B(2p) Mm‘vq‘vs‘(26.,66.) 17464446,) I] (5.21) M (5.19) Fm.’q.'A.(6B.)’ Pm.’q"A.(EE.,6B¢) MmO'qt’A‘,Bs(EP.) fl H (5.23) at E‘, m" [reference file] M"WWW"?(367') H (5.21) Fmi'qi'Ai’B', Pmi".'A°’B’(EE.) E=Bfl4 q‘p A. beference file] Figure 5.6: Schematic outline of the procedures leading to the reference representation and from the reference representation to the general transfer map. 148 5.5 Examples In this section, we will illustrate the profitable use of the method with several ex- amples. In order to evaluate speed and accuracy of the proposed approximation, we study a certain aberration coefficient of a quadrupole. Figure (5.7) shows the depen- dence of the expansion coefficient (xlzrza) as a function of the field B at the pole tip. Because functions like this can be closely approximated by polynomials, symplectic scaling (SYSCA) is very accurate. 0 10‘4/m -8 P10810(A) 3 —8 . —-10 . 4 —16 . -12 . /5 -24 . . . . _14 \. . . . 0 1 2 3 4 B/T 0 1 2 3 4 B/T Figure 5.7: left: (xlxxa) for a quadrupole as a function of the field at the pole tip. right: Error A of the approximation of (:clxxa) with different expansion orders for the reference representation at B =2T. Even at the border of the range in figure (5.7) the presented method is more accurate than the COSY standard integrator. Close to the value with which the reference file was produced, the accuracy increases drastically. The results in figure (5.7) were obtained by evaluating the symplectic reference representation to third, fourth, and fifth order. The accuracy can be further improved by increasing this 149 order which of course increases the computation time that has to be invested for creating the reference map in advance. This investment can be very much rewarding, especially when beamlines or spectrometers are being fitted or when system errors are analyzed so that maps of similar fringe fields are needed over and over again with only slightly different parameters. The SYSCA approximation is especially helpful in the design of a realistic system after approximate parameters of the elements have been obtained by neglecting fringe fields. These values can be used to create a reference file for symplectic scaling. In this way, a very high accuracy almost equivalent to accurate but time intensive numerical integration can be obtained. The time advantage of this method is illustrated in figure (5.8). Fringe fields do have noticeable effects already in first order. In the example of the A1200 [She92] isotope separator at the NSCL, the effect of the fringe fields on the calculated setting of the field strength is shown in figure (5.9). The fringe fields were described by Enge functions, and the Enge coefficients had been fitted to measured field data. Here the time advantage of the proposed approximation in the fit is three minutes versus two hours. As a measure of accuracy, we study the tilt angle 0 of the dispersive image plane and the opening aberration Co for various approximation methods. In the discussed device the coefficient (:rlaa) vanishes because of symmetry of the axial ray and anti symmetry of the dipole fields; therefore (mlaaa) is the relevant opening aberration, _ _M _ CL. aaa 9 — (ala)(.r]6) , Co —( I ). (5.24) Table (5.1) shows 9 and Co for various fringe—field models. The values of O with and without fringe fields differ by 0.5% for the first dispersive image plane in the A1200; the third order aberration, however, is completely wrong if fringe fields are 150 177.0 . ,F. 57.6. __ 142.0 . 46.2 . _, _, 107.0 . 34.9 . _ r—1 _— 71.6 1- —‘ 23.6 I- 36.3 . 12.3 . 1.0 l l l l l L L 7 1.0 J l L L l l L 0 1 2 3 4 5 6 7 0 1 2 3 4 5 62.4F _. 72.6 .. _ 50.1 . 58.3 . _ 37.9 . 44.0 . f— 25.6 . L___ 29.6 . _r- 13.3 .. 15.3 . l 1.0 4 l l l l l A 1.0 L l L l l J Li 0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7 Figure 5.8: Factor of time advantage of SYSCA to numerical integration with accu- racy of 10’”9 as a function of the expansion order. From top left to bottom right: dipole, quadrupole, hexapole, octupole. 151 r- Yj I 1 I T7 I I I I I I r I I TI -‘ 1.000 :- a an 5 aa a J ._ fl ‘ - 0.975 :— 0 ——7 Z x Z I 1 0.950 :— — 0.925 :— + COSY With Fringe Field J i 0 c1 COSY with SYSCA : o 9 o E x SCOFF 1 ' ° t x 0 TRANSPORT with Fringe Field _— h I L I l J 1 L l I J l l l l l l l 1 .fi 1 2 3 4 5 Figure 5.9: Relative deviation of predicted field settings with SCOFF and SYSCA from the correct settings for five quadrupoles. The standard fringe—field approxi- mation of TRANSPORT is given as a reference; the deviation is mainly due to the neglect of quadrupole fringe fields. O and Co with SCOF F approximation 80.8840° -65.96m O and Co with dipole fringe fields only 81.1696° —65.96m O and Co with quad fringe fields only 81.2694° -—682.68m O and Co with SYSCA approximation 81.2701° —687.10m O and Co with actual fringe fields ' 81.2702° —687.10m Table 5.1: Tilt angle and opening aberration for various fringe—field models. disregarded. This comparison also shows that quadrupole fringe fields, although often disregarded, can have effects which dominate over dipole fringe fields. Nonlinear effects can be seen by sending a cone of particles through the 7th order A1200 map. The images with SCOFF and SYSCA approximation are shown in figure (5.10). The maximum angle used is 15mrad. Note that due to the difference in (scale, the beam spot computed with SCOFF is only one tenth as big. Trusting SCOFF would lead to a loss of most of the beam. 152 m 1mm 0.1mm 1——1 - 1—-1 Figure 5.10: Beam spots with SYSCA (left) and SCOFF (right) approximation. The plot produced with the exactfringe fields can not be distinguished from the plot produced with SYSCA. Note the difference in scale. Also for long term tracking in storage rings, fringe fields are influential. In figure 5.11 typical tracking pictures are displayed. An example storage ring for 1GeV pro- tons was optimized for a big dynamical aperture. Eight particles were tracked with phase space coordinates which had :6 and y components in order not to avoid x—y coupling. The left picture was computed without fringe—field maps, whereas the right picture was computed with fringe-field maps. When fringe fields are neglected, one can often not trust the computed dynamical aperture. O ’ ... . 0". .' ... ...” C... . I - Q .. n ......'..o.. '- 0‘ . ' :0, ~ “~- ~ . a O O .... O - ~ ~... . .I .' ' 1’4”" "‘1 if f f. '- 3 ’ ’. °°°°° ~ ‘ ”(’{b ‘\\‘ '. {I ’0. ...“ —. t . ‘ ..- ,, I row — ' : I .0 ‘ i , — E .l/l ’ f ‘9 \ ‘ 0 .' O ‘ . ‘.‘ o ‘. o . . .. . .0 ’ ’ ‘ 0 ~ ". Q. ‘. t..oo ‘." t ‘ ‘ ~ . ... I 0' ' ’1 . . ..o g .0... ...-4; ~ : LT... 1:0. 2; . ‘. ~ 0 Figure 5.11: 500 turn tracking without (left) and with (right) fringe fields trough an 1Gev proton storage ring. The effort involved in generating a symplectic approximation is rewarded when 153 repetitive tracking is being performed. The example lattice of choice is the proposed PSR II Ring. The 9“ order 5000 turn tracking pictures are displayed in figure (5.12). ... . 0‘ ‘. -‘. ... .' 0. Q. a .V’f .... '. '0. 00* . . ‘0... . . o ‘0 o . ... ..Q' ' ‘ "‘...:'~O#: o I20mrad 2cm I—-—i 2cm I20mrad 1-—-1 Figure 5.12: 5000 turn tracking with fringe fields obtained by numerical integration (left), SYSCA (middle), and a non—symplectic fringe—field approximation (right). The initial position of the particle is (x, y) = (3cm, 3cm) with no initial inclinations a." and I y . The tracking was performed with the described standard numerical integration, SYSCA, and a non—symplectic fringe—field approximation obtained by low accuracy numerical integration. Non—symplectic tracking strongly violates the conservation of phase—space volumes. SYSCA yields more stable results than the numerical integra- tionlsince the limited accuracy of the numerical integrator slightly violates symplec- ticity. The corresponding 9“ order maps were produced with the SYSCA mode in COSY INFINITY in 30 minutes, whereas the standard numerical integration took 15 hours, and the non-symplectic approximation took 44 minutes on a VAX 4000 90 computer. It is worthwhile mentioning that only maps of entrance fringe fields have to be saved. The exit fringe field of an element is described by the reversed map. A reversed map is computed by inverting a map and changing the signs of the incoming and the 154 outgoing momenta: ( Mr“) _Ma-l My"l _II’Ib—l _MT-1 1 115,3) Mrev (Ir—a'yv‘by—TJSE) Computing the inverse map 03" from the SYSCA reference file is just as direct as computing the transfer map [If itself. The generating function is a mixed function of initial and final variables. To compute the transfer map, it has to be partial inverted with respect to the initial variables. Partial inversion with respect to the final variables is a very similar process and leads to the inverse map. For the Lie exponent, exp(— : P :) acting on the identity f leads to the inverse of exp(: P :)L If a representation approximates maps with pole—tip fields close to B then fields close to —B can be approximated, also. The reason for this advantage is that the change of the sign in a multipole of symmetry Cu corresponds to a rotation of that element by 1—39- degrees. 5.6 The Lie Exponent We search for a representation of the order n symplectic Taylor map M : B“ —> B“ of a Hamiltonian system in the form M :n Lo (e‘fln+l)‘5) (5.26) where again “=,,” indicates that the right and left side agree up to order 72. Also here we use the previously advocated notation of the identity map 2', not to be confused with a real vector. L is the linear map and the Lie exponent PM“) is a polynomial of orders from three to n + 1. With the following induction, we will proof that the 155 Lie exponent always exists and the induction proof itself will reveal a methode to compute this exponent. Since the Jacobi matrix L of L is symplectic, it has an inverse. As in chapter (1), we write N for the Jacobian of a map 1V . The exponential has to satisfy e:P(n+l);§:n iii—1 0 AI: 52+ N (5.27) with the nonlinear map A7. The right—hand side has a symplectic Jacobian I + N, where I symbolizes the identity matrix. Symplecticity requires (I +N1TJ21(I +N) =..-1 12., (5.28) Jng + NTJM + NTJng z.-. 0 . (5.29) The first order of this equation shows that Jngz is symmetric. The parts of exact order m of 117 and PM“) are denoted by 117", and Pm respectively. This implies that the potential problem given by the second order of equation (5.27) 2.69315: —J2d5P3 = N2 (5.30) has a solution P3. The next order yields the equation .1 1 :P4 : 5': N3 — §(: P3 :)25. (5.31) If this potential problem for P4 has a solution the next step to obtain an equation for P5 is obvious. To show that all the appearing potential problems have a solution, we start with the assumption that P), was previously found as a solution of the potential problem IPk 3 E=k_1 5+ N — €:P(k_l):5 . (5.32) The potential problem for P1,“ is ZPk+1 25:]; 2+ N — 8:P(k):2-5. . (5.33) 156 We call the right hand side F), and note that it has no orders lower than k. This problem has a solution Pk“ if .1ng), is symmetric to order k — 1. Let the symplectic matrix I + Em be the Jacobian of eiflk):2', then the symplectic condition requires similar to equation (5.29) that szE(k) + E‘kfl‘J... + E("’T.J2dE(’°) = 0 . (5.34) To ensure the symmetry of Jngk, one needs to show that J2d(N — Em) + (NT - 13””)sz =k-10 - (5-35) The equations (5.29) and (5.34) show that this condition is equivalent to NTJng =,,_l ETJ2dEU=> . (5.36) Since we have kept in mind that N is nonlinear and therefore the Jacobian N has no constant part and furthermore that N =k_2 EU‘), we realize that the right and left side of equation (5.36) agree up to order k — 1. 5.7 The Generating Function for the Linear Map For the task of representing fringe—field maps which are to be stored as Taylor expan- sions of the magnetic field, it is important that all the manipulations can be done with parameter dependent maps and parameter dependent Lie exponents. The linear map E then corresponds to a matrix whose elements are Taylor expansions with respect of the parameters. In the specific case of SYSCA, the parameter is 63, the relative deviation of the pole—tip field from a reference value. The corresponding matrix L is symplectic to the expansion order of 63. However, if a specific value for 63 is inserted the resulting matrix is in general not symplectic. Therefore a symplectic represen- tation has to be found that agrees with E to the expansion order of the parameter 157 63. The generating functions are a suitable approach. We are concerned with the six dimensional space of transverse momentum, energy, and the corresponding canonical variables. We will keep the arguments general, however, because the procedure is the same for any symplectic 2d dimensional space. The following matrices will all be either d or 2d dimensional. The symplectic matrix L will be written in terms of d A B 1;:(0 D). (5.37) The first d dimensions describe coordinates (f and the second of dimensions describe dimensional submatrices: the corresponding momenta 13’. J2.) therefore has the structure J2.) = ( _3 3 ) . (5.38) Because of symplecticity _ DT —BT L I: 42.?sz = ( _CT AT ) (5.39) and therefore -1 __ ADT — BCT BAT — ABT _ LL ‘ ( CDT — DCT DAT - CBT ‘ I ' (5'40) The four commonly used generating functions depend on the initial and final po- sitions (22/; and momenta 15', / f. The relationship between final and initial coordinates is established by the equations ((3.5!) = ( 35.,~35,)F1(15i,15})» (5.41) (til-.153) = ( 3*.» 35,)F2(I31.