UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A novel optimization platform and its applications to the TRIUMF energy recovery linac Gong, Chris 2015

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2016_february_gong_chris.pdf [ 8.84MB ]
JSON: 24-1.0221373.json
JSON-LD: 24-1.0221373-ld.json
RDF/XML (Pretty): 24-1.0221373-rdf.xml
RDF/JSON: 24-1.0221373-rdf.json
Turtle: 24-1.0221373-turtle.txt
N-Triples: 24-1.0221373-rdf-ntriples.txt
Original Record: 24-1.0221373-source.json
Full Text

Full Text

A Novel Optimization Platform andIts Applications to the TRIUMFEnergy Recovery LinacbyChris GongB.Sc., The University of British Columbia, 2008A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFDOCTOR OF PHILOSOPHYinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2015© Chris Gong 2015AbstractA novel software platform for global optimization was developed to create abaseline design for the TRIUMF Energy Recovery Linac (ERL). The plat-form is parallel capable, scalable, and allows flexible combinations of variousaccelerator tracking tools such as Madx and Free Electron Laser (FEL) toolssuch as Genesis. The TRIUMF machine includes simultaneously a two-passERL and a rare isotope line. Many parameters are coupled, including RFand the separator system which are shared for all three linac passes. Theglobal optimization platform can study dynamic relationships between dif-ferent processes, a practice not easily performed with piecewise optimization.The FEL induced energy spread, which grows by an order of magnitude afterdeceleration and increases the difficulty of beam disposal, creates a tradeoff,or Pareto front, between the gain and the dump energy spread. Anotherfront forms between energy recovery and final energy spread due to RF set-tings. The Pareto fronts give insights on how objectives are related and therepercussions of design decisions. Pareto relationships are presented, alongwith potential lattice solutions found by the optimization platform.Chris GongiiPrefaceThis dissertation is based on the creation of a computational apparatus andresulting data of the TRIUMF Energy Recovery Linac optimization. Thework is the effort of myself and Y.C. Chao.The computational apparatus is a multiobjective optimizer that worksin a massively parallel environment (WestGrid), and can integrate and linkmultiple modeling tools and codes, i.e. engines, for a flexible method ofcreating global models of particle accelerators. All codes for linking differ-ent engines together are performed by the optimizer without need for userintervention.The majority of the platform codebase deals with handling multiple en-gines, transition between engines, and distribution of work to massivelyparallel computing clusters. These portions of the code are wholly designedand implemented by myself. The genetic optimization algorithm used is animplementation of SPEA2 by I. Bazarov. The code is refactored to suit theneeds of the platform. The fundamental algorithm is unchanged.The modeling engine used to simulate the accelerating cavities is theEmpirical Model (EM). The idea for EM was proposed by Y.C. Chao. TheEM code was wholly written by myself, originally in Java for TRIUMF highlevel applications development. It was ported to C++, also by myself, foroptimization purposes. Interpolation tables used by EM were extracted byY.C. Chao.The setup of the optimization problem involving 62 free parameters, 13objectives, and 3 constraints, presented in chapters 3 and 4 are my ownwork, with guidance and discussions with Y.C. Chao.The optimization results and analysis presented in chapter 5, as wellas the analysis of the accelerator baseline solution presented in chapter 6,are my original work. This includes studying the tradeoffs between differentmachine parameters and presenting an accelerator baseline solution completewith layout coordinates of all optical components in TRIUMF standardformat.All images and tables used in this dissertation were produced by myself.Publications arising from the work presented in this dissertation:iiiPrefaceC. Gong and Y.C. Chao, “A New Platform for Global Optimization,”IPAC12, New Orleans, 2012.C. Gong and Y.C. Chao, “The TRIUMF Optimization Platform andApplication to the E-Linac Injector,” ICAP12, Rostock, 2012.C. Gong and Y.C. Chao, “A Novel Optimization Platform and Its Appli-cations to the TRIUMF Energy Recovery Linac,” ICAP15, Shanghai, 2015.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiList of Programs . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . xiiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 Description of the TRIUMF E-Linac and ERL . . . . . . . 93 Description of the Global Optimization Platform and En-gines Used for ERL Modeling . . . . . . . . . . . . . . . . . . 174 Modeling the ERL . . . . . . . . . . . . . . . . . . . . . . . . . 265 Optimization Results and Tradeoff Studies . . . . . . . . . . 446 ERL Baseline Design . . . . . . . . . . . . . . . . . . . . . . . 897 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100vTable of ContentsAppendicesA Notations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109B Software Design of the Optimization Platform . . . . . . . 111C List of Parameters and Constants in the ERL Optimization 148D ERL Major Components and Layout . . . . . . . . . . . . . 154viList of Tables1.1 Pareto dominance example . . . . . . . . . . . . . . . . . . . 74.1 ERL modeling and topology information . . . . . . . . . . . . 304.2 RIB modeling and topology information . . . . . . . . . . . . 334.3 Dipole modes in the 9-cell cavity . . . . . . . . . . . . . . . . 426.1 ERL Baseline parameters . . . . . . . . . . . . . . . . . . . . 90C.1 Optimization parameters - ERL beam . . . . . . . . . . . . . 148C.2 Optimization parameters - rare isotope beam . . . . . . . . . 149C.3 Optimization parameters - RF . . . . . . . . . . . . . . . . . 150C.4 Optimization parameters - optics . . . . . . . . . . . . . . . . 151C.5 Optimization parameters - undulator . . . . . . . . . . . . . . 153D.1 ERL beam parameters . . . . . . . . . . . . . . . . . . . . . . 154D.2 RIB parameters . . . . . . . . . . . . . . . . . . . . . . . . . . 154D.3 FEL parameters . . . . . . . . . . . . . . . . . . . . . . . . . 155D.4 RF requirements . . . . . . . . . . . . . . . . . . . . . . . . . 155D.5 ERL quadrupole requirements . . . . . . . . . . . . . . . . . . 156D.6 RIB quadrupole requirements . . . . . . . . . . . . . . . . . . 157D.7 ERL dipole requirements . . . . . . . . . . . . . . . . . . . . . 158D.8 ERL element coordinates . . . . . . . . . . . . . . . . . . . . 159viiList of Figures1.1 Pareto front example . . . . . . . . . . . . . . . . . . . . . . . 61.2 Pareto front example (forming new front) . . . . . . . . . . . 72.1 ERL/RIB operations . . . . . . . . . . . . . . . . . . . . . . . 122.2 FEL operating principle . . . . . . . . . . . . . . . . . . . . . 142.3 Energy recovery operating principle . . . . . . . . . . . . . . . 153.1 Longitudinal phase space of Genesis slices . . . . . . . . . . . 213.2 ASTRA vs Empirical Model tracking . . . . . . . . . . . . . . 233.3 Validation of Empirical Model tracking results against ASTRA 254.1 ERL layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . 314.2 ERL optimization topology . . . . . . . . . . . . . . . . . . . 324.3 Layout of the main linac . . . . . . . . . . . . . . . . . . . . . 334.4 Layout of the first arc transport . . . . . . . . . . . . . . . . . 344.5 M56 in the first transport arc . . . . . . . . . . . . . . . . . . 364.6 Longitudinal phase space manipulation in the first transportarc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.7 Layout of ERL dump transport . . . . . . . . . . . . . . . . . 384.8 Longitudinal phase space manipulation in the second trans-port arc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 394.9 Layout of RIB transport . . . . . . . . . . . . . . . . . . . . . 415.1 Beam parameters before and after arc 1 . . . . . . . . . . . . 465.2 Effects of phase on beam parameters . . . . . . . . . . . . . . 475.3 Effects of acceleration phase on peak current . . . . . . . . . 485.4 Effects of acceleration phase on gain . . . . . . . . . . . . . . 495.5 M56 vs peak current, earlier generation . . . . . . . . . . . . . 505.6 Effects of phase on longitudinal parameter . . . . . . . . . . . 515.7 Peak current vsM56, no longitudinal correlation in initial bunch 535.8 Energy loss through lasing . . . . . . . . . . . . . . . . . . . . 545.9 Increase in energy spread due to lasing . . . . . . . . . . . . . 55viiiList of Figures5.10 Effects of gain on beam disposal . . . . . . . . . . . . . . . . 565.11 Gain vs energy recovery . . . . . . . . . . . . . . . . . . . . . 575.12 Energy deviation after lasing . . . . . . . . . . . . . . . . . . 585.13 Tracking beam loss after lasing . . . . . . . . . . . . . . . . . 595.14 Effects of deceleration on energy recovery and beam disposal 605.15 Effects of phase change on energy recovery and beam disposal 615.16 Tradeoff between energy recovery and energy spread . . . . . 625.17 Tradeoff between energy recovery and EDBT beam size . . . 635.18 Effect of RF phases and energy recovery . . . . . . . . . . . . 655.19 RF acceleration and deceleration curves . . . . . . . . . . . . 675.20 Choosing RF phase and time-of-flight for energy matching . . 685.21 Gain vs peak current . . . . . . . . . . . . . . . . . . . . . . . 695.22 Undulator matching conditions from global optimization . . . 725.23 Undulator matching conditions (2D) . . . . . . . . . . . . . . 735.24 Undulator matching conditions from local optimization . . . . 745.25 Quads to match beam into undulator . . . . . . . . . . . . . . 755.26 Beam size after linac pass 1 . . . . . . . . . . . . . . . . . . . 765.27 Effect of EHATQ1 on vertical beam size . . . . . . . . . . . . 765.28 EHATQ1 vs max dump size . . . . . . . . . . . . . . . . . . . 775.29 Beam size control in arc 1 vs EDBT . . . . . . . . . . . . . . 785.30 M56 vs vertical beta function symmetry . . . . . . . . . . . . 795.31 Constraint on βy symmetry . . . . . . . . . . . . . . . . . . . 795.32 Vertical beta function symmetry vs peak current . . . . . . . 805.33 Higher order effects in the linac . . . . . . . . . . . . . . . . . 815.34 Peak current vs M56 . . . . . . . . . . . . . . . . . . . . . . . 825.35 Measuring strengths of T566 . . . . . . . . . . . . . . . . . . . 835.36 Effects of T566 on gain . . . . . . . . . . . . . . . . . . . . . . 845.37 Effects of higher order terms on energy recovery . . . . . . . . 855.38 Effects of higher order terms on energy spread . . . . . . . . . 865.39 Parameter evolution in the optimization population . . . . . . 876.1 Linac to undulator transport . . . . . . . . . . . . . . . . . . 926.2 Undulator to linac pass 2 transport . . . . . . . . . . . . . . . 926.3 EDBT transport . . . . . . . . . . . . . . . . . . . . . . . . . 936.4 EHAT transport . . . . . . . . . . . . . . . . . . . . . . . . . 936.5 Chicane in RLA mode . . . . . . . . . . . . . . . . . . . . . . 946.6 Removing chicane in RLA mode . . . . . . . . . . . . . . . . 95B.1 Optimization top level . . . . . . . . . . . . . . . . . . . . . . 114B.2 Layers of the optimization problem . . . . . . . . . . . . . . . 115ixList of FiguresB.3 Example topology . . . . . . . . . . . . . . . . . . . . . . . . 115B.4 Optimization platform class diagram . . . . . . . . . . . . . . 117B.5 Variator state chart . . . . . . . . . . . . . . . . . . . . . . . 119B.6 Selector state chart . . . . . . . . . . . . . . . . . . . . . . . . 120B.7 Evaluator execution flowchart . . . . . . . . . . . . . . . . . . 122B.8 Evaluator multithreading flowchart . . . . . . . . . . . . . . . 123B.9 Evaluator class diagrams . . . . . . . . . . . . . . . . . . . . . 124B.10 Population evaluation statechart . . . . . . . . . . . . . . . . 126B.11 Work distribution to computational clusters . . . . . . . . . . 127B.12 Process monitoring with monitoring thread . . . . . . . . . . 129B.13 Process monitoring without monitoring thread . . . . . . . . 130B.14 Minion flowchart . . . . . . . . . . . . . . . . . . . . . . . . . 131B.15 Optimization topology XML block example . . . . . . . . . . 136B.16 Optimization filesystem structure . . . . . . . . . . . . . . . . 140B.17 Optimization multithreading sequence diagram . . . . . . . . 142xList of Programs3.1 Python program for converting the Genesis output distribu-tion from binary to ASCII . . . . . . . . . . . . . . . . . . . . 223.2 Sample Empirical Model interpolation table for a solenoid . . 24xiAcknowledgementsI would like to thank the members of my doctorate committee, Yu-ChiuChao, TomMattison, Rick Baartman, Reiner Kruecken, and Mike Craddock,for providing valuable guidance during the course of my studies.This work benefited from insightful correspondences with Sven Reiche,Gabe Marcus, Michelle Shinn, and members of the TRIUMF Beam Dynam-ics Group and Accelerator Division.Thanks to the members of the UBC Physics Department and of TRIUMFfor providing a supportive environment for this work to take place.The formatting of this thesis is based on the UBC LATEX template byMichael Forbes.Lastly, a tremendous thank you to my supervisor Yu-Chiu Chao, whoseguidance was instrumental.xiiDedicationThis work is dedicated to my parents for their continued support, and toCatherine Kwan for her seemingly endless patience and empathy.xiiiChapter 1IntroductionThe goals of this research are:1. Design and implement a novel software optimization platform, capa-ble of multiobjective optimization in massively parallel computing sys-tems.2. Use the platform to study the underlying physics of the TRIUMFEnergy Recovery Linac (ERL).3. Use the platform to find a baseline solution to the ERL, including RFsettings and a complete coordinates layout of all optical components.The ERL is a future upgrade to the Electron linear accelerator (E-linac),a machine currently under development at TRIUMF. The questions of con-cern are, for a given ERL machine layout where attributes such as driftlengths and quadrupole strengths are variables, does a transport solutionexist which has good lasing, energy recovery, and beam disposal? How doessuch a machine work? What are the relationships between RF, gain, andbeam transport? This matter is complicated by the fact that currently, nosimulation software has the capacity to model all the relevant physics of anERL from start to end, therefore performing global optimization on such amachine is very difficult.The E-linac is the accelerator at the center of the Advanced Rare IsotopELaboratory (ARIEL), a TRIUMF facility designed to advance Canada’s ca-pabilities in science and technology. The E-linac is a 0.5 MW, 10 mA,continuous wave (CW) rare isotope beam (RIB) driver for the photofissionof actinide targets, with emphasis on neutron-rich isotopes. The linac sharescommonalities with Energy Recovery Linac (ERL) designs and can be up-graded to fulfill such a purpose. Further descriptions of ARIEL and ERLscan be found in chapter 2.The list of notations and conventions used in this thesis areshown in Appendix A.The optimization problem can be defined as follows: what is the underly-ing physics of the TRIUMF ERL? Given a minimalist layout of accelerator1Chapter 1. Introductionoptical elements, can a good transport solution be found (or exist), andcan this solution satisfy other functions of the ERL, namely, maximize las-ing and energy recovery, and coexist with existing RIB operations? Whatquadrupole strengths and positions, and RF parameters are required for theaforementioned solution?The study and design of accelerators is typically done with piecewiseoptimization. For example, find the optimal beam conditions to maximizelasing in the Free Electron Laser (FEL), then design the arcs to match thoseconditions, and then produce RF conditions necessary for the arc designs.This type of piecewise, or local optimization, does not explore the global so-lution space or show tradeoffs between important parameters. For example,is it acceptable if lasing is reduced slightly but has a large benefit for energyrecovery and beam disposal? This piecewise scheme arose from necessityand practicality. No accelerator modeling tool currently has the capabilitiesto model every physical process contained in an ERL. The software for mod-eling the FEL cannot model the transport arcs or accelerating cavities, andvice versa. At the same time, attempting to perform global optimizationby exploring a “grid” of different local optimizations quickly fails for largeproblems due to the large number of permutations needed. If only two freeparameters exist, we can scan them, but if a large number of parametersexists, scanning is impractical. The ERL optimization contains over 60 freeparameters. An exhaustive scan is unfeasible.We solve this problem by creating a software platform for optimizationwhich can model the global ERL, and can make use of the many modelingtools available to model any physics necessary. The software architecturaldesign of the optimization platform is provided in Appendix B. A summaryof the platform’s major functionalities:1. Able to perform global ERL optimization. We are primarily interestedin modeling the ERL from acceleration, to lasing, to deceleration, tobeam disposal. Although designed with the ERL in mind, the plat-form is generic and supports many optimization problems from outsideaccelerator physics. It has been used in several other optimizationproblems, including the E-linac injector [49, 50], CSR [47], and theVECC injector linac [37]. Global optimization also allows us to seetradeoffs and dynamic relationships of the machine.2. Interface with different modeling tools, or engines. The platform cancall MADX [27] for arc design, GENESIS [84] for lasing, and automat-ically transfer values from one to the other for continuous modeling,2Chapter 1. Introductione.g. use the MADX output beam parameters as the input beam pa-rameters for GENESIS. This mitigates the need for local optimizationsor permutations. The engines can be combined in serial, parallel, orany combination thereof. This defines an optimization problem topol-ogy, with each modeled section being a vertex. The ERL topology canbe seen in 4.2.3. New modeling engines are easily added to the framework if necessary.This increases the flexibility of the platform and its adaptability tonew problems.4. Parallel capable - reduces running time by taking advantage of com-putational clusters. The platform is designed to take advantage ofCanada’s WestGrid [8] and similar batch systems.5. Exception handling - by allowing different engines, a multitude of er-rors can occur, including engine hang-ups that can destroy a time-consuming run. Dealing with large parallel jobs can also produce net-working and file system issues. The optimization platform is smartenough to handle the multitude of exceptions that can occur. Thisrequired significant investments in good software design.6. The platform uses a flexible XML input scheme to define the list ofobjectives, constraints, and parameters, as well as the optimizationtopology. The scheme is generic and not tailored toward any onemodeling engine.7. Allows the user to inject code for data manipulation and post process-ing.Designing the platform to operate with arbitrary engines with multithread-ing capabilities was a significant software challenge. Multithreading andparallel processing especially increased the complexity of the software de-sign due to concurrency issues and threading exceptions.The optimization platform uses the genetic algorithm, a class of mul-tiobjective optimization algorithms. The advantages of genetic algorithmsare:ˆ Multiobjective - the ERL problem involves several very different ob-jectives: maximize lasing, optimal energy recovery, beam disposal,and beam transport. A multiobjective algorithm allows us to see thetradeoffs, or Pareto fronts, between competing objectives. Optimizing3Chapter 1. Introductionusing a single objective algorithm requires the different objectives tobe combined into a single objective, in an arbitrary way.ˆ Can find multiple solutions - this is critical in a multiobjective problem,as one solution may not be optimal in all dimensions. An ERL thatprovides the best gain may not have the best energy recovery. Multiplesolutions can form Pareto fronts, which allow us to compare differentoptions and see the tradeoffs.ˆ Global - genetic algorithms do not require specific search space ge-ometries or convexity. Note: the platform is “global” in two senses. Itallows the global modeling of the ERL, and the optimization algorithmcan search the global search space.ˆ No gradients required.ˆ No initial search point required.Mathematical definitions of optimization and Pareto fronts can be foundin the next section. An example population with the Pareto front defined isshown in Fig. 1.1.The particular variant of the genetic algorithm used is the StrengthPareto Evolutionary Algorithm (SPEA2) [98]. The algorithm works as fol-lows:1. Create a set of ERL designs. Each ERL design has a random combi-nation of parameters, such as quad gradients and RF phases. This isthe initial population.2. Each design in the population, or individual, is evaluated on how wellthey satisfy the design objectives and constraints (fitness). The worstdesigns are thrown out.3. New designs are created by tweaking the best of the old designs (evo-lution), and added to the population to replace the ones thrown away.This is the new generation.4. Go to step 2, until a predetermined number of generations have passed.SPEA2 performed well in a comparison of genetic algorithms and justi-fies this choice [99]. Genetic algorithms are also comparable in efficacy andrunning time to Particle Swarm Optimization, another popular class of pop-ulation based search algorithms, in both single and multiobjective problems4Multiobjective Optimization[54, 58, 68, 95]. More information on SPEA2 and genetic algorithms can befound in the next section.The implementation of the TRIUMF optimization platform follows thePlatform and Programming Language Independent Interface for Search Al-gorithms (PISA [19]) concept which, for flexibility purposes, divides theplatform into two portions, Variator, responsible for evolution, and Eval-uator, responsible for fitness testing, i.e. how well does a solution satisfythe objectives. Other precursors to the TRIUMF optimization platform areAlternate PISA (APISA [11, 13]) written by I. Bazarov and Yet AnotherPISA (YAPISA [70]) by G. Goh. The TRIUMF platform is the first toallow for multiple simulation engines in the same optimization setup, andextending the software to allow for parallel computing was a challenge andstep forward.More details of the software aspect of the optimization platform are inAppendix B. For the application of the platform to the TRIUMF ERL seechapters 3 and 4. ERL optimization results can be found in chapters 5 and6.Multiobjective OptimizationWe now provide the mathematical formulation of a multiobjective optimiza-tion problem. A generic optimization problem can be stated as, given someoptimization parameters (or free parameters)x = (xi)s.t. Li ≤ xi ≤ Ui(1.1)where constants Li and Ui are the lower and upper bounds on xi and i =1, ..., N , the goal is to optimize one or more objective functionsmin gj(x) (1.2)where j = 1, ..., P . At least one objective is required. The system is subjectto multiple constraintsFk(x)≤Ck (1.3)where k = 1, ..., Q and Ck are constants. Optimizing with zero constraintsis allowed.The space of all x which satisfies the constraints is the search space. Wedo not assume the search space is convex (for definition of convexity, see[18]), hence our choice of the genetic algorithm, a pseudo-random search5Multiobjective Optimizationmethod. In theory, if given enough time, this algorithm will always find theoptimum solutions. Note that the algorithm SPEA2 can generate individualsthat do not satisfy all constraints, but these points are heavily biased againstin fitness selection and should disappear in later iterations.A point S1 in search space is dominant over point S2 if for all objectivefunctions gj∀gj(S1) ≤ gj(S2) ∧ ∃gj(S1) < gj(S2) (1.4)and S1 is nondominant over S2 if∃gj(S1) > gj(S2). (1.5)Solutions to the problem 1.2 is defined as the Pareto front. A point isin the Pareto front if it is not dominated by any other point. The goal ofoptimization is to find the Pareto front.A B C D E F G Emittance x Cost c Figure 1.1: Example population in objective space. Each blue point is partof the Pareto front because they are not dominated by any other points.An example of a Pareto front is shown in Fig. 1.1. Given a problem withtwo objectives: 1) minimize the machine cost c and 2) minimize the beamemittance ε, the optimization platform generates a population of differentmachines. The genetic algorithm ranks the fitness of each machine based onPareto dominance (Eqn. 1.4), shown in Table 1.1. We define the fitness ofthe machine by how many other machines dominate it (the lower the fitnessscore, the more fit it is).6Multiobjective OptimizationTable 1.1: Pareto dominance of points shown in Fig. 1.1. The definitions ofdominance and nondominance are given by eqns. 1.4 and 1.5. The fitnessof the machine is equal to how many machines dominate it. The lower thefitness score, the more fit that machine is.Machine Dominated by Fitness scoreA - 0B - 0C B,D 2D - 0E D 1F - 0G A,B,C,D,E,F 6Points not dominated by any other points form the Pareto front, e.g. A-B-D-F. Each iteration, the algorithm generates new machines which canachieve a better front. An example is shown in Fig. 1.2. A new machine,H, is created by the optimization platform, forming a new front A-H-F.A B C D E F G Emittance x Cost c H Figure 1.2: Example population in objective space. New machines are gen-erated which can change the Pareto front.The dominance shown in Table 1.1 is also used by the genetic algorithmas weights. The most dominated machine, G, is the most likely to be tossed7Multiobjective Optimizationaway, whereas machines on the Pareto front are the most likely to be keptand used as parents for new machines.8Chapter 2Description of the TRIUMFE-Linac and ERLEnergy Recovery Linacs (ERL) are a class of particle accelerators with a longhistory [51] dating back to 1965 when M. Tigner first proposed the concept[93]. Tigner discussed the large power required by radio-frequency (RF)accelerators to produce beams of sufficient quality for collider experiments.A modest collider with 500 MeV energy for each of two beams at 100 mAcurrent requires more than 100 MW RF power, with 50 MW carried awayby each beam. The high power requirement, combined with the fact that fora collider, the majority of the beam do not interact (e.g. the Large HadronCollider has 20 collisions per two crossing 108 proton bunches [3]), meansthat most of the beam power is wasted. While such high RF power maybe achievable (even by 1960 standards), the bluntness and inefficiency (notto mention economic consequences) of the approach encouraged Tigner toenvision alternative machine designs with energy recovery capabilities. Oneof his designs uses bending magnets to turn the beam back into the same RFdevice used for acceleration where the beam energy is given back to the RFsystem. This is the underlying principle of ERLs today. Energy recoveryallows for two benefits: 1) capital savings resulting from using less RF powerand/or 2) higher powered operation due to less demand on the RF systemcompared to a equivalent single-pass machine.In the following decades several ERL experiments and proofs-of-concepttook place. The first successful demonstration of energy recovery was per-formed in 1975 at the Chalk River Nuclear Laboratories in Canada [87].An electron beam was passed through the same S-band (microwave) linearaccelerator twice. The beam was accelerated to 8 MeV after the first pass.After the second pass the beam energy can range from 3 to 16 MeV, de-pending on the separation between the accelerator and the magnetic system.The beam passes through the linac once, is reflected by a compact magneticsystem, and goes through the linac a second time in the opposite direction.For this reason the machine is named the Reflexotron.The invention of the free electron laser (FEL) by J. Madey in 1971 [69]9Chapter 2. Description of the TRIUMF E-Linac and ERLopened up interesting avenues of research for accelerator driven FELs. Theconcept of energy recovery is very attractive for such machines, since FELlasing only use ≈ 1% of the electron beam energy.In 1985 an ERL at the University of California at Santa Barbara becameoperational [42]. Unlike Tigner’s vision of ERL for collider experiments, the3 MeV machine was used as a light source, using the ERL to drive anFEL operating in the far-infrared range (initial lasing at 400 µm). Themachine demonstrated the versatility of ERL based sources in condensedmatter experiments [43]. The UCSB machine uses electrostatic rather thanRF based accelerating devices, but does incorporate energy recovery.FEL experiments took place on the Stanford superconducting FEL (SCA-FEL) in the late 1980s [90]. The accelerator uses a 50 MeV RF linac sur-rounded by a recirculation loop. The loop had a controllable path length,allowing the machine to be operated in two modes: energy doubling mode(two acceleration passes in the linac) and energy recovery mode (one accel-eration pass and one deceleration pass). The purpose of the energy doublingmode was to decrease the FEL photon wavelength, which is proportional tothe inverse of the square of the electron Lorentz factor. Photon wavelengthsdown to 0.5 µm were achieved [85]. The SCA-FEL successfully demonstratedenergy recovery, although the demonstration was limited and not performedsimultaneously with FEL lasing.Near the same time, another successful demonstration of energy recov-ery took place in the Los Alamos National Laboratory FEL [44]. Electronsare accelerated by a 20 MeV accelerator (divided into two sections), passesthrough a wiggler, then through a 20 MeV decelerator. The deceleratorand the accelerator are separate devices, but are connected to the same RFgenerators through resonant bridge couplers (RF power splitter). Energy re-covery down to 3.5 MeV was achieved. The experiment was performed withthe wiggler at 0.7% extraction efficiency, 100 mA average electron current,and 68% recovery.The first high powered FEL was the IR Demo at Jefferson Lab [77].The IR Demo first lased in 1999 with 1.72 kW average power at 3.1 µm[76], compared to the previous record of 11 W at Vanderbilt University in1990 [22]. The machine convincingly demonstrated energy recovery usingthe same cavities for acceleration and deceleration, referred to as same-cellenergy recovery (SCER), while lasing simultaneously. This is in comparisonto earlier experiments which energy recovered without lasing [90] or did notuse the same cavities for acceleration and deceleration [44]. IncorporatingSCER resulted in a 5× reduction in RF power, allowing up to 5 mA av-erage beam current at ≈ 48 MeV beam energy, compared to 1.1 mA when10The TRIUMF ERLthe machine was operated in single pass mode (without SCER). A secondERL-FEL was constructed at JLab, named the IR Upgrade. The machineachieved lasing in 2003 at 5.7 µm with 8.5 kW average power [14, 16]. Thiswas subsequently optimized to 10.5 kW by 2006 [17].An ERL based FEL was developed at the Japan Atomic Energy Agency(formerly Japan Atomic Energy Research Institute). The machine was ini-tially a linac-driven FEL without ERL capabilities. It uses a 15 MeV su-perconducting RF linac to produce infrared lasing at 24 µm with 0.1 kWaverage power [74]. The linac-FEL was designed and upgraded to an ERL-FEL with both energy recovery and lasing demonstrated in 2002 [52, 53].The upgraded machine achieved 98% energy recovery.Other energy recovery experiments and machines include a GeV scaleenergy recovery experiment at the Continuous Electron Beam AcceleratorFacility [20] and an FEL driven by a multipass racetrack microtron in Novosi-birsk [9]. These ERL projects highlight the efficiency of having energy re-covery capabilities and the versatility of an FEL light source.The TRIUMF ERLThe TRIUMF E-linac, currently under construction, is a driver for pro-ducing rare isotope beams (RIB) via photofission of actinide targets [63].The emphasis is on producing neutron-rich isotopes for studies on nuclearstructure physics and material properties.The intention is to upgrade the E-linac to an energy recovery linac (ERL)and use the E-linac as the driver of an infrared light source. This requiresthe addition of a recirculation lattice, a free electron laser (FEL), and a newERL gun to the main linac. The operation is as follows (refer to Fig. 2.1):1. Beam from the ERL injector is sent into the main linac and acceleratedfrom 7.5 MeV to 45 MeV.2. ERL beam transported to the recirculating loop via RF separator.Three beams are separated by the separator: the ERL beam, therare isotope beam, and the ERL dump beam. RIB is separated fromthe ERL beam by matching to different RF phases in the separator,which introduces different dipole kicks to transport the beam to theirrespective beamline. The accelerated ERL beam and the dump ERLbeam are separated via dispersive effects in the first separator dipole.3. The loop is divided into two arcs, each introducing a 180◦ bend to the11The TRIUMF ERLbeam. The first arc turns the beam antiparallel to the main linac. Thefirst arc is followed by a four-dipole chicane for bunch compression.4. The beam travels through the Free Electron Laser (FEL) undulator,where electron energy is converted to coherent radiation via electron-light interaction (Fig. 2.2).5. After lasing, the beam goes through the second arc where it is broughtback to the main linac out-of-phase for deceleration back to 7.5 MeV.Beam energy is returned to the RF system.6. The decelerated beam finally is dumped after the second linac pass.The TRIUMF ERL is designed to be a dual RIB and ERL machine.The two beams are produced in two separate guns but share the injectorand linac transport. RIB operates at 650 MHz and fills every second bucketof the 1300 MHz RF system. The ERL beam will occupy unused RF buckets(RF periods unoccupied by RIB), therefore does not lower RIB performanceand allows for simultaneous operation.A complete overview of ARIEL and the E-linac can be found at [63].linac merger separator undulator chicane ERL/RIB shares common injector ERL pass 1 RIB pass 1 RIB extraction ERL pass 2 ERL recirculation ERL beam compression lasing ERL recirculation Figure 2.1: Overview of the TRIUMF dual ERL/RIB.12Principle of Free Electron LaserPrinciple of Free Electron LaserA Free Electron Laser (FEL) is a high powered light source, opening thedoor to many diffraction and microscopy experiments [15, 55, 63].The FEL converts electron beam energy into coherent radiation. Thedevice consists of an optical cavity formed from two mirrors surroundingan undulator, a series of alternating dipole magnets. Electrons travel ina sinusoidal trajectory, producing initial light via synchrotron radiation.The light is trapped in the optical cavity. Past this initial seeding, furtherradiation is produced via electron-light interactions [59].The undulator dipole magnets are located in close proximity to eachother, and are constructed from permanent magnets [36]. Electromagnetsrequire greater space between each other for coil placement and are not used.The E-linac undulator is designed as a planar undulator.The fundamental laser wavelength is given by the undulator equation [24,59]λ =λu2γ2r(1 +K2)(2.1)Where λu is the spacing or period of the undulator magnets, γ is theLorentz factor of the bunch centroid, and K is the dimensionless undulatorparameter, dependent on the dipole magnetic field and spacing. K havetypical values on the order of unity [66]. Taking λu = 4 cm and K = 0.7,estimated from machines of similar energy [34], we find λ to be severalmicrons, in the short- to mid-wavelength infrared region. It is theoreticallyeasy to operate the FEL up to several hundred microns in the far-infraredregion, by tuning the electron beam energy or K, which can be tuned bychanging the magnetic gap size [78].13Principle of Free Electron LaserFigure 2.2: An FEL-oscillator converts beam energy into laser energy. Thebeam travels in a sinusoidal path in the undulator, a series of alternatingdipole magnets. The energy is captured in a cavity by two high-reflectivitymirrors and increases until saturation. A small amount of power is out-coupled for useful purposes. The beam after lasing, or exhaust beam, isdistorted by the lasing process and is no longer useful. It is transportedback to the linac for energy recovery.These parameters are suitable for an FEL of the low gain oscillator type,i.e. using synchrotron radiation to induce lasing and having a mirror systemto capture the light while allowing a small amount through. Mirrors ofhigh-reflectivity are relatively easy to design for infrared light. Using anoptical cavity to store energy negates the need for a very long undulator,which is required for SASE [24], an alternative FEL with a one-pass opticalsystem. A short undulator also causes less beam distortions, making thebeam suitable for energy recovery.Physical RequirementsAt the writing of this document, the scientific requirements are not well-defined; therefore parameters of both the FEL and ERL can only be ap-proximated. The undulator period and number of periods were chosen tobe 4 cm and 25, respectively, for a total undulator length of 1 m. Thesevalues were based off the Peking University FEL design [34], a machine withsimilar beam energy and parameters to the E-linac.14Principle of Energy RecoveryPrinciple of Energy RecoveryThe recirculation lattice brings the ERL beam back into the linac for a sec-ond pass. The path length of the lattice is adjusted such that the beamarrives for the second pass at a decelerating RF phase, i.e. 180◦ offset fromthe accelerating pass, where the beam energy is reduced to the initial in-jection energy. The advantages of energy recovery are 1) the deceleratedbeam is easy to dispose of and 2) the beam energy is transferred back to thelinac, which intensifies the RF field for acceleration of further bunches. Thenet beam loading in an ERL is close to zero, thus consuming only a modestamount of RF power [73]. An illustration of energy recovery is shown inFig. 2.3.Linac RF waveform at 1.3 GHz 180o for optimal energy recovery Beam on  linac pass 1 Beam on linac pass 2 Figure 2.3: Energy recovery requires timing the pass 1 (accelerating) andpass 2 (decelerating) RF phases. When the two phases are offset by exactly180◦, optimal energy recovery occurs and the linac operates as if no beamloading is present. Equivalently stated, no RF power is drawn for the beam.Physical RequirementsThe ERL adds a recirculation lattice and a second gun to the E-linac. Thishigh brightness gun shares the injection cryomodule and main linac with thephotofission beam. The gun is projected to operate at 130 MHz, as opposedto the 650 MHz RIB, for a RIB to ERL beam ratio of 5:1. The ERL beamoccupies empty RF buckets, so does not interfere with RIB operation orrequire RIB to operate sub-optimally. The bunch charge is 100 pC, resultingin 10 mA beam current and 0.5 MW power.The injector transport, merger, main linac, and separator systems areshared by RIB and ERL operations. The separation system consists of an RFseparator, followed by a drift and ending with a septum. The recirculation15Design Tradeoffslattice begins at the septum exit, turns the beam 180◦, goes through theFEL, followed by another 180◦ arc, for a full 360◦ turn back into the mergerand the main linac for deceleration. The merger system and the main linacaccommodate three beams: ERL injected beam, ERL recirculated beam,and RIB injected beam. Refer to Fig. 4.1 for the machine layout schematic.Beam loading vector prefers a deceleration phase exactly 180◦ offset fromthe acceleration phase [46, 63]. With no beam loading, the optimal Q-valueQL,opt for an ERL is QL,opt = f/(2∆f) [65], where f is the cavity fundamen-tal mode and ∆f is the cavity detuning driven primarily by microphonics.When the ERL is operated in 180◦ mode, the machine operates as a highQ machine with QL = QL,opt [82]. The two beam current vectors canceland the RF system resembles one with no beam loading. The JLab IRFELenergy recovery demonstration shows that RF power was reduced from 36kW with no energy recovery, to 16 kW with energy recovery, at 1.1 mAbeam current [72].Design TradeoffsThe TRIUMF ERL contains many design challenges. Some items of interestare, but not limited to,1. How to design the bunch compression system? How do the RF and arctransport parameters play off against each and the undulator gain?2. How does the arc time-of-flight (TOF) impact the RF decelerationphase and energy recovery?3. How does the RF deceleration phase impact energy spread?4. Can lasing cause issues leading to beam loss?Results from optimization detailing tradeoff issues can be found in chapter5.16Chapter 3Description of the GlobalOptimization Platform andEngines Used for ERLModelingThe software architecture of the platform follows the PISA interface [19].The platform is divided into two parts: Variator and Selector.Variator handles the creation and evaluation of new ERL designs in thepopulation. The evaluation of each ERL design is broken into sections,with each section being handled by a modeling engine. Each section is as-signed to a worker node in a parallel environment for execution. Variatorhandles all parallel computing work assignments and ensures work is per-formed smoothly through a series of process control and exception handlingmechanisms.Selector handles the optimization using the genetic algorithm SPEA2[98]. Given an optimization population, Selector tests for Pareto dominanceand stochastically chooses ERL designs to use as parents for the next it-eration. Steps taken by the algorithm were outlined previously in chapter1.A complete description of the design and implementation of the opti-mization platform is provided in Appendix B.Connecting Different Engines for Global ModelingEach engine used in the ERL optimization (detailed below) is treated bythe optimization platform as independent. They are modular in that noengine, or machine section, need to explicitly know about any other. TheC++ code of the platform automatically creates the necessary input filesfor these engines to run and read the outputs from them.Transitions between engines are handled by the optimization platform.17Engines Used for ModelingFor instance, MADX and the Empirical Model (EM) are often used in serial.EM is used to model a cavity, followed by a transport section modeledby MADX. The transitioning is quite involved. For example, units needto be converted (EM uses MeV, MADX uses GeV). EM uses the beamsize, which needs to be converted to Courant-Snyder parameters for MADX(see Appendix A for definition). This complicated process is hidden fromthe user. The platform automatically performs unit conversions, parametertranslations, and other niceties.For a detailed explanation of how the engines are incorporated into theplatform and how to define the engine landscape for an optimization prob-lem, see Appendix B.Engines Used for ModelingAn engine in the context of optimization refers to a modeling software. Dif-ferent engines are used to model different sections of the machine, allowingus to take advantage of the features of each engine. The different enginesused for optimization are listed in this section.MADXMADX [27] is an accelerator design code from CERN. The MADX Twissmodule is a convenient method to retrieve beam parameters and transportmap elements of the lattice. However, MADX does not have good capabili-ties in tracking longitudinal parameters and Twiss calculation is 4D trans-verse only. Sections that involve non-RF elements are modeled by DIMADand MADX in parallel, which mitigates each engine’s disadvantages. Theparallel MADX-DIMAD modeling can be used to cross-check the validity ofeach other’s output.DIMADDIMAD v2.9 [88] is a tracking code with a similar input format to MADX.DIMAD has many modules, providing easy and flexible access to beam infor-mation. (for example, beam size is readily available in the BEAM module,including contributions from betatron oscillations, dispersion, and higherorder terms). Both envelope and particle tracking are used. DIMAD’s MA-TRIX module is useful for displaying second order transfer map elementsfor an entire section, whereas MADX requires manual concatenation of mapelements for individual optical devices.18Engines Used for ModelingA drawback of DIMAD is that the floating point output behaviour is of-ten unpredictable. Output of the MACHINE command may replace floatingpoint values with the string ********** if the floating point requires morethan ten significant digits to display. The output can also produce run-on floating point values, i.e. no delimiters, making automatic text parsingdifficult or impossible.GenesisGenesis v1.3 [84] is an industry standard for simulating undulator radiation.In the ERL optimization Genesis is used to produce the gain of the FEL atthe resonance wavelength. Genesis requires an input wavelength λ. Beamparameters can shift the resonance, therefore the gain is not optimal if λ isoff-resonance. We use Genesis’ scan feature to look for λ which maximizesgain. The scan feature is compatible only with running Genesis in steady-state mode. The electron bunch length is ∼ 100 times longer than theradiation wavelength (σz = 100 µm vs λ = 1 µm), therefore the steady-statemode is justified.The beam envelope is calculated from the distribution from Genesis.The beam coordinates x, x′, y, y′, δ are translated directly from the Gen-esis distribution. The coordinate z is not easily obtained due the difficultyin translating from the ponderomotive phase θ (see Appendix A for defini-tion). Genesis models the bunch as a series of slices, each with a radiationwavelength thickness. Transition between slices is not modeled [71] andthis complicates the translation from θ to z. We choose to keep z constantthroughout the undulator and the justification is as follows.The bunch length after the undulator is assumed to be identical to theinitial bunch length. This can be justified by looking at the longitudinalphase space evolution. For a single particle in the distribution, expand tofirst order the phase θ and ξ ≡ γ/γr − 1, the energy deviation from thebunch mean energy γr [59]:θ(r) = θ0 + ǫθ1= θ0 +ǫξ0(sin θ0 − sinφ02kuξ0− r cosφ0)θ0(r) = 2kuξ0r + φ0(3.1)where r is the position of the bunch center within the undulator and φ0 is19Engines Used for Modelingthe initial phase. The expansion parameter is given byǫ =eE0K[JJ ]2γ2rmc2[JJ ] = J0(K24 + 2K2)− J1(K24 + 2K2) (3.2)where the undulator parameter K = 0.7 and J0 and J1 are Bessel functions.The radiation field E0 can be estimated from the energy density S and beamcross sectional area AE0 =√Scµ0 =√Pcµ0A=√Pcµ0πσxσy= 60 MV/m (3.3)With P = 8 MW, σx = σy = 0.5 mm, this leads to ǫ = 5×10−3 m−1. Fora large energy spread of 0.01, the change in phase ∆θ = θ(r = 1 m)− θ(r =0 m) is 2.6 at φ0 = 0, 2.8 at φ0 = π/2, 3.6 at φ0 = π, and 2.6 at φ0 = 2π.Converting θ to z using the conversion factor of λ/2π, the change in bunchlength is several radiation wavelengths, even less at the beam interior. Giventhat the bunch occupies several 100 radiation wavelengths, particles withinthe distribution experience very little longitudinal movement. Thereforechanges in z are considered negligible. Under this assumption, the envelopeparameters evolve in the undulator asσz,+ = σzσδ,+ = BσδVzδ,+ = BVzδεz,+ = Bεz(3.4)where the (+) denotes the end of the undulator, B is a scaling constantdetermined from the Genesis output distribution, σz is the bunch length, σδis the bunch energy spread, Vzδ is the covariance between the coordinates zand δ, and ε is the emittance.A distribution is created from the Genesis output and tracked in arc 2.Momentum tail is a concern and can result in particle loss in the arc dipoles.We create a distribution using the momentum profile of one Genesis slice.Fig. 3.1 shows that the momentum profile varies little across simulationslices, justifying this approach.20Engines Used for Modeling160 200 230 260 300slices [arbitrary units]0102030I[A]−15−10 −5 0 5 10 15θ[λ/2π]8687888990γ−15−10 −5 0 5 10 15θ[λ/2π]8687888990γ−15−10 −5 0 5 10 15θ[λ/2π]8687888990γ−15−10 −5 0 5 10 15θ[λ/2π]8687888990γ−15−10 −5 0 5 10 15θ[λ/2π]8687888990γFigure 3.1: Top left: Genesis simulation slices vs the current in each slice, intime-dependent mode, at the undulator exit. The five following plots (leftto right, top to bottom) are the longitudinal profile at the five vertical redlines (slices 160, 200, 230, 260, 300). Plotted are the particle Lorentz factorγ vs the ponderomotive phase θ. The phase spaces undergo near identicaldistortions in the undulator, thus justifying using the steady-state mode inGenesis.Genesis outputs a particle distribution in binary format. In order totrack the distribution through the return arc, the binary was converted intoASCII via the Python script 3.1.21Engines Used for ModelingProgram 3.1 Python program that converts the Genesis output distribu-tion from binary to ASCII.data = [] # each row of ’data’ is one slicebuf = NCOLS*NPARTS # number of values for each slice# Read all data for a slice at a time,# append each to the array.# Each value is 8 bytes.# The final array has dimensions [NSLICES,buf]with open(filename,mode=’r’) as fp:for d in iter(lambda: fp.read(8*buf),’’):c = struct.unpack(’d’*buf,d)data.append(c)# transform to numpy matrixdata = numpy.array(data)NSLICES = len(data)slicedata = data[slice_to_read-1,:]# NOTE: np.reshape differs from Matlab reshape.# Suppose x=1...10:# np.reshape(x,[2,5]) returns# [[1 2 3 4 5]# [6 7 8 9 10]]# Matlab reshape(x,2,5) returns# 1 3 5 7 9# 2 4 6 8 10slicedata2= numpy.reshape(slicedata,[NCOLS,NPARTS])gamma = slicedata2[0,:]phase = slicedata2[1,:] # ponderomotive phasex = slicedata2[2,:]y = slicedata2[3,:]px = slicedata2[4,:] # gamma*betapy = slicedata2[5,:] # gamma*beta22Engines Used for ModelingEmpirical ModelAt TRIUMF, tracking the electron beam through RF devices is typicallydone using ASTRA [45]. The Empirical Model (EM) was written (by thepresent author) as a replacement to ASTRA, and works by using transfermaps pre-generated from ASTRA data [28, 32, 35, 48]. EM uses particletracking; each particle is propagated using map elements interpolated fromparticle momentum, RF phase, and RF amplitude.Astra EM Figure 3.2: ASTRA is the accelerator physics standard for tracking in RFdevices. The Empirical Model as an alternative due to the large compu-tation requirements of optimization. The EM running time is an order ofmagnitude less than ASTRA due to the small number of interpolation slicescompared to the large number of steps required for ASTRA tracking.Fig. 3.2 illustrates how EM tracking is performed compared to ASTRA.ASTRA, when given electromagnetic field maps of beamline devices, per-forms Runge-Kutta type tracking. This requires many small steps througha device. EM breaks the device into slices with the physics of each slice rep-resented by a transfer map. The typical number of slices is on the order ofone, much less than the number of Runge-Kutta steps necessary and resultsin a large saving in computation cost. For the ERL cavity, the EM modelrunning time is one order of magnitude less than the equivalent ASTRAmodel.The recirculation time-of-flight (TOF) is important for ERL phase match-ing between the two linac passes. EM tracks using time instead of position,thus the TOF through the cavity is easily obtainable. For devices withenergy changes, such as cavities, calculating the TOF is tricky because itcannot be easily inferred from the path length and velocity, thus makingEM an important component in the model.EM was created with the intention of having ASTRA’s accuracy, butwith running time suitable for online tracking in TRIUMF high level appli-cations using Java [6]. The low running time also makes EM suitable for23Engines Used for Modelingglobal optimization and a C++ port was made with such a purpose. Fig.3.3 shows the validation of EM tracking results compared against ASTRA.A sample of an interpolation table is shown in Program 3.2. The tableentries are coefficients of a Taylor expansion. For each column heading, thesix digits after the ‘x’ are the exponents of the dependencies on the six beamcoordinates. The column 1x102000 for instance, is the map element thatrepresents the dependence of x on x0 and y20. It is a third order transportelement, with one derivative against x0 and two against y0. Each line in thetable represents one grid point.Program 3.2 Sample Empirical Model interpolation table for a solenoid.The solenoid transfer map is interpolated from two coordinates: the fieldamplitude and the particle momentum. Only a portion of the full table isshown.# Empirical Element DataMaxB(G) P(MeV) 3x000000 6x000000 1x1020001.6000E-02 5.0380E-01 7.9411E-01 5.1924E-01 -2.6807E-011.6000E-02 5.3529E-01 8.1728E-01 4.9543E-01 -2.5850E-011.6000E-02 5.6678E-01 8.3556E-01 4.7528E-01 -2.4887E-011.6000E-02 5.9827E-01 8.5181E-01 4.5536E-01 -2.3963E-011.6000E-02 6.2976E-01 8.6582E-01 4.3676E-01 -2.3072E-011.6000E-02 6.6124E-01 8.7789E-01 4.1937E-01 -2.2245E-01The interpolation algorithm assumes the data grid is non-uniform. Thisis useful for elements that require different fineness for different parameterregions. Interpolation is performed using a combination of inverse-distanceweighting and nearest neighbor, and in the case of a regular grid, shouldsimplify to linear interpolation. Program 3.2 shows interpolation against twocoordinates (momentum and B-field amplitude). However, the implementedalgorithm works for arbitrary interpolation coordinate dimensions. Linearinterpolation is performed. See the EM design document [48] for algorithmicdetails.24Engines Used for Modeling−6 −4 −2 0 2 4 6 8x[mm]−20−15−10−505101520px[keV/c]AstraEM−3 −2 −1 0 1 2 3y[mm]−30−20−100102030py[keV/c]−1.0 −0.5 0.0 0.5 1.0z[mm]−0.4−[MeV/c]−6 −4 −2 0 2 4 6x[mm]−15−10−5051015px[keV/c]−6 −4 −2 0 2 4 6y[mm]−40−2002040py[keV/c]−1.5−1.0−0.5 0.0 0.5 1.0z[mm]−0.4−0.3−0.2−[MeV/c]Figure 3.3: Tracking through a 1.3 m 9-cell cavity, followed by a 0.65 mdrift, and another 1.3 m cavity (blue is Astra, red is the Empirical Model).The left column represents the initial beam phase spaces. The right are thephase spaces after tracking.25Chapter 4Modeling the ERLThis chapter details the modeling of the ERL and decisions made in theoptimization setup.The machine is modeled starting from the main linac entrance. The in-jector is not modeled in the optimization due to 1) RIB injection transport iswell-defined at the writing of this document, with the gun and injection cry-omodule built, thus unlikely to be affected by the results of optimization, 2)simulating the already-built RIB injector only increases simulation runningtime, and 3) a preliminary ERL transport tune exists for the linac, thereforewe can take the beam conditions at the linac entrance as the ERL input.Modeling software was chosen and/or written to provide full 6D transportdescription. Another important factor is for the software to be suitably ac-curate in the TRIUMF energy regime (10 to 50 MeV). The modeled machineis shown in Fig. 4.1. The beam, for most of the modeling, is represented byfirst order envelope parameters, which proved to be sufficient.The ERL beam is modeled with the deceleration pass and transport tothe exhaust beam dump. RIB modeling tracks the high energy beam pastthe separator septum and two more quads. The rest of the RIB delivery lineis not modeled.A 3 mm maximum transverse RMS beam size is strictly enforced at allpoints of transport. These are the acceptable conditions to minimize beamloss in the existing E-linac and high energy beam delivery design and wechoose to impose them globally in the recirculation loop as well.The main linac layout design is concrete and not debatable. However,the cavity gradients, cavity phases, and strengths of the quadrupole tripletin the center of the linac are free parameters for optimization. Likewise,positions of the merger and separator elements adhere to existing designnotes, but strengths of quadrupole magnets in these areas can be optimized.The rest of the recirculation lattice, from exit of the separator septum tothe merger entrance, is not tied to existing designs, thus all elements in thissection are free parameters.26Objectives and ConstraintsObjectives and ConstraintsHere we define the objectives and constraints used in the global optimizationproblem. Top level criteria are used when possible.The optimization objectives are (see Appendix A for notations):1. Find viable electron transport for the ERL. Solution must also supportRIB transport.2. Maximize gain: maximum FEL lasing, therefore radiation power.3. Edmp = Ein = 7.5 MeV: energy recovery condition. If dump energyequals injection energy, the beam loading vectors cancel for maximumRF efficiency [64].4. σx,dmp = 7 mm, σy,dmp = 7 mm: at the ERL beam dump, we relax thebeam size condition. Instead, we want to blow up the beam to reduceradiation heating.5. ηx = ηx′ = 0 at end of both arcs: dispersion suppression for botharcs. This separates the lattice into dispersion independent sectionsfor better modularity and tuning. Note that here arc 1 includes theseparator and mirror separator.6. αx = αy = ηx′ = 0 at both arc centers: we try to look for a designwith arc symmetries in βx, βy, and ηx. Although not critical to ERLoperations, symmetries make tuning easier. The layout of the arcoptics is symmetric as well to accommodate these conditions. ηy iszero everywhere and does not need to be considered. There is overlapwith item 5.7. Maximize ERIB : ERL operations cannot interfere with RIB opera-tions. High RIB energy is desired for photofission.The design constraints are:1. σx ≤ 3 mm, σy ≤ 3 mm: minimize beam loss by restricting transversebeam size everywhere. The bunch is modeled as Gaussian thus re-stricting the sigmas is sufficient. In the case when the bunch deviatesfrom Gaussian see item 3.2. σx,EDBT ≤ 3 mm, σy,EDBT ≤ 3 mm: for beam disposal, the beam sizein the dump transport section EDBT must be constrained (exceptin the last drift and the beam dump). This is complicated by the27Design Considerationslarge energy spread obtained from lasing and amplified by deceleration,which is converted to beam size by the EDBT dipole (Fig. 4.7). Thisconstraint overlaps with item 1, but we list it again to emphasize itsimportance.3. beam loss ≤ 10−5 per meter (or 0.001% m−1): additional beam losscondition imposed in arc 2. Lasing creates a momentum tail whichcan cause beam loss in the arc dipoles. Particle tracking is used in arc2 to enforce this condition, in addition to item 1. The value 10−5 m−1is the criterion for the existing E-linac design and we reuse it for theERL recirculation loop.Free parameters in the optimization run are the RF amplitudes, RFphases, and lattice optics. The complete list of free parameters andtheir search ranges can be found in Appendix C. List of initialbeam parameters and constants are also listed in Appendix C.Design ConsiderationsThe major design consideration is the E-linac itself. While most of therecirculating lattice are free variables, the E-linac design is concrete, and allupgrades must be compatible with it. The following were considered whenchoosing the initial layout.ˆ Space for the machine is limited. The ERL is housed in the hall wherethe E-linac is assembled. The size of the hall cannot be changed.ˆ The size of the machine must match the linac length of 8.31 m, ex-tending from the beginning of cavity 1 to the exit of cavity 4.ˆ A four-dipole design is chosen for both arcs 1 and 2. The highernumber of dipoles provides greater tunability. Arcs with greater thanfour dipoles were not considered as they require larger space and cost.ˆ Totally symmetric arcs and chicane dispersion suppression: the first180◦ transport is inherently not symmetric because the separator causesthe beam to enter the arc at a 7◦ off-angle and non-zero dispersion.We choose to restore symmetry by adding a mirror separator (Fig.4.1), which is a dipole-quad-dipole system which mirrors the real sep-arator. The inclusion of the separator and mirror separator means thefour-dipole arc 1 has a combined bending angle less than 180◦. While28Machine Layoutsymmetry is not compulsory, it is a desired transport characteristicand lends intuitiveness in tuning.ˆ All eight arc dipoles have identical geometry for simplicity and practi-cality. The arc 1 dipoles have smaller bending angles (see above point),achieved by operating at a lower gradient than the arc 2 dipoles.ˆ Initial beam parameters are taken from existing E-linac transport de-signs.ˆ While in theory all physical parameters can be included as free param-eters for the optimization software to sort out, our baseline choicesmust be grounded in reality. For example, the upper limit of thequadrupole gradient is limited by magnet designs already used for theE-linac. This removes immediate spurious and infeasible solutions.ˆ A minimum spacing of 25 cm between optical elements. This providesspace for support and diagnostic equipment.ˆ A minimal set of constraints and objectives are chosen to demonstratethat local micromanagement is not necessary. Top level requirementsforce local parameters into place. For example, we do not specifymatching conditions at the undulator, only requiring that FEL gain ismaximized.Machine LayoutThe ERL component layout is shown in Fig. 4.1. The translation of machineto simulation topology is shown in Fig. 4.2. Note that the linac sections areshared between three beams: ERL pass 1, ERL pass 2, and RIB. Parametersin these sections are coupled and difficult to optimize individually withoutthe global optimization platform.ERL section information is listed in Table 4.1. Several sub-systems inthe ERL modeling make use of multiple engines. This allows the drawbacksof one engine to be offset by a different engine, and is made possible by thevariable engine combination feature of the optimization platform.A combination of MADX and DIMAD is used to model sections of themachine without energy changes. The MADX output files contain numeri-cal values with greater number of significant digits and behave better thanDIMAD for text parsing. For example, an ill-behaved design can producedispersion up to 1015. MADX displays this large value, whereas DIMAD29Machine Layoutdisplays the string **********. Although dispersions of such large magni-tude are not physically realistic, it is useful for optimization as it providesinformation on search direction.The section from separator to FEL matching is commonly referred to as‘a1’. Likewise, the section from arc 2 matching to linac pass 2 matching isreferred to as ‘a2’. The modeling information for the RIB beam is listed intable 4.2.Table 4.1: List of ERL sections and their modeling information. The sectionsare listed from upstream to downstream.Section Modeling Tool(s) DescriptionLIN11 EMLinac cavity 1, ERL pass 1,including inter-cavity driftLIN12 EM Linac cavity 2, ERL pass 1EABT1 MADX,DIMAD ERL pass 1LIN13 EMLinac cavity 3, ERL pass 1,including inter-cavity driftLIN14 EM Linac cavity 4, ERL pass 1SEP MADX,DIMAD RF separatorARC1 MADX,DIMAD First arcSEP2 MADX,DIMADMirror separator to cancelthe effects of the separatorCHI MADX,DIMAD Compression chicaneFELM MADX,DIMAD FEL matching sectionUND Genesis 1 m undulatorA2M MADX,DIMAD Arc 2 matching sectionARC2 MADX,DIMAD Return arcMERG MADX,DIMAD 3-dipole mergerLIN21 EMLinac cavity 1, ERL pass 2,including inter-cavity driftLIN22 EM Linac cavity 2, ERL pass 2EABT2 MADX,DIMAD ERL pass 2LIN23 EMLinac cavity 3, ERL pass 2,including inter-cavity driftLIN24 EM Linac cavity 4, ERL pass 2EDBT MADX,DIMAD ERL dump line30MachineLayoutFigure 4.1: ERL baseline layout. The modeling of the machine starts at the first pass of the linac. Also includedis the RIB beam which operates simultaneously with the ERL beam. The linac consists of four cavities arrangedin two cryomodules EACA and EACB, separated by a transfer section EABT. The photofission beam exits theseparator to the high energy transport EHAT. The accelerated ERL beam exists the separator into arc 1. Thedecelerated ERL beam is disposed in the dump transport EDBT.31MachineLayout(MADX) a1 (DIMAD) a1 (EM) FEL beam (Genesis) FEL light (MADX) a2 (DIMAD) a2 (MADX) EABT2 (DIMAD) EABT2 (EM) lin22 (EM) lin21 (EM) lin24 (EM) lin23 (MADX) EDBT (DIMAD) EDBT END (MADX) EABT1 (DIMAD) EABT1 (EM) lin12 (EM) lin11 (EM) lin14 (EM) lin13 ERL BEGIN Coupling: shared linac parameters (MADX) EABT3 (DIMAD) EABT3 (EM) lin32 (EM) lin31 (EM) lin34 (EM) lin33 Coupling: shared linac parameters (MADX) EHBT (DIMAD) EHBT RIB BEGIN END Figure 4.2: The ERL represented by the simulation tools used to model each section. The subsystems are combinedtogether to form the optimization problem topology. The optimization platform allows engines to be combinedin arbitrary order, parallel or serial. The accelerating elements are modeled by the Empirical Model, for whichMADX and DIMAD analytical models are less than ideal for E-linac energies. Allowing for an arbitrary topologyof different modeling tools is a key and novel feature of the optimization platform.32Machine LayoutTable 4.2: List of RIB sections and their modeling information. The sectionsare listed from upstream to downstream. The linac is shared with the ERLbeam.Section Modeling Tool(s) DescriptionLIN31 EMLinac cavity 1, RIB pass,including inter-cavity driftLIN32 EM Linac cavity 2, RIB passEABT3 MADX,DIMAD RIB passLIN33 EMLinac cavity 3, RIB pass,including inter-cavity driftLIN34 EM Linac cavity 4, RIB passEHAT MADX,DIMAD Transport to RIB targetLinac Pass 1The 1.3 GHz linac [63] is composed of four superconducting niobium cavitiesaccelerating the 100 MHz ERL beam to 45 MeV and the 650 MHz photofis-sion beam for RIB production to 50 MeV going into the separator. A quadtriplet (EABT) is inserted between cavities 2 and 3. The two beams operateconcurrently, occupying different RF buckets (RF periods with no beam).The linac layout is shown in Fig. 4.3.Linac length = 8.31 m Cavity 1 Cavity 2 Cavity 3 Cavity 4 1.28 m (same for all cavities) 1.89 m 0.65 m 0.65 m Figure 4.3: The main linac consists of four superconducting niobium cavitieswith a quad triplet in the middle. The linac is shared between ERL andRIB operations.Tracking starts at the entrance of the linac with Gaussian beams. RIBparameters are taken from current design notes, with 10 MeV injectionenergy. ERL parameters are taken from initial design estimates, with 7.533Machine LayoutMeV injection energy. Certain parameters, in particular the longitudinalparameters, are estimated from the Peking University FEL [34, 39, 67, 96],which at 35 MeV and 120 pC bunch charge are comparable to the envisionedTRIUMF ERL. The injector introduces zero or negligible dispersion into thelinac. Initial beam parameters are not variables in optimization so they donot affect existing injection tunes.A phase difference is allowed between ERL and RIB operations. Thisrepresents different time-of-arrival of the two different bunches at the linac.The relative phase differences between individual cavities are equivalent forthe two modes of operation.Transport to UndulatorThe linac exit beam is delivered to the FEL via the first arc transport. Theinitial recirculation beam is separated from the photofission beam via theseparator, an RF kicker and septum combination, both modeled as dipoles,followed by the four-dipole arc 1, a mirror separator system, compressionchicane, and a five quad matching section into the undulator. The completesection is shown in Fig. 4.4. Magnets in the recirculation loop conform toexisting TRIUMF designs [10, 80, 81].Figure 4.4: The first arc transport delivers the beam from the linac to theundulator. It includes the separator, arc1, mirror separator, chicane, andFEL matching section. The mirror separator corrects for distortions fromthe separator and restores symmetry to the system. Note that the RFSeparator 2 is a regular dipole and not an RF device.The separator introduces a bend angle of 9.62◦ and a non-zero horizontaldispersion ηx leading into arc 1. Dispersion symmetry is restored by the34Machine Layoutaddition of a mirror separator system after the fourth arc dipole. The mirrorseparator contains two dipoles designed to reflect the effects and bend anglesof the septum and RF separator. The mirror septum is a simple dipolecontaining one chamber as opposed to two of the actual septum. The mirrorRF separator is a basic non-RF dipole, but shortened to 30 cm as opposedto 80 cm for the actual RF separator, with the shortening designed to savespace and cost. The RF separator is a weak dipole and the shortening hasnegligible impact on its behaviour, with MADX simulations showing changesof 10−3 m and 10−6 in ηx and η′x, respectively.The arc dipoles and quadrupoles are designed to be totally symmetric.The bending angle is evenly divided amongst the four arc dipoles, which,together with the bending angles from the separator and mirror separator,create a 180◦ turn. Solutions with symmetries βx, βy, and ηx are desired(ηy is 0 everywhere due to the lack of vertical dipoles).A totally symmetric four-dipole chicane is used to compress the beamleading into the undulator [86]. The chicane dipoles are identical with rect-angular faces. Both the length of the chicane drifts and dipole bendingangle are variables in the optimization, allowing for a wide range of chicaneM56. No objectives or constraints are placed on chicane M56 or compressionfactor. The top level requirements of beam transport and undulator gainare used as selection criteria. A variety of M56 can be achieved by manipu-lating the arc optics. In general, arc 1 rotates the beam from expanding tocontracting and the chicane compresses the beam (Fig. 4.5). The combinedM56 of the two systems determines bunch length and peak current at theundulator (Fig. 4.6).35Machine LayoutFigure 4.5: M56 of the linac-to-undulator transport is determined by arc1 and the chicane. Due to the layout of the dipoles, arc 1 has a naturalM56 < 0. A particle with less energy (red) takes a shorter path in arc 1 thana particle with more energy (blue). The chicane has M56 > 0. A particlewith less energy (red) takes a longer path in the chicane than a particle withmore energy (blue). The total M56 determines bunch compression at theundulator.z w After linac After 180 bend Arc 1 At undulator chicane Figure 4.6: Longitudinal phase space manipulation in the first transport arc.Both the linac-to-undulator and undulator-to-linac transports are de-signed for dispersion suppression. While efficient lasing operation does notrequire dispersion suppression, it leads to modularity in the system. Non-zero dispersion couples the two arcs, adding potential complications to beamtuning. The compression chicane is a completely symmetric system with noquads, therefore is inherently dispersion-free. This reduces the dispersioncondition to purely the two transport arcs.36Machine LayoutUndulatorThe focus of the modeling is a one-pass tracking of the electron beam inthe undulator at saturation using Genesis. Our goal is to find how theundulator, at steady state, affects beam dynamics. The FEL startup regimeis not modeled, since the light-beam interaction is weakest at this stage. TheFEL at saturation imparts the greatest disturbance on the beam, and it isour interest whether this exhaust beam can be adequately transported to thedump. A complete study of the optical system requires engineering decisions,e.g. choice of reflectivity for the mirrors and output coupling, which are tooearly to adopt at this stage of the design without scientific requirements set.Rough values estimated from machines with similar parameters are used.Based on the PKU-FEL ([34]), a machine of similar scope, the undulatorparameter K is chosen to be 0.7, and the undulator period λu is chosen tobe 4 cm. The optical wavelength is estimated from the undulator equation([36, 59]):λ =λu2γ2r(1 +K2)≈4 µm (4.1)where the Lorentz factor γr is based off the elinac energy of 45 MeV. Theintra-cavity saturation power Psat used for optimization is 8 MW. The valueis taken from the Daresbury ERL Prototype ([57, 89, 92]), which has similarbeam parameters and also 4 micron optical wavelength. A theoretical esti-mate ([59] Eqn. 4.61) using a combined mirror reflectivity R = 0.99, beampower Pbeam = 0.5 MW, and number of undulator periods Nu = 25,Psat≈ 12Nu(1−R)Pbeam∼1 MW (4.2)supports the assertion of megawatt saturation power. Energy loss fromundulator radiation is also incorporated into the simulation and carried intothe return arc. The list of undulator parameters is shown in Table C.5.Transport to DumpThe return loop transport consists of a five quad matching section, arc 2,another matching section, and a three dipole merger. After the linac pass 2,the exhaust beam is disposed of in EDBT (Fig. 4.7). The four arc dipolesare geometrically identical to those of arc 1 for ease of construction anddesign, each bending the beam by 45◦. The totally symmetric design is alsodispersion free if the incoming dispersion from arc 1 is zero. The mergeris taken from existing design [30]. During deceleration, energy spread can37Machine Layoutgrow by a factor of 10 due to anti-damping, potentially complicating beamdisposal. This necessitates energy spread compression in the linac. Thelongitudinal phase space manipulation is shown in 4.8.linac last cavity dump transport EDBT separator dipole quadrupole Figure 4.7: The dump transport consists of a quad shared with the ac-celerated beams, a dipole, followed by a quad doublet. The large energyspread out of the linac is converted to beam size by the dipole, potentiallycomplicating disposal.38Lattice Time-of-Flightz w Before undulator After undulator Lasing Entering linac 2 Arc 2 Dump Linac Figure 4.8: Longitudinal phase space manipulation in the second transportarc. The undulator induces energy spread, which increases due to linac anti-damping, and can complicate beam transport. Energy spread is shaped inthe linac second pass for easier beam disposal.Special care is paid to tracking beam loss in the return arc. The undu-lator introduces a momentum tail in the bunch which can be converted tobeam loss by the arc dipoles. The tail can be seen in Fig. 5.12.No free parameters exist in the second pass of the linac. All RF param-eters and EABT quad parameters are coupled to the first pass parameters.Lattice Time-of-FlightThe lattice time-of-flight is an important variable for phase matching be-tween linac passes. The TOF depends onˆ TOF in linac pass 1ˆ TOF in recirculating latticeˆ Variable-length chicaneˆ Beam energyˆ Variable length drifts in lattice39Lattice Time-of-FlightThe last parameter is a length parameter inserted into the arcs whichallows the recirculation TOF to vary with a value between 0 and 1 RFwavelength, allowing for all possible pass 2 phases. No constraints are placedon the pass 2 phase. The phase is automatically determined by the top-levelenergy recovery objective Edmp = Ein, where Ein is injector energy andEdmp is the energy at the dump after recovery.Each transport section creates a time slip variable designating the TOFthrough that section. The method of calculating the slip depends on thetool used to model the section. The Empirical Model inherently tracksusing time, whereas the slip in Madx can be inferred by the path length s.A global variable keeps track of the cumulative time slip of the beam andconverts it to phase when needed.The four linac pass 1 cavity phases, φ11, φ12, φ13, and φ14, are allowedto vary independently, where φ is the phase when the beam is at the cavityentrance. For convenience, we define φ110, φ120, φ130, and φ140, which repre-sent the cavity phases back-propagated to the simulation beginning, i.e. thecavity phases when beam is at linac entrance (by this definition, φ11 = φ110).The second pass phases, φ21, φ22, φ23, and φ24 are calculated byφ21 = φ110 +∆t(to pass 2, cav 1 entrance)φ22 = φ120 +∆t(to pass 2, cav 2 entrance)φ23 = φ130 +∆t(to pass 2, cav 3 entrance)φ24 = φ140 +∆t(to pass 2, cav 4 entrance)(4.3)where ∆t is the globally cumulative time slip. A similar set of phases existfor the RIB beam, labeled φ31, φ32, φ33, and φ34. An optimization variable,dφERL-RIB ≡ φ31 − φ11, allows the ERL and RIB beams to be phasedindependently, but the relative phases between cavities are identical sinceboth beams share the same linac, e.g. φ330 − φ320 = φ130 − φ120.Rare Isotope Beam TransportA small section of RIB transport section EHAT is included in the modelingto show consistency between ERL and RIB operations. The first section ofEHAT is shown in Fig. 4.9. The layout of elements in EHAT conform tothe existing baseline design [31] and are not variables for optimization.40Lattice Time-of-FlightFigure 4.9: The rare isotope beam is delivered to the RIB targets throughthe EHAT line. RIB is separated from the ERL beam using the RF kickerand septum.Other ConsiderationsCoherent Synchrotron RadiationThe largest sources of coherent synchrotron radiation (CSR) in the ERLcome from the chicane and arc dipoles. CSR simulations were performedusing the tracking program CSRtrack [40, 41], with 45 MeV, 100 pC Gaus-sian bunches. The chicane was modeled with four dipoles using identicalgeometric attributes as the chicane from the optimization setup. Trackingshowed a relative energy loss of 4×10−4, a negligible value and thus CSRcalculations were excluded from the optimization setup.The arc dipoles were similarly modeled using 45◦ dipoles from the opti-mization setup. Four dipoles produced a relative energy loss of 10−4, andagain small enough to be neglected from optimization. Together this showsthat CSR is not an issue in the recirculation lattice and therefore not mod-eled.CSR imparts a z-dependent kick on the bunch which can increase thetransverse emittance in the bending plane x [21, 23]. The increase is pro-portional to the inverse bunch length 1/σz [38, 94]. This could be an issuein the compression chicane where σz can be compressed to sub-millimeterlevel.We estimate the magnitude of emittance increase in εx using CSRtrack.For a 45 MeV bunch, the emittance increase was shown to be 2% in thechicane. The increase in σz is also 2%. The increases are negligible and thusCSR is excluded from the model.41Lattice Time-of-FlightBeam BreakupBeam breakup (BBU) causes deflections in the beam when traveling througha structure where electrons give up energy to the deflecting mode. If QL, theloaded Q-value of the mode, is sufficiently high, the electromagnetic fieldsof the mode will increase until the beam scrapes an aperture, resulting inbeam loss. The threshold current Ith is given by [56, 83]:Ith =2c2e(R/QL)QLω1M12 sin(ω∆t)(4.4)where e is the electric charge, (R/QL)QL is the impedance, ω is the angularfrequency of the mode, and ∆t is the recirculation time. Extensive simula-tion based on E-Linac recirculation geometry and optics has led to an upperlimit for the characteristic impedance of any higher order modes (HOM)in the E Linac cavity of (R/Q)QL ≤ 107 Ohm [61]. This was shown toallow operation of the E-Linac ERL safely below the beam breakup (BBU)threshold for all conceivable geometry and optics.2-pass ERL BBU instability was modeled with BI [12]. E-linac HOMdata up to 4 GHz was used in the simulation. Table 4.3 shows the mostdamaging dipole modes. Simulations showed Ith ≈ 30 mA, which is greaterthan the maximum E-linac current of 10 mA. The four cavities are identicalto first order, thus for modeling purposes the HOM data was treated to beidentical for all cavities. In real life the frequencies could spread by 0.1% to1%, R/Q andQ by up to 10% [60]. In addition, the quads of the recirculationloop contain enough freedom to shapeM12 (or M34) and mitigate unwantedexcitations, with optimization results showing more than 2π range in phaseadvance, both horizontal and vertical.Table 4.3: List of the most damaging dipole modes in the E-linac 9-cellcavity. Data provided by P. Kolb [62].f (GHz) QL R/Q (Ohm)3.837065314 2.89e6 0.000167783.840685874 1.34e6 0.0005371933.844854002 0.884e6 0.0010686071.890036165 0.482e6 0.03573770842Lattice Time-of-FlightSpace ChargeA question is whether space charge (SC) is a concern for tracking, which canhave an effect for energy regimes less than 10 MeV (ERL beam starts at 7.5MeV) [75]. Tracking shows SC is not an issue. Fig. 3.3 compares trackingfrom the Empirical Model (without SC) with ASTRA (with SC), with nosignificant differences.43Chapter 5Optimization Results andTradeoff StudiesThis chapter details results from the TRIUMF ERL optimization. Ourchief goal is to understand the underlying physics of the ERL transport andtradeoffs between parameters. We detail the following:ˆ How does the RF and arc transport compress the bunch to maximizelasing at the undulator?ˆ How does lasing affect beam parameters? Does this cause problemsfor energy recovery or the control of energy spread? Can lasing leadto beam loss in the return transport? How does energy spread evolvedue to lasing and energy recovery, as this is an important parameterwhich can lead to beam loss?ˆ Competition between objectives, such as energy recovery, energy spread,and how they affect the RF.ˆ Effect of the recirculation loop time-of-flight on energy recovery.ˆ Understand the physics of optimal input beam parameters into theundulator as obtained by global optimization.ˆ Transverse phase space control and the impact of demanding arc sym-metry.ˆ How do higher order transport terms impact lasing?We also show the evolution of important ERL parameters as the opti-mization proceeds.Most plots shown in this chapter represent either the optimization pop-ulation, or a subset of the population. Each point in the plots should be in-terpreted as an instance of the ERL, i.e. a particular machine design. Thesepopulation plots are useful in illustrating the physics of the ERL, and how44Bunch Compressiondifferent design parameters play against each other. Often the plots showPareto fronts between certain parameters, illustrating design tradeoffs.Some plots show a percentage of the population in order to reduce theclutter of showing the full population.Bunch CompressionProper phase space manipulation is important to rectifying many beam dy-namics challenges [79]. Proper bunch compression is important for pro-ducing a large peak current to maximize gain. We start by developing arelationship between RF phase and gain. We link the gain to the beam pa-rameters at the end of the linac, then those beam parameters to the initialRF phase, and finally show that the RF phase is related to the gain. Beamsat the start and end of arc 1 are related byVzz,a1+ = Vzz,a1 + 2M56Vzδ +M256Vδδ,a1 (5.1)Where (+) indicates the arc end and Vij is the covariance between i andj. This arises out of Va1+ = MVa1MT [25]. Fig. 5.1 shows that the threeterms of Eqn. (5.1), related to 〈z2〉, 〈zδ〉, and 〈δ2〉, all contribute to thebunch length at the arc end. Since the bunch length impacts peak current,all three terms are important to gain.45Bunch Compression0.00 0.05 0.10 0.15 0.20 0.25 0.30Vzz,a1+[mm2][mm2]0.00 0.05 0.10 0.15 0.20 0.25 0.30Vzz,a1+[mm2][mm2]0.00 0.05 0.10 0.15 0.20 0.25 0.30Vzz,a1+[mm2][mm2]0.00 0.05 0.10 0.15 0.20 0.25 0.30Vzz,a1+[mm2][mm2]Figure 5.1: Top left: Eqn. (5.1) relates z, zδ, and δ at the arc entrance to z+,where (+) refers to the arc end. The three terms of the equation are labeledA, B, and C for convenience. The line of best fit is f(〈zz〉+) = 1.047〈zz〉++4.956× 10−10. The errors on the slope and y-intercept are 6.096× 10−2 and7.086 × 10−9, thus showing the data is consistent with a completely linearrelationship. Top right, bottom left, bottom right: neglecting the first,second, and third term from f(〈zz〉+), respectively. None of the three canproduce the desired relationship with z+, thus showing all three terms areimportant contributors to bunch length, and therefore gain. The top rightplot, although linear, results in the impossible situation of a point-beamcreating a finite-sized beam (point-beam occurs when the y-axis is 0).The dependence of beam parameters on the linac acceleration phase φ1is shown in Fig. 5.2. As shown in Fig. 5.1 the beam parameters after thelinac impacts peak current at the undulator. We therefore expect φ1 toimpact lasing. Fig. 5.3 and 5.4 illustrate the effect of φ1 on peak currentand gain.46Bunch Compression324 326 328 330 332 334 336 338 340 342φ1[deg]〈zδ〉a1[µm]324 326 328 330 332 334 336 338 340 342φ1[deg]δa1[10−3]Figure 5.2: Beam parameters at the arc entrance are affected by the linacacceleration phase φ1. Top: φ1 vs zδ correlation. Bottom: φ1 vs energyspread δ. Both have roughly linear relationships with φ1 and shows theimpact φ1 can have on longitudinal beam parameters. Points shown havecavity gradients > 18 MV/m.Interestingly, there is no optimal phase corresponding to optimal gain.The two plots of Fig. 5.3 shows overlapping phase regimes produce differentpeak current Ip, and therefore gain, at different M56 (M56 is for combinedarc 1 and chicane). The top plot is a slice in negative M56 and shows apositive slope in Ip as φ1 increases. The bottom plot is a slice in positiveM56 and shows a negative slope in Ip as φ1 increases. Note the preference47Bunch Compressionfor negative M56 slice which produces higher Ip. Since large positive M56’sare not useful for compression, they are biased against in the optimization(large positive values do not even exist in later iterations, and the bottomplot uses data from an earlier iteration).322 324 326 328 330 332 334φ1[deg]5560657075808590Ip[A]310 312 314 316 318 320 322 324 326 328φ1[deg][A]Figure 5.3: The acceleration phase φ1 affects beam parameters and thereforeFEL gain. Top: peak current Ip vs φ1 with M56 in [-.15, -.05] m. Bottom:same parameters with M56 in [.3, .305] m. Ip depends on both the linac andthe arc transport, thus slices in M56 are required to isolate the effects ofφ1. The difference in slope between the two plots corresponds to how M56is matched to φ1.48Bunch CompressionThis illustrates that the combination of φ1 and M56 are needed to prop-erly compress the bunch, and that the ERL compression scheme is a com-bination of massaging the bunch in the linac, and compression in the arcand chicane transport. The highest peak current occurs near the regionM56≈− 0.1 m. The direct effects of RF on the gain are described next.Fig. 5.4 directly shows the effects of φ1 on the gain. The data forms thetypical shape of the RF curve, demonstrating the important role of RF inshaping the bunch for compression and lasing.320 325 330 335 340φ1[deg] 5.4: FEL gain vs acceleration phase φ1. Slices made on Courant-Snyder parameters at the undulator entrance: βx ≤ 2 m, βy ≤ 2 m, −0.5 ≤αx ≤ 0.5, and 0.5 ≤ αx ≤ 1.5. These slices are centered around the optimaltransverse matching conditions for the undulator (Eqn. 5.12), in order toisolate longitudinal effects. The solution set encompasses physics of thelinac, arc transport, and lasing. No single simulation tool can provide allthe physics modeling necessary.The highest gain occurs when φ1 is several degrees before crest (335◦).The baseline ERL layout presented in the next chapter has a φ1 in thisregion, roughly 10◦ before crest. Since the gain is affected by both longitu-dinal and transverse parameters, cuts were made in the transverse parameterspace to isolate the effects of φ1. The cuts made are described in the Fig.5.4 caption.The data of Fig. 5.4 includes multiple engine tracking through the linac,arc 1, chicane, and undulator, and is difficult to achieve without the com-bined modeling capabilities of the global optimization platform.49Bunch CompressionInterestingly, at the beginning of optimization, M56 can reach large pos-itive and large negative values (Fig. 5.5), courtesy of the large degrees offreedom given to the arc and the chicane. Also note that the linac is giventhe full search range for phase, on either side of the RF crest. Why then,is only a negative M56≈ − 0.1 m capable of compressing the beam? Wecan imagine a bunch accelerating on one side of the RF crest resulting in alongitudinal profile that requires a negative M56 for compression. Then thebunch accelerating on the opposing side of the RF should require a positiveM56. Fig. 5.5 should contain two peaks to either side of M56 = 0.−0.4 −0.2 0.0 0.2 0.4M56[m]0102030405060708090Ip[A]Figure 5.5: M56 of arc 1 and chicane vs peak current in an earlier generation.A negative M56 is required for compression.The reason that only one sign of M56 is selected for compression isbecause a longitudinally expanding bunch is used as the starting bunch(αz ≈ −1.5) before acceleration. A final contracting bunch is required forpositive M56 compression. The RF has a large impact on the resultinglongitudinal space, but the initial bunch is expanding too quickly and theRF cannot drive the expanding bunch to a contracting bunch (Fig. 5.6).Therefore, no positive M56 is required for compression.50Bunch Compression310 315 320 325 330 335 340 345 350φ1[deg]−3.0−2.5−2.0−1.5−1.0−0.5αz,EHATFigure 5.6: Effects of phase on longitudinal parameter αz. The RF hasa large impact in shaping the bunch, shown by the large spread in αz.However, αz started with a large negative value and the RF cannot drive itto positive.To confirm that the RF cannot drive the bunch αz from negative topositive, we estimate the RF effects on energy compression. The initialERL beam was given a bunch length σz = 0.0003 m and relative energyspread of 0.019, which translates to an absolute energy spread σE = 0.14MeV at an initial beam energy of 7.5 MeV. The RF phase range σφ occupiedby this bunch isσφ =σzfc2π = 0.008 (5.2)where f = 1.3 GHz is the RF frequency and c is the speed of light. Supposethe bunch is accelerated at φ degrees before crest. The energy spread afteracceleration is estimated from the change in particle energy E:E = E0 cosφdE = −E0 sinφ dφ≈ −E0φ dφ= −(40 MeV)(±20◦π/180◦)(0.008)= ∓0.11 MeV(5.3)where E0 ≈ 40 MeV is the energy gain in the linac if the particle is acceler-ated on-crest (7.5 MeV to 45 MeV). φ was allowed to vary in the optimization51Bunch Compression20◦ to each side of the crest. We therefore use φ = ±20◦ because these arewhere the RF slope is greatest. The phase spread was estimated dφ ≈ σφ.dE = ∓0.11 MeV is therefore the change in σE . The difference in signrepresents which side of the RF crest the bunch is accelerated at, wherethe RF slope can increase or decrease σE . Adding this value to the startingσE = 0.14 MeV, the final σE can be 0.25 MeV or 0.03 MeV. Using the initialαz = −1.5, αz after acceleration can benew αz = old αznew σEold σE= −1.50.25 MeV0.14 MeV accelerating on right of RF crest0.03 MeV0.14 MeV accelerating on left of RF crest=−2.7 accelerating on right of RF crest−0.2 accelerating on left of RF crest(5.4)The new αz range of [-2.7, -0.2] agrees with Fig. 5.6, and shows that thelinac did not drive the bunch from expanding to contracting.We test this hypothesis by reoptimizing with an initial longitudinallyupright bunch (αz = 0). αz after the linac therefore depends only on theRF phase. The results are shown in Fig. 5.7. Depending on the phase withwhich the bunch is accelerated, the output bunch can be longitudinallycontracting or expanding. Two peaks in M56 are seen, one near -0.1 m andone near +0.1 m, corresponding to the compression of an expanding andcontracting bunch, respectively. This is important, because the ERL gunis not yet designed and the exact bunch parameters at the linac are notknown. This suggests that the arc design has enough freedom to producethe necessary M56 for the compression in either direction of αz and RFphase.52Bunch Compression−0.2 −0.1 0.0 0.1 0.2M56[m]202530354045505560Ip[A]Figure 5.7: M56 vs peak current, starting with no longitudinal correlationin the initial bunch (αz = 0).Why does the RF system only produce positively chirped beams? Wewould like the RF to have the flexibility in providing any chirp required. Fora particle at coordinate z and RF angular frequency ω, the correspondingRF phase is φ = φ0 +∆φ = φ0 + ωz/c, with the bunch centroid at φ0. Weestimate the chirp magnitude at 20◦ before crest (φ0 = −20◦):M65 =∂∂zE=∂∂z(E0 cosφ)= −E0 sinφ∂φ∂z= −E0ωcsinφ≈ −E0ωcsinφ0= −(40 MeV)2π(1.3 GHz)csin(−20◦)= 372 MeV/m= 50 m−1 (scaled by 7.5 MeV initial energy)(5.5)Full RF effects on the initial energy spread of δ = 0.02 used for optimization53Effects of Lasing on Beam Transportrequires a minimumz = δ/M65= 0.02/50 m−1= 0.0004 m(5.6)By comparison, the initial bunch length σz = 0.0003 m used for opti-mization is too short for the RF to have full range-of-effect.As a result, we recommend that when a more realistic longitudinal profileof the injected beam becomes available, alternate solutions which bettermatch the M56 and RF phase are not overlooked.Effects of Lasing on Beam TransportAt the undulator, the electron-light interaction changes the electron beamproperties, including energy loss and increase in energy spread, which canhave consequences for the return transport.Lasing reduces the beam energy (Fig. 5.8). However, even with maxi-mum beam-laser interaction, the energy lost represents 0.001 of the beam(45 MeV beam losing 40 keV). This is insignificant and not a factor forenergy recovery.0.0 0.1 0.2 0.3 0.4 0.5 0.6g−40−35−30−25−20−15−10−5∆E[keV]Figure 5.8: Energy change ∆E through lasing. At maximum lasing, the 45MeV beam loses 0.1% energy. 10% of the population is shown.The lasing process also induces a beam energy spread (Fig. 5.9). Weare concerned with whether this distorted beam can be adequately energy54Effects of Lasing on Beam Transportrecovered and transported. Fig. 5.10 shows how the energy spread δ at thebeam dump and its dependence on lasing (notice the magnitude of δ beforeand after deceleration, compared with Fig. 5.9). Both figures show pointsin the same slice of the deceleration phase φ2. The increase in δ due todeceleration is evident.The lasing process induces an energy spread directly correlated with thegain. Thus, high lasing complicates beam disposal. This is unavoidablesince maximizing the gain is a top level design requirement.0.0 0.1 0.2 0.3 0.4 0.5g3.δund+[10−3]Figure 5.9: Energy spread after the undulator δund+ increases with an in-crease in gain, potentially complicating beam disposal in the return arc.Points are taken from a slice in φ2 between 150◦ and 160◦.55Effects of Lasing on Beam Transport0.0 0.1 0.2 0.3 0.4 0.5g20212223242526272829δdmp[10−3]Figure 5.10: δ at the dump transport EDBT, as a function of gain. Lasinginduced energy spread increases the final energy spread at the dump. Thisfigure is a representation of all the physics in an ERL, ranging from lasing toenergy recovery, and is made possible by the optimization platform. Pointsare taken from a slice in φ2 between 150◦ and 160◦, the same points fromFig. 5.9 after deceleration.We now explore the relationship between gain and energy recovery.56Effects of Lasing on Beam Transport0.0 0.1 0.2 0.3 0.4 0.5g7.[MeV/c]Figure 5.11: Top: gain vs energy recovered momentum. Gain and energyrecovery are two high level objectives, and it is possible to satisfy both con-currently. The red line indicates the optimal energy recovered momentumfor RF efficiency, and it is possible to maximally lase while maintainingoptimal RF efficiency. 10% of the population is shown.Fig. 5.11 shows the two top level ERL objectives of high gain and optimalenergy recovery do not compete with each other. Optimal energy recoverydepends strongly on the phase change between acceleration and deceleration,which in turn depends on the recirculation loop time-of-flight.Energy recovery depends on the 0th order transport term (TOF), whichlasing does not change, explaining the lack of correlation between the twoobjectives.Lasing distorts the longitudinal phase space, creating a momentum tail.Fig. 5.12 shows the momentum tail in the skewed normal distribution fromGenesis.57Effects of Lasing on Beam Transport−20 −15 −10 −5 0 5 10 15δund+[10−3][arbitrary units]Figure 5.12: Energy deviation of particles after lasing. A slight tail can beobserved. Its effects need to be tracked in the second transport arc.If not controlled, the tail can become an issue in the return arc transportdue to the arc dipoles, where dispersive effects can cause beam loss. In arc2 we demand three constraints:max σx < 3 mmmax σy < 3 mmbeam loss < 10−5(5.7)The third condition is not required for arc 1 because the beam in arc 1 isGaussian and has no tail.The beam loss condition is for particles whose radius exceeds 10 mm.Fig. 5.13 shows the result of imposing all three constraints in optimization.The addition of constraining beam loss proves to be a tighter condition thanconstraining beam size only, particularly for individuals with large maximumσx or σy. Although the beam envelope proves small enough, beam lossresults from the momentum tail and requires tightened global beam size forlossless transport. Fig. 5.13 shows that beam loss tightens the beam sizeconstraints to σx < 2.9 mm and σy < 2.8 mm. Fortunately, the figure alsoshows that approaching the maximum allowed beam size is not necessary toprovide a valid transport solution, with the optimization platform able tofind solutions with σx < 2 mm and σy < 1.5 mm, although such stringencyis not required.58Effects of Deceleration on Energy Recovery and Energy Spread1.8 2.0 2.2 2.4 2.6 2.8 3.0max σx, arc 2 [mm]σy,arc2[mm]beam size constraintsbeam size+loss constraintsFigure 5.13: Maximum horizontal and vertical beam size in the return arc.Blue are solutions which satisfy the first two constraints of Eqn. (5.7) re-garding beam size. Red are solutions which satisfy all three constraints,including beam loss < 10−5. The red solutions are a subset of the blue solu-tions. Included in the plot are only solutions with high gain (g > 0.3). Highgain is an optimization objective and also introduces stronger beam distor-tions, thus making the beam more susceptible to loss through the momentumtail. Imposing beam loss does constrain the arc transport further, as canbe observed at the high beam size regions where only blue solutions exist.However, plenty of individuals exist which satisfy all constraints; thereforebeam loss should not be an issue for transport.Effects of Deceleration on Energy Recovery andEnergy SpreadThe deceleration pass in the linac is important for both energy recoveryand energy compression. The deceleration phase φ2 is determined by theacceleration phase φ1 and the loop recirculation time ∆t. ∆t is allowed tovary in the optimization over more than one RF period. We now illustratethe consequences of choosing a φ2.Fig. 5.14 shows the important effects of RF phase. Both momentumand energy spread are greatly affected by the deceleration phase. The closerthe bunch enters the linac on-crest, the greater deceleration it experiences,but also leads to a decrease in RF slope and therefore less control of energy59Effects of Deceleration on Energy Recovery and Energy Spreadspread.121416182022242628δdmp[10−3]789101112pdmp[MeV/c]130 135 140 145 150 155 160 165φ2[deg]Figure 5.14: RF deceleration phase φ2 vs energy recovered momentum pdmpand energy spread at dump δdmp. This figure shows the importance of RFphase matching. The greatest decrease in pdmp occurs at the RF crest,which also corresponds to a lack of control in δdmp due to less RF slope.The two plots suggest a negatively correlated relationship between δdmp andthe ability to energy recover. Showing points with a 180◦±10◦ phase changebetween acceleration and deceleration phases.Since the time of arrival at the deceleration pass depends on the phasechange ∆φ ≡ φ2−φ1, both pdmp and δdmp show similar dependences on ∆φ(Fig. 5.15). The pdmp data agrees with the theoretical prediction that 180◦is optimal for energy recovery to 7.5 MeV/c.60Effects of Deceleration on Energy Recovery and Energy Spread121416182022242628δdmp[10−3]78910111213pdmp[MeV/c]150 160 170 180 190 200 210∆φ [deg]Figure 5.15: The change in phase between the two linac passes ∆φ is veryimportant for energy recovery. The optimal ∆φ for energy recovery is 180◦,where the plot shows almost precise recovery to the desired 7.5 MeV/c, butalso making it bad for energy compression. The figure shows 20% of thepopulation.Fig. 5.16 shows the tradeoff between parameters during deceleration,forming a Pareto front between energy recovery and energy spread. For max-imum energy recovery, we wish for the final momentum to be 7.5 MeV/c.Correspondingly, the minimum achievable relative energy spread is 0.021.Interestingly, the figure also shows that it is theoretically possible to recovermore beam energy than the linac put in, with some individuals of the op-timization population achieving less than 7 MeV/c, less than the injectionmomentum of 7.5 MeV/c. If energy spread at the dump becomes a pressingissue, it is possible to sacrifice energy recovery to achieve smaller energyspread.61Effects of Deceleration on Energy Recovery and Energy Spread7.0 7.5 8.0 8.5 9.0 9.5 10.0pdmp[MeV/c]182022242628δdmp[10−3]Figure 5.16: A Pareto front forms between energy recovered momentumpdmp and energy spread at dump δdmp, defined by the two lines. For energyrecovery, we wish to move left towards the vertical red line, which shows thedesired 7.5 MeV/c, equal to the injection momentum. At this momentumonly a 0.021 minimum δdmp is achievable. To lower δdmp, we wish to move tothe right. The two parameters are fighting against each other because theyboth depend on the RF phase. Showing same set of points as Fig. 5.14.62Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recovery7 8 9 10 11 12pdmp[MeV/c]σx,max[mm]Figure 5.17: A Pareto front forms between final momentum and EDBT beamsize, represented by the red line. Maximizing energy recovery increases theenergy spread, which leads to an increase in beam size. Optimally, we wantto move towards the lower left.Since the effect of a large energy spread is beam loss via beam scrapingin EDBT, we directly look at the relationship between energy recovery andbeam size. Fig. 5.17 shows the energy recovered momentum pdmp andthe maximum horizontal beam size in EDBT σx,max. A clear bound canbe seen, where beam size is larger as we approach the optimal 7.5 MeV/cenergy recovery point. To obtain optimal energy recovery, the decelerationpass must be driven at close to the RF crest, where less RF slope existsfor energy compression, leading to higher energy spread and thus largerbeam size. This demonstrates that having the best energy recovery requiressacrifices in beam control, and vice versa.Tradeoff Between the Recirculation LoopTime-of-Flight and Energy RecoveryFig. 5.18 illustrates the effect of acceleration phase φ1 on energy recovery.Three slices were taken in the recirculation loop time-of-flight ∆t: ∆t thatbrings the bunch back to the linac with ≈ 180◦ phase change (red), < 180◦(blue), and > 180◦ (black). These represent bunches which arrive at thelinac with perfect 180◦ timing, early, and later, respectively. Interestingly,63Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recoverywhen the energy recovered momentum pdmp and φ1 are plotted for the threeslices, we find that all three slices cross 7.5 MeV/c, which is the injectormomentum and represents our objective of Edmp = Ein for energy recovery.64Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recovery310 315 320 325 330 335 340acceleration phase φ1[deg]120130140150160170180decelerationphaseφ2[deg]127.752 ns<∆t <127.755 ns127.770 ns<∆t <127.773 ns127.790 ns<∆t <127.793 ns310 315 320 325 330 335 340acceleration phase φ1[deg][MeV/c]127.752 ns<∆t <127.755 ns127.770 ns<∆t <127.773 ns127.790 ns<∆t <127.793 nsFigure 5.18: Top: correlation between linac acceleration phase φ1 and de-celeration phase φ2. Each set represents solutions where the recirculationlattice has a particular range of time-of-flight ∆t. The red solutions repre-sent machines with an optimal phase difference of ∆φ = φ2 − φ1 ≈ 180◦,whereas the blue solutions are machines less than 180◦ (beam arrive at sec-ond pass earlier than optimal) and black solutions are machines greater than180◦ (beam arrive at second pass later than optimal). Bottom: energy re-covery as a function of φ1, for different slices of ∆t. Note that for all threeslices, solutions exist to satisfy the energy recovery objective Edmp = Ein atcertain phases. The slice represented by red is phase independent and cor-responds (see left figure) to the theoretically optimal ∆φoptimal = 180◦. Theblue and black slices show that satisfying Edmp = Ein can still be achievedat other phase settings.65Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recovery∆t and its effect on energy recovery is shown in Fig. 5.11. The role ofφ1 becomes clear once ∆t is fixed. If ∆t transports the beam such that φ1 is180◦ phase flipped from the deceleration phase φ2, optimal energy recoveryis achieved regardless of the initial φ1 due to the exact cancellation of theRF acceleration and deceleration curves, as shown by the red solutions inFig. 5.18:V2 = −V1sinφ2 = − sinφ1φ2 = φ1 + 2π(n + 0.5)(5.8)where the phases are independent of the recirculation time-of-flight ∆t. Ifa non-180◦ ∆φ is used, Edmp = Ein is still possible at specific φ1, as illus-trated by the blue and black solutions. The reason is as follows: φ1 andφ2 are coupled by ∆t, which determines the offset between the accelerationand deceleration curve (Fig. 5.19). Perfect cancellation occurs if the curvesare maximally out-of-phase (180◦ offset), or if they deviate such that φ1 iscoupled to φ2 on the opposing side of the RF trough, shown by the diago-nally dashed line. The horizontally dashed line indicates ∆φ needed for thisparticular instance. For Edmp = Ein, ∆t is expected to be linearly relatedto φ1 byφ2′ − φ1 = 2πf∆t = 2dφ+ 2π(n + 0.5)= 2(1.5π − φ1) + 2π(n + 0.5)∆t =2(1.5π − φ1) + 2π(n + 0.5)2πf= −2/(2πf)φ1 + const= −0.00427 ns/deg× φ1 + const(5.9)where the prime indicates the opposite side of the crest/trough and 2dφ isthe difference between the primed and unprimed phases.66Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recovery315 320 325 330 335 340 345 350 355accelerating phase φ1[deg][arbitraryunits]135 140 145 150 155 160 165 170 175decelerating phase φ2[deg]φ1φ2Figure 5.19: The data in Fig. 5.18 can be explained by the matching oflinac phases. Top: acceleration and deceleration by φ1 and φ2. Momentumunits are arbitrary; p = 1 is maximum acceleration for φ1 and p = 0 is max-imum deceleration for φ2. Data points are extracted from Empirical Modelinterpolation tables, which are created from Astra tracking. Accelerationand deceleration are reflections of each other with 180◦ flip, with the cavityimparting maximum effects at 335◦/155◦. The configuration correspondingto the red solutions is pictured, where the two curves are perfectly out-of-phase and RF effects cancel at all phases for optimal energy recovery. If theloop produces ∆φ 6= 180◦, Edmp = Ein can be achieved only at the pointat the mirror opposite of the RF crest/trough, illustrated by the diagonallydashed line. In this scenario, ∆φ for energy matching depends on the howfar φ1 is off-crest. The bottom of Fig. 5.18 illustrates this; the black andblue solution sets are non-180◦ offsets, and crosses the optimal 7.5 MeV/conly at specific φ1.This suggests a correlation between φ1 and ∆t for optimal energy re-covery. Eqn. 5.9 is plotted against the solution set with Edmp = Ein inFig. 5.20. Edmp = Ein is possible at all RF phases if the recirculation looptime-of-flight is tuned correctly.67Tradeoff Between the Recirculation Loop Time-of-Flight and Energy Recovery324 326 328 330 332 334 336 338 340φ1[deg]−0.04−0.03−0.02−∆t[ns+127.770ns]Figure 5.20: Relationship between φ1 and ∆t for energy matching to 7.5MeV/c. Edmp = Ein is possible for all initial phases. The fit is given byEqn. 5.9, with the theoretical slope = −0.00427 ns/deg. Included are pointswith pdmp in [7.45, 7.55] MeV/c.Fig. 5.20 seems to suggest that we can optimally energy recover at non-180◦ phase changes.It is important to note that the optimization was carried out for beamdynamics, and does not include RF loading effects. A previous study [82]has shown that non-180◦ phase matching has serious consequences for theRF system, including the need for higher RF power, due to beam vectors notcancelling and therefore causing beam loading. Another interesting avenueto pursue is that the linac supports three passes: ERL acceleration anddeceleration, and RIB acceleration. 180◦ may be optimal for ERL operation,but may not be globally optimal when operating simultaneously with RIB.The amount of deviation from 180◦ that can be tolerated by the systemdepends on the selection of the loaded Q-value and other RF parameters,and requires further studies beyond the scope of this thesis.The optimization platform is capable of finding 180◦ solutions optimalfor RF, as illustrated by the red solution set in Fig. 5.18. We intentionallychose not to set a constraint to enforce exactly 180◦ phase change in theoptimization setup, as we wanted to test the capabilities of the platformwhen we give it minimal guidance, as well as explore the repercussions ofnon-optimal phase matching, which indeed was illustrated in the tradeoffbetween energy recovery and energy compression (Fig. 5.16). It is inter-esting that the platform can find unexpected dynamics (Fig. 5.20) within68Undulator Conditions for Maximizing Gainthe system and demonstrates the advantage of the global approach overpiecewise.Undulator Conditions for Maximizing GainWe study the effects of beam conditions on the undulator gain. The theo-retical gain in an undulator is given by [59]g = −4√2π2IpIAK2[JJ ]2(1 +K2/2)3/2N2u√λλuh(χ) (5.10)where IA = 17045 A is the constant Alfve´n current, K is the undulatorparameter, [JJ ] = J0(K2/(4 + 2K2))− J1(K2/(4 + 2K2)) is a combinationof Bessel functions of the first kind J0 and J1, Nu the number of undulatorperiods, λ the radiation wavelength, λu the undulator wavelength, and h(χ)Madey’s function with χ the scaled energy deviation. Of note is the lineardependence on the peak current Ip, which was reproduced by optimizationin Fig. 5.21. Global optimization shows a maximum bunch compression(therefore maximum peak current) down to 0.15 mm.30 40 50 60 70 80Ip[A]0.300.350.400.450.500.550.600.65gFigure 5.21: Gain increases as peak current Ip increases; this is consistentwith FEL theory (Eqn. 5.10) and requires compression by the chicane. Theupper bound on the gain has the predicted linear functional dependenceon Ip, denoted by the line. The spread in gain is caused by other beamparameters.69Undulator Conditions for Maximizing GainWe look at the transverse matching conditions for the undulator, givenin the reference [34]. Horizontally the undulator resembles a drift, so wewant the incoming beam to be focusing and form a symmetric waist atthe undulator center. Vertically, dipole edge focusing, when averaged overthe undulator period λu, creates a section of constant focusing strengthK/(√2γrλu) where γr is the Lorentz factor of the bunch centroid. Thusvertically we would like a coasting beam. The matching Courant-Snyderparameters areβx,match = ZR + (L2u/4ZR)αx,match = Lu/(2ZR)βy,match =√2γr/(Kku)≈ZRαy,match = 0(5.11)where ZR is the Rayleigh length, Lu is the undulator length, and ku is theundulator wavenumber. Definitions of α and β are presented in AppendixA. Qualitatively, β is related to the beam size in x or y, and α is related tothe beam tilt (in x or y phase space). Using ZR = 0.5 m, Lu = 1 m, thematching conditions can be estimated asβx,match = 1 mαx,match = 1βy,match = 0.5 mαy,match = 0(5.12)The matching conditions produced by global optimization can be seenin Fig. 5.22, which match closely with values given by Eqn. 5.12. 2D plotsshowing the same matching conditions are shown in Fig. 5.23. The globalresults can be contrasted with results from local (standalone) optimizationin Fig. 5.24.Although α and β values agree closely, the gain values suggests an ad-vantage of the global optimization scheme over local optimization. The gainproduced by global optimization (Fig. 5.22) is higher than the gain fromthe standalone undulator optimization (Fig. 5.24). Global optimizationdoes not impose a limit on the peak current Ip, which is obtained from RFand arc tracking.On the other hand, local optimization starts at the undulator withoutthe preceding RF and arc sections, thus we require an initial search range forIp. The search range must be estimated and can impose an arbitrary limit.In particular, the values used to produce Fig. 5.24 were based on the PekingUniversity FEL [34]. The TRIUMF arc resulted in a higher Ip, and if the70Undulator Conditions for Maximizing GainPKU number was used, we would have arbitrarily limited the performanceof the undulator.This is a demonstration of the superiority of the global scheme. If thearbitrary upper limit for Ip is too high, the upstream transport can neverproduce it. If the limit is low, the performance is suboptimal. Globaloptimization solves the problem by letting the system dynamically choosethe correct Ip for gain and transport.71Undulator Conditions for Maximizing GainFigure 5.22: Matching conditions of the undulator in horizontal x (top) andvertical y (bottom), as produced by global optimization. The results can becontrasted with the standalone undulator optimization (Fig. 5.24), whichmatch closely. The results also match with the theoretical conditions definedby Eqn. 5.11.72Undulator Conditions for Maximizing Gain0 1 2 3 4 5 6 7 8βx[m]−3 −2 −1 0 1 2 3 4αx0. 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0βy[m]−6 −5 −4 −3 −2 −1 0 1 2αy0. 5.23: Another view of the undulator matching conditions. The datamatch very well with theory (βx = 1 m, αx = 1, βy = 0.5 m, αy = 0).73Undulator Conditions for Maximizing GainFigure 5.24: Matching conditions of the undulator in horizontal x (left)and vertical y (right), as produced by local (standalone) optimization ofGenesis. The standalone results match closely with global optimizationresults (Fig. 5.22), which validates the global ERL model. One aspectof the global scheme not reproduced by the standalone is the maximumgain (0.6 for global and 0.3 for standalone). This is because we need toexplicitly specify the input peak current Ip in the standalone case. Sincethe rest of the beamline is not designed, the range of Ip is not known andcan only be estimated (in this case, underestimated, leading to lower gain).The global scheme implicitly produces the Ip range from tunable knobs (suchas RF phase and quad strengths) which we can directly control in practice.Therefore, the derived Ip is more on target with the machine design.The dependence on small β for high gain also sets the focusing conditionsof the quads in the FEL matching section. Fig. 5.25 shows strong gradientsin the two quads immediately upstream of the undulator. All quads on thelattice are given bipolar freedom and their strengths left to the optimizationplatform to choose.74Transverse DynamicsFigure 5.25: Achieving a small beam size requires strong focusing in theFEL matching section. Plotted is K1 for the matching section quads Q4and Q5, the two quads prior to the undulator. Global optimization choseQ5 to be horizontally focusing, supporting the theory that the best matchfor the undulator is a horizontally converging beam (Eqn. 5.12). Q4 ishorizontally defocusing to create a defocusing-focusing section (D0F0).Transverse DynamicsWe look at how transverse objectives such as symmetry can impact beamtransport. Symmetries in βx, βy, and ηx were casted as objectives by re-quiring αx = 0, αy = 0, and η′x = 0, respectively, in the centers of both arc1 and arc 2 (six objectives total). ηy = 0 everywhere and thus inherentlysymmetrical.The beam size after acceleration is much larger vertically than horizon-tally (Fig. 5.26). As a reminder, we require less than 3 mm RMS beam sizeat all points.75Transverse Dynamics0.95 1.00 1.05 1.10 1.15 1.20 1.25 1.30σx, after linac pass 1 [mm]2.552.602.652.702.752.802.852.902.95σy,afterlinacpass1[mm]Figure 5.26: Horizontal and vertical RMS beam size after linac pass 1. Thevertical size is significantly larger and requires immediate focusing. Pointswith g > 0.3 shown.This requires the first quad after the linac, EHATQ1, to be verticallyfocusing (negative K1). Without immediate vertical focusing, the beam sizecan easily increase above the acceptable limit of 3 mm (Fig. 5.27).−1.0 −0.8 −0.6 −0.4 −0.2K1,EHATQ1[m−2]2.702.752.802.852.902.953.00maxσy,linactoarc1center[mm]Figure 5.27: Effect of EHATQ1 on vertical beam size in the arc. The Y-axisshows the maximum vertical size from the linac exit to the center of the firstarc. Points with g > 0.3 shown.76Transverse DynamicsThe dump section EDBT consists of the quad EHATQ1, a horizontalbending dipole to kick the beam to the dump transport, followed by a quaddoublet (Fig. 4.7). Control of horizontal beam size σx is critical to properbeam disposal so as not to lose beam via beam scraping. As was shown pre-viously, the energy spread increases dramatically after deceleration, whichis converted to horizontal size. Therefore, other contributions to beam sizemust be minimized. EHATQ1 cannot be strongly horizontal defocusing.Otherwise, a large beam going into the dipole becomes even larger due toenergy spread, and easily violates the beam size condition. Fig. 5.28 showsthe effect of EHATQ1 on maximum EDBT σx. A bound forms (blue line), inwhich the more horizontal defocusing EHATQ1 is, the more difficult limitingbeam size becomes.−1.0 −0.8 −0.6 −0.4 −0.2K1,EHATQ1[m−2]σx,EDBT[mm]Figure 5.28: EHATQ1 vs max horizontal beam size in EDBT. Too muchhorizontal defocusing going into the dipole causes problems because we ex-pect the beam size to get even larger after the dipole due to the high energyspread. Points with g > 0.3 shown.The above shows competing criteria for the settings of the quad EHATQ1.1. The vertical beam size is large coming out of the linac first pass. Im-mediate beam size control is required. For vertical focusing, EHATQ1is preferred to have a negative K1.2. For proper dump transport, the beam cannot be too large horizontally.Otherwise, along with the increase in horizontal size from dispersion,77Transverse Dynamicsthe beam becomes difficult to control. For this reason, EHATQ1 ispreferred to have a positive K1.The competing objectives between arc beam size and dump beam sizeforms a Pareto front shown in Fig. 5.29. EHATQ1 is pulled in both di-rections and the optimization platform settles on EHATQ1 being close toneutral. Both Fig. 5.27 and 5.28 show that EHATQ1 is concentrated near0, with a slight negative lean, showing that the arc beam constraint is moreurgent.This tradeoff in beam size should not cause a problem for beam transportas there exist solutions in which the beam size is under 3 mm RMS for botharc 1 and EDBT.2.70 2.75 2.80 2.85 2.90 2.95 3.00max σy, linac to arc 1 center [mm]σx,EDBT[mm]Figure 5.29: Beam size control in arc 1 vs EDBT. Arc 1 vertical beam sizeforms a Pareto front with EDBT horizontal beam size due to the opposingrequirements on the shared quad EHATQ1. The front is denoted by theline. σx in EDBT and σy in arc 1 cannot be minimized simultaneously.Fortunately, this is not an issue because we only require both < 3 mm,not minimized. The plot shows many solutions which satisfy this condition.Points with g > 0.3 shown.The large vertical beam size out of the linac also has an effect on sym-metry. We represent the βy symmetry by αy = 0 at the arc center. Fig.5.30 shows that if we want M56 to be near the desired value of ≈ −0.1 mfor bunch compression, we have to take away from symmetry. To explainthis refer to the sample lattice in Fig. 5.31. Symmetry requires the points78Transverse DynamicsA and A′ to be equal. This is impossible, because if A′ was made to be aslarge as A, then it would immediately be made even larger by the verticallydefocusing chicane, and the beam size constraint would be violated.−0.08−0.06−0.04−0.020.00 0.02 0.04 0.06 0.08M56, linac to und [m]−8−7−6−5−4−3−2αy,arc1centerFigure 5.30: M56 vs vertical beta function symmetry. For optimal lasing,M56 ≈ −0.1 m. This opposes βy symmetry defined by αy = 0. The expla-nation of the mechanism is given by Fig. 5.31.020406080100βx(black),βy(red),[m]0 5 10 15 20s [m]chicaneAA'Figure 5.31: βy (red) symmetry is constrained because of the large initialvertical beam size at point A. For βy symmetry to hold, the point A′ mustincrease to the same size as A (reflection around the arc 1 center at s ≈ 7m). This large beam size is made larger by the vertically defocusing chicane,violating the beam size constraint. The chicane begins near s ≈ 15 m. Forcomparison, βx (black) does not exhibit the same issue.79Higher Order Transport EffectsThis tradeoff between M56 and βy symmetry leads to a Pareto frontbetween βy symmetry and peak current (Fig. 5.32). For high gain operation,vertical symmetry has to be broken. This is acceptable because symmetry isnot required for the machine to function (besides, breaking one out of foursymmetries is not bad). In principle, this can be fixed by tuning the injectortransport to create a different set of beam conditions at the arc.−8 −7 −6 −5 −4 −3 −2αy, arc 1 center202530354045505560Ip[A]Figure 5.32: A Pareto front forms between the vertical beta function sym-metry and peak current Ip at the undulator, shown by the line. For βysymmetry we would like αy = 0. This requires moving towards the right inthe plot, but doing so would decrease Ip and thus the gain. High gain andαy = 0 (top right in the plot) are not achievable simultaneously. This effectis caused by optical requirements explained in Fig. 5.31.Higher Order Transport EffectsHigher order terms can contribute to linac transport due to the nonlinearRF slope. Fig. 5.33 shows the bunch longitudinal phase space after accel-eration with two different energies, 7 MeV and 40 MeV. Neither case showsappreciable higher order effects, and this justifies using first order envelopetracking in the cavities.80Higher Order Transport Effects−1.0 −0.5 0.0 0.5 1.0z [mm]−20−1001020δ[10−3]−1.0 −0.5 0.0 0.5 1.0z [mm]−10−50510δ[10−3]Figure 5.33: Longitudinal phase space after acceleration in cavities 1 (left)and 4 (right). Higher order effects are minimal, therefore first order enve-lope tracking is sufficient. Tracking performed using the Empirical Model.Note: shown here are the phase spaces for a particular ERL design, not thepopulation.The effects of transport map element T566 is of concern due to dipolesin the beamline. In arc 1 and the chicane, we are concerned with bunchcompression. To compare first and second order contributions to bunchlength, define the metricr56 =∣∣∣∣T566δM56∣∣∣∣ (5.13)where δ is the energy spread out of the linac. Fig. 5.35 shows the effectof T566. Fig. 5.34 shows peak current Ip, or inversely, bunch length, as afunction of the linear term M56. Maximum Ip, and thus gain, occurs nearM56 ≈ −0.08 m, therefore we limit the analysis near this region. This avoidsthe region M56 = 0, where Eqn. 5.13 breaks down.81Higher Order Transport Effects−0.10 −0.08 −0.06 −0.04 −0.02 0.00 0.02arc 1 + chicane M56[m]405060708090Ip[A]Figure 5.34: Peak current Ip occurs near M56 ≈ −.08 m. The highest gaintherefore occurs in the same range. Points shown have g > 0.3.Fig. 5.35 shows that r56 is at maximum 0.1, thus linear transport is atleast an order of magnitude stronger than the second order term. Translatedinto bunch length, the second order contribution T566δ2 is at most 0.026 mm,roughly 20% of the bunch length at the undulator (optimization shows thesmallest bunch length is 0.15 mm). Although small, the contribution is notnegligible.82Higher Order Transport Effects0.06 0.07 0.08 0.09 0.10 0.11arc 1 r56[m]0.0160.0180.0200.0220.0240.0260.0280.030arc1T566δ2[mm]Figure 5.35: The horizontal axis shows the metric r56 as defined by Eqn.5.13. Small bunch length is important for achieving high gain (Fig. 5.21);the vertical axis shows the second order bunch length contribution T566δ2.In the M56 range [-.1,-.07] m where maximum gain occurs, the maximumvalue of r56 is 0.1, demonstrating that in the optimal gain regime, the linearterm dominates. T566δ2 contributes up to 20% of the bunch length at theundulator (optimization shows compression down to 0.15 mm). Thereforethe second order effect, while not dominating, is non-negligible. To isolatethe second order effects, points shown have −.1 m < M56 < −.07 m.Fig. 5.36 shows the direct effect of T566 on gain. The plot is boundedby an inverse relationship; large T566 leads to larger bunch length, thereforelower gain. Minimizing T566 is desirable, but unnecessary, since requestingmaximum gain automatically forces global optimization to minimize T566.83Higher Order Transport Effects1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0arc 1 and chicane T566[m] 5.36: Direct effect of T566 on gain. Although bunch length at theundulator is primarily determined byM56, T566 provides non-negligible con-tributions. The figure shows gain is bounded by an inverse relationship(albeit fuzzy due to the dependence of gain on other parameters). High T566is not accommodating to high gain. Therefore a criterion for high gain isminimizing T566. This is automatically enforced in the optimization by thegain maximizing objective. Points shown are the same as in Fig. 5.35.In arc 2 we are concerned with the effects of T566 and its effect on energyrecovery. Fig. 5.37 shows that neither the linear nor the second order termhas a dominating effect on energy recovery, i.e. optimal energy recovery ispossible for a range of transport schemes. The top figure shows the lineareffect. Energy recovery near M56 = 0 is possible, and Eqn. 5.13 is not auseful metric in this scenario.Energy recovery is not dependent on T566; therefore higher order trans-port in arc 2 is not a concern. This is in agreement with previous results,which showed that the primary role of arc 2 is timing the bunch to RF phase,the phase being the primary determinant of optimal energy recovery (Fig.5.11).84Higher Order Transport Effects−0.25−0.20−0.15−0.10−0.05 0.00 0.05arc 2M56[m][MeV/c]2.0 2.5 3.0 3.5 4.0arc 2 T566[m][MeV/c]Figure 5.37: The objective of arc 2 is optimal energy recovery. Unlike arc1, which needed a specific M56 to achieve high peak current for lasing, arc2 can have a range of M56 and still achieve optimal energy recovery. Top:M56 can cross 0, thus using the metric r56 from Eqn. 5.13 is not useful.Bottom: effects of T566 on energy recovery. Optimal energy recovery ispossible for all values of T566, thus T566 should not be a concern. The rangeof T566 values shown in the figure spans the entire range of T566 exploredby the optimization platform, i.e. possible combinations given the initialoptimization parameters. There is no indication of a parasitic T566 value,i.e. a T566 value which causes inadequate energy recovery.Higher order terms do have an effect on energy spread. Fig. 5.38 showsthe effects of arc 2 optics on energy spread δ in the dump line. The top85Higher Order Transport Effectsfigure displays the typical RF curve, demonstrating that M56 affects RFmatching. The bottom figure shows large T566 can negatively impact beamdisposal by mismatching the bunch to the RF curve, resulting in less thanoptimal energy compression.−0.25−0.20−0.15−0.10−0.05 0.00 0.05 0.10arc 2M56[m]1618202224262830δdmp[10−3]2.0 2.5 3.0 3.5 4.0arc 2 T566[m]1618202224262830δdmp[10−3]Figure 5.38: Arc 2 transport can shape the energy spread δ at the dump.M56 (top) and T566 (bottom) shapes how the bunch is matched with the RFcurve and therefore influence energy spread compression in the linac. Theline in the bottom figure shows that a large T566 value is correlated withlarger δ which can potentially cause beam disposal problems. Points withg > 0.2 are shown.86Evolution of the Optimization PopulationEvolution of the Optimization PopulationHere we show how the population evolved over the course of the optimizationrun. This is useful in determining whether the algorithm is working correctlyand finding better Pareto fronts. The state of the run is measured in numberof iterations, or generations.Fig. 5.39 shows the evolution of several key parameters: gain, energyrecovered momentum, maximum beam size in EDBT, and the dump energyspread. The run at three different generations are compared. Generation Ais near the start of the run, B some time after, and C is the latest generation,or the generation we terminated the run at.0.0 0.1 0.2 0.3 0.4 0.5 0.6g050010001500200025003000# of solutions in populationgen Agen Bgen C6 7 8 9 10 11 12pdmp[MeV/c]0100200300400500600# of solutions in populationgen Agen Bgen C2.0 2.5 3.0 3.5 4.0max σx,EDBT[mm]02004006008001000# of solutions in populationgen Agen Bgen C14 16 18 20 22 24 26 28 30δdmp[10−3]0100200300400500600700# of solutions in populationgen Agen Bgen CFigure 5.39: Parameter evolution in the optimization population. From topleft and clockwise, the parameters are gain, dump momentum, dump energyspread, and maximum horizontal beam size in EDBT.The plot of gain shows that the population starts with almost all lowgain individuals, and the population evolves towards maximizing the gain.The plot of maximum beam size σx,EDBT shows designs evolving towardsthe 3 mm boundary in order to satisfy the beam size constraint. In gener-ation A, many designs are larger than the acceptable 3 mm. By generation87Evolution of the Optimization PopulationC almost all are within the 3 mm limit.The energy spread δdmp moves toward smaller values in response to theEDBT beam size constraint.The plot of pdmp shows that the number of designs from ≈ 7 MeV/c to≈ 12 MeV/c are increasing. The number of designs near the optimal 7.5MeV/c point is increasing because this best satisfies the energy recoveryobjective. The number of designs above optimality is also increasing due tothe tradeoff between pdmp and δdmp, as Fig. 5.15 demonstrates. High pdmpdesigns also tend to have low δdmp, thus more likely to satisfy the σx,EDBTconstraint.The optimization run can be completed in 30 days running over 60 com-puting nodes. The total number of ERL evaluations is roughly one million.88Chapter 6ERL Baseline DesignThis chapter outlines a particular solution taken from the optimization pop-ulation. The solution obeys all beam transport constraints, taking into ac-count the effects of energy spread and beam scraping in EDBT. The choicewas primarily dictated by high gain, followed by energy recovery. Symme-try in transport functions are a nice-to-have and are present in the solutionchosen.This specific design of the ERL beamline is described. Specific infor-mation on beam parameters, layout, and coordinates table can be foundin Appendix D. The baseline comprises the following sections (refer to thelayout diagram Fig. 4.1):1. cryomodule EACA pass 1 (containing cavities 1 and 2)2. EABT pass 13. cryomodule EACB pass 1 (containing cavities 3 and 4)4. high energy transport EHAT, including RF separator and septum5. first arc ARC16. mirror septum SEP27. chicane CHI8. undulator matching FELM9. free electron laser FEL10. arc 2 matching A2M11. second arc ARC212. merger MERG13. cryomodule EACA pass 214. EABT pass 215. cryomodule EACB pass 216. ERL beam dump EDBT17. RIB transport EHAT (alternate path after septum)The purpose for these sections of the E-linac is to transport and accel-erate a 100 pC/bunch ERL beam through EACA-EABT1-EACB from 7.5MeV to ≈ 45 MeV. The beam is transported through the SEP-ARC1-CHI-SEP2-FELM sections to arrive at the undulator FEL. Lasing follows alongwith slight energy loss. Beam transports through A2M-ARC2-MERG to89Beam Transportarrive at the linac. Pass 2 EACA-EABT-EACB decelerates the beam to 7.5MeV and beam is disposed in EDBT.A separate 16 pC/bunch beam accelerates through EACA-EABT-EACB-SEP-EHAT to the RIB photofission target from 10 MeV to ≈ 50 MeV.Optimization Parameters for the BaselineTop level parameters are listed below. Gain is defined as (dP/dz)/P , whereP is the radiation power.Table 6.1: ERL baseline parameters.Parameter ValueGain 0.5 m−1Initial momentum 7.5 MeVEDBT momentum 7.7 MeVσx ≤ 3 mm everywhereσy ≤ 3 mm everywhereEDBT energy spread 0.029EDBT max σx 3.0 mmEDBT max σy 1.9 mmDump σx 5.5 mmDump σy 6.0 mmBeam loss ≤ 10−5The design has a gain of 0.5 m−1. This is near the top of the optimizationsearch space for gain and satisfies our “maximize lasing” objective.The EDBT max σx is within our constraints, demonstrating that energyspread is contained and should not be an issue.Beam TransportHere we show the baseline optics functions (see Appendix A for definitions)of the ERL recirculation lattice. The transport solution obeys beam sizeconstraints everywhere, including beam loss ≤ 10−5. A detailed latticelayout coordinates table is shown in Table D.8 of Appendix D.90Beam TransportOptics functions of EHAT to FEL is shown in Fig. 6.1. From FEL toMERG is shown in Fig. 6.2. EDBT is shown in Fig. 6.3.91BeamTransport−2.0−1.5−1.0−ηx(green)[m]020406080100βx(black),βy(red),[m]0 5 10 15 20 25s [m]Figure 6.1: Optical functions of the linac to undulator transport. Arc 1 is centered at s ≈ 7 m. βx and ηx showsgood symmetry as desired. βy does not show good symmetry, as explained previously by Fig. 5.31. The chicanestarts near s ≈ 15 m.−2.0−1.5−1.0−ηx(green)[m]020406080100βx(black),βy(red),[m]0 2 4 6 8 10 12 14s [m]Figure 6.2: Optical functions of the undulator to linac pass 2 transport. All three functions plotted exhibit signsof symmetry in arc 2 (centered at s ≈ 5.5 m).92BeamTransport−2.0−1.5−1.0−ηx(green)[m]020406080100βx(black),βy(red),[m]0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5s [m]Figure 6.3: Optical functions of EDBT transport. Although ηx looks contained compared to Fig. 6.1 and 6.2, itseffect on the beam size σx = ηxδ can be large due to large energy spread δ in this section.−2.0−1.5−1.0−ηx(green)[m]020406080100βx(black),βy(red),[m]0 1 2 3 4 5 6 7s [m]Figure 6.4: Optical functions of EHAT transport for RIB photofission beam. Optimization shows viable transportsolutions for the RIB photofission beam, illustrating that simultaneous ERL operation does not interfere with RIBoperation.93ERL Compatibility with RIBERL Compatibility with RIBERL operations do not interfere with RIB operations. The simultaneousRIB transport solution in EHAT is shown in Fig. 6.4.ERL Compatibility with Energy DoublingThe first phase of ERL construction may see the possibility of building a barerecirculation lattice without the undulator. We wish to know whether thelattice can be used in energy doubling mode, labeled Recirculating LinearAccelerator (RLA) mode, before the undulator is added and full FEL op-eration begin. This was not part of the original optimization requirements,thus no provisions were taken to include it in the optimization. In principle,RLA can be added to the optimization in a similar manner as RIB.RLA operates by transporting the bunch to the linac pass 2, but with0◦ phase change from pass 1, for acceleration in both passes. In contrast,ERL requires 180◦ phase change. Therefore, compatibility with energy dou-bling requires that the ERL recirculation lattice can provide half an RFwavelength λ of freedom in path length.We begin by examining whether the chicane can provide this freedom.We add 0.25λ ≈ 6 cm to each of the first and third drift lengths L (see Fig.6.5). This results in a horizontal offset of 14 cm in the second drift betweenthe ERL and RLA beam trajectories. This is too much for the beam pipeand thus the chicane cannot provide the path length difference.Figure 6.5: Chicane in ERL mode (straight line) and RLA mode (dashedline).We examine a second possibility of replacing the chicane with a drift.The path length reduction is 2×6 cm ≈ 0.5λ (see Fig. 6.6), which is roughly94ERL Compatibility with Energy Doublingthe amount required for the phase change. Further small changes can beperformed by tuning M56 of the arcs. Thus if RLA operation is desiredbefore the full ERL is operational, we recommend leaving out the chicane inthe first phase of the recirculation loop construction and replacing it with adrift.Figure 6.6: Effect of removing chicane on path length.95Chapter 7ConclusionsThis work resulted in the creation of a generic optimization platform andthe creation of a baseline ERL design. This is an exciting development for itnot only is the first ERL in Canada, but the first combined functions ERLand RIB accelerator in the world, and highlights Canadian achievements inaccelerator physics. The baseline design is important for two reasons:1. A complete list of elements and coordinates of the ERL lattice wasprovided. It lays the foundation for the upgrade of the TRIUMF E-linac to a light source. Clear relationships and tradeoffs were shownregarding this specific design. Empirical relationships, such as thatof energy recovery and dump energy spread were derived, and theywill be of great value to physicists when the full ERL design is un-derway. In addition, all optical elements were designed in accordancewith existing magnet designs. This simplifies the design and does notrequire new studies in alternative magnets. The baseline also con-forms to engineering constraints. Space is set for support technologyand diagnostic devices to be inserted.2. The wider importance of the work is that it offers a study of a generalclass of ERLs. The physical processes of gain, recirculation, energyrecovery, and beam disposal are present in all ERLs in the world.Even physically larger ERLs with more energetic beams encounterthe same issues, as can be observed in the Continuous Electron BeamAccelerator Facility (CEBAF) at Jefferson Lab. The interplay betweenthese varied objectives were empirically demonstrated in the first start-to-end ERL optimization, and can be applied to any machine with alinac driver, a recirculation loop with a free electron laser, and energyrecovery.Some dynamic relationships that arose from the start-to-end ERL studyare:1. Lasing increases energy spread, and therefore complicates beam dis-posal.96Chapter 7. Conclusions2. Tradeoff between energy recovery and energy spread. We cannot havethe best of both worlds.3. Energy recovery is affected by the interplay between RF phase andrecirculation time. We showed that the commonly held assumptionof a 180◦ phase change is not the only solution for optimal energyrecovery.4. Lasing complicates beam transport in the return arc, and can resultin beam loss. In the TRIUMF ERL layout, this beam loss can bemitigated by appropriate arc optics.5. ERL operation is compatible with RIB operation. The ERL upgradedoes not affect the performance of the existing RIB photofission beam.6. High demand on optics in certain areas, such as conflicting demandson EHAT Q1.Some of the listed items are conceptually known processes, but the signif-icance of the optimization results is that they showed quantitatively howmuch these constraints and limitations matter, and whether they can becircumvented or avoided in the TRIUMF machine. The question of whethera transport solution exists for the given E-linac ERL layout was answeredand presented. The lattice layout proved sufficient and demonstrated thatthe optimization platform is a viable method of accelerator study and design.Additionally, the optimization platform was demonstrated to be a valu-able computational tool for physicists, both as an optimization tool andas a tool for studying global dynamic relationships in machines that wasnot possible previously. The ability to add or change modeling engines forstudying different processes contributes to its versatility and has never beenperformed before in the field of accelerator physics. No precedence for thistype of tool existed. Either genetic optimizations were performed in specificsections, such as injector optimization in the case of APISA, or an engineneeded to be created to encompass all the necessary physics. The TRIUMFplatform allows for global optimization with only the choice of the rightmodeling engines for the different physical processes needed to be included,and many engines exist in accelerator physics, each catering to its niche.The optimization platform was built with rigorous software engineeringstandards with the intentions of versatility and reusability. Indeed, appli-cation to other problems is already on the way. Collaboration with S. De-choudhury from the Variable Energy Cyclotron Centre (VECC) in Calcuttais in progress to optimize an injector linac in India. The problem requires97Chapter 7. Conclusionstwo modeling engines, one is the accelerator code TRACK, and the otheris a custom written Fortran engine. The custom written engine calculatesdrift lengths from initial cavity phases, and then passes the drift lengthsto TRACK for tracking. Neither engines were used for the TRIUMF ERLoptimization, but were easily implemented into the platform for the VECCproblem. This would have been impossible without the extensibility of theplatform or its multi-engine feature, and mitigated the need for VECC towrite a custom optimizer or tracking tool.The VECC problem was set up within two hours and put on a compu-tational cluster. The ease of transiting the platform into a parallel environ-ment is a great advantage, as large optimization problems are computation-ally intensive. The software design includes flexible and powerful exceptionhandling mechanisms, which makes the platform start-and-forget, allowingthe physicist to carry on with other tasks without the need for constantintervention and monitoring.Another work in progress is the optimization of a compression chicanefor the Deutsches Elektronen-Synchrotron facility in Germany. The projectmodels energy loss in a chicane for a 1 GeV electron beam, using the engineCSRtrack, with a second custom built engine that automatically generatesa bunch distribution given initial envelope parameters. A third vertex usingMADX runs in parallel to the first two, to calculate the transport mapelements of the chicane. This is another demonstration of the versatility ofthe optimization platform.The creation of the Empirical Model for ERL linac modeling is an in-cidental benefit. The design of the software makes it very useful for appli-cations outside of the optimization platform, where its fast running time issuitable for online applications. A port of Empirical Model was integratedinto XAL, a Java framework for creating high level applications for acceler-ators. XAL was already used for the first rounds of E-linac commissioningand proved to be invaluable. This demonstrates value of the optimizationproject outside of optimization.The optimization was designed as a beam dynamics study. Not includedare RF effects such as beam loading. Certain results, such as changing thearc path length, and therefore time-of-flight, for energy recovery, need to bechecked for compatibility with the RF system. When the acceleration anddeceleration phases are not 180◦ apart, less than optimal klystron operationfollows. In principle, RF loading can be added as another vertex in theoptimization. While optimal RF operations preferring 180◦ phase difference,a non-180◦ phase difference could be better for beam dynamics and energycompression. The results can be interesting and warrants further study.98Chapter 7. ConclusionsA limitation of the platform is its inability to perform discrete optimiza-tion. For example, what is the optimal number of dipole magnets in an arc?The produced baseline design can claim to be a good design with four-dipolearcs, but cannot claim to be unequivocally better than all designs with three-dipole arcs. Although the choice of four dipoles is intuitively better thanthree dipoles and the decision is physically sound, future questions may ariseregarding discreteness that are not as easily answered.Genetic algorithms can only perform optimization on continuous vari-ables. Therefore, answering questions regarding discreteness require eitheralgorithmic changes, or new features implemented into the platform. Onesuch feature that was discussed but not realized was the “prototype” feature.This allows the comparison of top level parameters such as final energy andemittances, while allowing for different arc layouts, or prototypes. We feelthis is the next stage in computational optimization in the field of acceleratorphysics.99Bibliography[1] Boost C++ Libraries. http://www.boost.org/. Accessed: 2015-09-04.[2] Git. https://git-scm.com/. Accessed: 2015-09-04.[3] LHC Machine Outreach.http://lhc-machine-outreach.web.cern.ch/lhc-machine-outreach/collisions.htm. Accessed: 2015-12-07.[4] Madx Quadrupole documentation.http://madx.web.cern.ch/madx/madX/doc/usrguide/Introduction/quadrupole.html. Accessed: 2015-10-31.[5] NumPy. https://www.numpy.org. Accessed: 2015-09-05.[6] Open XAL. http://xaldev.sourceforge.net/. Accessed: 2015-08-29.[7] The XML C Parser and Toolkit of Gnome libxml.http://www.xmlsoft.org/. Accessed: 2015-09-04.[8] WestGrid. https://www.westgrid.ca. Accessed: 2015-09-03.[9] Antokhin, E.A. et al. First lasing at the high-power free electron laserat Siberian center for photochemistry research. Nucl. Instrum. MethodsPhys. Res., Sect. A, 528(1–2):15–18, 2004.[10] Baartman, R. The Buckley Quads: Hysteresis, Calibration.http://lin12/text/ELinac/Quads/Buckley/2013Buckley.pdf, Dec2013. Accessed: 2015-12-10.[11] Bazarov, I. APISA. http://www.lepp.cornell.edu/~ib38/apisa/.Accessed: 2015-09-04.[12] Bazarov, I. bi – Beam Instability BBU Code.http://www.lepp.cornell.edu/~ib38/bbucode/src. Accessed:2015-10-30.100Bibliography[13] Bazarov, I. and Sinclair, C. Multivariate Optimization of a High Bright-ness DC Gun Photoinjector. Phys. Rev. ST Accel. Beams, Mar 2005.[14] Behre, C. et al. First lasing of the IR upgrade FEL at Jefferson lab.Nucl. Instrum. Methods Phys. Res., Sect. A, 528(1–2):19–22, 2004.[15] Ben-Zvi, I. and Kayran, D. and Litvinenko, V. High Average PowerOptical FEL Amplifiers. Proceedings of The 27th International FreeElectron Laser Conference, Aug 2005.[16] Benson, S. et al. High power lasing in the IR Upgrade FEL at JeffersonLab. Proceedings of the 2004 Free Electron Laser Conference, Aug 2004.[17] Benson, S. et al. High power operation of the JLab IR FEL driveraccelerator. Proceedings of the 2007 Particle Accelerator Conference,Jun 2007.[18] Bertsimas, D. and Tsitsiklis, J.N. Introduction to Linear Optimization.Athena Scientific, 1997.[19] Bleuler S. et al. PISA – A Platform and Programming Language In-dependent Interface for Search Algorithms. Proceedings of the 2003Evolutionary Multi-Criterion Optimization Conference, pages 494–508,Apr 2003.[20] Bogacz, A. et al. CEBAF Energy Recovery Experiment. Proceedingsof the 2003 Particle Accelerator Conference, May 2003.[21] Borland, M. Design and performance simulations of the bunch compres-sor for the Advanced Photon Source Low-Energy Undulator Test Linefree-electron laser. Phys. Rev. ST Accel. Beams, 4(7):074201, 2001.[22] Brau, C.A. The Vanderbilt university free-electron laser center. Nucl.Instrum. Methods Phys. Res., Sect. A, 319(1–3):47–50, 1992.[23] Braun, H. H. et al. Emittance growth and energy loss due to coherentsynchrotron radiation in a bunch compressor. Phys. Rev. ST Accel.Beams, 3(12):124402, 2000.[24] Buts, V.A. and Lebedev, A.N. and Kurilko, V.I. The Theory of Coher-ent Radiation by Intense Electron Beams. Springer, 2006.[25] Carey, D. C. The Optics of Charged Particle Beams. Harwood, 1987.101Bibliography[26] Carey, D.C., Brown, K.L. and Rothacker F. Third-Order TRANSPORTwith MAD Input: A Computer Program for Designing Charged ParticleBeam Transport Systems. Technical report SLAC-R-530, Fermilab-Pub-98-310, UC-414, Oct 1998.[27] CERN. MAD – Methodical Accelerator Design.http://madx.web.cern.ch/madx, 2015. Accessed: 2015-07-06.[28] Chao, Y.C. Recipe for Empirical Model Propagation.http://trshare.triumf.ca/~chao/Chris/HLA/Recipe for Empirical Model Propagation.pdf, Jun 2011. Ac-cessed: 2015-12-10.[29] Chao, Y.C. Baseline Pass 1 ERL Design.http://trshare.triumf.ca/~chao/Chris/Recirculation/7MeV_New/MergeBaselinePass1ERLRecFit3.opt, 2012. Accessed:2015-10-20.[30] Chao, Y.C. Baseline Pass 2 ERL Design.http://trshare.triumf.ca/~chao/Chris/Recirculation/7MeV_New/MergeBaselinePass2ERLDmpFit3.opt, 2012. Accessed:2015-10-20.[31] Chao, Y.C. Baseline RIB Extraction Design.http://trshare.triumf.ca/~chao/Chris/Recirculation/7MeV_New/MergeBaselinePass1ERLExtFit3.opt, 2012. Accessed:2015-10-20.[32] Chao, Y.C. Empirical Model Interpolation Tables.http://trshare.triumf.ca/~chao/InterpolationTables/, 2013.Accessed: 2015-12-10.[33] Chao, Y.C. E Linac EMBT-EABT-EHAT Phase One Major Compo-nents and Layout. TRIUMF design note TRI-DN-12-03, Mar 2014.[34] Chao, Y.C. Wiggler Parameters and Matching Conditions.http://trshare.triumf.ca/~chao/Chris/FEL/Wiggler Parameters and Matching Conditions.pdf, Jun 2014. Ac-cessed: 2015-12-10.[35] Chao, Y.C. et al. Low-Beta Empirical Models used in Online Modelingand High Level Applications. Proceedings of the 2011 InternationalParticle Accelerator Conference, Sep 2011.102Bibliography[36] Clarke, J. A. The Science and Technology of Undulators and Wigglers.Oxford University Press, Sep 2004.[37] Dechoudhury, S. Beam dynamics design of 160 MHz Linac in Asymet-rical Alternate Phase Focusing (A-APF) mode. Unpublished note, Jul2015.[38] Di Mitri, S. and Venturini, M. CSR induced emittance growthand related design strategies. USPAS Accelerator Physics Coursehttp://uspas.fnal.gov/materials/15Rutgers/lecture_We9.pdf,2015. Accessed: 2015-12-10.[39] Ding, Y.T. et al. Study on the planar undulator scheme with focusingproperties for PKU-FEL. Proceedings of the 2004 Free Electron LaserConference, Aug 2004.[40] Dohlus, M. and Limberg, T. CSRtrack.http://www.desy.de/xfel-beam/csrtrack/. Accessed: 2015-10-30.[41] Dohlus, M and Limberg, T. CSRtrack: faster calculation of 3D CSReffects. Proceedings of the 2004 Free Electron Laser Conference, Aug2004.[42] Elias, L.R. and Hu, J. and Ramian, G. The UCSB electrostatic acceler-ator free electron laser: First operation. Nucl. Instrum. Methods Phys.Res., Sect. A, 237(1):203–206, 1985.[43] Elias, L.R. and Jaccarino, V. and Yen, W.M. Scientific research withthe UCSB free electron laser. Nucl. Instrum. Methods Phys. Res., Sect.A, 239(3):439–442, 1985.[44] Feldman, D.W. et al. Energy recovery in the Los Alamos free electronlaser. Nucl. Instrum. Methods Phys. Res., Sect. A, 259(1–2):26–30,1987.[45] Floettmann, K. A Space Charge Tracking Algorithm.http://www.desy.de/~mpyflo/. Accessed: 2015-08-29.[46] Garoby, R. Beam Loading in RF Cavities. In Frontiers of ParticleBeams: Intensity Limitations, pages 509–541. Springer, 1992.[47] Gong, C. CSR Optimization in Chicane. TRIUMF unpublished note,2013.103Bibliography[48] Gong, C. Empirical Model Design Document. TRIUMF unpublishednote, 2013.[49] Gong, C. and Chao, Y.C. A New Platform for Global Optimization.Proceedings of the 2012 International Particle Accelerator Conference,May 2012.[50] Gong, C. and Chao, Y.C. The TRIUMF Optimization Platform andApplication to the E-linac Injector. Proceedings of the 2012 Interna-tional Computational Accelerator Physics Conference, Aug 2012.[51] Hajima, R. Energy Recovery Linacs for Light Sources. Reviews ofAccelerator Science and Technology, 03(01):121–146, 2010.[52] Hajima, R. et al. Design of energy-recovery transport for the JAERIFEL driven by a superconducting linac. Nucl. Instrum. Methods Phys.Res., Sect. A, 445(1–3):384–388, 2000.[53] Hajima, R. et al. First demonstration of energy-recovery operation inthe JAERI superconducting linac for a high-power free-electron laser.Nucl. Instrum. Methods Phys. Res., Sect. A, 507(1–2):115–119, 2003.[54] Hassan, R. et al. A Comparison of Particle Swarm Optimization andthe Genetic Algorithm. Proceedings of the 46th Structures, StructuralDynamics, and Material Conference, Apr 2005.[55] Herman, W. Fourth Generation Light Sources. Proceedings of the 1997Particle Accelerator Conference, May 1997.[56] Hoffstaetter, G. H. and Bazarov, I. V. Beam-Breakup Instability The-ory for Energy Recovery Linacs. Phys. Rev. ST Accel. Beams, 7, May2004.[57] Holder, D. J. et al. The Status of the Daresbury Energy RecoveryLinac Prototype. Proceedings of the 2008 European Particle AcceleratorConference, Jun 2008.[58] Jones, K. Comparison of Genetic Algorithm and Particle Swarm Op-timisation. Proceedings of the 2005 International Conference on Com-puter Systems and Technologies – CompSysTech2005, 2005.[59] Kim, K.J. and Huang, Z.R. and Lindberg, R. Introduction to thePhysics of Free Electron Lasers. Jun 2010. Unpublished manuscript.104Bibliography[60] Kolb, P. HOM Cavity Design for the TRIUMF Elinac. Proceedings of2011 International Conference on RF Superconductivity, Jul 2011.[61] Kolb, P. HOM Considerations for the TRIUMF Elinac. Presented atthe HOM Workshop 2012, Jun 2012.[62] Kolb, P. HOM Data. Private communications, 2015.[63] Koscielniak, S. ARIEL and E-linac Conceptual Design Report.http://lin12.triumf.ca/text/ELinac/CDR.ps, 2008. Accessed:2015-12-10.[64] Laxdal, R. Private communications, Feb 2015.[65] Liepe, M. et al. Pushing the Limits: RF Field Control at High LoadedQ. Proceedings of the 2005 Particle Accelerator Conference, May 2005.[66] Liu, C.Y. and Krafft, G. and Wang, G.M. Performance Evaluation ofUndulator Radiation at CEBAF. Proceedings of the 2010 InternationalParticle Accelerator Conference, May 2010.[67] Lu, H. H. et al. Research on the undulator used for PKU-FEL. Pro-ceedings of the 2004 Free Electron Laser Conference, Aug 2004.[68] Lukes, Z. Multi-Objective Optimization of Wire Antennas: Genetic Al-gorithms Versus Particle Swarm Optimization. Radioengineering Jour-nal, Dec 2005.[69] Madey, J.M.J. Stimulated Emission of Bremsstrahlung in a PeriodicMagnetic Field. Journal of Applied Physics, 42(5):1906–1913, 1971.[70] Marchetto, M. et al. Beam Dynamics Optimization of the TRIUMFElinac Injector. Proceedings of the 2009 Particle Accelerator Confer-ence, May 2009.[71] Marcus, G. Private communications, Jun 2015.[72] Merminga, L. Energy Recovering Linac Issues. Presented at the 2002Electron Ion Collider Accelerator Workshop, Feb 2002.[73] Merminga, L. Energy Recovery Linacs. Proceedings of the 2007 ParticleAccelerator Conference, Jun 2007.[74] Minehara, E.J. et al. A 0.1 kW operation of the JAERI superconduct-ing RF linac-based FEL. Nucl. Instrum. Methods Phys. Res., Sect. A,429(1–3):9–11, 1999.105Bibliography[75] Muratori, B.D. et al. Space Charge Effects for the ERL PrototypeInjector Line at Daresbury Laboratory. Proceedings of the 2005 ParticleAccelerator Conference, May 2005.[76] Neil, G. R. et al. Sustained Kilowatt Lasing in a Free-Electron Laserwith Same-Cell Energy Recovery. Phys. Rev. Lett., 84:662–665, Jan2000.[77] Neil, G. R. et al. The Jefferson Lab Free Electron Laser Program.International Symposium on Infrared Free Electron Laser and its Ap-plication, Tokyo, Jan 2002. JLAB technical document jlab-acc-02-02.[78] Onuki, H. and Elleaume, P. Undulators, Wigglers and Their Applica-tions. CRC Press, 2003.[79] Piot, P. and Douglas, D.R. and Krafft, G.A. Longitudinal Phase SpaceManipulation in Energy Recovering Linac-Driven Free-Electron Lasers.Phys. Rev. ST Accel. Beams, Oct 2003.[80] Planche, T. ARIEL S34 Dipoles design note.http://beamphys.triumf.ca/~tplanche/designs/elinac_dipoles/design_notes/Strong34/stong34.pdf, Apr 2013.Accessed: 2015-12-10.[81] Planche, T. ARIEL ’Y30’ Dipoles design note.http://beamphys.triumf.ca/~tplanche/designs/elinac_dipoles/design_notes/Y30/Y30.pdf, Jul 2013. Accessed:2015-12-10.[82] Powers, T. and Tennant, C. Implications of Incomplete Energy Re-covery in SRF-Based Energy Recovery Linacs. Proceedings of the 41stAdvanced ICFA Beam Dynamics Workshop on Energy Recovery Linacs,May 2007.[83] Rand, R. E. Recirculating Electron Accelerators. Harwood, 1984.[84] Reiche, S. GENESIS 1.3. http://genesis.web.psi.ch/. Accessed:2015-12-04.[85] Rohatgi, R. et al. The SCA/FEL program: Operation in the infrared,visible and ultraviolet. Nucl. Instrum. Methods Phys. Res., Sect. A,272(1):32–36, 1988.106Bibliography[86] Satogata, T. More Lattice Optics and Insertions. USPAS AcceleratorPhysics Course http://www.toddsatogata.net/2011-USPAS/2011-06-17-Insertions.pdf, Jun 2011. Accessed: 2015-12-09.[87] Schriber, S.O. and Heighway, E.A. Double Pass Linear Accelerator –Reflexotron. IEEE Transactions on Nuclear Science, 22(3):1060–1064,Jun 1975.[88] Servranckx, R. et al. Users Guide to the Program DIMAD. Technicalreport SLAC-R-285, Jan 2004.[89] Smith, S.L. The Status of the Daresbury Energy Recovery Linac Proto-type (ERLP). Proceedings of the 41st Advanced ICFA Beam DynamicsWorkshop on Energy Recovery Linacs, May 2007.[90] Smith, T.I. et al. Development of the SCA/FEL for use in biomedicaland materials science experiments. Nucl. Instrum. Methods Phys. Res.,Sect. A, 259(1):1–7, 1987.[91] Springob, N. VariantCalc.http://boost-spirit.com/repository/applications/show_contents.php. Accessed: 2015-09-05.[92] Thompson, N. Introduction to Free-Electron Lasers.http://www.astec.stfc.ac.uk/ASTeC/Resources/PDF/Thompson_FELs.pdf, 1997. Accessed: 2015-12-10.[93] Tigner, M. A Possible Apparatus for Electron Clashing-Beam Experi-ments. Il Nuovo Cimento, 37:1228–1231, Jun 1965.[94] Venturini, M. CSR-induced emittance growth in achromats: Linear for-malism revisited. Nucl. Instrum. Methods Phys. Res., Sect. A, 794:109–112, May 2015.[95] Vu, V. T. A Comparison of Particle Swarm Optimization and Differ-ential Evolution. International Journal on Soft Computing, Aug 2012.[96] Wang, G.M. et al. Energy recovery transport design for PKU FEL.Proceedings of the 2007 Particle Accelerator Conference, Jun 2007.[97] Yang, L.Y. and Li, Y.J. and Guo, W.M. and Krinsky, S. MultiobjectiveOptimization of Dynamic Aperture. Phys. Rev. ST Accel. Beams, May2011.107[98] Zitzler, E. and Laumanns, M. and Thiele, L. SPEA2: Improving theStrength Pareto Evolutionary Algorithm.http://e-collection.library.ethz.ch/eserv/eth:24689/eth-24689-01.pdf, May 2001. Accessed: 2015-12-09.[99] Zitzler, E. et al. Comparison of Multiobjective Evolutionary Algo-rithms: Empirical Results. Evol. Comput., 8(2):173–195, Jun 2000.108Appendix ANotationsAccelerator physics notations used in this document conform to TRANS-PORT notations [26].A particle is defined by the six coordinates x, x′, y, y′, z, and δ, wherex′ = px/p is the x-angle defined by the ratios of electron momentum to bunchcentroid momentum, y′ = py/p is similarly the y-angle, and δ = (pz − p)/pis the energy deviation of the particle to the centroid. z = −ct is the timedifference compared to the time of the bunch centroid, scaled by the negativespeed of light −c. Positive z denotes that the particle arrives at a locationbefore the centroid, in the centroid frame. The six coordinates can alsobe labeled from 1 to 6, or xi where i = 1, ..., 6. The RMS values of thecoordinates are defined by σ, for example, σx ≡√〈x〉.The particle energy is represented by E = γmc2, where γ is the particle’sLorentz factor and m is the electron rest mass. The bunch centroid energyis represented by Er and correspondingly γr.RMS bunch Courant-Snyder (CS) parameters are defined in the usualmanner:εx,rms =√〈x2〉〈x′2〉 − 〈xx′〉2βx,rms = 〈x2〉/εx,rmsαx,rms = −〈xx′〉/εx,rms(A.1)where 〈〉 denotes the mean, ε is the emittance, β is related to the beam size,and α is related to the beam tilt. Similar equations follow for y and z.The transfer map of a beamline lattice follows typical accelerator con-vention. The second order transport is given by the Taylor expansionxi =∑jMijx0j +∑j,kTijkx0jx0k (A.2)where x0i denotes the coordinates before the map is applied, and xi after.The first order map elements are given by Mij = ∂xi/∂x0j , and the secondorder elements are given by Tijk = (1/2)∂2xi/∂x0j∂x0k. For example, themap element M12 denotes how a beamline section affects the x coordinategiven the x′ coordinate of the particle at the start of the section.109Appendix A. NotationsDispersion ηx is defined as, for an off-momentum particle with energydeviation δ, ηxδ is the x deviation of the off-momentum closed orbit fromthe reference orbit. Similarly for dispersion prime η′x, η′xδ is the x′ deviation.Identical definitions exist for y and y′.RF phases are defined using φij , where i is the linac pass number, and jis the cavity number. i can range from 1 to 3, with 1 and 2 representing theacceleration and deceleration passes of the ERL beam, and 3 representingthe acceleration pass of the rare isotope beam. j can range from 1 to 4,representing the four 9-cell cavities of the linac. If only one subscript isdisplayed, e.g., φ2, the subscript represents the index i, with the cavitynumber j assumed to be 1.The phase φij is the phase (in degrees◦) of the jth cavity when the bunchcentroid is at the jth cavity entrance during pass i. We choose to use thecavity entrance phase so as to conform to the Empirical Model convention(see chapter 3). Sometimes, φij0 is used, denoting the initial phase of thejth cavity at the beginning of machine modeling, for pass i.FEL literature [59] typically replaces the particle coordinate z with theponderomotive phase θ = (k + ku)r − wτ + const, where k and w are thewavenumber and angular frequency of the FEL radiation, ku is the wavenum-ber of the undulator periodicity, r is the position in the undulator, and τ isthe time an electron arrives at r. The other five coordinates remain identical.Quadrupole strengths are given in units of K1 in units of [m−2], definedby [4, 25]K1 =1Bρ∂By∂x(A.3)where Bρ = p/q is the beam rigidity and q is the electric charge, and ∂By/∂xis the quadrupole gradient.110Appendix BSoftware Design of theOptimization PlatformIntroductionThe goal of the TRIUMF optimization project is to study the TRIUMF ERLdesign, although the input format is flexible enough to be easily extended toother optimization problems. The code executes on parallel machines andcan easily switch between different physics engines.The TRIUMF optimization framework is based upon APISA [11], whichis a realization of PISA [19], an interface designed for the global optimiza-tion of multiple objectives based on pseudo-random sampling algorithms.The genetic algorithm SPEA2 [98] is chosen, which has precedence of priorapplications to accelerator physics [13, 97]. A genetic algorithm is a popula-tion based algorithm which solves the optimization problem using analoguesto processes from evolutionary biology:1. Generate random starting population, e.g., a set of ERL designs eachwith different design parameters.2. Evaluate fitness of each individual, or ERL design, in the population.An ERL design which better satisfies the objectives and constraints ismore fit.3. Stochastically choose individuals (weighed by fitness), mutating or re-combining them to form new individuals.4. Throw away some individuals (designs with worst fitness more likelyto be thrown away). Add new individuals added to the population.Repeat from step 2.SPEA2 uses two processes for creating new individuals: mutation, wherea copy of a parent is made, then each parameter of the child copy has achance to change, i.e. mutate. The second is binary crossover where twoparents are chosen, and each parameter of the child is randomly taken from111Program Descriptionone parent or the other. The specifics of the algorithm can be found in[98]. The algorithm ranks fitness based on Pareto dominance. An exampleis shown later in this section.The algorithm iterates until the terminating condition is satisfied, i.e.maximum number of generations (iterations) is reached.PISA conceptually separates the software into two components, Variatorand Selector. Each program is compiled to its own binary and executedseparately. Data is exchanged between the two through text files.Program DescriptionThe optimization platform is written in C++ in a 64-bit Linux environment.RequirementsThe TRIUMF optimization platform must have the following:1. Allow multiple engines to be included into the same optimization prob-lem. This allows for a global setup and discovers dynamic relationshipstypically not possible with local optimization.2. Parallel capable - moving from local to global optimization necessitateshigher computing resources.3. Easily extensible to add new engines.4. A generic input file capable of describing all optimization problemsdescribed by section B.5. Good exception handling mechanism in preparation for the multitudeof errors that can occur with the inclusion of different simulation en-gines.6. Handle unit conversions between different engines.7. Create python codebase for common post-processing operations, suchas extraction of energy and emittances.PrerequisitesThe input file for the optimization program is defined in XML format. Thelibxml2- [7] XML parser was used. All C++ projects using libxmlmust include the flags -lxml2 -lm in the linker options.112Program DescriptionThe Boost libraries [1] are required, in particular boost filesystem,boost spirit and boost regex.The following environmental variables are required:ˆ LD LIBRARY PATH - directory containing the compiled optimization li-braries.ˆ PYTHONPATH - directory containing the python post-processing code.This path is optional if the included python code is not used.Configuration ManagementThe development environment gongc.triumf.ca uses Git [2] for versioncontrol.All design documentation are also version controlled.DeliverablesIn build order:1. libDiagnostics.so2. libAccessibility.so3. libEvaluator.so4. spea2 (executable)5. Variator (executable)6. Minion (executable)7. OptimizationMain (executable)The executable OptimizationMain is the entry point for running opti-mization.Top Level LogicThe three major components of the optimization program are (Fig. B.1):ˆ Variator - the optimization base. Controls mutation/crossbreeding,giving information to Selector on the current generation. Passes infor-mation to Evaluator on what runs to make.ˆ Selector - part of the PISA framework that evaluates fitness of indi-viduals. The SPEA2 algorithm used here is almost unchanged fromthe original.113Program Descriptionˆ Evaluator - wrapper framework for the physics engines. Gets informa-tion from Variator on what to run, runs the engines, then sends theinformation back to Variator. Evaluator is completely new and doesnot appear in PISA or APISA.Evaluator (evaluates population) Variator (population control) Selector (algorithm) <communicates via text files> Minion (evaluates individual) Figure B.1: Optimization top level logic. As per the PISA principle, theplatform is divided into Variator (population control) and Selector (popu-lation evolution algorithm). Due to the large amount of code required forthe Variator side, two new modules were created to divide the work. TheVariator base is a skeleton state machine. Evaluation of the entire popula-tion is performed by Evaluator, which then evaluates each individual in thepopulation via Minion. Minion communicates with Evaluator via text filesbecause due to the distributed nature of the platform, they may not executeon the same node. Selector and Variator also communicates via text filesfor decoupling.Variator and Selector are not coupled to each other programmatically,but rather run as separate programs that communicate by reading/writingtext files. They run as state machines with the state codes passed via thetext files. Evaluator is a part of Variator, and thus could be called directly.We make a distinction between global and local problems. A globalproblem is the overall ERL optimization. Given a population of individu-als (i.e. ERL designs), we select the best individuals in the hope to max-imize/minimize some parameters. Each optimization problem is a singleglobal problem.114Program DescriptionFigure B.2: Layers of the optimization problem. The top layer is the globalproblem, in which the optimum individual(s) are selected from a population.The lower layer is how each section of the ERL is produced.Figure B.3: Example topology. Simulation engines can be linked togetherin parallel or serial for flexibility.A local problem is the production of a single individual. Figure B.2 illus-trates the relationship between global and local problems. Each individualis an instance of the model (in this case, an instance of the ERL design). Wedo this through multiple local problems; the output of each local problem ispiped to the next local problem to form a continuous and complete simula-tion. It is convenient to represent the local problem as a graph (Fig. B.3).Each process is a call to a simulation program, such ASTRA, PARMELA, or115Program Descriptiona custom written program. Evaluator links the output of one process to theinput of the next. The graph object in Fig. B.3 is referred to as a topology,and each process is a vertex. Local problems do not always have to repre-sent a physical section of a beamline. Example 1: GPT to GPT2ASTRAto ASTRA. Here the local problem GPT2ASTRA is a conversion program,which takes the output of GPT and massages it to a format suitable for thesubsequent ASTRA run.Selector is purely concerned with the global problem, Variator with bothglobal and local, and Evaluator purely with local. The top level architecturaldesign is shown in Fig. B.4.116ProgramDescriptionFigure B.4: Optimization platform class diagram.117Program DescriptionVariatorVariator is written entirely in C and C++. The initial Variator codebaseis Ivan Bazarov’s APISA [11] framework, which in turn is based on PISA.APISA was written entirely in C, with the Variator section consisting of allglobal variables and global functions. Additionally, APISA works only withASTRA and runs locally. The upgrade from APISA to the TRIUMF PISAframework consisted of the following stages:1. In the near term, encapsulate the APISA Variator in a namespace“Variator”. The code should be usable at this point, but does not fitperfectly within the object-oriented framework.2. In the long term, modify the APISA Variator to be object-oriented.The reason is the exponential increase in the scale of the software.Moving to a multithreading, multiple engine setup significantly in-creases the complexity of the software, which is why the modularity ofthe OO-framework is desired. This is a difficult and time-consuming.Good OO-design requires that each class serves an explicit purpose.The global functions of the APISA Variator, which are split into mul-tiple files, does not fit this requirement.Variator acts as an interface in the creation of new generations and newindividuals. It does not do the real work for either processes, but directsthe flow. Variator communicates with Selector to determine fitness, andcommunicates with Evaluator to produce offspring. It is a state machine(see Fig. B.5).118Program DescriptionFigure B.5: Variator state chart.SelectorSelector is the component that evaluates the fitness of each individual, andthen determines the individuals of the population that are selected as parentsfor the next generation. The fitness evaluation algorithm is SPEA2 [98].The Selector state chart is shown in Fig. B.6. APISA’s implementationof SPEA2 is used as the codebase. Refactoring was performed on the code-base. Otherwise, algorithmic changes were not done. All SPEA2 settingsare moved to the optimization input XML file. This makes the XML file aunified location for all optimization input settings.119Program DescriptionFigure B.6: Selector state chart.EvaluatorEvaluator controls the physics engines. For each individual of the popu-lation, Variator calls Evaluator to run the engines, after which numericalvalues such as final emittances are produced. Each individual is obtained120Program Descriptionvia a sequence of engine processes. Evaluator does not exist in APISA, sowill be written from scratch following rigorous object-oriented design princi-ples. Evaluator also handles the multithreading using the POSIX threadinglibrary.121ProgramDescriptionFigure B.7: Evaluator execution flowchart.122Program DescriptionFigure B.8: Evaluator multithreading flowchart.Fig. B.7 shows the process of evaluating each individual, and the inter-actions between Evaluator and the external files and physics engines. TheEvaluator state diagram (Fig. B.8) shows the algorithm Evaluator usesto handle parallel computing. To take advantage of parallel-processing, athread is created for each vertex to be executed. Each thread uses ssh toaccess a node, and executes the engine instructions. The number of threadsshould never exceed the number of available nodes.123Program DescriptionFigure B.9: Evaluator class diagrams.Evaluator consists of two pieces. The backbone is written in C++ andacts as a mediator between the frontend and the Variator, and as a wrapperfor the physics engine. The frontend allows user-injected python code, usefulfor post-processing. The class design of Evaluator is shown in B.9. Evaluatorcontains three major classes and interfaces: EvaluatorMain, INodeManager,and IEngineManager.ˆ EvaluatorMain - this singleton class is the interface between Evaluatorand Variator. Variator passes data of individuals to this class to beevaluated. For the sake of modularity, this is the only entry pointavailable to Evaluator.ˆ INodeManager - interface for classes that handles computing nodes.EvaluatorMain creates a NodeManager depending on the environment,e.g., WestgridNodeManager if the optimization is run on WestGrid.The NodeManager determines which vertex is assigned to which node.Reasons for different NodeManager for different environments include– Different environments have different ways of retrieving the listof nodes available to do work– Different filesystem: Westgrid nodes all share the same filesystem,whereas a cluster of individuals do not. This introduces limits onwhich vertex can be executed on which node. A vertex whichrequires output files of previous vertices must be executed on the124Program Descriptionsame filesystem as those previous vertices, to avoid copying filesbetween systems. See below for more information.The NodeManager also handles thread creation, monitoring, and ter-mination.ˆ IEngineManager: this is the interface between Evaluator, which treatsnodes and engines as abstract objects, and the actual nodes and en-gines which executes the vertices. The EngineManager logs into theworking node by ssh and runs the simulation engine.Assigning Jobs to a NodeWe distinguish between vertex-based assignment and individual-based as-signment:ˆ Vertex-based - used for both local and Westgrid executions. Verticesbelonging to the same individual can be executed by different nodes.This is designed for systems where nodes share the same filesystem.ˆ Individual-based - designed for systems where nodes do not share thesame filesystem, for instance, a network of different computers. Allvertices of the same individual must be executed on the same node.Since vertices might exchange information with each other, this avoidshaving to copy files and data between nodes.125ProgramDescriptionFigure B.10: Population evaluation statechart.126Program DescriptionThe optimization platform assigns each individual to a node for eval-uation. A node represents a single processor, which may contain multiplecores. Multiple threads are created, with each thread evaluating a vertex ofthe individual’s topology. The relationship between the number of threadsand the number of cores in the node can be many-to-one. The evaluatealgorithm is shown in Fig. B.10.Figure B.11: Work distribution to computational clusters. The platformcan easily take advantage of high performance computing centers such asWestGrid. An important feature of the platform is scalability: it can be runon both large clusters and also individual local machines.An important feature of the platform is scalability. It can be run onboth large clusters and also individual local machines. This is significantbecause not all optimization problems require high performance computing.Any optimization problem setup, including the platform binaries and theproblem description files, can be switched from WestGrid to a local machineand run without recompiling any portion of the code. The work distributionmodel is shown in Fig. B.11. On a local machine, all nodes simply refer tolocalhost.127Program DescriptionProcess MonitoringOne of the biggest additions to the TRIUMF platform is process oversight.The increase in software scope also resulted in an increase in exceptionproduction. Process monitoring refers to exception handling in the multi-threading, parallel computing scheme. When the main node assigns a jobto a worker node via ssh, the worker process needs to be monitored, i.e.we need to know when the worker node finishes the job, so we can analyzethe results and free the worker node for additional jobs. There are twomethods of process monitoring used in the optimization platform: processvs processless.1. Process-based monitoring - requires Environment=WESTGRID and Pro-cessLess=NO, or Environment=LOCAL. Each thread creates a monitorprocess. The steps are shown in Fig. B.12:(a) Create monitor on head node.(b) Monitor assigns job to worker via ssh.(c) Monitor does not return until job finishes.The advantage of such a system is the monitor takes care of trackingevery aspect of the job. The disadvantage is for large optimizations,there are many monitor processes on the head node, resulting in alarge memory overhead.2. Processless job monitoring - requires Environment=WESTGRID and Pro-cessLess=YES. The steps are shown in Fig. B.13:(a) Create monitor on head node.(b) Monitor assigns job to worker via background ssh, which returnsthe PID. Monitor returns immediately.(c) Thread polls PID on worker node to check if engine finished run-ning.This is the preferred way of running because it does not overload thehead node with monitoring threads. More details can be found belowB.128Program DescriptionMain thread Thread 1 Thread 2 Thread 3 Monitor 1 Monitor 2 Monitor 3 <creates POSIX threads> <fork, exec> Workers (many nodes) Engine 1 Engine 2 Engine 3 <ssh> Wait for engine to finish, then ssh returns Figure B.12: Process monitoring with monitoring thread. This is the idealmethod for small optimization problems which do not require many nodes.Threading is the natural method to monitor worker nodes and check theirstatus.129Program DescriptionMain thread Thread 1 Thread 2 Thread 3 Monitor 1 Monitor 2 Monitor 3 <fork, exec> Workers (many nodes) Astra 1 Astra 2 Astra 3 <ssh> Engine started via background ssh. Monitor process returns immediately Thread polls PID of engine to check for execution finish <creates POSIX threads> Figure B.13: Process monitoring without monitoring thread. This is idealfor large optimization problems which require many nodes.Minion: Local Population EvaluatorFor processless monitoring 2, the head node (where Variator runs) does notsubmit individuals to nodes to be run. Instead, the run is outsourced toMinion. Variator starts Minion on a given node, and passes individuals tobe evaluated. Minion then runs in process monitoring mode and evaluatesthe individuals on its local host node. Variator and Minions exchange datathrough text files. The interaction flowchart is shown in Fig. B.14.The purpose of such a scheme is that the POSIX threading system, whilepowerful, is not designed for extreme numbers of threads. Large optimiza-tion problems often make use of dozens of worker nodes. If one thread isassigned to each node, terrible performance is caused by the large numberof threads due to virtual memory limitations, context-switching overhead,and scheduler overhead. Globally, Variator circumvents this threading prob-lem by executing Minion as a detached process. Variator then uses a singlethread, the main thread, to monitor the status of all Minions. Locally,Minion is allowed to assign a small number of threads to take advantage ofparallel computing on a multicore node, but not too many to cause virtualmemory problems. This combination of multiprocessing and multithreadingis an efficient way of executing a genetic algorithm.To allow Minion to operate correctly, set the XML parameter Use-130Program DescriptionMinion=YES. Minion works as a detached process, so it is necessary for theexecuting environment it runs in (which is different from the environmentthat Variator is running in) to contain the correct variables. Set the variablesLD LIBRARY PATH and PYTHONPATH in the .bashrc file.Figure B.14: Minion flowchart.Custom Python CodeCustom code allows the user to write problem-specific code in a conve-nient manner. In the template directory of each local problem is a file,customcode.py, with the function RunCustomCode. A list of running vari-ables is passed to this function, which the user can use to calculate newvariables and write them to file. The user can also read from program out-put files to obtain values to assist in the calculation of new quantities.Many typical functionalities for engines exist in included python files,e.g., astra.py, madx.py, to calculate and extract common accelerator physics131Input and Output Formatsvariables such as emittances, Courant-Snyder functions, and transport mapelements. The location of these files should be stated in a global shellvariable PYTHONPATH. By default, the path is $INSTALLPATH/python, where$INSTALLPATH is the location of the optimization engine binaries.Input and Output FormatsThe TRIUMF optimization platform uses two types of input files:1. Optimization input (global) - input for the optimization program2. Physics engine input (local) - each modeling engine has its own inputfiles and format, e.g., MADX uses input files based on TRANSPORTnotationThe global input format is presented first, followed by the local inputformat.Global Input File Problem.xmlGlobal input is defined in the XML file Problem.xml. This file includes1. Problem description - parameters, constraints, objectives, and topol-ogy2. Settings - variables related to the execution of the optimization pro-gram, such as output folder path and seedThe XML file problem.xml includes everything needed to define the cur-rent optimization problem. The XML file must begin with the tag <Pro-blem>. The following blocks must be included inside the problem tag (inorder):1. SettingsList2. UnitsList3. DefaultUnits4. AutoNDParamLists5. ParamList6. NDParamList7. ConstraintList8. ObjectiveList9. Topology132Input and Output FormatsSettings InputThe <Settings> tag contains attributes related to the run:<SettingsList><Setting name="RunName" value="sample_run" /><Setting name="BaseDir" value="./" /><Setting name="ExchangeDir" value="exch" /><Setting name="TemplateDir" value="template" /><Setting name="OutputDir" value="output" /><Setting name="OS" value="UNIX" /><Setting name="Environment" value="WESTGRID" />...</SettingsList>All name and value pairs are case-sensitive.Units and Default Units InputThe <UnitsList> tag contains the list of available units. Information onunits is provided in section B.Template for Automatic Lists InputOften we want to retrieve the same set of variables at many points on thebeamline. For example, we might want x, px, y, py at the points ’start’,’middle’, and ’end’ of the lattice. We could create all the variables like such<NDParamList><NDParam name="x_start" ... /><NDParam name="px_start" ... /><NDParam name="y_start" ... /><NDParam name="py_start" ... /><NDParam name="x_middle" ... /><NDParam name="px_middle" ... /><NDParam name="y_middle" ... /><NDParam name="py_middle" ... /><NDParam name="x_end" ... /><NDParam name="px_end" ... /><NDParam name="y_end" ... /><NDParam name="py_end" ... />...</NDParamList>133Input and Output FormatsThe alternative is to add a block after UnitsList:<AutoNDParamLists><AutoListDefinition name="beam"><Item name="x" unit="m" /><Item name="px unit="MeV/c" /><Item name="y" unit="m" /><Item name="py" unit="MeV/c" /></AutoListDefinition></AutoNDParamLists>Later in the NDParamList section:<NDParamList><NDParam ... /><AutoList listname="beam" label="_start" /><AutoList listname="beam" label="_middle" /><AutoList listname="beam" label="_end" /><NDParam ... /></NDParamList>Optimization will automatically expand the AutoList block to create all thevariables necessary. Auto lists are used only for non-decision parameters.For decision parameters, the upper and lower values must be specified foreach, therefore they cannot be automated.Decision Parameters InputDecision parameters x (see Eqn. 1.1) are defined under the <ParamList>tag. These are the free parameters that the platform will try to optimize.Boundaries and units can be specified for each parameter.<ParamList><Param name="Energy" minvalue="40" maxvalue="50"unit="MeV" /><Param name="Q1" minvalue="0" maxvalue="2"unit="T" /><Param name="B1Length" minvalue="5" maxvalue="10"unit="cm" /></ParamList>The value of the unit attribute must be predefined in the <Units> block.134Input and Output FormatsNon-Decision Parameters InputNon-decision parameters y(x) are defined under the <NDParamList> tag.These are the implicit parameters that are not free parameters. They canbe (but do not have to be) used as constraints and objectives.<NDParamList><NDParam name="energyloss" unit="MeV" /><NDParam name="K" unit="dimensionless" /></NDParamList>In the above example, energy loss is typically not a free parameter pre-determined at the start of beamline design, but rather a result of latticesettings.Constraints InputConstraints are Fj(x) defined under equation 1.3. They are represented bythe <ConstraintList> tag. Example:<ConstraintList><Constraint param="energyloss"direction="LT" value="1.2" unit="MeV" /></ConstraintList>The direction attribute can be less than (LT) or greater than (GT).Objectives InputObjectives are gi(x) defined under equation 1.2. They are implemented bythe <Objective> block as follows:<ObjectiveList><Objective param="emitx_end" direction="minimize" /></ObjectiveList>The direction attribute can be minimize, maximize, or equals. If thedirection is equals, the Objective element must also contain an additionalvalue attribute, for example:<Objective param="x"direction="equals" value="0.5*{mm}" />135Input and Output FormatsThe above is equivalent tomin |x− 0.5 mm| (B.1)representing that we want to optimize a parameter such that it is as closeto a given value as possible. The value attribute is not used and can beomitted if the direction is maximize or minimize.Topology InputThe execution order of the engines, or local problems, is defined under the<Topology> tag.Merger (Parmela) Parameters x RF (Astra) Separator (MADX) Beam delivery (Astra) E0 ”RF ”Sep Quad strengths Figure B.15: Optimization topology XML block example. The platform basecontains information on all the optimization parameters. The parametersare given to an engine (vertex) when needed.A topology is a directed graph. Fig. B.15 illustrates a sample topology.The round vertex x, represents the input parameters. The square verticesrepresent local problems. To produce an individual, Variator selects startingvalues for the global optimization parameters x. The beam runs throughthe local problems in the defined order. The merger section is modeled first,using Parmela. Then the beam passes through the RF, and so on. Any localproblem can make use of x. For example, we can choose to optimize the RFphase in the cavities, hence the ASTRA input file for the RF local problemcan include the parameter phiRF . The illustrated topology is defined inXML as<Topology>136Input and Output Formats<Vertex name="Merger1" type="PARMELA"inputfolder="m1" prereqs="" timeout="60" /><Vertex name="RF1" type="ASTRA"inputfolder="rf1" prereqs="Merger1" timeout="120" /><Vertex name="S1" type="PARMELA"inputfolder="s1" prereqs="RF1" timeout="600" /><Vertex name="BD1" type="ASTRA"inputfolder="bd1" prereqs="RF1,S1" timeout="60" /></Topology>Each Vertex node must have a type attribute. This can be ASTRA,PARMELA, etc. The Process must also have an inputfolder attribute,which states the relative location of the local input files. For instance,the location of an ASTRA Process A1 must contain the ASTRA input fileastra.in.The prereqs attribute is a comma-delimited list of vertices that mustcomplete before the current vertex can execute. In the above sample, RF1will not start until Merger1 finishes. Merger1 and S1 can run in parallel,since they are independent from each other. NOTE: if vertex A2 uses A1 asa prerequisite, A1 must be listed before A2 in the XML.The timeout attribute denotes the maximum wall-time (in seconds) thevertex is allowed to run for. The engine run is cancelled when this limitis reached, and the vertex must be processed again for the particular indi-vidual. The timeout is a safety feature to prevent a software hang-up fromdestroying the entire optimization run.Often, one engine is used to produce a file that will be used as an inputto another engine. For instance, GPT can be used to produce a particle dis-tribution file, which can then be used as the input file for ASTRA. Copyingfiles can be performed as follows:<Vertex name="v1" type="GPT" inputfolder="gpt1"prereqs="" timeout="600" /><Vertex name="v2" type="ASTRA" inputfolder="astra1"prereqs="v1" timeout="120"><FileRequest from="v1"name="out/dist" newname="in/dist" /></Vertex>In the above code, when the optimization engine begins to execute the enginev2, it will first copy the file out/dist from the directory of the engine v1 intothe directory of the engine v2, and renames the file to in/dist. To make137Input and Output Formatssure the file exists, it is a good habit to put vertex v1 as a pre-requisite forthe vertex v2.Local Engine Input FilesThe vertex input files, or local input files, are files necessary for an engineto run. They are organized into template directories. An example enginemight have an input file that takes in energy as a parameter. The vertextemplate directory could therefore contain the input file with energy as aparameter:Run beam at ${energy} MeV.When this vertex is evaluated for a particular individual, a copy of thetemplate is made and parsed to contain:Run beam at 10.2 MeV.where 10.2 is the optimization generated value of the energy for the partic-ular individual.Another example is an RF section of the linac modeled by ASTRA. Wemust provide the optimization program the path to the ASTRA input file.In the topology definition in Problem.xml, we specified an inputfolder at-tribute for each Process, i.e. local problem. This input folder contains tem-plate files. As an example, the RF model has one template file, astra.in.Within the template file astra.in, we can find the lineMaxE(1)=${MaxE1}where MaxE1 is the name of a parameter (within x or y(x)). When Variatorrequests a new individual to be created, Evaluator makes a directory forthe individual and copies the ASTRA template file to the directory. If fora particular individual, a value of 10 MeV for MaxE1 is selected. Variatormakes the substitution in the copied fileMaxE(1)=10The copied Astra input file can now be read by the program ASTRA. Eval-uator then starts ASTRA for the current local problem with the newlygenerated input file.138Input and Output FormatsLocal Engine OutputThe output of a local problem adheres to the conventions of the physicsengine used in the problem. ASTRA problems create output files accordingto the convention described in the ASTRA manual. In order to read datafrom a local output, the global framework must first check the type of processof the local problem (whether ASTRA or PARMELA, etc.), then follow thefolder and file convention of that process type.When the process of a local problem finishes running, we want to extractuseful information from the output data. Each modeling engine implementsIEngineManager and a post-run function. Within this function, we canextract commonly used information from the output files. As an example,we have ASTRA designed to model the RF section. During post-run, we canextract the energy of the beam and the size of the beam. These extractedvalues belong to the set of non-decision parameters y(x), and are passedback to Variator to be used, for instance, to check whether a constraint onbeam size is satisfied.Custom Code for Local ProblemsCustom code allows the user to manipulate the optimization parameters tosuit their programming needs. This is convenient, for instance, if for anAstra vertex, the user wants to extract the emittance only of particles thatare in front of the beam center.Each vertex contains a file customcode.py (see the folder structure sec-tion). IEngineManager executes the following algorithm:EvaluateVertexRun vertex engineWrites parameter list to customcode.inRun customcode.pycustomcode.py reads from customcode.inManipulates/extracts parameter values as neededOutputs modified values to customcode.outReads customcode.outUpdates parameter values for this individualEndThe python code is executed by system(). Embedded Python was ruledout as an option since the Global Interpreter Lock (GIL) allows only onePython object to be accessed at a time. This effectively makes the opti-mization program single-threading. Custom codes are anticipated to have139Input and Output Formatsmuch I/O interaction with the file system, so can be very computationallycostly. Using System makes this process multi-processing, so bypassing thesingle-threading issue.Global Optimization OutputThe optimization platform writes the population to a history file his. For-mat of the file isLn 1: # commentLn 2: # commentLn 3: id1 x0 x1 ... y0 y1 ... g0 g1 ... F0 F1 ...Ln 4: id2 x0 x1 ... y0 y1 ... g0 g1 ... F0 F1 ......Ln N+2: idN x0 x1 ... y0 y1 ... g0 g1 ... F0 F1 ...The first two lines are comment lines. The following lines list the indi-vidual ID, x, y, F , and g values for each individual of the population. Ndenotes the size of the population.With each generation of the genetic algorithm, the history file is over-written with the new generation. An XML setting HistoryBackupIntervalcan be used to periodically save the history file.Filesystem StructureFigure B.16: Optimization filesystem structure.140Input and Output FormatsFig. B.16 illustrates the folder structure of an optimization problem:1. Problem directory - base directory (RunName in figure) for the currentoptimization run. The global XML input file is located inside. Thisdirectory is unique to each optimization problem.(a) Global input file Problem.xml.(b) optexch is the location for exchange files from Variator and Se-lector, and global output files. The name is carried over fromAPISA.(c) The templates directory contains the engine files that define theproblem topology. In the above sample, there are two engines:TestEngine01 and TestEngine02. Within the folder of each engineare:i. Engine executable.ii. Template input files for the engine. Here it is testcase.in.iii. A customcode.py file, which contains user-written pythoncode to manipulate program output.(d) The output directory stores run files for each individual. Resultsfor each individual are stored its own directory, named after theindividual ID. The contents of each individual directory mimicthat of the template directory, but with parameter values specificto that individual. It is best to make sure the output directory isempty at the beginning of an optimization run. Some engines cancheck for the existence of an output file as a completion condition.A pre-existing output file can cause Evaluator to terminate athread prematurely.2. Executables directory (not shown in figure) containing the compiledC++ binaries (see list B). By default /home/<user>/opt. This pathneeds to be added to the environmental variable LD LIBRARY PATH.3. Directory (not shown in figure) containing python files for common en-gine operations. For example, the file madx.py contains functions formanipulating MADX output. By default /home/<user>/opt/python.This path needs to be added to the environmental variable PYTHONPATH.141Multithreading and Parallel Processing Implementation DetailsMultithreading and Parallel ProcessingImplementation DetailsTo take advantage of high performance computation clusters such as West-Grid [8], the optimization program uses both multithreading and parallelprocessing in its algorithms. The current implementation of the platformrequires that the cluster supports parallel computing, not distributed com-puting, i.e. all processors share the same storage space. The interactionbetween threads and processes are detailed in Fig. B.17.Figure B.17: Optimization multithreading sequence diagram.The above components involve Variator and Evaluator. Selector is single-threaded. The major components of the scheme in Fig. B.17 are1. Main thread - the main Variator thread. Handles all job assignments.When a worker thread finishes, the thread uses the node occupied bythe previous thread, to assign a new job.2. Worker threads - when the main thread assigns a job, it creates aworker thread to execute this job. Multiple worker threads can run142Multithreading and Parallel Processing Implementation Detailsconcurrently. On Westgrid, the maximum number of worker threadsdepends on the number of nodes granted by the Westgrid scheduler.For local runs, the number of nodes can be changed in Problem.xml.When a worker thread is created, it calls fork() to create two pro-cesses. The worker thread (parent process) waits for the child to com-plete, then calls post-run functions (custom code).3. Child process - the child process is created by the worker threadthrough fork(). The child process calls exec() to run the physicsengine, and terminates right after. The child process may be termi-nated prematurely by the main thread if a timeout clause is violated.4. Vertex state - an internal representation that stores the current statusof the job evaluation.Executing an EngineA modeling engine such as ASTRA is executed when a child thread is as-signed a vertex to evaluate. The thread creates a child process using fork().The child process uses execvp() to start the physics engine. For a local run,the command looks likeastra input > outputInternally, Evaluator tokenizes the command before it is passed to execvp:argv = {"astra", "input", ">", "output"}execvp("astra", argv)For a WestGrid run, since we want to tunnel into another processor, thecommand is of the formssh -q nodename cd vertexdir; astra input > outputThis is converted intoargv = {"ssh", "nodename", "cd", ...}execvp("ssh", argv)The command string is tokenized via strtok r(), the thread-safe versionof strtok.143Multithreading and Parallel Processing Implementation DetailsEvaluation TimeoutSince the optimization platform allows for a large array of different modelingengines, many exceptions can occur. This can arise from input files withbad combinations of values or bugs with the engine code, and can causethe engine to hang-up. The platform uses a timeout counter to check forhang-ups.The main thread keeps a list of process IDs for physics engine processes.When a worker thread creates a physics engine process via exec(), theresulting pid is added to this list. The job assignment algorithm in themain thread will periodically check whether any process violates the timeoutcondition. If so, the main thread sends a SIGKILL signal to terminate theprocess. The worker thread which is waiting for the physics engine to finishwill also join. The purpose of the timeout is to reap any process that mighthave encountered an error during execution and hangs as a result.A job which timed out is given the internal status CANCELLED. If thecancellation occurs during the first execution of the job, the job will bere-assigned. If the cancellation occurs during the re-execution, the job iscancelled permanently. The associated individual is given default values forall its non-decision parameters.The timeout value (in seconds) of a vertex is defined in Problem.xml.It is up to the user to determine an upper bound for each vertex. If thetimeout is too short, then it is possible that the physics engine is runningproperly but is not given the time to finish executing.WestGridA major component of the optimization program is to take advantage ofparallel processing on the WestGrid cluster.Variator controls the distribution of nodes. When Variator assigns a jobto Evaluator, it passes a node ID for Evaluator to use. Evaluator then runsthe job on this node. Evaluator is in charge of tracking the status of thisone node only. If the node dies or fails to respond, Evaluator must returnan error code to Variator, and the job is run again.Python on WestGridThe python codebase uses Numpy [5]. To take advantage of Numpy in cus-tom code, please adhere to specific policies regarding Numpy for each West-Grid cluster. On the WestGrid system Orcinus, for instance, the shell com-mand module load python2.4/extramust be run (command can be added144Unitsto .bashrc) first. Otherwise, the python line import numpy returns an error.The cluster-specific policies can be found at http://www.westgrid.ca/sup-port.A Note on Vertices in Parallel and Updating IndividualValuesConsider for instance, two non-decision parameters A1, A2, and two vertices,V1, V2. V1 is charged with updating A1 and V2 for A2. V1 returns valuesA1=1, A2=NaN. The NaN is because V1 does not know about A2. Whenthe engine finishes running, Evaluator will update the individual with thenew values. The individual now has A1=1, A2=NaN. Then V2 finisheswith A1=NaN, A2=1. Evaluator updates the individual to the new values.Instead of A1=1, A2=1, the individual now has A1=NaN, A2=1. Thecorrect value for A1 returned by V1 is overwritten by the NaN from V2.This would occur if the update mechanism is indiscriminant.To prevent good values from being overwritten by bad ones, we have todistinguish which is good. As a convention, an Individual object should startwith random values for decision variables, and INFINITY for non-decisionvariables. The INFINITY value is used as a check when update parametervalues. If an engine returns INFINITY for a particular parameter, it meansthat parameter should updated by a different engine. Therefore, the startingvalues of the Individual object are integral to the operation.This means that any customcode function which returns INFINITYshould be truncated to some finite value (default 1030), so as to not confusebetween a legitimate INF vs a do-not-overwrite INF.UnitsSince the optimization framework is designed to optimize over multiple (per-haps different) physics engines, it has to adapt to the unit convention of eachengine when necessary. For example, Empirical Model uses MeV whereasMADX uses GeV. The optimization platform takes care of the frustrationof unit conversions between engines.Parameters and constraints are defined with a unit attribute in the XMLinput file. Functions related to units and conversion are stored in the opti-mization platform’s Accessibility library. Accessibility keeps a set of units,each with a defined numerical value. For example, the user can define lengthunits ”m”=1 and ”mm”=0.001. Then if a parameter is defined as145Units<Parameter name="x" min="2" max="5" unit="mm" />In internal units, the above parameter can vary from 0.002 to 0.005.When parsing template files, the parameter values must be convertedinto the external units that the engine uses. The template parsing func-tion in Evaluator must remember to call the necessary conversion functions.Each IEngineManager class contains its own set of default units relevant tothe particular engine. These default engine units are stored as settings inProblem.xml.Unit conversions occur in the following areas:ˆ Creating parameters and constraints from the XML file. Need to con-vert values defined in XML into internal units.ˆ Parsing template files. Convert parameter values from internal intoexternal engine units.ˆ Writing customcode.in file. Convert into external units.ˆ Reading from customcode.out file. Convert from external engineunits into internal units.The local problem (including customcode.py) assumes all units are inthe local engine units. The global problem (including the global output file)assumes global units.Unit DefinitionsThe list of available units are defined in the XML file under the <UnitsList>block. For instance:<UnitsList><UnitType name="length" /><Unit name="m" value="1" /><Unit name="cm" value="{m}/1e2" /><Unit name="mm" value="{cm}/10" /></UnitType></UnitsList>All the units and types used in the Parameters and Constraints XMLblocks must be defined here first. The value attribute for a unit can referenceother units, as long as the referenced units are defined first. In the above, ifthe cm unit is defined first, a parsing error would occur. The unit values are146Unitstranslated from XML strings into floating points by 1) using Boost Regex toparse any references to other units, then 2) using Boost Spirit6 to calculatethe value. The units calculator is based on the example VariantCalc [91].Default UnitsThe optimization engine has a list of default units (labeled ”generic”). Inaddition, each engine can have its set of default units. If, for example,the XML does not define a default length unit for ASTRA, then the genericlength will be used. All custom code inputs are in the units of the engine thatcreated them. When reading from custom code outputs are also assumed tobe in the engine default units.Defaults units are defined in XML as:<DefaultUnits><DefUnit engine="GENERIC" type="length"defaultunit="m" /><DefUnit engine="GENERIC" type="angle"defaultunit="deg" /><DefUnit engine="ASTRA" type="length"defaultunit="cm" /></DefaultUnits>Every unit type must have a unit defined under the label GENERIC.This is the unit that the optimization program uses if an engine specificunit cannot be found.When defining parameters in XML, make sure to set numerically stablevalues for units. Upper and lower limits should not be smaller than 10−7.E.g., instead of setting parameter x to be between 10−7 and 10−5 m, definea new unit “um”=m/1e6, and set parameter x to be between 0.1 and 10 µm.This avoids any rounding errors that can occur during computation.147Appendix CList of Parameters andConstants in the ERLOptimizationHere we detail free parameters and their search ranges, beam parameters,and constants, as implemented in the ERL optimization setup. The objec-tives and constraints in the optimization are listed in section 4.Initial Beam ParametersInitial ERL beam parameters are constants taken from E-linac baseline de-signs [29, 34], and are listed in Table C.1. RIB parameters are taken frombaseline designs [31] and ASTRA injector simulations provided by Y.C.Chao. They are listed in Table C.4.Table C.1: Initial ERL beam parameters for optimization.Parameter ValueBunch charge Q 100 pCHorizontal normalized emittance εx,n 10 µmVertical normalized emittance εy,n 10 µmLongitudinal normalized emittance εz,n 47 µmβx 4.787 mαx -1.64881βy 1.08843 mαy 0.91827Bunch length σz 0.0003 mEnergy spread δ 0.019Electron energy Ein 7.5 MeVηx 0 mContinued on next page148RF ParametersTable C.1 – continued from previous pageParameter Valueηx′ 0Table C.2: Initial rare isotope beam parameters for optimization.Parameter ValueBunch charge Q 16 pCHorizontal geometric emittance εx 0.3257 µmVertical geometric emittance εy 0.3257 µmLongitudinal geometric emittance εz 0.331 µmβx 4.8 mαx -1.6βy 1.1 mαy 0.92Bunch length σz 0.00026 mEnergy spread δ 0.04Electron energy Ein 10 MeVηx 0 mηx′ 0RF ParametersCavity amplitude and phase were optimization variables, with all four cavi-ties operating independently from each other. The initial phases of the RIBand ERL beams can be different to represent the independent phasing of theRIB and ERL injectors. The list of RF parameters is shown in Table C.3.The parameter ranges are the maximum ranges as defined by the EmpiricalModel interpolation tables. The RF phases span 20◦ to either side of thecrest.149Lattice ParametersTable C.3: RF optimization parameters.Parameter ValueFrequency 1.3 GHzCavity 1 amplitude [16,20] MV/mCavity 2 amplitude [16,20] MV/mCavity 3 amplitude [16,20] MV/mCavity 4 amplitude [16,20] MV/mCavity 1 phase [310◦,350◦]Cavity 2 phase [310◦,350◦]Cavity 3 phase [310◦,350◦]Cavity 4 phase [310◦,350◦]Lattice ParametersLattice optics are mostly free parameters.Drifts in the merger, linac, and separator are not variables due to existingdesign constraints. Certain quads in the arcs are given higher ranges thanothers due to the excessive demands placed on them. For example, sectionswith one quad between two dipoles are given extra range to compensate forpossible strong dipole effects (see Fig. 4.1). Other quads, while capable ofreaching the same ranges [10], are not given the freedom because we seeksolutions with low to moderate optical demands. Note that the first and lastquads of both arcs have a high positive limit. This is to offset the strongvertical focusing of dipoles to either side of the quads. The dump doublethas higher K1 (see Appendix A for definition) because the beam energy islower in this section. In addition, all quads in the recirculation loop areallowed to move, i.e. the drifts around the quads can change. Drifts areminimum 25 cm to reserve space for diagnostics and other equipment.150Lattice ParametersTable C.4: Lattice optics parameters, listed from upstream to downstream.Drifts are labeled ‘L’. Arc quads are labeled with the system AxyQz, wherex is the arc number, y is the arc drift section (1, 2, or 3), and z is the quadnumber, so A21Q3 is the third quad in the first drift of arc 2.Section Parameter ValueEABT EABTQ1,K1 [-15,15] m−2EABT EABTQ2,K1 [-15,15] m−2EABT EABTQ3,K1 [-15,15] m−2EHAT EHATQ1,K1 [-15,15] m−2EHAT EHATQ2,K1 [-15,15] m−2EHAT (RIB only) EHATQ3,K1 [-15,15] m−2EHAT (RIB only) EHATQ4,K1 [-15,15] m−2Arc 1 matching A1MQ1,K1 [-15,30] m−2Arc 1 matching drift L1 0.25 mArc 1 matching drift L2 0.25 mArc 1 A11Q1,K1 [-15,30] m−2Arc 1 A12Q1,K1 [-15,15] m−2Arc 1 A12Q2,K1 [-15,15] m−2Arc 1 A12Q3,K1 same as A12Q1Arc 1 A13Q1,K1 same as A11Q1Arc 1 drift L1 [0.25,0.45] mArc 1 drift L2 [0.25,0.38] mMirror separator MSEPQ1,K1 same as A1MQ1Chicane matching CHIQ1,K1 [-15,15] m−2Chicane matching CHIQ2,K1 [-15,15] m−2Chicane matching drift L1 [0.25,0.75] mChicane matching drift L2 [0.25,0.75] mChicane drift L1 [0.8,2.0] mFEL matching FELMQ1,K1 [-15,15] m−2FEL matching FELMQ2,K1 [-15,15] m−2FEL matching FELMQ3,K1 [-15,15] m−2FEL matching FELMQ4,K1 [-15,15] m−2FEL matching FELMQ5,K1 [-15,15] m−2FEL matching drift L1 [0.25,0.5] mFEL matching drift L2 [0.25,0.5] mArc 2 matching A2MQ1,K1 [-15,15] m−2Arc 2 matching A2MQ2,K1 [-15,15] m−2Continued on next page151Undulator ParametersTable C.4 – continued from previous pageSection Parameter ValueArc 2 matching A2MQ3,K1 [-15,15] m−2Arc 2 matching A2MQ4,K1 [-15,15] m−2Arc 2 matching A2MQ5,K1 [-15,15] m−2Arc 2 matching drift L1 [0.25,0.5] mArc 2 matching drift L2 [0.25,0.5] mArc 2 A21Q1,K1 [-15,30] m−2Arc 2 A22Q1,K1 [-15,15] m−2Arc 2 A22Q2,K1 [-15,15] m−2Arc 2 A22Q3,K1 same as A22Q1Arc 2 A23Q1,K1 same as A21Q1Arc 2 drift L1 [0.25,0.45] mArc 2 drift L2 [0.25,0.67] mMerger matching EMBQ0,K1 [-15,15] m−2Merger matching EMBQ1,K1 [-15,15] m−2Merger matching EMBQ3,K1 [-15,15] m−2Merger matching EMBQ5,K1 [-15,15] m−2Merger matching EMBQ7,K1 [-15,15] m−2Linac matching EMBTQ6,K1 [-15,15] m−2Linac matching EMBTQ7,K1 [-15,15] m−2EDBT MQD1,K1 [-60,60] m−2EDBT MQD2,K1 [-60,60] m−2- Arc dipoles edge angle [15◦,30◦]- Chicane dipoles bend angle [10◦,25◦]- Variable length drifts [0.0,0.058] mA series of variable length drifts are inserted into various points on thelattice to provide further freedom in choosing the recirculation time-of-flight.Undulator ParametersThe list of undulator parameters is shown in Table C.5. Some parameterswere derived previously in section 4.152Undulator ParametersTable C.5: List of undulator parameters.Parameter ValueUndulator parameter K 0.7Beam energy 45 MeVBeam power 0.5 MWBunch charge 100 pCTransverse emittance, normalized 10 µmLongitudinal emittance, normalized 80 keV-psRMS energy spread 0.004RMS bunch length 1 psUndulator period 4 cmNumber of undulator periods 25Undulator length 1 mRadiation wavelength 4 µmintracavity saturated power 8 MW153Appendix DERL Major Components andLayoutParameters and coordinates for the ERL baseline presented in chapter 6 areshown here. The format of this section follows that of the E-Linac PhaseOne design note [33].Beam and Laser RequirementsERL beam properties are listed in Table D.1. Rare isotope beam propertiesare listed in Table D.2.Table D.1: ERL beam parameters.Parameter ValueHorizontal normalized emittance εx,n 10 µmVertical normalized emittance εy,n 10 µmLongitudinal normalized emittance εz,n 47 µmBunch length σz 0.0003 mEnergy spread δ 0.019Electron energy before acceleration 7.5 MeVElectron energy after acceleration 45.81 MeVElectron energy after lasing 45.78 MeVElectron energy after deceleration 7.7 MeVTable D.2: RIB parameters.Parameter ValueHorizontal geometric emittance εx 0.3257 µmContinued on next page154RF RequirementsTable D.2 – continued from previous pageParameter ValueVertical geometric emittance εy 0.3257 µmLongitudinal geometric emittance εz 0.331 µmBunch length σz 0.00026 mEnergy spread δ 0.04Electron energy before acceleration 10 MeVElectron energy after acceleration 48.3 MeVFEL parameters are listed in Table D.3. Gain is defined as (dP/dz)/P ,where P is the radiation power. Undulator geometric parameters are listedin Table C.5 and are not repeated here.Table D.3: FEL parameters.Parameter ValueGain 0.5 m−1Laser wavelength 3.8 µmRF RequirementsRF requirements are listed in tableD.4 for the four 1.3 GHz cavities (EACAand EACB each contain two cavities). The cavities are phased indepen-dently. ERL and RIB injectors are also phase independently.Convention: phases labeled “initial” are defined with the bunch centroidat the entry point of EACA. This represents the phases of the cavities at thesame time. Phases labeled “entrance” are defined with the bunch centroidat the respective cavity entrance. Under this convention 335◦ is the crest.Table D.4: RF requirements.Cavity ERL phase (deg) RIB phase (deg)EACA:CAV1, initial 328.14 337.27Continued on next page155Optics RequirementsTable D.4 – continued from previous pageCavity ERL phase (deg) RIB phase (deg)EACA:CAV2, initial 164.32 173.45EACB:CAV3, initial 228.41 237.54EACB:CAV4, initial 57.73 66.86EACA:CAV1, entrance 328.14 337.27EACA:CAV2, entrance 332.64 340.79EACB:CAV3, entrance 337.67 345.57EACB:CAV4, entrance 333.04 340.89Optics RequirementsERL quad requirements are listed in Table D.5. All quads are assumed tobe bipolar. Quadrupole strengths are listed in KdL, where K = K1 givenin Appendix A and dL is the effective length.Table D.5: ERL quadrupole requirements.Quad KdL (G)EABT:QEABTQ1 93EABT:QEABTQ2 -677EABT:QEABTQ3 682EHAT:QEHATQ1 -40EHAT:QEHATQ2 -410ARC1:QA1M1 3474ARC1:QA111 4711ARC1:QA121 1890ARC1:QA122 -1315ARC1:QA123 1890ARC1:QA131 4711SEP2:QA1M1 3474SEP2:QEHATQ2 -410FELM:QFELM1 2062FELM:QFELM2 -803FELM:QFELM3 -468FELM:QFELM4 -1453Continued on next page156Optics RequirementsTable D.5 – continued from previous pageQuad KdL (G)FELM:QFELM5 2385A2M:QA2M1 1019A2M:QA2M2 -3406A2M:QA2M3 3370A2M:QA2M4 -2903A2M:QA2M5 1477ARC2:QA211 6563ARC2:QA221 3205ARC2:QA222 -2728ARC2:QA223 3205ARC2:QA231 6563MERG:QEMBQ0 2766MERG:QEMBQ1 -645MERG:QEMBQ3 -172MERG:QEMBQ5 -270MERG:QEMBQ7 1099MERG:QEMBTQ6 -3395MERG:QEMBTQ7 2124EDBT:QMQD1 521EDBT:QMQD2 -159RIB quad requirements are listed in Table D.6. Contains overlap withERL quads. All quads are assumed to be bipolar.Table D.6: RIB quadrupole requirements.Quad KdL (G)EABT:QEABTQ1 93EABT:QEABTQ2 -677EABT:QEABTQ3 682EHAT:QEHATQ1 -40EHAT:QEHATQ2 -410EHAT:QEHATQ3 1137EHAT:QEHATQ4 1012157Optics RequirementsTable D.7: ERL dipole requirements.DipolePath length(m)Bend angle(deg)Entry angle(deg)Exit angle(deg)BdL(G-cm)EHAT:BRF 0.5000 0.27 0.00 0.00 715EHAT:BBSTR 0.3000 0.99 -0.25 1.18 2633EHAT:BBSEP 0.5000 8.36 0.00 8.36 22306ARC1:BA11 0.2906 40.19 29.99 29.99 107183ARC1:BA12 0.2906 40.19 29.99 29.99 107183ARC1:BA13 0.2906 40.19 29.99 29.99 107183ARC1:BA14 0.2906 40.19 29.99 29.99 107183SEP2:BBSEP 0.5000 8.36 8.36 0.00 22306SEP2:BBSTR 0.3000 1.26 1.18 -0.25 3348CHI:BCHI1 0.3500 23.15 0.00 23.15 61747CHI:BCHI2 0.3500 -23.15 -23.15 0.00 -61747CHI:BCHI3 0.3500 -23.15 0.00 -23.15 -61747CHI:BCHI4 0.3500 23.15 23.15 0.00 61747ARC2:BA21 0.2906 45.00 29.99 29.99 119936ARC2:BA22 0.2906 45.00 29.99 29.99 119936ARC2:BA23 0.2906 45.00 29.99 29.99 119936ARC2:BA24 0.2906 45.00 29.99 29.99 119936MERG:BMBA1 0.1500 -4.61 -2.30 -2.30 -12281MERG:BMBA2 0.3000 9.22 4.61 4.61 24563MERG:BMA3 0.1500 -4.61 -2.30 -2.30 -12281EDBT:BBSTR 0.3000 6.00 0.00 6.00 2720ERL dipole requirements are listed in Table D.7. All eight arc dipoles(four per arc) have the same geometry. The four arc 1 dipoles operate atlower strengths than the arc 2 dipoles since they do not need to account forthe full 180◦ turn. Some of the turn is absorbed by the separator and mirrorseparator.All four chicane dipoles have the same geometry.158Installation Coordinates TableInstallation Coordinates TableWe provide the list of ERL transport elements and their layout in the TRI-UMF E-hall. The list begins with the first linac cavity. In all the following:ˆ Uses TRIUMF standard Cyclotron Center coordinates (Cyclotron Cen-ter is (0,0))ˆ S: Cumulative path length with S=0 at RIB injector cathodeˆ X: +:East; -:West in Cyclotron Center coordinatesˆ Y: +:North; -:South in Cyclotron Center coordinatesˆ Naming convention: elements beginning with ’O’ are drifts; ’Q’ arequads; ’B’ are dipolesˆ Coordinates refer to the end of the elementTable D.8: ERL element coordinates. Layout shown in Fig. 4.1.Element S (m) X (m) Y (m)EACA:START 9.1595 -36.8470 -3.0114EACA:CAV1 10.4395 -36.8470 -1.7314EACA:O12 11.0895 -36.8470 -1.0814EACA:CAV2 12.3695 -36.8470 0.1986EABT:ODN04 12.7752 -36.8470 0.6043EABT:QEABTQ1 12.8552 -36.8470 0.6843EABT:ODN05A 13.0895 -36.8470 0.9186EABT:BEABD 13.2395 -36.8470 1.0686EABT:ODN05A 13.4737 -36.8470 1.3028EABT:QEABTQ2 13.5537 -36.8470 1.3828EABT:ODN06 13.7737 -36.8470 1.6028EABT:QEABTQ3 13.8537 -36.8470 1.6828EABT:ODN07 14.2594 -36.8470 2.0885EACB:CAV3 15.5395 -36.8470 3.3686EACB:O34 16.1895 -36.8470 4.0186EACB:CAV4 17.4695 -36.8470 5.2986EHAT:ODN08 17.9009 -36.4156 5.2986EHAT:QEHATQ1 17.9809 -36.3356 5.2986EHAT:ODN09 18.2149 -36.1016 5.2986EHAT:BRF 18.7149 -35.6016 5.2974EHAT:BBSTR 19.0149 -35.3016 5.2934EHAT:ODN11 20.0149 -34.3016 5.2715Continued on next page159Installation Coordinates TableTable D.8 – continued from previous pageElement S (m) X (m) Y (m)EHAT:QEHATQ2 20.0949 -34.2216 5.2698EHAT:ODN12 21.0949 -33.2216 5.2479EHAT:BBSEP 21.5949 -32.7216 5.2005ARC1:OA1M1 21.8449 -32.4716 5.1588ARC1:QA1M1 21.9949 -32.3216 5.1337ARC1:OA1M2 22.2449 -32.0716 5.0919ARC1:BA11 22.5355 -31.7810 4.9508ARC1:OA111 22.9051 -31.4114 4.6685ARC1:QA111 23.0551 -31.2614 4.5539ARC1:OA112 23.3855 -30.9310 4.3015ARC1:BA12 23.6761 -30.6404 4.0341ARC1:OA121 23.9970 -30.3195 3.7133ARC1:QA121 24.1470 -30.1695 3.5633ARC1:OA122 24.5010 -29.8155 3.2093ARC1:QA122 24.6510 -29.6655 3.0593ARC1:OA123 25.0050 -29.3115 2.7053ARC1:QA123 25.1550 -29.1615 2.5553ARC1:OA124 25.4758 -28.8407 2.2344ARC1:BA13 25.7664 -28.5501 1.9670ARC1:OA131 26.0968 -28.2197 1.7147ARC1:QA131 26.2468 -28.0697 1.6001ARC1:OA132 26.6164 -27.7001 1.3177ARC1:BA14 26.9070 -27.4095 1.1766ARC1:OA1M2 27.1570 -27.1595 1.1348ARC1:QA1M1 27.3070 -27.0095 1.1098ARC1:OA1M1 27.5570 -41.0776 9.4208SEP2:BBSEP 28.0570 -41.1249 8.9235SEP2:ODN12 29.0570 -41.1469 7.9237SEP2:QEHATQ2 29.1370 -41.1486 7.8437SEP2:ODN11 30.1370 -41.1705 6.8440SEP2:BBSTR 30.4370 -41.1738 6.5440CHI:OCHI1 30.9195 -41.1738 6.0616CHI:QCHI1 31.0695 -41.1738 5.9116CHI:OCHI2 31.5879 -41.1738 5.3931CHI:QCHI2 31.7379 -41.1738 5.2431CHI:OCHI3 32.4875 -41.1738 4.4935Continued on next page160Installation Coordinates TableTable D.8 – continued from previous pageElement S (m) X (m) Y (m)CHI:BCHI1 32.8375 -41.1040 4.1529CHI:OCHI11 33.6607 -40.7804 3.3961CHI:BCHI2 34.0107 -40.7106 3.0555CHI:OCHI21 37.3048 -40.7106 -0.2386CHI:BCHI3 37.6548 -40.7804 -0.5792CHI:OCHI31 38.4780 -41.1040 -1.3360CHI:BCHI4 38.8280 -41.1738 -1.6766FELM:OFELM1 39.1439 -41.1738 -1.9924FELM:QFELM1 39.2939 -41.1738 -2.1424FELM:OFELM2 39.5912 -41.1738 -2.4398FELM:QFELM2 39.7412 -41.1738 -2.5898FELM:OFELM3 40.0386 -41.1738 -2.8871FELM:QFELM3 40.1886 -41.1738 -3.0371FELM:OFELM4 40.4859 -41.1738 -3.3345FELM:QFELM4 40.6359 -41.1738 -3.4845FELM:OFELM5 40.9333 -41.1738 -3.7818FELM:QFELM5 41.0833 -41.1738 -3.9318FELM:OFELM6 41.5960 -41.1738 -4.4445FEL:UND 42.5960 -41.1738 -5.4445A2M:OA2M1 42.9576 -41.5354 -5.4445A2M:QA2M1 43.1076 -41.6854 -5.4445A2M:OA2M2 43.4235 -42.0014 -5.4445A2M:QA2M2 43.5735 -42.1514 -5.4445A2M:OA2M3 43.8895 -42.4673 -5.4445A2M:QA2M3 44.0395 -42.6173 -5.4445A2M:OA2M4 44.3554 -42.9332 -5.4445A2M:QA2M4 44.5054 -43.0832 -5.4445A2M:OA2M5 44.8213 -43.3991 -5.4445A2M:QA2M5 44.9713 -43.5491 -5.4445A2M:OA2M6 45.3460 -43.9238 -5.4445ARC2:BA21 45.6366 -44.2144 -5.3362ARC2:OA211 46.0518 -44.6296 -5.0426ARC2:QA211 46.2018 -44.7796 -4.9365ARC2:OA212 46.4866 -45.0644 -4.7351ARC2:BA22 46.7772 -45.3550 -4.4735ARC2:OA221 47.1327 -45.7106 -4.1179Continued on next page161Installation Coordinates TableTable D.8 – continued from previous pageElement S (m) X (m) Y (m)ARC2:QA221 47.2827 -45.8606 -3.9679ARC2:OA222 47.8945 -46.4724 -3.3561ARC2:QA222 48.0445 -46.6224 -3.2061ARC2:OA223 48.6563 -47.2341 -2.5944ARC2:QA223 48.8063 -47.3841 -2.4444ARC2:OA224 49.1619 -47.7397 -2.0888ARC2:BA23 49.4525 -48.0303 -1.8272ARC2:OA231 49.7373 -48.3151 -1.6258ARC2:QA231 49.8873 -48.4651 -1.5197ARC2:OA232 50.3025 -48.8803 -1.2261ARC2:BA24 50.5931 -49.1709 -1.1177MERG:ODQE 50.9931 -49.5709 -1.1177MERG:QEMBQ0 51.1431 -49.7209 -1.1177MERG:OMQ 51.6231 -50.2009 -1.1177MERG:QEMBQ1 51.7731 -50.3509 -1.1177MERG:OMQ 52.2531 -50.8309 -1.1177MERG:QEMBQ3 52.4031 -50.9809 -1.1177MERG:OMQ 52.8831 -36.8470 -5.9045MERG:QEMBQ5 53.0331 -36.8470 -5.7545MERG:OMQ 53.5131 -36.8470 -5.2745MERG:QEMBQ7 53.6631 -36.8470 -5.1245MERG:ODD 53.7631 -36.8470 -5.0245MERG:BMBA1 53.9131 -36.8410 -4.8747MERG:OD1 54.0922 -36.8266 -4.6961MERG:BMBA2 54.3922 -36.8266 -4.3965MERG:OD1 54.5714 -36.8410 -4.2179MERG:BMA3 54.7214 -36.8470 -4.0681MERG:ODD2 54.7714 -36.8470 -4.0181MERG:ODN01 55.0094 -36.8470 -3.7801MERG:QEMBTQ6 55.0894 -36.8470 -3.7001MERG:ODN02 55.3103 -36.8470 -3.4791MERG:QEMBTQ7 55.3903 -36.8470 -3.3991MERG:ODN03 55.7780 -36.8470 -3.0114EACA:CAV1 57.0580 -36.8470 -1.7314EACA:O12 57.7080 -36.8470 -1.0814EACA:CAV2 58.9880 -36.8470 0.1986Continued on next page162Installation Coordinates TableTable D.8 – continued from previous pageElement S (m) X (m) Y (m)EABT:ODN04 59.3937 -36.8470 0.6043EABT:QEABTQ1 59.4737 -36.8470 0.6843EABT:ODN05A 59.7080 -36.8470 0.9186EABT:BEABD 59.8580 -36.8470 1.0686EABT:ODN05A 60.0923 -36.8470 1.3028EABT:QEABTQ2 60.1723 -36.8470 1.3828EABT:ODN06 60.3923 -36.8470 1.6028EABT:QEABTQ3 60.4723 -36.8470 1.6828EABT:ODN07 60.8780 -36.8470 2.0885EACB:CAV3 62.1580 -36.8470 3.3686EACB:O34 62.8080 -36.8470 4.0186EACB:CAV4 64.0880 -36.8470 5.2986EHAT:ODN08 64.5195 -36.8470 5.7300EHAT:QEHATQ1 64.5995 -36.8470 5.8100EHAT:ODN09 64.8334 -36.8470 6.0440EDBT:BRF 65.3334 -36.8470 6.5440EDBT:BBSTR 65.6334 -36.8627 6.8434EDBT:ODQ40 66.0334 -36.9045 7.2413EDBT:BBSTR2 66.3334 -36.9359 7.5396EDBT:ODQ 66.4334 -36.9463 7.6391EDBT:QMQD1 66.5134 -36.9547 7.7186EDBT:ODQ40 66.9134 -36.9965 8.1164EDBT:QMQD2 66.9934 -37.0049 8.1960EDBT:ODQ40 67.3934 -37.0467 8.5938EDBT:BDMP 67.3944 -37.0468 8.5948163


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items