THEORY AND MODELING OF INTENSE ION BEAMS AND DIAGNOSTIC MEASUREMENTS IN ACCELERATOR FRONT ENDS By Chun Yan Jonathan Wong A DISSERTATION Michigan State University in partial fulfillment of the requirements Submitted to for the degree of Physics – Doctor of Philosophy 2019 ABSTRACT THEORY AND MODELING OF INTENSE ION BEAMS AND DIAGNOSTIC MEASUREMENTS IN ACCELERATOR FRONT ENDS By Chun Yan Jonathan Wong Beam dynamics in front ends of ion accelerators is often subject to strong space charge and significant transverse coupling due to magnetic fields in the source and the use of focusing solenoids in the beam line. This thesis starts with an overview and a brief primer on beam dynamics to frame two research themes: simulation and analytic studies on beam dynamics, and improvements on transverse beam diagnostics. Each theme is first developed in a general treatment, followed by applications to the front end of the Facility for Rare Isotope Beams (FRIB). The first research theme includes: (1) numerical modelling of accelerator front ends, and (2) analytic studies of coupled transverse beam dynamics. In (1), a set of readily adaptable simulation tools based on the PIC code Warp is constructed to model the transport and manipulation of intense multi-species beams with high levels of detail. The tools are applied to investigate beam evolution in the FRIB front end including the species separation process. In (2), general analytic results on beams with rotational symmetry and the formalism of eigen-emittances are derived to elucidate non-axisymmetric initial conditions at the ion source and beam states advantageous for solenoid transport. The results are employed to construct LEBT/MEBT (low/medium energy beam transport) tuning schemes currently adopted at FRIB for transverse matching of strongly magnetized beams into the RFQ and cryomodules respectively. The second research theme presents improvements on transverse beam diagnostics including: (1) corrections in the conventional data analysis model of high-resolution phase space measurements with Allison scanners, and (2) error minimization in beam matrix measurements and enhanced tomography capabilities using beam profile monitors. All techniques and algorithms under this theme have been verified using measurements at FRIB. To my family iii ACKNOWLEDGEMENTS First and foremost, I would like to express my heartfelt gratitude to my advisor, Professor Steven Lund, for his unfailing guidance and support. He granted me the freedom to pursue a variety of different projects, and always made time with a single knock on his door to share his wisdom in physics as well as in life. As a scientist, his excellence and integrity set the standards I aspire to meet, and his enthusiasm is a constant source of encouragement. As a mentor and teacher, his patience, his vision for education, and his strong commitment to preparing all students to succeed, continue to inspire me to contribute to the accelerator physics community. I was very fortunate to have the joy of learning from him. It has been my great privilege to work with the Accelerator Physics (AP) Department at FRIB. I thank Professor Peter Ostroumov for maintaining an exciting research environment in which I was nurtured, and for giving me invaluable opportunities to participate in beam studies. His remarkable intuition and astute judgment played an important role in shaping the course of my research. I am grateful to Professor Yue Hao who always expressed genuine interest in my work and well-being. I greatly benefited from his immense knowledge and career advice. I am also indebted to many congenial collaborators in AP: Kei Fukushima for many fruitful discussions on beam simulations; Tomofumi Maruta for his physical insight in beam tuning and his generous supervision of my experimental studies; Alexander Plastun for his prowess in beam optics and his perceptive observation of the devils lying in the details; Takashi Yoshimoto for his intriguing work; Tong Zhang for the joint efforts in application development that brought my understanding of accelerator control system to a whole new level; and Qiang Zhao for sharing his practical wisdom and extensive experience. They taught me a lot and their work constantly fueled my interest and curiosity. I would like to take this opportunity to wish AP and the FRIB project resounding success in the years to come. I owe many thanks to Professor Steven Lund (chair), Professor Philip Duxbury, Professor Yue Hao, Professor Peter Ostroumov, Eduard Pozdeyev and Professor Vladimir Zelevinsky for serving iv on my thesis committee and providing invaluable feedback. Professor Yoshishige Yamazaki also shared thoughtful comments on countless occasions. I gratefully acknowledge Professors Martin Berz and Kyoko Makino for their courses and textbooks on beam physics. The unique techniques and approaches I learned have been very beneficial to my understanding and research work. I would like to thank Scott Cogan, Diego Omitto and Rebecca Shane for their technical support on Allison scanner measurements at FRIB; Guillaume Machicoane and Jeff Stetson for their insight on ECR ion sources and the associated beam transport; and Haitao Ren for his support on modeling of the front end. Special thanks must be given to Steven Lidia for many enlightening conversations on beam diagnostics that provided much-needed reality checks and refocused my work on the important problems. My research updates were regularly presented at a meeting that mainly consists of MSU graduate students in accelerator physics. I thank Daniel Alt, Michael Balcewicz, Crispin Contreras, Robert Hipple, Bryan Isherwood, Andrew LaJoie, Didi Luo, Derek Neben, Christopher Richard, and Safwan Shanab for their questions and feedback that helped me understand my work. I am also grateful to Harry He, Felix Marti and Alfonse Pham who were only too willing to provide their astute comments. I appreciate the assistance from Kim Crosslan at the Department of Physics and Astronomy that allowed me to meander through all the paperwork in graduate school. I am deeply indebted to the US Particle Accelerator School for the invaluable learning oppor- tunities and social interactions that have been immensely beneficial to my thesis research and my growth in the field of accelerator science. I am grateful to the staff Susan Winchester and Irina Novitski for making every session such a memorable and addicting experience. Last but not least, I owe a special debt to Professor Simon Yu, my master’s thesis advisor at the Chinese University of Hong Kong, for guiding me into the field of accelerator physics and laying the foundation in both knowledge and research mentality. v TABLE OF CONTENTS . . . . . . . . . . . . LIST OF TABLES . LIST OF FIGURES . CHAPTER 1 . . . . INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1 High Intensity Ion Accelerators . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Facility for Rare Isotope Beams (FRIB) . . . . . . . . . . . . . . . . . . . . . . . 1.3 Electron Cyclotron Resonance (ECR) Ion Sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Beam Dynamics Challenges in Accelerator Front Ends ix x 1 1 2 4 5 CHAPTER 2 A PRIMER ON BEAM PHYSICS . . . . . . . . . . . . . . . . . . . . . . 2.1 Particle Motion and Phase Space . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Symplectic Dynamics, Beam Matrices and Emittance Conservation 7 7 . . . . . . . . 10 Symplectic Condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2.2 Beam Matrices 2.2.3 Emittance Conservation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.3 2D Linear Phase Space Motion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.4 Axisymmetry and Conservation of Canonical Angular Momentum . . . . . . . . . 14 . . . . . . . . 3.1 Warp Code . 3.2 Beam Line Element Models 3.3 Simulations of FRIB Front End CHAPTER 3 SELF-CONSISTENT SIMULATIONS . . . . . . . . . . . . . . . . . . . . 16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 . . . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.1 Initial Conditions & Loading . . . . . . . . . . . . . . . . . . . . . . . . . 18 3.3.2 Lattice Design & Matching Methodology . . . . . . . . . . . . . . . . . . 19 3.3.3 Envelope Model for Multi-Species Beam Transport . . . . . . . . . . . . . 23 Simulation Cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3.4 3.3.5 Warp Simulation Results: Case Comparisons . . . . . . . . . . . . . . . . 26 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30 3.4 Further Work . . . . . . CHAPTER 4 THEORY OF TRANSVERSE COUPLED BEAM DYNAMICS, WITH 4.1 Beams with Rotational Symmetry . . . . . . . . . . . . . . . . . . . . . . . . . APPLICATIONS TO THE FRIB FRONT END . . . . . . . . . . . . . . . 34 . 34 4.1.1 Effect of Symmetry on Beam Moments . . . . . . . . . . . . . . . . . . . 34 4.1.2 Derivation of Moment Relations from Symmetry . . . . . . . . . . . . . . 37 4.1.3 Correspondence between Discrete and Continuous Rotational Symmetry . 39 . . . . . . . . . . . . . . . . . . . . 42 4.1.4 Example: Beams with C3 Symmetry . 42 . 44 4.2 Coupled Beam Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2.1 Axisymmetric Beams and Eigen-emittances . . . . . . . . . . . . . . . . . 46 4.2.2 Magnetized Beam . . 47 2nd Order Moments . . . . . . . . . . . . . . . . . . . . . . . 3rd Order Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.4.1 4.1.4.2 vi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2.3 Mirror-Symmetric Beam Line Does Not Preserve Axial-symmetry . . . . . 49 . 51 4.2.4 Linear Solenoid Transport 4.2.5 Quasi-axisymmetric Conditions . 52 4.3 Applications to FRIB Front End . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.3.1 LEBT: Before Species Selection . . . . . . . . . . . . . . . . . . . . . . . 55 4.3.2 LEBT: After Species Selection . . . . . . . . . . . . . . . . . . . . . . . . 56 4.3.3 MEBT and LS1 . . 61 4.3.3.1 Optics Design and Matching Recipe . . . . . . . . . . . . . . . . 62 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3.4 . . CHAPTER 5 TRANSVERSE PHASE-SPACE MEASUREMENTS WITH ALLISON- . . . . . . . 5.2.2 Introduction . TYPE EMITTANCE SCANNERS . . . . . . . . . . . . . . . . . . . . . . 65 5.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 5.2 Additional Geometric Features . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 5.2.1 Geometric Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 Longitudinal Asymmetry in E-dipole Placement (l1 (cid:44) l2) . . . . . 68 5.2.1.1 Finite Slit Thickness d (cid:38) Slit Width s . . . . . . . . . . . . . . . 68 5.2.1.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Four Cases . Particle Simulations . . . . . . . . . . . . . . . . . . . . . . . . 70 5.2.2.1 Particle Transmission Quantified . . . . . . . . . . . . . . . . . . 72 5.2.2.2 5.2.2.3 Analytic Results . . . . . . . . . . . . . . . . . . . . . . . . . . 75 Two Corrections in Data Analysis . . . . . . . . . . . . . . . . . . . . . . 78 5.2.3.1 Voltage-to-Angle Relation . . . . . . . . . . . . . . . . . . . . . 78 5.2.3.2 Voltage-Dependent Weight Compensation . . . . . . . . . . . . . 78 5.2.4 Analysis of Synthetic Data . . . . . . . . . . . . . . . . . . . . . . . . . . 82 5.3 Corrections for Degenerate Phase Space Measurements . . . . . . . . . . . . . . . 83 Example: FRIB Allison Scanner . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.4 5.4.1 FRIB Front End . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5.4.2 Quadrupole Scan Measurements . . . . . . . . . . . . . . . . . . . . . . . 90 Total Current Normalization Anomaly . . . . . . . . . . . . . . . . . . . . 93 5.4.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 5.5 Conclusion . . 5.2.3 . . . . . . . . CHAPTER 6 OPTIMIZED TRANSVERSE PHASE-SPACE MEASUREMENTS WITH 6.1 Review of Linear Algebra . BEAM PROFILE MONITORS . . . . . . . . . . . . . . . . . . . . . . . . 97 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97 . . . . . . . . . . . . . . . . . . . . . . . . 98 6.1.1 Least Squares Approximation 6.1.2 Singular Value Decomposition . . . . . . . . . . . . . . . . . . . . . . . . 99 6.1.3 Vector Norms and Matrix Norms . . . . . . . . . . . . . . . . . . . . . . . 99 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 6.1.4 Condition Number . 104 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 1st Order Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 2nd Order Moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 . . . . . . . . . . . . . . . . . . . . . . . 107 . 108 6.2.1 Notation . 6.2.2 6.2.3 6.2.4 A Limitation of Solenoid Scans 6.2 Quadrupole Scan and Solenoid Scan: Working Principles . . . . . . . . . . . . . 6.3 Quadrupole Scan and Solenoid Scan: Error Quantification and Minimization . . . . . . . vii 6.3.1 Error Quantification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3.2 Error Minimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 . 108 6.4 Spatial and Phase Space Tomography . . . . . . . . . . . . . . . . . . . . . . . . . 111 . 111 . 114 6.5 Example: FRIB Measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.6 Further Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120 Spatial Tomography . . . . . . . . . . . . . . . . . . . . . . . . . . . . Phase Space Tomography . . . . . . . . . . . . . . . . . . . . . . . . . . 6.4.1 6.4.2 . . . . . . . . . . . . . . . . APPENDICES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122 APPENDIX A EXTENSIONS ON ALLISON SCANNER ANALYSIS . . . . . . . . 123 APPENDIX B ALTERNATIVE PROOF OF RELATIONS GOVERNING 2ND ORDER MOMENTS OF CN BEAMS . . . . . . . . . . . . . . . . . 136 APPENDIX C GENERATING A 4D COUPLED DISTRIBUTION . . . . . . . . . . 144 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 BIBLIOGRAPHY . . . . . . . . . viii LIST OF TABLES Table 3.1: Beam Parameters of U34+ based on the VENUS-like source. For definitions, see also Sec. 3.3.3 and ??. Temperatures and envelope parameters are the same for all ion species (indices suppressed). Canonical angular momentum, radial normalized emittance, and kinetic energy/rigidity vary with ion species. . . 21 Table 4.1: Beam states entailed by the requirement of 2nd-order axisymmetry in the lower LEBT. (cid:88) and X denote true and false respectively. . . . . . . . . . . . . . . . . 60 Table 4.2: First set of quadrupole scan parameters . . . . . . . . . . . . . . . . . . . . . . 61 Table 5.1: Geometric parameters applied in examples corresponding to the four geometric models in Fig. 5.2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 Table 5.2: Analytic formulae corresponding to the geometric models in Fig. 5.2 for V0 ≥ 0. Results for V0 < 0 obey the same formulae with V0 → |V0| and the x-axis reversed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 . . . . . . . Table 5.3: Beam parameters of the three distributions in Fig. 5.8. . . . . . . . . . . . . . . 84 Table 5.4: Beam parameters of the two distribution in Fig. 5.10 compared against those of the waterbag distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 Table 5.5: Allison scanner geometries at FRIB front end, before and after modification in Feb 2018. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 Table 5.6: Beam parameters of corresponding to the phase space distribution at VQ3 = . . . . . . . 1000 V in Fig. 5.13, computed with and without correction terms. . 92 Table 5.7: Solution of Eq. (5.11) for the quadrupole scan in Sec. 5.4.2, with measurement results x f calculated using conventional and improved analysis method respectively. 93 Table 6.1: Quadrupole parameters for Scan 1 . . . . . . . . . . . . . . . . . . . . . . . . . 116 Table 6.2: Quadrupole parameters for Scan 2 . . . . . . . . . . . . . . . . . . . . . . . . . 117 Table 6.3: Beam moments and emittances corresponding to the phase space distributions plotted in Fig. 6.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Table 6.4: Comparison between beam moments of the reconstructed spatial distribution in Fig. 6.11 and the values obtained from wire scanner measurements. . . . . . . 119 ix LIST OF FIGURES Figure 1.1: High intensity ion accelerators (Image courtesy of Prof. Jie Wei at FRIB) [1] . . Figure 1.2: Schematic of FRIB superconducting driver linac. Note that part of the front end is located on ground level, a complete illustration of the front end is shown in Fig. 1.3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.3: Schematic of the FRIB front end with emphasis on the location of diagnostics devices. Abbreviations are listed below: BCM - beam current monitor; BPM - beam position monitor; FC - Faraday cup; MHB - multi-harmonic buncher; PM - profile monitor. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 1.4: Schematic of an electron cyclotron resonance (ECR) ion source (Image cour- tesy of Dr. Daniele Leitner at LBNL). . . . . . . . . . . . . . . . . . . . . . . Figure 1.5: Measured charge state distribution of a uranium beam from the VENUS ion source [9]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 2.1: Frenet-Serret coordinate system. Image adapted from Ref. [11] Figure 2.2: Ellipse in x-x(cid:48) phase space [17]. Figure 3.1: CST Studio model of S4 solenoids at FRIB and the calculated fields at different . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2 3 4 5 6 8 radii within the clear bore aperture radius R = 164.2/2 mm. . . . . . . . . . . 19 Figure 3.2: CST Studio model of a Q7 electrostatic quadrupole triplet at FRIB with dimensions and potential plots. . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Figure 3.3: CST Studio model of a Q7 electrostatic quadrupole triplet at FRIB. Dimen- sions and potential contours at various axial locations are shown. . . . . . . . . 21 Figure 3.4: CSS lattice design used in Warp simulations. . . . . . . . . . . . . . . . . . . . 22 Figure 3.5: Comparisons between numerical solutions of the envelope equation and Warp simulation results in the transport line upstream of the CSS. These correspond to simulation cases 3 and 4 described in Sec. 3.3.4. . . . . . . . . . . . . . . . 24 Figure 3.6: Accessible lattice functions β and α at the entrance of the first dipole when the Bpeak in both solenoids varied independently between 0 T and 1.5 T. Results are shown for different space charge neutralization factors. Observe that there are regions of initial Twiss parameters that are inaccessible, and the region of accessible Twiss parameters decreases with decreasing neutralization factor . . 25 x Figure 3.7: Warp simulation results of Uranium species separation and collimation in the beam evolution from the ECR source extraction through the midpoint of the charge selection system showing collimation to two target species for high intensity FRIB operation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 Figure 3.8: Rms envelops σx,y a) and normalized rms-emittance εx,y b) of the reference species U34+ in the CSS. c) shows the centroid of U33+ and x y-particle projections of both target species at the CSS mid-point for two cases. Solid curves are x-plane and dashed curves are y-plane. Colors flag Warp sim- ulation cases. Black solid lines denote MADX design values in a) and c). The axial midplane of the CSS is indicated in c) (z in FRIB coordinates with origin outside machine). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 Figure 3.9: Extraction geometry of the VENUS-like ECR ion source. The electrode in red has negative voltage to prevent backflow of electrons. (Image courtesy of Dr. Daniel Winklener.) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 3.10: Beam sizes in the post-extraction beam line downstream of the VENUS-like ECR at FRIB, obtained from Warp simulations with different homogeneous neutralization factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 Figure 4.1: Simulated spatial distribution of Ar9+ ions at the extraction plane of the Artemis [5] ECR ion source at FRIB. The inner circle indicates the extraction aperture. x-y and u-v denote two possible transverse coordinate systems. Subsequent analysis will show that, despite severe triangulation, the rms envelope and emittance of the beam is the same along any transverse direction. Simulation results courtesy of the code developed by Dr. Vladirmir Mironov at JINR-Dubna [25, 26]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Figure 4.2: A schematic of transverse focusing and bending elements in the FRIB front end. Note that Segment 2 and 5 are mirror-symmetric about the axial mid-plane. Lengths are not drawn to scale and all quadrupoles are normal quadrupoles (i.e. upright orientation). . . . . . . . . . . . . . . . . . . . . . . 55 Figure 4.3: An example of FRIB LEBT tuning results where the focusing strengths of 17 elements (10 quadrupoles and 7 solenoids) were optimized to achieve desired beam conditions at the RFQ entrance. (Image courtesy of Dr. Tomofumi Maruta at FRIB Accelerator Physics Department) . . . . . . . . . . . . . . . . 58 Figure 4.4: An example of FRIB LEBT tuning results where the tuning scheme introduced in Sec. 4.3.2 was applied to optimize elements strengths in three stages. (Image courtesy of Dr. Tomofumi Maruta at FRIB Accelerator Physics Department) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 . . . . xi Figure 5.1: Side-view of the FRIB Allison scanner. . . . . . . . . . . . . . . . . . . . . . . 69 Figure 5.2: Four geometric cases and their respective characteristic trajectories. Case I: l1 = l2 = l and d = 0; Case II: l1 (cid:44) l2 and d = 0; Case III: l1 = l2 = l and d (cid:44) 0; Case IV: l1 (cid:44) l2 and d (cid:44) 0. zi is the longitudinal position where the beam particles enter the device and where the x-x(cid:48) phase space is measured. Figure 5.3: Blue dots denote particles that are transmitted in the simulations described in Sec. 5.2.2.1. They map x-x(cid:48) phase space regions at the device entrance z = zi that contribute to the data point at x = 0 mm. Results are shown for all four geometric cases in Fig. 5.2 for (a) V0 = 500V and (b) V0 = 1000V respectively. Geometric parameters for each case are listed in Table 5.1. Black lines define the phase space grid with increments determined by the scanning step sizes. Dashed blue lines indicate the nominal position x = 0 and angle x(cid:48) = x(cid:48) ref corresponding to the data point. The red rectangles denote . . the device resolution predicted by the conventional treatment (i.e. Case I). . . 71 . 73 . . . . . Figure 5.4: Illustrates how the shrunk phase space area attains its shape in Case III of Fig. 5.3b. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 Figure 5.5: Angular transmission factor T(x(cid:48)) for Cases I to IV in Fig. 5.3b. . . . . . . . . . 75 Figure 5.6: Voltage-to-angle relation for the selected beam. The red curve (analytic) overlays the blue curve (ballistic simulation). . . . . . . . . . . . . . . . . . . . 79 Figure 5.7: Plots of a) angular resolution ∆x(cid:48), and b) angular correction factor as a ref. The red curve (analytic) overlays the blue curve (ballistic . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 function of x(cid:48) simulation). . . Figure 5.8: Uniform phase space projection of the KV distribution, compared against results of synthesized measurement data analyzed using conventional and improved analysis methods. The left plot shows a sampling (1 out of 1000 particles) used to represent the distribution. The middle and right plots are fine mesh density plots with the same color scale. . . . . . . . . . . . . . . . . 83 Figure 5.9: Phase space region contributing to the data point at x = 0 mm and V0 = 200 V (presentation format as detailed in Fig.5.3), where the Allison scanner geometry corresponds to Case IV of Table 5.1 except for increased slit width to s = 300 µm. Scanning step sizes in x and V0 are 0.5 mm and 20 V respectively. 84 xii Figure 5.10: Waterbag distribution (left column), compared against synthesized measure- ment data analyzed without (middle column) and with (right column) correc- tions for degenerate measurements. The upper row corresponds to intensity, whereas the lower row shows all non-zero intensity points in the same color to highlight extent errors. Black and white ellipses denote the ideal boundary of the waterbag distribution analyzed. In the left column sampled particles (1 out of 1000) are plotted. Measurement data are synthesized with s = 300 µm, and scanning step sizes of 0.5 mm and 10 V. . . . . . . . . . . . . . . . . . . . 86 Figure 5.11: Fraction of particles that fall outside the ellipse (orientation defined by Twiss parameters) with area πε in each of the three distributions in Fig. 5.10. The vertical dashed line corresponds to the sharp phase space edge of the waterbag distribution at ε/εrms = 6. The blue curve approaches ε/εrms = 6 if the lower range of the ordinate is extended. . . . . . . . . . . . . . . . . . . . . . . . . . 87 Figure 5.12: Upstream end of the low energy beam transport (LEBT) line of the FRIB front end. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 Figure 5.13: y-y(cid:48) phase-space projections measured at the FRIB front end during a quadrupole scan. VQ3 denotes the voltage at the last quadrupole upstream of the Allison scanner. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 . . . . . . . Figure 5.14: Projected beam currents from Allison scanner measurements in a quadrupole scan, compared against the beam current measured at the Faraday Cup slightly upstream (see Fig. 5.12). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 Figure 5.15: Screenshot of measured currents on the x-direction Allison scanner on the Artemis beam line of FRIB, plotted against the bias voltage of the Faraday cup on the scanner. The scanner position and voltage were fixed at a phase space coordinate with high density throughout the measurements. . . . . . . . . 95 Figure 6.1: Schematic of the beam line section containing the first profile monitor at the FRIB Front End. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 Figure 6.2: Schematic of a wire scanner at the FRIB front end. Image courtesy of Tomofumi Maruta at FRIB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 Figure 6.3: A diagonal wire only gives the total current in each diagonal because and each grid cell along a diagonal always contributes equally to a ray. The alternating diagonals are plotted using chessboard coloring on the right. . . . . . . . . . . 113 Figure 6.4: A wire intersecting a 2 : 1 grid at 45o is equivalent to the case of a wire intersecting an even grid at 63o. . . . . . . . . . . . . . . . . . . . . . . . . . . 113 xiii Figure 6.5: For a rotated grid, the wires are no longer vertical, horizontal or diagonal with respect to it. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 Figure 6.6: Spatial tomography of several synthesized distributions. . . . . . . . . . . . . . 115 Figure 6.7: Schematic of the beam line section containing the first profile monitor at the FRIB Front End. The profile monitor is located at the position designated by f and the beam was reconstructed at position i. Green blocks denote identical (Q7 type) electrostatic quadrupoles. . . . . . . . . . . . . . . . . . . . . . . . . 116 Figure 6.8: Projection angles on the initial phase space corresponding to quadrupole scan parameters in Table 6.1 and Table 6.2. . . . . . . . . . . . . . . . . . . . . . . 117 Figure 6.9: Histograms of the normalized x-emittance with random errors applied to the measurement results. The applied errors have a truncated Gaussian distribu- tion where 3σ = 1% for Scan 1 and 3σ = 10% for Scan 2. . . . . . . . . . . . . 118 Figure 6.10: Comparison between tomographic reconstruction of the x-x(cid:48) phase space distribution, and the measured distribution at the Allison scanner propagated to the same location. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118 Figure 6.11: Spatial tomography with a grid rotated by 20o clockwise with respect to the standard x-y coordinate system. . . . . . . . . . . . . . . . . . . . . . . . . . . 119 Figure A.1: Potential contours of the FRIB Allison scanner in Fig. 5.1, including details of the fringe structure in the vicinity of the entrance slit (right). . . . . . . . . . 125 Figure A.2: How (a) shifting a trajectory for fixed angle and (b) changing the angle for fixed initial position allow one to calculate quantities listed in Table 5.2. . . . . 126 Figure A.3: Reversing particle velocity allows one to check results when the E-dipole placement is asymmetric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 Figure A.4: Comparison between transmission ratio in Case III and IV when weight correction effect is large. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 Figure A.5: Beam region is defined by the red ellipse. . . . . . . . . . . . . . . . . . . . . . 130 Figure A.6: Histogram of current values at data points outside the defined beam region. . . . 130 Figure A.7: Non-zero data points after background subtraction. . . . . . . . . . . . . . . . . 131 Figure A.8: Non-zero data points of the beam distribution, before and after applying the island filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132 xiv Figure A.9: Beam parameters setup screen for the Allison scanner GUI. (Image courtesy of Zhang Tong from the Accelerator Physics Department at FRIB). . . . . . . . 134 Figure A.10:Noise correction screen for the Allison scanner GUI. (Image courtesy of Zhang Tong from the Accelerator Physics Department at FRIB). . . . . . . . . . 135 xv CHAPTER 1 INTRODUCTION This chapter overviews high intensity ion accelerators and introduces the Facility for Rare Isotope Beams (FRIB). Electron cyclotron resonance (ECR) ion sources and beam dynamics challenges at front ends of high intensity ion accelerators are discussed. 1.1 High Intensity Ion Accelerators A summary of high intensity ion accelerators operating or being planned around the globe is shown in Fig. 1.1. They serve a variety of purposes including: 1. The intensity frontier in high energy physics, which uses secondary or tertiary beams of neutrinos, muons and kaons (e.g. PIP-II, USA; J-PARC, Japan). 2. Nuclear physics (e.g. FRIB, USA; RIKEN, Japan; RISP, Korea). 3. Spallation neutron production (e.g. SNS, USA; ISIS, UK; ESS, Sweden). 4. Accelerator driven system (ADS) for nuclear waste transmutation and subcritical power generation (e.g. CADS, China; MYRRHA, Belgium). 5. Plasma physics with beam-driven warm dense matter or high-energy-density physics (e.g. FAIR, Germany; HIAF, China). A comprehensive overview as of 2014 can be found in Ref. [1]. A common metric of high intensity ion accelerators is the average beam power on target defined by (maximum beam energy × average beam current). Proton accelerators have already achieved beam power above 1 MW, and multiple MW to >10 MW facilities are presently under construction. For heavy ion accelerators, FRIB is set to advance the continuous-wave (CW) power frontier by two orders of magnitude when it reaches full operating specifications at 400 kW. The FRIB accelerator is described in the next section. 1 Figure 1.1: High intensity ion accelerators (Image courtesy of Prof. Jie Wei at FRIB) [1] 1.2 Facility for Rare Isotope Beams (FRIB) The Facility for Rare Isotope Beams (FRIB) [2] is currently being constructed at Michigan State University. FRIB consists of a CW superconducting radiofrequency (SRF) linear accelerator that aims to accelerate all stable isotopes to energies >200 MeV/u and attain up to 400 kW CW beam power on target. Secondary rare isotope beams generated at the target will be separated in-flight and transferred to experimental halls as fast, stopped or re-accelerated beams. Re-acceleration capabilities are enabled by a SRF ReAccelerator (ReA). A schematic of the FRIB driver linac is shown in Fig. 1.2. It is approximately 470 m in length and consists of over 300 superconducting quarter-wave and half-wave coaxial resonators distinguishable into four types that are suited to different beam velocities [3]. The linac has a folded “paper-clip” layout which facilitates beam collimation after stripping at the end of linac segment 1 to enhance acceleration efficiency. A novel liquid lithium stripper will be employed [4]. The warm front end section of FRIB is illustrated in Fig. 1.3. There will be two electron cyclotron resonance ion sources (ECRIS, see section below) at the front end: ECRIS-1, ARTEMIS 2 Figure 1.2: Schematic of FRIB superconducting driver linac. Note that part of the front end is located on ground level, a complete illustration of the front end is shown in Fig. 1.3 [5], is a normal conducting source that was commissioned in September 2016; ECRIS-2, VENUS- like [6], is a superconducting source that can operate at higher frequency to produce higher currents and higher charge states in heavy ions like uranium. ECRIS-2 will be commissioned in two years’ time. The beam undergoes three stages of acceleration in the front end (see Fig. 1.3): 1. Accelerated by a 20 kV to 30kV potential upon extraction from the ECR. 2. Accelerated to 12 keV/u by an electrostatic gap prior to charge selection. 3. Accelerated to 500 keV/u by the RFQ. Species selection is achieved by a charge selection system that consists of two 90-degree bends for generating sizable charge separation. The lattice includes six electrostatic quadrupoles and the two dipoles include a linear field gradient to enhance dispersion for species separation. After charge selection, the beam line has transverse focusing (magnetic solenoids and electrostatic quadrupoles) and bending elements (electrostatic) for transporting the beam to the RFQ with design parameters, as well as choppers and bunchers for longitudinal prebunching of the beam before it enters the RFQ. A variety of diagnostic devices are installed along the beam line for monitoring beam transport [7]. Wire profile monitors are the most abundant, and play a key role in trajectory corrections throughout the LEBT. There are also scintillator view screens and phase space measurement devices in the form of Allison-type scanners (operational) and pepper-pot emittance meters (not 3 Figure 1.3: Schematic of the FRIB front end with emphasis on the location of diagnostics devices. Abbreviations are listed below: BCM - beam current monitor; BPM - beam position monitor; FC - Faraday cup; MHB - multi-harmonic buncher; PM - profile monitor. yet operational). We will discuss improved diagnostic measurements with Allison scanners in Chapter 5, and with profile monitors in Chapter 6. 1.3 Electron Cyclotron Resonance (ECR) Ion Sources Electron cyclotron resonance (ECR) ion sources [8] are widely used in high intensity heavy ion accelerators because they can stably produce high currents of multiply charged ions with high phase space density from a large variety of stable elements. A schematic of an ECR ion source is shown in Fig. 1.4. The plasma is confined axially by two solenoids that constitute a magnetic bottle, and radially by the sextupole field. Microwaves heat the electrons in the plasma, which in turn collide with ions to produce a high density of high charge state ions that are continuous (DC beam) extracted. The beam extracted from an ECR ion source contains ions of many different charge states. An 4 Figure 1.4: Schematic of an electron cyclotron resonance (ECR) ion source (Image courtesy of Dr. Daniele Leitner at LBNL). example beam from a VENUS-like ECR tuned for producing maximum current of U33+ has charge states ranging from 22+ to 40+ as shown in Fig. 1.5. The magnetic fields of the ECR ion source also induce complication for the dynamics of the extracted beam. Firstly, there exists a strong axial magnetic field where the ions are extracted, which gives ions a large canonical angular momentum. Secondly, due to the sextupole in the source, the extracted density profile of the beam is markedly non-axisymmetric. 1.4 Beam Dynamics Challenges in Accelerator Front Ends Accelerator front ends must preserve beam quality from the source while eliminating unwanted ions. This is necessary to achieve maximum beam power while avoiding machine damage down- stream due to loss of ions in the distribution tail. This section overviews challenges in understanding beam dynamics in accelerator front ends. Ion beams in accelerator front ends have low energies and high intensities. Both factors enhance space-charge effects in the beam evolution. This is further complicated by electron neutralization caused by ionization of residue gas in the beam pipe. Beams are typically extracted at a high magnetic field, which imparts large canonical angular momenta onto beam ions when the 5 Figure 1.5: Measured charge state distribution of a uranium beam from the VENUS ion source [9]. downstream transport is axisymmetric. Thus the transverse phase space of the beam is correlated. There are also many types of lattice elements, which have nonlinear fields with overlapping fringes. Transverse-longitudinal coupling also occurs as a result of beam acceleration and bunching. For heavy ion accelerators which employ ECR-type ion sources, there are further complications due to characteristics of ECR beams. The multi-species beam has strong space-charge forces with intricate inter-species interactions that are partially electron neutralized. Next, beam manipulations are required for species selection. Furthermore, the initial phase space distribution of the beam is non-axisymmetric. Optimal tuning of the front end in the presence of these complexities is essential to preserving beam brightness (current / emittance2) from the source and suppressing downstream beam losses in delicate and expensive SRF cavities. This is critical in high power CW machines like FRIB. 6 CHAPTER 2 A PRIMER ON BEAM PHYSICS This chapter overviews key beam physics results relevant to subsequent chapters in this dissertation. 2.1 Particle Motion and Phase Space The relativistic Hamiltonian of a charged particle in an electromagnetic field is given by [10]: (cid:17)2 where m is the mass, q is the charge, φ is the scalar potential and −→ P is the canonical momentum which is related to the kinetic momentum −→p and vector potential −→ A by: −→ P − q A H = eφ + c m2c2 + (cid:114) (cid:16)−→ (2.1) −→ −→ P = −→p + q A (2.2) In an accelerator, particles evolve about a reference orbit (also called the design trajectory) that is often a plane curve. Therefore, it is convenient to employ the Frenet-Serret coordinate system with vanishing torsion to describe particle motion. The notation commonly employed in the accelerator community is illustrated in Fig. 2.1. s is the path length along the curve from some initial point whereas the local radius of curvature is denoted by ρ(s). After adopting the Frenet-Serret coordinate system, the description of particle motion can be further simplified by noting that the location of beam line elements are defined in s, not t, so it is desirable to use s as the independent variable. Moreover, particles typically only deviate slightly from the reference orbit. Normalizing by the longitudinal momentum and expressing phase space coordinates as differences from those of a reference particle that follows the reference orbit with the design energy, the coordinates become small quantities that enable simplifying approximations. The sequence of coordinate changes outlined above can be achieved via a series of canonical transformations [12]. The resulting Hamiltonian in the absence of applied electric fields and self 7 Figure 2.1: Frenet-Serret coordinate system. Image adapted from Ref. [11] (cid:18) fields is: H = − 1 + x ρ (cid:19)(cid:34)(cid:16)1 + δ 2(cid:17) (cid:18) + px − qAx p0 (cid:19)2 (cid:18) + py − qAy p0 (cid:19)2(cid:35)1/2 (cid:18) − 1 + (cid:19) qAs p0 x ρ + (1 + δ) (2.3) where phase space variables are (x, px, y, py, ∆z, δ). px and py are conjugate momenta to x and y, which are coordinates in the Frenet-Serret coordinate system after the transformations. Longitudinal phase space coordinates ∆z and δ express position and energy deviations from the reference particle as follows: ∆z = −v0(t − t0) p − p0 p0 δ = (2.4) (2.5) where p0 and v0 are the momentum and velocity of the reference particle respectively, and t0 is the time at which the reference particle arrives at position s. Transverse magnetic fields can be expressed solely with As, so Ax and Ay are taken to be zero. Then px and py are related to the 8 mechanical momentum and thus the slope via: −→p · ˆx p0 −→p · ˆy p0 px = py = ≈ x(cid:48)(1 + δ) ≡ dx ds (cid:48)(1 + δ) ≡ dy ≈ y ds (1 + δ) (1 + δ) (2.6) (2.7) Given a beam line element, its vector potential −→ A = As ˆs can be inserted into the Eq. (2.3) to derive Hamilton’s equations that govern the single-particle dynamics. It is a common approximation to assume linear motion, whereupon all terms higher than 2nd order in the Hamiltonian are discarded. Discarding terms in the Hamiltonian guarantees that the resulting linear dynamics is symplectic, which implies nice properties including the preservation of phase space volume that we exploit in Sec. 2.2. With the assumption of small momenta in the linearization of the equations of motion, the difference between x(cid:48), y(cid:48) and px, py is a higher order term that can be neglected. Therefore, it is customary to describe linear dynamics using the following formalism: ˆx(s = s f ) = M(s f |si) ˆx(s = si) (2.8) where the ˆx is the phase space vector (x, px, y, py, ∆z, δ) and the its initial (s = si) and final (s = s f ) values are related by a transfer matrix M. This dissertation mainly studies transverse motion. If the 2nd-order Hamiltonian does not couple transverse and longitudinal phase space coordinates, it can be separated into a transverse part and a longitudinal part. The transverse part will evolve independently as follows: (2.9) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)s=s f x x(cid:48) y y(cid:48) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)s=si x x(cid:48) y y(cid:48) = M⊥ where M⊥ is the 4 × 4 transfer matrix. Furthermore, if the x-x(cid:48) and y-y(cid:48) phase space are not 9 coupled, M⊥ becomes block diagonal, whereupon the dynamics can be further simplified into: x(cid:48) (cid:169)(cid:173)(cid:173)(cid:171) x (cid:169)(cid:173)(cid:173)(cid:171) y y(cid:48) (cid:170)(cid:174)(cid:174)(cid:172)s=s f (cid:170)(cid:174)(cid:174)(cid:172)s=s f x(cid:48) (cid:169)(cid:173)(cid:173)(cid:171) x (cid:169)(cid:173)(cid:173)(cid:171) y y(cid:48) (cid:170)(cid:174)(cid:174)(cid:172)s=si (cid:170)(cid:174)(cid:174)(cid:172)s=si = Mx = My where Mx and My are 2 × 2 matrices. Analytic forms of transfer maps for standard beam line elements such as dipoles, quadrupoles and solenoids are readily available in the literature [13]. 2.2 Symplectic Dynamics, Beam Matrices and Emittance Conservation (2.10) (2.11) (2.12) (2.13) (2.14) 2.2.1 Symplectic Condition Transfer matrices obey the symplectic condition [14]: MTSM = S where S is the block diagnonal symplectic matrix given by: S = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) s 0 0 s ... 0 0 · · · and (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) · · · 0 0 0 ... . . . s 0 0 =(cid:169)(cid:173)(cid:173)(cid:171)0 0 0 0 with s =(cid:169)(cid:173)(cid:173)(cid:171) 0 1 −1 0 (cid:170)(cid:174)(cid:174)(cid:172) 2.2.2 Beam Matrices The second order moments of the distribution of a beam can be naturally represented as the elements of a matrix. Such a matrix is called the beam matrix or sigma matrix, and is commonly denoted 10 by σ. The dimension of σ depends on the dimension of the phase space under question. In this dissertation, we often use the 4 × 4 sigma matrix which comprises the 2 × 2 sigma matrices as its blocks. For a beam with n particles indexed with subscript i that runs from 1 to n, define y(cid:48)(cid:17) i x(cid:48) y σ = x 1 n x x(cid:48) y y(cid:48) i=1 (cid:16) (cid:170)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)i n =(cid:169)(cid:173)(cid:173)(cid:171)σxx σx y σyy =(cid:169)(cid:173)(cid:173)(cid:171) (cid:104)yy(cid:105) n σT x y σyy (cid:104)yy(cid:48)(cid:105) (cid:104)xx(cid:105) = 1 n i=1 where: σxx =(cid:169)(cid:173)(cid:173)(cid:171) (cid:104)xx(cid:105) (cid:104)xx(cid:48)(cid:105) (cid:170)(cid:174)(cid:174)(cid:172) (cid:104)xx(cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) and 2.2.3 Emittance Conservation (cid:170)(cid:174)(cid:174)(cid:172) σx y =(cid:169)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:170)(cid:174)(cid:174)(cid:172) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:104)yy(cid:48)(cid:105) (cid:104)y(cid:48)y(cid:48)(cid:105) xi xi etc. (2.15) (2.16) (2.17) (2.18) Following the treatment of Ref. [15], invariants composed of 2nd-order moments of a beam can be extracted by invoking the symplecticity of the linear transfer map. Consider 4D phase space vector: (cid:16) x = x(cid:48) x y y(cid:48)(cid:17)T Denote x(s = s2) as x2, x(s = s1) as x1. For x2 = M⊥x1, σ2S = x2xT2S = M⊥x1xT1MT⊥S = M⊥σ1SM−1⊥ 11 which shows that σ1S and σ2S are similar matrices and have the same eigenvalues. Therefore, the eigenvalues of σS are invariants of the motion. When the transverse motion is decoupled, the proof above also applies individually to 2D phase spaces x-x(cid:48) and y-y(cid:48) where M⊥ should be replaced by Mx or My, and σ by σxx or σyy. The resulting conserved eigenvalues can be expressed as the familiar 2D emittances: (cid:113)(cid:104)xx(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) − (cid:104)xx(cid:48)(cid:105)2 (cid:113)(cid:104)yy(cid:105) (cid:104)y(cid:48)y(cid:48)(cid:105) − (cid:104)yy(cid:48)(cid:105)2 εx = εy = 2.3 2D Linear Phase Space Motion This section briefly discusses the special status of 2D phase space ellipse in linear motion and its connection to the Twiss parameters. The treatment is inspired by Ref. [16] and omits discussion of periodic motion because accelerator front ends are single-pass systems. Consider an ellipse in x-x(cid:48) phase space (see Fig. 2.2) defined by: γx2 + 2αxx(cid:48) + βx(cid:48)2 = ε (2.21) The equation has redundancy in the sense that the ellipse remains unchanged if one rescales both sides by the same factor. To remove the redundancy, the normalization βγ−α 2 = 1 is conventionally chosen because it renders the ellipse area πε. Hence, given ε, two out of the three parameters α, β and γ completely determines the orientation of the ellipse including its extent and critical points. In practice, α and β are normally the two parameters used because they naturally connects to an envelope description of a beam. The phase space ellipse at s = s1 can be written in matrix form as: (2.19) (2.20) (2.22) (2.23) where (cid:16) x α1 β1 x(cid:48)(cid:17)(cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 det(cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 α1 β1 x(cid:48) (cid:169)(cid:173)(cid:173)(cid:171) x (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) = ε (cid:170)(cid:174)(cid:174)(cid:172) = 1 12 Figure 2.2: Ellipse in x-x(cid:48) phase space [17]. from the normalization constraint. It can be shown that, under uncoupled linear dynamics, an ellipse in x-x(cid:48) phase space transforms into another ellipse with the same area but different orientation. Take the transfer map from s1 to s2 to be Mx, the ellipse at s = s2 becomes: (2.24) x α1 β1 x(cid:48) (cid:16) x x(cid:48)(cid:17) MT (cid:169)(cid:173)(cid:173)(cid:171) x (cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 (cid:170)(cid:174)(cid:174)(cid:172) = ε (cid:170)(cid:174)(cid:174)(cid:172) Mx (cid:169)(cid:173)(cid:173)(cid:171)γ2 α2 (cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 (cid:170)(cid:174)(cid:174)(cid:172) Mx (cid:170)(cid:174)(cid:174)(cid:172) = MT MT  det(cid:169)(cid:173)(cid:173)(cid:171)γ2 α2 (cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 (cid:170)(cid:174)(cid:174)(cid:172) = det (cid:170)(cid:174)(cid:174)(cid:172) Mx (cid:17) det(cid:169)(cid:173)(cid:173)(cid:171)γ1 α1 = det(cid:16)MT (cid:170)(cid:174)(cid:174)(cid:172) det(Mx) α1 β1 α1 β1 α1 β1 α2 β2 x x x α2 β2 = 1 where: is represents an ellipse with a different orientation but the same area because Since a phase space ellipse only changes orientation under uncoupled linear optics, it is the ideal geometric object in phase space for describing a beam. For a beam that is elliptical in phase 13 space, α and β are sufficient to describe the orientation and determine how it evolves under linear optics. α, β and γ, are called the Twiss parameters of a beam. For a beam that is not elliptical in phase space, there is no unique way to define its Twiss parameters. The commonly adopted strategy is to calculate its second order moments and use them to define an ellipse by asking: if a uniformly filled ellipse were to have these second order moments, what would its Twiss parameters be? The correspondence to beam moments under this perspective is summarized by: βx = γx = αx = (cid:104)xx(cid:105) εx (cid:104)x(cid:48)x(cid:48)(cid:105) εx − (cid:104)xx(cid:48)(cid:105) εx (2.25) (2.26) (2.27) where εx is the rms emittance defined in Eq. (2.19). The beam is then represented by such an rms phase space ellipse, with area εx,or 4εx depending on the convention. The rms phase space ellipse should behave just as another other ellipse in phase space, provided the optics is uncoupled and linear. 2.4 Axisymmetry and Conservation of Canonical Angular Momentum We know from Noether’s Theorem [14] that every continuous symmetry of the Lagrangian of a system gives rise to a conserved quantity. For a charged particle in an axisymmetric beam line, the conserved quantity is the canonical angular momentum which can be derived as follows. The Lagrangian for a relativistic charged particle in an electromagnetic field is given by: L = −mc2 γ − eφ + e−→ v · −→ A (2.28) is the velocity and γ is the Lorentz factor. In cylindrical coordinates, the Lagrangian v ≡ d−→x dt where −→ reads: L(r, (cid:219)r, θ, (cid:219)θ, z, (cid:219)z) = − (cid:112)1 − ((cid:219)r2 + r2 (cid:219)θ mc2 2 + (cid:219)z2)/c2 − qφ + q(cid:0)(cid:219)r Ar + r (cid:219)θ Aθ + (cid:219)zAz(cid:1) (2.29) where r is the radial coordinate and θ is the azimuthal angle. For an axisymmetric system, neither φ(r, z) nor −→ A(r, z) depends on θ. Therefore, θ becomes a cyclic coordinate of L, and the canonical 14 angular momentum: ∂L ∂ (cid:219)θ ≡ Pθ = γmr2 (cid:219)θ + qr Aθ (2.30) is a conserved quantity. Since the conservation of Pθ holds for every particle, for a single species mononenergetic beam: (cid:68) (cid:69) (cid:104)Pθ(cid:105) = γm r2 (cid:219)θ + q (cid:104)r Aθ(cid:105) = const (2.31) This relation places a strong constraint on transverse single-particle dynamics in axisymmetric focusing systems. In practice, asymmetries due to misalignments and field errors will induce some degree of violation to this constraint. 15 CHAPTER 3 SELF-CONSISTENT SIMULATIONS This chapter describes the development of a set of self-consistent simulation tools that are readily adaptable to simulations of accelerator front ends with high levels of detail. Simulations are crucial to understanding the beam dynamics because the effects of non-ideal applied fields, intricate lattice geometries and intense non-uniform space-charge cannot be described purely analytically. We applied these simulation tools to analysis the VENUS-like beam line at the FRIB front end at full power. Many effects were discovered with causes including beam magnetization, space charge, field nonlinearities and their interplay. 3.1 Warp Code Warp [18, 19] is a open source particle-in-cell (PIC) code well-suited for self-consistent mod- eling of high intensity ion beams. In additional to an extensive list of standard beam line elements such as dipoles, quadrupoles and solenoids; the code enables time-dependent fields and detailed specifications of conductor surfaces and applied fields from imported field maps. This capability, discussed in more detail in the section below, is particularly useful for modeling accelerator front ends given the variety of beam elements involved. Warp also allows simulations with an arbitrary number of particle species including multiple ions with multiple charge states, which is necessary in the modeling of beams extracted from an ECR ion source. Warp simulations can be run in various modes including: 1) full 3D; 2) axisymmetric (r − z); and 3) transverse slice (x − y). The code has a variety of electrostatic and electromagnetic field solvers (including capabilities like adaptive and static mesh refinement) and is built around an adaptable python interface. Assuming there is no axisymmetry, the user can choose between full 3D mode and slice mode depending on whether longtudinal self-field is important. The beam can be initialized with pseudoequilibrium distributions adapted for space-charge[20], or with distributions obtained from other simulations or experimental data. 16 We have written a set of wrapper modules in Python that allows one to conveniently: 1. Set up beam lines by specifying the locations of focusing elements, conductors and beam collimators at several levels of detail. 2. Generate initial beam conditions given the current and phase space parameters of each ion species. 3. Assign diagnostics and record simulation results. 4. Control global properties and parameters of the simulation. Transverse x-y slice simulations are employed that adjust timesteps of individual particles advanced to advance slice-to-slice for a specified lattice advance step. Self-fields are generated by solving Poisson’s equation on a discretized spatial mesh with the correct transverse conducting aperture of the accelerator lattice. Applied electric and magnetic fields (including longitudinal electric) are imported via gridded potential and field elements. This model neglects relatively weak longitudinal self-fields. However, it retains correct longitudinal applied forces from lattice elements. The field solutions are consistent with local bend radii ρ of the tight dipoles used for charge selection in the front end (see overview in Chapter 1). These modules constitute a versatile code framework for simulations of ion accelerator front ends using Warp. The framework is highly adaptable to different accelerator facilities and different operating points, which helps expand its utility within the community. While we have not employed this functionality in the studies below, Warp also supports parallelization to allow high-detail runs and parametric studies. 3.2 Beam Line Element Models The transverse spatial extent of low energy beams in accelerator front ends can be similar to the extent of the beam pipe and apertures of focusing elements. In fact, it is not uncommon to have ∼ 1% beam loss in front ends, which indicate scraping during beam transport. Therefore, the outer particles in the beam distribution can evolve beyond the good field radius of beam line elements and 17 experience substantial field nonlinearities. Furthermore, large transverse-to-longitudinal aspect ratios in focusing elements increase nonlinearities and cause fringe fields to leak out far beyond the longitudinal extent of the element, sometimes overlapping with those from adjacent elements. The Warp simulations properly model these effects. Fig. 3.1 shows the CST Studio[21] model of S4 solenoids employed in the front end at FRIB. One can observe that field nonlinearties grow significantly as one goes to larger radii. Fringe fields also extend up to 15 cm outside the solenoid on both sides, which is half as long as the coil length. The CST Studio model of Q7 electrostatic qaudrupoles are shown Fig. 3.2 and Fig. 3.3. While the field is relatively more uniform, the Q7 in the triplet assembly contains collimators that can intercept beam particles that might be lost on electrode structures. The Warp simulations employ particle scraping on the detailed aperture structures like the electrostatic quadrupoles. Both examples show that it is necessary to implement detailed field models to perform realistic simulations of accelerator front ends. 3.3 Simulations of FRIB Front End 3.3.1 Initial Conditions & Loading x =(cid:10)x2(cid:11)). The 35 kV extraction potential sets the kinetic energy of each species, and 2 Table 3.1 lists beam parameters used in the Warp simulations. An initial axisymmetric waterbag distribution adapted for space-charge [20] is injected with no centroid offset at the extraction point of the ECR. Each species’ envelope extent fills the puller electrode aperture (i.e. Rpuller = 2σx where σ the longitudinal velocity spread and radial (thermal component) emittance arise from a 3 eV ion temperature. Since the beam aperture radius at extraction is relatively small in the sextupole field of the ECR ion source, we assume nearly constant beam canonical angular momentum (cid:104)Pθ(cid:105) in the initial transport from the ECR. The average canonical angular momentum of each species is set birth Bbirth/2, where z(peak 1) Bz(r = 0)dz is the average axial B-field between the peaks of the solenoidal field of the ECR. Since the B-field at launch differs from Bbirth, for (cid:104)Pθ(cid:105) = const as specified, an by assuming all ions are born at the same magnetic field, i.e. (cid:104)Pθ(cid:105) = q(cid:10)r2(cid:11) Bbirth =Þ z(peak 2) 18 Figure 3.1: CST Studio model of S4 solenoids at FRIB and the calculated fields at different radii within the clear bore aperture radius R = 164.2/2 mm. angular velocity is injected upon each ion of the species with: Here, [Bρ] is the particle rigidity, and we assume(cid:10)r2(cid:11) birth launch Bbirth − Blaunch (cid:48) 0 = θ (cid:34) (cid:10)r2(cid:11) (cid:10)r2(cid:11) (cid:35) birth =(cid:10)r2(cid:11) /(2[Bρ]). (3.1) launch. Simulated beam species and currents for U operation are given in Ref. [9]. 3.3.2 Lattice Design & Matching Methodology The initial beam line from the ECR ion source can be divided into two distinct sections: the axisymmetric transfer line from the ECR to the charge selection system (CSS), and the CSS which produces large dispersion for species selection. 19 395.2 mm296 mm164.2 mm456.1 mmyokecoilhalf lengthRadial FieldAxial FieldCross Section (Side View)Cross Section (Perspective View) Figure 3.2: CST Studio model of a Q7 electrostatic quadrupole triplet at FRIB with dimensions and potential plots. For species selection with no loss in target species and preservation of beam quality, the CSS should: 1) Be linearly achromatic. 2) Generate large dispersion with small βx (x-betatron lattice function) at collimation point to maximize selection resolution with vertical slits. 3) Output the beam with well-controlled envelopes for downstream transport. Requirements (1) and (3) can be met in an idealized linear optic limit with a symmetric beam envelope in the CSS by exploiting mirror symmetry about the axial mid-point of the CSS (strength of each element equals that of its mirror counterpart) and setting D(cid:48) = αx = αy = 0 (D = bend- plane dispersion function and αj = −β(cid:48) j/2, j = x, y) at the axial mid-point. The axisymmetric 20 70 mm 90 mm200 mm239 mm96 mm100 mm70 mm150 mm130 mmBetweeen electrode & endplateEnd of electrodeMid-quad Figure 3.3: CST Studio model of a Q7 electrostatic quadrupole triplet at FRIB. Dimensions and potential contours at various axial locations are shown. Table 3.1: Beam Parameters of U34+ based on the VENUS-like source. For definitions, see also Sec. 3.3.3 and ??. Temperatures and envelope parameters are the same for all ion species (indices suppressed). Canonical angular momentum, radial normalized emittance, and kinetic energy/rigidity vary with ion species. Quantity Transverse Temperature Canonical Angular Momentum Radial Normalized Emittance Radial Envelope Radial Envelope Angle Kinetic Energy Rigidity Longitudinal Temperature T Symbol (cid:104)Pθ(cid:105) /(mc) rms √ ε nr 2σx σr = σ(cid:48) r = dσr/dz T(cid:107) Value at Launch 3 eV 0.305 mm-mrad 0.015 mm-mrad 2.82 mm 0 3 eV Ek (before/after ES Gap) [Bρ] (before/after ES Gap) 5.00 / 12.00 keV/u 7.13×10−2 / 1.10 ×10−1 Tm initial beam envelope conditions are adjusted to satisfy requirement (2) and achieve an “attractive” beam envelope devoid of large excursions. The CSS design procedure was carried out with the code MADX [22] using a hard-edge lattice and single-particle dynamics. Results are shown in Fig. 3.4. The lattice is achromatic with zero 21 in = 0) and final dispersion (Dout = D(cid:48) initial (Din = D(cid:48) out = 0), and the envelopes are symmetric about mid-plane with lattice functions satisfying βin = βout and αin = −αout to ideally recover axisymmetry on exit. The 1 m-long 90-degree bend dipoles have slanted poles that generate a focusing strength equivalent to a superimposed quadrupole with κ0 = 0.365 m−2. The three electrostatic quadrupoles (ESQs) in the triplet have axial length 20.7 cm with κ1 = 6.62 m−2, κ2 = -14.8 m−2 and κ3 = 7.80 m−2. The initial envelope has βin = 4.383 m and αin = 0.306, which (see below) are consistent with upstream matching from the ECR. Figure 3.4: CSS lattice design used in Warp simulations. The transport line from the ECR must be tuned to deliver the required lattice functions entering the CSS. We assume no initial centroid offset (Din = D(cid:48) = 0). Applying εx = εeff/2 (see x/εx to convert between lattice functions and envelope Sec. 3.3.3), βx = σ size, the excitation of the two upstream solenoids are tuned to match the required lattice functions. x/εx and αx = −σxσ(cid:48) 2 in 22 s (m)0.01.2.3.4.5.6.7.8.0.00.51.01.52.02.53.0Dx(m)βxβyDx69.0970.1371.1772.2173.2574.2990o dipole90o dipoleESQ TripletESQ Triplet 3.3.3 Envelope Model for Multi-Species Beam Transport Beam evolution in the axisymmetric transport line from the ECR to the 1st dipole of the CSS can be described by a multi-species envelope equation: (cid:48)(cid:48) r j = σ +  qjV(cid:48) 2Ek j (cid:48) r j + σ qjV(cid:48)(cid:48) 4Ek j σr j − Q js fs σr j + σ 2 rs 2 r j + (cid:18) qj Bz0 (cid:19)2 j /(cid:16) (cid:17)2 mj βjc . σr j 2mj βjc 2 + (cid:104)Pθ(cid:105)2 rms ε r j 3 r j (3.2) (3.3) σ s,species √ 2σx j = σ Details of this model can be found in Ref. [23]. In brief: j indexes the jth species; qj and mj are the charge and mass of the species; βj and Ek j are the relativistic factor and kinetic energy √ 2(cid:104)x2(cid:105)j is the radial sigma of the species charge profile; V and of the species; σr j = B0z are the on-axis applied potential (of accelerating gaps) and the axial magnetic field which both can vary in s; (cid:104)Pθ(cid:105)j = const is the canonical angular momentum of the species; and fj is the fractional electron neutralization factor ( fj = 0 corresponds to neutralized). The canonical angular momentum measured by (cid:104)Pθ(cid:105) j /(mj βjc) contributes to the beam’s expansion in the same manner as the thermal geometric emittance (cid:118)(cid:116) Tj rms r j = ε 2 j c2 R mj β rms nr j at injection with ε = βj εr j = const (conserved normalized emittance) set by the initial value emerging from the source. In Eq. (3.3), R = 4 mm is the ECR extraction aperture and the value of rms ε r j corresponds to a uniform density beam species filling the aperture. In this model 2 ε eff,j = ε rms r j 2 + (cid:104)Pθ(cid:105)2 (mj βjc)2 j (3.4) constitutes an effective phase-space volume of the beam. The envelope model also assumes Gaussian density profile species are maintained in the evolution. The species kinetic energy Ek is set consistently with the applied on-axis potential V of electrostatic gaps. The envelope equations and energy equations of all species constitute a system of ordinary differential equations (ODEs) that can be numerically integrated along with an equation for the species energy using the linear 23 field components extracted from the 3D element models [B0z(s) = Bz(r = 0, s) for solenoids and V(s) for electrostatic gaps]. The integration begins at the ECR extraction point and ends at the entrance of the first dipole. Figure 3.5 shows that the envelope equation achieves good agreement with self-consistent Warp simulations both with and without space charge. A python based envelope package was constructed within the Warp interface to integrate the coupled nonlinear envelope equations 3.2 and the energy gain equation. This allowed the same lattice setup to be used as the full simulations. Even with little numerical optimization, this envelope model is more than 100 times faster than relatively efficient Warp transverse slice simulations of the same section. This efficiency allows parametric exploration of improved matching. Successively finer scans of the two solenoid fields (minimum increments being 0.002 Tesla) are carried out until target envelope conditions are matched by hitting target envelope parameters entering the achromatic CSS. Figure 3.5: Comparisons between numerical solutions of the envelope equation and Warp simulation results in the transport line upstream of the CSS. These correspond to simulation cases 3 and 4 described in Sec. 3.3.4. A caveat is appropriate for the matching procedure above: the lattice design is not solely dictated by the CSS, because the design incoming lattice functions may be unattainable within the range of possible solenoid ex citations. Fig 3.6 summarizes accessible lattice functions at the entrance of the first dipole for different space charge neutralization values. Colored dots indicate points where 24 ECRCSSsolsolgrated gapOn-Axis Peak Field35kV50kV Figure 3.6: Accessible lattice functions β and α at the entrance of the first dipole when the Bpeak in both solenoids varied independently between 0 T and 1.5 T. Results are shown for different space charge neutralization factors. Observe that there are regions of initial Twiss parameters that are inaccessible, and the region of accessible Twiss parameters decreases with decreasing neutralization factor accessible solenoid excitations can be found for the two solenoids connecting the ECR to the CSS to achieve the indicated lattice function values βx = βy ≡ β and αx = αy ≡ α = β(cid:48)/2. If the CSS operating point is incompatible with the matching capabilities of the solenoids, an alternative CSS operating point with attainable initial lattice functions must be selected. 3.3.4 Simulation Cases Fields from lattice elements in the axisymmetric transport line are imported from r-z electrostatic and magnetostatic simulations with Poisson [24] and CST Studio [21]. Warp also allows hard edge elements, and for present purposes, all elements in the CSS are modeled as hard-edge equivalent here to simplify comparisons to the MADX design lattice. Simulations with space charge include a uniform neutralization factor of 75% ( f = 0.75) except inside the grated gap (protected by a guard electrode included in V) and the CSS lattice downstream of the 1st electrostatic quadrupole. To illustrate effects of space charge and canonical angular momentum in the CSS, we launch 25 beams with the same effective phase-space area (εeff) as calculated from Eq. (3.4), but with different relative compositions from thermal emittance and canonical angular momentum. A “magnetized” beam has phase space contributions consistent with Table 3.1. A “thermal” beam has the same phase-space area (emittance) but with zero canonical angular momentum ((cid:104)Pθ(cid:105)j = 0). That is, for the thermal beam is unmagnetized with increased thermal temperature Tj for the same effective phase-space area measured by εeff,j. The following cases were run: 1) Thermal beam, no space charge. 2) Thermal beam, with space charge. 3) Magnetized beam, no space charge. 4) Magnetized beam, with space charge. Note that case 4 is close to what is expected in the laboratory since canonical angular momentum dominates the effective phase space volume. 3.3.5 Warp Simulation Results: Case Comparisons Before presenting how the case comparisons of Sec. 3.3.4 illuminate different effects in the beam dynamics, it is helpful to observe the phenomena of species separation and collimation in the charge selection section as illustrated in Fig. 3.7. The Warp simulation results correspond to the realistic case 4 with charge neutralization factor f = 0.75. For the two-species transport (U33+ and U34+) planned for the full power operation at FRIB, the optics is tuned with a hypothetical U33.5+ ion as the reference. As seen in the upper panel of Fig. 3.7 showing beam species extents in the x-plane, only U33+ and U34+ (purple and dark green) species make it past the final 4-jaw collimator aperture. Numerous species, including high space-charge intensity oxygen, are lost in the dipole and on the slits incorporated within the following electrostatic quadrupole triplet structure. Numerous species are also lost on the gate valve structure. Neutralization ends at the start of the first electrostatic quadrupole which has strong electric sweeping fields. This is also the axial location 26 of a gate valve. Corresponding line-charge measures of space-charge intensity are given in the middle panel of Fig. 3.7. This shows that most cumulative space-charge intensity emerging from the ECR is scraped. Simple analytic estimates of the depression in local particle phase advance due to beam space-charge are included in the lower panel of the figure. Only the region downstream of the four jaw collimator is show since the estimates are based on single species beam models. The depression measure shown is σ/σ0 the ratio of space-charge depressed phase advance to the single particle phase advance (σ/σ0 = 1 for single particle dynamics with zero beam intensity and σ/σ0 → 0 at the transport limit with full space-charge depression). Results are consistent with a roughly 10% space-charge depression (σ/σ0 = 0.9). We note that this operating point would be representative of a highest intensity facility operation. While space-charge is not negligible post species separation, it is not large and will further drop post acceleration in the RFQ. Lower current beams presently being used in facility commissioning have lower space-charge intensity and single particle dynamics downstream of species separation will be a good approximation in these cases. Warp PIC simulations of the 4 cases in the CSS are contrasted to understand deviations from the MADX optical design due to space charge, non-constant emittance, beam canonical angular momentum, and nonlinear optical effects. In this case the four simulations are run through the 2nd dipole of the CSS. Identifying the source of deviations from the ideal design with reference to the four test cases outlined in Sec. 3.3.4 will help guide future optimization of the CSS. Figure 3.8 shows the rms beam envelope (a) and emittance (b) evolution of the reference species U34+, and the centroid offset of U33+ along with mid-plane xy phase space projections (c). Contrasting the simulation cases, we observe the following behavior: • Beam rotation due to canonical angular momentum. • In cases with space charge, non-reference species have centroid offset upon exiting the CSS. • RMS emittance growth in the CSS which is mostly reversed on exiting the CSS. • RMS emittance exchange between the x- and y-planes for magnetized beams. 27 Figure 3.7: Warp simulation results of Uranium species separation and collimation in the beam evolution from the ECR source extraction through the midpoint of the charge selection system showing collimation to two target species for high intensity FRIB operation. Distribution x y projections in Fig. 3.8c illustrate how a magnetized beam becomes tilted in real space at the CSS midplane due to canonical angular momentum ((cid:104)Pθ(cid:105) (cid:44) 0; dropping species subscripts henceforth) whereas a thermal beam ((cid:104)Pθ(cid:105)) does not. Large initial (cid:104)x y(cid:48)(cid:105) and (cid:104)x(cid:48)y(cid:105) moments at the start of the dipole in the case of a magnetized beam causes (cid:104)x y(cid:105) to evolve in the non-axisymmetric lattice which manifests as a significant distribution tilt in the middle of the CSS. This tilt suggest that more elaborate shaped electrodes for final beam collimation may be beneficial. However, effective tuning of such a system would require improved understanding of the initial state of the beam emerging from the ECR including the degree of magnetization (value of (cid:104)Pθ(cid:105)). 28 Figure 3.8: Rms envelops σx,y a) and normalized rms-emittance εx,y b) of the reference species U34+ in the CSS. c) shows the centroid of U33+ and x y-particle projections of both target species at the CSS mid-point for two cases. Solid curves are x-plane and dashed curves are y-plane. Colors flag Warp simulation cases. Black solid lines denote MADX design values in a) and c). The axial midplane of the CSS is indicated in c) (z in FRIB coordinates with origin outside machine). 29 U33+U34+Therm w/ SCMag w/ SCa)b)c)(mid CSS)(mid CSS)xyxy While space charge of the multi-species beam has little impact on the envelope evolution in x and y, insofar as the matching is carried out equivalently, the cases with and without space charge in Fig. 3.8c show that space charge creates a ∼mm shift in the (cid:104)x(cid:105) centroid of U33+. However, the canonical angular momentum has little impact on either the envelope extent or the dispersive shift of the centroid. The space charge induced centroid shift at mid CSS suggests that optimal adjustment of the four jaw collimator will vary with ECR current extracted. This further complicates the impact of magnetization induced tilt on the beam x y profile. Figure 3.8b shows significant emittance growth inside the 1st bending dipole that is largely removed on exiting the 2nd dipole of the achromatic CSS beamline for simulation cases 1) – 3). Slight asymmetries for case 2) (thermal beam with space-charge) are likely induced by nonlinear space-charge fields that vary within the CSS evolution upstream and downstream of the midplane. Effects from nonlinear terms in the dipole arising from large x/ρbend may also contribute to asymmetries. The emittance increase is stronger for thermal beams than for magnetized beams – regardless of space charge. For case 4) (magnetized beam with space-charge) there is significant emittance exchange between the x- and y-planes that results in significantly different final emittance values. For single particle dynamics beams in quadrupoles and dipoles, there is no coupling between planes regardless of beam rotation. For beams aligned with the x-y symmetry axes with charge density nearly constant on elliptical symmetry surfaces, space-charge only induces weak coupling. However, for a rotated x y beam profile, space-charge coupling between the planes is enhanced. More analysis is necessary to probe if this exchange comes with net growth due to nonlinearities but the effect appears largely exhange. Understanding effects like these and predicting the x y tilt angle of the beam will help guide more effective choices of apertures and slit positions for optimal species collimation at higher intensities anticipated for high-power FRIB operating modes. 3.4 Further Work Simulation studies of the FRIB front end can be improved by including the physics of ion extraction at the ECR. We have completed the first step by including the extraction electrodes 30 into the lattice as conductors. A schematic of the electrodes of the VENUS-like ECR is shown in Fig. 3.9. We have yet to study the beam dynamics with the extraction electrodes included, and we expect to see interesting effects because experiments have shown substantial changes in beam conditions downstream when the extraction electrode positions are adjusted. Nonlinearities in extraction also have maximal effect near the emitting plasma meniscus of the source. Proper resolution of these effects will likely require r-z or full 3D (to include sextupole distortions) of the beam in the extraction region. A further improvement will be the use of realistic initial beam conditions derived from ECR simulations. A group at JINR, Dubna, Russia has developed sophisticated simulations of ECR ion sources [25, 26]. We briefly interacted with Dr. Vladimir Mironov on a visit to FRIB. If such steady-state simulations were carried out for the ECRs at FRIB, they could be run to accumulate a population of extracted ions. Then, Warp simulations can be set up to birth randomly chosen extracted ions and model their extraction in 3D mode until steady state is reached. Results of such simulations (perhaps also needing proper plasma meniscus modeling) could reveal much of the ion beam extraction physics that is presently omitted. Another important piece of physics that is currently missing from the Warp simulations is self- consistency in the electron neutralization model. Fig. 3.10 shows that the beam evolution heavily depends on the degree of neutralization, so a clear understanding is essential to front end design and tuning. Not only do we lack a detailed model electron neutralization effects, even in the primitive model of a homogeneous neutralization factor, we are not certain what it should be. Tentative steps towards developing a simulation neutralization model were presented in the Ph.D thesis of Daniel Winklener [27]. Extensions of such complicated models may be necessary for more fully realistic modeling. A critical detail may be transitions between neutralized and unneutralized (swept by electrostatic fields) regions. Lack of uniformity in neutralization factors there may enhance the nonlinear component of space-charge forces. A final problem concerns the generation of 4D coupled phase space distributions. The ini- tialization of uncoupled distributions with high space charge intensity is treated comprehensively 31 Figure 3.9: Extraction geometry of the VENUS-like ECR ion source. The electrode in red has negative voltage to prevent backflow of electrons. (Image courtesy of Dr. Daniel Winklener.) Figure 3.10: Beam sizes in the post-extraction beam line downstream of the VENUS-like ECR at FRIB, obtained from Warp simulations with different homogeneous neutralization factors. 32 In Sec. 3.3.1, coupling is inserted into an uncoupled distribution by the addition in Ref. [20]. of an angular velocity as described by Eq. (3.1). Alternatively, we have derived in Appendix C a method of generating a 4D coupled distribution with any set of 2nd-order moments via linear transformation. It should be beneficial to study how these treatments connect with each other and whether they produce simulation artifacts. 33 CHAPTER 4 THEORY OF TRANSVERSE COUPLED BEAM DYNAMICS, WITH APPLICATIONS TO THE FRIB FRONT END At the time of this writing, the analysis in this chapter pertaining to beams with rotational symmetry is under preparation for submission for publication [28]. 4.1 Beams with Rotational Symmetry This section develops analytic tools to describe moments of beam distributions with discrete or continuous rotational symmetry. Their implications and possible applications are also discussed. 4.1.1 Effect of Symmetry on Beam Moments Define transverse beam moments (cid:68) xb1 x(cid:48)b2 yb3 y (cid:48)b4(cid:69) ≡ Þ ∞ −∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ −∞ F(x, x(cid:48), y, y(cid:48))xb1 x(cid:48)b2 yb3 y(cid:48)b4dxdx(cid:48)dydy(cid:48) −∞ −∞ F(x, x(cid:48), y, y(cid:48))dxdx(cid:48)dydy(cid:48) Þ ∞ −∞ −∞ −∞ −∞ where F(x, x(cid:48), y, y(cid:48)) is the distribution function in transverse phase space and b1, b2, b3, b4 ∈ Z≥0. The moment is said to be of k-th order with k = b1 + b2 + b3 + b4. Note that the moments are calculated assuming each variable has zero mean, i.e. all four first order moments vanish. This corresponds physically to the beam centroid following the design (reference) orbit and must be true for beams with rotational symmetry, as will be proved in subsequent arguments. For a beam with n-fold rotational symmetry (i.e. Cn), the distribution function F(x, x(cid:48), y, y(cid:48)) is invariant under a rotation by θ = 2π/n. This implies all transverse beam moments remain unchanged under the transformation: (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) x x(cid:48) y y(cid:48) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) cos θ 0 − sin θ 0 cos θ 0 (cid:55)→ 0 − sin θ 34 sin θ 0 cos θ 0 0 sin θ 0 cos θ (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) x x(cid:48) y y(cid:48) (4.1) For an axisymmetric beam (i.e. beam with continuous rotational symmetry), the invariance of F and beam moments holds for any rotation angle θ. This is the only piece of information that we need to extract from the rotational symmetry of the beam. No assumption on space-charge intensity is made or implied. The rest of this section explores its consequences and applications. Given a beam distribution F with rotational symmetry, its beam moments must satisfy certain constraints in order for them to be invariant under the corresponding rotational transformations. Since rotations amount to linear transformations among phase space coordinates, one expects the constraints to arise from systems of linear equations among beam moments of the same order. Formally, denote (x, x(cid:48), y, y(cid:48)) as (x1, x2, x3, x4) and given the rotation matrix: where θ is a rotation angle of the symmetry, the relation: R = 0 − sin θ − sin θ cos θ 0 0 cos θ 0 sin θ 0 cos θ 0 0 sin θ 0 cos θ (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) R3j xj(cid:170)(cid:174)(cid:172)b3(cid:169)(cid:173)(cid:171) 4 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) R1j xj(cid:170)(cid:174)(cid:172)b1(cid:169)(cid:173)(cid:171) 4 R2j xj(cid:170)(cid:174)(cid:172)b2(cid:169)(cid:173)(cid:171) 4 C(a1, a2, a3, a4)(cid:68)    = j=1 a1 a2 a3 a4 (cid:68) xb1 1 xb2 2 xb3 3 xb4 4 = (cid:69) (cid:42)(cid:169)(cid:173)(cid:171) 4 (cid:69) j=1 (cid:68) xb1 1 xb2 2 xb3 3 xb4 4 j=1 holds for every b1, b2, b3, b4 ∈ Z≥0. Eq. (4.3) can be rewritten in the form: j=1 (4.2) (4.3) (4.4) R4j xj(cid:170)(cid:174)(cid:172)b4(cid:43) (cid:69) xa1 1 xa2 2 xa3 3 xa4 4 where a1, a2, a3, a4 ∈ Z≥0 and a1 +a2 +a3 +a4 = b1 +b2 +b3 +b4. Each coefficient C(a1, a2, a3, a4) is a polynomial of matrix elements of R, which means it is either zero or a polynomial of degree k in sin θ and cos θ. Thus, the set of all k-th order beam moments in Eq. (4.4) for non-negative integers bj form a 35 linear system: (cid:68) (cid:68) (cid:68) (cid:68) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 3, x0 2, x0 1 x0 xk 4 xk−1 3, x0 2, x0 x1 1 4 xk−1 x0 3, x0 2, x1 1 4 ... 3, xk−1 1 x0 x0 2, x1 4 1 x0 x0 2, x0 3, xk 4 (cid:69) (cid:68) (cid:69) (cid:69) (cid:69) (cid:69) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:68) (cid:68) (cid:68) (cid:68) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 3, x0 2, x0 1 x0 xk 4 xk−1 3, x0 2, x0 x1 1 4 xk−1 x0 3, x0 2, x1 1 4 ... 3, xk−1 1 x0 x0 2, x1 4 1 x0 x0 2, x0 3, xk 4 (cid:69) (cid:68) (cid:69) (cid:69) (cid:69) (cid:69) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = A(θ) (4.5) that describes relations among beam moments resulting from invariance under rotation by angle θ. A(θ) is a function of θ because the coefficients C(a1, a2, a3, a4) in Eq. (4.4) depend on θ, and so the relations change with the order of rotational symmetry. While Eq. (4.5) is a general result that already contains the properties of beam moments under rotational symmetry, the expressions are complicated because most elements in A(θ) are non-zero. It is difficult to know whether or how the equations in Eq. (4.5) can be manipulated to obtain simple equality relations that are practically useful. One trick to achieving some simplification uses the with 0 ≤ m < n. Therefore, for a beam with Cn symmetry, fact that Eq. (4.5) holds for all θ = 2mπ n a sum over the set of θ that preserves symmetry can be taken to obtain: 1 x0 2, x0 3, x0 xk 4 xk−1 3, x0 2, x0 x1 1 4 xk−1 x0 3, x0 2, x1 1 4 ... 3, xk−1 x0 1 x0 2, x1 4 2, x0 x0 1 x0 3, xk 4 1 x0 2, x0 3, x0 xk 4 xk−1 3, x0 2, x0 x1 1 4 xk−1 x0 3, x0 2, x1 1 4 ... 3, xk−1 x0 1 x0 2, x1 4 2, x0 x0 1 x0 3, xk 4 (cid:18)2mπ n−1 (cid:68) (cid:68) (cid:68) (cid:69) (cid:69) (cid:69) (cid:69) (cid:69) (cid:69) (4.6) m=0 (cid:69) (cid:69) (cid:68) (cid:68) = 1 n (cid:69) A n (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) Each coefficient contains summations of polynomials P with two trigonometric variables: (cid:68) (cid:19) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:68) (cid:68) (cid:68) (cid:68) (cid:18)2mπ n (cid:19)(cid:19) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:69) n−1 m=0 (cid:18) (cid:18)2mπ (cid:19) n sin P ,cos which may be reducible. For k = 2, this method invokes the same identities as the alternative proof presented in Appendix B, and enables the derivation of simple relations among 2nd order 36 (4.7) beam moments. The k = 2 results suggest that simple equality relations exist for k > 2, but the method of Eq. (4.6) is difficult to generalize to k > 2 because the elements of A [i.e. coefficients C(a1, a2, a3, a4) of Eq. (4.4)] grow exceedingly complicated as the moment order k increases. We employ a more efficient, yet less intuitive, method to derive moment relations from symmetry in Sec. 4.1.2. 4.1.2 Derivation of Moment Relations from Symmetry Having developed the intuition on the relationship between rotational symmetry and beam moments in Sec. 4.1.1, we devised a more abstract and efficient method to derive clean expressions of the equality relations. The method was inspired by Ref. [16], and is readily applicable to moments of arbitrary order in beams with any rotational symmetry. Define two complex conjugate pairs composed of transverse phase space coordinates: (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) w w w(cid:48) w(cid:48) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 0 0 0 e−iθ Upon a coordinate rotation by θ,(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) w w w(cid:48) w(cid:48) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:55)→ eiθ 0 0 0 0 e−iθ 0 0 0 0 eiθ 0 In analogy to Sec. 4.1.1, we construct complex moments(cid:68) Z≥0 and a1 +a2 +a3 +a4 = k. The real and imaginary parts of(cid:68) (cid:48)a3 w(cid:48)a4(cid:69) (cid:55)→ ei(a1−a2+a3−a4)θ(cid:68) wa1 wa2 w (cid:68) wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69) where a1, a2, a3, a4 ∈ wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69) each comprises (cid:48)a3 w(cid:48)a4(cid:69) wa1 wa2 w (4.10) sums of physical k-th order beam moments. Upon a rotation by θ, the complex moment undergoes the following transformation in accordance with Eq. (4.9): (4.8) (4.9) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) w w w(cid:48) w(cid:48) ≡ (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) x + iy x − iy x(cid:48) + iy(cid:48) x(cid:48) − iy(cid:48) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 37 If rotation by θ is a symmetry of the beam, (cid:68) wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69) remains unchanged upon the (cid:68) transformation which gives: (cid:48)a3 w(cid:48)a4(cid:69) (cid:48)a3 w(cid:48)a4(cid:69) due to symmetry. For every 4-tuple(a1, a2, a3, a4) where ei(a1−a2+a3−a4)θ (cid:44) 1,(cid:68) = ei(a1−a2+a3−a4)θ(cid:68) wa1 wa2 w wa1 wa2 w Eq. (4.11) is the key equation that efficiently generates equality relations among beam moments (4.11) wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69) = 0 must follow; this gives: Re(cid:16)(cid:68) Im(cid:16)(cid:68) (cid:48)a3 w(cid:48)a4(cid:69)(cid:17) (cid:48)a3 w(cid:48)a4(cid:69)(cid:17) wa1 wa2 w wa1 wa2 w = 0 = 0 (4.12) (4.13) For a beam with n-fold rotational symmetry, each 4-tuple (a1, a2, a3, a4) satisfying the following conditions: Condition 4.1.1. a1 + a2 + a3 + a4 = k. Condition 4.1.2. (a1 + a3) − (a2 + a4) > 0. Condition 4.1.3. (a1 + a3) − (a2 + a4) (cid:44) 0 (mod n). gives two unique relations on the k-th order transverse moments of the beam, in the form of Eq. (4.12) and Eq. (4.13). Condition 4.1.1 is merely a restatement of the moment being k-th order, while condition 4.1.3 asserts ei(a1−a2+a3−a4)θ (cid:44) 1. Condition 4.1.2 removes redundant relations for the following reason. Suppose a 4-tuple (a1, a2, a3, a4) fulfills conditions 4.1.1 and 4.1.3 but not 4.1.2: (cid:68) (cid:68) (cid:48)a3 w(cid:48)a4(cid:69) (cid:48)a4 w(cid:48)a3(cid:69) wa1 wa2 w wa2 wa1 w Its complex conjugate: (a1 + a3) − (a2 + a4) < 0 (4.14) (a2 + a4) − (a1 + a3) > 0 (4.15) = 0 = 0 38 is equivalent to and contains the same information as: (cid:68) (cid:48)b3 w(cid:48)b4(cid:69) wb1 wb2 w (4.16) Therefore, the relations generated by (a1, a2, a3, a4), which violates condition 4.1.2, are already covered by another 4-tuple (b1, b2, b3, b4) that satisfies condition 4.1.2. (b1 + b3) − (b2 + b4) > 0 = 0 angle θ keeps(cid:68) For an axisymmetric beam (i.e. beam with continuous rotational symmetry), rotation by any wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69) unchanged. Therefore, ei(a1−a2+a3−a4)θ (cid:44) 1 can always be satisfied for some θ. In this case, conditions 4.1.1 and 4.1.2 suffice for a 4-tuple (a1, a2, a3, a4) to impose unique relations upon k-th order transverse moments. An immediate corollary of these analytic arguments is the fact that any beam with rotational symmetry must have zero means (i.e. first order centroid moments) in all transverse phase space coordinates. For moments of order k = 1, when n (cid:44) 1, (a1, a2, a3, a4) = (1,0,0,0) and (a1, a2, a3, a4) = (0,0,1,0) always satisfy conditions 4.1.1, 4.1.2 and 4.1.3. Therefore, (cid:104)w(cid:105) = 0 and (cid:104)w(cid:48)(cid:105) = 0 must hold, which imply (cid:104)x(cid:105) = (cid:104)y(cid:105) = 0 and (cid:104)x(cid:48)(cid:105) = (cid:104)y(cid:48)(cid:105) = 0 respectively. 4.1.3 Correspondence between Discrete and Continuous Rotational Symmetry Armed with a general method for deriving relations imposed upon beam moments due to the rotational symmetry possessed by the beam, it is illuminating to ask: how do continuous and discrete rotational symmetries manifest themselves differently in beam moment relations? This question was part of the motivation for developing the analytic tools in Sec. 4.1.2, and will lead to arguments that enable simple proofs of certain useful results. We begin by counting the number of equality relations governing k-th order transverse moments of an axisymmetric beam. We define: • ξ(k) as the number of 4-tuples (a1, a2, a3, a4) that satisfy conditions 4.1.1; • χ(k) as the number of 4-tuples (a1, a2, a3, a4) that satisfy conditions 4.1.1 and 4.1.2; and • ˆχ(k) as the number of 4-tuples (a1, a2, a3, a4) that satisfy conditions 4.1.1 and (a1 + a3) − (a2 + a4) = 0. 39 ξ(k) is the number of distinct k-th order moments, which applies equally to real moments wa1 wa2 w(cid:48)a3 w(cid:48)a4(cid:69). It is equivalent to the total number of ways to distribute k indistinguishable balls into 4 distinguishable boxes, which gives (cid:104)xa1 x(cid:48)a2 ya3 y(cid:48)a4(cid:105) and complex moments(cid:68) ξ(k) =(cid:169)(cid:173)(cid:173)(cid:171)k + 3 (cid:170)(cid:174)(cid:174)(cid:172) = k (k + 3)! 3!k! = (k + 3)(k + 2)(k + 1) 3 × 2 × 1 (4.17) To find χ(k), we use the fact that the number of 4-tuples satisfying (a1 + a3) − (a2 + a4) > 0 and (a1 + a3) − (a2 + a4) < 0 are equal. Therefore, χ(k) = 1 2 [ξ(k) − ˆχ(k)] (4.18) For odd k, ˆχ(k) = 0 because(a1+a3)−(a2+a4) = 0 cannot occur. For even k,(a1+a3)−(a2+a4) = 0 occurs whenever a1 + a3 = a2 + a4 = k/2. The number of 2-tuples (a1, a3) such that a1 + a3 = k2 equals k2 + 1, likewise for (a2, a4). Since (a1, a3) and (a2, a4) are independent, ˆχ(k) = 0 equals . Summarizing these results with Eq. (4.17) and Eq. (4.18), we obtain (cid:16) k2 + 1(cid:17)2 χ(k) = (k+3)(k+2)(k+1) (cid:20)(k+3)(k+2)(k+1) 3×2×1 3×2×1 1 2 1 2 −(cid:16) k2 + 1(cid:17)2(cid:21) = k(k+2)(2k+5) 24 if k is odd, if k is even. (4.19)  Since each vanishing complex moment generates two relations, the total number of equality relations governing k-th order transverse moments of an axisymmetric beam is 2 χ(k). 2 χ(k) is the maximum number of moment relations that can be imposed upon k-th order moments of a beam due to rotational symmetry. For beams with discrete rotational symmetry Cn, the number of relations ≤ 2 χ(k) because some 4-tuples (a1, a2, a3, a4) which satisfy condition 4.1.1 and condition 4.1.2 may violate condition 4.1.3. Let m(n, k) be the number of 4-tuples (a1, a2, a3, a4) that satisfy conditions 4.1.1 and 4.1.2, but not condition 4.1.3, for the given n. Thus the number of moment relations governing a beam with Cn symmetry is given by 2[ χ(k) − m(n, k)]. We discuss the implications of m(n, k) = 0 and m(n, k) (cid:44) 0 below. We introduce the term “k-th order axisymmetry” to denote the set of 2 χ(k) relations obeyed by k-th order moments of an axisymmetric beam. A beam is said to be “k-th order axisymmetric” 40 if it obeys k-th order axisymmetry. Note that k-th order axisymmetry merely concerns properties of k-th order moments, so a beam can be k-th order axisymmetric but not l-th order axisymmetric for any l (cid:44) k. k-th order axisymmetry is a much weaker notion than axisymmetry, i.e. continuous rotational symmetry, which implies k-th order axisymmetry for all k. When m(n, k) = 0, a Cn beam is k-th order axisymmetric, i.e. its k-th order moments obey the same set of relations as those of an axisymmetric beam. This is true despite the fact that Cn and axisymmetry are different rotational symmetries. This equivalence indicates that k-th order axisymmetry renders the Cn beam indistinguishable from an axisymmetric beam in terms of k-th order beam moments. As a result, the Cn beam will share all properties of an axisymmetric beam that only depend on k-th order moments. For example,(cid:68) xk(cid:69) of an axisymmetric beam remains the same regardless of the orientation of the transverse coordinate system, so the same must hold for the Cn beam as well. m(n, k) = 0 is quite common and occurs in more than half of all (n, k) pairs. Referring again to conditions 4.1.2 and 4.1.3, to ensure that both of them are always satisfied, (a1 + a3)−(a2 + a4) = 0 (mod n) cannot occur unless (a1 + a3)−(a2 + a4) = 0. It is easy to show that this condition always holds if k < n, or if k is even and n is odd, or vice versa. A corollary of the above statement is that, given n, k = n is the smallest k at which m(n, k) (cid:44) 0 occurs. This agrees with the heuristic picture where, the higher the degree of discrete rotational symmetry, the closer the beam is to being axisymmetric and the higher the order of moments it takes for deviations to start. Let S(Cn) and S(axisymmetry) be the sets of unsimplified relations derived from Eq. (4.11) governing n-th order moments of a Cn beam and an axisymmetric beam respectively. When k = n, m(n, n) = n + 1, so S(axisymmetry) contains 2(n + 1) more elements (i.e. moment relations) than S(Cn). With n being the lowest moment order at which Cn and axisymmetry are distinguishable, the difference between S(Cn) and S(axisymmetry) can be interpreted as a manifestation of the difference between Cn and axisymmetry in terms of beam moments. Such a perspective can allow one to employ beam moments to quantify how close a Cn beam is from being axisymmetric. The quantification may be useful if one would like to gauge the impact of 41 elements with Cn symmetry on the beam distribution, or if one would like to generate a Cn beam via maximum entropy methods with prior knowledge on beam moments. The general arguments developed in this subsection will be applied to analyze the 2nd and 3rd order moments of a beam with C3 symmetry in Sec. 4.1.4. 4.1.4 Example: Beams with C3 Symmetry To demonstrate the utility of the analytic tools developed in Sec. 4.1.1 and 4.1.2, we derive equality relations governing the moments of a beam with C3 symmetry, and compare them against those of an axisymmetric beam. Beams with C3 symmetry are directly relevant to ECR ion sources which contain a sextupole (C3) in additional to solenoids (axisymmetric). Simulation results of an ECR beam at the extraction plane is shown in Fig. 4.1 where the spatial distribution of ions extracted is visibly triangulated due to the sextupole field of the ECR. Note that the dominant triangulation is shown over the plasma chamber extent but is still manifest within the extraction aperture. 4.1.4.1 2nd Order Moments For k = 2, the complex moments satisfying conditions 4.1.1 to 4.1.3 can easily be found to be (cid:104)ww(cid:105), (cid:104)ww(cid:48)(cid:105) and (cid:104)w(cid:48)w(cid:48)(cid:105). The equality relations they impose on transverse beam moments are given as follows: (cid:10)ww (cid:10)w w (cid:48) (cid:104)ww(cid:105) = 0 (cid:48)(cid:11) = 0 (cid:48)(cid:11) = 0 (cid:104)x y(cid:105) = 0 (cid:104)x y (cid:104)x(cid:48) (cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:48)(cid:105) = 0 y y(cid:105) ⇒ ⇒ ⇒ (cid:104)xx(cid:105) = (cid:104)yy(cid:105) (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:48)(cid:105) (cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:48) y 42 Figure 4.1: Simulated spatial distribution of Ar9+ ions at the extraction plane of the Artemis [5] ECR ion source at FRIB. The inner circle indicates the extraction aperture. x-y and u-v denote two possible transverse coordinate systems. Subsequent analysis will show that, despite severe triangulation, the rms envelope and emittance of the beam is the same along any transverse direction. Simulation results courtesy of the code developed by Dr. Vladirmir Mironov at JINR-Dubna [25, 26]. Together these constraints give: (cid:104)xx(cid:105) = (cid:104)yy(cid:105) (cid:48)(cid:105) (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:48) (cid:48)(cid:105) (cid:104)x y(cid:105) = 0 (cid:104)x(cid:48) (cid:48)(cid:105) = 0 y (cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:104)x y y y(cid:105) (4.20) (4.21) (4.22) (4.23) (4.24) (4.25) An alternative proof of the fact that Eq. (4.20) to (4.25) holds for a beam with 3-fold rotational symmetry is provided in Appendix ??. That proof utilizes trigonometric identities and is applicable to any n-fold rotational symmetry where n > 2. Eq. (4.20) to (4.25) are the same set of relations that are obeyed by a beam with continuous 43 rotation symmetry. Therefore, following the argument in Sec. ??, we can immediately conclude that, not only are the 2nd-order beam moments in x-x(cid:48) and y-y(cid:48) identical, they must be identical in any transverse direction. This can also be proved by writing out moments in an arbitrary direction with u = x cos θ + y sin θ, but the approach is much more tedious than the indistinguishability argument and is not easily extendable to higher order moments. The results also imply that the azimuthal orientation of the ECR sextupole (see Fig. 4.1) has no effect on 2nd-order moments. Small errors (e.g. mechanical misalignment) may introduce a limited degree of violation to this useful idealization. 3rd Order Moments 4.1.4.2 For k = 3, we use Eq. (4.19) to obtain χ(3) = 10, which means 10 complex moments satisfy conditions 4.1.1 and 4.1.2. We also know that, since the moment order k = 3 is the same as n, 4 complex moments must violate condition 4.1.3, namely: (cid:104)www(cid:105) (cid:10)www (cid:10)ww (cid:10)w w w w (cid:48) (cid:48) (cid:48) (cid:48)(cid:11) (cid:48)(cid:11) (cid:48)(cid:11) They provide no information because e3iθ = 1 for θ = 2π/3. The 6 complex moments which also satisfy condition 4.1.3 give: (cid:104)www(cid:105) = 0 ⇒ = 0 ⇒ (cid:68) www(cid:48)(cid:69) (cid:48)(cid:11) = 0 ⇒ (cid:10)www (cid:10)ww (cid:48)(cid:11) = 0 ⇒ w(cid:48)(cid:69) (cid:68) w(cid:48)(cid:69) (cid:68) (cid:104)xxx(cid:105) = −(cid:104)x yy(cid:105) (cid:48)(cid:105) = (cid:104)x(cid:48) 2(cid:104)x yy (cid:104)xxx(cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:104)xx(cid:48)x(cid:48)(cid:105) = −(cid:104)x y = 0 ⇒ 2(cid:104)x(cid:48) (cid:48)(cid:105) = (cid:104)x y (cid:48) = 0 ⇒ (cid:104)x(cid:48)x(cid:48)x(cid:48)(cid:105) = −(cid:104)x(cid:48) yy(cid:105) − (cid:104)xxx(cid:48)(cid:105) yy(cid:105) (cid:48) (cid:48)(cid:105) (cid:48)(cid:105) − (cid:104)xx(cid:48)x(cid:48)(cid:105) (cid:48)(cid:105) (cid:48) ww (cid:48) yy w w w (cid:48) (cid:48) (cid:48) y y y y (cid:48)(cid:105) (cid:104)yyy(cid:105) = −(cid:104)xx y(cid:105) 2(cid:104)xx(cid:48) (cid:104)xx y (cid:104)x(cid:48)x(cid:48) 2(cid:104)x(cid:48)x y (cid:48) (cid:104)y (cid:48)(cid:105) − (cid:104)yyy y(cid:105) = (cid:104)xx y (cid:48)(cid:105) = −(cid:104)yyy (cid:48)(cid:105) (cid:48) (cid:48)(cid:105) y(cid:105) = −(cid:104)yy y (cid:48)(cid:105) = (cid:104)x(cid:48)x(cid:48) y(cid:105) − (cid:104)yy (cid:48)(cid:105) = −(cid:104)x(cid:48)x(cid:48) (cid:48)(cid:105) (cid:48) y y y (cid:48) (cid:48)(cid:105) y 44 Further combination of these relations show that they impose the following independent conditions on third order moments: (cid:104)xxx(cid:105) = −(cid:104)x yy(cid:105) (cid:104)yyy(cid:105) = −(cid:104)xx y(cid:105) (cid:48)(cid:105) = (cid:104)x(cid:48) yy(cid:105) = −(cid:104)xxx(cid:48)(cid:105) (cid:104)x yy (cid:104)xx(cid:48) (cid:48)(cid:105) = −(cid:104)yyy (cid:48)(cid:105) y(cid:105) = (cid:104)xx y (cid:104)x(cid:48) (cid:48)(cid:105) = (cid:104)x y (cid:48) (cid:48)(cid:105) = −(cid:104)xx(cid:48)x(cid:48)(cid:105) y (cid:104)xx(cid:48) (cid:48)(cid:105) = (cid:104)x(cid:48)x(cid:48) (cid:48)(cid:105) y(cid:105) = −(cid:104)yy (cid:104)x(cid:48)x(cid:48)x(cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:48) (cid:48) (cid:48)(cid:105) = −(cid:104)x(cid:48)x(cid:48) (cid:104)y (cid:48)(cid:105) (cid:48)(cid:105) yy y y y y (cid:48) y (cid:48) y y Conversely, all third order moments of an axisymmetric beam must vanish. The same conclusion applies to all odd-order moments of an axisymmetric beam as well. This can be seen immediately if one applies a rotation by θ = π. This difference in moment relations between a C3 beam and an axisymmetric beam yields interesting implications for beam transport. Whereas the approximation of linear optics is often made, beam transport inevitably involve nonlinear terms in the transfer map. Consider a beam with some specific rotational symmetry entering a beam line with the nonlinear map: 4 i=1 4 4 i=1 j=1 xl = Rli xi + i ) Tli j xi xj + O(x3 (4.26) where subscripts can all attain integer values from 1 to 4 with l, m denoting final phase space coordinates and i, j, k denoting initial phase space coordinates. Second order moments of the beam exiting the beam line are given by: 4 4 i=1 j=1 (cid:10)xi xj (cid:11) + 4 4 4 (cid:16) i=1 j=1 k=1 (cid:104)xl xm(cid:105) = RliRm j RliTm jk + RmiTl jk (cid:17)(cid:10)xi xj xk (cid:11) + O(cid:16)(cid:68) (cid:69)(cid:17) x4 i (4.27) 45 which have a dependence on initial moments higher than second order due to nonlinear terms in the map. For an axisymmetric beam, all of its 3rd order moments vanish, so the leading correction to linear optics arises from 4th order moments. For a C3 beam, however, 3rd order moments are probably nonzero and can affect the final second order moments. Hence, linear optics is a better approximation in the sense that residual terms neglected are high order when the incoming beam is axisymmetric than when it has C3 symmetry. Last but not least, we note that all 3rd order moments of a C2 beam also vanish. Thus the discussion in the paragraph above also holds if we substitute axisymmetry by C2, the most commonly occurring rotational symmetry in beam lines due to the ubiquity of quadrupoles. This fact may render linear optics a surprisingly poor assumption when it is applied to C3 beams. The assumption may be poor because second order terms in the transfer map act more strongly on the rms evolution of C3 beams than on that of axisymmetric or C2 beams, and the poorness may be surprising because most of our experience is built upon the less error-prone cases of C2 and axisymmetry. 4.2 Coupled Beam Dynamics This section presents several results on linear coupled beam dynamics in single particle dynam- ics (i.e. self fields are neglected). 4.2.1 Axisymmetric Beams and Eigen-emittances This subsection derives the relation between εx, (cid:104)x y(cid:48)(cid:105) and the eigen-emittances for an axisymmetric beam. An axisymmetric beam has sigma matrix with the following (generally non-zero) 2nd-order moment elements: (4.28) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) σ = (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:104)xx(cid:105) (cid:104)xx(cid:48)(cid:105) 0 (cid:104)x y(cid:48)(cid:105) 0 (cid:104)xx(cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) − (cid:104)x y(cid:48)(cid:105) − (cid:104)x y(cid:48)(cid:105) (cid:104)xx(cid:105) (cid:104)xx(cid:48)(cid:105) 0 (cid:104)x y(cid:48)(cid:105) 0 (cid:104)xx(cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) 46 From the eigenvalues of σS, one obtains the two eigen-emittances: ε1 = εx +(cid:10)x y ε2 = εx −(cid:10)x y (cid:48)(cid:11) (cid:48)(cid:11) where εx is the usual rms emittance defined by Eq. (2.19). It follows that: (cid:10)x y εx = (cid:48)(cid:11) = ε1 + ε2 (cid:16) ε1 − ε2 2 2 (cid:17) (4.29) (4.30) (4.31) (4.32) 4.2.2 Magnetized Beam Consider an axisymmetric beam in an axisymmetric beam line going from si to s f where Bz0(si) (cid:44) 0 and Bz0(s f ) = 0 At s = si, due to the presence of axial magnetic field, canonical coordinates ˜x and laboratory coordinates x are related by (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) K = 0 ˜x = Kx 0 0 1 1 −ks/2 0 0 0 0 0 ks/2 0 1 1 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) where ks = qB mv The sigma matrix written in terms of canonical coordinates is given by σi = ˜xi˜xT i = KxixT i KT (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) σi = 0 0 0 1 1 −ks/2 0 0 0 0 0 ks/2 0 1 1 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)xx(cid:105)i (cid:104)xx(cid:48)(cid:105)i 0 (cid:104)x y(cid:48)(cid:105)i (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:48)(cid:105)i 0 (cid:104)xx(cid:48)(cid:105)i (cid:104)x(cid:48)x(cid:48)(cid:105)i 0 1 0 ks/2 1 0 0 0 0 −ks/2 1 0 0 0 1 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 0 (cid:104)xx(cid:48)(cid:105)i (cid:104)x(cid:48)x(cid:48)(cid:105)i − (cid:104)x y(cid:48)(cid:105)i − (cid:104)x y(cid:48)(cid:105)i (cid:104)xx(cid:105)i (cid:104)xx(cid:48)(cid:105)i 0 47 At s = s f , σ f = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)xx(cid:105) f (cid:104)xx(cid:48)(cid:105) f 0 (cid:104)x y(cid:48)(cid:105) f (cid:104)xx(cid:48)(cid:105) f (cid:104)x(cid:48)x(cid:48)(cid:105) f − (cid:104)x y(cid:48)(cid:105) f 0 0 − (cid:104)x y(cid:48)(cid:105) f (cid:104)xx(cid:105) f (cid:104)xx(cid:48)(cid:105) f (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:104)x y(cid:48)(cid:105) f 0 (cid:104)xx(cid:48)(cid:105) f (cid:104)x(cid:48)x(cid:48)(cid:105) f (cid:10)x y (cid:48)(cid:11) f =(cid:10)x y (cid:48)(cid:11) ks 2 (cid:104)xx(cid:105)i + i From conservation of canonical angular momentum: (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) σ f = (cid:104)xx(cid:105) f (cid:104)xx(cid:48)(cid:105) f 0 (cid:104)x y(cid:48)(cid:105)i + ks2 (cid:104)xx(cid:105)i (cid:104)xx(cid:48)(cid:105) f (cid:104)x(cid:48)x(cid:48)(cid:105) f − (cid:104)x y(cid:48)(cid:105)i − ks2 (cid:104)xx(cid:105)i 0 0 − (cid:104)x y(cid:48)(cid:105)i − ks2 (cid:104)xx(cid:105)i (cid:104)xx(cid:105) f (cid:104)xx(cid:48)(cid:105) f (cid:104)x y(cid:48)(cid:105)i + ks2 (cid:104)xx(cid:105)i 0 (cid:104)xx(cid:48)(cid:105) f (cid:104)x(cid:48)x(cid:48)(cid:105) f We can solve for eigenvalues of σiS and σ f S to get eigen-emittances: (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 2 i; I,II =  2 i,x + ±(cid:16) k2 2 (cid:104)xx(cid:105)2 s ks (cid:104)xx(cid:105)i + 2(cid:10)x y (cid:48)(cid:11) i + 2ks (cid:104)xx(cid:105)i (cid:10)x y (cid:48)(cid:11) +(cid:10)x y (cid:48)(cid:11)2 i i + k2 i (cid:17)(cid:113) (cid:10)x y 2  i,x (cid:48)(cid:11) i s (cid:104)xx(cid:105)i /4 + ks (cid:104)xx(cid:105)i (cid:104)x y(cid:48)(cid:105)i +(cid:10)x y (cid:48)(cid:11)2 i ±  f ,x (cid:16) ks (cid:104)xx(cid:105)i + 2(cid:10)x y (cid:48)(cid:11) (cid:17) i 2 f ; I,II =  2 f ,x  + k2 4 (cid:104)xx(cid:105)2 s i + ks (cid:104)xx(cid:105)i Eigen-emittances are invariants, i.e. Therefore, the final x-emittance is given by: i; I,II =  f ; I,II 2 f ,x  =  2 i,x + k2 4 (cid:104)xx(cid:105)2 s i + ks (cid:104)xx(cid:105)i 48 (cid:10)x y (cid:48)(cid:11) i 4.2.3 Mirror-Symmetric Beam Line Does Not Preserve Axial-symmetry In a beam line consisting of quads and dipoles only, equations of motion in x and y are decoupled. Thus the transfer matrix is block diagonal, with the general form: (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) M = m11 m12 0 0 0 0 0 m33 m34 0 m43 m44 m21 m22 0 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) Suppose there is a symmetric beam line, i.e. M(s f |si) = M1M2 . . . Mn−1MnMnMn−1 . . . M2M1 where M(sm|si) = MnMn−1 . . . M2M1 M(s f |sm) = M1M2 . . . Mn−1Mn (cid:170)(cid:174)(cid:174)(cid:172) Proof: m21 m11 m21 m22 Where si, sm and s f stand for the start, mid- and end point respectively. Lemma: if Mx(sm|si) =(cid:169)(cid:173)(cid:173)(cid:171)m11 m12 (cid:170)(cid:174)(cid:174)(cid:172) then Mx(s f |sm) =(cid:169)(cid:173)(cid:173)(cid:171)m22 m12 (cid:169)(cid:173)(cid:173)(cid:171) x = Mx(sm|si)(cid:169)(cid:173)(cid:173)(cid:171) x (cid:170)(cid:174)(cid:174)(cid:172)m (cid:170)(cid:174)(cid:174)(cid:172)i =(cid:169)(cid:173)(cid:173)(cid:171) x Mx(s f |sm)(cid:169)(cid:173)(cid:173)(cid:171) x (cid:170)(cid:174)(cid:174)(cid:172)m By time-reversal: x(cid:48) (cid:170)(cid:174)(cid:174)(cid:172)i −x(cid:48) x(cid:48) −x(cid:48) 49 where (cid:170)(cid:174)(cid:174)(cid:172) 0 0 −1 0 0 −1 m21 m22 m21 m11 0 0 −1 x(cid:48) 0 0 −1 0 0 −1 (cid:169)(cid:173)(cid:173)(cid:171) x x(cid:48) =(cid:169)(cid:173)(cid:173)(cid:171)1 (cid:169)(cid:173)(cid:173)(cid:171) x Mx(s f |sm)(cid:169)(cid:173)(cid:173)(cid:171)1 (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172)m x (sm|si)(cid:169)(cid:173)(cid:173)(cid:171)1 Mx(s f |sm) =(cid:169)(cid:173)(cid:173)(cid:171)1 (cid:170)(cid:174)(cid:174)(cid:172) M−1 T =(cid:169)(cid:173)(cid:173)(cid:171)1 (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) =⇒(cid:101)Mx =(cid:169)(cid:173)(cid:173)(cid:171)m22 m12 Mx =(cid:169)(cid:173)(cid:173)(cid:171)m11 m12 My =(cid:169)(cid:173)(cid:173)(cid:171)m33 m34 (cid:170)(cid:174)(cid:174)(cid:172) =⇒(cid:101)My =(cid:169)(cid:173)(cid:173)(cid:171)m44 m34 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)m (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)i (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 0 0 1 0 −1 0 0 0 0 1 0 −1 0 (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) 0 0 0 1 0 −1 0 0 0 0 1 0 −1 0 0 0 = Mx ⊗ My m43 m44 m43 m33 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) M−1 0 0 0 (cid:101)M = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) f (cid:170)(cid:174)(cid:174)(cid:172)i (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172)i (cid:16)(cid:101)MxMx (cid:17) ⊗(cid:16)(cid:101)MyMy = (cid:17)(cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) 50 For (cid:104)x y(cid:105)i = 0, (cid:104)x(cid:48)y(cid:48)(cid:105)i = 0 and (cid:104)x y(cid:48)(cid:105)i = −(cid:104)x(cid:48)y(cid:105)i: (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) f (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) = 2(m11m22 + m12m21)m34m44 − 2(m33m44 + m34m43)m12m22 −4m12m22m33m43 + (m11m22 + m12m21)(m33m44 + m34m43) 4m11m21m34m44 − (m11m22 + m12m21)(m33m44 + m34m43) −2(m11m22 + m12m21)m33m43 + 2(m33m44 − m34m43)m11m21 (cid:104)x y (cid:48)(cid:105)i which do not obey any of (cid:104)x y(cid:105) f = 0, (cid:104)x(cid:48)y(cid:48)(cid:105) f = 0 or (cid:104)x y(cid:48)(cid:105) f = −(cid:104)x(cid:48)y(cid:105) f in general. In fact, (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) f (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 0 L/2 −L/2 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = only if Mx = My 4.2.4 Linear Solenoid Transport Section 4.1 proves that the 2nd-order moments of an axisymmetric beam must obey the following relations: (cid:104)xx(cid:105) = (cid:104)yy(cid:105) (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:48)(cid:105) (cid:48) (cid:104)x y(cid:105) = 0 (cid:48)(cid:105) = 0 (cid:104)x(cid:48) y (cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:104)x y y y(cid:105) We call a beam that fulfills these relations 2nd-order axisymmetric, which is a much weaker notion than axisymmetric because it does not impose any requirements on higher order moments. This 51 if a subsection shows that 2nd-order axisymmetry is preserved by linear solenoid transport, i.e. 2nd-order axisymmetric beam enters a linear solenoid, it will retain its 2nd-order axisymmetry upon exit. Note that in this context we are neglecting space charge forces which, depending on the distribution, may exert non-axixymmetric self-fields. An axisymmetric beam will still be axisymmetric after passing through a solenoid, which means that it obeys 2nd-order axisymmetry both upstream and downstream. Furthermore, when the solenoid is linear, the final beam moments are only functions of initial beam moments of the same order. Therefore, a linear solenoid has a transfer map that preserves 2nd-order axisymmetry. This is a general property of the transfer maps of linear solenoids that is agnostic to beam information other than 2nd-order moments. So as long as a beam has 2nd-order axisymmetry, it will be preserved by a linear solenoid. Lastly, we note that, for a linear solenoid, not only does 2nd-order axisymmetric upstream imply 2nd-order axisymmetric downstream, considerations from time reversal also dictates 2nd- order axisymmetric downstream implies 2nd-order axisymmetric upstream. 4.2.5 Quasi-axisymmetric Conditions As shown in Sec. 4.2.4, linear solenoid optics preserves 2nd-order axisymmetry, so it would be ideal for a beam to obey the corresponding relations when it enters solenoid transport. In this subsection, we define a weaker set of conditions which we call quasi-axisymmetry. A non-axisymmetric beam satisfying its conditions have properties similar to that of an axisymmetric beam in solenoid transport. The conditions are weaker because it does not require εx = εy, which is a necessary condition for an axisymmetric beam. This has practical importance when a beam with εx (cid:44) εy has to be matched into a solenoid transport line using only normal quadrupoles. 52 The quasi-axisymmetry conditions read: αx = αy βx = βy (cid:104)x y(cid:105) = 0 (cid:48)(cid:11) = 0 (cid:48)(cid:11) = −(cid:10)x(cid:48) (cid:10)x(cid:48) (cid:10)x y y (4.33) (4.34) (4.35) (4.36) (4.37) y(cid:11) In terms of 2D phase space projections, the first two conditions describe rms ellipses in x-x(cid:48) and y-y(cid:48) phase spaces which have the same orientation and sizes of the ratio εx : εy. Therefore, (cid:104)xx(cid:105) = (cid:10)xx(cid:48)(cid:11) = (cid:10)x(cid:48)x(cid:48)(cid:11) = εy εx εy εx εy εx (cid:104)yy(cid:105) (cid:10)yy (cid:10)y (cid:48)(cid:11) (cid:48)(cid:11) (cid:48) y (4.38) (4.39) (4.40) If we express moments of a quasi-axisymmetric beam in an arbitrary transverse direction in terms of the moments in x, x(cid:48) and y, y(cid:48), we obtain, for u = x cos θ + y sin θ: (cid:104)uu(cid:105) = (cid:104)xx(cid:105) cos2 = (cid:104)xx(cid:105) cos2 = (cid:104)xx(cid:105) cos2 = (cid:104)xx(cid:105) (cid:20) 1 + θ + 2 (cid:104)x y(cid:105) sin θ cos θ + (cid:104)yy(cid:105) sin2 θ + (cid:104)yy(cid:105) sin2 θ + (cid:104)xx(cid:105) εy εx sin2 (cid:18) εy θ sin2 (cid:19) (cid:21) − 1 θ θ εx θ (4.41) Applying Eq. (4.38), one can check that Eq. (4.41) correctly recovers (cid:104)uu(cid:105) = (cid:104)xx(cid:105) and (cid:104)uu(cid:105) = (cid:104)yy(cid:105) at θ = 0 and θ = π/2 respectively. Similarly, (cid:10)uu(cid:48)(cid:11) =(cid:10)xx(cid:48)(cid:11)(cid:20) (cid:10)u(cid:48)u(cid:48)(cid:11) =(cid:10)x(cid:48)x(cid:48)(cid:11)(cid:20) (cid:19) (cid:19) (cid:18) εy (cid:18) εy εx εx 1 + 53 (cid:21) (cid:21) 1 + − 1 sin2 θ − 1 sin2 θ (4.42) (4.43) The Twiss parameters and emittance along u are given by: (cid:113)(cid:104)uu(cid:105) (cid:104)u(cid:48)u(cid:48)(cid:105) − (cid:104)uu(cid:48)(cid:105)2 Applying Eq. (4.41) to Eq. (4.43), it can be shown that: εu = αu = βu = − (cid:104)uu(cid:48)(cid:105) εu (cid:104)uu(cid:105) εu εu = εx (cid:18) εy (cid:20) (cid:19) = εx +(cid:0)εy − εx(cid:1) sin2 1 + εx − 1 sin2 (cid:21) θ θ = αx − (cid:104)uu(cid:48)(cid:105) εu (cid:104)uu(cid:105) εu = αu = βu = − (cid:104)xx(cid:48)(cid:105) = (cid:104)xx(cid:105) εx εx = βx (4.44) (4.45) (4.46) (4.47) (4.48) (4.49) Eq. (4.48) and Eq. (4.49) states that αu and βu has no dependence on θ and always equal their counterparts in x. This means that, for a quasi-axisymmetric beam, despite εx (cid:44) εy, the Twiss parameters are the same in all transverse directions. On the other hand, Eq. (4.47) shows that εu changes with θ. εu attains its maximum and minimum value along x and y respectively if εx > εy, and vice versa. 4.3 Applications to FRIB Front End This section shows how the results developed in this chapter can be applied to the FRIB front end. Fig. 4.2 shows a schematic of the FRIB front end where only elements relevant to transverse dynamics are displayed. The beam line is also divided into segments (1) through (7) which will play different roles in beam tuning. Section 4.3.1 discusses how results from Sec. 4.1 clarify the beam dynamics in segments 1 and 2 and motivate investigation on an unresolved phenomenon. Section 4.3.2 employs arguments in Sec. 4.2 to elucidate attainable beam states in segments 3 to 6 and presents a tuning scheme that can be readily implemented in the laboratory for achieving optimal matching into the RFQ. Section 4.3.3 analyzes the intertwined problems of how the optics 54 Figure 4.2: A schematic of transverse focusing and bending elements in the FRIB front end. Note that Segment 2 and 5 are mirror-symmetric about the axial mid-plane. Lengths are not drawn to scale and all quadrupoles are normal quadrupoles (i.e. upright orientation). in Linac Segment 1 (following segment 7) should be designed and how the MEBT (segment 7) should match a non-axisymmetric beam into Linac Segment 1. Space charge is relatively weak after the mid-plane of the CSS (segment 2), hence it is neglected in the discussion on tuning in segments 3 to 7. 4.3.1 LEBT: Before Species Selection Using the analytic tools developed Sec. 4.1, we proved in Sec. 4.1.4 that a beam with C3 symmetry has the same emittance in any transverse direction. This result should hold in the presence of multi-species space charge, chromatic aberrations and radial field nonlinearities because none of these effects break the discrete rotational symmetry of the system. This conclusion quickly revealed a conundrum at FRIB because diagnostic measurements often show >10% difference between εx and εy at the mid-point and end-point of the charge selection system (i.e. segment 2). Prior to this analysis, the difference used to be attributed to the sextupole 55 of the ECR ion source that possibly results in a triangulated beam, symmetry arguments from Sec. 4.1 invalidated such an explanation and gave impetus to a more thorough inquiry. We have identified two possible reasons for εx (cid:44) εy to occur in the charge selection system, neither of which we understand fully: (A) the beam has εx (cid:44) εy entering segment 2; and (B) emittance evolution in the dipole or prior to charge separation. Further investigation is underway on both fronts to explain this phenomenon. If the possibility A is true, the results in Sec. 4.1 compels us to argue by transposition that axisymmetry is broken in segment 1. Thus we searched for plausible causes and concluded that an old solenoid with poorly designed current leads may create strong multipole moments in the downstream fringe fields. Possibility B may also be true because εx (cid:44) εy was observed in Warp simulation results in Sec. 3.3 due to an interplay between space charge and beam magnetization, albeit with a smaller percentage difference. Furthermore, the discussion at the end of Sec. 4.1 shows that the evolution of second order moments of a beam with C3 symmetry is likely more susceptible to higher order terms in the transfer map than that of an axisymmetric beam. Since the transfer map of the 90-degree combined function dipole may contain significant nonlinearities, the use of axisymmetric beams in previous simulations could induce large errors in modeling the dipole transport. This suggests that further simulation studies should include non-axisymmetric initial conditions with C3 symmetry to probe how this changes the beam dynamics in segment 2. 4.3.2 LEBT: After Species Selection After charge selection, the LEBT has to match the beam into the RFQ while avoiding beam loss. The design of the optics can be done in two steps: (A) obtain complete information on transverse coupled beam moments after segment 2 via diagnostic measurements; and (B) tune quadrupoles and solenoids via optimization to accomplish design goals. There are a total of 17 knobs available for tuning the beam, they comprises the focusing strengths of 10 quadrupoles and 7 solenoids. The vertical drop only has six knobs rather than twelve because only two of the four quadrupole triplets 56 are independent in its design as a mirror-symmetric achromat. Design goals include small beam sizes at locations with small apertures (e.g. within the buncher and chopper), as well as: αx = αy = αin βx = βy = βin εx = εy (4.50) (4.51) (4.52) at the RFQ entrance where αin and βin denote the design Twiss parameters entering the RFQ. The most direct scheme is to employ all 17 knobs as optimization parameters and construct a penalty function that incorporates all design goals as constraints. This approach was adopted by the Accelerator Physics Department at FRIB. An example of optics design 1 based on the direct scheme is plotted in Fig. 4.3. Note that there is strong emittance exchange in the lower LEBT solenoid transport (segment 6) which causes the beam to enter the RFQ with large emittance differences in x and y. Such emittance exchange is a common problem that occurs because the beam does not fulfill 2nd order axisymmetry in the solenoid transport line. We designed an easily implementable tuning scheme based on results in coupled beam dynamics developed in Sec. 4.2. This treatment delineates attainable and preservable beam states in each beam line segment, which allow us to split the optimization into separate stages with well-defined goals. In the ensuing analysis, it is important to distinguish between beam line segments that: (A) con- tain solenoids, which couple the dynamics in x and y; and (B) contain dipoles and/or quadrupoles, which does not couple the couple the dynamics in x and y. It is proved in Sec. 4.2.4 that linear 1provided by Dr. Tomofumi Maruta 57 Figure 4.3: An example of FRIB LEBT tuning results where the focusing strengths of 17 elements (10 quadrupoles and 7 solenoids) were optimized to achieve desired beam conditions at the RFQ entrance. (Image courtesy of Dr. Tomofumi Maruta at FRIB Accelerator Physics Department) solenoids preserve the set of 2nd-order axisymmetric conditions: (cid:104)xx(cid:105) = (cid:104)yy(cid:105) (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:48)(cid:105) (cid:48)(cid:105) (cid:48) (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:104)x y(cid:105) = 0 (cid:104)x(cid:48) (cid:48)(cid:105) = 0 y (cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:104)x y y y(cid:105) For dipoles and quadrupoles, εx and εy are conserved quantities under linear optics (see Chapter 2). We extract from these results three conditional statements on beam states and beam transport: Corollary 4.3.1. (2nd-order axisymmetry) → (εx = εy) 58 "Catchall" Optimization Corollary 4.3.2. In linear optics, for a beam line consisting of solenoids and drifts, (incoming beam is 2nd-order axisymmetric) ↔ (outgoing beam is 2nd-order axisymmetric) Corollary 4.3.3. In linear optics, for a beam line consisting of normal dipoles, quadrupoles and drifts, (incoming beam has εx = εy) ↔ (outgoing beam has εx = εy) We start at the end (RFQ entrance) and work backwards. To achieve Eq. (4.50) to (4.52) at the RFQ entrance, it is highly desirable for the beam to be 2nd-order axisymmetric in the lower LEBT. 2nd-order axisymmetry guarantees: αx = αy βx = βy εx = εy throughout the solenoid transport, which means the beam optics in identical in any transverse direction and reduces to an effective 1D (radial) design. Hence, matching into the RFQ entrance with Twiss parameters αin and βin becomes a problem of tuning 4 solenoids to meet 2 constraints. Not only it the goal easily attainable, the problem also affords additional freedom to accommodate other constraints, such as small beam size at the multi-harmonic buncher. The next step concerns how 2nd-order axisymmetry can be achieved in the lower LEBT. A summary of the problem can be found in Table 4.1. To have 2nd-order axisymmetry after segments 5 and 6 implies having εx = εy as well (see Corollary 4.3.1), and Corollary 4.3.3 in turn implies εx = εy before the vertical drop. Concerning the conditions after the charge selection system, the results in Sec. 4.2.3 show that it takes Mx = My for an uncoupled mirror symmetry beam line to preserve 2nd-order axisymmetry. With only 3 quadrupoles on each half, segment 2 cannot fulfill this condition even if the beam is 2nd-order axisymmetric before segment 2. Therefore, the beam generally lacks 2nd-order axisymmetry after segment 2. In practice, as we noted in Sec. 4.3.1, measurements show that the beam almost always fail to even fulfill εx = εy after the charge selection system. Remaining uncertainties in the design are denoted by question marks and will be resolved below. 59 After segment 2 (charge selection system) After segment 3 (upper LEBT quad.) After segment 4 (upper LEBT sol.) After segment 5 (vertical drop) After segment 6 (lower LEBT) εx = εy X ? (cid:88) (cid:88) (cid:88) 2nd-order axisymmetric X ? ? (cid:88) (cid:88) Table 4.1: Beam states entailed by the requirement of 2nd-order axisymmetry in the lower LEBT. (cid:88) and X denote true and false respectively. Two possible schemes to tune the LEBT are shown in Table 4.2. We proposed Scheme 2 which tunes for εx = εy in the upper LEBT and uses the vertical drop to achieve 2nd-order axisymmetry. Scheme 1, which tunes for 2nd-order axisymmetry in the upper LEBT, is also discussed because it was an adaptation of the FRIB baseline design. Arguments below will illustrate why Scheme 2 is preferable. Scheme 1 is unlikely to work for the following reasons. Firstly, a prerequisite of the scheme is εx = εy after the charge selection system, otherwise εx = εy cannot hold after segment 3 either (by corollary 4.3.3). However, measurement results have shown that this prerequisite almost always fails to hold. Furthermore, even if εx = εy is true after segment 2, tuning for 2nd-order axisymmetry after segment 3 is difficult to succeed. In general, after εx = εy removes one condition, the beam must satisfy five more conditions to attain 2nd-order axisymmetry. However, segment 3 only contains four quadrupoles, so the task employs a four-parameter optimization to meet five conditions, which is unlikely to produce favorable operating points even if we disregard aperture limits and other secondary considerations. Before explaining how Scheme 2 works, we can show by deduction that it is the only possible arrangement. Since εx (cid:44) εy after segment 2, it entails (by Corollary 4.3.3) εx (cid:44) εy after segment 3 , which in turn implies that 2nd-order axisymmetry cannot hold after segment 3 (by Corollary 4.3.1) and, consequently, after segment 4 (by Corollary 4.3.2). The tuning in Scheme 2 is performed in three stages: 1. Stage 1 (segment 3 and 4): attain εx = εy at the entrance of segment 5. Seven knobs to tune for one condition. 60 After segment 2 (charge selection system) After segment 3 (upper LEBT quad.) After segment 4 (upper LEBT sol.) After segment 5 (vertical drop) After segment 6 (lower LEBT) εx = εy (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) X (cid:88) (cid:88) (cid:88) (cid:88) εx = εy X X (cid:88) (cid:88) (cid:88) Scheme 1 2nd-order axisymmetric Scheme 2 2nd-order axisymmetric X X X (cid:88) (cid:88) Table 4.2: First set of quadrupole scan parameters 2. Stage 2 (segment 5): attain 2nd-order axisymmetry at the entrance of segment 6. Six knobs to tune for five conditions. 3. Stage 3 (segment 6): match beam into RFQ design parameters. Four knobs to tune for two conditions. Each optimization contains more parameters than constraints, so it is likely that a favorable solution will be found even with secondary constraints such as a small beam size at the multi-harmonic buncher upstream of the RFQ. Fig. 4.4 shows an example optics design based on Scheme 2 of Table 4.2. The optimization was implemented by Dr. Tomofumi Maruta using a linear optics envelope code FLAME[29]. While emittance exchange occurs in the solenoids of segment 3, as a result of the optimization, the solenoid strengths are arranged such that εx = εy when the exchange ends. There exists small emittance exchange in the lower LEBT because the optimization does not attain perfect 2nd-order axisymmetry after the vertical drop, but the exchange is much smaller than the direct optimization design in Fig. 4.3, and εx is very close to εy at the RFQ entrance. This scheme is presently employed for front end tuning in FRIB and has superseded previous methods. 4.3.3 MEBT and LS1 Despite tuning efforts in the LEBT, the beam often has different emittances in x and y (εx (cid:44) εy) downstream of the RFQ. With only normal quadrupoles in the MEBT, it is impossible to achieve 61 Figure 4.4: An example of FRIB LEBT tuning results where the tuning scheme introduced in Sec. 4.3.2 was applied to optimize elements strengths in three stages. (Image courtesy of Dr. Tomofumi Maruta at FRIB Accelerator Physics Department) εx = εy and hence axisymmetry at the LS1 entrance. This subsection addresses three intertwined problems under these circumstances: 1) how to design the optics for the solenoid transport line LS1; 2) how to define the acceptance condition; and 3) how to match a non-axisymmetric beam into the design optics. In the following treatment, we first propose a design and matching recipe, and then discuss the underlying beam dynamics and why it will work. 4.3.3.1 Optics Design and Matching Recipe Firstly, the transverse optics of the LS1 solenoid transport line should be designed assuming the incoming beam is axisymmetric. This reduces the design to a 1D problem (optics and beam parameters identical along any transverse direction due to axisymmetry) which is easy to handle and is how it was done in the conceptual and engineering design of FRIB. The 1D design defines αin and βin, the incoming beam parameters required for matching, and εA, the maximum acceptable incoming emittance (i.e. acceptance). 62 Optimization:Stage 1Stage 2Stage 3Coupling Coeff. Next, we use the quasi-axisymmetric conditions defined in Sec. 4.2.5 to match a beam with εx (cid:44) εy into an axisymmetric optics by requiring: αx =αy = αin βx =βy = βin (cid:104)x y(cid:105) = 0 (cid:48)(cid:11) = 0 (cid:48)(cid:11) = −(cid:10)x(cid:48) y(cid:11) (cid:10)x(cid:48) (cid:10)x y y This is achieved by tuning the MEBT quadrupoles. The acceptance condition is fulfilled if εx ≤ εA and εy ≤ εA. To explain why the recipe works, we first note the transfer matrix of a linear solenoid can be decomposed as [30]: Msol = MRMF where MR is a rotation matrix and MF is focusing in both transverse directions: (4.53) (4.54) (4.55) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) MR = MF = cos ks(cid:96) 0 − sin ks(cid:96) 0 sin ks(cid:96) 0 cos ks(cid:96) 0 sin ks(cid:96) 0 cos ks(cid:96) 0 0 − sin ks(cid:96) 0 cos ks(cid:96) −ks sin ks(cid:96) 1 sin ks(cid:96) ks cos ks(cid:96) 0 0 0 0 cos ks(cid:96) −ks sin ks(cid:96) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) cos ks(cid:96) 0 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) 0 0 1 sin ks(cid:96) ks cos ks(cid:96) ks = Bz0/2[Bρ], (cid:96) is the effective length of the solenoid and Bz0 is the effective field. When a quasi-axisymmetric beam enters a linear solenoid, since the incoming beam has the same Twiss parameters in every transverse direction, and the focusing force is axisymmetric, the beam should retain quasi-axisymmetry upon exiting the solenoid. We arrive at this conclusion just 63 by invoking symmetry, but Eq. (4.53) tells us more. In addition to focusing, the linear solenoid also azimuthally rotates the beam by the Larmor angle ks(cid:96) [13]. This is relevant to the dynamics because Eq. (4.47) shows that εu, the emittance measured along ˆu where u = cos θ + sin θ, is dependent on θ. This means the transverse directions along which the beam size is maximum and minimum rotate by angle ks(cid:96) after passage through the solenoid. The beam evolution can be described using the following picture. The axisymmetric optics design and acceptance together defines a cylindrical envelope that varies in size along the solenoid transport line. For a quasi-asymmetric beam that is matched to the optics design with εx = εA > εy, the beam is spatially an ellipse as it enters the solenoid transport. The ellipse rotates in each solenoid and its two points of maximum extent always touches the cylindrical envelope. Thus a matched quasi-axisymmetric beam within acceptance guarantees that the beam is spatially bounded by the design cylindrical envelope throughout the solenoid transport. 4.3.4 Further Work In practice, the LS1 is not purely a solenoid transport because there exists quadrupole kicks from the quarter-wave SRF cavities. While the quadrupole kicks can be treated as small perturbations, it will be interesting to study how the solenoid optics can be designed such that the effect from such perturbations can be minimized. Additionally, construction and alignment errors can introduce further (hopefully small) symme- try breaking terms. Likewise, space charge nonlinearities of a nonuniform and non-axisymmetric beam can muddle the overall picture. Nevertheless, considerations presented in this chapter pro- vide a clear ideal basis for optimal tuning whose robustness against reasonable error sources can be evaluated with detailed simulations. 64 CHAPTER 5 TRANSVERSE PHASE-SPACE MEASUREMENTS WITH ALLISON-TYPE EMITTANCE SCANNERS The content of this chapter and Appendix A is an extended version of the work published in Ref. [31]. 5.1 Introduction Allison-type emittance scanners [32], or Allison scanners for short, are widely used to efficiently measure projected 2D phase space distributions of low-energy beams. An Allison scanner (see Fig. 5.1) consists of an entrance slit-plate (slit width s), an aligned exit slit-plate (slit width s) with an integrated Faraday cup, and a bipolar-biased electric dipole (voltage ±V0) placed between the two slits. The scanner assembly is translated mechanically (typically in discrete steps with a stepper motor) to change the slit position, and the dipole voltage V0 is swept at each step to select transmittable angles by varying the bending strength. The scanner samples one grid point in phase space for each coordinate and “E-dipole” voltage setting, with the beam density taken to be proportional to the current collected by the Faraday cup. The entire projected distribution is measured by recording currents collected over a range of position and voltage values that samples the full phase-space projection. Analysis of Allison scanner data requires a voltage-to-angle conversion and the angular resolu- tion of the device. Idealized relations were derived in Ref. [32] and summarized in Ref. [33]. They constitute the conventional analysis model which assumes the scanner has an idealized geometry with uniform hard-edge dipole fields, thin slit-plates and longitudinally symmetric placement of the E-dipole. Extending upon the work in [34], Sec. 5.2 augments the conventional analysis model to account for two commonly occurring geometric features that can lead to significant corrections due to: 1) the effective longitudinal asymmetry in E-dipole placement between the slits, and 2) the finite 65 slit-plate thickness. Their effects on particle transmission are studied using both simulations and analytic expressions which show excellent agreement. Asymmetry in E-dipole placement alters the slope of the linear relation between dipole voltage and particle angle. Finite slit thickness causes the weight of data points to vary as a decreasing function of their voltage values. The reduction becomes significant at large voltages and when the slits’ thickness-to-width ratio (cid:38) 1. When the interval between data points is smaller than the device resolution, the current measured at each data point is a weighted sum of the actual current densities at the data point and its neighbors. Sec. 5.3 shows that a detailed accounting of the phase space region contributing to each data point allows one to deconvolve the effective blur induced by overlapping data points and resolve the beam distribution more accurately. Reconstruction procedures are illustrated and verified using numerically generated phase-space distributions so that improvements are clearly characterized. Section 5.4 applies the improved model to experimental data from the Allison scanner used in the front-end of the Facility for Rare Isotope Beams (FRIB) at Michigan State University. Raw measurement data first undergo a noise removal scheme, then the processed data are analyzed with both the conventional model and the improved model. Comparisons reveal that the improved model produces more consistent results among a series of scans with changing focusing upstream. Beam moments are modified significantly relative to conventional treatments, with differences growing with the angular divergence of the beam distribution. The Python code tools used to perform the data analysis and noise removal are described and made available [35]. The modular programs are readily adaptable to other Allison scanners via a change of geometric input parameters. Section 5.5 concludes this study with an outline of directions for future improvements. 5.2 Additional Geometric Features This section first describes how asymmetric E-dipole placement and finite slit thickness can be modeled in Allison scanners. Effects induced are investigated using simulations that probe particle transmission given the device geometry and E-dipole voltage. The simulation results are accurately 66 modeled by analytic expressions for typical choices of scanner parameters. The implications for data analysis are demonstrated through simulated measurements of ideal particle distributions chosen to highlight those effects. 5.2.1 Geometric Model With reference to a realistic Allison scanner (see Fig. 5.1), the geometric information relevant to particle transmission and measurement data analysis can be captured using a model with six parameters. A schematic of the model for measuring the x-x(cid:48) phase space is illustrated in Case IV of Fig. 5.2 where • L is the E-dipole plate length; • l1 is the distance between the entrance-slit and the E-dipole plate; • l2 is the distance between the exit-slit and the E-dipole plate; • g is the gap between the E-dipole plates; • s is the slit width; and • d is the slit thickness. The model assumes that particle transmission has no y-dependence, which is a good approximation if no particle approaches the vertical extent of the device. This requires the vertical extent of the scanner to be much larger than 2(|ycen| + yrms) where ycen and yrms are the beam centroid position and rms width in the vertical direction. The geometric model hinges upon a correct definition of l1, l2 and d, so it is important to identify what performs the function of a slit. Figure 5.1 shows the design for the FRIB Allison scanner [lidia2016overview]. On each side of the E-dipole there is a thick plate with an opening (see zoom insert) that has large-angled relief cuts and a straight channel of minimum width. The relief cuts allow using thicker plates that enable cooling and have higher mechanical stability. The 67 relief cuts do not intercept any ions that would have otherwise passed through the slits, since the opening angle (30◦ in Fig. 5.1) far exceeds the transverse angles of beam ions. Hence the slit, i.e. the structure that limits particle transmission, is only taken to be the minimum-width channel, with width s = 60 µm and thickness d = 254 µm in Fig. 5.1. The relief cut is immaterial to the slit-to-E-dipole distance l1 in Fig. 5.1, hence l1 = 2.01 mm. In contrast, l2, the distance between the dipole plate and the exit slit, must take into account the longitudinal extent of the relief cut because it opens towards the E-dipole. Therefore, l2 = (2.06 + 3.18 − 0.254)mm = 4.986 mm. Note that the relief cuts as illustrated in a real system in Fig. 5.1 induce significant asymmetry in E-dipole placement with l2 = 4.986 mm > l1 = 2.01 mm. Conventional treatments model an Allison scanner with thin slits and symmetry in E-dipole placement, which corresponds to d = 0 and l1 = l2. However, real devices often deviate from these assumptions as discussed below. 5.2.1.1 Longitudinal Asymmetry in E-dipole Placement (l1 (cid:44) l2) Effective longitudinal asymmetry in E-dipole placement commonly arises in Allison scanners, probably unintentionally, because relief cuts on thick entrance- and exit-plates have been made in the same direction[36, 37]. This asymmetry was also present in the early implementation of the FRIB Allison scanner where relief cuts on both plates face the incoming beam as shown in Fig. 5.1. Since the relief cuts have no effect on particle transmission, it effectively displaces the end slit but not the entrance slit to make l2 > l1 in the FRIB design. We find that such asymmetries affect the transmitted angles for a given V0. 5.2.1.2 Finite Slit Thickness d (cid:38) Slit Width s The slit thickness d can be neglected when it is much smaller than the slit width s. This idealization starts to break down as Allison scanner designs employ decreasing slit widths to improve resolution. Recent examples include slit widths of s ≤ 100 µm in plans presented for GSI FAIR[38] and s = 38 µm implemented at TRIUMF[39]. Fig. 5.1 details the slit plate presently employed in the 68 Figure 5.1: Side-view of the FRIB Allison scanner. Table 5.1: Geometric parameters applied in examples corresponding to the four geometric models in Fig. 5.2 Case L [mm] 71.85 I 71.85 II III 71.85 71.85 IV Symmetric Asymmetric Symmetric Asymmetric l2 [mm] 3.498 4.986 3.498 4.986 s [µm] d [µm] E-dipole Placement l1 [mm] 3.498 2.01 3.498 2.01 g [mm] 7.91 7.91 7.91 7.91 60 60 60 60 0 0 254 254 Slits Thin Thin Thick Thick FRIB scanner where s = 60 µm and the slit thickness (not including irrelevant longitudinal extent of the 30◦ relief cut) is 254 µm. This approximately 4 : 1 aspect ratio effectively produces a narrow channel, which can scrape particles that would have otherwise passed through a slit-plate with no thickness. 5.2.2 Four Cases To investigate the effects of asymmetry in dipole placement and finite slit thickness, four Allison scanner models (see Fig. 5.2) are specified using the parameters listed in Table 5.1. The parameters in Case IV correspond to the geometric model of the FRIB Allison scanner in Fig. 5.1, which has asymmetric E-dipole placement (l1 (cid:44) l2) and thick slits (d (cid:44) 0). To mimic how the analysis would have been done if one or both geometric features are omitted, the FRIB design is modified into examples of Cases I to III (see Table 5.1) with thin slits (d = 0) and/or symmetric E-dipole 69 (cid:1)(cid:2)(cid:3)(cid:2)(cid:4)(cid:5)(cid:6)(cid:7)(cid:2)(cid:8)(cid:8)(cid:9)(cid:2)(cid:3)(cid:2)(cid:4)(cid:7)(cid:5)(cid:10)(cid:11)(cid:2)(cid:8)(cid:8)(cid:12)(cid:5)(cid:13)(cid:7)(cid:2)(cid:8)(cid:8)(cid:12)(cid:5)(cid:13)(cid:14)(cid:2)(cid:8)(cid:8)(cid:1)(cid:2)(cid:15)(cid:5)(cid:7)(cid:10)(cid:2)(cid:8)(cid:8)(cid:16)(cid:17)(cid:1)(cid:18)(cid:2)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:1)(cid:24)(cid:2)(cid:25)(cid:21)(cid:23)(cid:22)(cid:24)(cid:2)(cid:2)(cid:2)(cid:19)(cid:13)(cid:26)(cid:2)(cid:3)(cid:2)(cid:14)(cid:13)(cid:2)(cid:8)(cid:16)(cid:17)(cid:1)(cid:18)(cid:2)(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:1)(cid:24)(cid:2)(cid:25)(cid:21)(cid:23)(cid:22)(cid:24)(cid:2)(cid:2)(cid:2)(cid:27)(cid:19)(cid:13)(cid:28)(cid:23)(cid:29)(cid:23)(cid:30)(cid:23)(cid:31)(cid:2)(cid:32)(cid:33)(cid:34)(cid:14)(cid:13)(cid:8)(cid:12)(cid:11)(cid:35)(cid:2)(cid:8)(cid:15)(cid:5)(cid:7)(cid:10)(cid:2)(cid:8)(cid:8)(cid:15)(cid:13)(cid:20)(cid:36)(cid:20)(cid:20)(cid:8)(cid:2)(cid:20)(cid:37)(cid:2)(cid:26)(cid:21)(cid:17)(cid:22)(cid:38)(cid:39)(cid:21)(cid:17)(cid:22)(cid:151)(cid:151)(cid:151) placement (l1 = l2). Note that Case I corresponds to the conventional analysis model. All four cases have identical E-dipole length L, gap g and total inter-slit distance L + l1 + l2; thus all differences in analysis results arise from the geometric features. Characteristic trajectories in Fig. 5.2 illustrate the range of angles that can be transmitted at each voltage setting due to the finite slit width. We define: ref: angle at which an ion would enter and exit the slits at the same x-position 1, • x(cid:48) • x(cid:48) • ∆x(cid:48) = x(cid:48)max − x(cid:48) max / min: maximum/minimum transmittable angle, min: angular resolution, where x(cid:48) ≡ dx/dz and x = 0 corresponds to the device center-line. The x(cid:48)max-trajectory touches the lower edge of the entrance slit and the upper edge of the exit slit, and vice versa for the x(cid:48) min-trajectory. When the slits have finite thickness, one has to determine whether the extreme trajectory should touch the corner on the upstream or downstream side of each slit by selecting the combination that is not scraped. This complicates analytic modeling and results in a longitudinal shift of the extreme trajectories as visualized in Case III and IV of Fig. 5.2. 5.2.2.1 Particle Simulations Transmission properties of an Allison scanner design can be investigated using particle simulations. Two types of simulations are employed, namely “ballistic” and “realistic”, depending on the model of the dipole field. Realistic simulations employ a numerically calculated fringe field model for the dipole and include non-paraxial effects due to energy change. Ballistic simulations take the Allison scanner as an ideal dipole with uniform hard-edge fields over the axial extent L, with free drifts and collimating apertures on either side. Energy change from crossing E-dipole potential lines is neglected. Given V0, the particle trajectory consists of two straight lines outside the dipole and a parabola inside, and is scraped if it hits either of the slit plates. Space charge and scattering effects 1see Appendix A.3.2 for a subtle distinction in Case IV 70 Figure 5.2: Four geometric cases and their respective characteristic trajectories. Case I: l1 = l2 = l and d = 0; Case II: l1 (cid:44) l2 and d = 0; Case III: l1 = l2 = l and d (cid:44) 0; Case IV: l1 (cid:44) l2 and d (cid:44) 0. zi is the longitudinal position where the beam particles enter the device and where the x-x(cid:48) phase space is measured. are neglected in this study. A condition is derived in Appendix A.1 to show that non-paraxial effects are typically tiny. Details of the realistic model are given in Appendix A.2. Most simulation results shown are obtained using ballistic simulations because they are much faster than realistic simulations and can be summarized using analytic formulae. Results in Sec. 5.2.3 show that the differences between ballistic and realistic simulations are small for typical scanner geometries. The device center-line is chosen as the x-axis of the coordinate system for simplicity of de- scription. To analyze particle transmission at a given V0 value, the incoming angle which sends a 71 particle to the same position at both slits, i.e. x(cid:48) ref, is first obtained numerically via root-finding. Then, assuming that the phase space distribution is uniform over the acceptance region, an ensem- ble of incident particles are generated under the following conditions. Transversely, the particles have a uniform distribution about x(cid:48) = x(cid:48) ref and x = 0 with sufficient widths to fully populate the acceptance. Longitudinally, all particles have the same initial position zi and lie on the plane where the particles enter the device (see Fig. 5.2). For slits with finite thickness, the particles are aligned with the outer plane of the entrance slit. Every particle is advanced with ballistic simulations and all that can pass through the exit slit are recorded. Using 40Ar9+ ions with kinetic energy E = 12 keV/u as test particles, the phase space region that can be transmitted at V0 = 500V and V0 = 1000V for the four geometric cases in Table 5.1 are plotted in Fig. 5.3 for ballistic simulations. The red rectangle represents the position resolution ∆xideal and angle resolution ∆x(cid:48) ideal predicted by the conventional treatment ([32, 33], also reproduced under Case I in Table 5.2), where (5.1) ∆xideal = s, ∆x(cid:48) ideal = 2s L + 2l . A comparison of Cases I and III (l1 = l2) against Cases II and IV (l1 (cid:44) l2) reveals that x(cid:48) ref changes as a result of longitudinal asymmetry, even though the parameters are chosen such that all cases have the same inter-slit distance. The effects of finite slit thickness are manifested in Cases III and IV where the phase space area transmitted shrinks in both dimensions because d ≈ 4s (cid:44) 0. Furthermore, the shrinkage is larger in Fig. 5.3a than in Fig. 5.3b - this indicates that this effect increases with V0, which corresponds to a larger selected reference angle x(cid:48) ref. An illustration of the details of how the shrunk phase space area attains its shape is given by Fig. 5.4. 5.2.2.2 Particle Transmission Quantified To quantify the transmission properties at a specified value of V0, assuming the phase space distribution is uniform about x = 0 and x(cid:48) = x(cid:48) ref, the density of transmitted particles in phase space 72 (a) V0 = 500V (b) V0 = 1000V Figure 5.3: Blue dots denote particles that are transmitted in the simulations described in Sec. 5.2.2.1. They map x-x(cid:48) phase space regions at the device entrance z = zi that contribute to the data point at x = 0 mm. Results are shown for all four geometric cases in Fig. 5.2 for (a) V0 = 500V and (b) V0 = 1000V respectively. Geometric parameters for each case are listed in Table 5.1. Black lines define the phase space grid with increments determined by the scanning step sizes. Dashed blue lines indicate the nominal position x = 0 and angle x(cid:48) = x(cid:48) ref corresponding to the data point. The red rectangles denote the device resolution predicted by the conventional treatment (i.e. Case I). can be defined as f(x, x(cid:48)) =  1 0 if (x, x(cid:48)) ∈ A if (x, x(cid:48)) (cid:60) A , (5.2) where A = {(x, x(cid:48))| particle at (x, x(cid:48)) is transmittable}. Graphically, in Fig. 5.3, f(x, x(cid:48)) = 1 inside the blue area flagging the simulated phase space particles that are transmittable and f(x, x(cid:48)) = 0 outside. Two quantities that distill useful information concerning the transmitted phase space area are the angular transmission factor: T(x(cid:48)) = 1 f( ˜x, x(cid:48))d ˜x, (5.3) Þ ∆xideal/2 ∆xideal −∆xideal/2 73 41.041.542.042.543.043.544.0Angle x′ [mrad]Case Id=0l1=l2Case IIId≠0l1=l2−0.050.000.05Position x [mm]42.543.043.544.044.545.045.5Angle x′ [mrad]Case IId=0l1≠l2−0.050.000.05Position x [mm]Case IVd≠0l1≠l283.584.084.585.085.586.086.5Angle x′ [mrad]Case Id=0l1=l2Case IIId≠0l1=l2 0.050.000.05Position x [mm]87.087.588.088.589.089.590.0Angle x′ [mrad]Case IId=0l1≠l2 0.050.000.05Position x [mm]Case IVd≠0l1≠l2 Figure 5.4: Illustrates how the shrunk phase space area attains its shape in Case III of Fig. 5.3b. and total transmission factor W = 1 ∆x(cid:48) ideal Þ x(cid:48) ref+∆x(cid:48) x(cid:48) ref−∆x(cid:48) ideal/2 ideal/2 T( ˜x(cid:48))d ˜x(cid:48) . (5.4) Here, ∆x(cid:48) and ∆x(cid:48) ideal are the ideal position and angular resolution of the device described by Eq. (5.1), and ˜x and ˜x(cid:48) are dummy variables used to denote integration over position and angle respectively. The angular transmission factor T(x(cid:48)) represents the fraction of particles entering the slit with angle x(cid:48) that is transmitted, assuming a uniform distribution in initial x over the slit width. T(x(cid:48)) can be interpreted as a normalized projection of the blue areas in Fig. 5.3 onto the coordinate x(cid:48) by integrating over x. T(x(cid:48)) corresponding to the four cases in Fig. 5.3b is plotted in Fig. 5.5. For Cases I and II with thin slits, T(x(cid:48)) = 1 when x(cid:48) = x(cid:48) ref and decreases linearly until T(x(cid:48)) = 0 at ideal. Thick slits in Cases III and IV diminish both the range of transmittable x(cid:48) and x(cid:48) = x(cid:48) the fraction of particles transmitted at a given x(cid:48). ref ± ∆x(cid:48) 74 xzllddxzxz Figure 5.5: Angular transmission factor T(x(cid:48)) for Cases I to IV in Fig. 5.3b. ideal defined by the device resolution. The total transmission factor W represents the area ratio between the transmittable phase space region and ∆xideal × ∆x(cid:48) It can be thought of as the ratio between the blue area and area of the red rectangle in any individual plot in Fig. 5.3. Alternatively, it is the ratio between the area under T(x(cid:48)) and the area enclosed by the two vertical dotted lines in Fig. 5.5. W always satisfies W ≤ 0.5 with W = 0.5 corresponding to the ideal geometry (Case I). This result means that within the rectangular phase space region defined by the resolution obtained conventionally, no more than half of the area is actually transmittable. The fact that W decreases with V0 in Cases III and IV leads to significant corrections in the data analysis. 5.2.2.3 Analytic Results Exploiting the simple geometry of particle trajectories in free drifts and ideal dipole fields, all quantities relevant to particle transmission defined in this section can be derived analytically following the steps outlined in Appendix A.3. Table 5.2 summarizes the formulae derived for all four cases. The formulae presented have been numerically verified. The linear coefficient in the voltage-to-angle relation giving x(cid:48) ref as a function of V0 is shown to change with asymmetry in dipole placement. Expressions for T(x(cid:48)) in Cases I to III agree with 75 −0.8−0.6−0.4−0.20.00.20.40.60.8Relative Angle x′−x′ref [mrad]0.00.20.40.60.81.0Angular Transmission T(x′−x′ref) Ɗx′ideal2−Ɗx′ideal2Δase IΔase IIΔase IIIΔase IV (cid:104)1 − (x(cid:48) Fig. 5.5 where T(x) is maximum at x(cid:48) = x(cid:48) refd/s)(cid:105) in Case III arise from effects resulting from the finite slit thickness. The factor is ref and declines linearly to 0 at x(cid:48) = x(cid:48) max/min. Factors of ≤ 1 and decreases linearly with increasing V0 since x(cid:48) the blue areas in Fig. 5.3 when d (cid:44) 0. ref ∝ V0. This effect causes the shrinkage of T(x(cid:48)) and W in Case IV are more complicated due to the coupling between asymmetric dipole placement and finite slit thickness. Details of the derivation can be found in Appendix A.3.2. However, in most relevant cases, T(x(cid:48)) and W of Case III and Case IV are very close except at large angles. This can be observed in Fig. 5.5 where the two plots almost overlap even at V0 = 1000 V (Fig. A.4 in Appendix A.3.2 highlights the differences when they are discernible). Therefore, the discussion concerning Case III typically applies to Case IV with only small modifications. The agreement between analytic formulae and simulation results, both ballistic and realistic, is demonstrated in the next section where we discuss how the additional geometric features lead to two corrections. 76 Table 5.2: Analytic formulae corresponding to the geometric models in Fig. 5.2 for V0 ≥ 0. Results for V0 < 0 obey the same formulae with V0 → |V0| and the x-axis reversed. I IV III II symmetric (l1 = l2) thin (d = 0) asymmetric (l1 (cid:44) l2) (cid:16)1 + l2−l1 thin (d = 0) qV0L gE L+l1+l2 (cid:17) 1 2 Case Dipole placement Slit thickness x(cid:48) ref x(cid:48)max − x(cid:48) ref x(cid:48) ref − x(cid:48) min ∆x(cid:48) T(x(cid:48)) for x(cid:48) ≥ x(cid:48) ref T(x(cid:48)) for x(cid:48) < x(cid:48) ref W 1 2 qV0L gE s L+2l s L+2l 2s L+2l x(cid:48)max−x(cid:48) x(cid:48)max−x(cid:48) ref x(cid:48)−x(cid:48) min x(cid:48) ref−x(cid:48) min 1 2 symmetric (l1 = l2) thick (d (cid:44) 0) (cid:19) (cid:18) 1 qV0L gE 2 s−x(cid:48) refd L+2l+d s−x(cid:48) refd L+2l+d (cid:18) 1 − x(cid:48) (cid:18) 1 − x(cid:48) (cid:19)2 1 − x(cid:48) 1 − x(cid:48) refd s refd s refd s refd s (cid:19) (cid:19) 2s L+2l+d x(cid:48)max−x(cid:48) x(cid:48)max−x(cid:48) ref x(cid:48)−x(cid:48) min x(cid:48) ref−x(cid:48) min 1 2 (cid:18) 1 2 (cid:17) asymmetric (l1 (cid:44) l2) thick (d (cid:44) 0) qV0L gE L+l1+l2 (cid:16)1 + l2−l1 (cid:16) (cid:18) s−x(cid:48) refd L+l1+l2+d s − L+2l1 L+2l2 x(cid:48) refd x(cid:48) 1 − L+l1+l2 refd L+2l2 s (cid:17) 1 2s (cid:19) L+l1+l2+d L+l1+l2+d see Appendix A.3.2 see Appendix A.3.2 see Appendix A.3.2 s L+l1+l2 L+l1+l2 s 2s L+l1+l2 x(cid:48)max−x(cid:48) x(cid:48)max−x(cid:48) ref x(cid:48)−x(cid:48) min x(cid:48) ref−x(cid:48) min 1 2 77 5.2.3 Two Corrections in Data Analysis Analytic and simulation results above reveal two significant corrections on the data analysis of Allison scanners relative to the conventional treatment that assumes symmetric dipole placement and thin slits [32, 33]. One correction arises from asymmetry in dipole placement and alters the proportionality constant in the linear relation between dipole voltage (V0) and selected angle (x(cid:48) ref). The other correction adjusts the weights of data points as a function of x(cid:48) ref to compensate for scraping due to the finite slit thickness. Both corrections can be applied using the quantities analytically expressed in Table 5.2. We illustrate these effects for the FRIB Allison scanner geometry shown in Fig. 5.1 for an 40Ar9+ ion with kinetic energy E = 12 keV/u. 5.2.3.1 Voltage-to-Angle Relation As shown in Cases II and IV of Table 5.2, when the dipole placement is asymmetric, the reference angle x(cid:48) ref has an additional factor [1 + (l2 − l1)/(L + l1 + l2)] relative to the symmetrical model results (Case I and III). For the FRIB Allison scanner, l2−l1 = 2.976 mm, so the relief cuts generate substantial effective asymmetry. In Fig. 5.6, x(cid:48) ref versus V0 is plotted for the symmetrical treatment (Table 5.2, Case I or III), improved analytic results (Table 5.2, Case II or IV), and both ballistic and realistic numerical simulations. Both types of simulations agree almost perfectly with analytic results to show a slope of 88.4 mrad/kV, which deviates from the symmetrical result 85.2 mrad/kV by roughly 4% as a consequence of the effective asymmetry. The close correspondence between realistic simulations and analytic results shows that fringe-fields and non-paraxial effects impose minimal changes on particle trajectories for typical device geometries and beam parameters, thus validating the analytic calculations as an accurate and efficient way to obtain x(cid:48) ref. 5.2.3.2 Voltage-Dependent Weight Compensation As mentioned in Sec. 5.2.2, thick slits cause W to decrease as a function of x(cid:48) ref. To understand implications of this effect, it is important to delineate what the blue area, red rectangle and black 78 Figure 5.6: Voltage-to-angle relation for the selected beam. The red curve (analytic) overlays the blue curve (ballistic simulation). rectangle represent in Fig. 5.3. Phase-space scans are typically performed by measuring the current at regular discrete steps in slit position (x) and voltage (V0). After voltage-to-angle conversion, the data are represented by a grid in phase-space where each grid represents a rectangle of size (position step) × (angle step). In general, grid dimensions will not coincide with the spatial and angular resolution of the device. In Fig. 5.3, the black rectangle denotes the phase space grid while the conventional device resolution defines the red rectangle. In this study, scanning steps exceed the respective device resolution in both position and angle, as is typical of scans taken in FRIB (see Sec. 5.4). Assuming a uniform distribution of particles within each grid, the actual beam current falling within each phase space grid cell should be calculated by Igrid = Ifc × Black Rectangle Area Red Rectangle Area × 1 W , (5.5) where Ifc is the Faraday cup current collected at the data point, and W is the ratio between blue area and red rectangular area defined in Sec. 5.2.2. With thin slits, Eq. (5.5) merely rescales all data points uniformly because W = 1/2 is a constant. Thus one does not need to distinguish between Igrid and Ifc when one calculates moments of the 79 02004006008001000Dipole Voltage V0 [V]0102030405060708090Angle x′ref [mrad]Conventional AnalysisSimulation (Ballistic)Simulation (Realistic)Analytic beam distribution in phase space with measured data, unless one also wants to compute the total beam current accurately to verify normalizations. With thick slits, W is not constant and decreases quadratically in (x(cid:48) ref × d/s) (see Table 5.2). Therefore, Ifc has to be rescaled differently depending on the angle x(cid:48) ref. The larger the slit thickness-to-width ratio d/s, the more important this effect becomes. This effect must be taken into account for accurate moment calculations. In Fig. 5.7, the angular resolution ∆x(cid:48) and 1/W are plotted as a function of x(cid:48) ref for the FRIB Allison scanner. Curves for conventional treatment (Table 5.2, Case I), improved analytic results (Table 5.2, Case II or IV), and numerical simulations are shown. Note again the minimal difference between ballistic and realistic simulation results. ref = 70 mrad, W−1 ≈ 4 as opposed to W−1 = 2 when x(cid:48) The correction introduced by the effect of finite slit thickness can be substantial. For example, in Fig. 5.7b, when x(cid:48) ref = 0 mrad. Therefore, if we apply the conventional model where all data points have equal weights, 0-mrad data points would wrongly weigh twice as much as ±70-mrad data points. Failure to correct the weights will distort beam measurements and moment evaluations with increasing amplitude as beam angular extent increases. In the laboratory, increasing the angular extent of the beam phase space distribution at the Allison scanner by adjusting the applied focusing lattice is often needed to obtain sufficient angular information. An order-of-magnitude inspection based on ideal device resolution in Eq. (5.1) reveals the reason. (L + 2l) typically ∼ 0.1 m, therefore ∆xideal ∆x(cid:48) ideal (cid:39) L + 2l 2 (cid:39) 1 20 mm mrad . (5.6) Since measurements usually take scanning steps exceeding the device resolution, Eq. (5.6) says that if a similar level of detail is desired in both dimensions, the distribution’s angular extent in mrad should be around 20 times its position counterpart in mm. This often forces one to employ relatively narrow scanner slit width s at the expense of reduced Faraday cup currents that lead to greater noise issues, particularly for heavy-ion beams whose currents are typically lower than that of proton beams. 80 (a) (b) Figure 5.7: Plots of a) angular resolution ∆x(cid:48), and b) angular correction factor as a function of x(cid:48) ref. The red curve (analytic) overlays the blue curve (ballistic simulation). 81 020406080100Angle x′ref [mrad]0.80.91.01.11.21.31.41.51.6Angular Re olution Ɗx′ [mrad]Conventional Analy i Simulation ΔBalli tic)Simulation ΔReali tic)Analytic020406080100Angle x′ ef [m ad]0123456T ansmission Facto Inve sed W−1Conventional AnalysisSimulation (Ballistic)Simulation (Realistic)Analytic 5.2.4 Analysis of Synthetic Data The effects of the two corrections on data analysis can be tested using synthetic measurement data as follows. After generating an idealized distribution of particles in phase space with given Twiss parameters and emittance, the transmission of each particle at each (x,V0) pair of the scan can be tested using ballistic simulations. The particle number transmitted at each (x,V0) setting is recorded to emulate the measurements that would have been taken by an Allison scanner. It is unnecessary to test every (x,V0) pair for a given particle as many pairs are ruled out based on the device resolution. The artificial data set thus synthesized is processed using both the conventional and improved analysis to analyze errors in beam moments and phase space distortions in comparison with the idealized distribution. Beam moments can be defined in terms of a continuous particle distribution function represent- ing the beam as: (cid:10)g(u, u(cid:48))(cid:11) ≡ Þ ∞ −∞ Þ ∞ Þ ∞ Þ ∞ −∞ f(u, u(cid:48))g(u, u(cid:48))dudu(cid:48) −∞ −∞ f(u, u(cid:48))dudu(cid:48) (5.7) (5.8) (5.9) where f(u, u(cid:48)) is the distribution function in u-u(cid:48) phase space, and g(u, u(cid:48)) = ua1u(cid:48)a2 with a1 + a2 being the order of the moment. Here, u represents either x or y. Beam moments can be equivalently calculated from the particle perspective as: for a distribution with N discrete particles, and within a gridded phase space as: g(ui, u(cid:48) i) 1 N N (cid:10)g(u, u(cid:48))(cid:11) = i=1n m (cid:10)g(u, u(cid:48))(cid:11) = m i=1n i=1 j=1 Ji j g(uj, u(cid:48) i) j=1 Ji j where Ji j is the current density in u-u(cid:48) phase space in the i, j cell. Equation (5.9) is equivalent to discretizing Eq. (5.7). A data set for the Allison scanner geometry in Case IV of Table 5.1 is generated using a 12 keV/u 40Ar9+ beam with KV distribution [20, 40]. The KV distribution has uniform projected 82 Figure 5.8: Uniform phase space projection of the KV distribution, compared against results of synthesized measurement data analyzed using conventional and improved analysis methods. The left plot shows a sampling (1 out of 1000 particles) used to represent the distribution. The middle and right plots are fine mesh density plots with the same color scale. density with an x-x(cid:48) phase space ellipse. This simple uniform projection makes it easy to visualize distribution distortions. The KV distribution and results from conventional and improved analysis methods are plotted in Fig. 5.8. Although the KV distribution has a uniform x-x(cid:48) projection within an elliptical boundary by construction, in the synthesized measurement, the beam density artificially diminishes strongly at larger angles under the conventional analysis, whereas the correct uniformity is accurately preserved under the improved analysis. Beam moments and the associated Twiss parameters αx, βx with rms emittance εrms are listed in Table 5.3. The results of the improved analysis match the simulated distribution very closely, whereas there is > 10% deviation in the Twiss parameters if the conventional analysis is used. The conventional analysis also underestimates the total beam current by 25%. These deviations in the conventional analysis will amplify if the angular extent of the KV distribution increases. 5.3 Corrections for Degenerate Phase Space Measurements In the previous section, all analysis was made assuming that the measurement step size of the scan always exceeds the device resolution. This assumption is illustrated in Fig. 5.3, where the red rectangle denoting the device resolution is completely encompassed by the black grid defined by the interval between data points. In this section, we discuss how to correct measurement results 83 −15−10−5051015Position x [mm]−80−60−40−20020406080Angle x′ [mrad]KV Distrib tion−10−50510Position x [mm]Conventional Analysis−10−50510Position x [mm]Improved Analysis0.00.20.40.60.81.00.00.20.40.60.81.0 Table 5.3: Beam parameters of the three distributions in Fig. 5.8. Parameters (cid:104)xx(cid:105) [mm2] (cid:104)xx(cid:48)(cid:105) [mm mrad] (cid:104)x(cid:48)x(cid:48)(cid:105) [mrad2] αx βx [m] εrms [mm mrad] Itotal [mA] Distribution KV 25.0 192.0 1600 -3.43 0.446 56.0 2.00 Conventional Analysis 20.9 156.5 1304 -3.00 0.400 52.1 1.49 Improved Analysis 25.0 192.3 1603 -3.43 0.447 56.0 2.00 Figure 5.9: Phase space region contributing to the data point at x = 0 mm and V0 = 200 V (presentation format as detailed in Fig.5.3), where the Allison scanner geometry corresponds to Case IV of Table 5.1 except for increased slit width to s = 300 µm. Scanning step sizes in x and V0 are 0.5 mm and 20 V respectively. when this assumption is no longer valid. Intuitively, when the device resolution exceeds the grid size, the current measured has contri- butions from multiple grid cells. Assuming this occurs only in the angle coordinate but not in the position coordinate, Fig. 5.9 shows an example where the phase space region extends over five grid cells vertically. The current measured at a data point is thus a weighted sum of the current densities from contributing grid cells. Assuming uniform current density in each grid cell, the weighting factor equals the overlapping area between the blue quadrilateral and the grid cell. Therefore, for each column in the phase space plot, if the array of angular values are numbered from 1 to n, a 84 Angle x' [mrad] system of linear equations can be written down (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) A1,1 ... An,1 · · · A1,n ... . . . · · · An,n (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) j1 ... jn (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) , I1 ... In (5.10) where In is the current measured at the n-th angle value, An,m is the area in m-th grid occupied by the phase space region contributing to In, and jn is the current density at the n-th grid. Fig. 5.9 shows several elements of the matrix A, where In corresponds to the current measured at V0 = 200 V in this case. Utilizing detailed information from Sec. 5.2 on the phase space region contributing to a data point, the matrix A can be obtained via transmission simulations or numerical integration of T(x(cid:48) − x(cid:48) ref) (see Table 5.2). The current densities j can thus be solved for as unknowns to deconvolve the information intermixed in current measurements. Corrections for degenerate measurements are tested using a synthesized data set generated with the methods described in Sec. 5.2.4 using a waterbag distribution. The waterbag distribution [20, 40] is applied corresponding to uniform transverse energy in phase space up to a sharp cutoff value, with self-field energy neglected. The Allison scanner geometry is the same as Case IV of Table 5.1 except s = 300 µm (i.e. wider slits), while the scanning step size in V0 is 10V. The uncorrected and corrected phase space distributions are plotted in Fig. 5.10. Beam parameters are listed in Table 5.4. Note that the uncorrected phase space distribution appears smoother than the corrected distribu- tion because the uncorrected current measured at each grid point takes contributions from multiple neighboring grid cells, thereby causing a blurring effect. In other words, while the grid size is determined by the scanning step size and can be made arbitrarily small, details finer than the device resolution are misleading unless corrections are applied to deconvolve overlapping measurements. Figure 5.10 shows that the uncorrected distribution also has ghost tails that leak out of the 100% beam ellipse on both sides of the angular extent of the distribution. The effect is illustrated clearly when all non-zero points are colored equally (second row in Fig. 5.10). This arises because currents 85 Figure 5.10: Waterbag distribution (left column), compared against synthesized measurement data analyzed without (middle column) and with (right column) corrections for degenerate measurements. The upper row corresponds to intensity, whereas the lower row shows all non-zero intensity points in the same color to highlight extent errors. Black and white ellipses denote the ideal boundary of the waterbag distribution analyzed. In the left column sampled particles (1 out of 1000) are plotted. Measurement data are synthesized with s = 300 µm, and scanning step sizes of 0.5 mm and 10 V. measured at points outside the ellipse also take contribution from interior grid points, so the extent of the leakage is determined by the device resolution. Fig. 5.11 shows that such ghost tails can cause significant overestimation of particles with large angles. Therefore, a correct characterization of degenerate measurements and device resolution is potentially important for studies of beam halo where particles outside the core can extend to large angles. The treatment above assumes the device resolution only exceeds the step size in x(cid:48) but not in x. This is quite common because the angular resolution limit is much more likely to be reached by the angular step size, rather than the same occurring in the position dimension, due to the ease of electronically sweeping and recording voltage steps relative to mechanical translation of the 86 Figure 5.11: Fraction of particles that fall outside the ellipse (orientation defined by Twiss parameters) with area πε in each of the three distributions in Fig. 5.10. The vertical dashed line corresponds to the sharp phase space edge of the waterbag distribution at ε/εrms = 6. The blue curve approaches ε/εrms = 6 if the lower range of the ordinate is extended. Table 5.4: Beam parameters of the two distribution in Fig. 5.10 compared against those of the waterbag distribution Conventional Phase Space Correction Analysis 9.03 11.91 27.48 -1.16 0.88 10.31 2.00 9.03 11.91 25.01 -1.30 0.99 9.17 2.00 Parameters (cid:104)xx(cid:105) [mm2] (cid:104)xx(cid:48)(cid:105) [mm mrad] (cid:104)x(cid:48)x(cid:48)(cid:105) [mrad2] αx βx [m] εrms [mm mrad] Itotal [mA] Waterbag Distribution 9.00 12.00 25.00 -1.33 1.00 9.00 2.00 87 0246810Relative Emi ance ε/εrms10−310−210−1100Frac ion of Par icles Ou side Area πεεmax in idealwa erbag dis ribu ionWa erbag Dis ribu ionw/o Degeneracy Correc ionsw/ Degeneracy Correc ions scanner (see also last paragraph of Sec. 5.2.3). If the step size in x is smaller than s, the method above can be extended to 2D whereby all points in the plot have to be solved at once. This can result in a very large matrix inversion problem unless contributions are limited to some level of nearest-neighbor zones. 5.4 Example: FRIB Allison Scanner This section applies the improved analysis model on measurements at the front end of the Facility for Rare Isotope Beams (FRIB). Sec. 5.4.1 describes the FRIB front end. Section 5.4.2 highlights the discrepancies between conventional and improved analysis results, and validates the improved treatment by demonstrating that it delivers more consistent results among a series of measurements performed in a quadrupole scan. The underlying Python code tools we employ for data analysis of Allison scanner measurements have been made publicly available in a git [41] software repository [35]. Tools include an adaptable code to implement the ideal and improved analysis methods described in Sec. 5.2, and algorithms to process measurement data for noise thresholding. These tools are self-contained, documented, and readily adaptable to a whole variety of applications. They are made available as a community resource. The version used in this study corresponds to November, 2018. Future algorithm advances will also be posted and documented in the repository. The algorithms have also been incorporated into a graphical user interface (GUI) developed by Tong Zhang that runs as an application in the FRIB control system. Information on this GUI is given in Appendix A.5. 5.4.1 FRIB Front End The front end of the Facility for Rare Isotope Beams (FRIB) [2] began commissioning in June 2017. A schematic of the upstream end of the low energy beam transport (LEBT) line is shown in Fig. 5.12. Beam ions are produced in an Artemis-type ECR source [5] with a 15 kV extraction voltage. The beam then traverses a short transport section with solenoid focusing and an electrostatic gap biased to accelerate the target ion species to 12 keV/u. Species are separated horizontally (in 88 Figure 5.12: Upstream end of the low energy beam transport (LEBT) line of the FRIB front end. x) via a large dispersion generated in a 90◦ magnetic dipole. This is followed by an electrostatic quadrupole triplet and a four-jaw collimator to scrape unwanted ion species in a high dispersion region. Two Allison scanners, one for each transverse direction, are located downstream of the four-jaw collimator. The collimator jaws are adjusted such that only the target species and possibly traces of contaminant ions whose mass-to-charge ratios M/Q differ from the target ion species by < 1% remain. The initial implementation of the Allison scanner had geometries corresponding to Case IV in Table 5.1. Motivated by this work, the direction of the relief cuts on the entrance slit plate was flipped in February 2018, thus greatly reducing the longitudinal asymmetry in E-dipole placement. Table 5.5 compares the geometries before and after the modification. All measurements in this section were obtained under the new device geometry, hence the corrections have negligible 89 ARTEMIS ECRQuadrupole Triplet4-Jaw CollimatorAllison ScannerFaraday CupSolenoids90o DipoleElectrostatic Gap Table 5.5: Allison scanner geometries at FRIB front end, before and after modification in Feb 2018. Design Before Symmetrization After Symmetrization L [mm] 71.85 71.85 l1 [mm] 2.01 4.938 l2 [mm] 4.986 4.986 g [mm] 7.91 7.91 s [µm] 60 60 d [µm] 254 254 contribution from asymmetry and mainly arise from thick slit effects. 5.4.2 Quadrupole Scan Measurements A set of six measurements were made on a 50 µA (measured separately by Faraday cup) Ar9+ ion beam using the y-plane Allison scanner, where VQ3, the voltage at the last quadrupole in the triplet upstream, varied from 0V to 5000V at 1000V intervals. Each scan was performed with y-spatial steps of 1 mm and voltage steps of 20 V. Raw data were processed using the noise removal procedures outlined in Appendix A.4. The resulting data set was analyzed using both the conventional and improved methods as outlined in Sec. 5.2. Three y-y(cid:48) phase-space projections (with corrections) at different VQ3 are shown in Fig. 5.13. To illustrate the differences between the conventional and improved analysis methods, we compute moments of the beam measured at VQ3 = 1000 V both with and without geometric corrections. The calculations with geometric corrections use the post-Feb 2018 geometry in Table 5.5, while the geometry’s symmetric (l = [l1 + l2]/2 = 4.962 mm on both ends) and thin-slit (d = 0 µm) counterpart gives results with both geometric features ignored. (cid:10)y Results are listed in Table 5.6. All second-order moments have centroids subtracted, i.e., 2(cid:11) ≡(cid:10)(y − (cid:104)y(cid:105))2(cid:11). We observe differences up to ∼10% in beam moments, thus confirming the importance of corrections for accurate measurements. While the example above shows the magnitude of the corrections, it does not provide evidence on which result is more valid. If the corrections introduced in the analysis constitute an “improved” treatment over the conventional one, a set of quadrupole scan measurements analyzed with the former should exhibit better consistency than those analyzed with the conventional methods. This check is performed below. 90 Figure 5.13: y-y(cid:48) phase-space projections measured at the FRIB front end during a quadrupole scan. VQ3 denotes the voltage at the last quadrupole upstream of the Allison scanner. 91 −10−50510152025Position y [mm]−40−2002040Angle y′ [mrad]VQ3=1000 Vεrms=20.1 mm mrad−100102030Position y [mm]−20−100102030Angle y′ [mrad]VQ3=3000 Vεrms=20.7 mm mrad−100102030Position y [mm]−20−100102030Angle y′ [mrad]VQ3=5000 Vεrms=20.0 mm mrad Table 5.6: Beam parameters of corresponding to the phase space distribution at VQ3 = 1000 V in Fig. 5.13, computed with and without correction terms. Conventional Phase Space Percentage Correction Difference Analysis Parameters (cid:104)y(cid:105) [mm] (cid:104)y(cid:48)(cid:105) [mrad] (cid:104)yy(cid:105) [mm2] (cid:104)yy(cid:48)(cid:105) [mm mrad] (cid:104)y(cid:48)y(cid:48)(cid:105) [mrad2] εrms [mm mrad] Itotal [µA] 7.02 4.36 8.84 17.8 76.9 19.0 36.8 7.09 4.64 9.76 19.6 80.8 20.1 39.6 1.0% 6.0% 9.4% 9.2% 4.8% 5.5% 7.1% Conclusive evidence of improvement arises from first order moments of the beam. Measurement results of a quadrupole scan and the linear transfer map at each setting can be used to form a system of linear equations: (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) M(VQ3 = V1) ... M(VQ3 = Vn) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) I2 ... I2 (cid:169)(cid:173)(cid:173)(cid:171)xi x0 (cid:170)(cid:174)(cid:174)(cid:172) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) xf(VQ3 = V1) ... xf(VQ3 = Vn) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) . (5.11) x(cid:48)(cid:17)T (cid:16) Here, VQ3 is the voltage of the last quad, I2 is the 2× 2 identity matrix, M(VQ3) is the 2× 2 transfer matrix consisting of the quadrupole at VQ3 and the quadrupole-to-scanner drift, and x = is the state vector with xi being the initial centroid position and angle entering the quad, x0 being the position offset and tilting angle of the Allison scanner with respect to the centerline of the quad, and xf(VQ3) being the measured centroid position and angle corresponding to the measurement at VQ3. x Equation (5.11) solves for the initial centroid position and angle of the beam, as well as the position offset and tilting angle of the Allison scanner. xf(VQ3) is computed differently depending on whether the conventional or improved analysis is applied. Results listed in Table 5.7 show that solutions from the improved analysis are much more realistic. Tilting angles of 33.6 mrad would not go undetected in mechanical alignment, nor would a beam with centroid angles (cid:104)x(cid:48)(cid:105) = 37.7 mrad transport through the downstream beam line without significant loss. Therefore, the experimental data support the validity of the improved analysis over the conventional treatment. 92 Table 5.7: Solution of Eq. (5.11) for the quadrupole scan in Sec. 5.4.2, with measurement results x f calculated using conventional and improved analysis method respectively. Analysis Conventional Improved Beam x(cid:48) i xi [mm] -5.9 -0.7 [mrad] 37.7 5.7 Scanner x(cid:48) 0 [mrad] -33.6 -1.2 x0 [mm] -18.0 3.0 The large differences between the two sets of solutions in Table 5.7 can be explained as follows. xf(VQ3) from the conventional analysis can be viewed as a perturbation to xf(VQ3) from the improved analysis. The sensitivity of the solutions xi and x0 to the perturbation can be quantified by the condition number [42] of the coefficient matrix in Eq. (5.11). Quadrupole scan parameters corresponding to larger condition numbers will be more sensitive to the perturbation. For the case presented in Table 5.7, the condition number is large. Choosing scan parameters corresponding to a well-conditioned coefficient matrix with a smaller condition number will be conducive to error minimization. Such issues have been explored in Sec. 6.3 5.4.3 Total Current Normalization Anomaly Despite the successful experimental verification described in Sec. 5.4.2, the quadrupole scan measurements also present a conundrum. Using Eq. (5.5) and detailed knowledge of the phase space area contributing to each data point summarized in Table 5.1, it is possible to calculate the current within each grid cell in an Allison scanner measurement. The summation of all grid cell current contributions should be consistent with the total beam current that can be accurately measured by the full-beam Faraday cup (see Fig. 5.12) located at approximately the same beamline location as the Allison scanner. Results for the quadrupole scan in Sec. 5.4.2 are shown in Fig. 5.14. The beam current is 50 µA according to direct Faraday cup measurements. However, the projected total currents from Allison scanner measurements are consistently ≈ 40 µA, which are 20% smaller than expected. Note that the improved analysis does make a small small correction, and the agreement is even worse if the conventional analysis is applied. But the discrepancies are still large. 93 Figure 5.14: Projected beam currents from Allison scanner measurements in a quadrupole scan, compared against the beam current measured at the Faraday Cup slightly upstream (see Fig. 5.12). Causes of the discrepancy are unknown and remain a subject of investigation. The bias voltage on the Allison scanner Faraday Cup, which has a great impact on the current reading, has been eliminated as a possible cause by the measurement results in Fig. 5.15. When ions hit the Faraday cup, secondary electrons may escape and amplify the measured current. Therefore, a bias voltage has to be applied to a Faraday cup to suppress such “ghost” currents. Fig. 5.15 shows that the measured current on the Faraday cup converges at bias voltages below -100 V, and FRIB usually employs a bias voltage of - 150 V. To amplify the current by 20%, the bias voltage has to be adjusted to roughly -25 V, which is clearly wrong because the curve has a large slope at that point which indicates that there exist strong ghost currents. Upcoming studies will compare the current discrepancy at beam currents ranging from 50 µA to 200 µA. If the discrepancy factor is not the same, the results will suggest there are effects at play other than an unknown scaling factor in the data processing. To the author’s knowledge, the literature contains no comparison between the projected total current from Allison scanner measurements and direct measurements with a Faraday cup. Under- standing reasons for this anomaly would help establish the underlying efficacy of precision Allison 94 010002000300040005000Quadrupole Voltage VQ3 [V]0102030405060Projected Beam Current [A]Current Measured at Faraday CupConventionalImproved Figure 5.15: Screenshot of measured currents on the x-direction Allison scanner on the Artemis beam line of FRIB, plotted against the bias voltage of the Faraday cup on the scanner. The scanner position and voltage were fixed at a phase space coordinate with high density throughout the measurements. scanner measurements. Current is in some sense the lowest order beam moment - the fact that it is not reproduced accurately by Allison scanner measurements is worrisome and unavoidably casts doubt on results on higher-order phase space moments. The author invites the community to collect data on this issue on other machines to potentially gain insight into the anomaly. 5.5 Conclusion We incorporated two important geometric features for Allison scanners into an improved model that extends the conventional analysis. Both effects were modeled with analytic formulae that have been verified by particle simulations. Asymmetric E-dipole placement alters the voltage-to-angle 95 relation for the selected beam, whereas slit thickness (cid:38) slit width introduces additional scraping and requires a data point weight correction that is quadratic in voltage. These effects increase with the angular extent of the beam distribution, and can significantly change results of phase space measurements relative to conventional analysis methods. Applying these corrections to Allison scanner data at the FRIB front end led to (∼10%) changes in beam moments; these corrections were crucial for obtaining consistent results in a quadrupole scan, which in turn provided experimental verification for the improved analysis methods. Detailed knowledge of the phase space area contributing to each data point also allows one to deconvolve overlapping measurements when the device resolution exceeds the spacing between data points. This was demonstrated using synthetic data sets. Python programs incorporating the improved analysis and noise removal for measurement data are made available [35]. They are readily applicable to any Allison scanner given its device geometry to implement accurate analytic models. Tools will be updated as continued improvements are made. An easy-to-use GUI was developed as an FRIB application and is described in Appendix A.5. 96 CHAPTER 6 OPTIMIZED TRANSVERSE PHASE-SPACE MEASUREMENTS WITH BEAM PROFILE MONITORS This chapter discusses how measurements from beam profile monitors can be utilized to obtain information on the transverse phase space of the beam consistent with linear particle dynamics. Linear dynamics over a limited transport length is a reasonable assumption in FRIB downstream of the charge selection system. Examples shown here primarily relate to wire-scanner profile monitors and tomography as employed in the FRIB front end (see Chapter 3 and 4). However, methods presented can be broadly applicable. Section 6.1 reviews results from linear algebra that are relevant to improving measurements from phase space diagnostics. Section 6.2 describes how quadrupole and solenoid scans can be used to measure the beam phase space moments. Statistical methods to quantify the errors in such measurements and an efficient scheme to minimize the errors are given in Sec. 6.3. In addition to obtaining the first and second order moments of the beam distribution, methods of tomography can be employed to reconstruct the phase space in detail. Techniques for tomographic reconstruction of the spatial density profile and 2D phase space are discussed in Sec. 6.4. At the timing of this writing, the analysis in this chapter is under preparation for submission for publication [43]. 6.1 Review of Linear Algebra This section reviews concepts and tools in linear algebra that are relevant to beam diagnostics using profile monitors. The treatment is based upon Ref. [42, 44, 45] which emphasize numerical aspects of the subject. We restrict our discussion to systems of linear equations Ax = b with full column rank, i.e. A ∈ Rm×n where m ≥ n and rank(A) = n. This corresponds to the case where there are at least as many equations as unknowns. In our applications, b consists of measurement data, x consists 97 of the unknowns (e.g. initial conditions upstream or spatial density in a grid) and A describes the linear relationship between measurement data and unknowns. In practice, a sufficient number of measurements is often made in beam diagnostics to over-constrain the system (i.e. make m > n) to allow better accounting for measurement uncertainties. 6.1.1 Least Squares Approximation Generally, an over-determined system of linear equations Ax = b has no exact solution. In that case, it is useful to find the least squares solution that minimizes a quadratic measure of deviation r(cid:124)r for r = Ax − b. The least squares solution x can be found as follows. Require r to be an extremum with ∇x r(cid:124)r = 0 Then: ∇x (cid:0)x(cid:124) (cid:124) A (cid:124)b + b(cid:124)b(cid:1) = 0 Ax − 2x(cid:124) A (cid:124) Ax − 2A 2A (cid:124) A (cid:124)b = 0 Ax = A (cid:124)b Hence, the least squares solution can be solved by inverting A(cid:124)Ax = A(cid:124)b, giving which is called the normal equation, and is defined as the pseudo-inverse of A. Here, B−1 denotes the inverse of a square matrix B and A† denotes the transpose of a (generally non-square) matrix A. For a square matrix B with full rank, the pseudo-inverse becomes the normal inverse, i.e. B† = B−1. 98 (cid:124) x =(cid:0)A A(cid:1)−1 A(cid:1)−1 A† ≡(cid:0)A (cid:124) (cid:124)b, A (cid:124) A , (6.1) (6.2) 6.1.2 Singular Value Decomposition Given a matrix A ∈ Rm×n where m ≥ n, it can always be factorized by the singular value decomposition (SVD) which decomposes A as follows: where U ∈ Rm×m and V ∈ Rn×n are both orthogonal matrices (i.e. UUT = UTU = I), and (6.3) (6.4) A = UΣVT (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) 0 σ1 0 σ2 ... ... 0 0 ... ... 0 0 . . . . . . . . . 0 0 ... . . . σn ... ... 0 0 Σ = ∈ Rm×n is a rectangular diagonal matrix. σi’s are called the singular values of A and they are typically arranged such that σ1 ≥ σ2 ≥ · · · ≥ σn. From AT AV = V ΣTΣ, we observe that V consists of the eigenvectors of AT A. Similarly, U consists of the eigenvectors of AAT. Since rank(AAT) = n, the last m − n columns of U has to be filled in to complete an orthonormal basis. Note that for square matrices, the set of singular values σi(A) and eigenvalues λi(A) are different in general and should not be confused. 6.1.3 Vector Norms and Matrix Norms A vector norm is a function (cid:107).(cid:107) : Rn → R with the three properties: 1. (cid:107)x(cid:107) = 0 =⇒ x = 0 (positive definite) 2. (cid:107)x + y(cid:107) ≤ (cid:107)x(cid:107) + (cid:107)y(cid:107) (satisfies triangle inequality) 3. (cid:107)αx(cid:107) = α (cid:107)x(cid:107) (absolutely homogeneous) 99 In colloquial terms, a vector norm measures the size of a vector. The family of vector p-norms of an n-dimensional vector x are defined as: (cid:32) n i=1 (cid:107)x(cid:107)p = (cid:33)1/p |xi|p . They include the most commonly used vector norms: (cid:113)|x1|2 + |x2|2 + · · · + |xn|2 (cid:107)x(cid:107)1 = |x1| + |x2| + · · · + |xn| (cid:107)x(cid:107)2 = (cid:107)x(cid:107)∞ = max(|x1| ,|x2| , . . . ,|xn|) L1-norm (Manhattan norm) L2-norm (Euclidean norm) L∞-norm (infinity norm) (6.5) (6.6) (6.7) (6.8) Unless otherwise specified, the length or magnitude of a vector in physics is usually synonymous with its Euclidean L2-norm. A matrix norm is a function (cid:107).(cid:107) : Rm×n → R that measures the size of a matrix. It has similar properties to a vector norm: 1. (cid:107) A(cid:107) = 0 =⇒ A = 0 (positive definite) 2. (cid:107) A + B(cid:107) ≤ (cid:107) A(cid:107) + (cid:107)B(cid:107) (satisfies triangle inequality) 3. (cid:107)αA(cid:107) = α (cid:107) A(cid:107) (absolutely homogeneous) A common family of matrix norms are called induced norms, whereby the matrix norm is defined in terms of vector norms as follows: (cid:107) A(cid:107) ≡ max (cid:107)x(cid:107)(cid:44)0 (cid:107) Ax(cid:107) (cid:107)x(cid:107) (6.9) Note that both norms on the R.H.S. are vector norms. The induced norm essentially finds, for a given vector norm, the maximum factor by which a matrix can increase such a vector norm. Two corollaries of the definition are: (cid:107) A(cid:107) = max (cid:107)x(cid:107)=1 (cid:107) Ax(cid:107) (cid:107) Ax(cid:107) ≤ (cid:107) A(cid:107) (cid:107)x(cid:107) ∀ (cid:107)x(cid:107) (cid:44) 0 100 (6.10) (6.11) We will use (cid:107) A(cid:107)2 to denote the induced norm of A corresponding to the L2 vector norm. This is the only matrix norm we will use in the rest of the chapter. As an example of the definition of induced norms, we prove two results below that will be useful in the discussion on condition number in Sec. 6.1.4. Firstly we prove (cid:107) A(cid:107)2 = max (cid:107)x(cid:107)2=1 (cid:107) Ax(cid:107)2 = σ1(A), (6.12) where σ1(A) denotes the largest singular value of A. Proof: (cid:107) Ax(cid:107)2 = = √ √ √ x(cid:124)A(cid:124)Ax x(cid:124)V Σ(cid:124)U(cid:124)UΣV(cid:124)x x(cid:124)V Σ(cid:124)ΣV(cid:124)x max (cid:107)x(cid:107)2=1 = √ (cid:107) Ax(cid:107)2 = max (cid:107)x(cid:107)2=1 = max (cid:107)V(cid:124)x(cid:107)2=1 x(cid:124)V Σ(cid:124)ΣV(cid:124)x √ x(cid:124)V Σ(cid:124)ΣV(cid:124)x (cid:112)y(cid:124)Σ(cid:124)Σy (y = V (cid:124)x) (∵ VV (cid:124) = I) = max (cid:107)y(cid:107)2=1 = max(σi) (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 Next, for the pseudo-inverse A† ≡ (A(cid:124)A)−1 A(cid:124) defined in Sec. 6.1.1, = 1 σn(A) (6.13) where σn(A) denotes the smallest singular value of A. 101 Proof: (cid:124) (cid:124) = V (cid:124) A ΣV (cid:124) A = V Σ A(cid:1)−1 (cid:0)A Σ(cid:1)−1 (cid:124)(cid:0)Σ Σ(cid:1)−1 (cid:124) = V(cid:0)Σ A(cid:1)−1 (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 (cid:13)(cid:13)(cid:13)V(cid:0)Σ Σ(cid:1)−1 (cid:16)(cid:0)Σ Σ(cid:1)−1 = σ1 A = V (cid:124) Σ Σ (cid:124) (cid:124) (cid:124) (cid:124) (cid:0)A (cid:124) U (cid:124) (cid:124) (cid:124)(cid:13)(cid:13)(cid:13)2 (cid:124)(cid:17) U Σ = 1 σn(A) 6.1.4 Condition Number When we solve Ax = b, the solution x will be changed by perturbations both in A and in b, which we denote as δA and δb respectively. The condition number of A, commonly denoted by κ(A), is a number that quantifies the sensitivity of the solution to such perturbations. In our context, b typically represents beam measurements with error δb, and A is constructed from linear transfer maps with uncertainties δA. We first show the results for a special case and follow with the general case. To illustrate the framework, we first consider the special case where δA = 0, A ∈ Rn×n and rank(A) = n. Let x and ˜x be the solutions to: Using: Ax = b A˜x = b + δb Ax = b (cid:107)b(cid:107) = (cid:107) Ax(cid:107) (cid:107)b(cid:107) ≤ (cid:107)x(cid:107) (cid:107) A(cid:107) (cid:107)x(cid:107) ≤ (cid:107) A(cid:107) 1 (cid:107)b(cid:107) 102 (6.14) (6.15) (6.16) δx =(cid:101)x − x (cid:13)(cid:13)(cid:13) (cid:13)(cid:13)(cid:13)A−1 (cid:107)δx(cid:107) ≤(cid:13)(cid:13)(cid:13)A−1(cid:13)(cid:13)(cid:13) (cid:107)δb(cid:107) δx = A−1 δb (cid:107)δx(cid:107) = δb we obtain where κ(A) is the condition number of A defined as: (cid:107)x(cid:107) ≤ κ(A)(cid:107)δb(cid:107) (cid:107)δx(cid:107) (cid:107)b(cid:107) (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 = σ1(A) σn(A) ≥ 1 κ(A) ≡ (cid:107) A(cid:107)2 (6.17) (6.18) (6.19) The relative change in the solution measured by (cid:107)δx(cid:107) /(cid:107)x(cid:107) is bounded by the relative change in the measurement data (cid:107)δb(cid:107) /(cid:107)b(cid:107) multiplied by κ(A). Therefore, in the context of diagnos- tic measurements, the smaller κ(A) is, the less sensitive the results are to relative errors in the measurements. The treatment can be extended to the general case where A ∈ Rm×n and there are perturbations to both A and b. From Sec. 6.1.1, finding the least squares solution to Ax = b is equivalent to solving the normal equation (cid:124) A Ax = A (cid:124)b. (6.20) Applying perturbations to both A and b and leaving only linear terms, the normal equation becomes (cid:124) δA (cid:124) Ax + A (cid:124) δAx + A Aδx = δA (cid:124) (cid:124)b + A δb. (6.21) It can be shown, analogously to Eq. (6.18), that the relative change in x is bounded by the relative changes in A and b multiplied by powers of κ(A): (cid:18) (cid:107)b(cid:107)2 (cid:107) Ax(cid:107)2 (cid:19) (cid:107)δb(cid:107)2 (cid:107)b(cid:107)2 + (cid:107)δA(cid:107)2 (cid:107) A(cid:107)2 (6.22) (cid:107)δx(cid:107)2 (cid:107)x(cid:107)2 ≤ κ(A)2 (cid:107)r(cid:107)2 (cid:107) Ax(cid:107)2 (cid:107)δA(cid:107)2 (cid:107) A(cid:107)2 + κ(A) 103 Proof: (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 (cid:107)(δb − δAx)(cid:107)2 (cid:19) = κ(A) are used (cid:107)δAx(cid:107)2 (cid:107) A(cid:107)2 (cid:107)x(cid:107)2 + (cid:107)δx(cid:107)2 = = (cid:107)δx(cid:107)2 (cid:107)x(cid:107)2 δA δA + (cid:124) (cid:124) (cid:124) (cid:124) (b − Ax) + A† (δb − δAx) δx =(cid:0)A A(cid:1)−1 (cid:13)(cid:13)(cid:13)(cid:0)A (cid:124) (b − Ax)(cid:13)(cid:13)(cid:13)2 (cid:13)(cid:13)(cid:13)A† (δb − δAx)(cid:13)(cid:13)(cid:13)2 A(cid:1)−1 ≤(cid:13)(cid:13)(cid:13)(cid:0)A A(cid:1)−1(cid:13)(cid:13)(cid:13)2 (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 (cid:107)(δb − δAx)(cid:107)2 (cid:13)(cid:13)δA (cid:124)(cid:13)(cid:13)2 (cid:107)(b − Ax)(cid:107)2 + (cid:13)(cid:13)(cid:13)(cid:0)A A(cid:1)−1(cid:13)(cid:13)(cid:13)2 (cid:13)(cid:13)δA (cid:124)(cid:13)(cid:13)2 (cid:107)(b − Ax)(cid:107)2 + (cid:19) (cid:18)(cid:107)δb(cid:107)2 (cid:107) A(cid:107)2 (cid:107) A(cid:107)2 (cid:107) A(cid:107)2 (cid:107) A(cid:107)2 ≤ κ(A)2 (cid:107)δA(cid:124)(cid:107)2 A(cid:1)−1(cid:13)(cid:13)(cid:13)2 (cid:107) A(cid:107)2 where(cid:13)(cid:13)(cid:13)(cid:0)A (cid:13)(cid:13)(cid:13)A†(cid:13)(cid:13)(cid:13)2 (cid:107) A(cid:107)2 (cid:18) (cid:124) (cid:18) (cid:107)δb(cid:107)2 (cid:18) (cid:107)b(cid:107)2 + κ(A) 2 = κ(A)2 and (cid:107) A(cid:107)2 (cid:107)δb(cid:107)2 + κ(A) (cid:107) A(cid:107)2 (cid:107)x(cid:107)2 ≤ κ(A)2 (cid:107)δA(cid:124)(cid:107)2 (cid:107) A(cid:107)2 (cid:107)x(cid:107)2 ≤ κ(A)2 (cid:107)δA(cid:107)2 (cid:107)r(cid:107)2 (cid:107) A(cid:107)2 (cid:107) Ax(cid:107)2 = κ(A)2 (cid:107)r(cid:107)2 (cid:107)δA(cid:107)2 (cid:107) Ax(cid:107)2 (cid:107) A(cid:107)2 (cid:107)r(cid:107)2 (cid:107) A(cid:107)2 + κ(A) + κ(A) (cid:107) Ax(cid:107)2 (cid:107) Ax(cid:107)2 (cid:107)δA(cid:107)2 (cid:107)x(cid:107)2 (cid:107) A(cid:107)2 (cid:107)x(cid:107)2 (cid:107)δb(cid:107)2 (cid:107)b(cid:107)2 (cid:107)δA(cid:107)2 (cid:107) A(cid:107)2 + (cid:124) (cid:107)r(cid:107)2 (cid:107) A(cid:107)2 (cid:107)δAx(cid:107)2 (cid:107) A(cid:107)2 (cid:107) A(cid:107)2 (cid:107) A(cid:107)2 + (cid:107) A(cid:107)2 (cid:19) (cid:19) + 6.2 Quadrupole Scan and Solenoid Scan: Working Principles This section describes a formulation to extract 1st and 2nd order moments of a beam at a location upstream of quadrupoles or solenoids via beam profile measurements with less information at a common downstream location taken under varying focusing strengths. This method is typically called a “quadrupole scan” for quadrupole optics, and “solenoid scan” for solenoid optics. One or more quadrupoles or solenoids may be varied. The method is routinely used to support beam trajectory correction and envelope matching under linear single particle dynamics. It is applicable to any measurements that can measure beam profiles including view screens, wire scanners and beam position monitors. Space charge forces, both linear and nonlinear, are neglected in this formulation because it relies on transfer maps which cannot include dependence on the beam distribution. 104 6.2.1 Notation The notation used in this section are defined as follows. Greek letters α, β, µ, ν denote 4D-phase space indices 1, 2, 3, 4 with x1 = x, x2 = x(cid:48), x3 = y, x4 = y(cid:48). M is the 4 × 4 transverse transfer matrix and (cid:104). . . (cid:105) denotes beam moments, i.e. statistical average over the beam distribution. Subscripts i and f denote the initial and final values with respect to the reference orbit. The final values are to be distinguished from the values measured on the profile monitor, which are denoted by the subscript m. 6.2.2 1st Order Moments The final position and initial phase space coordinates are related by: (cid:11) (cid:11) (cid:10)xβ (cid:10)xβ M1β M3β 4 4 β=1 β=1 (cid:104)x(cid:105) f = (cid:104)y(cid:105) f = i i (6.23) (6.24) If the coordinate system of the profile monitor is misaligned to that of the beam line, the measured x and y centroid positions are different from those with respect to the beam’s reference orbit: (cid:104)x(cid:105)m = (cid:104)x(cid:105) f + xshift = (cid:104)y(cid:105)m = (cid:104)y(cid:105) f + yshift = + xshift + yshift i i (6.25) (6.26) (cid:11) (cid:11) (cid:10)xβ (cid:10)xβ M1β M3β 4 4 β=1 β=1 Therefore, the transverse shifts resulting from misalignments can be treated as unknowns along with the initial first order moments. One obtains two equations from each measurement. With n 105 sets of measurement data, they can be solved from the system of equations (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) Row 1 of M1 1 0 Row 3 of M1 0 1 Row 1 of M2 1 0 Row 3 of M2 0 1 ... Row 1 of Mn 1 0 Row 3 of Mn 0 1 ... ... (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:104)x(cid:105)i (cid:104)x(cid:48)(cid:105)i (cid:104)y(cid:105)i (cid:104)y(cid:48)(cid:105)i xshift yshift = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x(cid:105)1 (cid:104)y(cid:105)1 (cid:104)x(cid:105)2 (cid:104)y(cid:105)2 ... (cid:104)x(cid:105)n (cid:104)y(cid:105)n (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (6.27) In the symbolism of Ax = b, b consists of the measured data and A is the coefficient matrix. The first 4 columns of A depend on the linear transfer map corresponding to each measurement, whereas the last two columns are (1,0,1,0, . . . ,1,0)(cid:124) and (0,1,0,1, . . . ,0,1)(cid:124) respectively. This treatment thus incorporates any shifts induced by misalignment. 6.2.3 2nd Order Moments Denoting: 4 β=1 x f ,α = Mαβxi,β (6.28) The second order moments at the measurement point are related to the initial second order moments as follows: (cid:10)(xα − (cid:104)xα(cid:105))(xµ −(cid:10)xµ (cid:11))(cid:11) (cid:42) 4 4 (cid:11)) Mαβ(xβ −(cid:10)xβ 4 4 (cid:10)(xβ −(cid:10)xβ (cid:11) whose results must cancel. The related centroid subtraction the shift affects both xµ and(cid:10)xµ Note that the equation holds regardless of whether the profile monitor is displaced, because f Mµν(xν − (cid:104)xν(cid:105)) (cid:11))(xµ −(cid:10)xµ (cid:43) (cid:11))(cid:11) i MαβMµν β=1 ν=1 = = β=1 ν=1 i 106 is omitted in the notation below for ease of presentation, in which case final values of 2nd-order spatial moments can be written as: β=1 4 4 4 β=1 ν=1 4 4 4 ν=1 β=1 ν=1 (cid:10)xβxν (cid:10)xβxν (cid:10)xβxν (cid:11) (cid:11) (cid:11) M1βM1ν M1βM3ν M3βM3ν (cid:104)xx(cid:105) f = (cid:104)x y(cid:105) f = (cid:104)yy(cid:105) f = i i i (6.29) (6.30) (6.31) Each profile measurement can give some or all of the above 2nd-order spatial moments. Then, the 10 initial 2nd-order moments can be solved via a system of linear equations Ax = b where: x =(cid:0)(cid:104)xx(cid:105)i,(cid:104)xx(cid:48)(cid:105)i,(cid:104)x y(cid:105)i,(cid:104)x y (cid:48)(cid:105)i,(cid:104)x(cid:48)x(cid:48)(cid:105)i,(cid:104)x(cid:48) y(cid:105)i,(cid:104)x(cid:48) (cid:48)(cid:105)i,(cid:104)yy(cid:105)i,(cid:104)yy (cid:48)(cid:105)i,(cid:104)y (cid:48) (cid:48)(cid:105)i y y (cid:1)(cid:124) and b with n measurements reads: b =(cid:0)(cid:104)xx(cid:105)1 ,(cid:104)x y(cid:105)1 ,(cid:104)yy(cid:105)1 , . . . ,(cid:104)xx(cid:105)n ,(cid:104)x y(cid:105)n ,(cid:104)yy(cid:105)n (cid:1)(cid:124) A is the coefficient matrix that depends on the transfer map for each measurement. 6.2.4 A Limitation of Solenoid Scans For an ideal, aligned solenoid, transport symmetries in linear optics mandate: M11M14 = M12M13 M31M34 = M32M33 M11M34 + M31M14 = M12M33 + M32M13 (6.32) (6.33) (6.34) The coefficient matrix A has 10 columns corresponding to the 10 initial 2nd order moments. However, the columns corresponding to (cid:104)x y(cid:48)(cid:105)i and (cid:104)x(cid:48)y(cid:105)i are always identical. Therefore, only the sum of two terms, (cid:104)x y(cid:48)(cid:105)i + (cid:104)x(cid:48)y(cid:105)i can be solved but not each term individually. 107 6.3 Quadrupole Scan and Solenoid Scan: Error Quantification and Mini- mization 6.3.1 Error Quantification In quadrupole and solenoid scans, the beam matrix is solved using a linear system of equations Ax = b. The error in the solution x can be quantified via statistical methods. There are many sources of errors in a quadrupole / solenoid scan including: misalignments, field deviations and measurement errors. These errors constitute perturbations to A and b in Ax = b. Given estimates of the magnitude of each error, one can apply random errors and solve for Ax = b in each case. Frequency count of the solutions give rise to a probability distribution of x. With the additional condition that unphysical answers are discarded, this treatment is equivalent to a Bayesian analysis with a uniform prior on the domain of physicality. 6.3.2 Error Minimization From the discussion on condition number in Sec. 6.1.4, it is clear that the error of quadrupole and solenoid scan results can be minimized by choosing a set of scan parameters that give a coefficient matrix A with as small a condition number as possible. The condition number κ(A) is defined in Eq. (6.19). This requirement also requires having transport without beam loss in each measurement. The use of condition number in error minimization of beam measurements is inspired by the work at the GSI Helmholtz Centre for Heavy Ion Research (GSI) [46, 47]. However, the papers by GSI did not appear to analyze how to choose scan parameters such that the condition number is minimized. Better choice of scan parameters can reduce the condition number and thereby decrease the relative uncertainty of results. A simple estimate can show that an exhaustive search over all possible sets of scan parameters is impractical. Suppose there are two knobs (e.g. focusing strenghs in a quadrupole doublet) where each knob can attain 10 values, thus giving 100 possible settings in total. To choose a set of scan 108 parameters for four measurements, the number of possible combinations equal: (cid:169)(cid:173)(cid:173)(cid:171)100 4 (cid:170)(cid:174)(cid:174)(cid:172) = 100! 4! × 96! ≈ 4 × 106 . For scans with more knobs and more measurements, the number can be orders of magnitude larger. It would be very computationally inefficient to build all possible coefficient matrices, compute their singular values, and select the one with the smallest condition number. Therefore, one has to rely on other ideas to efficiently obtain a set of scan parameters. The Paul Scherrer Institute (PSI) proposed choosing quadrupole parameters that correspond to discrete steps in particle phase advance between the measurement and reconstruction point [48]. However, the phase advance describes rotation in normal coordinates, whereas the actual rotation in phase space depends on the orientation of the invariant ellipses, which is different for each focusing setting and may not even exist (in a single-pass system, there is no need for a cell to satisfy the stability condition). Building upon PSI’s idea, we believe the viewpoint from tomography is the most direct and beneficial. Quadrupole settings alter the linear map from zi to z f and in effect provide a different projection of the initial phase space onto the 1D spatial measurements of the beam profile monitor. If we choose a set of measurements which correspond to a diverse range of projection angles, the corresponding matrix A should have a low condition number which reduces uncertainties in the solution. The projection angle corresponding to the linear map of a quadrupole transport line can be found as follows. The initial and final phase space coordinates are related by: Consider a line in x f -x(cid:48) is given by: f phase space defined by x f = x0 where x0 in a constant. Its image in xi-x(cid:48) i x(cid:48) f m21 m22 (cid:169)(cid:173)(cid:173)(cid:171)x f (cid:170)(cid:174)(cid:174)(cid:172) =(cid:169)(cid:173)(cid:173)(cid:171)m11 m12 (cid:169)(cid:173)(cid:173)(cid:171)xi (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:172) =(cid:169)(cid:173)(cid:173)(cid:171) m22 −m12 (cid:169)(cid:173)(cid:173)(cid:171)x0 (cid:169)(cid:173)(cid:173)(cid:171)xi (cid:170)(cid:174)(cid:174)(cid:172) −m21 m11 x(cid:48) f x(cid:48) i x(cid:48) i (cid:170)(cid:174)(cid:174)(cid:172) 109 (6.35) (6.36) Figure 6.1: Schematic of the beam line section containing the first profile monitor at the FRIB Front End. which are parametric equations of xi and x(cid:48) in xi-x(cid:48) is given by i with x(cid:48) f as the parameter. Hence the slope of the line i tan θ = −m11 m12 (6.37) where θ corresponds to the projection angle on the initial phase space. The argument is illustrated by Fig. 6.1. Here we employ normalized dimensionless coordinates in measuring x and x(cid:48) (e.g. normalization by 1 mm and 1 mrad respectively). We must add one caveat concerning the dimensionless coordinates in the aforementioned treatment. The projection angle is actually dependant on the choice of units for the phase space coordinates. For instance, in Fig. 6.1, if the unit of angle is changed from mrad to rad, the plot is compressed 1000-fold in y and the angle θ will decrease drastically. Thus we must ask: what units should we choose for the phase space coordinates? The most natural answer is a choice of units such that the beam has roughly equal sizes in the two dimensions. Heuristically, it often suffices to use the units adopted at a certain facility, e.g. mm for position and mrad for angles at FRIB, because the distribution typically has similar extent in both coordinates in these units. This is likely the reason why such units are chosen in the first place. In a more sophisticated treatment, two measurements with projection angles close to 0 and π/2 respectively, should be made first. They give a good estimate of the beam size in the original phase space and the associated projection angles will change minimally under a change of units. Then, 110 xix'ix'fxfMM-1θ the phase space units can be rescaled in accordance with the beam extent (measured by rms size or maximum extent), and the projection angle corresponding to each setting can be recalculated under the new units. Suppose mm and mrad are used as units of position and angle respectively. The new projection angles will be given by: α = x[mm] x(cid:48)[mrad] tan θnew = α tan θ (6.38) (6.39) The choice of settings for other measurements can now be selected as before, with a diverse spread of projection angles (under the rescaled coordinates) being the primary criterion. Techniques explored in this section and their connection to the condition number κ are presented in FRIB examples in Sec. 6.5. 6.4 Spatial and Phase Space Tomography 6.4.1 Spatial Tomography Spatial tomography harnesses information from beam profile monitors to reconstruct the beam’s image in x y space. One approach is the algebraic reconstruction technique (ART) which solves for densities on a grid as unknowns in a system of linear equations. A projection consists of many steps where each individual step is called a ray. A ray can be visualized as a line of finite thickness that samples the distribution and threads a path through the grid. For an m × n grid, the system of equations reads: where ρj is the discretized current density in the j-th grid cell, Ii is the measured current of the i-th ray, and A is the projection matrix with ai j being the area of the j-th grid cell that falls under the i-th ray. 111 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) a1,1 a2,1 ... ak,1 . . . a1,m×n . . . a2,m×n ... . . . ak,m×n ... (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) ρ1 ρ2 ... ρm×n (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) I1 I2 ... Ik (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) = (6.40) Figure 6.2: Schematic of a wire scanner at the FRIB front end. Image courtesy of Tomofumi Maruta at FRIB. In FRIB and several other facilities, the wires in beam profile monitors are oriented vertically, horizontally and diagonally. A schematic of a wire scanner device at FRIB is shown in Fig. 6.2; it consists of two diagonal wires and one vertical wire. If we attempt to apply ART to conduct spatial tomography with such a system, we quickly run into problems due to severe redundancy in the system of equations. For an m× n grid, it can be observed that horizontal and vertical wires only provide information on the total current in each row and column, thereby giving at most m and n equations respectively. This is true regardless of how dense the steps are in each projection. While a diagonal wire can couple information from different rows and columns, it is not effective and only provide no more than m + n − 1 equations. This can be seen by coloring the grid as one would a chessboard. There are m + n − 1 diagonals of squares that alternate in color, and the wire can only give information on the total current in each of the diagonals. The argument above is illustrated in Fig. 6.3. If a diagnostic station contains vertical, horizontal and two diagonal wires, for a m × n grid, the total number of equations < 4m + 4n − 2. The number is not equal to 4m + 4n − 2 due to further redundancies among the information provided by individual wires (e.g. the total current). In all cases of interest, 4m + 4n − 2 is much smaller than mn which is the number of equations required for solving the current densities in all grid cells. 112 xy Figure 6.3: A diagonal wire only gives the total current in each diagonal because and each grid cell along a diagonal always contributes equally to a ray. The alternating diagonals are plotted using chessboard coloring on the right. Figure 6.4: A wire intersecting a 2 : 1 grid at 45o is equivalent to the case of a wire intersecting an even grid at 63o. The most direct solution would be to change the physical orientation of the wires. However, this is costly to do, and, perhaps more importantly, the existing orientation is optimal for measuring beam moments (cid:104)xx(cid:105), (cid:104)x y(cid:105) and (cid:104)yy(cid:105). Beam moments should take precedence over spatial tomography. The other remedies are indirect ways of altering the wire angles. The first method is to use grids with uneven horizontal and vertical sizes. This amounts to a rescaling of one of the axis, and a diagonal that intersects an uneven grid at 45o is equivalent to the a wire intersecting a regular grid at a different angle. An example is shown in Fig. 6.4. However, this method only removes the redundant information from diagonal wires, not horizontal and vertical, and the tunability is low because a grid with vastly uneven sizes in the two directions is unlikely to be desirable. A more effective and implementable method is to rotate the grid on which the densities are to be solved by making a judicious choice of transformed coordinaces. The coordinate system is 113 133355224446rescaleuneven grid (1:2)even grid (1:1)45o63o Figure 6.5: For a rotated grid, the wires are no longer vertical, horizontal or diagonal with respect to it. an abstract construction anyway and does not have to be aligned with the axis of the beam line coordinates. This allows wires that are horizontal, vertical and diagonal with respect to the beam line to intersect the grid at angles of one’s choice. Fig. 6.5 shows a grid that is rotated by 25o with respect to the standard one. The ART technique for spatial tomography was tested using synthesized measurements from several numerically generated particle distributions. Both uniform and waterbag distributions are standard, while hollow distributions have been observed at FRIB and dual peak will be relevant when FRIB operates dual-species transport. The tests assumed the profile monitor has three wires: one vertical, one horizontal and one diagonal. The reconstruction results for both standard and rotated grids are plotted in Fig. 6.6. One can observe that the rotated grid produces more truthful reconstructions for the uniform and waterbag distributions. For hollowed and dual peak distributions, it may be beneficial to employ a variety of rotated coordinates and to combine their results. Optimized choices for such treatment will require further study. 6.4.2 Phase Space Tomography Unlike spatial tomography, phase space tomography can be conducted with as many projections as one desires. Hence, the reconstruction can be conducted using a standard grid, and the choice of focusing parameters between the profile monitor and reconstruction point can be chosen with the 114 Standard GridRotated Grid Figure 6.6: Spatial tomography of several synthesized distributions. 115 UniformWaterbagHollowDual PeakParticle DistributionStandard GridRotated Grid Figure 6.7: Schematic of the beam line section containing the first profile monitor at the FRIB Front End. The profile monitor is located at the position designated by f and the beam was reconstructed at position i. Green blocks denote identical (Q7 type) electrostatic quadrupoles. V1 V2 V3 V4 -2657 4513 -4295 300 -2657 4513 -4295 1300 -2657 4513 -4295 2300 -2657 4513 -4295 3300 -2657 4513 -4295 4300 -2657 4513 -4295 5300 Table 6.1: Quadrupole parameters for Scan 1 same reasoning as Sec. 6.3.2. Examples of phase space tomography at FRIB are shown in Sec. 6.5. 6.5 Example: FRIB Measurements To test the implementation of the methods developed in this chapter, two quadrupole scans were conducted using the first profile monitor at the FRIB Front End on 2019/01/18. A schematic of the relevant beam line section is shown in Fig. 6.7. There are two quadrupole doublets between the measurement and reconstruction points. The voltages applied to the four quadrupoles are listed in Table 6.1 for Scan 1 and Table 6.2 for Scan 2. The parameters in Table 6.1 is typical of how quadrupole scans were usually performed where only the focusing strength of the quadrupole immediately upstream of the profile monitor is varied. The parameters in Table 6.2 were chosen using the techniques described in Sec. 6.3.2 to sample the initial distribution from a wide range of projection angles. The projection angles on the initial phase space corresponding to each scan are shown in Fig. 6.8. Note that Scan 2 has a much wider spread of angles than Scan 1. The performance of Scan 1 and Scan 2 in both beam matrix measurements and phase space tomography are compared. The phase space distribution measured at the upstream Allison scanner 116 -500 3000 -4500 4000 -4000 4500 -4500 4000 -4500 -4000 V1 5500 4500 V2 -3000 -4500 V3 500 3500 V4 Table 6.2: Quadrupole parameters for Scan 2 -1500 4000 -4500 4000 Figure 6.8: Projection angles on the initial phase space corresponding to quadrupole scan parameters in Table 6.1 and Table 6.2. is propagated to the reconstruction location to serve as a benchmark. Hard edge equivalent transfer matrices are applied for the quadrupoles. As discussed in Sec. 6.2, beam moments are calculated from a system of linear equations where the coefficient matrix A depends on the choice of scan parameters. In the current example, the condition number κ(A) as defined in Sec. 6.1.4 equals 466 and 14.5 for Scan 1 and Scan 2 respectively. With this > 30-fold difference in the condition number, the errors in Scan 1 are expected to be much larger than those in Scan 2. This has been verified by applying random errors to measurement results and solving for the emittance in each case. The histograms corresponding to 10000 perturbed solutions in each case are plotted in Fig. 6.9, where the standard deviation of x-emittance εx is 0.029 mm-mrad and 0.005 mm-mrad for Scan 1 and Scan 2 respectively. Despite the fact that applied errors for Scan 1 are chosen to be 10 times smaller than those applied for Scan 2, the measurement errors in Scan 1 are still 6 times larger. This strongly illustrates the viability of the scan methodology presented in Sec. 6.3.2 to reduce uncertainties in quadrupole scans. 117 −4−2024xi [mm]−4−2024x′i [mrad]Measurement Set 1−4−2024xi [mm]−4−2024Measurement Set 2 (a) Scan 1 (b) Scan 2 Figure 6.9: Histograms of the normalized x-emittance with random errors applied to the measurement results. The applied errors have a truncated Gaussian distribution where 3σ = 1% for Scan 1 and 3σ = 10% for Scan 2. Figure 6.10: Comparison between tomographic reconstruction of the x-x(cid:48) phase space distribution, and the measured distribution at the Allison scanner propagated to the same location. The fact that the initial beam conditions are better reconstructed in Scan 2 than in Scan 1 also manifest itself in tomographic reconstruction. Fig. 6.10 compares the reconstructed phase space distribution from Scan 1 and Scan 2 against the distribution propagated from Allison scanner measurements. The tomography results from Scan 1 are very noisy, whereas the reconstructed phase space of Scan 2 agrees much more closely with Allison scanner measurements. The beam moments and emittances calculated from the reconstructed phase space distribution are listed in Table 6.3. The algorithms for spatial tomography with rotated grids presented in Sec. 6.4.1 was also tested as part of the experiment. Fig. 6.11 shows the reconstructed x-y spatial density distribution of the beam on a grid that is rotated by 20o. Beam moments of the reconstructed phase space are 118 0.000.020.040.060.080.100.120.140.16εx [mm-mrad]0100200300400500600700FrequencyAllison ScannerUnperturbed Solution0.060.070.080.090.10εx [mm-mrad]025050075010001250150017502000FrequencyAllison ScannerUnperturbed SolutionPhase Space Tomography from Measurem Set 1Phase Space Tomography from Measurem Set 2Propagated Distribution from Allison Scanner Measurements Tomography Tomography Scan 1 Scan 2 Projected from Allison Scanner (cid:104)xx(cid:105) [mm2] (cid:104)xx(cid:48)(cid:105) [mm-mrad] (cid:104)x(cid:48)x(cid:48)(cid:105) [mrad2] εx normalized [mm-mrad] 21.4 56.3 220.0 0.199 6.49 22.9 132.7 0.107 7.09 21.4 105.6 0.087 Table 6.3: Beam moments and emittances corresponding to the phase space distributions plotted in Fig. 6.3. Figure 6.11: Spatial tomography with a grid rotated by 20o clockwise with respect to the standard x-y coordinate system. Wire Scanner Measurements Rotated Grid Tomography (cid:104)xx(cid:105) [mm2] (cid:104)x y(cid:105) [mm2] (cid:104)yy(cid:105) [mm2] 44.2 24.7 44.6 46.9 24.6 51.3 Table 6.4: Comparison between beam moments of the reconstructed spatial distribution in Fig. 6.11 and the values obtained from wire scanner measurements. transformed back into x-y and compared against moments directly obtained from beam profile measurements with the wire scanner. Table 6.4 shows the results where the agreement is quite close, thus validating the algorithms as a potential expansion of the diagnostic capabilities of wire scanners at FRIB. 119 6.6 Further Work We have developed error minimization techniques and phase space tomography algorithms for quadrupole scans with profile monitors and tested them using the wire scanner data at FRIB. Further quadrupole scan experiments should be performed which can serve multiple purposes: 1. Test the method of making projects on rescaled phase space outlined at the end of Sec. 6.3.2. 2. Develop the above method into a general automated algorithm for online selection of optimal quad scan parameters. 3. Further testing and development of automated 2D phase space tomography as an additional output of quad scans. 4. Use the first profile monitor to benchmark beam moments and reconstructed phase space from quadrupole scans against Allison scanner measurements propagated to the same location. We would like to emphasize the significance of item 4, where the beam’s 2D phase space roughly half way between the Allison scanner and 1st profile monitors is obtained using distinct methods and compared. Good agreement in such a benchmark will lend confidence to the analysis methods of both devices as well as the transfer maps employed. Two extensions on the error minimization techniques will be highly beneficial. Firstly, the choice of optimal quad scan parameters in Sec. 6.3.2 only applies to 2D phase spaces x-x(cid:48) and y-y(cid:48). The four coupled moments are also measured by the quadrupole scan, and they are likely even more susceptible to measurement errors because they have one more unknown than 2D phase spaces. An understanding of error minimization for coupled moments should be the next step in advancing the physical model behind quadrupole scans. An even more difficult problem is how error minimization techniques can be applied to solenoid scans, where nine coupled moments are solved simultaneously (not ten, see Sec. 6.2.4). One plausible idea is to generalize the projection angle picture in 2D and apply it to 4D. For a given set of solenoid strengths, each measured moment correspond to some projection on the initial 4D phase 120 space, which can be represented by a point on a unit 3-sphere. However, with three angles instead of one, the idea of a diverse range of projection angles becomes much more fuzzy. Even if one associates a “spread” in projection angles by inter-point geodesic distances, it is not obvious how an optimal set can be chosen efficiently. It may be possible to reduce dimensionality by employing rotated frame Larmor variables or other choices of canonical variables in pure solenoid scans. That may enable the formulation of a procedure analogous to that developed for quadrupole scans. Much investigation remains to be done, and we hope such a pursuit may also provide a physical picture on why (cid:104)x y(cid:48)(cid:105) and (cid:104)x(cid:48)y(cid:105) cannot be measured separately in a solenoid scan. 121 APPENDICES 122 APPENDIX A EXTENSIONS ON ALLISON SCANNER ANALYSIS In this appendix, we include details associated with the Allison scanner analysis (Chapter 5) including: estimates on the range of validity of the paraxial approximation employed (A.1); a numerical model used to check analytic results (A.2); key steps in the derivation of analytic formulae (A.3); noise removal schemes for experimental measurement data (A.4); and implementation of the algorithms in the FRIB control system (A.5). A.1 Non-paraxial Effects: When They are Negligible This appendix estimates the magnitude of non-paraxial effects in an Allison scanner and provides a simple condition to determine whether they are negligible. Consider an idealized device geometry, i.e. symmetric E-dipole placement and thin slits, which corresponds to Case I in Fig. 5.2. Given V0, a reference particle with the corresponding angle x(cid:48) ref reaches maximum x−displacement at the axial center (z = l + L/2). The ratio between the axial kinetic energy decrement ∆E ≡ E − Emin at this point and the original axial kinetic energy E of the ion measures the strength of non-paraxial effects. Taking uniform hard-edge dipole fields and the paraxial approximation,  For x(z = 0) = 0 and x(cid:48)(z = 0) = x(cid:48) ref = 1 2 qV0L gE , the particle position at the center of the device is x(cid:48)(cid:48)(z) = − qV0 gE if l ≤ z ≤ L + l 0 otherwise . (cid:19) (cid:18)3L qV0 gE L qV0 gE L2 . 8 + l 1 x(z = l + L/2) = 2 ≈ 1 5 The energy change ∆E is maximum at z = l + L/2 where the particle has its closest approach to 123 the E-dipole. The fractional energy change is ∆E E = 1 E x(z = l + L/2) g/2 qV0 = (cid:19)2 (cid:18) L g (cid:19)2(cid:18) qV0 E . 2 5 For non-paraxial effects to be negligible, such energy gain should be much smaller than the kinetic energy of the beam, i.e. ∆E (cid:28) E, or (cid:19)2 2 5 (cid:19)2(cid:18) qV0 (cid:18) L (cid:19)2(cid:18) (cid:18)Q E g (cid:28) 1. (cid:19)2 40 A eV0 E0 [eV/u] (cid:28) 1. Taking the approximation L/g ≈ 10 for a typical plate spacing, this condition reduces to Here, for ions we take E0 = E/A to be the kinetic energy per atomic mass unit and q = Qe with Q being the charge state and e the elementary charge. Non-paraxial effects are negligible at V0 values that satisfy the above condition. For example, consider the condition applied to the E0 = 12 keV/u 40Ar9+ ion beam measured in the FRIB front end in Sec. 5.4. Even at V0 = 1000 V, the condition reduces to ≈ 0.014 (cid:28) 1, which is well satisfied. Thus it is not surprising that the realistic and ballistic simulations agree well for typical parameters. A.2 Realistic Model A Python [49] code is employed to numerically integrate particle equations of motion in a realistic x-z 2D field map of the device geometry that is generated from the electrostatic field code POISSON [24]. Ex = −∂φ/∂x and Ez = ∂φ/∂z electric field data are exported from POISSON onto a high-resolution, uniform x-z mesh (dx = dz = 0.2 mm), and then imported into the Python code. Fields at the particle position are calculated using bilinear interpolation from the gridded field data [50]. In Fig. A.1, the applied field potential φ of the FRIB Allison scanner is contoured showing enhanced detail near the entrance slit. Only half the geometry in the zoomed figure is contoured since φ(−x, z) = −φ(x, z). 124 Figure A.1: Potential contours of the FRIB Allison scanner in Fig. 5.1, including details of the fringe structure in the vicinity of the entrance slit (right). Particles are advanced with non-relativistic equations of motion d2x dt2 = d2z dt2 = q m q m Ex, Ez. The independent variable can be transformed exactly from time t to axial coordinate z giving (cid:16) qEx m x(cid:48)(cid:17) x(cid:48) t(cid:48) m − qEz m t(cid:48)3 − qEz t(cid:48)2  d dz where (cid:48) ≡ d/dz, and t(cid:48) = 1/vz.   x t x(cid:48) t(cid:48)  = (A.1) This numerical model includes full fringe field effects entering and exiting the dipole field region, as well as non-paraxial effects due to energy change as the particle crosses potential lines. Image charges, beam space charge and scattering effects are neglected. The state vector in Eq. (A.1) describing the particle trajectory is advanced using the ODE package within Scientific Python (SciPy)[51] for specified initial particle coordinate x, angle x(cid:48), and dipole voltage V0 (field data scaled). The code takes into account scraping on all boundaries. To solve for x(cid:48) max,min, note that the corresponding trajectory must touch the slits at two points (see Fig. 5.2). We employ a numerical root-finding procedure to solve for the initial x(cid:48) that connects the upstream point to the downstream point. x(cid:48) ref is solved analogously with the condition x = 0 at both ends. 125 z [cm]x [mm]Electric Dipole Plate (V = V0)Electric Dipole Plate (V = -V0)x [mm]z [mm]0-5-10-15510150123456701234560123456 (a) (b) Figure A.2: How (a) shifting a trajectory for fixed angle and (b) changing the angle for fixed initial position allow one to calculate quantities listed in Table 5.2. A.3 Analytic Results Appendix A.3.1 sketches key steps in the derivation of the expressions in Table 5.2 using the geometric models shown in Fig. 5.2. Complications involved in Case IV are outlined in Appendix A.3.2. A.3.1 Sketch of Derivation Using Geometric Models Table 5.2 is mostly derived based on two simple rules: 1) trajectories can be shifted in x-position; and 2) changing the initial x(cid:48) = x(cid:48) at z = zi by δx(cid:48) results in a displacement δx = (z − zi)δx(cid:48) i downstream. Subscripts i and f denote initial (entrance) and final (exit) locations respectively. In Fig. A.2, these principles are applied to Case III (symmetric E-dipole placement, thick slits) to illustrate several calculations. We observe that among a uniform spatial distribution of particles entering the slit with x(cid:48) ref, a fraction s1/s is collimated by the slit plate due to its thickness, where s1 = x(cid:48) ref), which manifests itself in Fig. 5.3 as a horizontal shrinkage of the blue area in Case III in comparison with Case I. in the transmission coefficient T(x(cid:48) refd. This explains the factor (cid:19) (cid:18) 1 − x(cid:48) refd s To calculate x(cid:48)max, we note that the corresponding trajectory touches the slit plates at the positions shown in Fig. A.2. Consider a particle with the same xi entering the slits with angle x(cid:48) ref. Applying the second rule, we obtain s2 = (x(cid:48)max − x(cid:48) refd as calculated ref)(d + 2l + L), while s1 = x(cid:48) 126 Thick slitlldds1Thick slitlldd (a) (b) Figure A.3: Reversing particle velocity allows one to check results when the E-dipole placement is asymmetric. above. Therefore, s2 = s − s1, ref)(d + 2l + L) = s − x(cid:48) refd, s − x(cid:48) refd d + 2l + L , x(cid:48)max − x(cid:48) (x(cid:48)max − x(cid:48) ref = which is one of the results in Table 5.2. Note that the derivation above assumes bending within the E-dipole, which does not hold when V0 = 0. In that case, instead of ∆x(cid:48) = 2s/(L + 2l + d), one can draw straight trajectories between furthest corners of the slits to see that ∆x(cid:48) = 2s/(L + 2l + 2d). Since d (cid:28) L + 2l, the difference is very small and will be neglected. (cid:44) −x(cid:48) f The calculations for asymmetric E-dipole placement were carried out using the same principles, but one must note that x(cid:48) due to the asymmetry. This makes calculations in Case IV i (asymmetric E-dipole placement, thick plate) much more complicated due to different scraping factors on the two ends. Rather than showing the tedious calculations, it is interesting to observe a succinct way to check the validity of the expressions in Case IV of Table 5.2 via direction reversal. Fig. A.3a shows Case IV trajectories and Fig. A.3b their direction-reversed counterparts, where the trajectories corresponding to maximum and minimum angles are interchanged depending on the direction of incoming particles. For the x(cid:48) = x(cid:48)max trajectory in the normal case, the final particle angle equals x(cid:48)max − qV0L min for the direction-reversed case. Therefore, we can gE , which equals −x(cid:48) 127 check whether (cid:19) (cid:18) L + 2l2 (cid:32) L + 2l∗ ˆL + 2 ˆL qV0L gE (cid:34)1 1 2 = − qV0L gE 2 (cid:32) (cid:33) s − x(cid:48) refd ˆL + d − 1 ˆL + d 1 + l∗ (cid:33)(cid:35) − qV0L gE s − L + 2l∗ 1 L + 2l∗ 2 2 remains constant and l∗ x(cid:48) refd where inter-slit distance ˆL = L + l1 + l2 = L + l∗ 2 = l1. All results for which the E-dipole placement is asymmetric have been verified with these procedures. 1 = l2, l∗ A.3.2 Results for Case IV Procedures sketched above are applied to Case IV to obtain, for l2 > l1: ref ≤ x(cid:48) for x(cid:48) for ˜x(cid:48) ≤ x(cid:48) < x(cid:48) ref for x(cid:48) < ˜x(cid:48) , 2(l2 − l1)d (L + l1 + l2 + 2d)(L + 2l2) (cid:21) ,  T = W = D1 c1 x(cid:48)max−x(cid:48) x(cid:48)max−x(cid:48) ref x(cid:48) ref−x(cid:48) ref− ˜x(cid:48) c2 + x(cid:48)− ˜x x(cid:48) x(cid:48) ref− ˜x(cid:48) c1 ˜x(cid:48)−x(cid:48) ˜x(cid:48)−x(cid:48) min L + l1 + l2 c2 2s2 , (cid:33) (cid:32) (cid:18) (cid:20) 1 − 1 − x(cid:48) refd s 1 − ˜x(cid:48)d (c1s)2 + (c2s)2 2(L + l1 + l2 + d) + (cid:19) s , , ˜x(cid:48) = x(cid:48) ref c1 = c2 = where D1 = ref − ˜x(cid:48)). ref only if l1 (cid:44) l2 and d (cid:44) 0, so the case ˜x(cid:48) ≤ x(cid:48) < x(cid:48) 2(c1 + c2)(x(cid:48) s Observe that ˜x(cid:48) (cid:44) x(cid:48) the E-dipole placement is asymmetric and slits are thick. ref is only relevant when The reason there exists a third region in T(x(cid:48)) is associated with a subtlety in the definition of x(cid:48) ref. Previously, x(cid:48) ref is defined as the angle at which particles enter and exit the slits at the same x-position, with no distinction on which side of the slits we mean when the slits are thick. 128 Figure A.4: Comparison between transmission ratio in Case III and IV when weight correction effect is large. Taking the inner facing side of the entrance slit as z = 0, x(cid:48) ref is defined as the angle at which particles attain the same position at z = 0 and z = ˆL where ˆL = L + l1 + l2 is the inter-slit distance. Then there is another angle ˜x(cid:48) with which the particles attain the same position at z = −d and z = ˆL + d. x(cid:48) ref = ˜x(cid:48) in Case III because the entering and exiting angles are the same. However, in Case IV, since the entering and exiting angles are different due to the asymmetry, x(cid:48) ref (cid:44) ˜x(cid:48). When the d/s ratio is large, the effects can be significant as shown in Fig. A.4 for d/s (cid:39) 4. A.4 Noise Removal Scheme The noise removal scheme for raw data is illustrated using an example measurement from the FRIB Front End (see Sec. 5.4). The scheme first defines the beam region and uses the data points outside the beam region to characterize the background. Then the background is subtracted and unrealistic islands in the distribution are filtered. Details of the procedure are discussed below. Such noise thresholding is crucial for heavy ion beam measurements because the total beam current is typically ∼50 µA; this leaves many data points with merely ∼ nA of collected currents which may only be ∼10 times larger than noise fluctuations. Algorithms implementing the noise removal scheme applied in this study are incorporated in the Python programs available in Ref. [35]. The noise removal scheme employed here has no capability to correct for ghost signals caused by particles that scatter after impacting the E-dipole plates; such ghost signals have been mitigated 129 −0.20−0.15−0.10−0.050.000.050.100.150.20x′−x′ref [mrad]0.0000.0250.0500.0750.1000.1250.1500.1750.200T(x′−x′ref) x′ref̃x′̃ase IIĨase IV Figure A.5: Beam region is defined by the red ellipse. Figure A.6: Histogram of current values at data points outside the defined beam region. by using E-dipole plates whose surfaces have a staircase profile to reduce grazing incidences [33, 52]. A.4.1 Specifying the Beam Region Figure A.5 shows an ellipse that surrounds the entire beam distribution, thereby defining a region that contains the beam. Such a region can be designated by user input or generated automatically. It should be self-consistent in the sense that, after background subtraction and island filter, no non-zero data point touches the boundary specifying the beam region. 130 −10−505101520Position x [mm]−40−2002040Angle x′ [mrad]−1.0−0.50.00.51.01.52.02.53.0−1.0−0.9−0.8−0.7−0.6Current I [nA]0100200300400500600700800Count+2σ−2σ Figure A.7: Non-zero data points after background subtraction. A.4.2 Background subtraction After designating the beam region, all exterior data points are used to characterize the background. In this example, the background noise has an average µ = −0.778 nA with standard deviation σ = 0.079 nA. Fig. A.6 shows a histogram of the measured currents in exterior data points, where almost all values lie below µ + 2σ. Therefore, a µ + 2σ cutoff is typically applied for interior data points, whereby any point with current < µ + 2σ is regarded as pure noise and assigned the value 0. Subsequently, all remaining non-zero data points have µ subtracted from them to correct for the background. Note that µ < 0 in this case, which is likely due to amplifier characteristics. Fig. A.7 shows all non-zero data points after applying background subtraction to the raw data in Fig. A.5. A.4.3 Island Filter The beam distribution after background subtraction often contains spurious islands that probably arise from noise that exceeds the cutoff value. Fig. A.8 zooms in on the upper right sector of the beam distribution in Fig. A.7 to show a number of islands before and after they are filtered. The island filter we employ is a modification of the median filter widely used in speech and image processing [53]. For each non-zero data point, the island filter examines neighbors in an n x n grid and counts the number of non-zero points. If the number < k, the data point in question is assigned a zero value after the filter has processed all data points. The filter can be applied 131 −10−505101520Position x [mm]−40−2002040Angle x′ [mrad] Figure A.8: Non-zero data points of the beam distribution, before and after applying the island filter. until it no longer has any effect on the distribution, which may require several iterations to enable layer-by-layer removal of clusters of islands. In Fig. A.8, n = 5 and k = 3, which means that if fewer than 3 data points amongst 24 neighboring points are non-zero, the data point at center is deemed part of a noise island and assigned the value zero. The filter takes two iterations to reach convergence in this case. A.5 Implementation in FRIB Control System The improved analysis of Allison scanner measurements described in this dissertation and Refs. [31, 34] have been incorporated into applications within the FRIB Control System. The initial application was made available to members of the Accelerator Physics Department in the form of a Jupyter notebook [54] which can be used both online and offline. The online mode provides the functionality of setting up and starting a scan using the notebook. Once the measurement data is ready, it is read from the control system, saved and analyzed with the algorithm. The data is saved in such a way that they can be readily analyzed by running the notebook offline. Tong Zhang from the Acceleartor Physics Department at FRIB incorporated the algorithms into a graphical user interface (GUI) that eases usage. The interface is set up to allow enhanced configurability while being straightforward to use with defaults for operators not familiar with 132 5.07.510.012.515.017.520.0Position x [mm]2025303540Angle x′ [mrad]Unfiltered5.07.510.012.515.017.520.0Position x [mm]Filtered physics issues. and prepares it for transition into a tool for the operators. The new software also communicates better with the control system and allows real-time update of the data and scanner status throughout a scan. Two screenshots of the software are shown in Fig. A.9 and Fig. A.10. 133 Figure A.9: Beam parameters setup screen for the Allison scanner GUI. (Image courtesy of Zhang Tong from the Accelerator Physics Department at FRIB). 134 Figure A.10: Noise correction screen for the Allison scanner GUI. (Image courtesy of Zhang Tong from the Accelerator Physics Department at FRIB). 135 APPENDIX B ALTERNATIVE PROOF OF RELATIONS GOVERNING 2ND ORDER MOMENTS OF Cn BEAMS This appendix presents an alternative, more direct, proof of the fact that, for all n > 2, the 2nd order moments of a beam with n-fold rotational symmetry (i.e. Cn) obey the same relations as those of an axisymmetric beam, i.e.: (cid:104)xx(cid:105) = (cid:104)yy(cid:105) (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:48)(cid:105) (cid:48) (cid:104)x y(cid:105) = 0 (cid:48)(cid:105) = 0 (cid:104)x(cid:48) y (cid:48)(cid:105) = −(cid:104)x(cid:48) (cid:104)x y y y(cid:105) (B.1) This is true regardless of the orientation of the coordinate axis on the transverse plane. Methods outlined here would become increasingly cumbersome when applied to higher-order moments. B.1 Spatial Density and Local Velocity Define F(r, θ, vr, vθ) as the beam distribution function in 4D transverse phase space where r and θ are the usual cylinderical-polar spatial coordinates and vr and vθ denote the corresponding angles in a cylindrical polar representation. The usual transverse particle angles (i.e. transverse velocity normalized by vz, x(cid:48) = vx/vz) are related to vr and vθ by: x(cid:48) = vr cos θ − vθ sin θ (cid:48) = vr sin θ + vθ cos θ y 136 A conveniently normalized spatial distribution function f(r, θ) withÞ 2π 0 Þ ∞ 0 f(r, θ)rdrdθ = 1 is defined as: f(r, θ) = Þ 2π 0 Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ −∞ −∞ 0 −∞ F(r, θ, vr, vθ)dvr dvθ −∞ F(r, θ, vr, vθ)rdrdθdvr dvθ The local average flow angle of particles within the distribution is denoted (cid:174)V(r, θ) = Vr(r, θ)ˆr + Vθ(r, θ) ˆθ (components vary locally as a function of r and θ) and is obtained by averaging over angular degrees of freedom: Þ ∞ Þ ∞ Þ 2π Þ ∞ −∞ 0 (cid:174)V(r, θ) = Þ ∞ Þ ∞ −∞ F(r, θ, vr, vθ)(vr ˆr(θ) + vθ 0 ˆθ(θ))dvr dvθ −∞ F(r, θ, vr, vθ)rdrdθdvr dvθ −∞ For a beam with n-fold rotational symmetry, its spatial density and local flow angle of the particle distribution must remain unchanged upon a coordinate rotation by angle φ = 2jπ/n ∀j ∈ Z, i.e. f(r, θ) = f(r, θ + Vr(r, θ) = Vr(r, θ + Vθ(r, θ) = Vθ(r, θ + 2jπ ) n 2jπ ) n 2jπ ) n (B.2) (B.3) (B.4) A sufficient, but not necessary, condition for these relations to hold is the even stronger statement that F(r, θ, vr, vθ) satisfies F(r, θ, vr, vθ) = F(r, θ + 2jπ n , vr, vθ) ∀j ∈ Z In that case, Eqs. (B.2)–(B.3) can be verified from: f(r, θ) = = Þ 2π Þ 2π 0 0 Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ −∞ −∞ −∞ −∞ 2jπ ) n 0 0 = f(r, θ + −∞ F(r, θ, vr, vθ)dvr dvθ −∞ F(r, θ, vr, vθ)rdrdθdvr dvθ −∞ F(r, θ + n , vr, vθ)dvr dvθ 2jπ n , vr, vθ)rdrdθdvr dvθ −∞ F(r, θ + 2jπ 137 and Vr(r, θ) = ˆr(θ) · (cid:174)V(r, θ) −∞ 0 Þ ∞ Þ ∞ Þ ∞ Þ ∞ Þ ∞ −∞ F(r, θ, vr, vθ)vr dvr dvθ Þ ∞ Þ ∞ −∞ F(r, θ, vr, vθ)rdrdθdvr dvθ Þ ∞ Þ ∞ Þ ∞ −∞ −∞ F(r, θ + 2nπ3 , vr, vθ)vr dvr dvθ −∞ F(r, θ + 2nπ3 , vr, vθ)rdrdθdvr dvθ −∞ 2nπ 3 ) −∞ Þ 2π Þ 2π 0 = = 0 0 = Vr(r, θ + (B.5) Similarly, Vθ(r, θ) = Vθ(r, θ + 2nπ3 ). B.2 Useful Trigonometric Identities In this section, several trigonometric identities are proved in anticipation of their utility in subsequent arguments. The first identity is: (cid:18) n j=1 (cid:18) n j=1 sin (cid:19) cos θ + 4jπ n = θ + 4jπ n = 0 ∀n > 2 (B.6) (cid:19) (cid:19)(cid:21) Identity 1 can be efficiently proven in a complex number representation with Euler’s theorem n eix = cos x + i sin x with i ≡ √−1: (cid:20) (cid:18) (cid:17) (cid:44) 1, this implies that: For exp(cid:16) 4πi exp j=1 θ + i n (cid:20) (cid:18) exp i θ + 4jπ n (cid:19)(cid:21) (cid:18)4πi n (cid:19) n (cid:19)(cid:21) j=1 4jπ n = 0 4jπ n = exp (cid:18) (cid:20) i θ + exp n j=1 The fact that both the real and imaginary parts of the L.H.S. must vanish individually gives us Eq. (B.6). Note that when n = 2, the proof does not work because exp(4πi/2) = 1. 138 Next we define sums: j=1 to prove two corollaries of Identity 1 as: (cid:19) (cid:18) θ + 2mπ n θ + sin j=1 S2 ≡ n S3 ≡ n S4 ≡ n j=1 cos (cid:18) cos2(cid:18) sin2(cid:18) θ + θ + (cid:19) (cid:19) (cid:19) 2jπ n 2jπ n 2jπ n S2 = 0 ∀n > 2 S3 = S4 ∀n > 2 (B.7) (B.8) (B.9) (B.10) (B.11) Identity 2 can be proven using the trigonometric identity sin x cos x = 1 2 sin(2x) to rewrite S2 as (cid:19) (cid:18) (cid:18) S2 ≡ n n j=1 1 2 = j=1 sin 2θ + 4jπ n (cid:19) (cid:18) sin (cid:19) cos θ + 2jπ n θ + 2jπ n It then follows from Identity 1 with θ → 2θ that S2 = 0 ∀n > 2. Similarly, Identity 3 can be proven using the trigonometric identities cos2 x = 1 2 + 1 2 cos(2x) 139 and sin2 x = 1 2 − 1 2 cos(2x) and adding in zero symmetrically: S3 ≡ n 1 2 j=1 n 2 + cos2(cid:18) n n n 2 + 0 2 − 1 sin2(cid:18) n n j=1 j=1 2 = = = = j=1 ≡ S4 θ + 2jπ n (cid:19) (cid:19) (cid:18) (cid:18) θ + 2jπ n cos 2θ + 4jπ n by Identity 1 4jπ n 2θ + cos (cid:19) (cid:19) by Identity 1 Showing that S3 = S4 ∀n > 2. B.3 Proof of Beam Moment Relations In this section, we prove the 2nd order axisymmetric beam moment relations in Eq. (B.1) hold for a beam with Cn symmetry by the following trick. For such a beam, all moments remain the same when the coordinate system rotates by 2jπ/n where j is an integer. By summing over j = 1 to j = n, the trigonometric identities derived in Sec. B.2 can be applied to simplify resulting trigonometric functions and show that the 2nd order axisymmetric moment constraints hold. All relations are proven for completeness, but with decreasing levels of detail in each successive sub-proof as the underlying manipulations are analogous. Proof of (cid:104)x y(cid:105) = 0: To prove (cid:104)x y(cid:105) = 0, we note that (cid:104)x y(cid:105) = Þ 2π Þ ∞ 0 0 r2 cos θ sin θ f(r, θ)rdrdθ 140 By symmetry, for all integer j, Using Eq. (B.12), we can write: (cid:19) (cid:19) )rdrdθ f(r, θ + 2jπ n f (r, θ) rdrdθ (B.12) (cid:19) f (r, θ) rdrdθ (cid:19) (cid:19) sin sin (cid:18) (cid:18) (cid:19) 2jπ n 2jπ n (cid:18) θ + θ + 2jπ n 2jπ n (cid:18) (cid:18) (cid:18) θ + θ + r2 cos r2 cos 0 0 Þ 2π Þ ∞ Þ 2π Þ ∞ Þ 2π Þ ∞ n Þ ∞ Þ 2π 0 0 0 0 r2 cos θ + 2jπ n sin θ + 2jπ n r2S2 f(r, θ)rdrdθ 0 0 by Identity 2 (cid:104)x y(cid:105) = = n(cid:104)x y(cid:105) = (cid:104)x y(cid:105) = j=1 1 n = 0 Proof of (cid:104)xx(cid:105) = (cid:104)yy(cid:105): (cid:104)xx(cid:105) = = (cid:104)yy(cid:105) = = 0 0 1 n Þ 2π Þ ∞ Þ ∞ Þ 2π Þ 2π Þ ∞ Þ ∞ Þ 2π 0 0 0 0 1 n 0 0 r2 cos2 θ f(r, θ)rdrdθ r2S3 f(r, θ)rdrdθ r2 sin2 θ f(r, θ)rdrdθ r2S4 f(r, θ)rdrdθ from Identity 3: S3 = S4 ∴ (cid:104)xx(cid:105) = (cid:104)yy(cid:105) 141 Proof of (cid:104)xx(cid:48)(cid:105) = (cid:104)yy(cid:48)(cid:105): 0 0 1 n Þ ∞ Þ 2π Þ ∞ Þ 2π Þ ∞ Þ 2π Þ ∞ Þ 2π 0 0 0 0 1 n 0 0 (cid:104)xx(cid:48)(cid:105) = = (cid:104)yy (cid:48)(cid:105) = = r cos θ [Vr(r, θ) cos θ − Vθ(r, θ) sin θ] f(r, θ)rdrdθ [S3Vr(r, θ) − S2Vθ(r, θ)] f(r, θ)r2drdθ r sin θ [Vr(r, θ) sin θ + Vθ(r, θ) cos θ] f(r, θ)rdrdθ [S4Vr(r, θ) + S2Vθ(r, θ)] f(r, θ)r2drdθ from Identity 3: from Identity 2: S3 = S4 S2 = 0 ∴ (cid:104)xx(cid:48)(cid:105) = (cid:104)yy (cid:48)(cid:105) Henceforth, we will suppress angular arguments of Vr(r, θ), Vθ(r, θ), and f(r, θ) to further abbreviate. Proof of (cid:104)x y(cid:48)(cid:105) = −(cid:104)x(cid:48)y(cid:105): (cid:104)x y (cid:48)(cid:105) = = (cid:104)x(cid:48) y(cid:105) = = 0 0 1 n Þ 2π Þ ∞ Þ 2π Þ ∞ Þ 2π Þ ∞ Þ ∞ Þ 2π 0 0 0 0 1 n 0 0 r cos θ [Vr sin θ + Vθ cos θ] f rdrdθ r [S2Vr − S3Vθ] f rdrdθ r sin θ [Vr cos θ − Vθ sin θ] f rdrdθ r [S2Vr + S4Vθ] f rdrdθ from Identity 3: from Identity 2: S3 = S4 S2 = 0 (cid:48)(cid:105) = −(cid:104)x(cid:48) y(cid:105) ∴ (cid:104)x y 142 Proof of (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y(cid:48)y(cid:48)(cid:105): (cid:104)x(cid:48)x(cid:48)(cid:105) = = (cid:48) (cid:104)y (cid:48)(cid:105) = y = 0 0 1 n Þ 2π Þ ∞ Þ ∞ Þ 2π Þ ∞ Þ 2π Þ ∞ Þ 2π 0 0 0 0 1 n 0 0 (cid:104) (cid:104) (cid:105) (cid:105) f rdrdθ f rdrdθ [Vr cos θ − Vθ sin θ]2 f rdrdθ r − 2S2VrVθ + S4V2 S3V2 θ [Vr sin θ + Vθ cos θ]2 f rdrdθ r + 2S2VrVθ + S3V2 S4V2 θ from Identity 3: from Identity 2: S3 = S4 S2 = 0 ∴ (cid:104)x(cid:48)x(cid:48)(cid:105) = (cid:104)y (cid:48) (cid:48)(cid:105) y Proof of (cid:104)x(cid:48)y(cid:48)(cid:105) = 0: (cid:104)x(cid:48) (cid:48)(cid:105) = y = Þ ∞ Þ 2π Þ 2π Þ ∞ 0 (cid:110) 0 1 n 0 0 (cid:104) [Vr cos θ − Vθ sin θ][Vr sin θ + Vθ cos θ] f rdrdθ (cid:105) + VrVθ [S3 − S4](cid:111) f rdrdθ S2 V2 r + V2 θ from Identity 3: from Identity 2: ∴ (cid:104)x(cid:48) y (cid:48)(cid:105) = 0 S3 = S4 S2 = 0 143 APPENDIX C GENERATING A 4D COUPLED DISTRIBUTION In this appendix, we present a method for generating an arbitrary 4D coupled distribution using a similarity transformation. Recall the 4D sigma matrix is given by: (cid:104)xx(cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:104)xx(cid:105) (cid:104)xx(cid:48)(cid:105) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) σ = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)x y(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)yy(cid:105) (cid:104)yy(cid:48)(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:104)yy(cid:48)(cid:105) (cid:104)y(cid:48)y(cid:48)(cid:105) (cid:16) y(cid:48)(cid:17). Since σ is a real symmetric matrix, ∃ orthogonal matrix R (i.e. =(cid:10)xx(cid:124)(cid:11) (C.1) where x(cid:124) = x(cid:48) R(cid:124) = R−1) such that: x y where λj, j = 1,4 denote the four eigenvalues of σ. The method can be described as follows. One first generates a 4D completely uncoupled distribution in the transformed coordinates ˜x which satisties: (cid:10)˜x˜x(cid:124)(cid:11) = Then, transforming ˜x in accordance with: λ1 0 0 0 0 λ2 0 0 0 0 λ3 0 0 0 0 λ4 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) R(cid:124) σR = λ1 0 0 0 0 λ2 0 0 0 0 λ3 0 0 0 0 λ4 (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) ˜x (cid:55)→ x = R˜x 144 (C.2) (C.3) (C.4) will result in a 4D coupled distribution with the specified moments. Proof. σ =(cid:10)xx(cid:124)(cid:11) (cid:104)xx(cid:48)(cid:105) (cid:104)x(cid:48)x(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) 0 0 0 λ2 0 0 λ3 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) (cid:104)x y(cid:48)(cid:105) (cid:104)x(cid:48)y(cid:48)(cid:105) (cid:104)yy(cid:48)(cid:105) (cid:104)y(cid:48)y(cid:48)(cid:105) (cid:104)x y(cid:105) (cid:104)x(cid:48)y(cid:105) (cid:104)yy(cid:105) (cid:104)yy(cid:48)(cid:105) 0 0 0 (cid:170)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:174)(cid:172) λ4 R(cid:124) = (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) (cid:104)xx(cid:105) (cid:104)xx(cid:48)(cid:105) (cid:104)x y(cid:105) (cid:104)x y(cid:48)(cid:105) (cid:169)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:173)(cid:171) λ1 0 0 0 = R = R(cid:10)˜x˜x(cid:124)(cid:11) R(cid:124) =(cid:10)R˜x˜x(cid:124)R(cid:124)(cid:11) Note that no assumption is made on the form of the decoupled distribution in ˜x, insofar as the 2nd-order moments satisfy Eq. (C.3). However, how distributions in ˜x should be interpreted after the transformation may warrant careful considerations (see the discussion at the end of Sec. 3.4). This formulation can be applied to generate a wide variety of coupled distributions with spec- ified 2nd-order moments to simulate magnetized beams emerging from ECR sources. Additional symmetry arguments (see Chapter 4) together with physics considerations may also be exploited to reduce the number of free parameters. 145 BIBLIOGRAPHY 146 BIBLIOGRAPHY [1] [2] J. Wei, “The Very High Intensity Future”, in Proceedings, 5th International Particle Acceler- ator Conference (IPAC 2014), Dresden, Germany, June 15-20, 2014 (2014), MOYBA01. J. Wei et al., “The FRIB Superconducting Linac - Status and Plans”, in Proceedings, 28th International Linear Accelerator Conference (LINAC16), East Lansing, USA, September 25-30, 2016 (2017), pp. 1–6. [3] T. Xu et al., “FRIB Cryomodule Design and Production”, in Proceedings, 28th International Linear Accelerator Conference (LINAC16), East Lansing, Michigan, September 25-30, 2016 (2017), WE2A02. [4] F. Marti, P. Guetschow, Y. Momozaki, J. Nolen, and C. Reed, “Development of a Liquid Lithium Charge Stripper for FRIB”, in Proceedings, 13th Heavy Ion Accelerator Technology Conference (HIAT2015), Yokohama, Japan, September 7-11, 2015 (2016), TUA1I01. [5] G Machicoane, D Cole, J Ottarson, J Stetson, and P Zavodszky, “ARTEMIS-B: A room- temperature test electron cyclotron resonance ion source for the National Superconducting Cyclotron Laboratory at Michigan State University”, Review of Scientific Instruments 77, 03A322 (2006). [6] D Leitner, C. Lyneis, S. Abbott, D Collins, R. Dwinell, M. Galloway, M Leitner, and D. Todd, “Next generation ECR ion sources: First results of the superconducting 28 GHz ECRIS–VENUS”, Nuclear Instruments and Methods in Physics Research B 235, 486–493 (2005). [7] S. Lidia et al., “Overview of Beam Diagnostic Systems for FRIB”, in Proceedings, 4th Inter- national Beam Instrumentation Conference (IBIC2015), Melbourne, Australia, September 13-17, 2015 (2016), MOPB071. [8] R Geller, Electron cyclotron resonance ion sources and ECR plasmas (IOP, Bristol, 1996). [9] D. Leitner, M. L. Galloway, T. J. Loew, C. M. Lyneis, I Castro Rodriguez, and D. S. Todd, “High intensity production of high and medium charge state uranium and other heavy ion beams with VENUS”, Review of Scientific Instruments 79, 02C710 (2008). [10] J. D. Jackson, Classical Electrodynamics, 3rd Edition (July 1998), p. 832. [11] R. D. Ruth, “Single particle dynamics and nonlinear resonances in circular accelerators”, in Lecture Notes in Physics (Springer, 1986), pp. 37–63. [12] S. Y. Lee, Accelerator Physics (WSP, 2019). 147 [13] H. Wiedemann, Particle accelerator physics; 4th ed. (Springer, Berlin, 2015). [14] H. Goldstein, C. Poole, and J. Safko, Classical mechanics (2002). [15] A. J. Dragt, F. Neri, and G. Rangarajan, “General moment invariants for linear hamiltonian systems”, Phys. Rev. A 45, 2572–2585 (1992). [16] M. Berz, K. Makino, and W. Wan, An introduction to beam physics, Series in high energy physics, cosmology, and gravitation (CRC Pr., Boca Raton, 2015). [17] J. Perl, J. Shin, J. Schumann, B. Faddegon, and H. Paganetti, “TOPAS: An innovative proton Monte Carlo platform for research and clinical applications”, Medical Physics 39, 6818 (2012). [18] Warp, http://warp.lbl.gov/. [19] D. P. Grote, A. Friedman, J.-L. Vay, and I. Haber, “The WARP Code: Modeling High Intensity Ion Beams”, in AIP Conference Proceedings, Vol. 749, 1 (2005), pp. 55–58. [20] S. M. Lund, T. Kikuchi, and R. C. Davidson, “Generation of initial kinetic distributions for simulation of long-pulse charged particle beams with high space-charge intensity”, Phys. Rev. ST Accel. Beams 12, 114801 (2009). [21] CST Studio Suite, https://www.cst.com/products/csts2. [22] MAD, http://mad.web.cern.ch/mad/. [23] S. Lund, E. Pozdeyev, H. Ren, and Q. Zhao, “RMS Emittance Measures for Solenoid Trans- port and Facitlity for Rare Isotope Beams Front-End Simulations”, in Proceedings, 54th ICFA Advanced Beam Dynamics Workshop on High-Intensity and High-Brightness Hadron Beams (HB2014), East Lansing, Michigan, November 10-14, 2014 (2015), MOPAB17. [24] M. Menzel and H. K. Stokes, User’s guide for the POISSON/SUPERFISH group of codes, tech. rep. LA-UR-87-115 (Los Alamos National Lab., NM (United States), 1987). [25] V. Mironov, S. Bogomolov, A. Bondarchenko, A. Efremov, and V. Loginov, “Numerical model of electron cyclotron resonance ion source”, Phys. Rev. ST Accel. Beams 18, 123401 (2015). [26] V. Mironov, S. Bogomolov, A. Bondarchenko, A. Efremov, and V. Loginov, “Numerical simulations of gas mixing effect in electron cyclotron resonance ion sources”, Phys. Rev. Accel. Beams 20, 013402 (2017). [27] D. Winklehner, “Ion beam extraction from electron cyclotron resonance ion sources and the subsequent low energy beam transport”, PhD thesis (Michigan State University, 2013). 148 [28] J. C. Wong and S. M. Lund, “Beams with Rotational Symmetry”, Phys. Rev. Accel. Beams (forthcoming). [29] Z. He, M. Davidsaver, K. Fukushima, D. Maxwell, G. Shen, Y. Zhang, and Q. Zhao, “De- velopment Status of FRIB On-line Model Based Beam Commissioning Application”, in Proceedings, 28th International Linear Accelerator Conference (LINAC16): East Lansing, Michigan, September 25-30, 2016 (2017), MOPRC015. [30] A. W. Chao and M. Tigner, Handbook of Accelerator Physics and Engineering (3rd Printing) (World Scientific Publishing Co, 1999). [31] J. C. Wong, S. M. Lund, and T. Maruta, “High resolution phase space measurements with Allison-type emittance scanners”, Phys. Rev. Accel. Beams 22, 072801 (2019). [32] P. W. Allison, J. D. Sherman, and D. B. Holtkamp, “An emittance scanner for intense low-energy ion beams”, IEEE Transactions on Nuclear Science 30, 2204–2206 (1983). [33] M. Stockli, R. Welton, R Keller, and M Leitner, “Emittance studies with an Allison scanner”, Review of Scientific Instruments 77, 03B706 (2006). [34] C. Y. J. Wong et al., “Physics Model of an Allison Phase-Space Scanner, with Application to the FRIB Front End”, in Proceedings, 6th International Beam Instrumentation Conference (IBIC2017), Grand Rapids, USA, August 20-24, 2017 (2018), pp. 72–76. [35] https://github.com/AllisonScanner/Analysis. [36] M. Stockli, W Blokland, T Gorlov, B Han, C Long, T Pennisi, and S Assadi, “Low-energy emittance studies with the new SNS Allison emittance scanner”, in Proceedings, 23rd Particle Accelerator Conference (PAC’09), Vancouver, Canada, May 4-8, 2009 (2009), pp. 3974– 3976. [37] R. D’Arcy, M. Alvarez, J. Gaynier, L. Prost, V. Scarpine, and A. Shemyakin, “Characterisation of the PXIE Allison-type emittance scanner”, Nuclear Inst. and Methods in Physics Research, A A815, 7–17 (2016). [38] C. Ullmann, R. Berezov, J. FIls, R. Hollinger, V. Ivanova, O. Kester, W. Vinzenz, N. Chauvin, and O. Delferrière, “Status and Computer Simulations for the Front End of the Proton Injector for Fair”, in Proceedings, 5th International Particle Accelerator Conference (IPAC 2014), Dresden, Germany, June 15-20, 2014 (2014), pp. 1120–1122. [39] A. Laxdal, F. Ames, R. Baartman, W. Rawnsley, A. Sen, V. Verzilov, G. Waters, R. Hari- wal, M. Kownacki, and R. Paris, “Allison Scanner Emittance Diagnostic Development at TRIUMF”, in Proceedings, 27th Linear Accelerator Conference (LINAC2014), Geneva, Switzerland, August 31-September 5, 2014 (2014), pp. 829–833. 149 [40] M. Reiser, Theory and design of charged particle beams (John Wiley & Sons, Hoboken, 2008). [41] https://git-scm.com. [42] J. W. Demmel, Applied numerical linear algebra, Vol. 56 (Siam, 1997). [43] J. C. Wong and S. M. Lund, “Error Minimization in Transverse Phase Space Measurements Using Quadrupole and Solenoid Scans”, Proceedings, 4th North American Particle Acceler- ator Conference (NAPAC2019): Lansing, Michigan, September 1-6, 2019 (forthcoming). [44] G. H. Golub and C. F. Van Loan, Matrix computations, Vol. 3 (JHU press, 2012). [45] L. N. Trefethen and D. Bau III, Numerical linear algebra, Vol. 50 (Siam, 1997). [46] C Xiao, L Groening, P Gerhard, M Maier, S Mickat, and H Vormann, “Measurement of the transverse four-dimensional beam rms-emittance of an intense uranium beam at 11.4 MeV/u”, Nuclear Instruments and Methods in Physics Research Section A 820, 14–22 (2016). [47] C Xiao, M Maier, X. Du, P Gerhard, L Groening, S Mickat, and H Vormann, “Rotating system for four-dimensional transverse rms-emittance measurements”, Phys. Rev. Accel. Beams 19, 072802 (2016). [48] E. Prat and M. Aiba, “Four-dimensional transverse beam matrix measurement using the multiple-quadrupole scan technique”, Phys. Rev. ST Accel. Beams 17, 052801 (2014). [49] Python, https://www.python.org/. [50] C. K. Birdsall and A. B. Langdon, Plasma physics via computer simulation (CRC press, 2004). [51] SciPy, https://www.scipy.org/. [52] M. Stockli, M Leitner, D. Moehs, R Keller, and R. Welton, “Beam dumping ghost signals in electric sweep scanners”, in AIP Conference Proceedings, Vol. 763, 1 (AIP, 2005), pp. 145– 158. [53] J. W. Tukey, Exploratory data analysis, Vol. 2 (Reading, Mass., 1977). [54] T. Kluyver, B. Ragan-Kelley, F. Pérez, B. E. Granger, M. Bussonnier, J. Frederic, K. Kel- ley, J. B. Hamrick, J. Grout, S. Corlay, et al., “Jupyter Notebooks-a publishing format for reproducible computational workflows.”, in ELPUB (2016), pp. 87–90. 150