NONLINEAR DYNAMICS OF CHARGED PARTICLES: ANALYSIS, COMPUTATION AND SIMULATION By Robert Hipple A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Physics - Doctor of Philosophy 2019 NONLINEAR DYNAMICS OF CHARGED PARTICLES: ANALYSIS, COMPUTATION ABSTRACT AND SIMULATION By Robert Hipple This is a multi-tier thesis, centered on modeling the dynamics of charged particles in storage rings and a time of flight mass spectrometer. Topics include: • High order storage ring tracking with fringe fields (Chapters 2,3) • Low-level code benchmarking (Chapters 4,5) • Second order electrostatic element analysis and modeling (Chapters 6) • Oak Ridge Isomer Separator and Spectrometer (Chapter 7) There is no shortage of simulation software in beam physics. Faced with the myriad of options, a researcher will typically select one or two reliable codes, learn them well, and stand by them. It is the very nature of computer simulation that different codes have various strengths and differences, relating to the level of relevant physics modeled and numerical methods applied. This situation underscores the importance of benchmarking multiple codes against one another. This philosophy holds especially true in accelerator modeling, where nonlinearities and the sheer complexity of the simulated scenarios make analytic verification of results difficult, if not impossible. A methodology to explore whether the results of a simulation represent real physics is to run the same test across multiple codes and compare the results. While pursuing my graduate studies at MSU, I have been able to implement this multi-code strategy in several contexts. Simulating the COSY-J¨ulich storage ring in MAD8 and Zgoubi, the dynamic aperture seemed acceptable. However, upon porting the lattice to COSY Infinity, we discovered the surprising result that the lattice as specified was unstable when fringe fields were taken into account. This result would not have been apparent had we limited our simulations to only the original two codes. A similar analysis was performed on the HESR storage ring to be constructed at FAIR at GSI. There it was discovered that the nonlinearities which limit dynamic aperture were not visible while using impulse approximation modes in the codes MAD8 and MADX. Porting the lattice to COSY INFINITY, we were able to model the nonlinearities to a high degree of accuracy. Often when codes disagree in their predictions, the cause of the discrepancy is not immediately apparent. Such was the case when comparing the codes LISE++ and COSY INFINITY in their modeling of spherical and cylindrical electrostatic deflectors. Disagreement in the second order aberration initiated a full analytic investigation into the nature of these aberrations. New physics clarified by that work is also presented in this dissertation for the first time. Finally, an analysis of the Oak Ridge Isomer/Isotope Spectrometer/Separator (ORISS), was performed. ORISS is an electrostatic multiply reflecting time-of-flight (MRTOF) mass separator that was built by the University Radioactive Ion Beam Consortium (UNIRIB) [41]. MRTOF devices are typically used at very low particle intensities due to strong Coulomb repulsion at the particle turning points. Questions were opened on whether the device could operate effectively in a high resolution mass separation mode at high particle intensities. Simulations performed on the iCER High Performance Computing Cluster at MSU show that the device can be operated effectively in this mode at high intensities if voltages are adjusted properly. This was the first analysis of the device to take the effects of intense space charge into consideration and the results are presented in this dissertation for the first time. Dedicated to my wife, Gina iv ACKNOWLEDGMENTS I would like to acknowledge those without whom this Dissertation could not have been completed! Many thanks, in alphabetical order, to Dr. Daniel Bazin, Dr. Phil Duxbury, Dr. Steven Lund, Dr. Peter Ostroumov, Dr. Mike Syphers, Dr. Oleg Tarasov, and Dr. Kirsten Tollefson for all of their help. Also, a special thanks to Dr. Lund for editing the early versions. v TABLE OF CONTENTS LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix Chapter 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 Orbits . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 2 Fringe Field Analysis at COSY J¨ulich . . . . . . . . . . . . . . . 2.1 The COSY Storage Ring . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Adding Fringe Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Fixing an Unstable Lattice Model . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Side by Side Comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Symplectic vs. Nonsymplectic Tracking . . . . . . . . . . . . . . . . . . . . . 2.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 3 The High Energy Storage Ring HESR (COSY Infinity, MAD8, MADX) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 MAD8, MADX and Thin Elements . . . . . . . . . . . . . . . . . . . . . . . 3.2 Tracking and Fringe Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 4 Benchmarking COSY Infinity, ZGOUBI, MAD8, GICOSY . . 4.1 Modeling Dipoles With Edge Angles . . . . . . . . . . . . . . . . . . . . . . 4.2 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 5 Comparative Tracking . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Model for Benchmarking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Test 1: Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Test 2: COSY INFINITY . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Test 3: Zgoubi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Test 4: MADX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 6 Electrostatic Deflectors for LISE++ . . . . . . . . . . . . . . . . 6.1 Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Comparison of Analytical Results, COSY Infinity . . . . . . . . . . . . . . . 6.2.1 First Order Electrostatic Cylindrical Deflector . . . . . . . . . . . . . 6.2.2 First Order Electrostatic Spherical Deflector . . . . . . . . . . . . . . 6.3 Lagrangian Analysis for Comparison With Runge-Kutta . . . . . . . . . . . 6.4 Term by Term Comparison of Analysis and RK . . . . . . . . . . . . . . . . 1 3 11 16 16 21 24 27 36 45 46 48 53 56 57 57 70 71 71 73 73 74 76 77 78 79 83 83 87 90 94 vi 6.5 Independent Check Following Berz,Makino,Wan . . . . . . . . . . . . . . . . 110 6.6 BMW to Second Order . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 6.7 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Chapter 7 Oak Ridge Isomer Spectrometer and Separator (ORISS) . . . 121 7.1 Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121 7.2 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130 7.3 Overview of Operational Tuning . . . . . . . . . . . . . . . . . . . . . . . . . 132 7.4 Particle-In-Cell (PIC) Simulations . . . . . . . . . . . . . . . . . . . . . . . . 140 7.5 Space Charge Effects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 7.6 Tuning for Mass Separation at High Space-Charge Intensity . . . . . . . . . 156 7.7 Simulation of Realistic Distributions . . . . . . . . . . . . . . . . . . . . . . 172 7.8 . . . . . . . . . . . . . . . . . . . . . . . . . 178 7.9 Summary of ORISS investigations . . . . . . . . . . . . . . . . . . . . . . . . 182 ISOLDE MR-TOF Simulations Chapter 8 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 Chapter 9 Talks and Publications . . . . . . . . . . . . . . . . . . . . . . . . . 187 APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 189 APPENDIX A Fringe Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 190 APPENDIX B Second Order Beam Dynamics . . . . . . . . . . . . . . . . . . . . 196 APPENDIX C Second Order Analysis of Electrostatic Elements (Newtonian) . . 206 APPENDIX D Derivation of the Envelope Equation . . . . . . . . . . . . . . . . 238 APPENDIX E Paraxial Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . 245 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 249 vii LIST OF TABLES Table 1.1: List of codes used in this thesis. . . . . . . . . . . . . . . . . . . . . . . . 2 Table 4.1: x and y tunes of the test lattice as calculated by various codes. . . . . . . 59 Table 4.2: Higher order coefficients for the test lattice for the various codes, with and without fringe fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 Table 7.1: Electrode settings in volts. Each column refers to an electrode specified in Fig. 7.9, each row is a bias configuration referred to in the text. . . . . . . 142 viii LIST OF FIGURES Figure 1.1: Initially parallel trajectories in a uniform magnetic field. . . . . . . . . . Figure 1.2: First order trajectories of multiple initially parallel particles through a uniform magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.3: Third order trajectories of multiple initially parallel particles through a uniform magnetic field. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.4: Linearly increasing field. . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.5: Phase space diagram of two adjacent simple harmonic oscillators. . . . . 3 4 5 5 7 Figure 1.6: Twiss parameters for a phase space ellipse. . . . . . . . . . . . . . . . . . 10 Figure 1.7: An intuitive illustration of resonance. . . . . . . . . . . . . . . . . . . . . 14 Figure 1.8: Integer resonance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Figure 1.9: Half-integer resonance. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 Figure 2.1: Simplified Schematic of COSY J¨ulich Synchrotron. . . . . . . . . . . . . 17 Figure 2.2: COSY Infinity tracking, no fringe fields. . . . . . . . . . . . . . . . . . . 18 Figure 2.3: Zgoubi tracking, no fringe fields. . . . . . . . . . . . . . . . . . . . . . . . 19 Figure 2.4: MAD8 tracking picture, no fringe fields. . . . . . . . . . . . . . . . . . . 19 Figure 2.5: COSY Infinity order 9 tracking, symplectic, no fringe fields, x-x’ phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . space. Figure 2.6: COSY Infinity order 9 tracking, symplectic, no fringe fields, y-y’ phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . space. 20 20 Figure 2.7: COSY INFINITY default fringe field profile. . . . . . . . . . . . . . . . . 22 Figure 2.8: Zgoubi quadrupole fringe field profile (from [24]). . . . . . . . . . . . . . 23 Figure 2.9: Telescope segments for COSY J¨ulich. . . . . . . . . . . . . . . . . . . . . 24 Figure 2.10: Magnification requirements for COSY J¨ulich telescope segments. . . . . . 25 ix Figure 2.11: 9th order x-x’ tracking, nonlinear response. . . . . . . . . . . . . . . . . . 28 Figure 2.12: 9th order y-y’ tracking, nonlinear response and beam loss starting at . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . around 6 cm. 29 Figure 2.13: Detail of region where beam begins to deteriorate (x=6-6.5 cm.) . . . . . 30 Figure 2.14: Zgoubi xa tracking with fringe fields. . . . . . . . . . . . . . . . . . . . . 31 Figure 2.15: Zgoubi yb tracking with fringe fields. . . . . . . . . . . . . . . . . . . . . 32 Figure 2.16: MAD8 phase space in x-x’ using fringe field model outlined. . . . . . . . 35 Figure 2.17: MAD8 phase space in y-y’ using fringe field model outlined. . . . . . . . 36 Figure 2.18: COSY INFINITY nonsymplectic tracking. . . . . . . . . . . . . . . . . . 39 Figure 2.19: COSY INFINITY symplectic tracking. . . . . . . . . . . . . . . . . . . . 40 Figure 2.20: Representation of Thick Quadrupoles by Thin Lenses, SSC-N-33, Aug ’85 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [30]. 42 Figure 2.21: MAD8 tracking with thick element nonsymplectic tracking. . . . . . . . . 44 Figure 2.22: MAD8 tracking with thin lens elements giving symplectic tracking. . . . 44 Figure 2.23: ZGOUBI tracking with thick element nonsymplectic tracking. . . . . . . 45 Figure 3.1: The High Energy Storage Ring (HESR). . . . . . . . . . . . . . . . . . . 46 Figure 3.2: Relative error in (x—x) as a function of slices for MADX bend element. . 50 Figure 3.3: Relative error in (x—a) as a function of slices for MADX bend element. . 50 Figure 3.4: Relative error in (a—x) as a function of slices for MADX bend element. . 51 Figure 3.5: Relative error in (x—x) as a function of slices for MADX quadrupole. . . 51 Figure 3.6: Relative error in (x—a) as a function of slices for MADX quadrupole. . . 52 Figure 3.7: Relative error in (a—x) as a function of slices for MADX quadrupole. . . 52 Figure 3.8: 9th order symplectic tracking, 30k turns, with no fringe fields, 3 Gev . . . . . . . . . . . . . . . . . . . . . . . antiprotons (5 secs CPU time). 54 x Figure 3.9: 9th order symplectic tracking, 30k turns, with fringe fields. . . . . . . . . 54 Figure 3.10: 9th order symplectic tracking, 30k turns, with fringe fields, near dynamic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . aperture. Figure 3.11: MAD8 tracking, 128 slices, 30,000 turns, 17 minutes CPU time. No non- linearities observed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 4.1: Parallel-faced dipole with corresponding sector magnet of identical deflec- . . . . . . . . . . . . . . . . . . . . . . . . . . tion overlaid for reference. 55 55 58 Figure 4.2: A parallel-faced dipole is equivalent to a drift of length R sin φ. . . . . . 59 Figure 4.3: The effect of reducing aperture on several nonlinear map elements. . . . 60 Figure 4.4: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (ignoring fringe fields). . . . . . . . Figure 4.5: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (with fringe fields). . . . . . . . . . 66 66 Figure 4.6: ZGOUBI y-b tracking pictures dominated by the nonsymplectic behavior. 67 Figure 4.7: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (with fringe fields). . . . . . . . . . 68 Figure 4.8: MAD8 tracking pictures with fringe fields. . . . . . . . . . . . . . . . . . 69 Figure 5.1: (x,x’) coordinate of an off-axis cosine-like trajectory in a dipole. . . . . . 72 Figure 5.2: (x,a) coordinate tracking by Runge-Kutta of a cosine-like trajectory through a dipole at 7-70% radius (courtesy Eremey Valetov, MSU). . . . . . . . . Figure 5.3: (x,a) coordinate tracking by COSY INFINITY of a cosine-like trajectory in a dipole from 7-70% radius (courtesy Eremey Valetov, MSU). . . . . . Figure 5.4: (x,a) ZGOUBI tracking of a cosine-like trajectory through a dipole at 7-70% radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.5: ZGOUBI locations of particles after first 8 turns. Each branch is a suc- cessive turn. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.6: (x,a) MAD8 tracking of a cosine-like trajectory through a dipole at 7-70% radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 74 75 75 76 xi Figure 5.7: (x,a) MADX tracking of a cosine-like trajectory through a dipole at 7-70% radius. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 5.8: MADX locations of particles after first 4 turns. Each branch is a successive turn. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.1: x as a function of deflection angle. Figure 6.2: Plot of x − (x|x)x0 = (x|xx)x2 0. . . . . . . . . . . . . . . . . . . . . . . . Figure 6.3: Analytically calculated value of (x|xx). . . . . . . . . . . . . . . . . . . . Figure 6.4: x = (x|a)a0 + (x|aa)a2 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6.5: RK result, (x|aa) = x − (x|a)a0. . . . . . . . . . . . . . . . . . . . . . . . 76 77 96 97 97 98 99 Figure 6.6: Analytic solution of (x—aa). . . . . . . . . . . . . . . . . . . . . . . . . . 100 Figure 6.7: RK calculation of (x|ax). . . . . . . . . . . . . . . . . . . . . . . . . . . . 101 Figure 6.8: Analytic solution of (x|xa). . . . . . . . . . . . . . . . . . . . . . . . . . 102 . . . . . . . . . . . . . . . 103 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 Figure 6.9: RK plot of a as a function of deflection angle. Figure 6.10: RK plot of (a|aa). Figure 6.11: (a|aa) by taking the derivative of (x|aa). . . . . . . . . . . . . . . . . . . 105 Figure 6.12: RK calculation of (a|xx). . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 6.13: Analytic calculation of (a|aa). . . . . . . . . . . . . . . . . . . . . . . . . 106 Figure 6.14: RK calculation of (a|xa). . . . . . . . . . . . . . . . . . . . . . . . . . . . 107 Figure 6.15: Analytic Solution of (a|xa). . . . . . . . . . . . . . . . . . . . . . . . . . 107 Figure 6.16: RK (A—AA). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 Figure 6.17: Analytic A—AA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Figure 6.18: RK (A—XX). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 Figure 6.19: Analytic (A—XX). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Figure 7.1: Thomson’s device (1913). . . . . . . . . . . . . . . . . . . . . . . . . . . 122 xii Figure 7.2: Aston’s device (1920s). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 Figure 7.3: A magnetic mass spectrograph requires a physically large device to ac- commodate a wide mass range. A mass spectrometer requires a large field range for a wide mass range, and low hysteresis for a sample rate. . 124 Figure 7.4: Time of flight detector with long path length. . . . . . . . . . . . . . . . 124 Figure 7.5: Time of Flight detector with multiple reflections. . . . . . . . . . . . . . 125 Figure 7.6: ORISS during construction (Drawing courtesy LSU). . . . . . . . . . . . 126 Figure 7.7: Oriss within vacuum jacket (Drawing courtesy LSU). . . . . . . . . . . . 126 Figure 7.8: Oriss showing pertinent components. . . . . . . . . . . . . . . . . . . . . 126 Figure 7.9: Electrostatic mirror geometry showing the numbering of each electrode. . 127 Figure 7.10: Mirror electrode geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Figure 7.11: Conical lens geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 Figure 7.12: On-axis potential of ring of charge. . . . . . . . . . . . . . . . . . . . . . 133 Figure 7.13: On axis field of a ring of charge. . . . . . . . . . . . . . . . . . . . . . . . 134 Figure 7.14: Axial field off axis to second order for various values or r. . . . . . . . . 135 Figure 7.15: Off-axis electric field Er of a ring of charge. . . . . . . . . . . . . . . . . 136 Figure 7.16: Trajectory of a particle revealing both focusing and defocusing of a ring of charge. z=0 m is the center of the device, z=-.6 m is the location of the mirror. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136 Figure 7.17: Radial orbit of off-axis particle. z = 0 m is the center of the device, z=-.6 m is the location of the mirror. . . . . . . . . . . . . . . . . . . . . 137 Figure 7.18: Envelope of a stable orbit between two rings of charge. The mirrors are located at ±.06 m, and the particle starts at z = 0 m,r = .003 m. . . . . 137 Figure 7.19: Transverse phase space stability. . . . . . . . . . . . . . . . . . . . . . . . 138 Figure 7.20: Potential of a ring of charge (x is radial, y is axial). . . . . . . . . . . . . 138 Figure 7.21: Example of unstable orbit that would be lost. . . . . . . . . . . . . . . . 139 xiii Figure 7.22: Plot showing agreement between 238U+ ion repulsion in Warp and analytic solution via Coulomb’s law. . . . . . . . . . . . . . . . . . . . . . . . . . 144 Figure 7.23: Plot of bunch centroid separation (lower curve) and RMSz (upper curve) vs. time. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 Figure 7.24: Typical potential with no additional transverse focusing. . . . . . . . . . 146 Figure 7.25: Optimal energy focusing. . . . . . . . . . . . . . . . . . . . . . . . . . . . 147 Figure 7.26: Energy isochronous, with no additional transverse focusing. . . . . . . . 148 Figure 7.27: Same potential with slight transverse focusing. . . . . . . . . . . . . . . . 149 Figure 7.28: Point to parallel. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 Figure 7.29: Mass separation. Solid curve is distance between centroids of ion species, dashed curve is RMS radius. . . . . . . . . . . . . . . . . . . . . . . . . . 150 Figure 7.30: Initial configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 Figure 7.31: Separation after 1 reflection. . . . . . . . . . . . . . . . . . . . . . . . . . 151 Figure 7.32: Bunch centroid separation compared with bunch RMS radius, 10000 par- ticles, electrode settings from row A of table 7.1. Due to space charge, expansion of the bunch dominates any bunch separation. . . . . . . . . . 152 Figure 7.33: Doubling the mirror voltage to 200V as in row B of 7.1 increases the bunch separation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152 Figure 7.34: Mirror voltage of 400V, row C of table 7.1. . . . . . . . . . . . . . . . . . 153 Figure 7.35: Mirror voltage 800V, row D of table 7.1. . . . . . . . . . . . . . . . . . . 153 Figure 7.36: Electrode settings row D, now 20000 particles. . . . . . . . . . . . . . . . 154 Figure 7.37: Electrode settings row D, now 40000 particles. . . . . . . . . . . . . . . . 154 Figure 7.38: Electrode settings row D, now 80000 particles. . . . . . . . . . . . . . . . 155 Figure 7.39: Minimal potential configuration. . . . . . . . . . . . . . . . . . . . . . . . 156 Figure 7.40: Energy isochronous condition as a function of beam energy. . . . . . . . 158 xiv Figure 7.41: Energy Focusing : Focusing strength F= 0kV (n=10e3, E=3.4kV, V=10kV, row G). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158 Figure 7.42: Transverse trajectory : Focusing strength E3= 0 kV (n=10e3, E=3.4kV, V=10kV, row G). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159 Figure 7.43: Transverse and Longitudinal focusing E3=1 kV, row H. . . . . . . . . . . 160 Figure 7.44: Transverse and Longitudinal focusing E3=2 kV, row I. . . . . . . . . . . 161 Figure 7.45: Transverse and Longitudinal focusing E3=3 kV, row J. . . . . . . . . . . 162 Figure 7.46: Transverse and Longitudinal focusing E3=4 kV, row K. . . . . . . . . . . 163 Figure 7.47: Focusing strength E3= 0 kV, with space charge, 106 particles, row G. . . 163 Figure 7.48: Focusing strength E3= 1 kV, with space charge, 106 particles, row H. . . 164 Figure 7.49: Focusing strength E3= 4 kV, with space charge, 106 particles, row K. . . 164 Figure 7.50: Focusing strength E3= 8 kV, with space charge, 106 particles, row N. . . 164 Figure 7.51: Focusing strength E3= 10 kV, with space charge, 106 particles, row P. . 165 Figure 7.52: Tracking picture of 106 particles over 20 orbits, row P. . . . . . . . . . . 165 Figure 7.53: Initial configuration, 106 particles at z = 0, 50%U 238, 50%N 237. . . . . . 166 Figure 7.54: Particle bunch reaches the turning point, note the spread due to optics and space charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 Figure 7.55: Bunch returns to center, focusing element constrains the expansion due to space charge. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168 Figure 7.56: Successful separation after several reflections, with intense space charge. 169 Figure 7.57: Particle bunch initial condition. Note the larger number of tracer partices. 170 Figure 7.58: Bunch reaches the turning point, but the optics has distorted the shape of the bunch. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170 Figure 7.59: In an untuned system with space charge, the particles do not separate. . 171 Figure 7.60: Schematic of FRIB ReA (Reaccelerator) showing the Beam Cooler Buncher.172 xv Figure 7.61: Detailed diagram of ReA platform. The BCB (not shown) is directly below EBIT. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173 Figure 7.62: ECB longitudinal kinetic energy distribution. . . . . . . . . . . . . . . . 174 Figure 7.63: ECB transverse distribution. . . . . . . . . . . . . . . . . . . . . . . . . . 175 Figure 7.64: BCB equivalent distribution with no emittance. . . . . . . . . . . . . . . 175 Figure 7.65: Same as Fig. 7.64, but including emittance. . . . . . . . . . . . . . . . . 176 Figure 7.66: Same as Fig. 7.65, including space charge. . . . . . . . . . . . . . . . . . 176 Figure 7.67: with space charge (100k particles), with velocity dispersion, 20 reflections. 177 Figure 7.68: Typical potential of ISOLTRAP MRTOF. . . . . . . . . . . . . . . . . . 179 Figure 7.69: IM potential implemented in ORISS using the electrode settings in row G of table 7.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180 Figure 7.70: Time of flight spectrum of A = 53 nuclides delivered from ISOLDE (53Cr+,53Ca+, and the reference ion 39K+. . . . . . . . . . . . . . . . . 180 Figure 7.71: TOF spectrum of 53Cr and 53Ca, 5000 ions. . . . . . . . . . . . . . . . . 181 Figure A.1: Ideal dipole between two oppositely directed, unbounded current sheets. 191 Figure A.2: Field of oppositely directled current slices. . . . . . . . . . . . . . . . . . 191 Figure A.3: Top view of exit fringe field. . . . . . . . . . . . . . . . . . . . . . . . . . 194 Figure B.1: An excerpt from Brown’s seminal paper [5], showing the traditional coor- dinate system. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198 Figure B.2: Excerpt from Wiedemann . . . . . . . . . . . . . . . . . . . . . . . . . . 203 Figure D.1: Rectangular region of phase space. . . . . . . . . . . . . . . . . . . . . . 239 xvi Chapter 1 Introduction This dissertation is a document of how I came to learn beam physics. The novice quickly realizes that the subject is not “distilled”; that is to say there are very few textbooks that are both comprehensive and accessible to the beginner. I personally acquired many books on the subject, and am constantly referring to all of them to gain a better understanding of a particular topic. One strategy I employed was to learn to use as many different beam physics codes as I could lay my hands on. Each code lends a different perspective to what it means to “do beam physics”. It is my hope that this dissertation distills some of the knowledge that I have acquired in a manner that benefits the reader. A major portion of this thesis is the analysis of physical problems with several simulation codes, and comparing the results. For ease of reference, table 1.1 lists all of the codes which contributed to this work. 1 Code Type Language Source COSY Infinity Taylor series based transfer map Fortran with propri- etary scripting lan- guage cosyinfinity.org ZGOUBI Ray tracing fortran zgoubi.sourceforge.io MAD8 MADX 2nd order map transfer Fortran mad8.web.cern.ch 2nd order transfer map, symplectic in- tegration C,C++,Fortran madx.web.cern.ch GICOSY 5th order map transfer Fortran web-docs.gsi.de/ we- ick/gicosy/ Warp Particle in cell (PIC) Python, Fortran warp.lbl.gov LISE++ various C++ lise.nscl.msu.edu Table 1.1: List of codes used in this thesis. 2 1.1 Orbits Beam physics is the study of charged particles moving through electromagnetic fields. The simplest magnetic field is a uniform magnetic field, and as any schoolchild knows, the path of a charged particle in a uniform magnetic field is circular in a plane transverse to the field. But even this simple case has subtle properties - imagine two identical charged particles launched in parallel through a uniform magnetic field. Their trajectories will be as shown in Fig. 1.1. Figure 1.1: Initially parallel trajectories in a uniform magnetic field. The two particles move with initially parallel trajectories, but their paths will eventally cross. It’s easy to work out a trigonometric formula for their separation x as a function of angle θ: x = x0 cos θ + (cid:113) R2 − x2 0 sin2 θ − R. To first order, this separation is simply x ≈ x0 cos θ, and the distance between the two par- ticles varies harmonically. In this approximation, their trajectories are seen to converge near 3 θ x X0 θ = π/2. In fact, if we plot multiple such parallel trajectories in first order approximation, we see the following pattern (Fig. 1.2): Figure 1.2: First order trajectories of multiple initially parallel particles through a uniform magnetic field. This convergence of all trajectories is called “weak focusing”. Much of the language of optics has been co-opted by the beam physics community to the point where beam physics is often referred to as “Charged Particle Optics”. Such focusing is only approximately perfect. If we plot the same trajectories to third order, we find that the convergence is not to a single point, but is subject to aberrations similar to those found in optics (Fig. 1.3). 4 DIPOLE TEST 00 BEAM PLOT XX.XX.XXX YY:YY:YY X-MAX 1.400 m Y-MAX 0.280 m 0.500 m SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC Figure 1.3: Third order trajectories of multiple initially parallel particles through a uniform magnetic field. A slightly more complicated field is a single component, linearly increasing transverse magnetic field (Fig. 1.4). Figure 1.4: Linearly increasing field. This type of field also gives rise to simple harmonic motion. The simple harmonic motion seen in first order beam physics is often referred to as “Betatron motion” since the Betatron was the first machine for which this motion was studied in depth. The most common fields used in beam dynamics are constant fields (“dipoles”) and linearly increasing fields (“quadrupoles”). Since these elementary fields give rise to simple harmonic motion, it is seen that much of beam physics is nothing more than the study of harmonic oscillation! Simple harmonic motion (SHM) is studied most readily in phase space, and beam physics is no exception. 5 DIPOLE TEST 00 BEAM PLOT XX.XX.XXX YY:YY:YY X-MAX 1.400 m Y-MAX 0.280 m 0.500 m SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC SECTOR MAGNETIC Let us begin with the equations of SHM in postion, angle form: x(cid:48)(0) ω x(t) = x(0) cos(ωt) + x(cid:48)(t) = −ωx(0) sin(ωt) + x(cid:48)(0) cos(ωt). sin(ωt) Squaring and adding these equations yields: x2(0) + x(cid:48)2(0) ω2 = x2(t) + x(cid:48)2(t) ω2 = constant, which is usually rewritten in the form ωx2 + 1 ω x(cid:48)2 = (Area of phase space ellipse)/π ≡ “emittance ”(). Geometrically, this is the equation of an ellipse. The quantity , which is the area of said ellipse divided by π, is a constant of the motion. As a particle travels through a system consisting of various magnetic elements, the restoring force acting on the particle varies rapidly. The particle will be traveling through a given linear field, executing SHM, and then suddenly pass into another element with a different restoring force. What is the general behavior of the particle? Here is a simple experiment performed with the software COSY Infinity [34] - we set up a system consisting of two nonoverlapping linear elements with different field strengths. Then we track a particle through the system many times and plot the phase space coordinates of the particle. The result is the following x-x’ phase space diagram (Fig. 1.5): 6 Figure 1.5: Phase space diagram of two adjacent simple harmonic oscillators. It is perhaps counterintuitive that two dissimilar harmonic oscillators acting sequentially will give rise to a simple phase space ellipse. However, this “hybrid” ellipse has the inter- esting property that at the turning points (the points of maximum x), there is a nonzero momentum. This is not the case for a single harmonic oscillator, for which the turning points are characterized by zero momentum. Because an ellipse has an eccentricity, it can also have an orientation (a circle does not have an orientation). What causes the orientation of an ellipse? To see why this is so, we need to investigate a more general form of the harmonic oscillator equation: ¨x + k(t)x = 0. This is called Hill’s equation, and the restoring force is now a function of time. It is intuitively easy to visualize a qualitative solution of this equation. One can imaging a spring executing SHM and then the restoring force of the spring is gradually changed over time. Hill’s equation admits a solution of the form x(t) = u(t)ev(t), where u(t) is called the “amplitude function” and v(t) is called the “phase function”. A 7 special case of Hill’s equation arises when the restoring force k(t) is a periodic function of time [1]. In this case the amplitude function u(t) is also periodic, with the same periodicity as k(t). This is called Floquet’s theorem. Solutions of Hill’s equation are typically written in the form x(s) =(cid:112)β(s) cos(ψ(s) + φ) (cid:114)  β(s) x(cid:48)(s) = (α(s) cos(ψ(s) + φ) + sin(ψ(s) + φ)). The amplitude function u(t) has been replaced with a new function called the “beta function” β(s), and a scaling factor . These equations are not written with time as the independent variable. Since the dynamics occurs while particles are traveling along the circumference of the storage ring, we use the variable s to simply mean the arc length along the path of the reference particle following the centerline of the machine. Also, in the equation for the derivative, we have defined the variable α = −β 2 . If we consider s = 0, then ψ(s) = 0 and we derive the following two relations: x(0)(cid:112)β(0) (cid:114) (cid:20) β(0)  cos φ = sin φ = (cid:21) x(cid:48)(0) − α(0)x(0) β(0) 8 Squaring these and adding allows us to eliminate φ. This yields the following two equations: (cos ψ(s) + α(0) sin ψ(s))x(0) +(cid:112)β(s)β(0) sin ψ(s)x(cid:48)(0) x(s) = β(s) β(0) x(cid:48)(s) = −(1 − α(s)α(0) sin ψ(s) + (α(s) − α(0)) cos ψ(s) x(0) (cid:112)β(s)β(0) (cid:115) (cid:115) + (cos ψ(s) − α(s) sin(ψ)x(cid:48)(0) β(0) β(s) These can be written more conveniently in matrix notation: x(cid:48)(s)  x(s)  = (cid:112)β(s)β(0) sin ψ(s) (cid:114) β(s) (cos ψ(s) − α(s) sin(ψ) β(0)   x(0) x(cid:48)(0)   (cid:114) − (1−α(s)α(0) sin ψ(s)+(α(s)−α(0)) cos ψ(s) β(s) β(0)(cos ψ(s) + α(0) sin ψ(s)) √ β(s)β(0) We are interested in the special case where the β function is periodic. (Recall that this is true in a storage ring by Floquet’s theorem, because the restoring force is periodic). Therefore β(s) = β(0), α(s) = α(0), and the equation of motion for the particle simplifies to:  x(s) x(cid:48)(s)  = cos ψ(s) + α sin ψ(s) −γ sin ψ(s)   x(0) x(cid:48)(0)  β sin ψ(s) cos ψ(s) − α sin ψ(s) We’ve introduced a third and final variable, γ = 1 + α2 β 9 The quantities α, β, and γ are referred to as the Twiss parameters for an ellipse. To see these in context, refer to the following diagram (Fig. 1.6) from Wiedemann [2]: Figure 1.6: Twiss parameters for a phase space ellipse. α is a measure of how the ellipse is rotated. If α = 0, the motion of the particle is simple harmonic motion such as a single spring. This can be seen by setting α = 0 in equation (1.1). If we know the Twiss parameters at some initial point, we can calculate them at any subsequent point provided we have the first order transfer matrix in hand. The equation m2 11 −2m11m12 m2 12 −m11m21 m11m22 + m12m21 m22m12 m2 21 −2m22m21 m2 22  =   β α γ  ,   β0 α0 γ0 where β0, α0 and γ0 are the initial values of the Twiss parameters and the mij are the elements of the first order transfer map of the system, describes the evolution of the Twiss parameters. For example , the transfer matrix of a quadrupole of length s and field strength 10 k = e p0 dB dx is given by [2]m11 m12  = m21 m22  cos s −√ √ k √ k k sin s  . k 1√ k √ sin s √ k cos s We find, for example, that the evolution of the beta function through the quadrupole is given by [2] β(s) = β0 cos2 s √ k − 2α0 1√ k √ k sin s √ cos s k + γ0 sin2 s √ k. 1.2 Stability To aid the intuition, we have exploited the similarity between constant and linear field elements and harmonic oscillators. Besides fields with a restoring force that drives the particle back toward equilibrium, we also find electromagnetic elements with repulsive fields corresonding to k(t) < 0. We find, therefore, that a typical system will consist of both restoring and repulsive elements. The question arises regarding the stability of the system. Will a particle traveling through the system remain in the system, or will its orbital amplitude grow without bound? This is the question of stability. Let the first order behavior of a system be given by the matrix M = a b  . c d If a particle passes through the system n times, the evolution of the coordinates of the 11 particle is given by  = M n x y  . x0 a0 If the coordinates of the particle remain bounded, the system is stable as n → ∞. The conditions of stability can be discovered by looking at the eigenvectors and eigenvalues of M . We express the initial coordinates of the particle in terms of the eigenvectors (cid:126)v1 and (cid:126)v2 of M : In this way, the motion of the particle through the system n times is given by a0 x0  = A(cid:126)v1 + B(cid:126)v2.  = AM n (cid:126)v1 + BM n (cid:126)v2 = Aλn x0 a0 M n 1 (cid:126)v1 + Bλn 2 (cid:126)v2. For the coordinates of the particle to remain bounded we see that λn 2 must remain bounded as n → ∞. The conditions under which this is so can be found by looking at the 1 and λn solutions of the eigenvalue characteristic equation [2], namely (cid:115)(cid:18)Trace(M) (cid:19)2 − 1. λ± = Trace(M ) 2 ± 2 We see that the eigenvalues depend upon the Trace of M . If Trace of M is greater than 2, λ± > 1 and therefore λn → ∞ as n → ∞. In this case the system is unstable. 12 On the other hand, if Trace(M ) is less than 2, then the quantity (cid:115)(cid:18)Trace(M) (cid:19)2 − 1 ± 2 is a pure imaginary, and the two eigenvalues are therefore conjugate. Since the determinant of M = 1, we also have λ+λ∗− = 1. Hence λ+λ∗− = ((cid:60)λ+ + (cid:61)λ−)((cid:60)λ+ − (cid:61)λ−) = ((cid:60)λ+)2 + ((cid:61)λ−)2 = 1 and the eigenvalues have modulus 1. This means that the phase space coordinates of the particle do not grow or shrink - they only rotate. Hence if |Tr(M )| < 2 the system is stable. Stability of the transfer matrix alone is not the whole picture. There is another form of stability. There are two independent forms of periodic motion in a storage ring: the betatron motion of the particle about the reference orbit, and the orbital motion of the particle around the ring. A simple analogy is a wobbly ceiling fan (Fig. 1.7). 13 Figure 1.7: An intuitive illustration of resonance. As the fan turns on its axis, we often observe an up/down motion of the individual blades. The blades oscillate vertically with a given frequency, while the fan rotates with an independent frequency. If these two motions are commensurate, the fan will begin to resonate and will wobble more and more violently. On the other hand, if these two motions are not commensurate, then on average there will not be any net transfer of energy from the orbital motion to the wobbling motion. The motion of the fan will not be smooth, but there will be no tendency for the fan to oscillate out of control. A similar situation exists in storage rings. The orbit of the particles around the ring corresponds to the rotation of the fan about the axis, while the betatron motion of the particles about the reference orbit corresponds to the wobbling of the fan blades about their equilibrium point. Clearly the important quantity is the commensurability between the orbital motion and the betatron motion. The ratio of betatron frequency to orbital frequency is called tune. If there is a single perturbation in the ring, say at θ = 0, then resonance will only occur if the betatron oscillation is an integer multiple of the orbital motion. This is called an integer resonance (Fig. 1.8). 14 Figure 1.8: Integer resonance. In a similar way, if there are two diametrically opposed perturbations in the ring, say at θ = 0 and θ = π/2, this is called a half-integer resonance (Fig. 1.9). Figure 1.9: Half-integer resonance. 15 Chapter 2 Fringe Field Analysis at COSY J¨ulich Presently there is a great deal of interest in the use of storage rings to search for an electric dipole moment (EDM) in hadrons [6]. This requires utilizing the storage ring as a precision measuring device [7]. A significant challenge in understanding the detailed behavior of storage rings comes from a careful study of the axial variation of optical focusing elements [8]. Various tracking codes differ in their ability to model such behavior. Some codes model elements as posessing a uniform field within, with zero field outside (hard edge) while others allow the fields to physically vary along the element. 2.1 The COSY Storage Ring A major player in the search for hadron EDMs is the COoler SYnchrotron (COSY) at Forschungzentrum J¨ulich. We modeled a simplified version of this lattice with three popular particle tracking codes - COSY Infinity [11], Zgoubi [10] and MAD8 [9]. To establish a baseline, we began with a simplified hard-edge model which only incorpo- rated bending and focusing elements. There are 16 sextupoles present in the actual lattice which were not included in the model. A schematic of the lattice can be found in Fig. 2.1. 16 Figure 2.1: Simplified Schematic of COSY J¨ulich Synchrotron. After implementing the storage ring elements into the three different codes, we confirm that the first order 4 × 4 transfer matrices of a single orbit in the ring are identical: COSY    0.000 0.00 −0.977 −1.078 0.035 −0.984 0.000 0.00 0.000 −0.517 −10.90 .0652 −0.558 0.000 0.000 0.000 0.000 0.00 −0.977 −1.078 .0352 −0.984 0.000 0.00 0.000 −0.518 −10.90 .0652 −0.558 0.000 0.000 0.000 0.000 0.00 −0.977 −1.079 .0352 −0.984 0.000 0.00 0.000 −0.518 −10.90 .0652 −0.558 0.000 0.000 0.000 17 Infinity: Zgoubi MAD8    . This establishes a baseline agreement between the three codes COSY Infinity, Zgoubi, and MAD8. Next we compare the computed dynamic aperture of the system as computed by the three codes. Figure 2.2: COSY Infinity tracking, no fringe fields. Here we have initial particles released at 1 cm intervals from the reference orbit and tracked for 1000 turns. There is no limit on the dynamic aperture for this case (the beampipe at COSY J¨ulich is only 4.5 cm). The dynamic aperture as determined by ZGOUBI for the same lattice is shown in Fig. 2.3: 18 Figure 2.3: Zgoubi tracking, no fringe fields. Here again, we have no apparent limits on the dynamic aperture. Finally, Fig. 2.4 shows the tracking portrait produced by MAD8 for the same lattice: Figure 2.4: MAD8 tracking picture, no fringe fields. All three codes show unlimited dynamic aperture. With COSY Infinity we can increase 19 the accuracy of the tracking. For example, here is tracking in both x (Fig. 2.5) and y (Fig. 2.6) directions using a 9th order transfer map: Figure 2.5: COSY Infinity order 9 tracking, symplectic, no fringe fields, x-x’ phase space. Figure 2.6: COSY Infinity order 9 tracking, symplectic, no fringe fields, y-y’ phase space. Again we find no hint of a limitation. Based on the results so far, one would conclude that the system is remarkably stable. However, we have only studied the behavior of the 20 system in the absence of fringe fields. Each element has been assumed to have an ideal field which begins and ends abruptly at the physical boundaries of the element. As shown in Appendix A, such elements are unphysical. 2.2 Adding Fringe Fields COSY INFINITY has the ability to globally “turn on” fringe fields with a single command [11]. Every element in the system is given a default set of fringe fields by the issue of the FR3 (fringe fields type 3) command. Each dipole and quadrupole can be modeled with these default fringe fields. The standard method of implementing the fringe fields associated with the FR3 command is by means of the Enge function [11]: By(s) = 1 + ea1+a2(s/D)+a3(s/D)2+a4(s/D)4+a5(s/D)5+a6(s/D)6 B0 A typical set of values for the as are a1 = 0.478959 a2 = 1.911289 a3 = −1.185953 a4 = 1.630554 a5 = −1.082657 a6 = 0.318111 These are the default values used by COSY INFINITY to model the fringe field of a dipole. 21 Using these values gives rise to the fringe field profile shown in Fig. 2.7: Figure 2.7: COSY INFINITY default fringe field profile. with a similar profile for quadrupoles. After we apply fringe fields to the elements of the storage ring, the first task is to check the single turn transfer matrix. Here is the transfer matrix for COSY J¨ulich with fringe fields turned on:  -0.9739104 1.954368 0.00 0.00 0.01832738 -1.063567  0.00 0.00 0.00 −0.7993542 −6.705731 0.05219644 −0.8131372 0.00 0.00 0.00 Notice that the motion in the horizontal direction is now unstable (|Trace| = 2.0374774 > 2). To confirm this result, we can turn fringe fields on in ZGOUBI and compare the transfer matrices. Implementation of fringe fields in Zgoubi is similar to that of COSY INFINITY in 22 that we need to specify the same Enge coefficients used by COSY Infinity [12]. In the case of Zgoubi each element must manually be given both entry and exit Enge coefficients. In addition, Zgoubi requires you to specify how accurately you wish to integrate through the fringe fields. The extent of the integration zone must be specified, how extended are you considering your fringe field, and how small should your step size be throughout the fringe field. These parameters are illustrated in Fig. 2.8. Figure 2.8: Zgoubi quadrupole fringe field profile (from [24]). Here is the corresponding single turn transfer matrix in Zgoubi:  -0.961831 2.08049 .02259699 -1.08985 0.00 0.00 0.00 0.00  0.00 0.00 0.00 0.00 −0.209283 17.4232 −.05576599 −0.135606 Note that (|Trace| = 2.051681 > 2) with rough agreement with the COSY Infinity results 23 given above, providing independent confirmation that this fringe field configuration renders the COSY J¨ulich storage ring unstable in the horizontal direction. 2.3 Fixing an Unstable Lattice Model Without a stable transfer matrix, tracking is not possible and we are unable to proceed. Therefore, we must find a way to return the lattice to stability. There are certain design principles which must not be violated, and also there are symmetries in the lattice which need to be respected (Fig. 2.9) Figure 2.9: Telescope segments for COSY J¨ulich. 24 Figure 2.10: Magnification requirements for COSY J¨ulich telescope segments. The quadrupoles with matching fill patterns must remain equal in strength. This limits the number of degrees of freedom. Both straight sections must also be telescopes with unit magnification (Fig. 2.10). Using the MINPACK LMDIF optimization algorithm built into COSY INFINITY [13], we find a set of quadrupole strengths that respect the symmetry of the system, have the proper magnification, and restore stability to the ring. Here are the scaling multipliers for the default design quadrupole amplitudes we found 25 for each of these quadrupole groups: λ1 = 1.022183439914847 λ2 = 1.019331666908882 λ3 = 1.009453575838502 λ4 = 1.019948781502361 λ5 = 0.999410894910437 λ6 = 1.019427238706680 λ7 = 1.031603150073639 λ8 = 1.013424638586861 This shows that stability retaining proper optics can be effected by an adjustment of only few percent in the quad strengths. With these adjustments, we are able to recreate the original design transfer matrix in both COSY INFINITY: −0.9774877 −1.078547 0.03521565 −0.9841743 0 0 0 0 0 0 0 0 −0.5176308 −10.90340 0.06520659 −0.5583641   , 26 and ZGOUBI:   . −0.9774877 −1.078548 0.03521565 −0.9841743 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 −0.5176308 −10.90340 0.06520659 −0.5583641 Note that these agree to a high precision. 2.4 Side by Side Comparison With these new quadrupole strengths in hand, we can now run side by side comparisons with Zgoubi, which also has full fringe field modeling for quadrupoles and dipoles. Fig. 2.11 and Fig. 2.12 show the results of tracking with COSY INFINITY. 27 Figure 2.11: 9th order x-x’ tracking, nonlinear response. 28 Figure 2.12: 9th order y-y’ tracking, nonlinear response and beam loss starting at around 6 cm. Notice the appearance of resonances in the x-x’ motion. The dynamic aperture is limited in this case. Fig. 2.13 is a closeup of the phase space region where the beam begins to deteriorate. 29 Figure 2.13: Detail of region where beam begins to deteriorate (x=6-6.5 cm.) Next we perform the same tracking experiment in Zgoubi (Fig. 2.14 and Fig. 2.15). Note the horizontal shift consistent with fringe field effects in the dipoles. The loss of beam at ≈ 6 cm predicted by COSY Infinity is also suggested, taking lack of symplecticity in Zgoubi into account. 30 Figure 2.14: Zgoubi xa tracking with fringe fields. 31 Figure 2.15: Zgoubi yb tracking with fringe fields. MAD8 also has provisions for fringe fields, but only in the thick bending elements. How- ever, we can add approximate fringe fields to the thin-element model using the “Fringe Integral (FINT)” techniques developed by Brown [5]. We must provide the value of FINT as a parameter to the software. FINT enters the equations of motion in the manner illustrated 32 by Eq. 2.1.  =   x a y b   x0 a0 y0 b0  +  g2 ρ cos2 β I1 0 0 0  where I1 = (cid:90) (cid:90) (B0 y − By)dσdσ 0 0 0 1 0 1 ρ tan β 1 0 0 1 0 0 0 0 − 1 ρ tan βv 1 and (cid:18) β − g ρ 1 + sin2 β cos3 β (cid:19) I2 with I2 = (cid:90) ∞ −∞ By(B0 − By) B2 0 dσ tan βv = tan (2.1) The value of the FINT computed by Mathematica I2 = 0.0369458. As a check, we can calculate the vertical portion of the transfer matrix by hand for a simple 15 degree dipole, 1 m radius, no edge angles, 9 cm full aperture with COSY Infinity default fringe fields: 1.008708359 0.26179939 .06681663 1.008708359  0.00 .066816630 33 Hand Calculation MAD8  0.96592583 0.25881905 −0.25881905 0.96592583 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.0087084 0.26179939 1.0087084  COSY Infinity “FR3” Zgoubi 0.9659126 0.2589177 −0.2588190 0.9659126 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.008727 0.2622809 0.00 0.06683583 1.008727  .977473 .258496 −.151348 .978166 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00351 .262081 0.00 .03506525 1.00566   We can compare the MAD8 thin element tracking to the COSY Infinity and Zgoubi tracking. MAD8 has a bending element, but it is not good for tracking. Better to use a thin element dipole, and an edge angle, and add the fringe fields. This model yields: Hand Calculation: MAD8  1.032769823 .036767837  1.8117336 1.032769823 1.0000000 1.8117333 0.00000000 0.00000000 0.69388939E − 17 1.0000000 0.00000000 0.00000000 0.00000000 0.00000000 0.96683044 0.00000000 1.8325957 0.00000000 −0.035599178 0.96683044  34 COSY Infinity 1.000000 0.0000000 0.0000000 0.0000000 1.000000 1.811733 0.0000000 0.0000000 0.0000000 0.9655335 −0.03696677 0.0000000 0.0000000 0.0000000 1.832596 0.9655335   Corresponding MAD8 are shown in Fig. 2.16 and Fig. 2.17. Even with the fringe fields added in, none of the aperture-limitiing nonlinear dynamics are visible. The tracking order of MAD8 is simply not high enough. Figure 2.16: MAD8 phase space in x-x’ using fringe field model outlined. 35 -0.100-0.0500.00.0500.100x (m)Table name = TRACKTABLECandidate, Fringe, Order 6, 1k turns, [1-9cm],& xaWin32 version 8.51/1520/08/14 14.41.59-0.020-0.015-0.010-0.0050.00.0050.0100.0150.020px/p0 Figure 2.17: MAD8 phase space in y-y’ using fringe field model outlined. 2.5 Symplectic vs. Nonsymplectic Tracking One aspect of the tracking we have overlooked thus far is the notion of symplecticity. First we present the mathematical background. Propagation of a particle through a Hamiltonian system in particle optical coordinates can be described as follows: (xf , af , yf , bf , tf , Ef ) = f (xi, ai, yi, bi, ti, Ei) where f is an invertible mapping called the flow of the system, and we denote x(cid:48) = a and y(cid:48) = b. It can be shown [3] that the flow of a Hamiltonian system is symplectic. What precisely does this mean? To say that the flow f is symplectic is to say that the Jacobian of f satisfies the condition Jac(f ) · J · Jac(f )t = J 36 -0.100-0.0500.00.0500.100y (m)Table name = TRACKTABLECandidate, Fringe, Order 6, 1k turns, [1-9cm], ybWin32 version 8.51/1520/08/14 14.43.50-0.020-0.015-0.010-0.0050.00.0050.0100.0150.020py/p0 We can use this fact to derive many useful properties. As a simple example, consider the free particle Hamiltonian, H = p2/2m. In this case, Hamilton’s equations of motion are ˙x = ∂H ∂p = p m , ˙p = −∂H ∂x = 0 The solution for this system can be expressed in matrix form: x(t)  = 1 t/m  x0  p(t) 0 1 p0 The matrix is the flow of the system, and it is symplectic:  1   0 0 1 −1 0  1 t/m  =  0 0 1 1 −1 0  t/m 1 This is true in general. The Jacobian matrix for Hamiltonian flow in one coordinate and its conjugate momentum is  ∂Q ∂p ∂P ∂p . ∂q ∂P ∂q  ∂Q  ∂Q ∂q ∂Q ∂p Imposing the requirement of symplecticity gives  ∂Q ∂q ∂Q ∂p   0 1 −1 0   ∂Q ∂q ∂P ∂q ∂P ∂q ∂P ∂p  = ∂Q ∂p ∂P ∂p  =  0 1 −1 0  ∂Q ∂p ∂Q ∂p ∂Q ∂q ∂Q ∂q ∂Q ∂q ∂Q ∂p ∂P ∂q ∂p − ∂P ∂p − ∂P ∂p ∂P ∂P ∂q ∂q − ∂P ∂q − ∂P ∂p ∂P The diagonal entries are trivial identities, but the off-diagonal entries yield many interesting results. 37 First, we expand the off diagonal entries of the Jacobian in a Taylor series. (cid:88) (cid:88) Let Q = (Q|qapb)qapb and P = (P|qcpd)qcpd a,b c,d Taking the required partial derivatives a,b c,d (Q|qapb)aqa−1pb (cid:88) (cid:88) (P|qcpd)qcdpd−1 (cid:88) (cid:88) C,D (P|qC pD)CqC−1pD (Q|qApB)qABpB−1 A,B ∂Q ∂q ∂P ∂p ∂P ∂q ∂Q ∂p = = = = and forming the required products gives (cid:88) (Q|qapb)(P|qcpd)adqa+c−1pb+d−1 − (cid:88) (P|qC pD)(Q|qApB )BCqA+C−1pB+D−1 = 1. a,b,c,d A,B,C,D Equating powers of q and p in this expression gives many conditions imposed by symplec- ticity. For example, for the q0p0 terms then we find (Q|p)(P|p) − (P|q)(Q|p) = 1 which is the well-known result that the determinant of a Hamiltonian is unity. But looking at higher orders, we find a specific relation of the form 2(Q|q2)(P|p) − 2(P|q2)(Q|p) + (Q|q)(P|qp) − (P|q)(Q|qp) = 0. 38 The correctness of this result is confirmed by simulations with COSY INFINITY. If a tracking code does not preserve symplecticity, this means that phase space area is not conserved as required by Hamiltonian mechanics. This is indicated by a “spiraling” of the particle orbits. For example, Fig. 2.18 shows nonsymplectic tracking by COSY INFINITY through the COSY J¨uelich storage ring at ForschungZentrum J¨ulich. The characteristic spiralling due to the nonphysical phase space damping apparent and large for higher ampli- tudes. For comparison, symplectic tracking through the same configuration shows that the spiraling is gone, and the orbit is a perfect ellipse (Fig. 2.19). Figure 2.18: COSY INFINITY nonsymplectic tracking. 39 Figure 2.19: COSY INFINITY symplectic tracking. Codes like MAD8 [28] and MADX [27] use thin elements to preserve symplecticity. A quadrupole can be approximated by a thin element. The bending nature of the quadrupole is simply   1 √ k sin L −√ 0 k ≈ −LK 1. The thickness of the quadrupole is modeled by drifts of length √ k L 2 / √ L k 2 L 2 tan ≈ L 2 . Putting these elements together yields a thin element quadrupole, 1 L/2   1 0 1 0 −kL 1 0 1  1 L/2  =  , 1 − L2K/2 L(1 − L2K/4) −kL 1 − L2K/2 which is a symplectic matrix. 40 In fact, from the definition of symplecticity, we can make a very general statement about thin elements. Again, consider the Jacobian: ∂Q ∂q ∂P ∂q  ∂Q ∂p ∂P ∂p “Thin element” means simply that ∂Q/∂q = 1 and ∂Q/∂P = 0, in other words, there is no change in the position of the particle - only a change in the angle.  1 ∂P ∂q  0 1 We easily verify that this matrix is symplectic: 1 ∂P ∂q   0 0 1 1 −1 0   1 ∂P ∂q  = 0 1  0 1 −1 0  MAD8 does not use symplectic tracking by default. However, if we use thin elements, we can perform symplectic tracking with MAD8, up to errors in numerical precision. A quadrupole can be most accurately replaced with thin elements according to the following scheme developed by Ricard Talman [30] summarized in Fig. 2.20. Each row in the figure represents an increasingly accurate thin lens approximation for a thick quadrupole. 41 Figure 2.20: Representation of Thick Quadrupoles by Thin Lenses, SSC-N-33, Aug ’85 [30]. The simplest scheme is shown in the top row of the figure, whereby the quadrupole is simply replaced with a thin lens at the very center. The bulk of the quadrupole is replaced with drifts. The fourth option has been shown to be the most accurate. We replace each quadrupole of the lattice with this equivalent thin lens configuration. COSY J¨ulich also uses parallel-faced dipole magnets for bending. Because of the edge angles, these elements focus in the y-direction. Therefore, for thin-element tracking, every rectangular dipole must be replaced with a drift preceded and followed by two edge-focusing elements of the form 42 1 L/2 0 1 L/2 0  0 0 0 1 0 0   0 0 0 1 0 0 1 1 L/2 0 0 0 0 1 0 1 0 0 tan(ψ) 1   0 0 0 0 0 0 1 0 0  0 0 0 1 L/2 0 1 When these replacements are done correctly, the first-order thin-element model of a single- turn transfer map should match the original: Compare to the original design transfer matrix received from ForschungZentrum-J¨ulich:   −0.977 −1.053 0.00 .035 −0.985 0.00 0.00 0.00 .064 −0.610 0.00 0.00 0.00 0.00 −0.977 −1.078 0.035 −0.984 0.00 0.00 0.00 −0.518 −10.9 .065 −0.558 0.00 0.00 0.00   With standard (thick) elements, MAD8 is seen to be nonsymplectic (Fig. 2.21) with charac- teristic phase space damping. 43 Figure 2.21: MAD8 tracking with thick element nonsymplectic tracking. However, when we replace the lattice with thin elements, we can perform symplectic tracking. The resultant (symplectic) tracking pictures are shown in Fig. 2.22: Figure 2.22: MAD8 tracking with thin lens elements giving symplectic tracking. ZGOUBI does not have the capability to perform symplectic tracking. (Fig. 2.23) shows 44 -0.05-0.03-0.010.010.030.05x (m)Table name = TRACKTABLEBaseline, TRANSPORT, Order 2, 20k turns, [X,A]Win32 version 8.51/1518/06/14 07.21.22-0.010-0.008-0.006-0.004-0.0020.00.0020.0040.0060.0080.010px/p0 the phase space tracking picture from ZGOUBI for the COSY J¨uelich system. Note the unphysical damping error is limited. Figure 2.23: ZGOUBI tracking with thick element nonsymplectic tracking. 2.6 Summary The apparent instability in our model of the Cosy J¨uelich storage ring when fringe fields were taken into account was unexpected, given that the facility has been in operation for decades. The first thought was that some error in the coding was the cause of this issue. It was not until the exact same result was obtained in a wholly independent manner that we began to suspect the design settings. This incident was the first of many demonstrating that simulating a system across multiple codes is a worthwhile investment of the time and effort involved. This work enabled more reliable, consistent modeling. 45 -.04-.020.00.020.04-.008-.006-.004-.0020.00.0020.0040.0060.008Zgoubi|Zpop 18-Aug-14 * ZGOUBI COSY Julich 00 Y’ (rad) vs. Y (m) Mi-ma H/V: -5.356E-02 5.328E-02/ -9.283E-03 9.203E-03 Part# 1- 40000 (*) ; Lmnt# 1; pass# 1- 2001; 2001 Chapter 3 The High Energy Storage Ring HESR (COSY Infinity, MAD8, MADX) Figure 3.1: The High Energy Storage Ring (HESR). The High Energy Storage Ring (HESR), part of FAIR at GSI in Darmstadt, Germany, is a facility dedicated to the field of high energy antiproton physics (Fig. 3.1) [14]. The HESR will be operated in two different modes - a high luminosity mode and a high resolution mode. The experimental requirements for HESR dictate its design characteristics. It must be a low impedance ring with large dynamic aperture, and it must be very stable in terms of beam parameters for high precision experiments. At the interaction point, the beam size must be 46 very small, which entails space charge considerations. Like COSY-J¨ulich, the HESR is designed as a racetrack shaped storage ring. But the circumference of the ring and number of elements is greatly increased. It has a magnetic bending power of 50 Tm and a circumference of 574 m, including two 132 m long straight sec- tions. The straight sections have large gaps for the installation of instrumentation, detectors, and electron cooling. For a simulation, we perform a side-by-side comparison of HESR simulations between MAD8 and MADX. The first step is to match the first-order transfer matrices. The lattice was ported into COSY INFINITY v9.1, MAD8 v8.51 Win32 and MAD-X 5.02.00 (64 bit, Linux). Here are the resulting first-order transfer matrices in the x− a directions (note a=x’ and b=y’): MAD8 v8.51 Win32: −1.0456753 −1.7327760 0.18033015 −0.65749692  MAD-X 5.02.00 (64 bit, Linux) −1.0456753 −1.7327760 0.18033015 −.65749692  COSY INFINITY v9.1 and discrepancy  −1.046011 −1.730853 .1803385 −.6576032   3.2093353e − 4 1.1110129e − 3 .46301816e − 5 1.6161722e − 4  Not surprisingly, the two versions of MAD agree. But there is a slight discrepancy between 47 these codes and COSY INFINITY. Can we account for this discrepancy? And which model is correct? Since the lattices in this simulation are composed strictly of dipoles and quadrupoles, we can look closely at the representation of these elements in the three codes. For a simple sector bend of 45◦, the transfer matrix is calculated analytically to be Hand Calculation cos (s/R) − 1 R sin (s/R)  =   . .98982144 4.1886283 −4.8353571e − 3 .98982144 R sin (s/R) cos (s/R) The same element in COSY infinity obtains COSY INFINITY and relative error  .98982144 4.1886283 −.48353571e − 2e − 2 .98982144   . 0 0 0 0 3.1 MAD8, MADX and Thin Elements Unlike MAD8, MADX does not work with “thick” elements. Instead, one takes a thick element and slices it into a variable number of thin elements. Here we slice the same 45◦ bending element into various sized slices: MADX (1 slice) and relative error  .98980413 4.1814717 −.48518298e − 2 .98980413   1.7490427e − 5 1.7085863e − 3 −3.4065529e − 3 1.7490427e − 5  48 MADX (128 slices) and relative error  .98982144 4.1886287 −.48353571e − 2 .98982144  0 −9.5496657e − 8  0 0 We can perform a similar analysis for a thick quadrupole. Here is the analytically calculated transfer matrix for the defocusing direction of a quadrupole:  cosh −√ (cid:17) (cid:16)√ (cid:16)√ ks (cid:17) k sinh k cosh ks Hand calculation and COSY Infinity  = (cid:17) 1.064308 .6128072  .2166274 1.064308 1√ k sinh (k) (cid:16)√ and the calculated transfer matrix for MADX with 1 and 128 slices: MADX and relative error (slices=1) 1.06363 .619089  .2121 1.06363 6.370336e − 4 −1.025086e − 2  2.089948e − 2 6.370336e − 4 MADX and relative error (slices=128) 1.0643076 .61280682  .21662735 1.0643076 9.3957791e − 8 6.527338239e − 7  9.3957791e − 8 0 Working with 128 slices is extreme, and computationally costly. Is there an optimum value? In figure 3.2we plot the relative error vs the number of slices on a log-log plot for the various dipole and quadrupole slice approximations: 49 Figure 3.2: Relative error in (x—x) as a function of slices for MADX bend element. Figure 3.3: Relative error in (x—a) as a function of slices for MADX bend element. 50 Figure 3.4: Relative error in (a—x) as a function of slices for MADX bend element. Figure 3.5: Relative error in (x—x) as a function of slices for MADX quadrupole. 51 Figure 3.6: Relative error in (x—a) as a function of slices for MADX quadrupole. Figure 3.7: Relative error in (a—x) as a function of slices for MADX quadrupole. 52 3.2 Tracking and Fringe Fields Since computer run time was not a primary concern, we chose to run our simulations with 128 slices. We want to estimate the dynamic aperture of the system using both COSY Infinity and MADX. With COSY Infinity we can turn on the fringe fields. Here is the transfer matrix of the HESR lattice hard-edge elements: COSY Infinity FR0  −1.04611 −1.730853 −1.730853 −0.6576032  and here it is with fringe fields turned on: FR3  −1.061066 −1.462613 .1675836 −.7114451  Unlike the situation for COSY J¨ulich analyzed in Chapter 2, the lattice remains stable when fringe fields are taken into consideration. There is substantial change, but the system is still stable and ready for tracking. Fig. 3.8 shows the result of tracking the HESR field to 9th order, for 30,000 turns, with fringe fields turned off. Not surprisingly, the dynamic aperture is unlimited. Next we turn on fringe fields, to determine the effect on dynamic aperture. Fig. (3.9) shows tracking of the same system, this time with fringe fields turned on. A clear limitation of the dynamic aperture is present around at around 4.5 cm. This is the design aperture of the system. 53 Figure 3.8: 9th order symplectic tracking, 30k turns, with no fringe fields, 3 Gev antiprotons (5 secs CPU time). Figure 3.9: 9th order symplectic tracking, 30k turns, with fringe fields. Fig. 3.10 presents a detailed view of the region near the limit of the dyamic aperture. The nonlinear dynamics of the system is apparent. 54 Figure 3.10: 9th order symplectic tracking, 30k turns, with fringe fields, near dynamic aperture. Now we track the same system with MADX, representing the thick elements by 128 slices. Fig. 3.11 shows the resulting phase space. After 30k turns, and 17 minutes (≈ 20× longer run time), none of the nonlinearities which limit the dynamic aperture are visible. Figure 3.11: MAD8 tracking, 128 slices, 30,000 turns, 17 minutes CPU time. No nonlinear- ities observed. 55 3.3 Summary In this chapter we essentially repeated the same analysis for the GSI HESR as we did for Cosy J¨uelich, with satisfactory results. However, the comparison between MAD8 and its successor MADX revealed some interesting challenges.There are major functional changes between MAD8 and MADX, most notably the loss of any “thick” elements. The replacement of thick elements with thin elements is a nontrivial exercise, as shown by Talman [30]. In addition, the codebase of MADX is currently under heavy development, and this can lead to instabilities. Issues relating to the viability of the modeling can open issues on whether instabilities evident in code models properly represent the physical ring. 56 Chapter 4 Benchmarking COSY Infinity, ZGOUBI, MAD8, GICOSY Of utmost importance in this regard is the capability of tracking codes to treat all nonlinear effects arising from the detailed field distributions present in the system, not the least of which are fringe fields. In previously published work [17][18], (also Chapters 2 and 3 of this thesis), we performed parallel tests of various tracking codes in order to compare and contrast the results. In this study [57] , we continue this line of research and extend the scope to analyze dipoles with edge angles. 4.1 Modeling Dipoles With Edge Angles A common replacement for bending magnets of the sector variety are dipoles with edge angles. These differ from traditional sector magnets in that the reference orbit enters and exits at an angle with respect to the face of the magnet (Fig. 4.1). 57 Figure 4.1: Parallel-faced dipole with corresponding sector magnet of identical deflection overlaid for reference. This difference results in a change in the optics, especially with respect to fringe fields. The transfer map of a parallel-faced dipole must be independent of the original horizontal position xi since the motion is invariant with respect to change in this direction. This implies that many matrix elements vanish. On the other hand, under the presence of fringe fields, the motion depends on the vertical position y, and several new matrix elements appear. We will show that fringe fields in conjunction with edge angles dramatically limit the dynamic aperture. It can be shown [19] that the horizontal transfer map of a parallel-face dipole of deflection angle φ and deflection radius R is exactly a drift: 1 R sin φ  0 1 The “length of the drift”, given by (x|a) = R sin φ, requires some interpretation which can be made clear by referring to the following figure (4.2): 58 Figure 4.2: A parallel-faced dipole is equivalent to a drift of length R sin φ. This suggests a simple test to confirm the accuracy with which a code implements edge angles. It must always be the case that (x|x) = 1, even when fringe fields are taken into consideration. We chose to test this assertion with the codes COSY INFINITY v9.1 [34], ZGOUBI v6.0.1 [24], MAD8 v8.51 [28] and GICOSY [25]. The tracking lattice used consisted of four 90 degree dipoles of radius 10 meters, each having entrance and exit angles of 22.5 degrees. The dipoles have a 10 cm full aperture and are separated by 5 m drifts. Implementing this lattice with the aforementioned codes, we find the tunes given in Table 4.1. Code COSY INF. (no fringe fields) COSY INF. (with fringe fields) ZGOUBI (no fringe fields) ZGOUBI (with fringe fields) GICOSY (fringe) MAD8 (fringe) Tune y .8655516 .8534879 .86572914 .83195955 .853462867 .917089883 Tune x .8671905 .86714278 .86719261 .86818895 .867160222 .943335438 Table 4.1: x and y tunes of the test lattice as calculated by various codes. For those codes that support nonlinear elements, a comparison of map coefficients for the cases with and without fringe fields is given in Table 4.2. When fringe fields are taken 59 COSY INF. (w/o ff) -.3330669e-15 2.000000 .5 -.1110223e-15 -.5551115e-15 1.000000 .000000 .000000 .000000 COSY INF. (w/ ff) .13552497e-9 2.000008 .47113316 -.33329235e-8 .24161645e-8 1.000008 .5886810 .4861485 -2.1406605 ZGOUBI GICOSY (w/ ff) (w/ ff) -8.799e-2 -2.777411e-6 6.372e-2 1.999998 .4711359 1.16 -8.928427e-8 N/A -3.462375e-6 N/A .9999951 N/A -4.485e-2 .5886785 .244 .4861427 -2.140659 4.22 (x—xx) (x—xa) (x—aa) (x—xxx) (x—xxa) (x—xaa) (x—yy) (x—yb) (x—bb) Table 4.2: Higher order coefficients for the test lattice for the various codes, with and without fringe fields. into account, these coefficients are dependent upon the aperture of the elements. Starting with the initial aperture of 10 cm, we can reduce the aperture and watch the behavior of the various coefficients (Fig. 4.3). Figure 4.3: The effect of reducing aperture on several nonlinear map elements. The transfer matrix for a small 10 degree bend with 1 meter deflection radius as calculated 60      to 9th order by COSY INFINITY is 1.000000000000000 .1736481776669303  0 1.000000000000000 If we increase the deflection angle to 90 degrees, we get 1.000000000000000 1.000000000000000 −.2220446049250313e − 15 .9999999999999999 Fringe fields should not affect these outcomes. COSY INFINITY provides a convenient method for generating a default set of fringe fields which are representative of realistic values. These fringe fields are specified by giving the Enge coefficients [26]. Repeating the test with fringe fields we find that for the 10 degree deflection angle: .9999999999997435 .1712151523389203 −.1779412557232258e − 13 1.000000000000256 (4.1) (4.2) corresponding to a change of less than 10−10 in the (x|x) element. Increasing the angle to 90 degrees, we find:  1.000000000020302 .9422624094620093 .3590956618514349e − 12 .9999999999807103 with an equally small variation from the previous. From this we conclude that even under the most stringent of conditions, COSY INFINITY maintains its fidelity to sub-parts per billion. ZGOUBI also implements fringe fields, and can use the same Enge coefficients as 61 those specified in COSY INFINITY. We can therefore run the same tests and compare the outcomes. First, without fringe fields, here is the transfer matrix of a 10 degree bend angle with 1 meter deflection radius as calculated by ZGOUBI:  1.00000 .173648 −4.135581e − 12 1.00000 These are the reasonable values one would expect. Also, with a 90 degree deflection angle, ZGOUBI still performs well:  1.00001 1.00001 1.065148e − 9 .999992  But these results are neglecting fringe fields. With fringe fields enabled ZGOUBI calculates the following transfer matrix:  1.00067 .173290 −3.426946e − 12 1.00056 with a corresponding discrepancy of .67e − 4. Although this is not a large change for a single pass system, storage rings are intended to be used for many turns. Such a slight error in a single pass can cause substantial discrepancies in long term behavior. Furthermore, increasing the deflection angle to 90 degrees we find:    (4.3) (4.4)  1.01332 1.00273 9.364740e − 10 .986917 62 and now the error in the (x|x) term has increased to 1.332%. GICOSY can also be subjected to these tests, with the additional capability to perform 3rd order tracking for comparison with COSY INFINITY’s higher order tracking. Here is the first-order transfer map for a parallel-faced bending element with deflection angle of 10 degrees and a 1 meter deflection radius, using the same Enge coefficients as were used previously for COSY INFINITY and ZGOUBI (compare with (Eq. 4.1) and (Eq. 4.3)):    .9999999 .1729761 −2.814809e − 7 1.000000 1.000000 .1729761 −3.157162e − 8 1.000000    With GICOSY’s option to perform third order calculations we obtain: When the deflection angle is increased to 90 degrees, all else remaining the same, we get a noticeable deterioration in accuracy (compare with Eq. 4.2): .9999999 .9838112 −5.716222e − 8 1.000000 ZGOUBI has the ability to output second-order transfer maps. These maps are estimated by launching 61 particles at random within a region of phase space, and computing the resulting map. We can take this opportunity to compare these maps, computed by integration, with those of COSY INFINITY, which is a pure transfer map code. Here are the second and third order elements from COSY INFINITY for a 90 degree deflection angle, 1 meter deflection 63 radius, without fringe fields: (x|xx) = −.3330669e − 15 (x|xa) = 2.000000 (x|aa) = .5 (x|xxx) = −.1110223e − 15 (x|xxa) = −.5551115e − 15 (x|xaa) = 1.000000 (x|yy) = .000000 (x|yb) = .000000 (x|bb) = .000000 Notice the null values for the (x|yy), (x|yb), (x|bb) terms. Obviously a constant vertical magnetic field will not influence the vertical trajectory of a particle in the midplane. However, enabling fringe fields changes this as can be seen from the same transfer map, this time with 64 fringe fields enabled: (x|xx) = .13552497e − 9 (x|xa) = 2.000008 (x|aa) = .47113316 (x|xxx) = −.33329235e − 8 (x|xxa) = .24161645e − 8 (x|xaa) = 1.000008 (x|yy) = .5886810 (x|yb) = .4861485 (x|bb) = −2.1406605 and here are the corresponding elements from GICOSY: (x|xx) = −2.777411e − 6 (x|xa) = 1.999998 (x|aa) = .4711359 (x|xxx) = −8.928427e − 8 (x|xxa) = −3.462375e − 6 (x|xaa) = .9999951 (x|yy) = .5886785 (x|yb) = .4861427 (x|bb) = −2.140659 65 The effects of the (x|yy), (x|yb), and (x|bb) terms can be clearly demonstrated by comparing the tracking pictures without fringe fields (Fig. 4.4) and with fringe fields (Fig. 4.7): Figure 4.4: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (ignoring fringe fields). Figure 4.5: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (with fringe fields). The trajectories are at 1 cm increments from the reference orbit, out to 30 cm. The limitation on the dynamic aperture due solely to effect of fringe fields is readily apparent. ZGOUBI can also do tracking with fringe fields, which can be specified with the same Enge coefficients as those used by COSY INFINITY (Fig. 4.6). Because of the nonsymplectic behavior of the tracking, the y-b phase space orbits are seen to rapidly decay, making it impossible to see any of the effects due to the fringe fields. The nonlinear elements available 66 for output from ZGOUBI are: (x|xx) = −9.437e − 10 (x|xa) = .993 (x|aa) = 1.26 (x|yy) = .814 (x|yb) = −.253 (x|bb) = .747 This feature of ZGOUBI does not seem to give good agreement with the transfer maps calculated by COSY INFINITY and GICOSY. There are several parameters which impact the accuracy of fringe field modeling in ZGOUBI and this issue is a subject of ongoing investigation. Nevertheless, can ZGOUBI detect the limitation in dynamic aperture due to the fringe fields? Fig. 4.6 is the tracking picture in y-b phase space. Figure 4.6: ZGOUBI y-b tracking pictures dominated by the nonsymplectic behavior. For comparison purposes, Fig. 4.7 is a tracking picture with COSY INFINITY running at an artificially low order of 3 and with symplectic tracking disabled. 67 Figure 4.7: y-b tracking in a simple storage ring composed of four 90 degree sector magnets with 22.5 degree edge angles (with fringe fields). MAD8 v8.51 was also used for tracking through the test lattice, although it has limited ability to simulate fringe fields. MAD8 has two ways of modeling a parallel-faced dipole; SBEND and RBEND. SBEND is a sector bend, but this can be given parallel faces by adding appropriate edge angles. RBEND is nominally parallel-faced. Repeating the 10 degree parallel-faced dipole test with 1 meter deflection radius, the MAD8 SBEND element produces the following transfer matrix: 1.0000000 .17364818  0.0000000 1.0000000 and with a 90 degree deflection angle:1.0000000 1.0000000  0.0000000 1.0000000 68 Now with the RBEND command, first 10 degree deflection:  .17364818 .13877788e − 16 1.0000000  1.0000000 1.0000000 1.0000000  and 90 degrees: 0.0000000 1.0000000 It is worthy of mention that the (x|a) element is not identically 0 for the 10 degree RBEND in contrast to the SBEND. The results of tracking in MAD8 to determine the dynamic aperture are negative. MAD8 does not use Enge coefficients to model the fringe fields; rather, it uses the fringe integral [28]. Calculating the fringe integral for the fringe fields used in the above tests for COSY INFINITY and ZGOUBI, we find the y-b tracking picture seen in Fig. 4.8. Figure 4.8: MAD8 tracking pictures with fringe fields. As a final note, we attempted to run similar tests for MADX 5.02.00 (64-bit linux) [27], but were unable to change the edge angles in this version of the code. In conclusion, the vertical aperture is limited exclusively by the fringe field effects, which introduce non- vertical field components. Ignoring these effects (as in MAD), or computing them only 69 approximately (as in ZGOUBI, or by forcing COSY INFINITY to run at unusually low order), yields inaccurate results. 4.2 Summary In this chapter, we compare the behavior of various tracking codes with respect to their accuracy in handling edge angles, especially in the presence of fringe fields. In addition, the relation of the various map elements as a function of aperture was explored. As far as the prediction of dynamic aperture as fringe fields are taken into account, it is clear that the higher order ability of COSY Infinity to handle these cases is unmatched in accuracy with respect to the other codes. Results show that it is important to accurately represent the fringe fields of magnets together with higher order symplectic tracking to viably represent a physical lattice. 70 Chapter 5 Comparative Tracking Recently a set of benchmarks for tracking software based on high precision analytical es- timates for spin tracking has been proposed [35]. Our group at Michigan State University is actively engaged in implementing these benchmark tests in a variety of codes. Here we present some results of this effort. 5.1 Model for Benchmarking A straightforward test, easy to implement in any code, is off-reference motion through a dipole. The motion of a cosine-like ray through such a constant field is diagrammed in Fig. 5.1. 71 Figure 5.1: (x,x’) coordinate of an off-axis cosine-like trajectory in a dipole. The solid circle is the reference orbit of radius R, and the location of the test particle is represented by the dashed circle and angle θ. The dashed circle is offset from the reference orbit by the quantity x0. The 16 radial lines in the circle of the reference orbit denote the locations where the coordinates of the test particles will be plotted. A tracking code will calculate the position of the displaced particle by the particle optical coordinates x and a which can be derived from geometric considerations: (cid:113) x = x0 cos θ + x(cid:48) = x0 sin θ − x2 0 sin θ cos θ 0 sin2 θ (R2 − x2 (cid:113) 0 sin2 θ − R R2 − x2 Here x(cid:48) is the derivative with respect to polar angle, which can be used to obtain the respective particle optical coordinate for slope or relative momentum. The main appeal of this test lies in the fact that regardless of the launching position x0, the orbit will return to exactly this point after one revolution through the lattice. Any deviation from this shows a limitation of the code, either due to numerics, or more importantly, through approximations 72 θ x X0 in the dynamics. 5.2 Test 1: Runge-Kutta The first test is a comparison of a Runge-Kutta integrator written in-house at MSU with COSY INFINITY at 19th order. Figure 5.2 shows the results of Runge-Kutta tracking (10000 turns @ ≈ 300 seconds CPU time). This is in excellent agreement with the analytic solution all the way out to .9R. Figure 5.2: (x,a) coordinate tracking by Runge-Kutta of a cosine-like trajectory through a dipole at 7-70% radius (courtesy Eremey Valetov, MSU). 5.3 Test 2: COSY INFINITY Figure 5.3 is a plot of the same system tracked to 19th order by COSY INFINITY for 10000 turns. What is apparent is the loss of fidelity at approximately .9R. This is to be expected, as COSY INFINITY tracking is based on a Taylor expansion about the reference orbit, and there will always be some radius beyond which a Taylor expansion will lose its accuracy. On 73 the other hand, a Runge-Kutta integrator suffers no such limitation. We therefore expect the accuracy of an RK integrator to be independent of R. Nonetheless, even at .9R, agreement of COSY INFINITY with the analytic solution is on the order of 10−5. Additionally, COSY INFINITY is an order of magnitude faster; this plot took only ≈ 30s. Figure 5.3: (x,a) coordinate tracking by COSY INFINITY of a cosine-like trajectory in a dipole from 7-70% radius (courtesy Eremey Valetov, MSU). 5.4 Test 3: Zgoubi Similar benchmarking was performed using ZGOUBI (Fig. 5.4). This tracking picture is substantially different from that of COSY INFINITY; it is worthwhile to investigate why. Fig. 5.5 shows the result of eight single turns through the dipole. The initial locations of the particles are along the horizontal. The successive turns form branches which spiral away from the horizontal. The increasing curvature of the particle configuration causes the phase space orbits to fill out and become solid lines. Another interesting observation is the behavior of the first 5 particles from the center. These particles do not share in the spiraling motion; 74 thus we suspect a software-related discrepancy between the handling of particles at x < 5cm and those for which x > 5cm. Figure 5.4: (x,a) ZGOUBI tracking of a cosine-like trajectory through a dipole at 7-70% radius. Figure 5.5: ZGOUBI locations of particles after first 8 turns. Each branch is a successive turn. To see whether particles indeed return to their starting position, Fig. 5.6 shows the same tracking test for MAD8. This tracking uses thick elements, which yield qualitatively similar orbits to ZGOUBI. The nonlinear features apparent in the ZGOUBI tracking picture of Fig. 5.4 are not apparent in the MAD8 tracking. 75 Figure 5.6: (x,a) MAD8 tracking of a cosine-like trajectory through a dipole at 7-70% radius. 5.5 Test 4: MADX Finally, we repeat the same benchmark test with MADX, which does not use thick elements for tracking. Fig. 5.7 shows the tracking picture. The asymmetrical nature of the orbits is lacking, which is due to approximations in the Hamiltonian MADX makes to allow the use of thin lens kicks for thick elements. Again the orbits are filled, as we found in ZGOUBI. If we study the first few turns and plot the location of the particles after each turn, we observe a similar “spiraling” as present in ZGOUBI. Figure 5.7: (x,a) MADX tracking of a cosine-like trajectory through a dipole at 7-70% radius. 76 Figure 5.8: MADX locations of particles after first 4 turns. Each branch is a successive turn. 5.6 Summary Nonphysical “stress tests” are often good ways to evaluate the accuracy of a particular code. While it is not likely that dipoles of radius 1 m will be used in any physical sense, such a simulation can point out weaknesses in an algorithm which may not be readily apparent when modelling an actual machine. Clarification of such limits can help reveal subtle errors. For example, some tracking errors which are obvious to observe when the model is pushed beyond the paraxial limit, may not be noticeable in a realistic simulation until after many thousands of orbits. In the case of ZGOUBI, we the difference treatment of orbits at the radial extents of 5 cm boundary would not have been apparent under normal circumstances. 77 Chapter 6 Electrostatic Deflectors for LISE++ The program LISE++ [68] has been developed to calculate the transmission and yields of fragments produced and collected in a spectrometer. It was desired to add second order electrostatic deflectors (cylindrical and spherical types) to the code. Upon review of the literature, the required matrix elements were found to have been computed by Hermann Wollnik in 1965 [69]. Relevant matrix elements were coded, but after initial tests, some questions about their correctness arose. The strategy for verifying the correctness of these calculations was to derive them in several independent ways, and compare the results with various independent codes. What follows is a brief summary of that analysis which demonstrates that the implementation in LISE++ is correct. 78 6.1 Analysis Rather than follow the more traditional approach via equations of motion (see appendices B and C), we follow Wollnik and Glaser[58] in using the index of refraction for an electrostatic field φ in curvilinear coordinates u,w: (cid:104) ee0 (cid:105)1/2 . φ[(1 + u)2 + u(cid:48)2] G = U where e0 is the fundamental charge, and e is the number of such charges. In the original paper, only the results are tabulated, so we present the explicit derivation here. First expand G via a Taylor series in u and u(cid:48): G = G00 + G10u + G01u(cid:48) + G20u2 + G11uu(cid:48) + G02u(cid:48)2 + G30u3 + G21u2u(cid:48) + G12uu(cid:48)2 + G03u(cid:48)3 where u(cid:48) = du/dw, w being the path length of the reference orbit. Upon substitution of this form of G into the Euler-Lagrange equations of motion (cid:19) (cid:18) ∂G ∂u(cid:48) − ∂G ∂u = 0, d dw we find that ∂G (cid:19) (cid:18) ∂G ∂u(cid:48) = G01 + G11u + 2G02u(cid:48) + G21u2 + 2G12uu(cid:48) + 3G03u(cid:48)2, ∂u(cid:48) (cid:19) ∂G ∂u − ∂G ∂u = G11u(cid:48) + 2G02u(cid:48)(cid:48) + 2G21uu(cid:48) + 2G12u(cid:48)2 + 2G12uu(cid:48)(cid:48) + 6G03u(cid:48)u(cid:48)(cid:48), = G10 + 2G20u + G11u(cid:48) + 3G30u2 + 2G21uu(cid:48) + G12u(cid:48)2, = 2G02u(cid:48)(cid:48) − 2G20u − G10 + 2G12uu(cid:48)(cid:48) + G12u(cid:48)2 − 3G30u2 + 6G03u(cid:48)u(cid:48)(cid:48) = 0. d dw (cid:18) ∂G ∂u(cid:48) d dw 79 The coefficients which need to be evaluated are G10, G20, G02, G12, G30andG03. Setting C = ee0/U and S = (1 + u)2 + u(cid:48)2, we abbreviate G = C1/2φ1/2S1/2. Again, φ is the potential. For example, in a spherical deflector the Taylor expansion for a (cid:20)1 2 (cid:21) 1/u potential is φ(u) = r0E0 − u + u2 − u3 To evaluate the coefficients of G, we need to precalculate a few values: r0E0 , 2 φ(0) = φ(cid:48)(0) = −r0E0, φ(cid:48)(cid:48)(0) = 2r0E0, φ(cid:48)(cid:48)(cid:48)(0) = −6r0E0, s(0) = 1, ∂s(0) = 2, ∂u ∂2s(0) ∂u2 = 2, ∂s(0) ∂u(cid:48) = 0, ∂2s(0) ∂u(cid:48)2 = 2, ee0E0r0 2U0 = 1. This gives (cid:104) φ−1/2φ(cid:48)S1/2 + φ1/2S−1/2S(cid:48)(cid:105) = 0, G10 = ∂G ∂u = C1/2 2 80 ∂2G 1 ∂u2 2 C1/2 G20 = = + 4 φ−1/2φ(cid:48)S−1/2S(cid:48) − 1 1 2 2 C1/2 [−1 2 [−1 2 φ−3/2φ(cid:48)2S1/2 + φ−1/2φ(cid:48)(cid:48)S1/2 + φ−1/2φ(cid:48)S−1/2S(cid:48) 1 2 φ1/2S−3/2S(cid:48)2 + φ1/2S−1/2S(cid:48)(cid:48)] φ−3/2φ(cid:48)2S1/2 + φ−1/2φ(cid:48)(cid:48)S1/2 + φ−1/2φ(cid:48)S−1/2S(cid:48) = − 1 2 = −1 2 4 φ1/2S−3/2S(cid:48)2 + φ1/2S−1/2S(cid:48)(cid:48)] , G02 = G12 = , 1 2 1 2 ∂2G ∂(cid:48)2u(cid:48) = ∂3G ∂u∂u(cid:48)2 = −1, ∂ ∂u [−1 2 φ−3/2φ(cid:48)2S1/2 + φ−1/2φ(cid:48)(cid:48)S1/2 + φ−1/2φ(cid:48)S−1/2S(cid:48) − 1 2 φ1/2S−3/2S(cid:48)2 G03 = C1/2 6 + φ1/2S−1/2S(cid:48)(cid:48)] C1/2 = φ−3/2φ(cid:48)2S−1/2S(cid:48) φ−5/2φ(cid:48)3S1/2 − φ−3/2φ(cid:48)φ(cid:48)(cid:48)S1/2 − 1 4 φ−1/2φ(cid:48)(cid:48)S−1/2S(cid:48) 3 [ 4 φ−3/2φ(cid:48)φ(cid:48)(cid:48)S1/2 + φ−1/2φ(cid:48)(cid:48)(cid:48)S1/2 + φ−3/2φ(cid:48)2S−1/2S(cid:48) + φ−1/2φ(cid:48)(cid:48)S−1/2S(cid:48) − 1 2 φ−1/2φ(cid:48)S−3/2S(cid:48)2 + φ1/2S−5/2S(cid:48)3 − φ1/2S−3/2S(cid:48)S(cid:48)(cid:48) 3 4 φ−1/2φ(cid:48)S−1/2S(cid:48)(cid:48) − 1 φ1/2S−3/2S(cid:48)S(cid:48)(cid:48) + φ1/2S−1/2S(cid:48)(cid:48)(cid:48)]. 2 1 2 6 − 1 2 − 1 2 − 1 4 1 2 + φ−1/2φ(cid:48)S−3/2S(cid:48)2 + φ−1/2φ(cid:48)S−1/2S(cid:48)(cid:48) This is a very large expression, so before we verify it, let us determine exactly which of these coefficients are relevant to the (x|xa) element, for example. Going back to the equation of 81 motion 2G02u(cid:48)(cid:48) − 2G20u − G10 + 2G12uu(cid:48)(cid:48) + G12u(cid:48)2 − 3G30u2 + 6G03u(cid:48)u(cid:48)(cid:48) = 0, 2G02u(cid:48)(cid:48) − 2G20u − G10 = −2G12uu(cid:48)(cid:48) − G12u(cid:48)2 + 3G30u2 − 6G03u(cid:48)u(cid:48)(cid:48), We use variation of parameters to solve this. The homogeneous equation is 2G02u(cid:48)(cid:48) − 2G20u − G10 = 0. The idea is to substitute u = A cos(kθ) + B sin(kθ), u(cid:48) = −kA sin(kθ) + kB cos(kθ), u(cid:48)(cid:48) = −k2A cos(kθ) − k2B sin(kθ) = −k2u, into the expression or equivalently into This gives 2¯guu(cid:48)(cid:48) + ¯gu(cid:48)2 + 3F30u2, [3F30 − 2k2¯g]u2 + ¯gu(cid:48)2. [3F30 − 2k2¯g][A2 cos2(kθ) + 2AB cos(kθ) sin(kθ) + B2 sin2(kθ)]+ ¯gk2[A2 sin2(kθ) − 2AB cos(kθ) sin(kθ) + B2 cos2(kθ)]. 82 We are only interested in the coefficients of the AB terms, 6[F30 − k2¯g] cos(kθ) sin(kθ) The variation of parameters integral looks like this: (cid:20) (cid:20) 6(F30 − k2¯g) k 6(F30 − k2¯g) (cid:90) (cid:90) (cid:90) (cid:90) (cid:18)w (cid:21) = = (cid:21) sin(kθ) cos2(kθ) sin(kθ)dθ − cos(kθ) sin2(kθ) cos(kθ)dθ sin(w) k2 6(F30 − k2¯g) (cid:20) k2 6(F30 − k2¯g) 3k2 sin2(w) cos(w)dw − sin(w) (cid:19) (cid:18)w 0 cos2(w) sin(w)dw − cos(w) (cid:19)(cid:21) (cid:104)− sin(w)(cos3(w) − 1) − cos(w) sin3(w) (cid:105) − cos(w) sin3(w) cos3(w) 3 3 0 = = 2(F30 − k2¯g) k2 sin(kθ)(1 − cos(kθ)), which matches the (x|xa) from Wollnik exactly. The other terms are derived in a similar fashion, and have been implemented into LISE++ consistently. 6.2 Comparison of Analytical Results, COSY Infinity 6.2.1 First Order Electrostatic Cylindrical Deflector Next we compare the electrostatic dipole models of COSY Infinity 9.1 and LISE++ 10.0.8. We want to compare both cylindrical and spherical cases up to order 2. The first test is a proton with kinetic energy of 1 eV (so that it is nonrelativistic), radius 1 m, through a deflection angle of θ = 20 degrees. 83 Analytically, first order transfer matrix is (x|x) (a|x)  =  cos −√ √ 2θ √ 2θ 2 sin (x|a) (a|a) 1√ 2 √ sin √ 2θ cos 2θ  =  .880607466 .335060062 −.670120124 .880607466  Also, we want to compare the results of Wollnik [69]. u(w) = u2c1 + α2s1/k1 + F10(1 − c1)/k2 1 The notation requires explanation. F10 = 0 for particles that are on-energy, k1 = √ 2 for a cylinder, c1 = cos k1w and s1 = sin k1w. Not surprisingly, we find match to first order. Setting up this configuration in COSY Infinity we get (x|x) (a|x)  =  .880607466 .335060062 −.670120123 .880607466  (x|a) (a|a) and with LISE++ we get (x|x) (a|x)  = (x|a) (a|a)  .88062 .33506 −.67003 .88062  Next we want to look at the second order transfer maps. Using Wollnik’s approach, we start with u(w) = H0 + Huu + H(cid:48) uu(cid:48) + Huuu2 + Huu(cid:48)uu(cid:48) + Hu(cid:48)u(cid:48)u(cid:48)2 Each of the coefficients are very general expressions covering a wide variety of configurations. Since we are only interested in on-energy particles in a cylindrical deflector, they simplify 84 somewhat. We have: H0 = 0 Hu = cos H(cid:48) 1√ u = 2 Huu = −1 6 √ 2w √ 2w √ 7 sin2 sin (cid:104) (cid:105) √ 2w 2w + 1 − cos H0, Hu, and Hu(cid:48) represent the first order coefficients we expect. The second order coefficient Huu is what we need to confirm. According to COSY INFINITY in the test case, this term should be Huu=-.28185099384. This value matches the analytical calculation. LISE++ gives Huu=-2.8187e-4. This is in agreement considering that LISE++ works in mm-mrad. The next term to check is Hu(cid:48)u(cid:48). According to Wollnik this should be (cid:104) Hu(cid:48)u(cid:48) = 1 12 √ 2w − 8(1 − cos √ 2w) 7 sin2 = .051381097 (cid:105) whereas COSY Infinity gives Hu(cid:48)u(cid:48)=.05138109687 and LISE++ gives Hu(cid:48)u(cid:48)=5.1381e-5 which agrees because LISE++ works in mrads. Finally, we calculate Huu(cid:48) which analytically is √ Huu(cid:48) = − 7 3 2 sin √ 2w(1 − cos √ 2w) + 2 ∗ 1√ 2 √ sin 2w = .57677822790 Here we have corrected a factor of 2. COSY Infinity gives Huu(cid:48)=.57677822805 and LISE++ gives Huu(cid:48)= 5.7678e-4, which is a good agreement (The factor of 2 in the second term of the analytic expression is missing in Wollnik). 85 Next we want to look into the u(cid:48) terms. These appear in Wollnik’s work as u(cid:48)(w) = D(cid:48) uu + D(cid:48) u(cid:48)u(cid:48) + D(cid:48) uuu2 + D(cid:48) uu(cid:48)uu(cid:48) + D(cid:48) u(cid:48)u(cid:48)u(cid:48)u(cid:48) It simply remains to evaluate these coefficients √ u = −s∗p = −p sin(pw) = − D(cid:48) 2 sin √ 2w which evaluates to D(cid:48) gives D(cid:48) u=-.670120124. COSY Infinity gives D(cid:48) u=-.67012012325 and LISE++ u=-.67003 and hence all three agree. √ D(cid:48) √ uu = −sin 2w 2 3 (cid:104) √ 2w + 1 8 cos (cid:105) Wollnik evaluates to D(cid:48) D(cid:48) uu =-.3083909478703871 and LISE++ evaluates to D(cid:48) uu =-.8985037333705609 and COSY evaluates to uu =-8.9869e-4. Wollnik and LISE++ agree, but not COSY. D(cid:48) √ u(cid:48)u(cid:48) = − 4 3 2 sin √ 2w(1 − cos √ 2w) This evaluates to D(cid:48) D(cid:48) u(cid:48)u(cid:48) =-.3483946185202147 and LISE++ gives D(cid:48) u(cid:48)u(cid:48) =-.05333822632965834 and COSY evaluates to u(cid:48)u(cid:48) =-5.3344e05. Again, LISE++ agrees with Wollnik, but not COSY. Finally, the cross term. Duu(cid:48) = −4 3 (2 sin2 √ 2w + cos √ 2w − 1) which evaluates to Duu(cid:48)=-.4395579292132237 and LISE++ gives Duu(cid:48)=-4.3961e-04 but 86 COSY gives Duu(cid:48)=-.1098894819385779. 6.2.2 First Order Electrostatic Spherical Deflector Next we want to perform the same comparisons on the spherical electrostatic deflector models of COSY Infinity 9.1 and LISE++ 10.0.8. Analytically, first order transfer matrix is (x|x) (a|x)  =  cos φ sin φ − sin φ cos φ (x|a) (a|a)  =  0.9396926210 0.3420201430 −0.3420201430 0.9396926210  Also, we want to compare the results of Wollnik. u(w) = u2c1 + α2s1/k1 + F10(1 − c1)/k2 1 In the present case, F10 = 0 for particles that are on-energy, k1 = 1 for a spherical deflector, c1 = cos k1w and s1 = sin k1w. Not surprisingly, we find a match to first order. Setting up this configuration in COSY Infinity we get (x|x) (a|x)   0.9396926209 0.3420201433 −0.3420201426 0.9396926209  (x|a) (a|a) and with LISE++ we get(x|x) (a|x)  = (x|a) (a|a)  0.93971 0.34202 −0.34193 0.93971  Turning our attention to the electrostatic spherical deflector, we return to the horizontal 87 displacement u(w): u(w) = H0 + Huu + H(cid:48) uu(cid:48) + Huuu2 + Huu(cid:48)uu(cid:48) + Hu(cid:48)u(cid:48)u(cid:48)2 Again, we are only interested in on-energy particles. H0 = 0 Hu = cos w H(cid:48) u = sin w Huu = − sin2 w The H0, Hu, and Hu(cid:48) are the first order coefficients we expect. The second order coefficient Huu is what we need to confirm. According to COSY INFINITY, this term should be Huu=-0.1169777782, which matches the analytical calculation. LISE++ gives Huu=-0.0001170000. This is in agreement consid- ering that LISE++ works in mm-mrad. The next term to check is Hαα. According to Wollnik this should be Hu(cid:48)u(cid:48) = sin w2 + cos w − 1 = .056670399 while COSY Infinity gives Hu(cid:48)u(cid:48)=0.0566703992 and LISE++ gives Hu(cid:48)u(cid:48)=0.0000056671 which again agrees because LISE works in mrads. Finally, we calculate Huu(cid:48) which analytically is Huu(cid:48) = −2 sin w(1 − cos w) + 2 sin w = .64278761 88 (Again, the factor of 2 in the second term is absent in Wollnik) COSY Infinity gives 0.6427876098 and LISE++ gives 0.0006427900 , which is again a good agreement. Next we want to look into the α terms. These appear in Wollnik’s work as u(cid:48)(w) = D(cid:48) uu + D(cid:48) u(cid:48)u(cid:48)u(cid:48) + D(cid:48) uuu2 + D(cid:48) uu(cid:48)uu(cid:48) + D(cid:48) u(cid:48)u(cid:48)u(cid:48)u(cid:48) It simply remains to evaluate these coefficients D(cid:48) u = −s∗ = − sin(w) = − sin w which evaluates to D(cid:48) gives D(cid:48) u=-0.34193 and hence all three agree. u= -.342020143. COSY Infinity gives D(cid:48) u=-0.3420201426 and LISE++ D(cid:48) uu = − sin w cos w = −0.3213938050 COSY is identically 0 in this term, but LISE++ evaluates to D(cid:48) uu =-0.0003215800. Wollnik and LISE++ agree, but not COSY. D(cid:48) u(cid:48)u(cid:48) = − sin w(1 − cos w) = −.020626338 COSY evaluates to D(cid:48) u(cid:48)u(cid:48)=0.3420201433 and LISE++ gives D(cid:48) u(cid:48)u(cid:48) =-0.0000206320. Again, LISE++ agrees with Wollnik, but not COSY. 89 Finally, the cross term. D(cid:48) uu(cid:48) = −2 sin2 w − cos w + 1 = −0.1736481780 LISE++ gives D(cid:48) uu(cid:48)=0 for this element. At the time, it was not known whether there was a software issue in COSY INFINITY or the uu(cid:48)=-0.0001737000 but COSY is identically D(cid:48) analysis. As an independent check, we checked these results against a Runge-Kutta code, the results of which are presented in the next section. 6.3 Lagrangian Analysis for Comparison With Runge- Kutta To confirm the correctness of the analytical results, we need an independent check. A Runge- Kutta (RK) code was chosen for validation. This analysis follows the general plan of attack by Wollnik. While Wollnik used the Euler-Lagrange equations resulting from treating the electric field as an optical system with a varying index of refraction, for comparison with RK code we work with the Lagrangian directly to provide a somewhat independent method of calculation. L = 1 2 m(r2 ˙θ2 + ˙r2) − φ Instead of using r, we will use r = r0(1 + u): L = 1 2 m(r2 0(1 + u)2 ˙θ2 + r2 0 ˙u2) − φ = mr2 0[(1 + u)2 ˙θ2 + ˙u2] − φ 1 2 90 The Euler-Lagrange equations of motion are given by: (cid:18)∂L (cid:19) − ∂L ∂u = 0. d dt ∂ ˙u To this end, we evaluate and (cid:18) ∂L (cid:19) ∂ ˙u d dt = mr2 0 ¨u ∂L ∂u = mr2 0(1 + u) ˙θ2 − ∂φ ∂u . The equation of motion obtained is mr2 0 ¨u − mr2 0(1 + u) ˙θ2 + ∂φ ∂u = 0. Using the cyclic variable θ: (cid:18) ∂L (cid:19) ∂ ˙θ d dt = d dt (mr2 0(1 + u)2 ˙θ) = 0, or Hence, or mr2 0(1 + u)2 ˙θ = constant ≡ l. m2r4 0(1 + u)4 ˙θ2 = l2 mr2 0(1 + u) ˙θ2 = l2 mr2 0(1 + u)3 , 91 which can be substituted into the equation of motion 6.1 yielding mr2 0 ¨u − l2 mr2 0(1 + u)3 + ∂φ ∂u = 0. Turning our attention to the potential, we have for a spherical potential φ = −k r = − k r0(1 + u) ,with k a constant, and the field is therefore ee0E = − 1 r0 ∂φ ∂u = − k r2 0(1 + u)2 . To follow Wollnik’s nomenclature, we use k = ee0E0r2 0. The equation of motion becomes mr2 0 ¨u − l2 mr2 0(1 + u)3 + ee0E0r0 (1 + u)2 = 0. We will expand the (1 + u) terms to second order in u so first rewrite mr2 0 ¨u − l2 mr2 0 (1 + u)−3 + ee0E0r0(1 + u)−2 = 0. The expansions are as follows: and (1 + u)−3 ≈ 1 − 3u + 6u2 (1 + u)−2 ≈ 1 − 2u + 3u2. 92 Using these expansions, the equation of motion 6.1 takes the form mr2 0 ¨u − l2 mr2 0 (1 − 3u + 6u2) + ee0E0r0(1 − 2u + 3u2) = 0. Combining like terms, (cid:32) mr2 0 ¨u + − 2ee0E0r0 3l2 mr2 0 (cid:33) (cid:32) u + − l2 mr2 0 (cid:33) (cid:32) + ee0E0r0 = − 3ee0E0r0 6l2 mr2 0 (cid:33) u2, and using l2 = m2v2r2 0 = mee0E0r3 0 in this expansion we get or Because it follows that mr2 0 ¨u + ee0E0r0u = 3ee0E0r0u2 ¨u + ee0E0 mr0 u = 3ee0E0 mr0 u2. mv2 r0 = ee0E0, v2 r2 0 = ee0E0 mr0 = ω2 and the equation of motion simplifies to ¨u + ω2u = 3ω2u2. (6.1) This is not the type of equation that can be solved by variation of parameters. The right hand side is a function of u, and not a function of the independent variable. So what we do in 93 this case is an approximation. We assume that u is the general solution. This approximation makes the right hand side a function of the independent variable, and that allows us to use variation of parameters to obtain an approximate solution. Solving the homogeneous equation: ¨u + ω2u = 0 u(t) = A cos ωt + B sin ωt gives with A = u(0) and B = 1 ω ˙u(0) To get the second order solutions, we need to take this homogeneous solution and sub- stitute this into the inhomogeneity term giving 3ω2(A cos ωt + B sin ωt)2. Multiplying this out, we get 3ω2A2 cos2 ωt + 6ω2AB cos ωt sin ωt + 3ω2B2 sin2 ωt. 6.4 Term by Term Comparison of Analysis and RK First, work with the term containing AB. This becomes the integral (cid:90) t sin ωt (cid:20) 1 ω cos ωˆt(6ω2AB cos ωˆt sin ωˆt)dˆt − cos ωt 0 0 94 (cid:90) t (cid:21) sin ωˆt(6ω2AB cos ωˆt sin ωˆt)dˆt , where ˆt is a dummy variable of integration. Evaluating this integral gives: 2AB sin ωt(1 − cos ωt), and the entire solution is (because B = 1/ω) (u(0)|u(0) ˙u(0)) = sin ωt(1 − cos ωt). 2 ω This result does not match Wollnik, because Wollnik is relativistic. It matches RK numerical results, because that is nonrelativistic. Next evaluate the (u(0)|u(0)u(0)) case. This consists of the AA term when we substitute the homogeneous solution into the inhomogeneity. sin ωt cos ωˆt(3ω2AA cos ωˆt cos ωˆt)dˆt − cos ωt 0 0 sin ωˆt(3ω2AA cos ωˆt cos ωˆt)dˆt The solution is (u(0)|u(0)u(0)) = sin2 ωt + 1 − cos ωt, which again matches the nonrelativistic RK. Finally, the BB part of the solution: sin ωt cos ωˆt(3ω2BB sin ωˆt sin ωˆt)dˆt − cos ωt 0 0 sin ωˆt(3ω2BB sin ωˆt sin ωˆt)dˆt (cid:90) t (cid:90) t (cid:20) 1 ω (cid:20) 1 ω (cid:90) t (cid:90) t (cid:21) (cid:21) And the solution for this is (B = 1/ω) (u(0)| ˙u(0) ˙u(0)) = 1 ω2 (cid:104) (cid:105) 2(1 − cos ωt) − sin2 ωt . 95 These match the nonrelativistic RK simulation. Fig 6.1 shows the output from the RK simulation used to verify the results obtained. Figure 6.1: x as a function of deflection angle. This is x as a function of the electrostatic element deflection angle, or x = (x|x)x0 + (x|a)a0 + (x|xx)x2 0 + (x|ax)a0x0 + (x|aa)a2 0. If we set the initial condition a0 = 0 then we have only x = (x|x)x0 + (x|xx)x2 0, and we can check the accuracy of the (x|xx). Fig. 6.2 is a plot of x − (x|x)x0. 96 Figure 6.2: Plot of x − (x|x)x0 = (x|xx)x2 0. Compare this to Fig. 6.3, a plot of the analytically calculated value of (x|xx). Figure 6.3: Analytically calculated value of (x|xx). 97 According to the RK, when the deflection angle is 45 degrees, (x|xx) = .793 and according to the analytic formula, (x|xx) = .79289. Next we want to look at the case where x0 = 0 or x = (x|a)a0 + (x|aa)a2 0, shown in Fig. 6.4. Figure 6.4: x = (x|a)a0 + (x|aa)a2 0. Fig. 6.5 shows the output from the RK integrator, (x|aa) = x − (x|a)a0. 98 Figure 6.5: RK result, (x|aa) = x − (x|a)a0. and we must compare this to the analytic solution: (cid:105) 2(1 − cos ωt) − sin2 ωt (cid:104) (x|aa) = 1 ω2 which is shown in Fig. 6.6 99 Figure 6.6: Analytic solution of (x—aa). Finally, we need to look at the cross term (x|xa). Since this requires both x0 and a0 to be nonzero, we need to subtract from x all other terms: x − (x|x)x0 − (x|a)a0 − (x|xx)x2 0 − (x|aa)a2 0 = (x|xa) ∗ a0 ∗ a0 Fig. 6.7 is the result from the RK: 100 1234561.(cid:180)10(cid:45)82.(cid:180)10(cid:45)83.(cid:180)10(cid:45)84.(cid:180)10(cid:45)8XAAAnalytic Figure 6.7: RK calculation of (x|ax). which is to be compared to the analytic solution sin ωt(1 − cos ωt) 2 ω shown in Fig. 6.8: 101 Figure 6.8: Analytic solution of (x|xa). Now let us consider the (a|∗) terms. Fig. 6.9 is the (a|∗) and shows the output of the RK as a function of deflection angle: 102 123456(cid:45)2.(cid:180)10(cid:45)8(cid:45)1.(cid:180)10(cid:45)81.(cid:180)10(cid:45)82.(cid:180)10(cid:45)8XXAAnalytic Figure 6.9: RK plot of a as a function of deflection angle. As a Taylor expansion this looks like a = (a|x)x0 + (a|a)a0 + (a|xx)x2 0 + (a|ax)a0x0 + (a|aa)a2 0 If we set x0 = 0, then a = (a|a)a0 + (a|aa)a2 0 So we can look at the RK calculated plot of the (a|aa) term, plotted in Fig. 6.10 103 Figure 6.10: RK plot of (a|aa). Calculating the α term is done by evaluating u(cid:48)(1 − u). The idea is to pick the term you are interested in, and recreate it. We are interested in the BB term. I need the BB term of du/dw, and add that to the BB terms of du dw u The BB term of du/dw is simply the derivative of (cid:105) 2(1 − cos ωt) − sin2 ωt (cid:104) (x|aa) = 1 ω2 which is sin ωt(1 − cos ωt) 2 ω 104 and the BB term of (cid:18) d (cid:19) sin ωt sin ωt dt ω ω du dw u = = 1 ω cos ωt sin ωt Subtracting these gives which matches Wollnik. sin ωt ω (2 − 3 cos ωt) Now this all works fine if you are trying to match Wollnik, but it does not look like the RK test. To match the RK test simply take the derivative of the corresponding (x|∗) term. In this case we want (x|aa) = d dt 2 sin ωt ω (1 − cos ωt) which matches the RK as shown by Fig.6.11: Figure 6.11: (a|aa) by taking the derivative of (x|aa). 105 123456(cid:45)0.0002(cid:45)0.00010.00010.0002AAAAnalytic Next look at (a|xx). Fig. 6.12 shows the result of the RK calculation: Figure 6.12: RK calculation of (a|xx). and Fig. 6.13 is the analytic calculation Figure 6.13: Analytic calculation of (a|aa). 106 123456(cid:45)0.00015(cid:45)0.00010(cid:45)0.000050.000050.000100.00015AAAAnalytic Finally, we want to look at the (a|xa) term. First, Fig. 6.14 shows the RK simulation, Figure 6.14: RK calculation of (a|xa). and here is the analytic calculation of (a|xa) in Fig. 6.15 Elements (a|xx), (a|xa), (a, aa) Figure 6.15: Analytic Solution of (a|xa). for spherical electrostatic deflector (a|aa) = 2 sin ωt ω (1 − cos ωt) 107 123456(cid:45)0.0004(cid:45)0.0003(cid:45)0.0002(cid:45)0.00010.00010.0002AXAAnalytic (cid:104) (cid:105) (a|xa) = 2 sin2 ωt − cos2 ωt + cos ωt (a|xx) = ω sin ωt(2 cos ωt + 1) Finally, we need to change coordinates from derivatives with respect to t or ωt to derivatives with respect to arc length. This means α(w) = du dw r0 r = du dw 1 1 + u ≈ u(cid:48)(1 − u) = u(cid:48) = u(cid:48)u To do this for (α|αα) we need all terms with αα in them. The first term (u(cid:48)) will be (u(cid:48)|uu(cid:48)). But then we need to subtract from that the product (u(cid:48)|u(cid:48))(u|u(cid:48)) which is (α|αα) = 2 sin ωt ω (1 − cos ωt) − 1 ω cos ωt sin ωt = sin ωt ω (2 − 3 cos ωt) Fig. 6.16 shows what (α|αα) looks like from the equation: and Fig. 6.17 shows the graph of Figure 6.16: RK (A—AA). the derived formula: Next consider the (α|uu) matrix element. Now we need (α|uu) = (u(cid:48)|uu) − (u(cid:48)|u)(u|u) 108 Figure 6.17: Analytic A—AA. which is simply (α|uu) = ω sin ωt(2 cos ωt + 1) + ω cos ωt sin ωt = ω sin ωt(3 cos ωt + 1) This settled the issue, and the spherical electrostatic deflector code in COSY INFINITY Figure 6.18: RK (A—XX). has since been modified [59] 109 123456(cid:45)3.(cid:180)10(cid:45)8(cid:45)2.(cid:180)10(cid:45)8(cid:45)1.(cid:180)10(cid:45)81.(cid:180)10(cid:45)82.(cid:180)10(cid:45)83.(cid:180)10(cid:45)8AAAAnalytic Figure 6.19: Analytic (A—XX). 6.5 Independent Check Following Berz,Makino,Wan As an independent check of the accuracy of these second order matrix elements, it is instruc- tive to derive them from the traditional equations of motion. We start with Berz, Makino, Wan (BMW) Eq. 3.22 [19]. a(cid:48) = d ds (cid:18)px (cid:19) p0 = (1 + hx) = (1 + hx) + h ps p0 1 + η 1 + η0 γ(v) γ(v0) p0 Ex ps χe0 ZeEx p0v0 p0 ps + h ps p0 . In this equation, m is the rest mass of both the reference particle and the test particle, v0 is the speed of the reference particle, v is the speed of the test particle, Ex is the electric field at the location of the test particle, h is the instantaneous curvature of the path of the reference particle, χe0 = Ze/p0v0 is the electric rigidity of the reference particle, s is the 110 123456(cid:45)2.(cid:180)10(cid:45)8(cid:45)1.(cid:180)10(cid:45)81.(cid:180)10(cid:45)82.(cid:180)10(cid:45)8AXXAnalytic arc length along the reference path, p0 is a normalization constant, and px and ps are the components of the test particle momentum in the directions normal and tangential to the path of the reference particle, respectively. Because p0 is a constant, we can move it out of the derivative and cancel leaving dpx ds = (1 + hx) = (1 + hx) = (1 + hx) = (1 + hx) γ(v) γ(v0) γ(v) γ(v0) ZeEx v0 p0 ps + hps ZeEx p0 v0 γ(v)mvs + hps p0 vs + hps ZeEx γ(v0)mv0 ZeEx + hps = ZeExt(cid:48) + hps. vs In the last step we applied from page 62 of BMW vs = dL dt = dL ds ds dt = (1 + hx) t(cid:48) . The hps term can be cast in a different form by virtue of hps = hγ(v)mvs = γ(v)m(1 + hx)h/t(cid:48). The γ(v) does not cancel from this term, therefore this is a relativistic equation of motion requiring a relativistic Lagrangian. First let us consider the nonrelativistic form of the BMW equation of motion, i.e., = ZeExt(cid:48) + m(1 + hx)h/t(cid:48). dpx ds 111 The equation of motion follows from the nonrelativistic Lagrangian, L = 1 2 m[x(cid:48)2 + (1 + hx)2]/t(cid:48)2 − V. To form the Euler-Lagrange equations of motion when path length is the independent vari- able, we need to consider the action integral. From Lanczos [70], we see that upon replacing time t with a parameter s, the action integral becomes (cid:90) t2 (cid:90) s2 Ldt = Lt(cid:48)ds t1 s1 which leads to the path-dependent form of the Euler-Lagrange equations of motion: (cid:18) ∂Lt(cid:48) (cid:19) ∂x(cid:48) d ds − ∂Lt(cid:48) ∂x = 0. There are two ways to eliminate t(cid:48) from the equations of motion. The first way (as in BMW), is to use conservation of energy and solve for t(cid:48). This gives E − ZeV = T = = 1 2 m(1 + hx)2/t(cid:48)2 mx(cid:48)2/t(cid:48)2 + m[x(cid:48)2 + (1 + hx)2]/t(cid:48)2, 1 2 1 2 or t(cid:48) = (cid:115) m 2 [x(cid:48) + (1 + hx)2] E − ZeV . Substituting this expression for t(cid:48) into Lt(cid:48) and simplifying gives (cid:113) m 2 [x(cid:48)2 + (1 + hx)2][E − 2ZeV ] √ E − ZeV . Lt(cid:48) = 112 Setting E = 0 allows us to simplify this expression to (cid:113)−2mZeV [x(cid:48)2 + (1 + hx)2]. Lt(cid:48) = The second way (as in Wollnik [69]) is to use Jacobi’s method, and eliminate t(cid:48) from the Lagrangian itself, forming the “modified Lagrangian” ¯L as follows: ¯L = Lt(cid:48) − ptt(cid:48) = (L − pt)t(cid:48) = 2T t(cid:48) = 2(E − ZeV )t(cid:48) = −2ZeV t(cid:48) The new action integral becomes (1 + hx)2 + x(cid:48)2ds (cid:90) s2 s1 (cid:113) (cid:90) s2 (cid:114) m (cid:90) s2 (cid:90) s2 (cid:113)−2mZeV [x(cid:48)2 + (1 + hx)2] 2T t(cid:48)ds = (cid:113) √ 2mT s1 (1 + hx)2 + x(cid:48)2ds 2T s1 s1 2T ¯Lt(cid:48)ds = = = This explicitly demonstrates that these two methods are equivalent in that they lead to identical equations of motion. 6.6 BMW to Second Order We start with BMW Eq. 3.22 [?]. a(cid:48) = d ds d ds (cid:18) px (cid:18) px p0 p0 1 + η 1 + η0 γ(v) γ(v0) p0 Ex ps χe0 ZeEx p0v0 + h ps p0 p0 ps + h ps p0 (cid:19) (cid:19) = (1 + hx) = (1 + hx) 113 At this point we make the explicit assumption that p0 is a constant. Then we can move it out of the derivative and cancel leaving dpx ds = (1 + hx) = (1 + hx) = (1 + hx) = (1 + hx) γ(v) γ(v0) γ(v) γ(v0) ZeEx v0 p0 ps + hps ZeEx p0 v0 γ(v)mvs + hps ZeEx γ(v0)mv0 ZeEx vs p0 vs + hps + hps. Finally, we assume nonrelativistic physics so that we can write mvs = ps: dpx ds = (1 + hx) mZeEx ps + hps. We need an expression for ps. We use conservation of energy U U − ZeV = K = p2 2m = x + p2 p2 s 2m . Where V is the potential. Solving for ps gives (cid:113) ps = 2m(U − ZeV ) − p2 x. So this is the equation of motion. Dividing the answer by p0 will give a(cid:48) as in (cid:113) + h 2m(U − ZeV ) − p2 x. a(cid:48) = dpx dps = (cid:112)2m(U − ZeV ) − p2 (1 + hx)mZeEx x 114 We want to solve this equation explictly for an electrostatic spherical deflector. That means we need Ex and V . A particle with circular motion at radius R = 1/h will have k = − 2U0R Ze where U0 is the total energy of such a particle. This will not be the same as U , the total energy of a test particle. We have Ex = − k (R + x)2 and V = k R + x . Now we need the first order approximation to the right hand side. The field and potential to first order are V = −Ex0R (cid:18) (cid:19) x2 R2 (cid:19) x , → 1 − x R + (cid:18) 1 − 2 R Ex = −Ex0 Ex = −Ex0, ∂Ex ∂x = 2Ex0 R , V = −Ex0R, ∂V ∂x = Ex0. Another assumption that we make is that the energy of the particle is the same as the energy of a particle on the reference orbit. Since that orbit would be circular, the energy is U = − 1 2 ZeEx0R. How does this modify the expression for the tangential component of momentum ps? Using ps = 2mZeEx0 this modifies the equation of motion to (cid:104) dpx ds = (1 + hx)mZeEx (cid:16) R 2 − x (cid:17) − p2 x (cid:105)1/2 2mZeEx0 (cid:20) (cid:18) R 2 (cid:21)1/2 (cid:19) − p2 x − x (cid:20) + h 2mZeEx0 115 , (6.2) (cid:19) − x (cid:18) R 2 − p2 x (cid:21)1/2 . (6.3) And to first order, Eq. 6.2 becomes ps = [2mZeEx0R]1/2 − (cid:20)mZeEx0 (cid:21)1/2 R x. (6.4) Now we evaluate the first derivative of the right hand side of Eq 6.3 in several steps, after making the substitution (1 + hx)mZeE(cid:48) x hmZeEx A1/2 + x + 2AmZeE(cid:48) 2AhmZeEx 2AA1/2 2AA1/2 2AhmZeEx + 2AmZeE(cid:48) 2AA1/2 −2m2Z2e2E2 x0 + 4m2Z2e2E2 2 (cid:21) x A = (cid:20) A1/2 − x (cid:19) 2mZeEx0 (cid:18) R − p2 − (1 + hx)mZeExA(cid:48) AhA(cid:48) − mZeExA(cid:48) 2AA1/2 x − mZeExA(cid:48) + AhA(cid:48) x0 − 2m2Z2e2E2 2(mZeEx0R)(mZeEx0R)1/2 2AA1/2 2AA1/2 + −[mZeEx0R]1/2 = = = = = −4m2Z2e2E2 x0 2(mZeEx0R)(mZeEx0R)1/2 Hence we have shown that dpx ds = − p0 R2 x or d ds 116 (6.5) ; hA(cid:48) 2A1/2 + x0 x0 − 2m2Z2e2E2 = −−2mU R2 = − p0 R2 . R2 (cid:18) px (cid:19) p0 = − x R2 . What are the Field, Potential, and transverse momentum to second order? (cid:18) (cid:19) (cid:18) V = −Ex0R Ex = −Ex0 1 − x R 1 − 2x R + x2 R2 − x3 (cid:19) R3 (cid:20) + (cid:18)1 2 3x2 R2 − x R (cid:21)1/2 (cid:19) − p2 x + x2 R2 px = 2mZeEx0R From here we need to evaluate the second derivatives. For convenience, I repeat the first derivatives here: hmZeEx A1/2 + (1 + hx)mZeE(cid:48) x A1/2 − (1 + hx)mZeExA(cid:48) 2AA1/2 hA(cid:48) 2A1/2 + 117 First term: Second term: Third term: Fourth term: hmZeE(cid:48) A1/2 x − hmZeExA(cid:48) 2AA1/2 hmZeE(cid:48) A1/2 x mZeE(cid:48)(cid:48) A1/2 x xA(cid:48) − mZeE(cid:48) 2AA1/2 + −hmZeExA(cid:48) 2AA1/2 xA(cid:48) − mZeE(cid:48) 2AA1/2 − mZeExA(cid:48)(cid:48) 2AA1/2 3mZeExA(cid:48)2 4AAA1/2 + hA(cid:48)(cid:48) 2A1/2 − hA(cid:48)2 4AA1/2 Now we want the denominators to match up. 4hmZeE(cid:48) 4AAA1/2 xAA − 2hmZeExAA(cid:48) 4AAA1/2 −2hmZeExAA(cid:48) 4AAA1/2 4hmZeE(cid:48) 4AAA1/2 xAA xAA(cid:48) 4AAA1/2 + xAA 4mZeE(cid:48)(cid:48) 4AAA1/2 xAA(cid:48) − 2mZeE(cid:48) 4AAA1/2 2hAAA(cid:48)(cid:48) 4AAA1/2 − 2mZeE(cid:48) − 2mZeExAA(cid:48)(cid:48) − hAA(cid:48)2 4AAA1/2 4AAA1/2 3mZeExA(cid:48)2 4AAA1/2 + Now that all the fractions have common denominators, I need to evaluate everything. 4hmZe(E(cid:48) x)(A)(A) − 2hmZe(Ex)(A)(A(cid:48)) (cid:18)2Ex0 (cid:19) R =4hmZe (mZeEx0R)2 − 2hmZe(−Ex0)(mZeEx0R)(−2mZeEx0) =4m3Z3e3E3 x0 118 4hmZeE(cid:48) =4hmZe(E(cid:48) xAA + 4mZeE(cid:48)(cid:48) x)(A)(A) + 4mZe(E(cid:48)(cid:48) xAA − 2mZeE(cid:48) xAA(cid:48) (cid:18) (cid:19) x)(A)(A) − 2mZe(E(cid:48) =4hmZe (mZeEx0R)2 + 4mZe x)(A)(A(cid:48)) (mZeEx0R)2 −6Ex0 R2 (cid:18)2Ex0 (cid:19) (cid:19) (cid:18)2Ex0 R (mZeEx0R)(−2mZeEx0) − 2mZe R = − 8m3Z3e3E3 x0 − 2hmZe(Ex)(A)(A(cid:48)) − 2mZe(E(cid:48) = − 2hmZe(−Ex0)(mZeEx0R)(−2mZeEx0) − 2mZe x)(A)(A(cid:48)) − 2mZe(Ex)(A)(A(cid:48)(cid:48)) + 3mZe(Ex)(A(cid:48))2 (mZeEx0R)(−2mZeEx0) (cid:18)2Ex0 (cid:19) R (cid:18)4mZeEx0 (cid:19) R + 3mZe(−Ex0)(−2mZeEx0)2 − 2mZe(−Ex0)(−2mZeEx0R) = − 24m3Z3e3E3 x0 2hAAA(cid:48)(cid:48) − hAA(cid:48)2 =2h(A)2(A(cid:48)(cid:48)) − h(A)(A(cid:48))2 (cid:18)4mZeEx0 (cid:19) − h(mZeEx0R)(−2mZeEx0)2 R =2h(mZeEx0R)2 =4m3Z3e3E3 x0 Combining all the fractions x0 − 8m3Z3e3E3 x0 − 24m3Z3e3E3 x0 + 4m3Z3e3E3 x0 4m3Z3e3E3 = − 24m3Z3e3E3 4A2A1/2 x0 = −6(mZeEx0R)1/2 R3 = −6p0 R3 Hence we have shown that dpx dx = − p0 R2 x − 3p0 R3 x2 or (cid:18) px (cid:19) p0 d ds = − x R2 − 3 x2 R3 119 This is the second order equation of motion for an electrostatic spherical deflector, which agrees with Wollnik Eq. 6.1. 6.7 Summary This chapter encapsulates the main idea behind this dissertation, namely that a seemingly innocuous discrepancy between two codes can give rise to an in-depth investigation with numerous benefits. First, any existing problems with the code are ferreted out and resolved, benefitting the user base and increasing the confidence with which the code can be used. Second, the investigation itself can be very instructive and informative for those performing the analysis. No less than four methods of analyzing the same physics (and finally demon- strating their equivalence) were undertaken. None of this would have occurred had only one code been used. This careful analysis uncovered implementation errors in COSY Infinity which were subsequently corrected [59]. 120 Chapter 7 Oak Ridge Isomer Spectrometer and Separator (ORISS) The Oak Ridge Isotope/Isomer Spectrometer and Separator (ORISS) is an electrostatic multiply reflecting time-of-flight (MRTOF) mass separator that was built around 2010 by the University Radioactive Ion Beam Consortium (UNIRIB) [41]. The device was never fully commissioned due to the Oak Ridge group losing funding, and ended up at Michigan State University for use at the Facility for Rare Isotopes and Beams (FRIB). MRTOF devices are typically used at very low particle intensities due to strong Coulomb repulsion at the particle turning points. It is the intent of this research to determine whether this device can be adapted to function under the high intensity modes which will be prevalent at FRIB. 7.1 Background The earliest mass separator was Thomson’s Parabola Spectrograph [61], shown in Fig. 7.1. “Spectrograph” means that the presence of multiple ion species are recorded simultaneously by a single detector, as in a photograph. This is in constrast to a mass spectrometer which is designed to measure a single ion species. The ion beam to be analyzed was collimated and passed through parallel electric and magnetic fields, which caused the ions to be deflected in mutually perpendicular directions. An ion would be deflected by the electric field by 121 Figure 7.1: Thomson’s device (1913). an amount y ∝ q perpendicular to the first by x ∝ q mv2 . This same ion will be deflected by the magnetic field in a direction mv . Eliminating v between these two equations yields y ∝ m q x2. In other words, particles impacting the output graph along a parabola will have the same mass to charge ratio, which can then be calculated. Aston [62] pointed out that any spread in the velocity of the ions would cause the width of the parabolas to increase, hindering an accuate measurement. His proposal (Fig. 7.2) was allow the ions to pass through an electric field P1 − P2 first. A spread in velocities would cause a divergence of the beam as they exited slit D. This divergence would be compensated by the weak focusing as the beam passes through the magnetic field represented by the circular magnetic poles. As a result, all ions of a given q/m ratio would be focused onto a single spot on the graph G. In the 1940’s and 1950’s, mass spectrometers based on magnetic sectors came into common use (Fig. 7.3. The idea is very simple, particles with differing masses will be detected at differing radii. However, there are some limitations to this technique [51]. A magnetic mass spectrograph requires a physically large device to accommodate a wide mass range. A mass spectrometer, 122 Figure 7.2: Aston’s device (1920s). which is a device that measures only a single q/m ratio at a time, requires a large field range for a wide mass range. Also, in order to accomodate a rapid sampling rate, a low hysteresis is required. MRTOF devices do not suffer from these limitations. Furthermore, an added benefit of MRTOF is elegance of concept; the system can be compact and is relatively easy to fabricate. The basic principle of MRTOF operation is simple. Consider two ions with identical kinetic energies, K1 = K2 = K, such as can be prepared by launching identically charged particles through a common voltage drop. Then for ions travelling in free space, K1 =K2 = K, 1 2 m1v2 1 = 1 2 m2v2 2, m1v2 1 = m2v2 2. Clearly if m1 (cid:54)= m2, then v1 (cid:54)= v2, and the distance between the ions (∆d) will increase in time. This is the meaning behind the designation “Time of Flight” separator. If ∆d is the 123 Figure 7.3: A magnetic mass spectrograph requires a physically large device to accommodate a wide mass range. A mass spectrometer requires a large field range for a wide mass range, and low hysteresis for a sample rate. Figure 7.4: Time of flight detector with long path length. separation between the ions, we have √ ∆d = d2 − d1 = t(v2 − v1) = t 2K (cid:18) 1√ m2 (cid:19) , − 1√ m1 from which we see that we can increase the magnitude of ∆d by increasing the value of t, the length of time the ions fly. One way to increase the flight time of the ions is by increasing the path length (Fig. 7.4). An alternative approach is to place electrostatic mirrors at each end of the flight path, causing the ions to be reflected multiple times between the 124 Figure 7.5: Time of Flight detector with multiple reflections. mirrors (Fig. 7.5). This is the technique used in “Multiply-Reflected Time of Flight” (MRTOF) separators. The figure of merit describing the efficiency of a mass separator is “mass deviation”. Mass deviation can be written as dm m , but in the literature of mass spectroscopy it is traditional to take the reciprocal of this quantity called “mass resolution” or m dm. From D = vt, and getting v from the kinetic energy, (cid:18)2K (cid:19) m = t2. D2 Taking differentials, which leads to the simple relation (cid:18)2K (cid:19) D2 2t dt, dm = m dm = t 2dt , again showing that increasing time of flight increases mass resolution. Fig. 7.6 is a photograph of ORISS resting on a workbench at Louisiana State University during construction. (All photographs and drawings are courtesy of Louisiana State University). Fig. 7.7 shows a mechanical drawing of the device within its vacuum jacket, and Fig. 7.8 is a representation showing the two mirror electrode regions at each end of the central drift. The device is about 130 cm from end to end. During injection, the potential of the left mirror is momentarily reduced, allowing the particles to enter the system. After injection, the original potentials are restored, and the end electrodes act as mirrors. Fig. 7.9 shows a detailed depiction of 125 Figure 7.6: ORISS during construction (Drawing courtesy LSU). Figure 7.7: Oriss within vacuum jacket (Drawing courtesy LSU). Figure 7.8: Oriss showing pertinent components. 126 Figure 7.9: Electrostatic mirror geometry showing the numbering of each electrode. Figure 7.10: Mirror electrode geometry. the left electrostatic mirror. The right mirror geometry is symmetric. There are 8 ring electrodes, and a conical electrostatic lens. All electrodes within each mirror, the conical lens, and the central drift region can each be independently biased. The ring electrodes (see Fig. 7.10) have a large aperture and can be approximately modeled as biased rings of diameter 7 cm. The conical lenses (see Fig. 7.11) have a smaller aperture diameter of 5 cm. This is the smallest aperture in the system, and therefore limits radial excursion of the particles. These conical lenses can be removed, and replaced with the large-aperture variety [65]. Nevertheless, we acted on the assumption that our beam would always be restricted to 127 Figure 7.11: Conical lens geometry. 2.5 cm radial excursion. The principle of time of flight mass spectrometry depends critically on the initial ion kinetic energies being identical. This is because we want any difference in velocity to be due to a difference in mass alone. Since there will always be some initial distribution of energies ∆E, the device should ideally be operated to ensure that ions with slightly different initial kinetic energies have isochronous orbits. Since an ion with a higher kinetic energy will have a higher velocity, this requires that the higher energy particle travel a slightly longer path into the reflecting mirror. This leads to a condition of energy isochrony. A derivation, based on Ref. [42], follows. The energy of a nonrelativistic ion moving through a free drift region is purely kinetic: with constant velocity E = K = mv2. 1 2 (cid:114) 2 m √ E. v = The time required for an ion with some reference energy E to travel the drift length L is 128 therefore (cid:114) m 2E . tE = L/v = L For off-energy ions the time required is tE+∆E = L (cid:114) m 2(E + ∆E) . The difference between these two times is, to first order, ∆tDrift ≡ tE+∆E − tE = L m 2(E + ∆E) − L (cid:114) (cid:114) m 2E (cid:114)m 2 ≈ −L 1 2 1 E3/2 ∆E. Next, we perform an analysis for ions traveling through a constant axial electric field Ez = −a = constant, which crudely represents the electrostatic mirror. For a particle of energy E, the turning point in the field will be at z = E/ea, and we have: E = v = tE = tE+∆E = mv2 + eaz, √ 1 2 (cid:114) 2 (cid:90) dz (cid:114)m m v(z) 2 0 E − eaz, (cid:114) m (cid:90) E+∆E = 2 ea ∆tφ = tE+∆E − tE = (cid:90) E/ea 0 dz√ E − eaz √ dz E + ∆E − eaz (cid:112)2m(E + ∆E) √ = , 2mE ea (cid:112)2m(E + ∆E) (cid:114) m √ ea ≈ 1 ea 2mE ea 2 = − ∆E√ E . ea 129 For energy isochrony, we require condition the cancellation of these time differences, that is or which simplifies to (cid:114) m 2 −L ∆tφ + ∆tDrift = 0, (cid:114) m 2 1 ea ∆E√ E = 0. 1 2 1 E3/2 ∆E + E ea = L 2 . 7.2 Literature Review Historically, MRTOF for mass spectrometry has been developed for use at very low intensity and analyzed as single particle devices. A good resource for understanding the historical development of time-of-flight mass analyzers is the review by Hermann Wollnik[36]. In that review, the earliest mention of such a device is 1948 [37], constructed by A.E. Cameron and D.F. Eggers. This device used 3 electrodes, at voltages +300 V, -500 V, 0 V, and had a pulse width of 5 microseconds. It was not a multiple-reflection type, but rather used a single 3.47 m tube. Particles were introduced into one end, and the time of flight to reach the other side was recorded. The first TOF analyzer which took energy isochrony explicitly into account was created in 1955 by Wiley and McLaren [40]. This device was designed to exploit a new type of ion gun. Wiley lists three advantages for a TOF mass spectrometer and one main disadvantage. Speed is one advantage. This is critical for isotopes which may have a short half-life. A second 130 advantage is that the entire source spectrum can be analyzed, even if the source conditions vary appreciably. The third advantage is that operation relies on biased electrodes, rather than a high mechanical precision. The devices can therefore be relatively cheap and light. On the other hand, the main limitation of a TOF mass spectrometer is limited resolution arising from an ever-present initial spatial and kinetic energy distribution. Following this Wiley and Mclaren device, no improvements were made in TOF mass spectroscopy until 1973, when Mamyrin et al. sent ions into a single electrostatic mirror [38], allowing for an effective doubling of the length of a compact device. Mamyrin describes the advantage of his device as “focusing of the time-of-flight of the ion packets on the basis of energy”. The electrostatic mirror was a flat electric field created by a grid, a backplane, and potential divider rings. In his review of the device, Wollnik points out that the grid is not innocuous and acts as many small lenses. Also the grid has finite transparency. Hence the next development in the progression of time of flight mass spectroscopy was extended to be a grid-free device. This relieves the distortion caused by the grid, but it produces a fringe field outside of the mirror region and complicates device tuning. At this point Wollnik references his own paper with the statement:“The most accurate method to obtain an optimized inhomogeneous field distribution Ez(z) seems to have been shown in the late 1980s”. The reference is to a paper “A Transversally and Longitudinally Fo- cusing Time-of-Flight Mass Spectrometer” [39]. There Wollnik describes optical constraints for focusing and energy isochrony. The next step in the evolution is the reuse of the same channel multiple times for high resolution in a compact device. That is the multiple in “multiply” reflected TOF analyzers. There are several different geometries to achieve this. The one we are interested in utilizes 131 coaxial mirrors. A similar device was sent to Comet 67P/Churyumov-Gerasimenko on the Rosetta mission [52]. This exploited the compact and light system for a demanding space application. ORISS is a grid-free, coaxial multiply reflecting time of flight mass separator and spectrometer [53]. 7.3 Overview of Operational Tuning A simple first model of ORISS is as a superposition of coaxial rings of charge. We can therefore use the well-known on-axis potential of a thin ring of charge Q and radius R at axial position z=0 [60] φ(r = 0, z) = 1 4π0 Q√ R2 + z2 from which we can calculate the off-axis values (to second order) of Ez and Er. Fig. 7.12 is a picture of this potential. The on-axis field Ez = −∂φ/∂z is Ez(r = 0, z) = 1 z 4π0 [R2 + z2]3/2 which is shown in Fig. 7.13. To first order, the value of Ez off-axis is the same as the value of Ez on-axis. To get the value of Ez(r, z) for points off the axis to second order, we apply the Laplace equation for an axisymmetric potential and obtain (see Appendix E): Ez(r, z) ≈ −φ(cid:48) + φ(cid:48)(cid:48)(cid:48)r2. 1 4 (7.1) 132 Figure 7.12: On-axis potential of ring of charge. Here the derivatives are evaluated with r = 0 and the primes denote differentiation with respect to z. Evaluating needed derivatives: z , φ(cid:48) = − Q 4π0 φ(cid:48)(cid:48)(cid:48) = − Q 4π0 (cid:20) [R2 + z2]3/2 − 9z [R2 + z2]5/2 (cid:21) . + 15z3 [R2 + z2]7/2 and inserting into Eq. 7.1 we obtain Ez to second order: Ez = −φ(cid:48) + 1 4 φ(cid:48)(cid:48)(cid:48)r2 z = Q 4π0 [R2 + z2]3/2 (cid:20) + Q 4π0 1 4 9z [R2 + z2]5/2 − 15z3 [R2 + z2]7/2 (cid:21) r2. The second order r2 term is superimposed onto the existing on-axis field. In Fig 7.14 the axial field is plotted for a few values of r0. Curvature of the on-axis potential also gives to 133 01-2-1012V/V0z/R Figure 7.13: On axis field of a ring of charge. leading order a linear radial field ∝ r. We see that Er = −∂φ ∂r ≈ 1 2 φ(cid:48)(cid:48)r = −1 2 1 [R2 + z2]3/2 − 3z2 [R2 + z2]5/2 (cid:18) (cid:19) r (7.2) Fig. 7.15 is a plot of the linear field Er slightly off the reference orbit. Note that the field is strongly focusing near the plane of the ring at z = 0. We can take advantage of this focusing region if we set the turning point of the orbit close enough to the ring. In Fig. 7.16, we plot a particle approaching the mirror from the right, with an initial transverse displacement of .1 mm. First the particle is deflected away from the axis, then, as it reaches the focusing region, it is strongly focused toward the axis. There is a minimum energy below which a particle will not be focused. This is shown in Fig. 7.17, where we show a particle approaching the mirror with an initial radius of 3 mm, but the energy is too low to reach the focusing region. This particle is therefore deflected from the axis and will be radially lost. The energy must be high enough to reach the focusing region, but not so high that it is overfocused and exits the system. Between these two extrema, we find the envelope of a stable orbit between 134 0-2-1012V/V0z/R Figure 7.14: Axial field off axis to second order for various values or r. two rings of charge (see Fig. 7.18). For a given energy within the stable range, we find a limit (see Fig. 7.19) on the transverse acceptance (here at about 6 mm for the example in Fig. 7.18). What is the nature of this limit? If we look at a plot that simultaneously shows the potential in the radial and transverse coordinates (see Fig. 7.20), we find that the radial field Er(r) goes as r, whereas the axial field Ez(r) goes as r2. Hence the axial field of the mirror will always eventually dominate the transverse field at large enough off-axis transverse distances. The unstable orbit depicted in Fig. 7.21 illustrates this eventual dominance of the axial field Ez over the radial field Er. That’s as far as we can take the simple ring focusing model. To make further progress we must introduce separate focusing elements to allow simultaneous tuning for transverse and energy-isochronous focusing. From the literature, we can set any of the middle electrodes (see Fig. 7.9) to a negative value to form an einzel lens. It can be shown [44] that any einzel lens is focusing. 135 -101-2-1012V/V0z/Rr=.0r=.25r=.5 Figure 7.15: Off-axis electric field Er of a ring of charge. Figure 7.16: Trajectory of a particle revealing both focusing and defocusing of a ring of charge. z=0 m is the center of the device, z=-.6 m is the location of the mirror. 136 -0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05-5-4-3-2-1012345V/V0z/Rr=.1 Figure 7.17: Radial orbit of off-axis particle. z = 0 m is the center of the device, z=-.6 m is the location of the mirror. Figure 7.18: Envelope of a stable orbit between two rings of charge. The mirrors are located at ±.06 m, and the particle starts at z = 0 m,r = .003 m. 137 Figure 7.19: Transverse phase space stability. Figure 7.20: Potential of a ring of charge (x is radial, y is axial). 138 Figure 7.21: Example of unstable orbit that would be lost. 139 7.4 Particle-In-Cell (PIC) Simulations The software used in these simulations is part of the Warp code[54]. Warp is a well-known, well-maintained and long-standing code that uses the particle-in-cell (PIC) method for simu- lating large numbers of interacting charged particles. The Lorentz equation of motion is used to advance macroparticles with a physical q/m ratio. Each macroparticle can represent many real particles, and the ratio of real particles to macroparticles is called the macroparticle’s “weight”. The electric fields, both those applied by external conductors and those which are beam-generated, are calculated on a spatial mesh. Conducting structures are placed onto the mesh. The macroparticles, which are represented as a density on the same mesh, are advanced in time by the leap-frog method which is second order in time. Use of this low order advance is desirable since only one solution of Poisson’s equation is required per time step. However, this requires use of smaller timesteps. What makes working with scientific simulations challenging is the ease with which com- puters can give the researcher a misleading sense that they are doing actual physics, when in reality the simulations my be creating inaccurate pictures unless run with proper convergence. A critical component (one which is often overlooked or simply bypassed) is benchmarking. For example, before an investigator can model thousands of charged particles moving under the influence of complex electrostatic fields, how well can the software model situations in which the analytic solution is available? We have already discussed this with regards to the potentials. Knowing that Warp can simulate known potentials (and knowing what situations to avoid), increases the level of confidence with which more advanced potentials can be simulated. Before attempting to model the full ORISS geometry with Warp, a few benchmark tests 140 were designed and implemented. The first test was to place a simple charged particle within the potential of a charged sphere of 1 cm radius and compare numerical and analytic poten- tials. This simple model was chosen because the result is easy to verify. It was discovered that the simulation grid had to be extraordinarily large to match the analytic solution. This is because the presence of the grid boundary altered the shape of the potential at distances of only a few centimeters. Once the potential matched the analytic solution, the next step was to test the propagation of a single particle (not macroparticle) within this potential. A simple energy conservation test was used whereby the particle was allowed to propagate some distance through the potential, and then its direction of motion was reversed. The particle should retrace its steps to the initial position and indeed this was the case, to a high degree of precision. The next simulation test was the ring of charge. Although a single ring of charge had already been found to be unsuitable for ORISS modeling, it still formed a beneficial test with at least a partially analytic solution. Both the shape of the potential (to second order) and the trajectory of particles within the potential were tested with satisfactory results. Once the initial benchmark tests had been passed, it was necessary to load the actual configuration of ORISS conductors into Warp. This is quite easy. Warp has a wide variety of tools to load conducting geometric shapes from which the ORISS conductors can be constructed. A model that matched the geometric specifications of ORISS was coded and set up to allow for arbitrary biases. Referring to Table 7.1, we see the various electrode settings which will be referred to in several places in this document. The table also uses Fig. 7.9 to identify the various electrodes. Once the model had been set up, it was necessary to determine the required size of the mesh. ORISS is a long, thin, cylindrically symmetric device. Because the potential outside of the device is not relevant to the simulation, a 141 Electrode Settings E1 +0 A +0 B +0 C +0 D +0 E +0 F +0 G +0 H +0 I +0 J +0 K +0 L M +0 +0 N +0 O P +0 Q -5000 R +0 E2 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 -5000 +0 E3 +0 +0 +0 +0 +0 -72000 +0 -1000 -2000 -3000 -4000 -6000 -7000 -8000 -9000 -10000 -5000 +0 E4 +50 +50 +50 +50 -9000 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 E8 E7 +50 +50 +50 +50 E5 +50 +50 +50 +50 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 E6 +100 +50 +200 +50 +400 +50 +50 +800 +3000 +3000 +4500 +8000 +24000 +36000 +10000 +10000 +10000 +10000 +10000 +10000 +10000 +10000 +10000 +10000 +6000 +2000 +2000 +7000 +14425 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 Table 7.1: Electrode settings in volts. Each column refers to an electrode specified in Fig. 7.9, each row is a bias configuration referred to in the text. cylinder of 2 m length and radius 10 cm was used. This mesh is large enough to enclose the conductor, but not too large to cause an undue amount of processing effort for solving for the potential. Warp has the ability to use refined meshes in areas requiring greater detail, so we used a finer mesh in the area near the axis. The possibility was considered of leveraging the axial mirror symmetry of ORISS and use a mesh that was reflected at the midpoint. However, simulations have been carried out in (r, z) using mesh refinement giving high speed and efficiency and we wanted to preserve the full geometry to allow eventual injection and extraction gating. So the reflection symmetry was not pursued. The on-axis potential using the smaller mesh, as well as the potential gen- erated by the maximal gridsize (7m x 7m) were compared and shown to be satisfactory in the 142 region of interest near the axis. Two meshes were used; a “coarse” mesh for rough simula- tions and a “fine” mesh for more detailed runs. The coarse mesh was 1.5625 mm/gridpoint in the transverse direction (10cm), and 1.953 mm/gridpoint in the longitudinal direction (2 m). The fine mesh was 24.4 µ m/gridpoint in the transverse direction, and .4883 mm/gridpoint longitudinal. It is also important to benchmark the modeling of space charge. A simple test in this regard is to take two singly charged 238U+ ions, and place them at rest a distance 1 mm from each other and calculate at which rate they will separate. A 238U+ ion has a mass of 3.952e-25 kg. The initial acceleration of each ion is determined from Coulomb’s law: a = q2 m4π0r2 = 583.8m/s2 This initial acceleration will drop as the particles become substantially separated, but it sets the scale. Setting up this scenario in Warp reveals that the simulation is accurate as shown in Fig. 7.22. This test also helps visualize the scale at which space charge becomes relevant. 143 Figure 7.22: Plot showing agreement between 238U+ ion repulsion in Warp and analytic solution via Coulomb’s law. 7.5 Space Charge Effects FRIB will be a CW SRF linac with high 400kW power. Rare isotopes will be separated by a multistage fragment separator and be transported to a variety of systems including an RF cooler buncher in which up to 108 ions will be trapped and cooled. These bunches will be the input for ORISS. Is it even possible to run this many particles through ORISS? This was a question that had to be answered before even considering whether ORISS could be used for delicate measurements such as mass spectrometry. The results of tracking two bunches containing ions under the influence of space charge through ORISS in a typical potential configuration is shown in Fig. 7.23. The bunches were initially combined; the first bunch consisted of 238U+ ions and the second contained 237Np+ ions. The lower curve shows the separation ∆z between the centroids of the two bunches as a function of time, while the upper curve shows the evolution of the axial RMS extent σz of the 238U+ bunch as a function of time. Note that the RMS extent σz grows more 144 Figure 7.23: Plot of bunch centroid separation (lower curve) and RMSz (upper curve) vs. time. rapidly than the separation of the bunches due to time of flight mass separation. Therefore space charge is inhibiting this operating point from behaving as a separator. Because space charge gives an energy spread to the particles, we require ∆z to be independent of ∆E, the energy deviation, at the end of a full orbit. To operate with space charge, we need to find a potential configuration that is both energy isochronous and point-to-point focusing. To explore this possibility, simulations are started with low voltage electrode settings as shown in row A of Table 7.1. The shape of the on-axis potential is reproduced in Fig. 7.24. Given these settings, we sought a kinetic energy that was energy-isochronous. This was found by testing a range of central reference energies and then launching a bunch with ±5% spread 145 Figure 7.24: Typical potential with no additional transverse focusing. of energies. The objective was to find a reference energy that resulted in a longitudinal focal point at the midpoint of the system (z = 0). For these electrode settings, such a reference energy was found to be 37.63 eV (5523 m/s). As can be seen from Fig. 7.25, all particles in a range about this central energy return to a focus at the midpoint. The second criterion to be met was transverse parallel to point focusing. A parallel set of rays must be focused to a point at the center of the device. As it happens these electrode settings do not focus sufficiently, as shown in Fig. 7.26. An additional focusing influence is required. From the earlier discussion regarding off-axis fields (Eq. 7.2), we have Er = φ(cid:48)(cid:48)r 1 2 so a focusing lens can be created by adding a very slight φ(cid:48)(cid:48) < 0 near the turning point. This can be accomplished by setting electrode E3=-4V. The resulting potential is shown in Fig. 7.27. This produces negative curvature as we move away from the mirror. Starting 146 Figure 7.25: Optimal energy focusing. from the maximum potential, as we move away in the positive direction, the slope of the potential becomes more negative giving φ(cid:48)(cid:48)(z) < 0 and a resulting Er that is focusing. As we continue to move away from the peak, we reach a point of inflection where the slope begins to increase, hence in this region Er > 0 and is therefore defocusing. Repeating the analysis which produced Fig. 7.26 with this additional focusing influence gives the parallel-to-point focusing we require as shown in Fig. 7.28. The successful separation of 238U+ and 237Np+ after only 5 reflections is illustrated in Fig. 7.29. Figs. 7.30 and 7.31 show two screenshots from a movie of the same configuration. Next we increase the initial particle density by increasing the number of particles while keeping the initial bunch volume constant. By increasing the voltage of the mirrors, which increases the velocity of the bunch and reduces the time spent at the turning points, we can achieve separation before the space charge has a chance to fully manifest itself. This technique will work if the different species are present in similar quantities. The separation of the species has the effect of reducing the density of the bunch, which in turn reduces the effect of the space charge. Fig. 7.32 shows an example 147 Figure 7.26: Energy isochronous, with no additional transverse focusing. where the mirror voltage is set so low that the space charge effects dominate mass separation. This negative effect can be reversed by increasing the mirror voltage (Figs. 7.33, 7.34, 7.35). However, as the space charge continues to increase by increasing the number of particles, the intensity of the space charge grows to dominate this effect. The next plots (Fig. 7.38 shows what happens as we increase the number of particles to 100k and model space charge consistently. Unfortunately, for this electrode voltage configuration, increasing the mirror voltage is no longer sufficient to compensate. We need to adapt the focusing potentials to compensate the force of the space charge repulsion. The next step was to try 105 particles, and see if we could keep the particles contained. The first working electrode settings are shown in row E of Table 7.1. The first objective is to find a particle energy that is energy isochronous. This corresponds to a velocity of 42424 meters/second. However, this arrangement is not parallel to point and transverse focusing, in fact, it is not even stable. So it was necessary to reduce this velocity by 20%. This algorithm allowed us to find electrode settings that worked all the way up to 106 148 Figure 7.27: Same potential with slight transverse focusing. particles, but at that time we ran into another problem. The time it took to simulate > 106 particles began to increase dramatically. It was not known at the time if this was due to simple scaling, or if some other influence was increasing the run time. So at this point we began to run macroparticles. Using 106 particles with a macroparticle weight of 100 gets us to the goal of 108 particles. Using 106 macroparticles with a weight of 100 corresponds to 108 physical particles. With electrode settings given in row F of Table 7.1, this configuration is energy isochronous at vref = 113715. As before, in order to make this arrangement focusing and stable we need to reduce this velocity by 20%. Success! This is a very coarse operating point, which only contains the beam, and should be treated as hypothetical, especially because of the large electrode voltages used (which may exceed breakdown). Nonetheless, it serves to show what factors need to be taken into consideration when working with such intense space charge. Next, in Sec. 7.6 we construct more reasonable electrode settings which can handle this high amount of space charge intensity more elegantly. 149 Figure 7.28: Point to parallel. Figure 7.29: Mass separation. Solid curve is distance between centroids of ion species, dashed curve is RMS radius. 150 Figure 7.30: Initial configuration. Figure 7.31: Separation after 1 reflection. 151 Figure 7.32: Bunch centroid separation compared with bunch RMS radius, 10000 particles, electrode settings from row A of table 7.1. Due to space charge, expansion of the bunch dominates any bunch separation. Figure 7.33: Doubling the mirror voltage to 200V as in row B of 7.1 increases the bunch separation. 152 Figure 7.34: Mirror voltage of 400V, row C of table 7.1. Figure 7.35: Mirror voltage 800V, row D of table 7.1. 153 Figure 7.36: Electrode settings row D, now 20000 particles. Figure 7.37: Electrode settings row D, now 40000 particles. 154 Figure 7.38: Electrode settings row D, now 80000 particles. 155 7.6 Tuning for Mass Separation at High Space-Charge Intensity Being able to contain the beam is only part of the challenge. If the optics are not optimal, nonlinearities can dominate the separation process due to the large path length needed for separation. To explore this, while keeping the parameter space manageable, we settle on a simple form of potential using only a single electrode for a mirror and an additional electrode for a focusing lens. Such a potential is shown in Fig. 7.39 for voltages given in Table 7.1. There are several “knobs” we can adjust. Figure 7.39: Minimal potential configuration. • n - Number of physical particles per bunch • E - bunch Energy • E3 - voltage of focusing lens • E8 - mirror voltage 156 All other electrode voltages are grounded for simplicity. Even with this simplification, this is a large parameter space. The first objective is to narrow it down, and make it more manageable. Initially, we choose a very low particle density (n = 103) so space charge is not strong. Due to the requirements of energy isochrony, we must have energy focusing at the midpoint of the device. This in turn is dependent upon the energy of the beam relative to the electrode strength of the mirror. To understand this, refer to Fig. 7.40. The curve is the shape of the mirror potential in the proximity of the turning point, with the potential φ(z) satisfying φ(cid:48)(cid:48)(z) < 0 at the z of the reference particle turning point. Indicated are two separate possible turning points, for different energy levels of the reference particle. The higher energy reference particle will travel a longer path length. The slope of the potential is exaggerated in the figure, to illustrate that for a gentler slope, a more energetic particle will travel deeper into the potential as compared to a particle with the reference energy. Conversely, a lower energy particle will travel less deeply into the potential. In this way, choosing the energy of the reference particle has the effect of “tuning” the energy isochronous condition. The energy focusing of the system is illustrated by Fig. 7.41. The horizontal axis is time (or simulation timesteps) and the vertical axis is z, the longitudinal position of the bunch. The bunch starts at the center of the system (z = 0 cm) and moves toward the reference mirror location (z = .6 cm). At the turning point (around t=13000 timesteps), the particle begins to move back toward z = 0. A distribution of energies is launched. The particles of various energies focus around t=15000 timesteps. Due to aberrations, the focusing is not linear, but approximate. The ideal would be for the longitudinal beam waist to be at z = 0, the system center. With the voltage of the focusing lens E3 set to 0 V (row G of table7.1, the transverse trajectories of the particles are as depicted in Fig. 7.42, showing inadequate transverse fo- 157 Figure 7.40: Energy isochronous condition as a function of beam energy. cusing. Clearly we would like to increase the transverse focusing in order to obtain parallel to point. But what will happen to the longitudinal focusing as we apply the transverse lensing? We examine what happens in Figs. 7.41 and 7.42. In Fig. 7.43 we see the effect Figure 7.41: Energy Focusing : Focusing strength F= 0kV (n=10e3, E=3.4kV, V=10kV, row G). of changing only the voltage of the focusing electrode to 1 kV. The longitudinal focal point has moved somewhat closer to z = 0, which is fortuitous, and increased transverse focusing is evident, albeit still inadequate. As we continue to increase the voltage E3 on the focusing 158 Figure 7.42: Transverse trajectory : Focusing strength E3= 0 kV (n=10e3, E=3.4kV, V=10kV, row G). lens, the longitudinal focusing moves towards z = 0 (see Figs. 7.44, 7.45), a desired result. At E3=4 kV we find that we have passed the optimum. We cannot (for this electrode con- figuration) reach optimal longitudinal and transverse focusing simultaneously, so we must be satisfied with a compromise. Using this electrode voltage configuration, we now increase the density of physical particles by two orders of magnitude (to n = 106 particles) and observe the effects of space charge on the optics. The results are shown in Fig. 7.47. Fortunately, energy focusing remains fairly stable to preserve the isochronous condition. Note the large change in transverse scale. As before, we want to increase the voltage E3 of the focusing electrode to try for point to parallel focusing. Fig. 7.48 shows how the focusing lens has the desired countering effect on the radial spread of the beam. As we reach a focusing electrode setting of 10 kV, the optics are such that the beam will stay radially confined over many turns, even with a high degree of space charge. Stable electrode settings for up to 108 particles have been found, but do they work to separate particles? Figs. 7.53, 7.54, 7.55 and 7.56 present a few screenshots from the author’s first movie which successfully showed separation in a case with intense space charge (106 particles). This result demonstrates that if the optics are tuned in such a way that energy isochrony is maintained, and first 159 Figure 7.43: Transverse and Longitudinal focusing E3=1 kV, row H. order optics with space-charge are kept from expanding beyond the limits of the system, then separation can occur even in the presence of intense space-charge. However, it must be emphasized that this process is sensitive to the optics. Figs. 7.57, 7.58, 7.59 show a similar run with non-optimal optical settings. The quality of the movie is slightly improved; more tracer particles are used. 160 Figure 7.44: Transverse and Longitudinal focusing E3=2 kV, row I. 161 Figure 7.45: Transverse and Longitudinal focusing E3=3 kV, row J. 162 Figure 7.46: Transverse and Longitudinal focusing E3=4 kV, row K. Figure 7.47: Focusing strength E3= 0 kV, with space charge, 106 particles, row G. 163 Figure 7.48: Focusing strength E3= 1 kV, with space charge, 106 particles, row H. Figure 7.49: Focusing strength E3= 4 kV, with space charge, 106 particles, row K. Figure 7.50: Focusing strength E3= 8 kV, with space charge, 106 particles, row N. 164 Figure 7.51: Focusing strength E3= 10 kV, with space charge, 106 particles, row P. Figure 7.52: Tracking picture of 106 particles over 20 orbits, row P. 165 Figure 7.53: Initial configuration, 106 particles at z = 0, 50%U 238, 50%N 237. 166 Figure 7.54: Particle bunch reaches the turning point, note the spread due to optics and space charge. 167 Figure 7.55: Bunch returns to center, focusing element constrains the expansion due to space charge. 168 Figure 7.56: Successful separation after several reflections, with intense space charge. 169 Figure 7.57: Particle bunch initial condition. Note the larger number of tracer partices. Figure 7.58: Bunch reaches the turning point, but the optics has distorted the shape of the bunch. 170 Figure 7.59: In an untuned system with space charge, the particles do not separate. 171 7.7 Simulation of Realistic Distributions Up to this point, all simulations have used initial conditions that were somewhat arbitrarily formulated since the primary objective was understanding the optics of the ORISS device. For simulations to be complete, we must more accurately represent the type of distribution injected into ORISS. Fig. 7.60 shows a schematic of the FRIB Re-Accelerator (ReA) including the Beam Cooler Buncher (BCB) which cools and bunches the ions injected into ORISS. The idea is that ORISS will reside between the BCB and the EBIT (Electron Beam Ion Figure 7.60: Schematic of FRIB ReA (Reaccelerator) showing the Beam Cooler Buncher. Trap). For a sense of perspective, a diagram of the ReA platform is shown in Fig. 7.61. Realistic simulations of the BCB ion accumulation have been carried out by Dr. Ryan Ringle [55]. This data can be applied to better represent initial ions to be injected into ORISS. An example simulation accumulated the coordinates of 105 K+ ions with a mean energy of 4.887 keV. Note there was a high energy tail on the distribution (see Fig. 7.62). For simulation purposes, this was converted to a Gaussian by cutting the distribution at the 172 Figure 7.61: Detailed diagram of ReA platform. The BCB (not shown) is directly below EBIT. maximum, and reflecting the lower half of the distribution about the cut. The axial bunch length is .0268 m which translates into an RMS bunch length of 172 ns. Fig. 7.63 shows the transverse distribution. The transverse RMS emittance of this distribution is rms =.85 π-mm-mrad. What effects might we expect from bunches with nonzero temperature? From the nonrelativistic, electrostatic paraxial equation in cylindrical coordinates [44], (cid:20)φ(cid:48)(cid:48) (cid:21) R(cid:48)(cid:48) = − R + 16 2 rms R3 + K R 4φ Here R is the envelope radius, and K is the dimensionless perveance. Because the emittance term ∼ 1/R3 relative to 1/R for the perveance term, we expect the modest thermal effects to have a weaker impact on beam defocusing than space-charge. We test this hypothesis by advancing this distribution through a typical parallel to point electrode configuration. Fig. 7.64 shows the result of a single reflection of the BCB configura- tion space distribution with the emittance-related velocity spreads artificially set to 0. Also, 173 Figure 7.62: ECB longitudinal kinetic energy distribution. we take n = 1000 particles so that the density is low enough that there is insignificant space charge defocusing. There is no focusing electrode and the mirror electrode E8 is 14425 V (row R of table 7.1). The focusing is provided by the mirror electrode alone. Fig. 7.65 repeats the test, but this time using the BCB derived Gaussian-equivalent velocity distribution. The transverse spread of the beam is apparent. Next in Fig. 7.65 we increase the density of the particles by two orders of magnitude (to n = 105), which increases space charge sufficiently to have a pronounced effect. However, the shape of the beam remains essentially unchanged. We can conclude that a nonzero emittance associated more realistic BCB-injected distribu- tions will not impact the operation of ORISS to a significant degree. Finally, in Fig. 7.67, we advance through 20 reflections, and find that the ions are still contained within the aperture of the system, with modest losses. 174 Figure 7.63: ECB transverse distribution. Figure 7.64: BCB equivalent distribution with no emittance. 175 Figure 7.65: Same as Fig. 7.64, but including emittance. Figure 7.66: Same as Fig. 7.65, including space charge. 176 Figure 7.67: with space charge (100k particles), with velocity dispersion, 20 reflections. 177 7.8 ISOLDE MR-TOF Simulations In this section we simulate the ISOLDE ISOLTRAP MRTOF using Warp to test the relia- bility of our modeling. ISOLDE (Isotope Separator On Line DEtector) is a facility located at CERN on the PS Booster [56]. Part of this facility is the ISOLTRAP tandem Penning trap mass spectrometer to which a multi-reflection time-of-flight mass separator (ISOLTRAP MRTOF) was added in 2010. For comparison purposes, the length of the ISOLTRAP MRTOF (IM) is only 80 cm long, compared with ORISS’ 130 cm. IM only has 5 aperture electrodes compared with ORISS’ 8. A third difference is the technique by which injection and extraction is effected. In ORISS, the mirror potentials are temporarily reduced to low values to permit the passage of ion bunches. IM injects bunches with a high energy, and then slows them by temporarily increasing the potential of the midsection drift region. Since we are not modeling injection or extraction, this difference is not a concern. A typical potential implemented in IM is depicted in Fig. 7.68 [48]. We can mimic this potential in ORISS with the electrode settings shown in row G of table 7.1. which yield the potential shown in Fig. 7.69. According to [50], there are three principal uses for an MRTOF. First, it can be used for the analysis of a beam of unknown composition. The beam enters the MRTOF, and a time spectrum is produced by which the various components of the beam can be identified. Second, the beam can be used as a separator to purify beams. After the contaminant species have been separated, they can be removed by a kicker in the extraction beamline. This can produce pure samples required for Penning trap measurements, for example. A third use for the MRTOF is for mass spectrometry of nuclides that have a short half-life. All of these use cases are valuable at FRIB as well. 178 Figure 7.68: Typical potential of ISOLTRAP MRTOF. An experiment which we apply to test our Warp model of an MRTOF based on ORISS is the first measurement by IM of 53Ca+ by comparing the TOF to a couple of known substances, namely 39K+ and 53Cr+. A plot showing the result of the IM experiment is shown in Fig. 7.70 [49]. According to Ref. [49], the mass excess (actual mass - mass number) of 53Ca is -29387.8(43.3) keV/c2 where the quantity in parentheses is the sta- tistical uncertainty. So actual mass = mass excess + 53u = [-29387.8+(53)931494095.4 ]eV/c2=4.936915767e10 eV/c2 = 8.80098015e − 26 kg = mass of 53Ca. The mass of 53Cr from [49] 8.79112959e-26 kg, this is lighter which is consistent. Fig. 7.71 presents the results from our simulation, using the ORISS MRTOF configuration with the IM bias shown in Fig. 7.69 1000 particles and no space charge. The successful separation seen here is based on approximately 20 reflections (10 orbits). These results are consistent with the IM results. First, the time of separation needs to be only 1 ns for a Bradbury-Nielsen gate to effect separation [66]. Second, the IM experiment was run with extremely low particle counts, justifying running such a small number of particles. 179 Figure 7.69: IM potential implemented in ORISS using the electrode settings in row G of table 7.1. Figure 7.70: Time of flight spectrum of A = 53 nuclides delivered from ISOLDE (53Cr+,53Ca+, and the reference ion 39K+. 180 Figure 7.71: TOF spectrum of 53Cr and 53Ca, 5000 ions. 181 7.9 Summary of ORISS investigations Time of Flight mass spectrometry is a well-established technique of beam analysis, but it is usually used under conditions of low current. Some empirical and analytic studies including space charge were undertaken in the early 2000s ([63], [64]), which treated the MRTOF as a potential well with potential walls linear in the axial direction. This analysis considered only the longitudinal evolution of the bunch size, and utilized Simion, which has only an approximate model for space charge [67]. We were fortunate to have received early on precise design drawings of ORISS which al- lowed us to create a highly accurate model. The settings of the electrodes were not available, and we proceeded to study the system from first principles. Realizing that ORISS was similar to a system of coaxial rings of charge, this simple and analytically tractable model was explored. The first and second order transverse and longitundinal single-particle optics of a symmetric two-mirror system using only two coaxial biased rings of charge were studied. It is well known that this model is unstable, and soon it became apparent that more sophisicated optics were required. Because of the large number of parameters, the investigation was limited to systems with reflective symmetry about the midpoint. Each side was comprised of a single mirror, and a single focusing lens. From the literature, the biases of the focusing lenses were chosen to have the opposite polarity as the mirror electrodes. The criteria were twofold; system was to be energy isochronous and parallel-to-point focusing. Operating points were initially found for systems with low particle counts. To ensure that our model was sound, a simulation based on the measurement of 53Ca+ by the MRTOF at CERN/ISOLDE was conducted successfully. 182 Then, higher particle counts were used to increase the influence of space charge on the dynamics. Utilizing proportionately higher electrode biases was helpful here; by increasing the longitudinal velocity of the particles, the transverse spread of the beam due to space charge was suppressed. An operating point capable of containing a beam up to 108 particles was found empir- ically. This was considered only as a proof of principle, and work began immediately to duplicate these results with lower electrode biases. A heuristic development cycle evolved, comprised of increasing the particle counts, increasing the electrode biases, and maintaining energy isochrony and parallel to point focusing. It was discovered that unless these criteria were adhered to, the nonlinearities induced by space charge overwhelmed the longitudinal separation of particles by mass alone. The transverse and longitudinal optics must be tuned to keep these nonlinearities in check in order that mass separation can be achieved. 183 Chapter 8 Conclusions The instability of the COSY J¨ulich storage ring under COSY Infinity was at first thought to be due to a coding error. After numerous checks failed to uncover any obvious mistakes, an investment was made in procuring and mastering the code ZGOUBI. Only after this second code produced the same instability did it become apparent that we were seeing an actual physical phenomenon. This set the course for investigation, mainly searching for a stable configuration under the constraints presented. But the advantages went beyond this single problem. ZGOUBI became another tool for research, always available to provide an independent perspective on a problem. When it came time to match the application of fringe fields between COSY INFINITY and ZGOUBI, it became necessary to study both methods of implementation, which were quite different. The impact of symplectic tracking on the fidelity of a simluation also became evident during this investigation. As a researcher, seeing widely varying approaches to the same physical problem enhances one’s understanding greatly, an additional benefit. Applying the same multi-simulation methodology to the GSI HESR led to useful results in a much shorter time. The learning curve had already been surmounted, and lessons learned from previous work were brought to bear on a new, much larger, much more complex system. Extending the multi-simulation paradigm to include MAD8 and MADX brought to light the severe limitations that ignoring fringe fields presents. Benchmarking need not be connected to an actual scientific endeavor. As shown in Chapter 184 4, it is useful to compare the performance of various codes in a series of simple, well-defined tests. There we showed the various fidelities of COSY INFINITY, ZGOUBI, MAD8, and GICOSY when simulating simple systems such as dipoles with edge angles, fringe effects in elements of varying aperture, and the dynamic apertures of simple storage rings. In Chapter 5, we demonstrate an additional benefit of crosscode benchmarking - revealing coding errors. Code testing is an integral part of development, and too often it is difficult to keep up the required never-ending stream of testing. This is even more true if the code is under active expansion, and used by many institutions. Cross-checking codes against each other as part of the scientific methodology helps the developer test his or her code. The investigator becomes a de facto code tester, uncovering bugs and reporting them back to the developer. The benefits of crosscode testing is not limited to testing codes alone, however. In Chapter 6, we compare code results to theory. The benefits to this line of research are many. If a simulation code produces results that do not agree with theory, what to do? Clearly, one must check the theory and the code. Reproducibility is part of the very fabric of science, and when theory does not match experiment (even the simulated kind), this produces the very best kind of science. Even in Chapter 7, where we are modeling an existing system with already proven code, we benefitted by comparing our results with a simple Runge-Kutta advanced-based code, to boost understanding and confidence in results. This was critical due to the fact that such simulations under conditions of intense space charge had never been attempted. Prior to running full simulations with multiparticle bunches, a substantial amount of time and effort was devoted to running simpler simulations with small numbers of particles. In this way, an intuition for the behavior of the system was developed. As the models grew more sophisticated, new results were compared with baseline results already established. Knowing how the system behaved without space charge, one could clearly see which new 185 effects were due to space charge alone. Any unexpected phenomena were scrutinized against the backdrop of expected results. Furthermore, it was important to make sure that any results could be easily reproduced by other researchers. This meant that any simulations had to be kept straightforward, with as few ad hoc conditions as possible. An operating point capable of containing a beam up to 108 particles was empirically found, but was considered only as a proof of principle. Work began immediately to duplicate these results with more realistic electrode biases. A heuristic development cycle evolved, comprised of increasing the particle counts, increasing the electrode biases, and maintaining energy isochrony and parallel to point focusing. It was discovered that unless these criteria were adhered to, the nonlinearities induced by space charge overwhelmed the longitudinal separation of particles which could be acheived by mass difference alone. The transverse and longitudinal optics must be tuned to keep these nonlinearities in check in order that mass separation can be achieved. In conclusion, the benefits of running multiple codes against a physical problem are numerous. Predominant among these benefits are cross-checking, as well as an enhanced physical understanding of the processes involved. This Thesis has demonstrated numerous examples of both of these. As stated in the introduction to this document, a researcher will typically invest in one or two reliable codes, spend time learning them well, and stand by them. In this document, we have demonstrated that only by subjecting a scientific endeavor to a wide variety of codes can one expose every facet of the physics behind the simulation. It is well worth the time and effort invested. 186 Chapter 9 Talks and Publications Published Works • “A Comparison of Storage Ring Modeling with COSY Infinity, ZGOUBI, and MAD8” – Proceedings, Ninth Intl Conference on Charged Particle Optics, Brno, Czech Re- public, 2014 • ”Implementation of Benchmarks for Precision Tracking in Storage Rings” – Proceedings, ICAP 2015, Shanghai,China • ”Comparison of Tracking Codes for the Determination of Dynamic Aperture in Storage Rings” – Proc. 7th Int. Particle Accelerator Conf. (IPAC’16), Busan, S. Korea, 2016 Talks • Implementation of Benchmarks for Precision Tracking in Storage Rings – Talk, ICAP 2015, Shanghai, China • High-Order Modeling of Fringe Field Effects in Storage Rings – Talk, CPO9 (Charged Particle Optics 9) 2014, Brno, Czech Republic 187 • *Nonlinear Dynamics & Resonances in the g-2 Ring – Talk, g-2 Collaboration Meeting, Fermilab, April 2016 *co-author 188 APPENDICES 189 APPENDIX A Fringe Fields “Constant” and “linearly increasing” fields are idealizations. Every physical element has a beginning and an end, and the fields do not begin and end abruptly. Within a constant magnetic field, a particle will travel in a perfect circle, but when the field begins to weaken at the boundaries, the curvature of the particle’s trajectory will decrease. Often this alteration in the motion significantly alters the motion of the particle from the designed intent. Another effect arises from one of Maxwell’s static field equations, in two dimensions within a vacuum, gives ∇ × (cid:126)B = 0 ∂Bx ∂y = ∂By ∂z . From this we see that if the Bx component of a field varies in the y direction, for example, the By component will likewise vary in the x direction. This additional field component will affect particle motion in ways that might not have been intended. These effects are referred to collectively at “fringe effects”. Fringe effects are difficult to calculate analytically, but their importance makes it im- perative to study their effects. There is a simple case which can be solved in closed form, namely the fringe fields of an ideal dipole. 190 An ideal dipole can be modelled with two identical, infinitely large current sheets, with oppositely directed currents (Fig. A.1). Figure A.1: Ideal dipole between two oppositely directed, unbounded current sheets. Outside the sheets, the fields exactly cancel. Between the sheets, they reinforce to form a perfectly constant dipole field. The problem is we want a particle to be able to pass through this constant dipole field. This requires a gap in the sheets. We can introduce such a gap by removing a section of the sheet, but there is an equivalent way. Rather than remove a slice of a sheet, we can superimpose a slice of oppositely directed current onto the sheet. This will cancel the effect of both currents, and mimic the effect of removing a slice. These slices will introduce their own field which is shown in the following Poisson [4] output: This field is nonzero outside of the sheets, and nonuniform between the sheets. Figure A.2: Field of oppositely directled current slices. When this field is added to the idealization of the infinite current sheets, we can clearly see 191 the nature of the fringe fields. Since the model is so simple, it can be described analytically as follows. From the Biot-Savart law [60] we have, for one of the sheets, (cid:90) (cid:126)B = µ0 4π (cid:126)K((cid:126)r(cid:48)) × ˆR S R2 da(cid:48) = µ0K 4π (cid:90) ∞ (cid:90) y0/2 −(z − z(cid:48))ˆy + (y − y(cid:48))ˆz dy(cid:48)dx(cid:48) (cid:90) y0/2 −y0/2 µ0K 4π −∞ −y0/2 (cid:2)−(z − z(cid:48))ˆy + (y − y(cid:48))ˆz(cid:3) dy(cid:48)(cid:90) ∞ R3 dx(cid:48) R3 −∞ The ˆy and ˆz integrals can be evaluated separately. (cid:90) ∞ −∞ dx(cid:48) R3 = = (cid:90) ∞ −∞ −µ0K 2π (z − z(cid:48))ˆy (cid:90) y0/2 −y0/2 µ0K 4π dx(cid:48) (cid:2)(x − x(cid:48))2 + (y − y(cid:48))2 + (z − z(cid:48))2(cid:3)3/2 (cid:90) y0/2 dy(cid:48) (y − y(cid:48))2 + (z − z(cid:48))2 = −y0/2 (y − y(cid:48))dy (y − y(cid:48))2 + (z − z(cid:48))2 = µ0K 4π We can summarize these results by (cid:19) (cid:18)y − y(cid:48) z − z(cid:48) ˆy + (cid:104) µ0K 4π ln (cid:126)B = µ0K 2π arctan = 2 (y − y(cid:48))2 + (z − z(cid:48))2 −y0/2 (cid:16) arctan z − z(cid:48) µ0K 2π (cid:19)(cid:12)(cid:12)(cid:12)(cid:12)y0/2 (cid:18) y − y(cid:48) (y − y(cid:48))2 + (z − z(cid:48))2(cid:17)(cid:12)(cid:12)(cid:12)(cid:12)y0/2 (cid:12)(cid:12)(cid:12)(cid:12)y0/2 (y − y(cid:48))2 + (z − z(cid:48))2(cid:105) ˆz ln −y0/2 −y0/2 There are two ways in which these fringe fields alter particle motion from the design trajec- tory. First, as the particle approaches what will eventually become a uniform magnetic field, it traverses a region where the field is growing from null to its final value. This will cause the path of the particle to curve - slightly at first, but with increasing curvature. Finally, when the particle has reached the uniform field, the curvature of the path will be maximal and constant. 192 The second way in which a fringe field can affect particle motion arises from the fact that there is a longitudinal component to the field off the midplane. This will affect the oblique motion of the particle off the midplane. Let’s discuss each of these effects in turn. From s = rθ and mv2/r = qvBtrans we find dθ = qBtrans p ds Integrating we find the total angular deflection is Btransds fringe (cid:90) θ = q p Hence the total angular deflection depends on the integral of the magnetic field, rather than its functional value. We can therefore replace a fringe field with a constant field of appropriate length to cause the exact same angular deflection. This appropriate length is called the effective length and is defined as that length of constant field that will yield the exact same angular deflection as a varying field. To analyze the variation in the trajectory which occurs off the midplane, refer to Fig. A.3: 193 Figure A.3: Top view of exit fringe field. In the direction σ, we assume that the dipole field decreases linearly as By = B0 − B0 σ l 194 Because the vertical component is changing we have from ∇ × (cid:126)B = 0 ∂By ∂σ = ∂Bσ ∂y . So we will have a Bσ component off the miplane: Bσ ≈ ∂Bσ ∂y y = ∂By ∂σ y = −B0 l y y. l Finally, this gives So Btrans = Bσ sin α = − B0 sin α (cid:90) l/ cos α θ = q p 0 (cid:18) −B0 sin α l (cid:19) yds = − q p B0y tan α. (A.1) In words, this means that an off-midplane particle with an oblique trajectory will feel a vertical deflection due to fringe fields. 195 APPENDIX B Second Order Beam Dynamics Linear dynamics forms the foundation of Beam Physics. The challenge comes from higher order dynamics. A pioneering attempt to investigate higher-order beam physics was under- taken by Karl Brown [5] who developed the following program. First, the differential equation for the trajectory of a charged particle in a static magnetic field with midplane symmetry is derived. The solution of this equation is assumed to be a Taylor series solution about a central trajectory (the reference orbit). Terms to second order in the initial conditions are retained. The coefficients of this Taylor series expansion are seen to be solutions to a harmonic oscillator equation. The first order coefficients are nothing more than sines and cosines. The coefficients of the second order terms are also harmonic oscillators, but these are driven. The Green’s function technique is used to solve these nonhomogeneous equations. Finally, we discover that we can express the second-order results in terms of the first-order coefficients. We begin with the traditional equation of motion ˙(cid:126)p = e((cid:126)v × (cid:126)B) 196 Letting (cid:126)T be the position vector of the particle, we can rewrite this equation without time: (cid:32) (cid:33) (cid:32) d dT d (cid:126)T dT p = e d (cid:126)T dT × (cid:126)B (cid:33) or (cid:32) (cid:33) d (cid:126)T dT × (cid:126)B d2 (cid:126)T dT 2 + d (cid:126)T dT p dp dT = e (B.1) But the right hand side of (Eq. B.1) is perpendicular to d (cid:126)T /dT , as is d2 (cid:126)T /dT 2. So dp/dT = 0 which gives: d2 (cid:126)T dT 2 = e p d (cid:126)T dT × (cid:126)B (cid:32) (cid:33) (B.2) We want to transform equation (Eq. B.2) into particle optical coordinates as shown in Fig. B.1. This coordinate system is the traditional one used in beam physics. The starting point is to pick an arbitrary trajectory and call that the reference trajectory. In the figure this is called s. While this trajectory is arbitrary, it is usually chosen as the design trajec- tory. For example, if there is a quadrupole in the system, the reference trajectory will go straight along the axis of the quadrupole. The reference trajectory is always assumed to lie in a plane. For most applications, this is not a limitation. Each point of the reference orbit has a coordinate triad affixed to it. If there is curvature at a point, the unit vector normal to that point away from the curvature is called ˆx and is one of the coordinates. The tangent to s is another coordinate, called ˆs. Finally, since the reference path is assumed always to be in a plane, the third coordinate y is perpendicular to the plane. It is these coordinates x, y, s that we use to specify the location of the particle. Also, since the reference orbit can be curved, we will have need for the curvature h, defined by the reciprocal of the radius of curvature. 197 Figure B.1: An excerpt from Brown’s seminal paper [5], showing the traditional coordinate system. 198 We transform the equation of motion (Eq. B.2) into these new coordinates. The first step is to change the independent variable from T , which is the arc length along the path of the particle in question, to s, the arc length along the reference orbit. This gives (cid:126)T(cid:48)(cid:48) − 1 2 (cid:126)T(cid:48) T(cid:48)2 dT(cid:48)2 ds = e p T(cid:48)( (cid:126)T(cid:48) × (cid:126)B) To bring in the x and y coordinates, we make the following substitutions: 1 2 = x(cid:48)x(cid:48)(cid:48) + y(cid:48)y(cid:48)(cid:48) + (1 + hx)(hx(cid:48) + h(cid:48)x) T(cid:48)2 = x(cid:48)2 + y(cid:48)2 + (1 + hx)2 (h = local curvature) dT(cid:48)2 ds (cid:126)T(cid:48) = ˆxx(cid:48) + ˆyy(cid:48) + (1 + hx)ˆx (cid:126)T(cid:48)(cid:48) = ˆx(x(cid:48)(cid:48) − h(1 + hx)) + ˆyy(cid:48)(cid:48) + ˆs(2hx(cid:48) + h(cid:48)x) After these substitutions are made, we get the following very complicated-looking vector formula (in component form): ˆx component: x(cid:48)(cid:48) − h(1 + hx) − x(cid:48) T(cid:48)2 (x(cid:48)x(cid:48)(cid:48) + y(cid:48)y(cid:48)(cid:48) + (1 + hx)(hx(cid:48) + h(cid:48)x) = T(cid:48)(y(cid:48)Bs − (1 + hx)By) e p ˆy component: y(cid:48)(cid:48) − y(cid:48) T(cid:48)2 (x(cid:48)x(cid:48)(cid:48) + y(cid:48)y(cid:48)(cid:48) + (1 + hx)(hx(cid:48) + h(cid:48)x)) = T(cid:48)((1 + hx)Bx − x(cid:48)Bs) e p It is important to point out that these are still exact equations. They are nothing more than the original Lorentz equation of motion expressed in the curvilinear coordinates employed in beam physics. 199 To proceed, we multiply the expressions out, and expand the 1/T(cid:48)2 in a Taylor series. Retaining terms to second order only (this is the first instance of approximation), we find: x(cid:48)(cid:48) − h(1 + hx) − x(cid:48)(hx(cid:48) + h(cid:48)x) = y(cid:48)(cid:48) − y(cid:48)(hx(cid:48) + h(cid:48)x) = T(cid:48)(y(cid:48)Bs − (1 + hx)By) T(cid:48)((1 + hx)Bx − x(cid:48)Bs)) e p e p (B.3) So far we have said nothing about the magnetic field Bx, By, Bs. Neglecting the beam itself, this field is current free, therefore (cid:126)B = ∇φ, for some magnetic potential φ. The Laplacian in our beam physics coordinate system is: (cid:18) (cid:19) ∇2φ = 1 1 + hx ∂ ∂x (1 + hx) ∂φ ∂x + ∂2φ ∂y2 + 1 1 + hx ∂ ∂s (cid:18) 1 1 + hx (cid:19) ∂φ ∂s = 0 If we assume that the magnetic scalar potential has odd symmetry about the midplane, then we can make the expansion φ(x, y, s) = ∞(cid:88) ∞(cid:88) m=0 n=0 A2m+1,n xn n! y2m+1 (2m + 1)! . Unfortunately, when we substitute this into the Laplacian, the result is not very illuminating - −A2m+3,n = A(cid:48)(cid:48) 2m+1,n + nhA(cid:48)(cid:48) 2m+1,n−1 −nh(cid:48)A(cid:48) 2m+1,n−1 + A2m+1,n+2 + (3n + 1)hA2m+1,n+1+ n(3n − 1)h2A2m+1,n + N (n − 1)2h3A2m+1,n−1 + 3nhA2m+3,n−1 +3n(n − 1)h2A2m+3,n−2 + n(n − 1)(n − 2)h3A2m+3,n−3. 200 Nonetheless, we will only require the first relation with m = n = 0. −A3,0 = A(cid:48)(cid:48) 1,0 + A1,2 + hA1,1 This is because we restrict ourselves to second order only: φ(x, y, s) = (A10 + A11x + 1 2! A12x2)y + A30 y3 3! And we see we have all the coefficients that we need. Turning our attention now to the magnetic field, By(x, y, s) = ∂φ ∂y = A10 + A11x + 1 2! x2 + 1 2! A3,0y2. Because of the odd symmetry of the potential with respect to the midplane, on the midplane we have only By. Hence By(x, 0, s) = A10 + A11x + 1 2! A12x2. The Aij are functions of s alone. So we are free to set x = 0. This means that A10 = By(0, 0, s) , A11 = ∂By ∂x , A12 = ∂2By ∂x2 To avoid writing derivatives, we adopt the following notation: By(x, 0, s) = By(0, 0, s)(1 − nhx + βh2x2). 201 Using these values instead of Aij, and setting By(0, 0, s) = hp0/e we can write (−nh2y + 2βh3xy) (h − nh2x + βh3x2 − 1 2 (h(cid:48)y − (n(cid:48)h2 + 2nhh(cid:48) + hh(cid:48))xy) (h(cid:48)(cid:48) − nh3 + 2βh3)h2) (cid:16)p0 (cid:16)p0 (cid:16)p0 e e (cid:17) (cid:17) (cid:17) Bx(x, y, s) = By(x, y, s) = Bs(x, y, s) = e These are the second-order expressions for the magnetic field that we will use. We plug these magnetic fields into the 2nd order equations of motion B.3 to obtain x(cid:48)(cid:48) + (1 − n)h2x = gδ + (2n − 1 − β)h2h3x2 + h(cid:48)xx(cid:48) + x(cid:48)2 (h(cid:48)(cid:48) − nh3 + 2βh3)y2 1 2 hy(cid:48)2 − hδ2 +h(cid:48)yy(cid:48) − 1 2 y(cid:48)(cid:48) + nh2y = 2(β − n)h3xy + h(cid:48)xy(cid:48) − h(cid:48)x(cid:48)y + hx(cid:48)y(cid:48) + nh2yδ. +(2 − n)h2xδ + 1 2 Note: in the above we used the expansions (cid:113) T(cid:48) = x(cid:48)2 + y(cid:48)2 + (1 + hx)2 and p0 p = p0 p0(1 + δ) = 1 − δ + δ2. This is a very important result - these are nothing more than the equations of harmonic oscillation - with a wide variety of driving terms! In fact, from Wiedemann [2], the same equations of motion are worked out to third order (Fig. B.2), and are nothing more than the equation of motion for a harmonic oscillator, with a wide variety of driving terms. Continuing with the second order analysis, the x(s) motion 202 Figure B.2: Excerpt from Wiedemann 203 of a particle can be expanded as +(x|x2 x = (x|x0)x0 + (x|x(cid:48) 0)x0x(cid:48) 0 + (x|x0x(cid:48) 0)x2 (cid:48)2 (cid:48)2 0δ)x(cid:48) 0 + (x|x(cid:48) +(x|x 0 )x 0)y0y(cid:48) 0 + (x|y0y(cid:48) +(x|y0)y2 0)x(cid:48) 0 + (x|δ)δ 0 + (x|x0δ)x0δ 0δ + (x|δ2)δ2 (cid:48)2 0 + (x|y (cid:48)2 0 )y 0 We plug this into Brown’s 2nd order equation of motion: x(cid:48)(cid:48) + (1 − n)h2x = hδ + (2n − 1 − β)h2h3x2 + h(cid:48)xx(cid:48) + +(2 − n)h2xδ + 1 2 x(cid:48)2 (h(cid:48)(cid:48) − nh3 + 2βh3)y2 1 2 hy(cid:48)2 − hδ2 +h(cid:48)yy(cid:48) − 1 2 We equate like terms. For an example, if we equate x0 terms we get: (x|x2 0)(cid:48)(cid:48) + (1 − n)h2(x|x2 0) = (2n − 1 − β)h3(x|x0)2 + h(cid:48)(x|x0)(x|x0)(cid:48) + (cid:48)2 h(x|x0) 1 2 To simplify matters, we assume n, β, and h are constant: (x|x2 0)(cid:48)(cid:48) + (1 − n)h2(x|x2 0) = (2n − 1 − β)h3(x|x0) + (cid:48)2 h(x|x0) 1 2 (B.4) We find the particular solution with the Green’s function, (cid:90) t (cid:20) (x|x2 0) = 0 (cid:21) G(t, τ )dτ. (2n − 1 − β)h3(x|x0)2 + (cid:48)2 h(x|x0) 1 2 204 substitutions and obtain (x|x2 0) =(cid:90) t (cid:20) = + 1 3 h 6 This gives (cid:21) h(cos kxτ ) (cid:48)2 1 2 (2n − 1 − β)h3(cos kxτ )2 + 0 = (2n − 1 − β)h3 cos2 kxτ G(t, τ )dτ + 1 2 hk2 x [sinkxt + (1 − cos kxt)] + (2n − 1 − β) (cid:2)2(1 − cos kxt) − sin2 kxt(cid:3) . (cid:90) t 0 h3 k2 x G(t, τ )dτ (cid:90) t 0 1 6 k2 xh sin2 kxτ G(t, τ )dτ Here G(t, τ ) is the Green’s function of B.4 From the first order solutions, we can make the (x|x0) = cos(kxt) and (x|x0)(cid:48) = −kx sin(kxt), (cid:104) (cid:105) 2(1 − cos kxt) − sin2 kxt . (x|x2 0) = 1 3 (2n − 1 − β) h3 k2 x [sin kxt + (1 − cos kxt)] + h 6 205 APPENDIX C Second Order Analysis of Electrostatic Elements (Newtonian) Relativistic case A separate method of derivation is given in Berz, Makino, Wan (BMW) [19]. We also have a method of solution so we can pinpoint the discrepancy between (a|xx) and (a|ax). The equations of motion we are interested in analyzing are x(cid:48) = a(1 + hx) a(cid:48) = (1 + hx) , p0 ps (cid:20) 1 + η 1 + η0 Ex χe0 p0 ps + b Bs χm0 p0 ps − By χm0 (cid:21) + h ps p0 (C.1) We need to work with Eq. C.1, which is quite complicated. So at this time we perform some simplifications. Let us call(cid:20) 1 + η 1 + η0 (cid:21) = M. Ex χe0 p0 ps + b Bs χm0 p0 ps − By χm0 Then we can write a(cid:48) = (1 + hx)M + h ps p0 . (C.2) We need the first order Taylor expansion of Eq. C.2. This involves substantial computation. 206 The first order coefficient with respect to x is given by ∂a(cid:48) ∂x = hM + (1 + hx) ∂M ∂x + h ∂ ∂x (cid:18) ps (cid:19) p0 . First, let us evaluate ∂/∂x(ps/p0). From equation (3.18) of BMW,[19] we have (cid:18) η(2 + η) η0(2 + η0) (cid:19)−1/2 − a2 − b2 = (N − a2 − b2)−1/2. p0 ps = We evaluate (cid:18) ps (cid:19) p0 d dx = (N − a2 − b2)−1/2 dN dx 1 2 . Now from page 68 of BMW, [19], (cid:19)(cid:18) (cid:18) N = (cid:19) (cid:18) (cid:19) 1 + δ − ZeEx0 η0mc2 x η0 = 1 + 1 + 1 + 2 + η0 η0 (cid:18) 1 δ − ZeEx0 (2 + η0)mc2 x (cid:19) δ − ZeEx0 mc2 + 1 η0 x + η0 2 + η0 δ2 − 2eZEx0 (2 + η0)mc2 xδ 2 + η0 2 + η0 + Z2e2E2 x0 η0(2 + η0)m2c4 x2, and so (cid:18) 1 2 + η0 + 1 η0 (cid:19) dN dx = −ZeEx0 mc2 We will also need N (0) = 1 and − 2eEx0 (2 + η0)mc2 δ + 2 Z2e2E2 x0 η0(2 + η0)m2c4 x. (cid:18) 1 2 + η0 (cid:19) . + 1 η0 dN (0) dx = −ZeEx0 mc2 207 Next we need to evaluate dM/dx with (cid:20) 1 + η 1 + η0 M = Ex χe0 p0 ps + Bs χm0 p0 ps b − By0 χm0 (cid:21) (1 + nbx) Here M is seen to consist of three terms. So it makes sense to split it into three components and work with each separately, that is, MI = 1 + η 1 + η0 Ex χe0 p0 ps p0 ps b MII = Bs χm0 MIII = − By0 χm0 (1 + nbx) First, work with MI using relations from page 68 of BMW [19]: MI = 1 + η0 1 + η0 δ − ZeEx0 (1 + η0)mc2 x (1 + nex)(N − a2 − b2)−1/2 Differentiating MI with respect to x: ∂MI ∂x = − ZeEx0 (1 + η0)mc2 −Ex0 χe0 (1 + nex)(N − a2 − b2)−1/2 η0 1 + η0 η0 δ − ZeEx0 (1 + η0)mc2 x δ − ZeEx0 (1 + η0)mc2 x 1 + 1 + η0 ne(N − a2 − b2)−1/2 (1 + nex)(N − a2 − b2)−3/2 ∂N ∂x Next we work with MII: MII = Bs χm0 p0 ps b = Bs χm0 (N − a2 − b2)−1/2b 208 (cid:20) (cid:20) (cid:20) − 1 2 + 1 + (cid:21) −Ex0 χe0 (cid:21) −Ex0 (cid:21) −Ex0 χe0 χe0 and Similarly, ∂MII ∂x = −1 2 Bs χm0 ∂N ∂x . (N − a2 − b2)−3/2b (cid:19) (cid:18) ∂MIII ∂x = d dx −By0 χb0 (1 + nbx) = − By0 χm0 nb. This completes all of the derivatives with respect to x. We now must evaluate these deriva- tives at x = a = b = 0. Going back to the very beginning, we have a(cid:48) = (1 + hx)M + h ps p0 , da(cid:48) dx = hM + (1 + hx) dM dx + h d dx (cid:18) ps (cid:19) p0 . All of these quantities need to be evaluated at x = a = b = 0. What is M (0)? (cid:20) (cid:21) Ex0 χe0 MI(0) = 1 + η0 1 + η0 δ − ZeEx0 (1 + η0)mc2 x (1 + nex)(N − a2 − b2)−1/2 = −Ex0 χe0 because N (0) = 1. Continuing with MII (0), MII(0) Bs χm0 p0 ps = Bs χm0 (N − a2 − b2)−1/2b = Bs χm0 b = 0, and finally (cid:18) By0 (cid:19) χm0 MIII(0) = − (1 + nbx) = − By0 χm0 209 So we have MI(0) + MII(0) + MIII (0) = M (0) = −Ex0 χe0 + 0 − By0 χm0 = −h. Repeating for emphasis, M (0) = −h. Next, we need to evaluate ∂M/∂x at x = a = b = 0. There are three parts to this: ∂MI ∂x = − ZeEx0 (1 + η0)mc2 −Ex0 χe0 (1 + nex)(N − a2 − b2)−1/2 (cid:20) (cid:20) − 1 2 + 1 + η0 1 + η0 η0 δ − ZeEx0 (1 + η0)mc2 x δ − ZeEx0 (1 + η0)mc2 x 1 + 1 + η0 ne(N − a2 − b2)−1/2 (1 + nex)(N − a2 − b2)−3/2 dN dx (cid:21) −Ex0 (cid:21) −Ex0 χe0 χe0 to evaluate this, recall that N (0) = 1 and ∂N (0) ∂x = −2ZeEx0 mc2 1 + η0 η0(2 + η0) . So and giving, ∂MI(0) ∂x = − ZeEx0 (1 + η0)mc2 −Ex0 χe0 −Ex0 χe0 + ne − Ex0 χe0 ZeEx0 mc2 1 + η0 η0(2 + η0) , ∂MII(0) ∂x = −1 2 Bs χm0 (N − a2 − b2)−3/2b ∂N ∂x = 0, ∂MIII(0) ∂x = − By0 χm0 nb. 210 Putting these all together we get ∂M (0) ∂x = − ZeEx0 (1 + η0)mc2 = ZeEx0 (1 + η0)mc2 + −Ex0 χe0 Ex0 χe0 −Ex0 χe0 ne − Ex0 χe0 ne − Ex0 χe0 ZeEx0 mc2 − Ex0 χe0 ZeEx0 mc2 1 + η0 η0(2 + η0) + 0 − By0 χm0 nb 1 + η0 η0(2 + η0) − By0 χm0 nb. We are trying to evaluate ∂a(cid:48) ∂x = hM (0) + ∂M (0) ∂x + h ∂ ∂x (cid:18) ps (cid:19) p0 and so far we have M (0) and ∂M (0)/∂x. Last we need (cid:18) ps (cid:19) p0 ∂ ∂x = −h ZeEx0 mc2 1 + η0 η0(2 + η0) So we have the result ∂a(cid:48) ∂x = = −h2 − Ex0 χe0 = −h2 − Ex0 χe0 = −h2 − Ex0 χe0 ne − By0 χm0 ne − By0 χm0 ne − By0 χm0 nb + nb + nb + ZeEx0 mc2 ZeEx0 mc2 ZeEx0 mc2 (cid:20) Ex0 (cid:20) Ex0 χe0 (cid:18) 1 (cid:18) (cid:20)Ex0 1 + η0 χe0 1 + η0 − 1 + η0 η0(2 + η0) −1 (cid:18) −1 − h 1 + η0 η0(2 + η0) 1 + η0 (cid:21) − h (cid:21) − h η0(1 + η0)(2 + η0) η0(2 + η0) η0(2 + η0) χe0 (1 + η0)2 (cid:19) (cid:19) (cid:19) (cid:21) which matches Eq. 4.2 in BMW precisely. Next we want to go to second order, and only for the electrostatic case. The first step is to use a spherical field: Ex = − χe0R (R + x)2 211 Using the same nomenclature as before, and introducing a new symbol Q = Ze mc2 , We also want to evaluate in advance 1 + η 1 + η0 = 1 − QV 1 + η0 , (cid:18) 1 − QV η0 N = η(2 + η) η0(2 + η0) = 1 + η 1 + η0 ∂ ∂x (cid:19)(cid:18) 1 − QV 2 + η0 = − Q 1 + η0 ∂V ∂x (cid:19) = 1 − 2Q , (cid:18) 1 + η0 (cid:19) η0(2 + η0) V + Q2V 2 η0(2 + η0) , = − ∂N ∂x 2Q η0(2 + η0) (1 − η0 − QV ) ∂V ∂x , A problem with these expressions is that η0 is in the denominator. This will be addressed when we analyze the non-relativistic limit. (N − a2)−1/2 + h(N − a2)1/2 (cid:19)(cid:18)−Ex (cid:19) a(cid:48) = (1 + hx) = (1 + hx) (cid:18) −Ex 1 + η 1 + η0 χe0 1 − QV 1 + η0 (N − a2)−1/2 + h(N − a2)1/2 χe0 212 (cid:18) (cid:18) (cid:18) 1 − QV 1 + η0 + (1 + hx) + (1 + hx) + (1 + hx) 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 h(N − a2)−1/2 ∂N ∂x + 1 2 (N − a2)−1/2 (cid:19) (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18) ∂ (cid:19)(cid:18) (cid:19)(cid:18) Ex χe0 −Ex χx0 ∂x (cid:19) (N − a2)−1/2 (N − a2)−1/2 − 1 2 (N − a2)−3/2 ∂N ∂x ∂a(cid:48) ∂x = (cid:18) h (cid:19) (cid:19)(cid:18)−Ex χe0 − Q ∂V ∂x χe0 213 (cid:18) 1 − QV 1 + η0 (cid:19) (cid:19)(cid:18) ∂ ∂x −Ex χe0 (N − a2)−1/2 (N − a2)−3/2 ∂N ∂x (N − a2)−3/2 dN dx (N − a2)−3/2 ∂N ∂x (N − a2)−3/2 ∂N ∂x (N − a2)−5/2 ∂N ∂x − 3 2 − 1 2 (N − a2)−3/2 ∂2N − 1 ∂x2 2 h(N 2 − a2)−1/2 ∂2N ∂x2 1 2 ∂2a(cid:48) ∂x2 = (cid:18) h − Q (cid:18) (cid:18) ∂V 1 + η0 ∂x 1 − QV 1 + η0 ∂V ∂x − Q − Q 1 + η0 + h + h + (1 + hx) + (1 + hx) + (1 + hx) + h 1 − QV 1 + η0 (N − a2)−1/2 1 + η0 − Q 1 + η0 ∂V ∂x − Q 1 + η0 (N − a2)−3/2 ∂N ∂x χe ∂x χe0 (cid:19) (cid:19) (cid:19) − 1 2 − 1 2 (N − a2)−1/2 χe0 ∂2V ∂x2 (N − a2)−1/2 + h (N − a2)−1/2 (cid:19)(cid:18) ∂ χe0 −Ex χe0 ∂V ∂x −Ex χe0 (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18)−Ex (cid:19)(cid:18) (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18)−Ex (cid:19)(cid:18) ∂ (cid:19) (cid:19)(cid:18) (cid:19)(cid:18)−Ex (cid:19) (cid:19)(cid:18) ∂ (cid:19)(cid:18) ∂2 (cid:19) (cid:19)(cid:18) ∂ (cid:19)(cid:18) 1 (cid:19) (cid:19)(cid:18) (cid:19) (cid:19)(cid:18) Ex (cid:19)(cid:18) (cid:19)(cid:18) (cid:19) (cid:19)(cid:18) ∂ (cid:19)(cid:18) (cid:19)(cid:18) Ex (cid:19)(cid:18) (cid:19) (cid:19)(cid:18) (cid:19)(cid:18) Ex (cid:19)(cid:18) Ex − 1 2 − 1 2 χe0 ∂V ∂x ∂x − 1 2 (cid:19) χe0 (N − a2)−1/2 −Ex χe0 Ex χe0 Ex χe0 (N − a2)−3/2 ∂N ∂x (N − a2)−1/2 (cid:19) 2 Ex χe0 ∂V ∂x (cid:19) ∂x2 ∂x χe0 χe0 ∂x ∂x x − Q 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 − Q 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 (N − a2)−1/2 χe0 h(N 2 − a2)−3/2 ∂N ∂x + (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:18) (cid:19) (cid:18) (cid:18) + (1 + hx) + (1 + hx) + (1 + hx) + (1 + hx) + (1 + hx) + (1 + hx) + (1 + hx) (cid:18) 1 (cid:19)(cid:18) − 1 2 2 + h 1 − QV 1 + η0 214 = ∂a(cid:48) ∂a∂x (cid:18) h 1 − QV 1 + η0 + (1 + hx) + (1 + hx) (cid:18) + (1 + hx) − ah −1 2 (cid:19) (cid:19)(cid:18)−Ex χe0 ∂V ∂x − QV 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 (cid:18) (cid:18) (cid:18) (cid:19) χe0 (N 2 − a2)−3/2 ∂N ∂x a(N − a2)−3/2 (cid:19) a(N − a2)−3/2 (cid:19) a(N − a2)−3/2 (cid:18) (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18) ∂ (cid:19) (cid:19)(cid:18)−Ex χe0 −Ex χe0 ∂x a −3 2 (N − a2)−5/2 ∂N ∂x and as a check let us evaluate ∂a(cid:48)/∂x∂da. (−1/2)(N − a2)−3/2(−2a) (cid:19)(cid:18)−Ex (cid:19)(cid:18) ∂ (cid:19)(cid:18) (cid:19)(cid:18) Ex χe0 −Ex χx0 ∂x (cid:19) (−1/2)(N − a2)−3/2(−2a) (cid:19) (−1/2)(N − a2)−3/2(−2a) (cid:19) (−5/2)(N − a2)−3/2(−2a) −1 2 ∂N ∂x = ∂2a(cid:48) ∂x∂a (cid:18) h 1 − QV 1 + η0 + (1 + hx) + (1 + hx) + (1 + hx) (cid:19) (cid:19)(cid:18)−Ex χe0 ∂V ∂x − Q 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 (cid:18) (cid:18) (cid:18) χe0 (−1/2)h(N − a2)−3/2(−2a) + 1 2 ∂N ∂x 215 cleaning it up a little bit = ∂2a(cid:48) ∂x∂a (cid:18) h 1 − QV 1 + η0 + (1 + hx) (cid:19) (cid:19)(cid:18)−Ex χe0 − Q ∂V ∂x + (1 + hx) + (1 + hx) 1 + η0 1 − QV 1 + η0 1 − QV 1 + η0 ah(N − a2)−3/2 ∂N ∂x (cid:18) (cid:18) (cid:18) + 1 2 a(N − a2)−3/2 (cid:19) a(N − a2)−3/2 (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18) ∂ (cid:19) (cid:19)(cid:18) Ex ∂x a(N − a2)−3/2 χe0 −Ex χx0 (3/2)a(N − a2)−5/2 ∂N ∂x χe0 (cid:18)1 (cid:19) (N − a2)−1/2(−2a) ∂a(cid:48) ∂a = (1 + hx) (cid:18) = (1 + hx) (cid:18) 1 − QV 1 + η0 1 − QV 1 + η0 (cid:19) χe0 (cid:19)(cid:18)−Ex (cid:19)(cid:18) (cid:19) (cid:19)(cid:18)−Ex (cid:19)(cid:18)−Ex χe0 1 − QV 1 + η0 χe0 ∂2a(cid:48) ∂a∂a = (1 + hx) (cid:18) 2 (N − a2)−3/2(−2a) + h −1 2 a(N − a2)−3/2 − ah(N − a2)−1/2 (cid:19) [(N − a2)−3/2 + 3a2(N − a2)−5/2] − [h(N − a2)−1/2 + ha2(N − a2)−3/2] Next we derive the same derivatives for the nonrelativistic case. 216 Nonrelativistic case Before we can work with the nonrelativistic case, we need to find the correct nonrelativistic forms of x(cid:48) = a(1 + hx) (cid:34) p2 p2 0 (cid:35)−1/2 , − a2 − b2 (cid:34) a(cid:48) = (1 + hx) 1 + η 1 + η0 Ex χe0 − a2 − b2 p2 p2 0 (cid:35)−1/2 (cid:34) p2 p2 0 + h (cid:35)1/2 . − a2 − b2 For nonrelativistic velocities, we have It can also be shown that 1 + η 1 + η0 = γ γ0 = 1 η(1 + η) η0(1 + η0) is the ratio of the test particle kinetic energy to the reference particle kinetic energy and remains such in the limit of nonrelativistic velocities [19]. Using these results we can write the nonrelativistic equations of motion as (cid:20) K (cid:20) K K0 Ex χe0 K0 (cid:21)−1/2 (cid:21)−1/2 − a2 x(cid:48) = a(1 + hx) − a2 a(cid:48) = (1 + hx) (cid:21)1/2 − a2 (cid:20) K K0 + h (C.3) We need K as a function of x. This can be done by considering the potential energy. The function − χe0R R + x 217 is a function with the property that its negative gradient is equal to the electric field at x. We can add any constant to this function without changing that property. Therefore let us add a constant such that when x = 0, the value of the potential is 0. This means that V = −qχe0R R + x + qχe0 = qχe0 x R + x . So the ratio of the kinetic energies which appear in Eq. C.3 takes the form K K0 = K0 − V K0 = 1 − V K0 = 1 − qχe0x K0(R + x) = 1 − qχe0x 1 2qχe0(R + x) = 1 − 2x R + x = R − x R + x . Finally we arrive at the nonrelativistic equations of motion: (cid:20)R − x (cid:20)R − x R + x R + x Ex χe0 − a2 (cid:21)−1/2 (cid:21)−1/2 − a2 x(cid:48) = a(1 + hx) a(cid:48) = (1 + hx) (cid:20)R − x R + x + h (cid:21)1/2 − a2 For the spherical case, the equations become (cid:20)R − x R + x −R (cid:21)−1/2 (cid:20)R − x − a2 (R + x)2 R + x x(cid:48) = a(1 + hx) a(cid:48) = (1 + hx) (cid:21)−1/2 (cid:20)R − x R + x + h (cid:21)1/2 − a2 − a2 To 0th order, these equations are simply x(cid:48) = a a(cid:48) = 0 218 To first order, we need a(cid:48) = h −R − a2 (R + x)2 R + x −R (cid:20) R − x (cid:21)−1/2 (cid:20)R − x + (1 + hx)(−2) (cid:20)R + x −R (R + x)2 (−1 (cid:21)−1/2 −2R − a2 (cid:21)−1/2 (cid:20)R − x (cid:20) R + x (R + x)2 − a2 1 2 = −h + (1 + hx) (R + x)3 ) 2 R + x R + x + h R + x R (R + x)2 (cid:21)−1/2 (cid:21)−3/2 − a2 − a2 R R R + x (cid:20)R − x (cid:20)R + x (cid:21)−1/2 R + x R + x R (R + x)2 . + 2(1 + hx) (R + x)3 − (1 + hx) (cid:20)R + x R + x (R + x)2 − a2 (cid:21)−1/2 (cid:21)−3/2 −2R − a2 − a2 (R + x)2 R (R + x)2 − h Setting a = x = 0 x(cid:48) = a a(cid:48) = − 1 R2 which agrees with the spherical case. In order to derive the nonrelativistic equations of motion for the cylindrical case, we recall the equations of motion Eq. C.3 (cid:20) K (cid:20) K K0 Ex χe0 K0 (cid:21)−1/2 (cid:21)−1/2 − a2 x(cid:48) = a(1 + hx) − a2 a(cid:48) = (1 + hx) (cid:21)1/2 . − a2 (cid:20) K K0 + h 219 The electric field in this case is Ex = − χe0 R + x . As before, we also need to calculate the kinetic energies. First, the potential energy as a function of x is Then and V = qχe0 ln(R + x) − χe0 ln R. K = K0 − V = 1 2 qχe0 − (qχe0 ln(R + x) − qχe0 ln R) K K0 = 1 − 2 ln R + x R . This gives the equations of motion for a cylindrical deflector (cid:21)−1/2 (cid:20) x(cid:48) = a(1 + hx) a(cid:48) = (1 + hx) R + x − a2 (cid:20) 1 − 2 ln −1 R + x R 1 − 2 ln R + x R − a2 . (cid:21)−1/2 (cid:20) + h 1 − 2 ln R + x R − a2 (cid:21)1/2 . To compare with known results we need to linearize these equations of motion. When we linearize the equation for x(cid:48) we get x(cid:48) = a. 220 To linearize the a(cid:48) equation of motion we need ∂a(cid:48) ∂x = −1 R + x h (cid:20) 1 − 2 ln (cid:21)−1/2 − a2 R + x R (cid:20) (cid:21)−1/2 (cid:21)−3/2(cid:18) −2 − a2 R + x (cid:19) + (1 + hx) + (1 + hx) (cid:18)1 (cid:19)(cid:20) + h 2 1 (cid:18) (R + x)2 −1 R + x 1 − 2 ln −1 2 − a2 R + x 1 − 2 ln (cid:19)(cid:20) R 1 − 2 ln R + x (cid:21)−1/2 −2 R . R + x R + x − a2 R This is evaluated at x = a = 0 to get a(cid:48) =1 − 2 R2 x which agrees with the well-known answer. A more general way to show the correctness of this result is to use an electric field of the form (cid:16) (cid:17) . 1 − n x R Ex = −Ex0 We need the potential energy of this field V giving Ex = − ∂V ∂x with (cid:18) x − n x2 2R (cid:19) . V = qEx0 Since this is 0 at the reference orbit, we have (cid:18) x − n (cid:19) . x2 2R K = K0 − V = K0 − qEx0 221 K0, the kinetic energy on the reference orbit is K0 = 1 2 qREx0 and therefore K K0 (cid:16) x − n x2 2R (cid:17) 1 2 qREx0 = 1 − qEx0 (cid:16) (cid:17) . x − n x2 2R R = 1 − 2 So this more general equation of motion becomes (cid:20) K (cid:20) K K0 Ex χe0 K0 (cid:21)−1/2 (cid:21)−1/2 − a2 x(cid:48) = a(1 + hx) − a2 a(cid:48) = (1 + hx) (cid:21)1/2 . − a2 (cid:20) K K0 + h or (cid:20) (cid:18) (cid:19) x(cid:48) = a(1 + hx) a(cid:48) = −(1 + hx) (cid:18) (cid:20) + h 1 − 2 R x − n 1 − 2 R 1 − n (cid:16) − a2 (cid:18) x2 2R 1 − 2 R (cid:17)(cid:20) (cid:21)1/2 x R − a2 . (cid:19) Ex0 χe0 x − n x2 2R (cid:21)−1/2 x − n x2 2R (cid:21)−1/2 (cid:19) − a2 Linearizing the equations of motion we get x(cid:48) = a 222 and evaluating ∂a(cid:48) ∂x = (cid:16) 1 − n − h(1/R) − (1 + hx)(1/R) − (1 + hx)(1/R) (cid:19)1 − 2 (cid:18)1 2 + h −1/2 −1/2 (cid:17) x − n x2 2R − a2 (cid:17) x R (cid:17)1 − 2 (cid:16) (cid:17)1 − 2 (cid:16)− n (cid:17)(cid:18) (cid:16) (cid:16) (cid:17) 1 − n x R x − n x2 R 2R R x − n x2 2R − a2 R R 2R x − n x2 (cid:17) (cid:16) (cid:19)1 − 2 (cid:16) −1/2(cid:18) −1 2 − a2 − a2 R (cid:19)(cid:16) (cid:17) 1 − n x R − 2 R at x = a = 0 to get a(cid:48) = −3 − n R2 which agrees with the well known general solution. Second order solution −3/2(cid:18) − 2 R (cid:19)(cid:16) (cid:17) 1 − n x R We want to consider second order electrostatic fields of the form. (cid:18) 1 − n1 Ex = −Ex0 x R − n2 x2 R2 (cid:19) 223 We need the potential energy V of this field. (cid:18) x − n1 V = qEx0 Since this is 0 at the reference orbit, we have K = K0 − V = K0 − qEx0 And the kinetic energy on the reference orbit is (cid:19) x2 2R − n2 x3 3R2 (cid:18) x − n1 (cid:19) x2 2R − n2 x3 3R2 K0 = 1 2 qREx0 and therefore = 1 − qEx0 K K0 (cid:16) x − n1 2R − n2 x2 x3 3R2 1 2 qREx0 (cid:17) (cid:18) x − n1 = 1 − 2 R (cid:19) x2 2R − n2 x3 3R2 Again repeating the equations of motion Eq. C.3 for ease of reference, (cid:20)K0 (cid:20) K K Ex χe0 K0 (cid:21)−1/2 (cid:21)−1/2 − a2 (cid:21)1/2 − a2 (cid:20) K K0 + h x(cid:48) = a(1 + hx) − a2 a(cid:48) = (1 + hx) 224 and substituting K, K0 and Ex into these equations we get (cid:21)−1/2 (cid:19) − a2 x2 2R − n2 x3 3R2 (cid:20) (cid:18) x − n1 1 − 2 (cid:18) R 1 − n1 x(cid:48) = a(1 + hx) a(cid:48) = −h(1 + hx) (cid:18) x − n1 (cid:20) + h 1 − 2 R x2 2R − n2 (cid:19)(cid:20) (cid:19) x2 2R x2 − n2 R2 x3 3R2 x R − n2 − a2 (cid:21)−1/2 − a2 (cid:19) (cid:18) x − n1 (cid:21)1/2 x3 3R2 1 − 2 R To keep these expressions manageable, let (cid:18) x − n1 (cid:19) − a2 x2 2R − n2 x3 3R2 S = 1 − 2 R and T = 1 − n1 x R − n2 x2 R2 We will also need some related expressions T (0) = 1 −n1 R (0) = dT dx d2T dx2 (0) = −2n2 R2 S(0) = 1 −2 R 2n1 R2 (0) = dS dx d2S dx2 (0) = dS da d2S dxda d2S dada (0) = 0 = 0 = −2 225 Utilizing these simplifications, the equations of motion take the abbreviated form x(cid:48) = a(1 + hx)S−1/2 a(cid:48) = −h(1 + hx)T S−1/2 + hS1/2. In order to expand these equations to second order, we need to evaluate first order derivatives: x(cid:48) = a ∂a(cid:48) ∂x = −h2T S−1/2 − h(1 + hx) S−1/2 + ∂T ∂x h(1 + hx)T S−3/2 ∂S ∂x 1 2 + hS−1/2 ∂S ∂x 1 2 As a check, we verify that this matches the known first order solution. (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)0 da(cid:48) dx x(cid:48) = a = −h2(3 − n1) To continue, we need to the solutions of the first order equations of motion. x(cid:48) = a a(cid:48) = −h2(3 − n1)x This is a well-known result so we can write down the solutions (x|x) = (a|a) = cos (x|a) = R√ 3 − n1 sin √ 3 − n1 R (a|x) = − (cid:16)(cid:112)3 − n1s/R (cid:17) (cid:16)(cid:112)3 − n1s/R (cid:17) (cid:16)(cid:112)3 − n1s/R sin (cid:17) 226 Next we move to evaluate the second order derivatives: ∂a(cid:48) ∂x = −h2T S−1/2 − h(1 + hx) S−1/2 + ∂T ∂x h 2 (1 + hx)T S−3/2 ∂S ∂x + h 2 S−1/2 ∂S ∂x (cid:18) (cid:19) S−1/2 − h2T −1 2 S−1/2 − h(1 + hx) ∂2a(cid:48) ∂x2 = − h2 ∂T ∂x − h2 ∂T ∂x T S−3/2 ∂S ∂x + + h 2 + h2 2 (1 + hx)T S−3/2 ∂2S h 2 ∂x2 + (cid:19) (cid:18) −1 2 (cid:18) ∂2T S−3/2 ∂S ∂x ∂x2 S−1/2 − h(1 + hx) (cid:19)2 (cid:18) ∂S (cid:18) S−3/2 ∂S (cid:19) ∂x S−3/2 ∂T ∂x −1 2 h 2 h 2 ∂x ∂T ∂x + h 2 S−1/2 ∂2S ∂x2 (1 + hx) + (1 + hx)T S−3/2 ∂S (cid:19) ∂x S−5/2 (cid:18) ∂S (cid:19)2 dx −3 2 When this is evaluated at x = a = 0 we get Kaxx = 0 227 We also need the mixed partial derivative ∂a(cid:48) ∂x = −h2T S−1/2 − h(1 + hx) S−1/2 + dT dx h 2 (1 + hx)T S−3/2 ∂S ∂x + h 2 S−1/2 ∂S ∂x ∂2a(cid:48) ∂x∂a (cid:19) (cid:18) = − h2T −1 2 − h(1 + hx) S−3/2 ∂S (cid:18) (cid:19) ∂a (cid:19) (cid:18) −1 2 ∂T ∂x −3 2 + (1 + hx)T (cid:19) (cid:18) h 2 −1 2 h 2 S−3/2 ∂S ∂a ∂S ∂x S−3/2 ∂S ∂a S−5/2 ∂S ∂a S−1/2 ∂2S ∂x∂a ∂S ∂x h 2 + + (1 + hx)T S−3/2 ∂2S ∂x∂a h 2 Because all derivatives with respect to a are 0, we have ∂2a(cid:48) ∂x∂a = Kaxa = 0 a(cid:48) = −h(1 + hx)T S−1/2 + hS1/2 = −h(1 + hx)T ∂a(cid:48) ∂a ∂a(cid:48) ∂a∂a = (cid:19) (cid:18) −1 2 (cid:18) (cid:18) S−3/2 ∂S ∂a (cid:19) (cid:19)(cid:18) (cid:19) −3 2 + h 2 S−1/2 ∂S ∂a (cid:19)2 (cid:18) ∂S ∂a S−5/2 S−3/2 ∂2S ∂a2 S−1/2 ∂2S ∂a2 h 2 + − h(1 + hx)T − h(1 + hx)T −1 2 −1 2 S−3/2 ∂S ∂a (cid:19) (cid:18) + h 2 −1 2 ∂a(cid:48) ∂a∂a = Kaaa = −2h These derivatives, when evaluated at x = a = 0, give us the equation of motion to second 228 order. The idea is to follow BMW [19]. The equation of motion is x(cid:48) = Kxxx + Kxaa + Kxxxx2 + Kxaxax + Kxaaa2 a(cid:48) = Kaxx + Kaaa + Kaxxx2 + Kaaxax + Kaaaa2 To match the notation of BMW, we set x  Kxx Kxa a Kax Kaa  (cid:126)r = ˆM = , and (cid:126)N2 = Kxxxx2 + Kxxaax + Kxaaa2 Kaxxx2 + Kaxaax + Kaaaa2  so we can rewrite the equation of motion as (cid:126)r (cid:48) = ˆM(cid:126)r + (cid:126)N2 Looking only at the first order part, we have a solvable equation of motion x(cid:48) = (x|x)xi + (x|a)ai a(cid:48) = (a|x)xi + (a|a)ai 229 To match the notation of BMW we write (x|x) (a|x)  and (cid:126)ri =  xi ai (x|a) (a|a) ˆL = So the first order solution is The second order solution is ˆL(cid:126)ri x = (x|x)xi + (x|a)ai + (x|xx)x2 a = (a|x)xi + (a|a)ai + (a|xx)x2 i + (x|ax)xiai + (x|aa)a2 i + (a|ax)xiai + (a|aa)a2 i i again, to match the notation of BMW, write R2x  = R2a (x|xx)x2 (a|xx)x2 (cid:126)R2 = i + (x|ax)xiai + (x|aa)a2 i + (a|ax)xiai + (a|aa)a2 i i  Now we substitute this solution into the equation of motion. Since there are two components 230 to the equation of motion, we can do one component at a time. First x: i + (x|ax)(cid:48)xiai + (x|aa)(cid:48)a2 i i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] i + (x|ax)xiai + (x|aa)a2 i ]2 i ]× i + (x|ax)xiai + (x|aa)a2 (x|x)(cid:48)xi + (x|a)(cid:48)ai + (x|xx)(cid:48)x2 = Kxx[(x|x)xi + (x|a)ai + (x|xx)x2 + Kxa[(a|x)xi + (a|a)ai + (a|xx)x2 + Kxxx[(x|x)xi + (x|a)ai + (x|xx)x2 + Kxxa[(x|x)xi + (x|a)ai + (x|xx)x2 [(a|x)xi + (a|a)ai + (a|xx)x2 + Kxaa[(a|x)xi + (a|a)ai + (a|xx)x2 i + (a|ax)xiai + (a|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ]2 If we perform the multiplications in the Kxxx, Kxxa, Kxaa terms, retaining only the second order terms, and combining like terms, i i + (x|ax)(cid:48)xiai + (x|aa)(cid:48)a2 (x|x)(cid:48)xi + (x|a)(cid:48)ai + (x|xx)(cid:48)x2 = Kxx[(x|x)xi + (x|a)ai + (x|xx)x2 + Kxa[(a|x)xi + (a|a)ai + (a|xx)x2 + [Kxxx(x|x)2 + Kxxa(x|x)(a|x) + Kxaa(a|a)2]x2 + [2Kxxx(x|x)(x|a) + Kxxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kxaa(a|x)(a|a)]xiai + [Kxxx(x|a)2 + Kxxa(x|a)(a|a) + Kxaa(a|a)2]a2 i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] i i Taking a look at this rather complicated expression, we see some things that can be rec- ognized. We see ˆL(cid:126)ri. This is multipled by the first row of ˆM . Also we recognize R2x and R2a. But the last three terms are new. We give a new name to these, Q2. But Q2 must be a vector. So we should write (cid:126)Q2, and the portion that we see in this expression is the 231 first component, namely Q2x. For what follows we should write down the second component here. We need to recall the second of the equations of motion, the equation for a(cid:48): a(cid:48) = Kaxx + Kaaa + Kaxxx2 + Kaaxax + Kaaaa2 i + (a|ax)(cid:48)xiai + (a|aa)(cid:48)a2 i i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] i + (x|ax)xiai + (x|aa)a2 i ]2 i ]× i + (x|ax)xiai + (x|aa)a2 (a|x)(cid:48)xi + (a|a)(cid:48)ai + (a|xx)(cid:48)x2 = Kax[(x|x)xi + (x|a)ai + (x|xx)x2 + Kaa[(a|x)xi + (a|a)ai + (a|xx)x2 + Kaxx[(x|x)xi + (x|a)ai + (x|xx)x2 + Kaax[(x|x)xi + (x|a)ai + (x|xx)x2 [(a|x)xi + (a|a)ai + (a|xx)x2 + Kaaa[(a|x)xi + (a|a)ai + (a|xx)x2 i + (a|ax)xiai + (a|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ]2 If we perform the multiplications in the Kaxx, Kaax, Kaaa terms, retaining only the second order terms, and combining like terms, i i + (a|ax)(cid:48)xiai + (a|aa)(cid:48)a2 (a|x)(cid:48)xi + (a|a)(cid:48)ai + (a|xx)(cid:48)x2 = Kax[(x|x)xi + (x|a)ai + (x|xx)x2 + Kaa[(a|x)xi + (a|a)ai + (a|xx)x2 + [Kaxx(x|x)2 + Kaxa(x|x)(a|x) + Kaaa(a|a)2]x2 + [2Kaxx(x|x)(x|a) + Kaxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kaaa(a|x)(a|a)]xiai + [Kaxx(x|a)2 + Kaxa(x|a)(a|a) + Kaaa(a|a)2]a2 i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] i i 232 It is hard to see, but the linear part corresponding to xi and ai is (x|a)(cid:48) (a|a)(cid:48) (a|x)(cid:48) (x|x)(cid:48) (x|x)(cid:48) (a|x)(cid:48)  xi ai (x|a)(cid:48) (a|a)(cid:48) or  xi ai  =  = Kxx Kxa Kxx Kxa Kax Kaa Kax Kaa   (a|x) (x|x) (x|x) (a|x)   (x|a) (a|a) (x|a) (a|a) which agrees with BMW [19] equation (4.36). In that notation, it is written d ˆL ds = ˆM ˆL The expressions corresponding to x2 i , xiai and a2 i (the nonlinear part) are i i + (x|ax)(cid:48)xiai + (x|aa)(cid:48)a2 i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] (x|xx)(cid:48)x2 = Kxx[(x|xx)x2 + Kxa[(a|xx)x2 + [Kxxx(x|x)2 + Kxxa(x|x)(a|x) + Kxaa(a|a)2]x2 + [2Kxxx(x|x)(x|a) + Kxxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kxaa(a|x)(a|a)]xiai + [Kxxx(x|a)2 + Kxxa(x|a)(a|a) + Kxaa(a|a)2]a2 i i We already have names for these quantities, they are R(cid:48) 2x = KxxR2x + KxaR2a + Q2x 233 likewise, looking at the a component of the equation of motion, i + (a|ax)(cid:48)xiai + (a|aa)(cid:48)a2 i i + (x|ax)xiai + (x|aa)a2 i ] i + (a|ax)xiai + (a|aa)a2 i ] (a|xx)(cid:48)x2 = Kax[(x|xx)x2 + Kaa[(a|xx)x2 + [Kaxx(x|x)2 + Kaxa(x|x)(a|x) + Kaaa(a|a)2]x2 + [2Kaxx(x|x)(x|a) + Kaxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kaaa(a|x)(a|a)]xiai + [Kaxx(x|a)2 + Kaxa(x|a)(a|a) + Kaaa(a|a)2]a2 i i or R(cid:48) 2a = KaxR2x + KaaRaa + Q2a These can be combined into a single equation d (cid:126)R2 ds = ˆM (cid:126)R2 + (cid:126)Q2 If (cid:126)Q2 = 0, then we have the same equation as before. So ˆL is the solution of the linear part. We need to find a particular solution to this equation. We try a particular solution of the form (cid:126)R2 = ˆL (cid:126)T 234 Trying this solution gives = ˆM (cid:126)R2 + (cid:126)Q2 d (cid:126)R2 ds d ˆL (cid:126)T = ˆM ˆL (cid:126)T + (cid:126)Q2 ds d ˆL ds (cid:126)T = ˆM ˆL (cid:126)T + (cid:126)Q2 + ˆL + ˆM ˆL (cid:126)T = ˆM ˆL (cid:126)T + (cid:126)Q2 ˆL (cid:126)T ds (cid:126)T ds ˆL (cid:126)T ds = (cid:126)Q2 So we have a way of evaluating the vector (cid:126)T (cid:90) (cid:126)T = ˆL−1 (cid:126)Q2ds has the form (cid:126)R2 =  It is important to spell out exactly the link between T and (a|xa), (a|xx). Recall that (cid:126)R2 R2x  = R2a (x|xx)x2 (a|xx)x2 i + (x|ax)xiai + (x|aa)a2 i + (a|ax)xiai + (a|aa)a2 i i R2a contains (a|xx) and (a|ax) directly, in the x2 term and a xiai term in (cid:126)T . That means there must be x2 i and xiai terms. So there must be an x2 i i and xiai components to (cid:126)Q2. And we can see that there are: Q2x = + [Kxxx(x|x)2 + Kxxa(x|x)(a|x) + Kxaa(a|a)2]x2 + [2Kxxx(x|x)(x|a) + Kxxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kxaa(a|x)(a|a)]xiai + [Kxxx(x|a)2 + Kxxa(x|a)(a|a) + Kxaa(a|a)2]a2 i i 235 and Q2a = + [Kaxx(x|x)2 + Kaxa(x|x)(a|x) + Kaaa(a|a)2]x2 + [2Kaxx(x|x)(x|a) + Kaxa[(x|x)(a|a) + (x|a)(a|x)] + 2Kaaa(a|x)(a|a)]xiai + [Kaxx(x|a)2 + Kaxa(x|a)(a|a) + Kaaa(a|a)2]a2 i i Getting very close. Let’s try to get (a|xx) first. I need the x2 i terms of Q2x and Q2a: Q2x = Kxxx(x|x)2 + Kxxa(x|x)(a|x) + Kxaa(a|a)2 Q2a = Kaxx(x|x)2 + Kaxa(x|x)(a|x) + Kaaa(a|a)2 These need to be evaluated. Q2x = 0 this is because Kxxx = Kxax = Kxaa = 0 cos2((cid:112)3 − n1s/R) Q2a = − 2 R Next step is to evaluate  cos(cid:0)√ 3 − n1s/R(cid:1) 3 − n1s/R(cid:1) sin(cid:0)√ 3−n1 R − R√  3 − n1s/R(cid:1) sin(cid:0)√ 3 − n1s/R(cid:1) cos(cid:0)√ 3−n1   0 Q2a ˆL−1 (cid:126)Q2 = √ Now performing the integrations: (cid:90) (cid:18) (cid:90) T2x = T2a = − R√ sin 3 − n1 (cid:16)(cid:112)3 − n1s/R cos (cid:17)(cid:19) Q2ads (cid:16)(cid:112)3 − n1s/R (cid:17) Q2ads 236 Here is the evaluation of T2x and T2a: (cid:18)1 (cid:18) 3 (cid:19) cos3(cid:16)(cid:112)3 − n1s/R (cid:17) − 1 sin3(cid:16)(cid:112)3 − n1s/R (cid:17) − 1 (cid:16)(cid:112)3 − n1s/R sin 3 (cid:17)(cid:19) 3 T2x = T2a = −2R 3 − n1 −2√ 3 − n1 Here we finally we arrive at the result (a|xx) = The x2 i term of R2a = (a|x)T2x + (a|a)T2a or the x2 √ 3 − n1 R − cos i terms of R2a = (cid:16)(cid:112)3 − n1s/R (cid:17) (cid:16)(cid:112)3 − n1s/R sin (cid:17) 2R (cid:18) 3 − n1 2√ 3 − n1 sin (cid:19) (cid:18)1 cos3(cid:16)(cid:112)3 − n1s/R (cid:17) − 1 (cid:16)(cid:112)3 − n1s/R (cid:17) − 1 sin3(cid:16)(cid:112)3 − n1s/R 3 3 (cid:17)(cid:19) 3 We have already done the work for: (x|xx) = The x2 cos((cid:112)3 − n1s/R) i term of R2x = (x|x)T2x + (x|a)T2a = cos3((cid:112)3 − n1s/R) − 1 (cid:18)1 (cid:19) + R√ 3 − n1 R (−2) 3 − n1 sin(cid:112)3 − n1s/R) 3 1√ 3 − n1 (−2) (cid:18) sin((cid:112)3 − n1s/R − 1 3 sin3((cid:112)3 − n1s/R) (cid:19) 3 237 APPENDIX D Derivation of the Envelope Equation Beams with initially rectangular phase space configurations form piecewise linear envelopes as they move through thin elements. Beams with initially elliptical phase space configurations moving through linear focusing elements retain their elliptical shapes, and form smooth envelopes. To see this, consider a beam with an initially rectangular phase space configuration passing through a drift. Represented in phase space, Fig. D.1 At any point along the path of the beam, the “envelope” is the maximum value of x. So initially the maximum value of x is ±x0. The question is how does the envelope evolve over time? This would be a function of the value of x(cid:48) 0 which each particle happens to have initially. xL = x0 + x(cid:48) 0L. By the time the particle reaches L, the largest value of xL will be reached by the particle posessing the largest values of x0 and x(cid:48) 0. This will be the point at the upper-right corner of the initial distribution. As z grows from 0 to L, this maximum will grow linearly. Hence the envelope with grow linearly from the initial position to the final position. Now how do things change when the initial rectangular configuration is replaced with an initial elliptical configuration? At the initial point, the transverse envelope of the beam is 238 Figure D.1: Rectangular region of phase space. 239 just the semimajor axis of the ellipse. What happens to this ellipse after it passes through a drift? How does the envelope change as a result? For this we need the general equation for an ellipse, and then an equation for the envelope. If a beam has the shape of an ellipse,then the envelope is the edge of that ellipse. As the ellipse evolves in phase space, so too does the envelope. Taking a and b to be the phase space radii, or x2 a2 + x(cid:48)2 b2 = 1, x2 + b a a b x(cid:48)2 = ab = ≡ . area π 240 Generally, we want the ellipse to have an arbitrary rotation by angle θ.  cos θ sin θ − sin θ cos θ  =   x x(cid:48)  x cos θ + x(cid:48) sin θ −x sin θ + x(cid:48) cos θ  Substituting these transformed quantitiies into the original equation for the ellipse gives b a + (x2 cos2 θ + 2xx(cid:48) cos θ sin θ + x(cid:48)2 sin2 θ) a b (x2 sin2 θ + 2xx(cid:48) sin θ cos θ + x(cid:48)2 cos2 θ) = , and collecting like terms gives (cid:18) b (cid:19) cos2 θ + sin2 θ a b a x2 + 2 (cid:19) (cid:18) b a − a b xx(cid:48) cos θ sin θ + (cid:18) b a sin2 θ + cos2 θ a b (cid:19) x(cid:48)2. 241 After making the well-known substitutions (cid:19) = γ sin2 θ a b (cid:18) b (cid:18) b (cid:18) b a cos2 θ + (cid:19) a − a b xx(cid:48) cos θ sin θ = α (cid:19) sin2 θ + cos2 θ a b = β, a we arrive at γx2 + 2αxx(cid:48) + βx(cid:48)2 = . In addition, it can be shown that βγ − α2 = 1 [2], which allows us to eliminate γ from the equation of the ellipse yielding x2 + (xα + x(cid:48)β)2 = β. From Liouville’s theorem [2], the value of  is a constant (being the area of the phase space ellipse divided by π). Therefore it follows that d = 0. Hence or d = ∂ ∂x dx + ∂ ∂x(cid:48) dx(cid:48) = 0 = − dx da ∂ ∂x(cid:48) ∂ ∂x . Taking the partial deriviative of this equation with respect to x(cid:48) gives xα + x(cid:48)β = 0 242 at the envelope; hence Rα + R(cid:48)β = 0. From the original equation this means that or, equivalently, R2 = β R =(cid:112)β As the ellipse is transformed throughout the system, the envelope of the ellipse varies as β. In addition, we have which yields Rα + R(cid:48)β = α(cid:112)β + R(cid:48)β = 0 R(cid:48) = −α(cid:112)β/β = −α(cid:112)/β We can invert these equations to get β = R2/, α = −RR(cid:48)/ Inserting these quantities into the equation of the ellipse, x2 + (x(−RR(cid:48)/) + x(cid:48)R2/)2 = R2/ or 2x2/R2 + (xR(cid:48) − x(cid:48)R)2 = 2, 243 and after differentiating by z and simplifying gives (R(cid:48)x − Rx(cid:48)) R(cid:48)(cid:48)x − Rx(cid:48) + 2 R3 (cid:20) (cid:21) = 0. The prefactor is never zero, so we must have Finally, we have so R(cid:48)(cid:48)x − Rx(cid:48)(cid:48) − 2 R3 x = 0. x(cid:48)(cid:48) + k2(z)x = 0, R(cid:48)(cid:48) + R(cid:48)k2(z) − 2 R3 = 0. which is called “the envelope equation”. 244 APPENDIX E Paraxial Equation Laminar beam Derivation of the paraxial equation. Following Reiser [44], the potential of an axially sym- metric charge distribution is a solution of the Laplace equation. In cylindrical coordinates this is 1 r ∂ ∂r We attempt a solution of the form (cid:18) (cid:19) r ∂φ ∂r + ∂2φ ∂z2 = 0. ∞(cid:88) ν=0 φ2ν(z)r2ν (E.1) φ(r, z) = where we only need even powers due to the requirement that the field on the axis r = 0 must vanish. Substituting this proposed solution into Eq. E.1 yields Shifting the indices in the first term yields ∞(cid:88) 4ν2φ2ν(z)r2ν−1 + ν=1 ν=0 ∞(cid:88) ∞(cid:88) 4(ν + 1)2φ2(ν+1)(z)r2ν + ν=0 245 φ(cid:48)(cid:48) 2ν(z)r2ν = 0. ∞(cid:88) ν=0 φ(cid:48)(cid:48) 2ν(z)r2ν = 0, or (cid:104) 4(ν + 1)2φ2ν+2 + φ(cid:48)(cid:48) 2ν (cid:105) ∞(cid:88) ν=0 r2ν = 0. This is a recursion relation that allows us to write or 4φ2(z) + φ(cid:48)(cid:48) 0 = 0, φ2(z) = −1 4 φ(cid:48)(cid:48) 0. This is an important result. It means that higher order terms in this expansion can be calculated solely from the value of the potential on the axis φ0. This on-axis potential is often known exactly. We have to second order: φ(r, z) ≈ φ(0, z) − 1 4 ∂2φ(0, z) ∂z2 r2 = φ0 − 1 4 φ(cid:48)(cid:48) 0r2, and Er ≈ −∂φ ∂r ≈ φ(cid:48)(cid:48) 0r 2 . The equation of motion can therefore be approximated by m¨r = qEr = φ(cid:48)(cid:48) 0. qr 2 Transforming from t to z as the independent variable gives r(cid:48)(cid:48)v2 z + r(cid:48)vzv(cid:48) z = φ(cid:48)(cid:48), qr 2m 246 or r(cid:48)(cid:48) + r(cid:48) vzv(cid:48) z v2 z = φ(cid:48)(cid:48). qr 2mv2 z (E.2) If we make the assumption that vz >> vr, then kinetic energy ≈ 1 2 mv2 z . If we make the further simplifying assumption that kinetic energy + potential energy = 0, we can write mv2 z = −qφ0. 1 2 (E.3) Differentiating both sides with respect to z gives mvzv(cid:48) z = −qφ(cid:48) 0 Eq. E.3 and Eq. E.4 can be used to simplify Eq. E.2: r(cid:48)(cid:48) + φ(cid:48) 0 2φ0 r(cid:48) + φ(cid:48)(cid:48) 0 4φ0 = 0. which is the nonrelativistic envelope equation. 247 BIBLIOGRAPHY 248 BIBLIOGRAPHY [1] Cronin, Jane. Differential Equations: Introduction and Qualitative Theory. New York: M. Dekker, 1994. [2] Wiedemann, Helmut. Particle Accelerator Physics, 4th ed. Cham: Springer International Publishing, 2015. [3] Berz, Martin, Kyoko Makino, and Khodr Shamseddine. Advances in Imaging and Elec- tron Physics, 108: Modern Map Methods in Particle Beam Physics. Burlington: Elsevier, 1999. [4] Menzel, M.T, and H.K Stokes. User’s Guide for the Poisson. Washington, D.C: United States. Dept. of Energy, 1987. [5] K. L. Brown. A First and Second Order Matrix Theory for the Design of Beam Transport Systems and Charged Particle Spectrometers. Adv. Part. Phys. 1, 71 (1968). [6] Andreas Lehrach, Bernd Lorentz, William Morse, Nikolai Nikolaev, Frank Rathmann, arxiv:1201.5773v1. 27 Jan 2012. [7] Andreas Lehrach, private communication. [8] Hermann Wollnik, Optics of Charged Particles. Academic Press, 1987. [9] Hans Grote, F. Christoph Iselin, The MAD Program (Methodical Accelerator Design), Version 8.13/8, User’s Reference Manual. CERN/SL/90-13 (AP) (Rev. 4). [10] M´eot, S. Valero, ZGOUBI Users’ Guide - version 4.3. CEA DSM DAPNIA-02-395 (2013). [11] M. Berz, K. Makino, COSY INFINITY 9.1 Beam Physics Manual. MSUHEP 060804- rev (2013). [12] Enge, Deflecting Magnets in Focusing of Charged Particles. ed. A. Septier, Academic Press (1967). [13] M.L. Shashikant, M. Berz, B. Erdelyi, COSY Infinity’s EXPO Symplectic Tracking for LHC. Institute of Physics CS 175 (2004) 299-306. [14] FAIR Baseline Technical Report. GSI (2006). 249 [15] Muon g-2 Technical Design Report. 2015, arXiv:1501.06858. [16] JEDI Collaboration, http://collaborations.fz-juelich.de/ikp/jedi. [17] R. Hipple, M. Berz, Microscopy and Microanalysis 21 Suppl. 4 (2015). [18] R. Hipple and M. Berz,Implementation of Benchmarks for Precision Tracking in Storage Rings, in Proceedings of ICAP 2015, Shanghai, China, (2015) MODBC3. [19] M. Berz, K. Makino, W. Wan, An Introduction to Beam Physics. CRC Press (2015). [20] M. Berz, K. Makino, COSY INFINITY 9.2 Beam Physics Manual, MSUHEP-151103 (2015). [21] U. Bechstedt, et al., Nucl. Instrum. Methods Phys. Res. B 113 (1996) 26-29. [22] B. Lorentz, HESR Linear Lattice Design in Proceedings of EPAC 2008, Genoa, (2008) MOPC112. [23] E.M. Metodiev et al., Nucl. Instrum. Methods Phys. Res., Sect. A 797, 311 (2015). [24] F. M´eot, S. Valero, ZGOUBI Users’ Guide, CEA DSM DAPNIA-02-395 (2013). [25] H. Weick, GICOSY User’s Reference Manual, CERN/SL/90-13 (AP) (Rev. 4). [26] H. Enge in Focusing of Charged Particles, A Septier, Ed. New York, NY, USA: Academic Press, 1967, pp. 203-264. [27] H. Grote, et al., The MAD-X Program (Methodical Accelerator Design) Version 5.03.00, User’s Reference Manual http://mad.web.cern.ch/mad/. [28] H. Grote, F. Christoph Iselin, The MAD Program (Methodical Accelerator Design), Version 8 http://mad.web.cern.ch/mad/. [29] H. Wollnik, M. Berz, Relations Between Elements of Transfer Matrices Due to the Condition of Symplecticity. Nuclear Instruments and Methods in Physics Research A, 238 (1985). [30] R. Talman, Representation of Thick Quadrupoles by Thin Lenses, SSC-N-33, (1985). [31] Andreas Lehrach, Searching Electric Dipole Moments in Storage Rings. Proceedings of ICAP 2015, Shanghai, China, (2015). [32] J. Grange, Muon g-2 Technical Design Report, (2015), arXiv:1501.06858. 250 [33] U. Bechstedt, Nucl. Instrum. Methods Phys. Res. B 113 (1996) 26-29. [34] M. Berz, K. Makino, COSY INFINITY 9.2 Beam Physics Manual, MSUHEP-151103 (2015). [35] E.M. Metodiev et al., Nucl. Instrum. Methods Phys. Res., Sect. A 797, 311 (2015). [36] Hermann Wollnik, History of Mass Measurements in Time-of-Flight Mass Analyzers, International Journal of Mass Spectrometry, 349-350 (2013) 38-46. [37] A.E. Cameron, D.F. Eggers, An Ion Velocitron, Review of Scientific Instruments 19, (1948) 605. [38] B. A. Mamyrin, The Mass-Reflectron, A New Nonmagnetic Time-of-Flight Mass Spec- trometer With High Resolution, Zh. Eksp. Teor. Fiz. 64, (1973) 82-89. [39] R. Kutscher, R. Grix, G. Li and H. Wollnik, A Transversally and Longitudinally Focusing Time-of-Flight Mass Spectrometer, International Journal of Mass Spectrometry and Ion Processes, 103 (1991) 117-128. [40] R. Wiley, Time-of-Flight Mass Spectrometer with Improved Resolution, Review of Sci- entific Instruments 26, 1150 (1955). [41] A. Piechaczek, MTOF- A High Resolution Isobar Separator for studies of Exotic De- cays in Proceedings of the Fourth International Conference of Fission and Properties of Neutron-Rich Nuclei, Sanibel Island, USA, (2007). [42] Robert J. Cotter, Time-of-Flight Mass Spectrometry - Instrumentation and Applications in Biological Research, American Chemical Society, Chapter 2 (1997) [43] Stanley Humphries, Charged Particle Beams, John Wiley and Sons, (1990). [44] Martin Reiser, Theory and Design of Charged Particle Beams, 2nd Edition, Wiley-VCH, (2008). [45] S. Kreim, Recent Exploits of the ISOLTRAP Mass Spectrometer, Nuclear Instruments and Methods in Physics Research B 317 (2013) 492–500. [46] R.N. Wolf, On-line Separation of Short-Lived Nuclei by a Multi-Reflection Time-of- Flight Device, Nuclear Instruments and Methods in Physics Research A 686 (2012) 82–90. [47] R.N. Wolf, M. Eritt, G. Marx, L. Schweikhard, Hyperfine Interactions 199 (2011) 115. [48] F. Wienholtz, International Journal of Mass Spectrometry 421 (2017) 285–293. 251 [49] F. Wienholtz, Masses of Exotic Calcium Isotopes Pin Down Nuclear Forces, Nature, Volume 498, (2013) 346. [50] R.N. Wolf, ISOLTRAP’s Multi-Reflection Time-of-Flight Mass Separator/Spectrometer, International Journal of Mass Spectrometry, 349–350 (2013) 123–133. [51] Yavor, Mikhail, Optics of Charged Particle Analyzers. Amsterdam: Academic Press, (2009). [52] http://sci.esa.int/rosetta/. [53] Nuclear Instruments and Methods in Physics Research B 266 (2008) 4510–4514. [54] http://warp.lbl.gov. [55] Dr. Ryan Ringle, Facility for Rare Isotope Beams, private communication. [56] https://en.wikipedia.org/wiki/On-Line Isotope Mass Separator. [57] R. Hipple, M. Berz, Comparison of Tracking Codes for the Determination of Dynamic Aperture in Storage Rings, in Proc. 7th Int. Particle Accelerator Conf. (IPAC’16), Busan, Korea, May 2016, paper WEPOY053, pp. 3114–3116. [58] Glaser W. Elektronen- und Ionenoptik. In: Fl¨ugge S. (eds) Optics of Corpuscles / Ko- rpuskularoptik. Encyclopedia of Physics / Handbuch der Physik, vol 6. Springer, Berlin, Heidelberg (1956). [59] E. Valetov, M. Berz, MSUHEP-180212, Michigan State University (2018). [60] Jackson, J D. Classical Electrodynamics, 2nd ed.. New York: J. Wiley, (1975). [61] Gross, J¨urgen, Mass Spectrometry, Third Edition, Springer (2013). [62] Aston, F.W., The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 38:228, (1920) 707-714. [63] Pedersen, Strasser, Amarant Physical Review A, Volume 65 (2002) 042704. [64] D Strasser, Geyer, Pedersen, Phys. Rev. Lett. 89 (2002) 28. [65] Ken Carter, private communication. [66] Peter Ostroumov, Facility for Rare Isotope Beams, private communication. [67] Manura, David, and David A. Dahl. Simion 8.0 User Manual. Ringoes, NJ: Scientific Instrument Services, (2007). 252 [68] http://lise.nscl.msu.edu/lise.html [69] Hermann Wollnik, Second Order Approximation of the Three Dimensional Trajectories of Charged Particles in Deflecting Electrostatic and Magnetic Fields Nuclear Instruments and methods 34 (1965) 213-221 [70] Lanczos, C., The Variational Principles of Mechanics, 7th ed.Dover Publications(2012) 253