Dr Shinji Machida

Dr Suzie Sheehy

Dr Yoshihiro Ishi

Pr Pr Chris

Mori Yoshiharu

Dr Scott Berg

Pr Antoine Philippe Cerfon

Piot

Yann Dutheil

Karim Hamdi

Charles Kitegi

Nathalie Bouet

Dean Hidas Silvia Verdu Andres

Nick Tsoupas

Nick Simos

Mei Bai

Sylvaine Harold Kirk

François Méot

Steve Peggs

Dejan

Kentaro Mihara

Kai Shih

Irina Petrushina

Hui Yuan

Pierre Wu

Dhananjay Ravikumar I Korysko

Wish

you all a very successful completion of your PhD. Lastly, I would like to thank my family for all their love and encouragement. Words cannot express how grateful I am to my mother, Aicha Haj Tahar, for guiding me as a person and for all the sacrifices she made on my behalf. To my older brother, Mejed Haj Tahar, thank you for being another huge part in my foundation, for all of the advice and wise words you've provided me and for

Abstract

High power proton accelerators allow providing, by spallation reaction, the neutron fluxes necessary in the synthesis of fissile material, starting from Uranium 238 or Thorium 232. This is the basis of the concept of sub-critical operation of a reactor, for energy production or nuclear waste transmutation, with the objective of achieving cleaner, safer and more efficient process than today's technologies allow.

Designing, building and operating a proton accelerator in the 500-1000 MeV energy range, CW regime, MW power class still remains a challenge nowadays. There is a limited number of installations at present achieving beam characteristics in that class, e.g., PSI in Villigen, 590 MeV CW beam from a cyclotron, SNS in Oakland, 1 GeV pulsed beam from a linear accelerator, in addition to projects as the ESS in Europe, a 5 MW beam from a linear accelerator.

Furthermore, coupling an accelerator to a sub-critical nuclear reactor is a challenging proposition: some of the key issues/requirements are the design of a spallation target to withstand high power densities as well as ensure the safety of the installation.

These two domains are the grounds of the PhD work: the focus is on the high power ring methods in the frame of the KURRI FFAG collaboration in Japan: upgrade of the installation towards high intensity is crucial to demonstrate the high beam power capability of FFAG. Thus, modeling of the beam dynamics and benchmarking of different codes was undertaken to validate the simulation results. Experimental results revealed some major losses that need to be understood and eventually overcome. By developing analytical models that account for the field defects, one identified major sources of imperfection in the design of scaling FFAG that explain the important tune variations resulting in the crossing of several betatron resonances. A new formula is derived to compute the tunes and properties established that characterize the effect of the field imperfections on the transverse beam dynamics. The results obtained allow to develop a correction scheme to minimize the tune variations of the FFAG. This is the cornerstone of a new fixed tune non-scaling FFAG that represents a potential candidate for high power applications.

As part of the developments towards high power at the KURRI FFAG, beam dynamics studies have to account for space charge effects. In that framework, models have been installed in the tracking code ZGOUBI to account for the self-interaction of the particles in the accelerator. Application to the FFAG studies is shown.

Finally, one focused on the ADSR concept as a candidate to solve the problem of nuclear waste. In order to establish the accelerator requirements, one compared the performance of ADSR with other conventional critical reactors by means of the levelized cost of energy. A general comparison between the different accelerator technologies that can satisfy these requirements is finally presented.

In summary, the main drawback of the ADSR technology is the high Levelized Cost Of Energy compared to other advanced reactor concepts that do not employ an accelerator. Nowadays, this is a show-stopper for any industrial application aiming at producing energy (without dealing with the waste problem). Besides, the reactor is not intrinsically safer than critical reactor concepts, given the complexity of managing the target interface between the accelerator and the reactor core.

To the memory of my father, Hedi Haj Tahar. Unfortunately he could not see this process through to its completion. However, he offered all the support to make it possible. I am eternally grateful. q particle charge (same species) r i particle trajectory bending radius

x particle horizontal coordinate in the Frenet-Serret frame y particle vertical coordinate in the Frenet-Serret frame z longitudinal position of the particle relative to the center of the bunch (X, Y ) cartesian coordinates in the lab frame You may delay, but time will not, and lost time is never found again.

-B. Franklin.

INTRODUCTION

In recent years, the concept of Accelerator Driven Sub-critical Reactor (ADSR) has gained more interest worldwide as a potential candidate to solve the problem of nuclear waste. This concept, the focus of the present dissertation, requires the coupling of a high power proton accelerator with a sub-critical core. Among the many options that can serve as a proton driver, the Fixed Field Alternating Gradient (FFAG) accelerators. The work presented in this thesis is dedicated to demonstrate the high power capability of FFAGs.

Context

Although the nuclear energy has one of the lowest impacts on the environment of any energy source, the general public is very sceptical on this, because nuclear is in general associated with accidents or with High Level Waste (HLW). Nuclear wastes, in particular the spent nuclear fuel, have been a major concern worldwide since the early days of civilian nuclear power: In the mid 50s, when early commercial nuclear power reactors were built in the US, it was assumed that the spent fuel would be reprocessed within a few years after discharge to recover uranium and plutonium using technologies that were developed for the Manhattan project. For this particular reason, spent fuel pools were designed and sized to accept the waste of up to three cores, no more. However, soon, this solution was doomed to failure as one facet of US non proliferation policy. In the 70s, most scientists, regulators and proponents of the nuclear energy realized that high-level waste is a major issue that needs to be solved in order for the nuclear energy to become a long-term option in the future. Discussions became very ethical as illustrated in the 1995 NEA report: "those who generate the wastes should take responsibility, and provide the resources, for the management of these materials in a way which will not impose undue burdens on future generations". This is also the central message of the encyclica letter(LAUDATO SI', § 184 [START_REF]ENCYCLICAL LETTER LAUDATO SI' OF THE HOLY FATHER FRANCIS[END_REF]): "In the face of possible risks to the environment which may affect the common good now and in the future, decisions must be made based on a comparison of the risks and benefits foreseen for the various possible alternatives ... Some projects, if insufficiently studied, can profoundly affect the quality of life of an area due to very different factors such as unforeseen noise pollution, the shrinking of visual horizons, the loss of cultural values, or the effects of nuclear energy use. The culture of consumerism, which prioritizes short-term gain and private interest, can make it easy to rubber-stamp authorizations or to conceal information.". The problem became political rather than technical and, in 1982, the Nuclear Waste Policy Act (NWPA) was passed in the US which established a comprehensive national program for the safe, permanent disposal of highly radioactive wastes: according to this act, the Department of Energy (DOE) has the responsibility to site, build, and operate a deep geologic repository for the disposal of HLW and spent nuclear fuel [START_REF] Epa | Summary of the Nuclear Waste Policy Act[END_REF]. A Nuclear Waste Fund was enacted which required the operators of nuclear plants to pay a tenth of a cent per kilowatt-hour to the government in return for DOE taking responsibility for spent nuclear fuel. However, prior to the waste disposal, the spent fuel remains the responsibility of the utility. Several amendments to the NWPA were added: as of 1987, amendment directed the DOE to investigate only one potential site for HLW disposal: the Yucca mountain in Nevada. In 2008, two years before it was expected to open, the Obama administration closed out funding for the project which has already cost $12 billion. Although there is no technical, scientific or safety basis for this decision, which we will not further discuss given its political facet, it is important to mention that no alternative solution was proposed in the US since then. Besides, the quantities of civilian waste generated continue to pile-up and now they exceed the storage capacity of the Yucca mountain cancelled plant. This means that already a second repository is needed although the first one is not built yet. In the meantime, the nuclear utilities are required to keep their waste on-site and as soon as the DOE failed to meet its obligations, they started suing it. On May 2014, a federal court ruled that DOE's current waste plans were too vague to allow a reasonable estimate to be calculated. Therefore, DOE was obliged to stop collecting fees until provisions are made to collect nuclear waste. An excellent overview of the political aspects of civilian nuclear waste disposal can be found in [START_REF] Holt | Civilian nuclear waste disposal[END_REF].

In order to review the different waste management options, one first needs to be specific about the waste definition.

Definition of waste

The definition of nuclear waste depends on the strategy of the country: if the spent fuel is not reprocessed, as is the case in USA and Sweden, then the entire fuel is considered as High Level Waste (HLW). If, on the other hand, the used fuel is reprocessed such as is the case in France, UK, Russia and Japan, then the waste is treated as a resource not as a waste. Reprocessing is a technology that consists in separating some long-lived radioactive elements from the used fuel, allowing the uranium and plutonium to be recycled. This technology is problematic given that it enhances the possibility for nuclear proliferation. Besides, the cost of fabricating the recovered uranium from reprocessing into fresh fuel is more than that of the fresh uranium from ore. This justifies why many countries, even after opting for this technology, are stockpiling the recovered uranium rather than using it, awaiting for the price of natural uranium to increase.

In order to evaluate the performance of a nuclear fuel cycle, one key parameter is the fuel utilization. Fig. 1.1 shows the highly inefficient use of the fuel resources in a typical Once-Through fuel Cycle (OTC): the used fuel typically contains ∼ 3% of fission products and minor actinides. These elements are poisonous and make the control of the reactor kinetics difficult, hence it is no longer practical to continue using the fuel. In the United States, Light Water Reactors (LWR) are widely employed for generating electricity based on the OTC strategy. A typical mass flow of Pressurized Water Reactor (PWR) is summarized in Table 1. [START_REF]ENCYCLICAL LETTER LAUDATO SI' OF THE HOLY FATHER FRANCIS[END_REF] where one can see that only 5% of the initial uranium resources are used. These estimates also show that, approximately, 2000 tons of nuclear waste are produced per year in the US (for a total power capacity of 100 nuclear reactors i.e. ∼ 100 GWe). Now, in order to assess the quantities of waste produced, one needs to look at the path of minor actinides in the fuel cycle. This is illustrated in Fig. 1.2: the two fuel cycles, uranium and thorium, are shown in red and green respectively. The formation of the Minor Actinides (MA) results from neutron captures, β-decay as well as α-decay. The elements above curium are usually neglected, given that they are produced in

Options for nuclear waste disposal

There exists several options for nuclear waste disposal:

• Ocean disposal: from 1946 through 1993, this was an option for nuclear waste disposal. There exists several damping locations in the sea. However, it is almost impossible to obtain good estimates about the quantities. Since 1993, this method is banned by international treaties. Nowadays, the annual production of 2000 tons of spent fuel in the US is approximately equivalent to the nuclear fallout of 1000 Hiroshima bombs. So, this method cannot be envisioned because its environmental footprint would be disastrous.

• Shipping into outer space: although this method is feasible (economically as well), the main issue is the risk of failure during the launching process of the rocket. This method is therefore disregarded.

• Geologic repository: This solution requires more time for the spent fuel to become stable than any empire or state has ever lasted. However, it is envisaged in most proposed solutions due to its economic competitiveness.

• Reprocessing: as discussed in the previous paragraph, several countries opted for this strategy as part of their plans for closing the fuel cycle. The main issue remains the high cost of the technology.

• Advanced reactor concepts closing the fuel cycle (transmutation is part of this). Mostly fast reactors are considered to reduce the total radiotoxicity of the nuclear waste. However they cannot accommodate large fraction of minor actinides. This justifies the use of sub-critical reactors, hence the need for an external source of neutrons. We discuss this concept in the next section.

1.4 Accelerator Driven Sub-critical Reactor (ADSR)

Overview

The accelerator driven sub-critical reactor is a hybrid technique combining a particle accelerator with a subcritical core. Most ADSR proposals assume proton accelerators delivering a continuous beam with an energy around 1GeV. The protons strike a heavy metal target, in solid or liquid state, ejecting neutrons among other particles in the forward direction with a lower energy than the incident particle has. These neutrons will also strike other nuclides, ejecting other lower energy particles. The cascade continues until all the energy is spent.

For a 1 GeV proton beam impinging on a lead target, the cascade length is around 1 m. This process is referred to as the spallation process. The spallation neutrons are then introduced into the sub-critical core to induce further nuclear reactions. The reactor core could be thermal or fast. The accelerator could be linear or circular. All these aspects will be discussed in detail in chapter 2. A conceptual diagram of ADSR is shown in Fig. 1.3.

We distinguish three potential missions of ADSR:

• The first mission is to eliminate the most offending isotopes prior to storage in a geological repository: bombardment with a large number of neutrons is a way to transmute radiotoxic elements (fission products, actinides) into less dangerous species. This concept, originally due to H. Takahashi [START_REF] Tkahashi | Actinide Transmutation by the Spallation Process[END_REF], is referred to in the literature as Accelerator Transmutation of Waste.

• Another role for ADS is in a breed-and-burn fuel cycle using fertile resource feed. The concept has been discussed since many years mostly for breeding fissile material from fertile thorium and uranium. These materials are provided in the fuel, or in a breeder blanket surrounding the core or in both. This idea originally developed in 1993 when Carlo Rubbia suggested to use ADSR as an energy amplifier, i.e. to produce more energy with fission than the energy needed to power the accelerator.

• The third mission is that of a fissile pre-breeder for critical reactor fleet. This is seen as a solution for the uranium scarcity problem. Production of fissile material occurs by neutron irradiation of fertile material, particularly Uranium-238 and Thorium-232. The feed material contains no natural fissile isotopes. This process can be summarized by the following chain of events:

Concepts of ADSR

Different concepts of accelerator driven systems exist in the literature. The definition of the concept relies on the basic criteria of the sub-critical blanket:

• Type of fuel: Uranium-Plutonium, Thorium-Uranium, MA, etc. • Neutron spectrum: comparing the ratio of fission to absorption, it can be seen that a fast spectrum has an advantage in the transmutation. However, a thermal spectrum may achieve the same performance if neutron fluxes are comparable to those of fast reactors, i.e 10 15 n/cm 2 /s.

• The effective neutron multiplication factor k ef f is the average number of neutrons from one fission that cause another fission. This parameter is fixed by the geometry of the blanket and the type of fuel in use. If k ef f <1, the system is called "sub-critical": it cannot sustain any chain reaction. For an ADS, this is a requirement, and it will be discussed later on in this document. But, for the moment, note that this factor can range from 0 to 1 and should be carefully chosen as a compromise between safety issues, technological limitations and cost of the system.

• Type of fuel (solid/liquid): a wide experience is acquired on solid fuel, although liquid fuels, such as molten salt, present the advantage of online processing. This would reduce the inventory of fission products, control corrosion and improve neutron economy by removing fission products with high neutron absorption cross-section, especially xenon. However, there are many challenges there.

• Power of the installation: depending on the goals foreseen (burning of nuclear waste, energy production, both,...), the reactor power, and thus the beam power required are determined.

Now, we will present some examples of the most important ADS concepts in the world, since the early 90s:

PHOENIX Concept (Brookhaven National Laboratory), 1991

The PHOENIX concept [START_REF] Vantuyle | The PHOENIX Concept[END_REF] proposed a large proton accelerator to drive a sub-critical reactor core which is mainly composed of MAs (Np, Am, Cm). This concept assumes a linear accelerator that can accelerate a proton beam up to 1.6 GeV with an average intensity of 104 mA. After passing through a beam expander, the protons will hit up to eight target modules rendering a total power of 3600 M W th , as indicated in (Fig. 1.4). Each target module would have fuel containing MAs oxide. These are cooled using liquid sodium. The PHOENIX machine was designed to play a major role within a waste-reduction plan called the Clean Use of Reactor Energy (CURE) proposed by the Westinghouse Hanford Corporation.

Molten Salt ADS Concept (Los Alamos National Laboratory), 1992

Los Alamos proposed several ADS concepts in the 1990s [7] [8]. This diversity is explained by the desire to meet the objectives of various civil and military programs.The initial project was aimed at burning MAs and The PHOENIX concept [START_REF] Vantuyle | The PHOENIX Concept[END_REF].

fission products 99 T c and 129 I in order to reduce the period of isolation of the waste from 10 000 years to a few hundred years. The most ambitious proposal consisted in a molten salt accelerator-driven waste burner proposed in Ref. [7] and dedicated to the destruction of transuranic elements. The system uses a thermal flux of neutrons in the 10 16 n/cm 2 /s range in a sub-critical reactor (k ef f = 0.96). The reactor core concept is composed of fluoride salt flowing through a lattice of graphite channels (= moderators). The expected efficiency for conversion of thermal to electric power is 44 %, and the efficiency of the accelerator is 45%. The latter operates at an energy of 1 GeV and a beam current of 49.3 mA.

Rubbia Energy Amplifier (CERN), 1993

In 1993, CERN proposed a hybrid system based on thermal neutrons, whose aim was to produce, from Thorium, clean energy (less radioactive waste) while being immune from proliferation risks [START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF]. This system was called "Energy Amplifier (EA)" in that the beam power is amplified by the system. The EA operates indefinitely in a closed fuel cycle, namely, the discharge of the fuel load with the exception of short lived fission fragments is re-injected in the blanket with the addition of natural Thorium to compensate for the fuel burn-up. The EA is a breed-and-burn application of an ADS.

In the first variant using thermal neutrons, the flux was limited to 10 14 n/cm 2 /s to avoid an indesirable increase in reactivity due to 233 P a decay (in case of beam shutdown), and also to avoid premature destruction of 233 P a through captures before decaying to 233 U . The EA module consists of a 200 MWth unit (k ef f = 0.93) with a dedicated 800 MeV proton accelerator of 6.25 mA. In that case, a cyclotron was envisaged. The burn-up increase which is typically obtained in a PWR thanks to enrichment is here guaranteed by increasing the beam current to compensate for the decrease of k ef f due to the build-up of fission products (no online processing was envisaged here). The maximum current achievable with a cyclotron posed a serious limitation to this proposal.

To overcome a number of limitations, the design evolved to opt for a fast spectrum blanket. In this case, the concentration ratio 233 U -to-233 P a at the equilibrium is much higher than in a thermal spectrum. Due to the lower capture cross sections of fission products in a fast spectrum, k ef f varies less in this later iteration of the EA design. In Ref. [START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF], k ef f = 0.98 is selected which means an energy gain of the order of 100 to 150 [START_REF] Rubbia | A high gain energy amplifier operated with fast neutrons[END_REF]. The accelerator itself is a 1 GeV cyclotron capable of delivering a beam intensity up to 20 mA.

Accelerator Transmutation of Waste (U.S. Department of Energy), 1999

In the ATW concept [START_REF]A Roadmap for Developing Accelerator Transmutation of Waste (ATW) Technology: A Report to Congress[END_REF] [START_REF] Yang | [END_REF], hazardous components of the nuclear waste, in particular minor actinides, will be destroyed to a high degree. The excess neutrons generated by the accelerator and subsequent fissions are used in a blanket-reflector containing the LLFP to be transmuted. Because of the excess neutrons, front end reprocessing is not required: the ATW system requires only the removal of the Zirconium cladding and the Uranium from nuclear reactor spent fuel. The separated Uranium is sent to a permanent storage while Pu and General layout of the Energy Amplifier complex [START_REF] Rubbia | Conceptual Design of a Fast Neutron Operated High Power Energy Amplifier[END_REF].

higher actinides are completely eliminated. Thus, the problem of siting a geological repository is not eliminated with an ATW, but it does provide a number of advantages:

• If the separations are sufficiently well done, the remnant waste could be stored in near-surface sites instead of geologic facilities.

• The Pu is never separated from the spent fuel which would eliminate the risk for proliferation issues. On the other hand, an ATW system operating at the same power level than LWRs could burn the waste of four LWRs, if its operating life were the same.

In conclusion, hybrid systems offer promising prospects due to a favourable neutronic balance and a large sub-criticality level that allow handling large amounts of MA. Confirmation that these properties are effectively exploitable has led many countries to develop ADSR as part of their program to close the fuel cycle: India uses ADS as a way to support the Thorium fuel cycle, the Japanese ADSR program is designed to transmute transuranics in the frame of the OMEGA project, the Chinese Academy of Science (CAS) is also developing an ADS program since 2011 and, in Europe, MYRRHA is the leading project for ADSR. At Brookhaven National Laboratory (BNL), the interest on this technology has lead to an inter-directorate collaboration between two departments, the Nuclear Science and Tecnology Department (NSTD) and the Collider Accelerator Department (C-AD). This collaboration took place over the period May 2013-May 2015 to discuss about the possibility of ADSR and was generally considered as a discussion forum producing many documents and proposals. A summary of this collaboration can be found in [START_REF] Méot | A SUMMARY OF AN INTER-DIRECTORATE C-AD / NSTD ACCELERATOR DRIVEN SUBCRITICAL REACTOR[END_REF].

Next, we will focus on an ADSR research facility: the Kyoto University Critical Assembly (KUCA), a feasibility evaluation of ADSR.

ADSR in KURRI

At Kyoto University Research Reactor Institute (KURRI), the world's first injection of spallation neutrons generated by high energy protons into a reactor core was accomplished on March 4, 2009 [START_REF] Pyeon | First Injection of Spallation Neutrons Generated by High-Energy Protons into the Kyoto University Critical Assembly[END_REF]. The KUCA is used as the sub-critical core for this experiment. An FFAG synchrotron accelerates the proton beam up to 100-150 MeV which is delivered to a tungsten (W) target (80 mm diameter × 10 mm thickness) located outside the critical assembly for safety reasons. In this configuration, a neutron guide is employed to lead the high energy neutrons into the center of the core. Several upgrades took place in the facility over the period 2009-2016 and the latest schematic of the installation is shown in Fig. 1.6. A 3D view of the KURRI-FFAG accelerator complex connected to KUCA is shown in Fig. 1.7.

An FFAG is a fixed field ring accelerator with strong focusing: the latter is achieved by alternation of the positive and negative bending field enabling strong focusing in both planes. Given the property of a fixed field, the FFAG combines properties of cyclotron and synchrotron. There are two different types of FFAGs: scaling FFAG where the number of betatron oscillations remains constant during the acceleration, and non-scaling FFAG where the number of betatron oscillations varies with the energy. The first type was proposed independently in the mid-1950s by T. Ohkawa in Japan, K. Symon in the United States and A. Kolomensky in Russia1 . Then, followed two periods of the FFAG development: the first, in Wisconsin, by MURA (Midwestern Universities Research Association) where three electron machines were constructed and tested in the period 1956-1967. Then followed a quiet period of nearly 30 years with limited interest in the FFAG technology due to several challenges (see Ref. [START_REF] Symon | MURA days[END_REF] for details). FFAGs were finally revived in the late 90s in Japan in the frame of the neutrino factory and medical applications R&D: a proton prototype was built and operated in 1999 by Mori's group at KEK (500 keV machine, 2.5 m diameter). In addition, more than 6 other scaling FFAGs were built in Japan in the subsequent years. Characterization and experiments are still carried on.

In 1997, Mills [START_REF] Mills | Proc. 4th Int. Conf. Physics Potential and Development of µ + µ -Colliders[END_REF] and Johnstone realized that the rapid crossing of the betatron resonances in a non-scaling FFAG mitigates their effect on beam emittance. The linear non-scaling FFAG concept arose and a demonstrator of this technology was successfully built and commissioned in Daresbury, UK, in 2010 [START_REF] Machida | Acceleration in the linear non-scaling fixed-field alternating-gradient accelerator EMMA[END_REF]. The activity in the US and Canada has mostly focused on the non-scaling concept [START_REF] Berg | Review of Current FFAG Lattice Studies in North America[END_REF]. However, the focus of this dissertation is on the scaling type in the frame of the KURRI FFAG collaboration in Japan: the main advantage of this concept is that the slow (following from the low energy gain per turn) crossing of betatron resonances is avoided, since the betatron tunes are fixed. Betatron resonances are a major concern in circular accelerators since the particle traverses many times the same magnet so that any error or field defect may be encountered many times. Therefore, having a fixed tune machine of the scaling type yields large dynamic aperture allowing large beams to be accelerated with minimum losses.

In order to meet the demand for high beam intensity, the FFAG complex in KURRI was upgraded to use a H -linac to inject the beam using charge exchange injection with negative hydrogen ion beams [START_REF] Mori | Present status of FFAG proton accelerators at KURRI[END_REF] in replacement of the earlier, double FFAG stage injector system (see Fig 1.8). The ADSR-FFAG accelerator is shown in Fig. 1.8: it consists of a scaling FFAG type that has 12 Defocusing-Focusing-Defocusing (DFD) dipole-triplet cells. At present, the number of protons per bunch is up to 3 × 10 9 corresponding to a 10 nA average current injected from the linac. The detailed parameters of the ADSR in- The present beam intensity is sufficient for the ADSR experiment (I < 1 nA to obtain P th,c < 10 W ). However the important beam losses need to be explained. The main reasons for the poor overall transmission of the machine (0.25% see table 1.3) are the dominant losses at injection as well as the crossing of a series of harmful betatron resonances. The latter represents a major problem that needs to be understood and eventually overcome: the major property of the scaling FFAG is that the number of betatron oscillations remains constant. However, as will be discussed in Chapter 3 of the present dissertation, certain types of field imperfections are unavoidable and must be accounted for. Given the complexity of the field profile in these machines (all multipoles are present in the field expansion ∝ R k ), a simulation campaign was established to benchmark different simulation codes as part of the intensity upgrade. In chapter 3, some aspects of the beam dynamics in fixed field ring accelerators are presented. Of particular interest is the case of the scaling FFAG that is discussed in detail with focus on the mathematical correctness of the model assumed. Analysis of the field profile and its derivatives revealed a major defect in the design of the magnets following which the average field index exhibits radial and azimuthal variations.

In order to explain the effects of these imperfections on the beam dynamics of the FFAG, analytical models were developed in chapter 4. It is shown that the Symon formula of the tunes as well as the scaling property of the orbits are violated in presence of azimuthal variations of the average field index of the magnets. Using the Bogoliubov method of averages, one derived a new formula for the tunes that does account for any field defect. This allowed to generalize the concept of a scaling FFAG to a non-scaling one which satisfies well defined tune variations. Besides, a general property was established which shows that, under some conditions, the number of betatron oscillations is a monotonic function of the energy. Exploiting the above results has lead us to establish a correction scheme to fix the tunes of the KURRI FFAG and avoid the problem of resonance crossing. The main idea of this original concept is to alternate the difference of the average field index of the focusing and defocusing magnets in the FFAG. This is the cornerstone of a new concept of a fixed tune non-scaling FFAG that we introduce and discuss in the last part of chapter 4. One important outcome of this study is to show that the two cardinal conditions of a scaling FFAG are sufficient but non-necessary conditions to obtain a fixed tune machine. In other words, the condition of scaling of the orbits at every azimuthal position is non-necessary and can be relaxed into an easier one to implement where the scaling of the orbits is achieved in an average sense: the idea being to implement alternating defects that alternate the monotonic behavior of the phase advance per cell. Study of the Dynamic Acceptance (DA) showed that for stronger alternating defects, the DA becomes smaller. This is essentially due to the case when the average field index of the D-magnet is larger than that of the F-magnet. A method was finally established for the DA calculation that allows to calculate the continuous data over the entire domain of stability by a simple barycentric interpolation through a discrete set of values.

One main outcome of this study was to show that trim coils must be considered in the preliminary design phase of the FFAG. This will enable a control of the field index of the magnets and therefore the tunes.

The main goal being to demonstrate the FFAG capability to produce high intensity beams, suitable for ADSR, one focused the study in chapter 5 on the space charge effects: an overview of the basics is first presented. Implementation of a space charge module in the tracking code ZGOUBI is then shown. This is based on the thin lens approximation. Some of the results are presented and application to the FFAG case is shown. As one discusses in the last part of chapter 5, the key property of the constant tune of the scaling FFAG is no longer valid in presence of space charge. Instead, one established a way to counteract the tune depression induced by the space charge effects by perturbing the average field indices of the magnets. This concept mainly works for uniform beam distributions where the tune spread is neglected. However, for other cases, the method is useful to move the core distribution in the tune diagram in a way to avoid harmful resonance crossings, as is the case in the KURRI FFAG. Another important result of these studies was to show that the tune spread in the new concept of the fixed tune non-scaling FFAG is similar to the ideal scaling FFAG case. This new concept is thus a promising candidate to demonstrate the high beam power capability in FFAG given its flexibility.

Chapter 2 presents the basic principles of ADSR. To begin with, one focused on the ADSR concept as a candidate to solve the problem of nuclear waste. In order to establish the accelerator requirements, one compared the performance of ADSR with other conventional critical reactors by means of the levelized cost of energy. A general comparison between the different accelerator technologies that can satisfy these requirements is also presented.

The wind scrambles and thunders over hills with a voice far below what we can hear. Whalesong, birdsongs boom and twitter. Sea, air, everything's a chaos of signals and even those we've named veer and fall in pieces under our neat labels. Waveshow to speak of the structure of waves when all disperses and there's nothing fixed to tell? -Philip Holmes, Background Noise Chapter 2

Accelerator-Reactor coupling

In this chapter, we first present some of the basic principles of the ADSR concept in view of better understanding the physics of sub-critical multiplying systems and how they differ from the conventional critical reactors. In particular, one focuses on the question of safety of ADSR. Unlike generally over-stated, ADSR are not intrinsically safer than conventional reactors: although super-criticality accidents are eliminated, overpower transients are a serious issue which can only be remediated for deeply sub-critical systems, therefore the economics of the system become questionable. In the second part of this chapter, we use the concept of the levelized cost of energy in order to compare an ADSR facility with other advanced reactor concepts. One assumes that the main mission of an ADSR is that of a waste burner, showing that dealing with the nuclear waste requires a supplemental cost of energy in the range 20-70 %. In the third and last part of this chapter, one defines the beam requirements of an accelerator suited for an ADSR and compares the different accelerator technologies that can fulfil this mission. The emphasis is made on the multi-accelerator scheme which resolves some of the key challenges and enhances the safety of the installation. A comparison of the costs of existing high power accelerators as well as the associated challenges is finally shown.

Basic principles of ADSR

Sub-criticality level

The level of sub-criticality, also referred to as the effective multiplication factor k ef f , is an average quantity that describes the way the neutrons multiply in the blanket. This is defined as the ratio of the number of neutrons produced by fission in one neutron generation to the number of neutrons lost through absorption or leakage in the preceding generation. k ef f writes in the following way:

k ef f = νΣ f Σ a P P + L (2.1)
where P is the production rate, L the leakage rate, Σ f the macroscopic fission cross section, Σ a the macroscopic absorption cross section and ν the number of neutrons released per fission.

In a conventional critical reactor, k ef f must remain close to 1 in order to ensure the steady state character of the reactor operation. This dictates the fuel type and content that can be handled.

Importance of the source neutrons

Although the external source of neutrons does not modify the k ef f of the subcritical blanket, it is important to realize that the multiplication of the spallation neutrons may be different from that of the neutrons that are born in the fuel elements from fission or other processes. Therefore, it is important to make this distinction between k ef f and k s where k s represents the source multiplication factor. This also suggests that all the neutrons are not valuable in the same way for transmutation, as expected. Several questions need to be answered:

• How is the kinetics of an ADSR different from the kinetics behavior of a critical reactor?

• How is the safety of an ADSR affected by the way the first generation neutrons (from spallation reaction) multiply in the core? Some of the content that we present next can be found in Ref [START_REF] Haj Tahar | BEAM STABILITY ANALYSIS FOR ADS[END_REF].

CHAPTER 2. ACCELERATOR-REACTOR COUPLING

Accelerator-Reactor coupling

If N 0 is the number of primary neutrons injected in the core per proton, from the spallation source (N 0 is an average quantity), and k i is the multiplication factor of the generation i, then the number of prompt fission neutrons after multiplication is:

N = N 0 (1 + k 1 + k 1 k 2 + k 1 k 2 k 3 + ... + k 1 k 2 ...k p + ...) -N 0 = N 0 k 1 (1 + k 2 + k 2 k 3 + k 2 k 3 k 4 + ... + k 2 k 3 ...k p + ...) (2.
2)

The first neutron generations have a multiplication factor k i that differs from k ef f , the multiplication factor of the blanket: the distinction is important since the spectrum of the fission neutrons is highly dependent on the spectrum of the neutrons from which they originated. Each fission produces ν ∼ 2.5 neutrons. Thus, the total number of fissions per proton is:

N f = N ν = N 0 ν k 1 S (2.3)
where

S = 1 + k 2 + k 2 k 3 + k 2 k 3 k 4 + ... + k 2 k 3 ...k p + ... (2.4)
The thermal power of the core is thus given by:

P th.c = E f (M eV )I(A) N 0 ν k 1 S + P dh (2.5)
where E f is the fission energy, I is the beam current and P dh is the power of the decay heat: the decay of the short-lived radioisotopes created in fission such as the fission products and the actinides continues even if the accelerator is shut-down, i.e (I = 0). This can represent up to 7% of the blanket thermal power, immediately after shutdown. For simplicity, we will assume here that P dh ∼ 0.

Another assumption in this analysis is that all the spallation neutrons injected in the core per second multiply in the same way. In other words, all series (defined by Eq. (2.4)) which originate from a primary proton converge in the same way: this is generally not the case and is mainly dependent on the time structure of the proton beam. Also, no reactivity feedback effects are considered here.

Safety analysis

A major advantage of ADS is its enhanced safety: if the sub-criticality level (or k ef f ) is far below 1 (this is generally chosen in order to accommodate for any possible positive reactivity insertion), super-criticality accidents are eliminated. However, it follows from Eq. (2.5) that overpower transients can be a serious issue: if the proton beam power (N 0 I) and/or the neutron source amplification (k 1 S) change, this implies that the thermal power of the core will follow as well. Therefore, it is crucial to understand the effect of the first generations of neutrons on the reactor power. Fig 2 .1 illustrates the importance of the first neutron generation k 1 which, for a given power, can lower the k ef f and thus increase the margin from criticality in an ADS. To have a better understanding of the source amplification, we will calculate the amplification term S defined in Eq. (2.4) above using a simplified model1 : Without any loss of generality, we define the sequence b i as follows:

k i+1 = (1 + b i )k ef f , 1 � i � p -1 k i = k ef f , i > p (2.6)
The sequence b i takes into account the importance of the first generation neutrons which have a multiplication factor different from k ef f up to the generation number p. p here can be made as big as possible. So, for the moment, the only assumption made is that it takes p generations before k i converges to k ef f . Here we assume that S = 1 1k ef f so that

P th ∝ k1 1 -k ef f It follows, S = 1 + k 2 + k 2 k 3 + k 2 k 3 k 4 + ... + k 2 k 3 ...k p + ... = 1 + h (1 + b 1 )k ef f + (1 + b 1 )(1 + b 2 )k 2 ef f + ... + (1 + b 1 )(1 + b 2 )...(1 + b p-1 )k p-1 ef f i + h (1 + b 1 )(1 + b 2 )...(1 + b p-1 )k p ef f + (1 + b 1 )(1 + b 2 )...(1 + b p-1 )k p+1 ef f + ... i = 1 + p-1 X i=1   i Y j=1 (1 + b j )   k i ef f +   p-1 Y j=1 (1 + b j )   k p ef f 1 -k ef f (2.7)
Then, by adding and subtracting,

S 0 = ∞ P i=0 k i ef f = 1 1 -k ef f = 1 + p-1 P i=1 k i ef f + k p ef f 1 -k ef f
, one obtains:

S = 1 1 -k ef f + p-1 X i=1   i Y j=1 (1 + b j ) -1   k i ef f +   p-1 Y j=1 (1 + b j ) -1   k p ef f 1 -k ef f = S 0 + δS (2.8)
The perturbation introduced is contained in the term δS: the first term of the latter represents the impact of the perturbation up to the generation number p. The last term accounts for the additional impact of the perturbation that is carried out by the later generations that converged to k ef f . Since the fission process is a chain reaction, the memory of the system is somehow conserved. Now, the thermal power of the core rewrites:

P th.c = E f (M eV )I(A) N 0 ν k 1 S 0 + E f (M eV )I(A) N 0 ν k 1 δS + P dh = P 0 + δP + P dh (2.9)
where P 0 represents the contribution of the neutrons which undergo fission reactions with a multiplication factor k i = k ef f while δP represents the excess or deficiency of the neutrons which undergo fission reactions with a multiplication factor k i 6 = k ef f . One should point out here that k ef f is an average quantity. In a critical reactor where there is no external source of neutrons, k ef f should remain close to 1 although the neutrons multiply in different ways depending on their spectrum and location in the core. In an ADS where the source is the driving term, the fact that the neutrons multiply in different ways makes the analysis far more complicated and intricate: the history of the spallation neutrons in the core which are the source of all forthcoming neutrons has an impact such that, if the multiplication factor k 1 of the first generation neutrons multiplies by 2, the thermal power of the core will double as well. This shows the importance of tailoring the spectrum of the target to maximize its efficiency. Although the spallation neutrons represent a small fraction of the total neutrons in the core, it allows more degrees of freedom to play with, that cannot be otherwise allowed in a critical reactor, since it would jeopardize its safety. Now, let us assume that b j = 1 j and inject it into Eq. (2.8). This form of b j was chosen in order to fit an example of the evolution of the multiplication factor as a function of the generation number (Fig. 1 in [START_REF] Nifenecker | Basics of accelerator driven subcritical reactors[END_REF]). It follows: as a function of the sub-criticality level: one can observe that, for deeply sub-critical systems, the effect is less severe. Also, note that the most important effect of the perturbation is carried out by the later neutron generations that converged to k ef f (solid vs dashed lines). In order to understand this result, the idea is to study the convergence rate of the series

S = 1 1 -k ef f + p-1 X i=1 i.k i ef f + (p -1) k p ef f 1 -k ef f ; p � 2 (2.
S 0 = ∞ P i=0 k i ef f = 1/(1 -k ef f ).

Convergence rate and the concept of neutron multiplication time

The convergence rate of the order p is defined as the relative error of the sum S 0 by estimating only the first (p-1) terms of this series: the error is down to 0.6% after 100 neutron generations.

ǫ = S p S 0 = 1 1 -k ef f - 1 -k p ef f 1 -k ef f ! × (1 -k ef f ) = k p ef f (
As the multiplication factor k ef f becomes higher and higher (k ef f ∼ 1), more time is needed for the system to converge to the final state level where S 0 = 1/(1 -k ef f ). Therefore the fission process continues for a long time after the injection of the spallation neutrons because the contribution of forthcoming generations becomes large. Thus, perturbing the first neutron generations will have an important impact mainly for slowlyconverging systems, i.e for nearly critical systems as shown in Fig. 2.2.

Based on the results above, we define the multiplication time T m as the time it takes for the spallation neutrons injected into the core to multiply and reach the amplification factor S 0 within a specified accuracy. It is a measure of the duration of the lasting fission process induced by the injection of the spallation neutrons into the core. Then, according to Eq. (2.11), it results

T m = ln(ǫ) ln(k ef f ) T g ∝ 1 | ln(k ef f )| ; k ef f � 1 (2.12)
where T g is the mean generation time, i.e the time between two consecutive neutron generations in the core. This result illustrates the importance of k ef f on the safety of an ADS: if k ef f is very close to 1, the reliability of the accelerator is no longer problematic since the fission process continues for a long time after the beam is off. However, this poses a serious safety issue and cannot be tolerated.

Point kinetics equations for ADSR

We write the coupled kinetics equations for a point reactor model with one group of delayed 2 neutrons:

         dn(t) dt = ρ -β ∧ n + λ.C + S dC(t) dt = β.n ∧ -λ.C (2.13)
where n is the neutron population, ρ is the reactivity (ρ = (k ef f -1)/k ef f ), β the fraction of delayed neutrons (assuming only one group), ∧ the generation time (∧ = l/k ef f ), l the neutron lifetime (or average life expectancy of neutrons in the reactor), λ the decay constant for precursor decay, C the density of the precursors and S is the external source of spallation neutrons. Temperature feedback effects are ignored. More information about the latter can be found in [START_REF] Nifenecker | Basics of accelerator driven subcritical reactors[END_REF].

Transmutation characteristics

The two major indicators of the transmutation performance are the transmutation speed and the transmutation yield. In a simplified model, suppose N 0 nuclei are irradiated with a time independent neutron flux Φ. Then the evolution of these nuclei includes two terms:

1. A term of destruction due to natural decay (λ nat ) and to the transmutation through neutron absorption with an absorption cross section σ.

2.

A term of production p supposed to be constant here.

The evolution of these nuclei is thus given by:

dN dt = -(λ nat + σ.Φ) N + p (2.14)
Solving this equation and assuming λ nat ≪ σ.Φ yields:

N (t) = N ∞ + (N 0 -N ∞ ). exp(-σΦt) ; N ∞ = p σΦ (2.15)
which shows that the initial inventory will evolve to an equilibrium with a time constant τ = 1/σΦ. Three different regimes can be distinguished:

1. N 0 > N ∞ : burner

2 delayed here refers to the fact that some of the neutrons are produced by daughter nuclides of precursors coming out of fission

2. N 0 = N ∞ : equilibrium 3. N 0 < N ∞ : breeder
Now, the transmutation performance can be evaluated via the transmutation speed or its yield (both linked): Transmutation speed : as deduced from Eq. (2.15), one can determine the time required to eliminate half the initial inventory (in the absence of any other competing reaction, i.e. p = 0): T 1/2 = ln 2/(σ.Φ) and is of the order of several years. Transmutation yield : this is defined as the ratio of the quantity transmuted to the initial one:

Y = ΔN N 0 = 1 - N ∞ N 0 (1 -exp(-σ.Φ.Δt)) (2.16)
In reality the amount of nuclides in the core evolves according to a coupled set of equations where all isotopes j 6 = i can decay into isotope i according to the Bateman equation:

dN i dt = -(λ i + σ i Φ) N i + X j6 =i (λ j R j→i + σ i P j→i Φ) N j (2.17)
where R j→i and P j→i are the branching ratio of the natural neutronic reaction respectively.

Comparison of ADSR with other advanced reactor concepts

The preceding section discusses the main aspects of ADSR and focuses on the differences with conventional critical reactors. Essentially, it shows that, although super-criticality accidents are eliminated, overpower transient effects can be a serious issue due to the coupling of the reactor core with the external source of neutrons. Aside from that, an ADSR facility is subject to the same safety requirements as any other conventional reactor facility. The question that naturally arises is the following: do the benefits of deploying an ADSR facility outweigh the challenges, and how can we evaluate its competitiveness with other advanced reactor concepts? Although there is no straightforward answer to this question, one needs to identify the main benefit of adding an expensive accelerator to a nuclear reactor facility. Following an inter-directorate C-AD/NSTD (Collider Accelerator Department/Nuclear Science and Technology Department) collaboration over the period May 2013-May 2015 [START_REF] Méot | A SUMMARY OF AN INTER-DIRECTORATE C-AD / NSTD ACCELERATOR DRIVEN SUBCRITICAL REACTOR[END_REF], it was accentuated that the main identifiable benefit of an ADSR is in a nuclear waste transmutation mission. In other words, remediating the problem of nuclear waste can be a sound motivation for ADSR.

From that perspective, one decided to compare the costs of energy produced in different technologies by internalising all the possible external costs, including the costs of the waste disposal.

Externalities

As defined in [START_REF]External costs, Research results on socio-environmental damages due to electricity and transport[END_REF] "An external cost, also known as an externality, arises when the social or economic activities of one group of persons have an impact on another group and when that impact is not fully accounted, or compensated for, by the first group." Air pollution is a typical example of this. In the matter of nuclear energy, external costs stem from the environmental impacts of waste or emissions.

In the context of sustainable energy goals, internalising the external costs is crucial. For instance, in the absence of taxation on carbon emissions, there is no incentive for the energy producers to switch to other sources with lower emissions. Therefore, the comparison of the different energy sources should be made using levelized costs. This means that the average costs of the energy generation over the lifetime of the project, including operation, fuel costs, decommissioning as well as waste disposal have to be included.

The Levelized Cost Of Energy (LCOE)

Cost reduction is the common goal of all power plants: for instance, in nuclear energy, this can be achieved by improving the efficiency of the plants, their reliability as well as their lifetime. Among the many questions to be answered is whether the advances of the fuel cycles contribute to reducing the costs of nuclear energy. Thus, a parameter that is commonly used to evaluate the economics of the back-end options of different types of power plants is the levelised cost of electricity generation (LCOE) [START_REF]The Economics of the Back End of the Nuclear Fuel Cycle[END_REF]. This is a particularly useful tool to compare the cost of electricity generation with different technologies. The LCOE ($/MWh) is defined as follows:

LCOE = P t Investment t + O&M t + Fuel t + Carbon t + Decommissioning t (1 + r) t P t Electricity t (1 + r) t (2.18)
where the subscript t denotes the year in which the electricity is produced and the expenses are made. The rest is defined in the following way:

• Electricity t : Amount of electricity produced in the year t;

• r: Annual discount rate;

• Investment t : Investment cost in the power plant in the year t;

• O&M t : Operations and maintenance cost in the year t;

• Fuel t : Fuel costs in year t (both front-end and back-end);

• Carbon t : Carbon cost in the year t;

• Decommissionning t : Decommissionning cost in the year t.

The denominator represents the discounted cash flow from selling the electricity at a constant price (assumption).

Fig. 2.4 shows a comparison of the LCOE in four countries: Although these studies have been carried out in different countries, and therefore reflect their national policies and are specific to their methods and strategies employed, it can be seen that the costs of nuclear energy are substantially cheaper than the renewables in four countries. At 7% discount rate, it is comparable with gas and cheaper than coal. System costs represent the total costs above plant-level costs to supply electricity at a given load and given level of security of supply [START_REF]The Economics of Nuclear Power[END_REF]. A discount rate of 7% is used.

Costs of Spent Nuclear Fuel (SNF) management

The parameter that should be employed to quantify the cost of the SNF management of an ADSR must include the cost of the electricity generated with the facility as well as the costs of the waste disposal. Besides, the back end component of the LCOE can be defined: in that case the above definition is summed over all back-end facilities. The latter are determined according to the strategy selected for the waste disposal: for a direct disposal, interim storage facilities and a deep geological repository are required. This is defined as the Once Through (or open) fuel Cycle (OTC). Fig. 2.5 shows a pie chart of the LCOE for the OTC which breaks down into its four main components: the investment (or capital cost), the reactor (non-fuel) Operating and Maintenance (O&M) costs, as well as the Back-end fuel cycle and Front-end fuel cycle costs: the front-end breaks down into the cost of raw uranium and the cost of its manufacture, i.e conversion, enrichment and fabrication, while the back-end breaks down into the cost of interim above ground storage of the spent fuel and the cost of a deep geological repository for ultimate disposal. The nuclear fuel cycle (front-end+back-end) accounts for ∼ 13 % of the total power generation costs which shows that the cost of nuclear energy is largely driven by the capital cost. For an advanced fuel cycle, several plants are required for reprocessing, fuel fabrication, waste conditioning as well as an underground repository.

The LCOE for ADSR

Several studies to assess the different waste management strategies lead to a similar conclusion [START_REF]The Future of Nuclear Power[END_REF]: "technical improvements to the waste management strategies in the once-through fuel cycle are potentially avail- 6 % discount rate is assumed. The details of the model assumed can be found in [START_REF]The Future of the Nuclear Fuel Cycle[END_REF].

able that could yield benefits at least as large as those claimed for advanced fuel cycles featuring waste partitioning and transmutation, and with fewer short-term risks. The most that can reasonably be expected of partitioning and transmutation schemes is to reduce the inventory of actinides in geologic repositories by perhaps two orders of magnitude."

In the following analysis, we will not compare the different fuel cycle options. We refer the interested reader to the 2014 comprehensive nuclear fuel cycle Evaluation & Screening report [START_REF] Wigeland | Nuclear fuel cycle evaluation and screening-Final report: Appendix B, Comprehensive set of fuel cycle options[END_REF] that was chartered by the US Department of Energy. Our focus will be on the ADSR option in order to evaluate the impact of adding an accelerator on the LCOE. From Eq. (2.18), one can identify the main parameters that change in comparison with a fast critical reactor installation:

• The Investment increases by x % depending on the accelerator type (cyclotron, linac, FFAG,...), size (beam final energy) as well as the beam intensity required and therefore k ef f (lim k ef f →1 x = 0): for instance, the cost of the 1.4 MW SNS linac, suitable for ADSR pre-industrial demonstrator represents roughly 10 % increase in the investment cost.

• The Fuel cost is different given that the fuel needs to be reprocessed. The front-end fuel cycle is the main reason of the high cost of advanced fuel cycle strategies given that nuclear reprocessing is not a competitive process due to the low price of natural uranium (scarcity of natural uranium is not yet a problem).

In the best case scenario, one can assume that the most offending radioisotopes in the spent fuel are eliminated yielding a reduction in the volume of the waste disposal facility. This is reflected as a potential saving in the cost of the back-end fuel.

• The Operation & Maintenance cost: these are generally assumed to be 20 % larger for a fast critical reactor as compared to a light water reactor [START_REF]The Future of the Nuclear Fuel Cycle[END_REF]. One further assumes that adding an accelerator will raise the O&M costs by roughly x %.

• The amount of electricity produced with an ADSR is lower than with a critical reactor that has the same thermal power output. This is due to the performance of the accelerator which does impact two key parameters of the ADSR installation: the amount of electricity drained by the accelerator lowers the overall thermal efficiency of the installation while the accelerator reliability lowers the capacity factor of the plant. In what follows we evaluate these two parameters.

Energy efficiency of ADSR

The goal of this section is to evaluate the corrected thermal efficiency of an ADSR. Some of the results presented next can be found in Refs [26][37].

CHAPTER 2. ACCELERATOR-REACTOR COUPLING From Eq. (2.5), it follows that:

P acc,grid P th,c = E p E f ν N 0 1 k 1 S 1 η acc (2.19)
where P acc,grid is the power drained from the grid by the accelerator and η acc is the wall-plug efficiency of the accelerator (one assume that P dh ∼ 0). Next, we compute the corrected thermal efficiency of an ADSR: this is essentially the net useful power output (a fraction of the electrical power is fed back to power the accelerator) divided by the total thermal power input:

η th = P el,c -P acc,grid P th,c = η th0 - E p E f ν N 0 1 k 1 S 1 η acc = η th0 -A 1 k 1 S 1 η acc (2.20)
where η th0 is the uncorrected thermal efficiency of the installation.

For a 1 GeV proton beam impinging on a lead target, N 0 ≈ 20 so that A ≈ 0.6.

The first spallation neutrons are the source of all forthcoming neutrons. Therefore, tailoring their spectrum in order to increase their importance vis-à-vis the fission process is important. If well mastered, this is very useful as will be discussed later on in this chapter. First, we assume that k i = k ef f for all i. Then, Eq.(2.20) becomes:

η th = η th0 -A 1 -k ef f k ef f 1 η acc (2.21)
where k ef f is the effective multiplication factor. In most ADS proposals k ef f is chosen in the range [0.95:0.98] and has to accommodate for any possible accident, i.e for any possible positive reactivity insertion during the operation of the reactor, including the fuel loading stage. Fig. 2.6 below shows a contour plot of the corrected thermal efficiency of an ADSR assuming an uncorrected reactor efficiency, η th0 ≈ 40%, consistent with a fast reactor concept. When k ef f is close to 1, the impact of η acc is negligible: the reactor is very near the critical state, and the power drained by the accelerator to sustain the core thermal power is negligible. The same result is shown in Fig. 2.7. Comments:

• To obtain the system's average efficiency, one normally averages the output over a period of time and divide it by the average input over the same period of time. However, this is not entirely true for a complex system such as ADS, especially when system failure occurs. Beam trips longer than few seconds induce a reactor shut-down which requires several hours to restart the reactor during which the facility is still running in the background.

• Reliability conflicts with efficiency: for instance, ADS reliability improves through redundancy. The main challenge is to improve the integrated efficiency over a long period of time with a fault-tolerant design.

• The previous calculation assumes that k i = k ef f for all i. In reality, given that the spallation neutron spectrum is decoupled from the fission neutron spectrum, different neutron generations have different multiplication factors. The history of the spallation neutrons in the core is the most important parameter for an ADS, because these are the source of all forthcoming neutrons. This can be turned into advantage provided one can tailor their spectrum. For instance, suppose that k i = k ef f for all i > 1. Then, from Eq (2.20),

η th = η th0 -A 1 -k ef f k 1 η acc (2.22) Fig 2.
8 illustrates the impact of the first neutron generation on the overall efficiency of an ADS. It can be seen that k 1 plays the same role as η acc . In other words, improving the importance of the first generation neutrons (by optimizing the design of the target) has the same energy efficiency impact as improving the wall-plug efficiency of the accelerator.

Accelerator efficiency

One can identify three main steps of the energy transformation from input to output:

1. AC current is the most efficient way to deliver the electrical power. However, most of the electronic devices require DC power to function. Thus the 1st transformation requires an AC to DC converter. We refer to the conversion efficiency as η DC . This is typically in the range 80 -90%. For PSI and SNS, η DC = 90%.

2. DC to RF power conversion. This transformation takes place in a RF amplification system such as triodes, RF tubes or klystrons. The most commonly used for linacs are klystrons for frequencies above about 300 MHz [START_REF] Wangler | RF linear accelerators[END_REF]. We refer to the conversion efficiency as η RF . Klystrons typically achieve 40% for pulsed operation and almost 60% for CW operation, suitable for ADS. Only a few high power klystrons offer 65%+ efficiency [START_REF] Lingwood | Klystrons", Proton Driver Efficiency Workshop[END_REF].

For the PSI main ring cyclotron, η RF ≈ 64% while for the SNS superconducting linac, η RF ≈ 30%. At this stage, one could also include the magnet power conversion efficiency. However, it is best to leave this for later. 3. RF power to beam power. We refer to its efficiency as η beam . In a superconducting cavity, nearly all of the RF power goes to the beam. A parameter that measures the effectiveness per unit power loss for delivering energy to the particle is the shunt impedance, Z = V 2 0 /P , where V 0 is the peak RF voltage in the gap and P is the average RF power loss, so that η beam ∝ Z. Superconducting technology is more suitable for a high intensity accelerator. However, the situation is non-trivial in the presence of beam loading. For a superconducting cavity, the efficiency of the cryogenic system needs to be included as well.

In conclusion, the wall-plug efficiency of an accelerator is given by:

η acc = η DC η RF η beam η other (2.23)
where η other is the additional contribution of the other sources of power consumption aside from the RF cavities, and is defined by:

η other = AC AC + magnet + cryo + ... (2.24)
For instance, an additional source of power consumption for the PSI cyclotron comes from the magnets, while for the SNS superconducting linac, the cryogenic system lowers the overall efficiency of the linac. Table 2.1 shows the current data for PSI and SNS which are the world leading accelerators in terms of average beam power produced. These data are not intended to make any comparison between the different technologies, but solely used here to identify the major sources of power losses that can be further reduced. For [START_REF] Seidel | Cyclotrons for high-intensity beams[END_REF] and for the SNS linac (P beam ≈ 1M W ) [START_REF] Galambos | [END_REF], as of 2010. The two transformations where there is room for improvement are the DC to RF power and the RF to beam power transformations.

accelerators (the beam requirements will be established in the next section), it can be shown that, with the current state of the art, the projected accelerator efficiency lies in the range [20% : 30%]. This represents [25% : 40%] downgrade of the uncorrected thermal efficiency of the power plant (assuming k ef f = 0.95).

Capacity factor of ADSR

In order to have a more accurate picture of the corrected thermal efficiency of an ADSR it is important to include the beam trip rate of the accelerator in the analysis. For that reason, one introduces the operational efficiency η oe of an ADSR facility: this is defined as the ratio between the total electricity produced by the plant during a period of time compared to the total potential electricity that could have been produced if the plant operated without an accelerator, i.e without any beam trips. It results that η oe can be expressed in the following way:

η oe = C C 0 (2.25)
where C and C 0 are the capacity factors of the power plant with and without the accelerator over the same period of time: the capacity factor is defined as the ratio of the plant's actual output over a period of time to its potential output if it operated continuously over the same period of time. As shown in Table 2 In a simplified model, let's assume that the accelerator has N well separated beam trips during the operation period T long enough to induce a reactor shut-down lasting Δt bt . Then the corrected capacity factor of the plant writes:

C C 0 = 1 - N × Δt bt T (2.26)
Δt bt is the time required to restore the accelerator beam and reach the nominal power of the reactor after shut-down: this is a multi-dimensional problem, therefore very difficult to estimate. However, it is vital for ADSR. The MYRRHA concept estimates that Δt bt ∼ 2h [START_REF] Abderrahim | [END_REF].

It is clear that, due to beam trips, η oe ≤ 1. From all the above, it results that the net useful electricity of an ADSR facility produced over a period of time t is given by:

P el = η th η oe P th,c = C C 0 η th η th0 P el0 (2.27)
where P el0 is the uncorrected electrical power of the reactor, i.e without taking into account the accelerator.

Based on all the above, one can calculate the LCOE of an ADSR in the following way:

LCOE LCOE min = C 0 C × η th0 η th (2.28)
where LCOE min is the LCOE equivalent to a fast critical reactor. Fig. 2.9 illustrates the impact of deploying an accelerator driven system on the LCOE where one neglected the effect of beam trips: if k ef f ∈ [0.95 : 0.97], then the achievable supplemental cost of ADSR lies in the range [20% : 70%]. For this reason, an ADSR facility whose unique mission is to produce energy cannot compete with any other advanced reactor concept. However, if one were to deal with the nuclear waste and reduce its radio-toxicity to low levels, then the above represents the price to pay. Note that the system costs, i.e. the total costs above the plant level to supply the electricity at a given load and given levels of security of supply [START_REF]The Economics of Nuclear Power[END_REF], are not considered here. These costs will depend on the reliability of the accelerator. Although the cost of deploying an ADSR system is high compared to other advanced reactor concepts, in order to evaluate its economic impact on the entire nuclear reactor fleet, the idea is to determine how many ADSR are needed to burn all the nuclear waste produced? Based on the industrial model of the EUROTRANS european research program, which assumes a maximum thermal power of the core P th,c ≈ 400M W th , a technical and economic assessment showed that 18 ADSR will be necessary to obtain a transmutation capability equivalent to a fast nuclear reactor fleet of 60 GWe (∼ nuclear power capacity in France). Note that the transmutation rate is limited by the thermal power of the core.

Conclusion

As mentioned in [START_REF]The Future of the Nuclear Fuel Cycle[END_REF], "the most important fact to keep in mind in considering any estimate of the cost of alternative fuel cycles is the high degree of uncertainty about key components of each cycle". This is just to remind that the technology of nuclear waste disposal is highly immature since no facility has already been built to meet all the requirements. However, in the best case scenario, deploying the ADSR technology will increase the LCOE by at least 20 %. Although the cost of the fuel cycle is a very important factor, one has to keep in mind that some other factors include the safety of the installation as well as the political facets (proliferation risks, etc).

Beam requirements and comparison of accelerator technologies

In order to assess the beam requirements of the accelerator, one has to look at the performance of the spallation target and its compatibility with the steady state character of the reactor operation. Monte Carlo simulations of the spallation process are thus necessary in order to determine the multiplicity of the target and thus the cost of producing the external neutrons. Some of the content that we discuss next can be found in Refs. [START_REF] Heidet | Accelerator-Reactor Coupling for Energy Production in Advanced Nuclear Fuel Cycles[END_REF] [38].

Proton energy

The neutron multiplicity and spectrum depend on the incident proton energy, the target nuclei as well as the engineering design: An interesting quantity, shown in Fig. 2.10, is the neutron cost, defined as the number of produced neutrons per incident particle normalized to the unit beam energy. It shows that the neutron production is less efficient above 1.3 GeV in the case of a lead target. There are two main reasons for this: first, pion production dominates at the cost of neutron production above 1 GeV. The second reason is the proton range which is typically limited to the length of the target, i.e. 55 cm at 1 GeV. This implies that the target is inefficiently used for higher energy particles. The target radius plays the same role as illustrated in Fig. 2.11: at lower radii, the number of escaped neutrons increases more steeply with the increasing energy. A saturation effect is observed at higher radii where the particle range is well contained within the target. Thus, a target having dimensions beyond that of the proton range does not offer much gain. However, as explained earlier, all the neutrons that escape from the target are not valuable in the same way for transmutation. Therefore, one has to look at the spectrum of the spallation target and how it compares to that of the fission process. A comparison between the spallation spectrum and the blanket spectrum is illustrated in Fig. 2.12: Although fast neutrons are more favorable for the transmutation process due to a higher fission/absorption cross section ratio, the tail of very high energy cascade neutrons (E n > 20M eV ) is very difficult to shield. These neutrons will likely pass through the bottom shield of the core and end up in the surrounding environment [START_REF] Meot | [END_REF]. This tail extends up to the incoming proton energy and can represent a non-negligible fraction of the entire spectrum. It results that the optimum proton energy lies in the range [0.5:1] GeV which is a trade-off between cost optimization, target heating as well as the shielding issue.

Beam intensity

Once the thermal power of the core is fixed, the required beam intensity can be determined by means of Eq. (2.5). However, it is important to note that there are two operational modes of ADSR: one in which the nominal power of the accelerator is fixed (which is optimized for better efficiency). When the blanket is depleted, this implies that the thermal power of the core decreases. The other operational mode consists in adjusting the accelerator power to follow the variation of the k e f f of the blanket: for instance, when the blanket is composed of pure thorium at the start of the fuel cycle, the build-up of fissile material (which implies an increase of k ef f ) does require a decrease of the accelerator power to avoid any criticality accident. If we assume that P acc is fixed, then one obtains:

I(A) = P th,c (M W th ) E f (M eV ) N 0 ν k ef f 1 -k ef f (2.29)
Assuming k ef f = 0.95 to avoid any criticality accident, N 0 = 20 which corresponds to 1 GeV proton beam impinging on a lead target and ν = 2.5, this yields: I(mA) = 3.29 × 10 -2 × P th,c (M W ). As shown in Fig. 2.13, the average beam current required can range from 10 mA for a modular reactor to 100 mA for an industrial 1 GWe reactor.

Beam time structure

The beam time structure is a crucial parameter to assess in order to maintain the steady state character of the reactor operation. Studying the impact of the time structure of the proton beam on the kinetics of the reactor allows to identify the required beam parameters as well as understand the start-up process of the reactor. Besides, this is important to understand the impact of the accelerator reliability which is identified in [START_REF] Abderrahim | Accelerator and target technology for accelerator-driven transmutation and energy production[END_REF] as the major challenge of ADSR: "ten to one-hundred fold improvement in long-duration beam trip rates relative to those achieved in routine operation of existing high power proton accelerators is necessary to meet industrial-scale ADS application requirements". In our analysis, we use a simplified model to evaluate the kinetics evolution of the neutronics of the ADSR [START_REF] Tahar | Effect OF THE BEAM TIME STRUCTURE ON THE NEUTRONICS OF AN ACCELERATOR DRIVEN SUBCRITICAL REACTOR[END_REF]: To begin with, we assume that p + 1 neutron bunches are injected into the core when the accelerator beam is turned on as illustrated in Fig. by t b . Besides, one assumes that the time structure of the neutron bunches is identical to that of the proton beam. Thus, S, the external source of spallation neutrons, writes in the following way:

S(t) = p X i=0 N s δ (t -it b ) (2.30)
where N s is the total number of spallation neutrons per bunch. Now, we introduce t r , the time elapsed since the arrival of the last bunch of neutrons, indexed p. Then, one can write:

     t = pt b + t r 0 � t r < t b 0 < t b =⇒ p = E t t b (2.31)
where E is the integer part function. We then solve the coupled point kinetics equations (2.13) for the last bunch of neutrons by means of Laplace transforms. This yields, after some simplifications:

n(t r ) = A 1 exp - t r τ 1 + A 2 exp - t r τ 2 (2.32)
where referred to as the prompt drop (green) characterized by the time constant τ 1 . The second term (blue) with much longer time characteristics τ 2 represents the contribution of the delayed neutrons. This is an illustrative example and these quantities depend on the fuel composition, the neutron spectrum as well as the sub-criticality level. If we neglect the contribution of the delayed neutrons, i.e. β ≈ 0 (this is what justifies the need for the external source of neutrons), then one calculates the total number of neutrons available in the core at any given time:

τ 1 = ∧ β -ρ ≈ l 1 -k ; τ 2 = ρ -β λρ A 1 = N s ; A 2 = -N s β ρ -β
N (t r ) N s = exp(- t r τ 1 ) + exp(- t r + t b τ 1 ) + ... + exp(- t r + pt b τ 1 ) = exp - t r τ 1 × 1 -exp -(p + 1) t b τ 1 1 -exp - t b τ 1 (2.33)
which shows that the impact of any beam interruption is an exponential decay with a time constant τ 1 . Eq.(2.33) can also be expressed as an explicit function of t,

N (t) N s = exp    - t -E( t t b )t b τ 1    × 1 -exp - E( t t b ) + 1 t b τ 1 1 -exp(- t b τ 1 ) (2.34)

Repetition rate

One particularly important parameter that intervenes in the above equation is the bunch repetition rate (≡ t b ). In order to evaluate its impact on the neutronics, we define the equivalence of two beams in terms of average beam current in the following way:

( N s t b is equivalent to ( f N s f t b (2.35)
where f is a time structure scale factor If we apply this equivalence in Eq. (2.34) above (by replacing (N s , t b ) by (f N s , f t b ), one obtains:

N (t) N s = f exp   - t -E( t f t b )f t b τ 1    × 1 -exp - E( t f t b ) + 1 f t b τ 1 1 -exp(- f t b τ 1 )
(2.36) Fig. 2.16 shows the evolution of the neutron population N (t)/N s for various values of the scale factor f : as can be seen, lowering f and thus increasing the bunch repetition rate is better for ADSR since it reduces the overall variation of the neutron population. The average state, when convergence is reached, is independent of the beam time structure. Only the average beam power of the accelerator determines the average thermal power of the blanket. To show this, let's calculate the average as well as the amplitude of variation of the neutron population. First, one introduces (2.35) into Eq. (2.33). In the limit where t → ∞, i.e when p → ∞, the steady state is reached, and N (t) becomes periodic in time with period f t b . This yields:

✵�✵✵✵✁ ✵�✵✵✵✂ ✵�✵✵✵✄ ✵�✵✵✵☎ ✵�✵✵✵✆ ❚✝✞✟✭s✮ ✁✶�✵ ✁✶�✆ ✂✵�✵ ✂✵�✆ ◆✟✠✡☛☞✌ ❝☞✠✌✡s ❢ ❂ ✍ ❢ ❂ ✎✏✑ ❢ ❂ ✎✏✒✑
hN (t r )i N s = 1 f t b N s Z f t b 0 f N s exp - t τ 1 1 -exp - f t b τ 1 dt = τ 1 t b = l t b 1 1 -k (2.37) N max N s = f 1 -exp(- f t b τ 1 ) (2.38) N min N s = f exp(- f t b τ 1 ) 1 -exp(- f t b τ 1 ) (2.39)
As shown in Fig. 2.17, the best choice for ADSR is a CW beam with a high repetition rate. As a general conclusion, the proton bunch separation time t b must be lower than the period of the rapidly decaying transient τ 1 . Table 2.3 shows the required repetition rate of the accelerator assuming the same number of protons per bunch (N ppb ∼ 10 10 ) and an overall variation of the neutron population limited to 1%. The condition writes in the following way: Reactor type variation due to a short interruption of the proton beam (� 1 ms) will be small. A full thermo-mechanical analysis of the design is thus needed to establish the number of beam interruptions allowed as well as their duration. Thus, the beam repetition rate can be reduced.

N max -hN i hN i ≤ ǫ (2.
P th,c [MW] I [mA] N ppb F [MHz] MYRRHA 80

Beam distribution

Another parameter of interest for ADSR is the beam distribution on the target: in order to calculate the impact of the latter, we use the same hypothesis as in the previous section, i.e. we neglect the contribution of the delayed neutrons. Then,

N (t) = n s (t n ) × (t n -t n-1 ) × exp(- t -t n τ 1 ) + n s (t n-1 ) × (t n-1 -t n-2 ) × exp(- t -t n-1 τ 1 ) + ... + n s (t 1 ) × (t 1 -t 0 ) × exp(- t -t 1 τ 1 ) (2.41)
where n s (t i ) is the rate of spallation neutrons ( neutrons s ) between t i and t i-1 . In the limit where t i -t i-1 → 0, this yields 4 ,

N (t) = Z tn t 0 n s (τ ) exp(- t -τ τ 1 )dτ ; t � t n � t 0 (2.42)
where t 0 (resp. t n ) represents the time at which the first (resp. last) bunch of spallation neutrons arrived to the core.

Since n s is a periodic function of time, it can be decomposed into a Fourier series:

n s (τ ) = a 0 + +∞ X n=1 a n cos[nωτ ] + b n sin[nωτ ] (2.43)
Injecting (2.43) into (2.42) and integrating it yields the expression of the evolution of N with any arbitrary periodic function of n s . Note that the speed of the exponential decay is insensitive to the time structure of the proton beam and depends only on the intrinsic properties of the reactor, i.e. the sub-criticality level and the neutron lifetime. Thus the beam time structure does not play any role in remediating the problem of beam trips. However, the beam radial distribution can increase/lower the target multiplicity as shown in Fig. 2.18: increasing the beam size decreases the neutron multiplicity for the present configuration: this can be explained by recalling the expression of the probability that a non-elastic nuclear collision occurs in the target:

P n = 1 -exp - R λ (2.44)
where R is the range of the proton (g/cm 2 ) and λ is the non-elastic nuclear collision mean free path. The latter is almost independent on the energy of the proton beam (above 100 MeV) [START_REF] Watanabe | Neutronics of pulsed spallation neutron sources[END_REF] so that the only parameter that plays is the range of the incoming particle: if the target is under-sized, i.e. the size of the target is less than the range of the incident proton beam, then, increasing the beam size will have an effect such as more outer particles will escape from the system before any nuclear collision process occurs.

Comparison of different accelerator technologies

There is no general criterion on how to provide the external source of neutrons. With the current state of the art, a spallation source is the most viable option to provide the high neutron flux to the core and minimize its cost. However, many requirements make the accelerator design challenging. In order to evaluate the best accelerator technology for ADSR, one needs to summarize the list of parameters/requirements to fulfil its mission. Analysis of multi-beam accelerator driven system [START_REF] Ludewig | Study of Multi-Beam Accelerator Driven Thorium Reactor[END_REF] showed the merit of using a multi-beam target system:

Beam property

1. "The fractional change in energy deposited in the fuel following a beam trip is significantly reduced for the cases using multi-beam target arrangements." This helps reduce the thermal stress and fatigue failures which relaxes the requirement of the number of allowed beam trips.

"

The production of spallation neutrons at multiple sites within the blanket reduces the peak/average energy deposition. This ensure that the transmutation of the fissile and fertile materials takes place more uniformly."

Another important merit of a multi-beam configuration is from the viewpoint of safety: if the accelerator beam power is fractioned into N beamlets, then the impact of overpower transient effects (see Eq. (2.9)) will be reduced since the likelihood that different neutrons sources (at different locations in the core) undergo the same perturbation is diminished.

Since the early 1990s, the obvious candidates that were considered to address the mission of ADSR are cyclotrons and linacs. However, due to the recent revival of interest in using the FFAG technology, a question that naturally occurs is what sort of improvement can the FFAGs bring. Several reviews summarize the projects and studies on fixed field ring methods for ADSR [START_REF] Calabretta | Cyclotrons and FFAG Accelerators as Drivers for ADS[END_REF] [START_REF] Méot | High Power from Fixed-Field Rings in the Accelerator-Driven Sub-Critical Reactor Application[END_REF].

In our approach, we analyze how each of the technologies can meet the requirements of ADSR:

• Energy: Linacs and FFAGs achieve the GeV energy range. For cyclotrons, high energies require separated sectors in order to allow space for the insertion of the RF cavities. However, the orbit separation becomes problematic for CW operation at relativistic energies. Besides, in order to avoid the radial π-mode stop band resonance, the radial betatron frequency ν r must not exceed N s /2 where N s is the number of sectors [START_REF] Livingood | Principles of cyclic Particle Accelerators[END_REF]. Given that ν r ≈ γ, this sets the fundamental energy limit in isochronous cyclotrons. Large number of sectors are needed for relativistic energies. For instance, for the PSI 590 MeV cyclotron, N s = 8. An ultimate 1 GeV proposed PSI cyclotron has N s = 12 (scaling of PSI cyclotron to 10 M W [START_REF] Th | The feasibility of high power cyclotrons[END_REF]). In summary, all three technologies can meet the energy requirement for ADSR (< 1 GeV ).

• Intensity: For ADSR missions requiring a single accelerator with a beam power above 10 M W , at present, only linacs are applicable [START_REF] Abderrahim | Accelerator and target technology for accelerator-driven transmutation and energy production[END_REF]. Besides, due to the RF power consumption, the CW operation can only be achieved with a superconducting technology. In reality though, and due to the reasons mentioned above as well as the challenges in designing a target to withstand several 10 M W of beam power [START_REF] Heidet | Accelerator-Reactor Coupling for Energy Production in Advanced Nuclear Fuel Cycles[END_REF], a multi-beam configuration is preferred. In the 5-10 M W range, cyclotrons are believed to be at their limits. A scaling law by W. Joho [START_REF] Joho | High intensity problems in cyclotrons[END_REF] states that, with constant losses at the extraction electrode, the maximum attainable current scales as:

I max ∝ 1 N 3 t (2.45)
where N t is the number of turns in the cyclotron. Thus, increasing the energy gain per turn is important for high intensity cyclotrons as verified in the PSI machine. FFAG accelerators potentially support repetition rates in the kHz range. Modulation of the RF frequency of 1 kHz or more is possible (unlike the magnetic field, which precludes the use of synchrotrons for ADSR).

A proof of principle machine was built in 1999 at KEK [START_REF] Aiba | Development of a FFAG proton synchrotron[END_REF] which demonstrated rapid acceleration compatible with 1 kHz repetition rate. Another reason for the track of high beam intensity in FFAG is the large Dynamic Aperture which potentially allows large beam size to reduce the space charge effects. In other words, the beam intensity that can be reached, which is a function of the number of particles per bunch (thus of the DA) and the repetition rate, makes FFAG a possible candidate as a proton driver. However, much efforts are still needed to bring these technologies into the M W -level. Recent advances in the non-scaling FFAG can bring this technology into the M W level [START_REF] Johnstone | Isochronous High Intensity Non-scaling FFAG Proton Drivers[END_REF].

• Reliability: This is a common issue of all accelerators. There is no other application that requires the same level of reliability that an ADSR requires. Although the maximum number of beam trips allowed depends on the application as well as the design, this number is obtained from thermo-mechanical considerations on the spallation target as well as the sub-critical assembly. The net useful electricity produced by an ADSR as well as the transmutation rate are limited by the reactor shut-down periods. Fig. 2. [START_REF] Berg | Review of Current FFAG Lattice Studies in North America[END_REF] shows the beam trip rates for existing high power accelerators. Many sources of failure explain the important trip rates as shown in Fig 2 .20: for instance, the high voltage breakdowns due to plasma discharges in the electrostatic devices accounts for 27 %.

The requirement for ADSR is to limit the number of beam trips exceeding few seconds to almost nothing (less than 10 per 3 months operation exceeding 3 seconds). However, as can be observed, with the current state of the art, an ADSR facility will not operate. This justifies much of the efforts put to develop faulttolerant accelerators with redundant system components and even the technology choice. For instance, the MYRRHA project opted for the linac technology in order to obtain a modular situation that is well suited to tackle the availability issue [53]. In case of a multiple accelerator scheme, the previous argument is no longer valid and one needs to evaluate the cost of each technology. • Cost: All the above requirements can somehow be relaxed if one were able to build as many accelerators as needed. In reality though, the cost is the main limitation to deploy the ADSR technology. Comparison of construction cost between a linear and a circular accelerator is not obvious given its dependence on several parameters: magnets, RF cavities, cryostats, beam instrumentation, safety systems, buildings, etc. However, approximated criterion can be derived to compare the linear and circular machines: in [START_REF] Ruggiero | Cost Comparison of Linear and Circular Accelerators[END_REF], it is shown that a circular accelerator is more economical than a linear accelerator provided the number of revolutions for the entire acceleration cycle is larger than a threshold number. From Eq. (2.45), it also appears that in order to aim for high intensity, the number of revolutions cannot be too high. Thus an adequate comparison of the costs is highly dependent on the intensity effects. A general finding is that, for low intensity machines, a circular accelerator is more economical than a linear one (number of turns can be made as large as possible provided some conditions are fulfilled). For very high intensity (several 10 mA), the linear accelerator is more economical. Thus the main idea is to find the threshold value for which the economics of both technologies are the same. Although we will not pursue this goal in the present dissertation, we introduce a quantity that is particularly useful for the analysis: the Levelized Cost Of Beam (LCOB) in analogy with the LCOE:

LCOB = P t Investment t + O&M t + Decommissioning t (1 + r) t P t (BeamPower t × Reliability × time) (2.46)
where the reliability is a weighted function that depends on the form and duration allowed for a beam trip to occur. For instance, ADSR application can tolerate beam trips of very short duration (< 1 s). However, longer beam trips are particularly detrimental because they induce a reactor shut-down.

To conclude, one compares the costs of the two existing most powerful accelerators in the world:

-the SNS linac at Oak Ridge which delivers a proton beam up to 1.4 MW. Its estimated cost as of 2016 is 730 M$ (∼ 600 M$ as of 2006 [START_REF] Mason | SPALLATION NEUTRON SOURCE PROJECT COMPLETION REPORT[END_REF]).

-the HIPA cyclotron at Paul Scherrer Institute which routinely delivers a proton beam of 1.3 MW. Its estimated cost as of 2016 is 120 M$ [START_REF] Seidel | Cyclotrons and superconducting linacs as high intensity driver accelerators[END_REF].

The cyclotrons are thus a natural and cost effective choice for high power applications at the MW-level.

• Beam losses: a very efficient extraction process is mandatory for all technologies; In fixed field rings, the difficulty arises due to the small turn separation. For instance, PSI HIPA uses an electrostatic channel which allows an overall transmission of 99.98 %. The impact of this parameter is important for the accelerator maintenance. The general criterion is that the activation caused by uncontrolled beam losses be limited to the level of 1 W/m to allow hands-on maintenance.

Conclusion

As shown earlier, the main drawback of the ADSR technology is the high LCOE cost compared to other advanced reactor concepts that do not employ an accelerator. Nowadays, this is a show-stopper for any industrial application aiming at producing energy (without dealing with the waste problem). Besides, the reactor is not intrinsically safer than critical reactor concepts, given the complexity of managing the target interface between the accelerator and the reactor core (overpower transients is one example). However, many of the main challenges (high power, safety, reliability) of the installation can be met by using a multi-beam target configuration where the beam power is split between several accelerators. This reduces the power densities on the target, improves the transmutation characteristics which becomes more uniform and ensures a better safety of the installation. In this configuration, cyclotrons have the potential to achieve the requirements, since this technology provides a cheaper cost than any superconducting linac.

The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the exact application of these laws leads to equations much too complicated to be soluble.

-Paul Dirac.

Chapter 3

Beam dynamics studies in Fixed Field Ring Accelerators (FFRA)

In this chapter, we present some aspects of the transverse and longitudinal beam dynamics in FFRA: we start from the basic principles in order to derive the linearized particle equations of motion for FFRA, essentially FFAGs and cyclotrons are considered here. This is a simple extension of a previous work dedicated to linear lattices that we generalized by including the bending terms to ensure its correctness for FFAG lattice. The space charge term (contribution of the internal coulomb forces of the beam) is contained as well, although it will be discussed in detail in another chapter. The emphasis is on the scaling FFAG type in the framework of a collaboration that has been on-going with the goal of better understanding the properties of the 150 MeV scaling FFAG at KURRI in Japan, and progress towards high intensity operation [22]. Some results of the benchmarking work between different codes are presented. Analysis of certain types of field imperfections revealed some interesting features about this machine which required the development of analytical models to better understand the simulation as well as the experimental results. Field imperfections create important change in the number of betatron oscillations of the particle in the accelerator. Therefore, one major question to be answered is whether one can find a way to remediate the significant tune changes in the FFAG. This chapter presents an introduction to the beam stability analysis that we will present in the next chapter.

Beam dynamics in FFAG accelerators

In a circular accelerator, the lattice is essentially dominated by several main magnets, focusing and defocusing to ensure the bending and the focusing of the beam. These magnets are separated by drift spaces used for insertion devices or accelerating cavities to ensure the basic functions of an accelerator. The motion of a charged particle is described as a series of oscillations around a fixed momentum orbit: This is commonly known as betatron motion. There is a great variety of fixed reference orbits. However, all fixed orbits that are considered in this thesis are flat and lie in the median plane of the accelerator. For a circular accelerator, we designate them as closed orbits. There is also a great variety of coordinate systems that can be used in order to describe the motion (global/local, cartesian/cylindrical/spherical,...). In what follows, we opt for the most commonly used, the Frenet-Serret coordinate system moving along the closed orbit: as shown in Fig. 3.1, the vector s is tangent to the reference orbit (shown in red). The transverse motion lies in the plane (x,y): x is in the median plane of the accelerator, is orthogonal to s and points to the center of the curve, while y is such that (x,s,y) form an orthonormal basis. The relative position of the charged particles with respect to the bunch centroid is denoted by z and referred to as the axial coordinate of the bunch. In general, we do not distinguish between s and z if the trajectory is a straight line. However, in the paraxial approximation, this remains valid for circular accelerators.

Equations of transverse motion of the particles including bending and space charge terms

On-momentum transverse equation of motion

This is an extension of previous work [START_REF] Lund | Transverse equilibrium distribution[END_REF]. We first write the equation of motion in the local curvilinear coordinate system (x,s,y) following the reference orbit by applying the relativistic form of Newton's second law and including the centrifugal force to account for the rotation movement:

Transverse restoring force = Centrifugal force + Electromagnetic force

(3.1)
or in equivalent notation:

F r = F c + F em (3.2)
A correct description of the particle motion should invoke relativity. Although Newton's second law is not valid in the relativistic case, it can be re-written in a covariant relativistic form, which ensures that the equation is valid in all inertial frames [START_REF] Leinaas | Classical Mechanics and Electrodynamic[END_REF]. Thus, the Lorentz force is entirely in accordance with special relativity. However, such is not the case in rotating frames which are non-inertial. For a rotating frame, the geometry is non-Euclidean and general relativity should be invoked. However, we will not do the treatment here. We simply introduce a fictitious force, the centrifugal force, according to Newton's definition (although non-covariant), to account for the rotation of the reference frame with respect to the global reference frame. A more elaborate calculation (that avoids the use of the centrifugal force) using the transformation from cartesian to curvilinear coordinate system yields the same result for the linearised particle equations of motion. 1 Now, for particle i in the bunch, we calculate each term of the above equation in the transverse plane (x,y) separately, using the path variable s rather than the time as the independent variable 2 (s is a monotonic function as is time to ensure the correctness of the derivation):

F r = d dt (mγ i v ⊥i ) = mv zi d ds (γ i v zi d ds x ⊥i ) (3.3) = mγ i v 2 zi d 2 ds 2 x ⊥i + mv zi ( d ds x ⊥i ) d ds (γ i v zi ) ≃ mγ b β 2 b c 2 x ′′ ⊥i + mβ b c 2 (γ b β b ) ′ x ′ ⊥i Thus, F r = mγ b β 2 b c 2 x ′′ ⊥i + (γ b β b ) ′ (γ b β b ) x ′ ⊥i (3.4) γ i in Eq. (3.
3) is the relativistic gamma factor defined by:

γ i = 1 q 1 -β 2 i = 1 q 1 -(β 2 i,x + β 2 i,y + β 2 i,z ) (3.5)
However, in the paraxial approximation, the transverse moments of the particle can be neglected compared to its axial moment (β i,x,y ≪ β i,z ) and the bunch mean axial Lorentz factor β b can be defined:

β b = 1 N N X i=1 v i .u z c ; N = number of particles per bunch (3.6)
Thus, in the derivation above, the particle gamma factor γ i was replaced with the bunch gamma factor γ b . Next, we write:

F c = mγ i v 2 i r i u x ≃ mγ b β 2 b c 2 r i u x (3.7)
Also (see appendix 5.6),

F em ≃ qE a ⊥ + qβ b cs × B a ⊥ + qB a s β b cx ′ ⊥i × s -q 1 γ 2 b ∂φ ∂x ⊥ (3.8)
where φ is the beam self potential generated uniquely by the internal forces of the beam. Now, inserting Equations (3.4), (3.7) and (3.8) into (3.2), one obtains:

x ′′ ⊥i + (γ b β b ) ′ (γ b β b ) x ′ ⊥i = 1 r i u x + q mγ b β 2 b c 2 E a ⊥ + q mγ b β b c s × B a ⊥ + qB a z mγ b β b c x ′ ⊥ × s - q mγ 3 b β 2 b c 2 ∂φ ∂x ⊥ (3.9)
which is the general particle equation of motion in the transverse plane.

•Next we write, r i = ρ + x, where ρ is the bending radius of the equilibrium orbit. Using a Taylor expansion to the first order, this gives:

1 r i = 1 ρ + x = 1 ρ (1 + x ρ ) -1 ≃ 1 ρ (1 - x ρ ) = 1 ρ - x ρ 2 (3.10)
Substituting Eq (3.10) into (3.9), and projecting onto the horizontal plane (x), one obtains:

x ′′ + (γ b β b ) ′ (γ b β b ) x ′ + 1 ρ 2 x = 1 ρ + q mγ b β 2 b c 2 E a x - q mγ b β b c B a y + qB a z mγ b β b c y ′ - q mγ 3 b β 2 b c 2 ∂φ ∂x (3.11)
The first order Taylor expansion of the field around the median plane yields:

B a y ≃ B a y0 + ∂B a y ∂x x = B a y0 (1 + ρ B a y0 ∂B a y ∂x x ρ ) = B a y0 (1 -n x ρ ) = mγ b β b c q 1 ρ (1 -n x ρ ) (3.12)
where n = -ρ B a y0 ∂B a y ∂x is the field index.

Next, we substitute Eq (3.12) into Eq (3.11), the x-transverse equation of motion becomes:

x ′′ + (γ b β b ) ′ (γ b β b ) x ′ + 1 -n ρ 2 x = q mγ b β 2 b c 2 E a x + qB a z mγ b β b c y ′ - q mγ 3 b β 2 b c 2 ∂φ ∂x (3.13)
It shows that the longitudinal component of the magnetic field introduces coupling.

•Now, projecting Eq (3.9) onto the vertical plane yields:

y ′′ + (γ b β b ) ′ (γ b β b ) y ′ = q mγ b β 2 b c 2 E a y + q mγ b β b c B a x - qB a z mγ b β b c x ′ - q mγ 3 b β 2 b c 2 ∂φ ∂y (3.14)
Also, in the vacuum aperture of the accelerator, the Maxwell equations must satisfy:

     ∇ . B a = 0 ∇ × B a = 0 (3.15)
which yields

∂B a x ∂y - ∂B a y ∂x = 0 =⇒ B a x = - mγ b β b c q n ρ 2 y (3.16)
so that the y-transverse particle equation of motion becomes:

y ′′ + (γ b β b ) ′ (γ b β b ) y ′ + n ρ 2 y = q mγ b β 2 b c 2 E a y - qB a z mγ b β b c x ′ - q mγ 3 b β 2 b c 2 ∂φ ∂y (3.17)
NB: The beam potential φ generated from the charge distribution will be treated in detail when the problem of the space charge forces is tackled (chapter 4). For the moment one shall bear in mind that, for a single particle theory, the contribution of the last term in the transverse equations of motion is removed.

Off-momentum equation of motion

The preceding analysis applies to particles having the same momenta. Here we investigate the effect of a momentum deviation Δp: We assume that the beam is an ensemble of N particles having different momenta, within a certain interval called momentum spread of the beam Δp/p. This is equivalent to replacing βγmc = p by p + Δp in the preceding equation. Starting from Eq (3.11), one obtains:

x ′′ + (γ b β b ) ′ (γ b β b ) x ′ + 1 -n ρ 2 x = δ ρ + q mγ b β 2 b c 2 1 1 + δ E a x + qB a z mγ b β b c 1 1 + δ y ′ - q mγ 3 b β 2 b c 2 1 (1 + δ) 2 ∂φ ∂x (3.18)
where δ = Δp/p and 1/ρ = qB a y /p. Solving this equation yields a solution of the homogeneous equation plus a particular solution of the inhomogeneous equation. The actual particle coordinate is thus a superposition of the betatron oscillation plus an additional contribution due to momentum errors. The off-momentum solution can be found by using the parameter variation method.

On solving the equations of transverse motion

In the accelerator physics literature, it is generally assumed that, when the damping term is removed, the above equations are of the form of a Hill's equation, for which general solutions exist and can be written in a closed form, using Floquet theory. The main assumption that allows to obtain such a result is that the coefficients of the above equation are periodic in time (or s). This gives rise to an amplitude function A(s) = p β(s) that varies periodically as well [START_REF] Courant | Theory of the alternating-gradient synchrotron[END_REF] (β in this context refers to the betatron functions). In reality, this may be a good approximation in certain cases, for instance a synchrotron where the accelerated orbit coincides with the closed orbit, but it may not be a good approximation for other types of circular machines such as FFAGs where the accelerated orbit moves outwards and therefore the coefficients of the above equation are no longer periodic in s. In what follows, we use the method of Wentzel-Kramers-Brillouin (WKB) in order to find approximate solutions of the equation of motion.

The WKB approximation

The previous equations of motion (Eq.3.13 & Eq.3.17), when there is neither any applied transverse electric field nor any longitudinal magnetic field components, and neglecting the space charge term become:

         x ′′ + (γ b β b ) ′ (γ b β b ) x ′ + 1 -n ρ 2 x = 0 y ′′ + (γ b β b ) ′ (γ b β b ) y ′ + n ρ 2 y = 0 (3.19)
The above equations can be written in the standard form:

x ′′ + p(s)x ′ + q(s)x = 0 (3.20)
This is a 2nd order linear differential equation with variable coefficients, for which, in general, a closed form solution is not known. We will first seek to remove the first order derivative in the above equation in order to simplify the analysis. In what follows, we use a method that has been developed in [START_REF] Hetrick | Dynamics of Nuclear Reactors[END_REF] (in the context of nuclear reactors) in order to express approximate analytical solutions of the point-reactor differential equations when the reactivity is an arbitrary specified function of time: Without any loss of generality, we assume that a solution of Eq (3.20) can be written in the form: x(s) = u(s)v(s). Thus Eq (3.20) becomes:

uv ′′ + (2u ′ + pu)v ′ + (u ′′ + pu ′ + qu)v = 0 (3.21)
and one only has to choose the function u(s) such that 2u ′ + pu = 0, i.e,

u(s) = exp - 1 2 Z s s 0 p(h)dh (3.22)
Now, the first order derivative is removed by making the Liouville-Green transformation

x(s) = v(s). exp - 1 2 R s s 0 p(h)dh
, which yields:

v ′′ + q(s) - 1 2 dp ds - p 2 4 v = 0 (3.23)
Now, back to Eq (3.19), the above formalism applies by simply replacing

p(s) = (γ b β b ) ′ (γ b β b )
which yields:

x(s) = v(s). exp - 1 2 Z s s 0 p(h)dh = v(s). exp - 1 2 Z s s 0 (γ b β b ) ′ (γ b β b ) dh = C √ γ b β b v(s) (3.24)
and v(s) is a solution of Eq (3.23) above. The above treatment shows why, in general, we refer to the first derivative term in Eq (3.19) as the "damping term". Next, we rewrite Eq (3.23) as:

v ′′ -h 2 (s)v = 0 ; h 2 (s) ∈ R ± (3.25)
where

h 2 (s) = - 1 -n ρ 2 + 1 2 d ds (γ b β b ) ′ (γ b β b ) + 1 4 (γ b β b ) ′ (γ b β b ) 2 ; horizontal plane h 2 (s) = - n ρ 2 + 1 2 d ds (γ b β b ) ′ (γ b β b ) + 1 4 (γ b β b ) ′ (γ b β b ) 2 ; vertical plane (3.26)
We now seek an approximate solution of Eq (3.25). We postulate: v = A 0 exp [ϕ(s)] where ϕ is a phase term (ϕ(s) ∈ C). Substitution into Eq (3.25) yields a second order non-linear differential equation:

d 2 ϕ ds 2 + dϕ ds 2 -h(s) 2 = 0 (3.27)
If dϕ ds is slowly changing (adiabatic change), an iterative solution can be found. The first iteration is:

dϕ ds ≃ h(s) (3.28)
The second iteration is:

dϕ ds 2 = h 2 - d 2 ϕ ds 2 ≃ h 2 - dh ds (3.29)
Now, we treat dh ds as small and we make the distinction between two cases (focusing/defocusing), according to the sign of h 2 (s):

If h 2 (s) > 0: dϕ ds ≃ h 2 - dh ds 1/2 ≃ h - 1 2h dh ds (3.30) =⇒ ϕ ≃ Z hds - 1 2 ln |h| (3.31) v = A 0 exp[ϕ(s)] ≃ A 0 p |h| exp Z hds (3.32)
and a complete approximate solution of Eq. (3.19) is given by:

x(s) = 1 p |h| 1 √ β b γ b A exp �R hds + B exp � - R hds (3.33)
where A and B are constants and h is defined in Eq (3.26) above.

If h 2 (s) < 0:

Repeating the same steps with h imaginary leads to:

x(s) = 1 p |h| 1 √ β b γ b C exp � i R |h|ds + D exp � -i R |h|ds (3.34)
•At the classical turning point, i.e when h = 0, the amplitude function diverges and the above solution is no longer valid. Such a condition is fulfilled at the point where the magnetic field changes its sign, i.e when the gradient is alternated. However, this difficulty can be handled by finding the higher order terms of the WKB approximation.

•The main approximation that we made in this analysis was to assume that the phase ϕ(s) was slowly changing which implies that h 2 (s) defined in Eq. (3.26) above is slowly changing as well. This adiabatic change becomes less and less valid when the energy gain per turn becomes important, since then, the change in γ b β b becomes important. This is generally a requirement to limit the space charge effects and also for fast crossing of the resonances.

•The amplitude of the oscillations as well as its frequencies are strongly dependent on the form of the magnetic field (to the first order in this analysis, the bending radius and the field see Eq (3.26)). As a consequence, any defect in the design of the magnets may be a source of imperfections and should be corrected.

In conclusion, this shows that, for slowly varying parameters of the linearized particle equations of motion, the damping of the transverse particle oscillations is given by:

x ∝ 1 p |h| 1 √ β b γ b (3.35) x ′ ∝ p |h| 1 √ β b γ b (3.36)
so that

ǫ rms = h x 2 x ′2 - xx ′ 2 i 1/2 ∝ 1 β b γ b (3.37)
This result will be discussed in more detail when investigating the damping law in FFAGs.

Finally, let's remind the Lagrange formula for the integration of the second order inhomogeneous linear differential equation (Eq. 3.18):

x ′′ + p(s)x ′ + q(s)x = f (s) (3.38)
Knowing a fundamental system of solutions x 1 , x 2 of the homogeneous Eq. (3.20), one can obtain a solution of (3.38) by means of the formula [START_REF] Dieudonné | Calcul infinitésimal[END_REF]:

x p (s) = Z s s 0 x 1 (t)x 2 (s) -x 2 (t)x 1 (s) W (t) f (t)dt (3.39)
where W(s) is the Wronskian of the two fundamental solutions defined by,

W (t) = x 1 (t)x ′ 2 (t) -x 2 (t)x ′ 1 (t) (3.40)

Longitudinal beam dynamics

The longitudinal fields usually exhibit a non-linear behaviour. They play two major roles: accelerate the beam in the longitudinal direction and possibly providing phase focusing (synchrotron acceleration).

We first start by deriving the longitudinal equation of motion: projecting Eq. (3.2) in the longitudinal plane yields [START_REF] Lund | Transverse equilibrium distribution[END_REF],

dp dt = d dt (mγv z ) = qE a z -q � v x B a y -v y B a x -q ∂φ ∂z i (3.41)
Now injecting Eq. (3.12) and (3.16) into (3.41), one obtains:

dp dt = qE a z -mγ b β b c v x ρ + (v y y -v x x) n ρ 2 -q ∂φ ∂z i (3.42)
which shows that the acceleration of the beam obeys a complicated equation that couples both planes and includes the field index n to the first order (without even considering the space charge term). One can conjecture that particles with different transverse momenta v x and v y in the bunch can have a seriously different behavior throughout the acceleration (the bunch length could increase for instance).

A remedy to this problem that was first discovered by Veksler in 1944 was to introduce the principle of phase stability: in order to ensure this principle, the requirements are [START_REF] Livingood | Principles of cyclic Particle Accelerators[END_REF]:

• The peak voltage across the accelerating gap must be somewhat greater than would otherwise be needed;

• The particle orbits must be such that a change in energy and momentum must be accompanied by a change in the period of revolution;

• The variation of the oscillator's RF frequency shall be made sufficiently slow (this is generally the case when the particles approach relativistic velocities);

Transition energy

The momentum compaction factor α is defined as:

α = dC/C dp/p = dR/R dp/p (3.43)
where C is the orbit circumference and R is the equivalent radius of the orbit defined as R = C/2π. Thus, α represents the fractional change of the orbit length as a function of the fractional change in momentum. Furthermore, it can be shown that [START_REF] Livingood | Principles of cyclic Particle Accelerators[END_REF],

dτ τ = - df f = 1 α - 1 γ 2 dp p = 1 γ 2 tr - 1 γ 2 dp p = -η dp p (3.44)
where γ tr is the transition energy of the particle, τ the revolution period, f the revolution frequency and η the slip factor. Thus, if γ < γ tr , an increase in momentum results in an increase of the revolution frequency and vice-versa. The particles that arrive earlier than the synchronous phase will experience a lower accelerating field and the particles that arrive later will experience a larger accelerating field. This maintains the particles oscillating around the synchronous phase and can only be achieved if the beam is provided to the leading edge of the wave as shown in Fig. 3.2, if the oscillator frequency is slow enough as explained earlier and if η > 0. However, the RF transverse electric fields also act on the beam near the edges of the accelerating gap (see Fig. 3.3), so that the particles that are focused longitudinally will experience transverse defocusing forces [START_REF] Wangler | RF linear accelerators[END_REF]. This effect can be understood and estimated using the transverse particle equations of motion (Eqs. (3.13) and (3.17) derived earlier) if we have an estimate of the transverse component of the electric fields E a x and E a y . For instance, if we consider the average effect of E a

x and E a y , one can use Eq (3.39) to estimate the closed orbit distortion due to the RF cavity.

Longitudinal equation of motion

It appears from the previous analysis that the relative phase advance of the particle with respect to the synchronous one is the relevant quantity to use in order to derive the equation of motion in the longitudinal plane. In order to do so, we derive two first order equations, the first for the energy gain after the cavity crossing, the second for the phase change relative to the synchronous particle. The derivation steps are mostly inspired from [START_REF] Reiser | Theory and Design of Charged Particle Beams[END_REF]:

For the first equation, we start from Eq.(3.42), and knowing that dp/dt = dE/ds (E being the kinetic energy of the particle) and assuming that the longitudinal electric field can be written in the simple form

E a zi = E m cosϕ, one obtains, dE ds = qE m cos ϕ -q ∂φ ∂z i (3.45) dE s ds = qE m cos ϕ s (3.46)
where we assume that the space charge force is zero at the center of the bunch (where the synchronous particle is located) and where we neglect the contribution of the transverse velocities to the change in the longitudinal momentum. E m represents the peak electric field while ϕ is the phase of the particle with respect to the peak and is defined by:

ϕ = ω rf t - ω rf s v 0 (3.47)
This yields,

d ds (ΔE) = dE ds - dE s ds = qE m (cos ϕ -cos ϕ s ) -q ∂φ ∂z i (3.48)
For the second equation of the phase difference, one can write:

d ds (Δϕ) = d ds (ϕ -ϕ s ) = - ω rf c 1 β - 1 β s (3.49)
Using the relation between β and γ and making a Taylor expansion, Eq (3.49) becomes:

d ds (Δϕ) = - ω rf c Δγ β 3 s γ 3 s = - 2π λ ΔE mc 2 1 β 3 s γ 3 s (3.50)
where λ = 2πc/ω rf is the synchrotron wavelength of the wave. Extracting ΔE from Eq (3.50) and substituting in (3.48) yields a second order non-linear differential equation for longitudinal motion:

d ds β 3 s γ 3 s d ds (Δϕ) = - 2π λ qE m mc 2 (cos ϕ -cos ϕ s ) - 2π λ q mc 2 ∂φ ∂z i (3.51)
Now, we assume that the energy change occurs slowly enough to consider that the factor β 3 s γ 3 s is almost constant during one phase oscillation period. We also neglect the space charge term here ∂φ ∂z ∼ 0. It results,

d 2 ϕ ds 2 + 2π λ qE m β 3 s γ 3 s mc 2 (cos ϕ -cos ϕ s ) = 0 (3.52)
For small amplitude oscillations around the synchronous phase, one has,

cos ϕ -cos ϕ s = cos(ϕ s + Δϕ) -cos ϕ s ≈ -Δϕ sin ϕ s (3.53)
so that Eq (3.52) reduces to a harmonic oscillation

d 2 ϕ ds 2 - 2π λ qE m sin ϕ s β 3 s γ 3 s mc 2 Δϕ = 0 (3.54)
Thus, stable oscillations will occur only if ϕ s < 0.

Now, we will transform the above equations into a more simplified form in order to find the points of the separatrix, the limit of stable trajectories in the longitudinal plane (and to be exploited later on for the code benchmarking): the derivation steps are in part inspired from [START_REF] Wangler | RF linear accelerators[END_REF]: To begin with, we introduce the following notation:

ǫ = ΔE mc 2 ; A = 2π β 3 s γ 3 s λ ; B = qE m mc 2 (3.55)
Equations (3.48) and (3.50) yield:

     ǫ ′ = dǫ ds = B (cos φ -cosφ s ) φ ′ = dφ ds = -Aǫ =⇒ φ ′′ = -AB (cos φ -cosφ s )
which, after integration yields:

Aǫ 2 2 + B (sin φ -φ cos φ s ) = H ǫ,φ (3.56)
where H ǫ,φ is a constant of integration identified as the Hamiltonian. Note that the previous result is based on the assumption that the acceleration rate is small. Furthermore, it can be shown that points of the separatrix must satisfy:

H ǫ,φ = -B (sin φ s -φ s cos φ s ) (3.57)

Beam dynamics in scaling FFAG

All of the above remains valid for scaling FFAG. The only difference is, now, we assume that the magnetic field obeys a well defined law, that we generally refer to as the scaling law. The scaling law guarantees that the number of betatron oscillations remain constant [START_REF] Symon | Fixed-Field Alternating-Gradient Particle Accelerators[END_REF]. However, it is susceptible to change due to imperfections. Therefore, in our analysis, we decide to take the magnetic field imperfections into account. In this section, we first show the origin of the scaling law for this type of accelerators. Then we summarize some of the benchmarking work undertaken to validate some of the simulation results. We also discuss about the main assumptions for the calculation of the scaling factor from the field map. The emphasis is made on the mathematical correctness of the calculation. We conclude this section by investigating the damping law in FFAGs in general, and scaling FFAGs in particular. Comparison of the simulation results that account for the non-linearities of the field with the analytical results based on the linearization of the equations of motion showed good agreement and helps validate the modelling work of the FFAG lattice.

We first re-write Eq (3.19) by dropping the damping term (fixed momentum) and making the transformation s = θR where θ is the azimuthal angle with the periodicity of the machine i.e. 2π/N , N being the number of sectors, and R the equivalent radius for each closed orbit, i.e. R = S/2π and S is the integrated path length over one closed orbit. This is the same transformation that can be found in the original paper by Symon in 1956 [START_REF] Symon | Fixed-Field Alternating-Gradient Particle Accelerators[END_REF]. Thus, we obtain:

         d 2 x dθ 2 + R 2 ρ 2 (R, θ) [1 -n] x = 0 d 2 y dθ 2 + R 2 ρ 2 (R, θ) n y = 0 (3.58)
Therefore, two sufficient conditions to have both the vertical and horizontal betatron oscillations constant with respect to the momentum, thus the forcing terms in the above equations constant, are:

∂n ∂p θ=const = 0 (3.59) ∂ ∂p R ρ θ=const = 0 (3.60)
These conditions are called the cardinal conditions of a scaling FFAG [START_REF] Machida | Fixed Field Alternating Gradient[END_REF]: the first one expresses the constancy of the field index with respect to the momentum while the second one expresses the similarity between the different closed orbits. In other words, closed orbits corresponding to different energies are homothetic with respect to the centre of the machine, therefore explaining the scaling property of this type of accelerator. However, one fundamental hidden assumption in this analysis is that the number of betatron oscillations of the accelerated orbit is the same as that obtained from the closed orbits. As was shown in Eq. (3.25), the forcing term of the equation of motion of the accelerated orbit differs from that of the closed orbits due to the presence of the terms containing the change in momentum of the particles, namely (p ′ /p) ′ and (p ′ /p) 2 . However, for small acceleration rates, this effect can be neglected.

If we define the chromaticity as the variation of the betatron tune with the momentum deviation, then a scaling FFAG is essentially a zero-chromatic FFAG. More detailed analysis of these properties will be tackled when developing the hard edge model. It results from both cardinal conditions that the field index averaged over the azimuth of the closed orbits is constant as well (independent of the radius). This is also defined as the scaling factor k:

k = R B ave dB ave dR ≈ R B ave dB ave dx = - R ρ n = const (3.61)
Integrating Eq. (3.61) yields:

B ave (R) = B 0 × ( R R 0 ) k (3.62)
where B 0 is the vertical component of the magnetic field3 averaged in azimuth for a given radius R 0 . This formula describes the radial dependence of the azimuthally averaged field over the entire width of the magnet.

In conclusion, for a radial sector type FFAG, the vertical component of the magnetic field in the median plane can be written in the general form:

B(R, θ) = B 0 × ( R R 0 ) k × F (θ) (3.63)
where F is the flutter function and R is the equivalent radius of the closed orbit as specified earlier. However,as will be shown below, this formula remains valid if we define R as the orbiting radius seen from the centre of the circular accelerator. The main assumption of this definition is that the function B(R,θ) is separable in radial and azimuthal coordinates, which will be discussed later.

If k is fixed, then R 0 can be chosen arbitrarily: to show this, let's assume that k is fixed and does not have any radial dependence. Then for an arbitrary radius R 1 (Eq 3.62) gives:

B 1 (R 1 ) = B 0 × ( R 1 R 0 ) k (3.64)
and injecting (Eq 3.64) into (Eq 3.62), one obtains:

B ave (R) = B 1 × ( R R 1 ) k (3.65)
which proves the assertion above. From this, one can obtain the value of the average scaling factor of the KURRI 150 MeV FFAG, by computing the average z-component of the magnetic field in the median plane for different radii and by means of (Eq 3.66): in principle, if the magnets obey to the exact scaling law, k has to remain constant. However, due to imperfections, one can observe some variations in k.

k(R) = log( B ave (R) B 0 ) log( R R 0 ) (3.66)

KURRI 150 MeV scaling FFAG and benchmarking work

Following the FFAG14 workshop held at BNL [73], a simulation campaign was established to benchmark several simulation codes. The main objective was to provide reliable modeling tools for FFAG type of accelerators and to better explain the results of the experiments at the KURRI 150 MeV scaling FFAG [START_REF] Sheehy | Progress on simulation of fixed field alternating gradient accelerators[END_REF]. Due to the nonlinearities of such a machine, the beam dynamics are intricate to simulate and require careful modeling. The longitudinal and transverse beam dynamics in FFAGs differ significantly from the one in fixed orbit machines such as synchrotrons. A more detailed discussion about the correct modeling in FFAG can be found in [START_REF] Berg | Correct Tracking in FFAGs[END_REF].

In what follows, we show some results of the benchmarking work: all the results that we present show good agreement either with the theory or with other codes. However, this required careful modelling and iterations before convergence. This is a crucial step to validate the models that we have built and will discuss in the following chapter.

Design of the injection line into the KURRI FFAG:

The first benchmarking test was carried out for the injection line into the main ring. The injection line transports the beam from the linac to a matching point inside the ring where a stripping foil has been inserted for the injection process. The orbit suitable for acceleration is the first closed orbit shown in Fig. 3.4. The injected beam has to match the properties of this reference orbit at a given azimuth that in this case is the stripping foil, located where the injection orbit and the closed orbit intersect. An easy and straightforward solution to find the injection line is to track the beam backwards from the matching point. This approach and the necessary computing tools were developed previously to design the injection line into the 800MeV/amu high power cyclotron [START_REF] Tahar | DESIGN OF THE INJECTION INTO THE 800 MeV/amu HIGH POWER CY-CLOTRON[END_REF]. [START_REF] Méot | Zgoubi Users Guide[END_REF] show good agreement (green and thick black curves superposed). The EARLIETIMES results were provided by Y. Ishi. The color density plot represents the magnetic field associated to the magnets composing the lattice.

Transverse beam dynamics comparison

The main property of this type of accelerator is that, in the ideal case, the number of betatron oscillations per turn, also referred to as the tune, is energy-independent (this is the so-called zero-chromaticity condition). This results from the r k dependence of the field and avoids any resonance crossing during acceleration. Therefore, one important parameter for the benchmarking is the evolution of the tune over the acceleration cycle. Comparison of the tracking results from several codes using OPERA 3D fieldmaps of the dipole triplet cell show a good agreement (see Fig. 3.5). However, this demonstrated that the tune variations are non-negligible (of the order of 15% in the vertical plane and 4% in the horizontal plane). Note that the tune measurements are not in good agreement with this result (see Fig. 3.6) which is mainly due to the difference between the simulated and the real field of the magnets. This is investigated in detail in next chapter where a stability analysis of FFAG lattices to field imperfections is elaborated. The details of the measurement as well as the methods and tools that were used to characterize the 150 MeV KURRI FFAG are discussed in [START_REF] Sheehy | Characterization techniques for fixed field alternating gradient accelerators and beam studies using the 150 MeV proton FFAG[END_REF].

Longitudinal beam dynamics comparison

In a scaling FFAG, the synchronous particle orbit spirals outwards with increasing energy. Therefore, the RF phase of the cavity must evolve in such a way as to follow the change in the revolution frequency of the particle.

The revolution period of the particle is given by:

τ rev = 2πR βc = τ 0 p p 0 -k/(1+k) E E 0 (3.67)
As shown in Fig. 3.7, the revolution frequency of the particle was calculated using the tracking in the fieldmap and the RF phase law implemented accordingly. Differentiating the particle momentum with respect to the radius, one can show that γ tr = √ 1 + k. Thus, for the KURRI 150 MeV FFAG, γ tr = 2.93. Given that the maximum kinetic energy of the protons yields γ = 1.16, this guarantees that η > 0 so that there is no crossing of the transition energy throughout the acceleration. i) No acceleration As shown earlier, the longitudinal particle trajectories for small acceleration rates are described by Eq.(3.56). As a first benchmarking test of the phase space trajectories in the longitudinal plane, we assumed no energy gain in the RF cavity, therefore φ s = -π/2. A contour plot of H ǫ,φ for the KURRI 150 MeV FFAG at injection is shown in Fig 3.8 (a). The separatrix, in dashed lines, shows the maximum excursion in phase space that any particle can have and still be captured in the bunch. Particles outside the separatrix are lost. Using the tracking code ZGOUBI, one calculated the stable phase space trajectories (delimited by the separatrix) by varying the initial particle coordinates in phase space. The results are shown in Fig 3.8 (b) and are in good agreement with the analytical results. Also, it is important to note that at the neighbourhood of the synchronous phase, the motion reduces to a harmonic oscillation (see Eq. 3.54) and the trajectories are ellipses. The conditions of the simulation are as follows: E inj = 11M eV , V rf = 4kV , τ rev = 0.625µs.

ii) With acceleration Comparison of the longitudinal phase space shows good agreement between ZGOUBI and SCODE (see Fig. 3.9). In the case of small synchrotron oscillations (the initial phase is close to the synchronous one), one can observe some differences essentially related to the initial conditions which are not the same in both codes. However, this is sufficient and useful when comparing multi-particle tracking results with slightly different settings.

Multi-particle tracking

Tracking particle distributions throughout the acceleration cycle and using the 3D field map of the dipole triplet cell revealed important emittance jumps in the horizontal plane at some specific locations (see Fig. 3.10). No such behavior is obtained in the vertical plane. These results were confirmed by both SCODE and ZGOUBI (although the initial conditions are different, rather by choice (the normalized emittance at injection is ∼ 7 times higher in ZGOUBI. Also, the particle distributions are not the same (Waterbag vs KV)), the results exhibit the same behavior). In an effort to explain the origin of that effect, several tests were carried out. Particle tracking using 2D mid-plane field map finally revealed that the effect disappears. This is still an open problem: one cannot assert that the quality of the 2D field map is better than that of the 3D field map. However, it is found that the observed effect is sensitive to the order of the interpolation polynomial used for the calculation of the field and its derivatives as well as the integration step size and the field map mesh size. All the above parameters were finely tuned to give convergent results as part of the benchmarking effort. As a consequence, it has been deemed necessary to compute new field maps of the FFAG dipole cell with finer mesh size in order to overcome and further investigate the origin of the problem. Further investigation will be needed to check the codes calculation in particular. The stable fixed point lies at φ = φ s ≡ -π/2[2π] (green point). (b) The initial RF phase is 30 degree and the synchrotron oscillations are much smaller around the synchronous phase ϕs = -30 deg. The phase oscillations are susceptible to the initial conditions which are not exactly the same in the two codes. However, the differences are less than 0.3 % so that the agreement is quite good.

Scaling factor analysis

The main objective of the following study is to investigate the origin of the tune variations in the KURRI FFAG.

Using the TOSCA field map ("TOSCA map f810 d1020.dat") of the KURRI machine, we calculated the average magnetic field ( R Bdθ/ R dθ) for one sector. A plot of the vertical component of the magnetic field for different radii is shown in (Fig. 3.11) below. The average scaling factor is then deduced by applying Eq (3.66) where R 0 is the reference radius that in our case was chosen to be R 0 = 445cm (see Fig. 3.12). These results show about 18% variation on the average scaling factor from injection to extraction (similar errors for the tunes were observed). In order to understand the evolution of the scaling factor with the radius, one can apply the previous formula to evaluate the local scaling factor as a function of the radius for different azimuthal positions. The results are shown in (Fig. 3.13): for the different azimuthal angles below 4 degrees (namely: 0, 1.5 and 3 degrees), the local scaling factor (which is obtained by using the formula:

k(R, θ = θ 0 ) = log( B(R, θ = θ 0 ) B 0 ) log( R R 0 )
) is constant and varies from 7.4 to 7.6. This corresponds to the central region of the main focusing magnet (θ = 0) where the magnetic field is constant (See Fig. 3.11). However, for the azimuthal positions above, the variations of the scaling factor are independent, in particular in the fringe field region. For the D-magnet, the average value of k is significantly higher than 7.6: this will be discussed in more detail later on in this chapter.

So, the assumption that the average scaling factor is constant does not apply in the FFAG ring. This could explain the strong tune variations that were obtained earlier (and will be discussed later on). However, is the assumption that R 0 can be chosen arbitrarily in Eq. (3.62) still valid? In order to answer this, let's rewrite Eq. (3.62) in the form

B ave (R) = B 0 × ( R R 0 ) k 0 (R) (3.68)
where now, we allow the radial variations of the scaling factor, k 0 (R).

If we choose our reference radius to be R 1 instead of R 0 , then Eq. (3.68) rewrites:

B ave (R) = B 1 × ( R R 1 ) k 1 (R) (3.69)
If the assumption were to remain true, then k 0 (R) = k 1 (R). However, we will derive an explicit transformation from (R 0 , B 0 , k 0 (R)) to (R 1 , B 1 , k 1 (R)):

On the one hand, one has,

B ave (R) = B 0 × ( R R 0 ) k 0 (R) = B 1 × ( R R 1 ) k 1 (R) (3.70)
On the other hand,

B 1 = B 0 × ( R 1 R 0 ) k 0 (R 1 ) (3.71)
By injecting (Eq 3.71) into the first term of (Eq 3.70), one obtain:

B ave (R) = B 1 × ( R 0 R 1 ) k 0 (R1) × ( R R 0 ) k 0 (R) (3.72)
Now, taking the Logarithmic of these expressions, one obtains:

log(B ave (R)) = log(B 1 ) + k 1 (R) log( R R 1 ) = log(B 1 ) + k 0 (R 1 ) log( R 0 R 1 ) + k 0 (R) log( R R 0 ) (3.73) so that, k 1 (R) = log( R R 0 ) log( R R 1 ) × k 0 (R) + log( R 0 R 1 ) log( R R 1 ) × k 0 (R 1 ) (3.74)
It results from this equation that:

• k 0 (R) = k 0 (R 1 ) = const =⇒ k 1 (R) = k 0 (R)
, which was proved in the first paragraph.

• R 0 = R 1 =⇒ k 0 (R 0 ) = k 1 (R 1 ) which is trivial. • R = R 0 =⇒ k 1 (R 0 ) = k 0 (R 1 )
which can be easily verified from (Fig. 3.

14).

And more generally:

∀(i, j) ∈ N × N, k i (R j ) = k j (R i )
where R i and R j are reference radii for k i and k j respectively.

In order to check this result, we compute the average scaling factor by varying the reference radius R 0 in the previous Equation (3.68). The results are shown in Fig. 3.14. Average scaling factor for different reference radii: Eq. (3.74) above was applied to transform (R0 = 445cm, B0, k0(R)) into (R1 = 465cm, B1, k1(R)) which gave perfect agreement with the result obtained from the field map. The discontinuity point simply corresponds to the case where R = R1 for which the formula is not defined. (all the other possible transformations in this plot were also checked and gave perfect agreement as expected).

In conclusion, k(R), as defined earlier, is a relative quantity and cannot be defined in an absolute way, unless it is constant. So a reference radius in the field map has to be well defined with respect to which the scaling factor can be determined. Yet how to calculate k(R) in a reliable way? Also, is there a reference radius that can be chosen in order to obtain consistent results? The next section answers this question.

A more accurate way to compute the average scaling factor

Let's use the original (correct) definition of the scaling factor:

k = dB/B dR/R (3. 75
)
where B is the field averaged in azimuth. This form has the merit that it does not have any assumption as to the form of the scaling factor, whether it is constant or is a function of R. Thus it is more rigorous to apply in order to compute k. From the field map, one applied Eq. (3.76) below:

k(R i ) = B i+1 -B i B i R i+1 -R i R i (3.76)
where R i+1 -R i is defined by the mesh size (here R i+1 -R i = 1 cm). The result is shown in Fig. (3.15). The correct scaling factor (curve in black) was obtained from Eq. (3.76) where no assumption is made as to its form: it can be seen that each one of the other curves which were obtained from the definition specified in Eq. (3.66) is valid only in the vicinity of the reference radius R 0 where the variations of k are still not too important. Note the oscillations of k which are due to the granularity of the field map (radial step size of 1 cm).

Conclusion:

The radial scaling law ) is only valid in the vicinity of the reference radius R 0 , simply because it is not a correct general solution of Eq. (3.75). Thus, when defining this equation, one has to choose R 0 carefully. However, as will be shown later, this form is particularly useful to simplify the analysis when writing the equations of motion in the vicinity of each of the closed orbits.

B(R) = B 0 × ( R R 0 ) k(R

The average scaling factors for the F and D magnets

The definition of k ( Eq.(3.75)) was applied again in order to compute the exact values of k for the focusing and defocusing magnets individually (the average value of the field B is taken over the width of each of these magnets). The results are shown in Fig. (3.16). For the F-magnet, the variations of k are weak, however for the D-magnet, one observes that the scaling factor changes significantly around an average value k ≈ 9.1. In order to understand the origin of these strong variations, mainly in the defocusing magnet, we compute the map of the average field index of the two magnets. As can be seen in Fig. 3.17, the main source of discrepancy is observed in the interaction region between the magnets and seems to mostly affect the defocusing one. This effect is mostly dominant in the injection region. In the following section, we develop an analytical model that takes into account the scaling imperfections of the magnetic field as a function of the radius. The next chapter will focus on the analysis of the effects of the azimuthal imperfections of the magnetic field on the beam dynamics.

Analytical solution of the magnetic field to account for scaling imperfections

The average field index k, is defined by Eq. (3.75). Based on the previous calculation of the scaling factor, it turns out that k is a function of R so that the previous form of the magnetic field (Eq. (3.62)) is not quite accurate to describe the magnets. In order to obtain the new form of the magnetic field, while keeping the same definition of the scaling factor, let's assume that k can be fitted with an n-order polynomial which is a function of the radius: k has the form

k(R) = n X i=0 a i R R 0 i (3.77)
Then by equating Eq. (3.75) and Eq. (3.77), one obtains:

dB B = n X i=0 a i R i 0 R i-1 dR = a 0 dR R + n X i=1 a i R i 0 R i-1 dR
which gives after integration: The new scaling law takes into account the R-dependence of the scaling factor and thus is more accurate to describe the magnet as can be seen from Fig. 3.18. The final parameters of the fitting are: n=3, a 0 = 7.6, a 1 = 4.74 × 10 -3 , a 2 = -5.39 × 10 -5 , a 3 = 9.04 × 10 -8 and R 0 = 480 cm.

ln( B B 0 ) = a 0 ln( R R 0 ) + n X i=1 a i R i 0 Z R R 0 R i-1 dR (3.78) B(R) = B 0 exp " a 0 ln( R R 0 ) + n X i=1 a i R i 0 (R i -R i 0 ) i # so that the general form of the magnetic field becomes, B(R, θ) = B(R)F (θ), i.e, B(R, θ) = B 0 R R 0 a 0 × exp n X i=1 a i R i -R i 0 i × R i 0 ! × F (θ) (3

Flutter function

If B(R,θ) were separable in radial and azimuthal coordinates, then one should obtain:

F (θ) = B(R, θ) hBi (R) (3.80)
Yet, when computing F(θ) for various radii (see Fig. 3.19), one observed that F has also a radial dependence. This proves that B(R,θ) is not a separable function.

In conclusion, the general form of the magnetic field is:

B(R, θ) = hBi (R) × F (R, θ) (3.81)
F (R) = hF (R, θ)i = 1 (3.82)
In order to better understand the radial dependence of the flutter, we start from Eq. (3.79) where we evaluate the field along a closed orbit (instead of a fixed radius orbit). We introduce the transformation: R(θ) = hRi × (1 + r(θ)) where r(θ) is a periodic function describing the periodic scalloping around the fixed radius orbit hRi (of the order of 10 %). This yields:

B(R(θ), θ) = B 0 hRi R 0 a 0 × [1 + r(θ)] a 0 × exp n X i=1 a i R(θ) i -R i 0 i × R i 0 ! × F (θ) (3.83)
The above equation is only valid under the assumption that all orbits scale with respect to the center of the ring. Otherwise r(θ) needs to be re-calculated for each orbit. •If the scaling factor is constant (i.e. R-independent), then a i = 0 ∀i ≥ 1 and Eq. (3.83) can be written in the form:

B(R(θ), θ) = B(hRi) × F 1 (θ) F 1 (θ) = [1 + r(θ)] a 0 × F (θ) (3.84)
This shows that the flutter function of the equilibrium orbit is constant and differs from that of the fixed radius orbit through a simple formula. This is only true if there are no scaling imperfections.

Turn separation in a scaling FFAG

In this section, we will derive the formula of the turn separation as a function of the energy, which will be useful for later analysis. We start from the scaling law:

r 2 = r 1 p 2 p 1 1 1 + k (3.85)
which is equivalent to:

ln r 2 r 1 = 1 1 + k ln p 2 p 1 (3.86) Injecting, pc = βγm 0 c 2 = 1 - 1 γ 2 1/2 γm 0 c 2 = � γ 2 -1 1/2 m 0 c 2 (3.87)
into 3.86, one obtains:

ln r 2 r 1 = 1 2(1 + k) ln γ 2 2 -1 γ 2 1 -1 = 1 2(1 + k) ln E 2 2 -m 2 0 c 4 E 2 1 -m 2 0 c 4 (3.88)
where E i refers to the total energy of the particle. We write: E 2 = E 1 + ΔE such that ΔE << E 1 (the energy gain per turn is in general negligible compared to the relativistic mass). Thus,

E 2 2 -m 2 0 c 4 E 2 1 -m 2 0 c 4 = (E 1 + ΔE) 2 -m 2 0 c 4 E 2 1 -m 2 0 c 4 ≈ 1 + 2E 1 E 2 1 -m 2 0 c 4 ΔE (3.89)
Injecting this into Eq. 3.88, and applying the exponential, one obtains:

Δr = r 2 -r 1 = r 1    1 + 2E 1 E 2 1 -m 2 0 c 4 ΔE 1 2(1 + k) -1    (3.90)
or, in a more general form:

Δr(E) = r(E)    1 + 2E E 2 -m 2 0 c 4 ΔE 1 2(1 + k) -1    (3.91) ≃ r(E) 1 + k E.ΔE E 2 -m 2 0 c 4 (3.92)
where r is the radius of the closed orbit for which the total energy of the particle is E (can be obtained from Eq.3.85 above). Fig. 3.20 shows the good agreement of the turn separation calculation, using the tracking code ZGOUBI or the analytical formula 3.91.

Off-momentum solution for scaling FFAG

• We will investigate the particle change of orbit due to a momentum error in a scaling FFAG: We write:

r = r 0 + x ; p 0 = qB 0 ρ 0 and B = B 0 ( r r 0 ) k . Thus, p = p 0 r r 0 k+1 = p 0 r 0 + x r 0 k+1 = p 0 1 + x r 0 k+1 ; x ≪ r 0 ≃ p 0 1 + (k + 1) x r 0 (3.93)
Thus,

δp p 0 = p -p 0 p 0 ≃ k + 1 r 0 x (3.94)
In conclusion, for a scaling FFAG, the orbit offset due to momentum change is given by:

x = D δp p 0 = r 0 k + 1 δp p 0 (3.95)
This shows that in order to minimize the beam extent due to momentum spread, it is preferable to have the scaling factor k as large as possible, within the transverse optics stability limits. Also, based on the previous result, it seems likely that there will be an interplay between the dispersion and the space charge effects: D ∝ 1/(k + 1) ≈ 1/ν 2 x . This would even imply that a fully depressed beam (ν x = 0) would have a dispersion function that goes to infinity.

Damping law in FFAGs

In order to investigate the damping law in FFAGs, we compare the solution of the non-linear equation of motion from ZGOUBI tracking with the solution of the linearized particle equations of motion based on the WKB approximation. We start from Eq.(3.26) where we neglect the changes in γ b β b . Substituting the expression of the field index n into it yields,

|h 2 (s)| = |n| ρ 2 = k(R) ρR = k(R) R 2 µ(R) ; µ(R) = R ρ ; vertical plane |h 2 (s)| = |1 -n| ρ 2 = k(R)µ(R) + µ 2 (R) R 2 ; horizontal plane (3.96)
For a scaling FFAG, it was shown earlier that the average scaling factor as well as the index of similarity of the orbits are constants. Thus, the only parameter that changes in the above expressions of h 2 is the radius of the accelerated orbit (for a synchrotron R is fixed). From Eq.(3.35) and (3.36), it results that, for a scaling FFAG, the damping of the amplitude of particle oscillations in the transverse plane is given by:

x ∝ √ R √ β b γ b (3.97) x ′ ∝ 1 √ R 1 √ β b γ b (3.98)
so that,

r x = x 2 1/2 ∝ √ R √ β b γ b (3.99) r ′ x = x ′2 1/2 ∝ 1 √ R 1 √ β b γ b (3.100) ǫ rms = h x 2 x ′2 - xx ′ 2 i 1/2 ∝ 1 β b γ b (3.101)
We generate a 2D mid-plane field-map for a scaling FFAG based on the analytical expression of the field and compare the results of tracking with the ZGOUBI "FFAG" analytical magnetic field model available in ZGOUBI. Both results agree. Therefore, in what follows, we only display the results of the analytical model. The first test that was carried out was for a single particle trajectory in the vertical plane. The results shown in Fig. 3.21 showed excellent agreement with those predicted by the WKB approximation. The comparison was carried out for the beam edge radii and the rms-emittances in the transverse plane as well (see Figs. [START_REF] Holt | Civilian nuclear waste disposal[END_REF].22 and 3.23 respectively) and the results are again in good agreement.

In conclusion, solving the linear equation of motion provides an accurate solution to the non-linear beam dynamics in scaling FFAGs.

Remark: if the orbits do not scale, the index of similarity of the orbits as well as the scaling factor change with the radius and have to be taken into account in evaluating the damping of the particle position and angle in the phase space. In that case, the damping in the horizontal and vertical plane is no longer the same due to the presence of the centrifugal force term (1/ρ 2 ) in the horizontal plane which changes with the radius.

Conclusion

The benchmarking work performed so far helped validate the simulation tools as well as establish the required level of accuracy in the simulation for good convergence of the calculated results. It is in particular found that the approximate solution of the linearized equations of motion is in good agreement with that of the non-linear equations based on Taylor expansions of the field and its derivatives up to the 5th order. Furthermore, analysis of the tune variations showed that the latter are non-negligible so that the Symon formula for the tunes does not apply. Thus, the poor overall transmission (see Table 1.3) may result from the crossing of the betatron resonances. In order to explain the origin of these defects, analysis of the scaling factor were undertaken which revealed some radial and azimuthal imperfections of the field. Some definitions were introduced to account for these imperfections and a detailed beam stability analysis is presented in the next chapter. I am a great believer in the simplicity of things and as you probably know I am inclined to hang on to broad & simple ideas like grim death until evidence is too strong for my tenacity.

-Sir Ernest Rutherford.

Chapter 4

Beam stability analysis and advanced FFAG concept

In the previous chapter, we have shown that certain types of field imperfections that violate the scaling law of FFAGs are important and need to be taken into account in the simulation. In particular, scaling imperfections may explain the important tune variations that were observed experimentally. In this chapter, we develop analytical models that allow an understanding of the effects of radial and azimuthal field defects on the beam dynamics of scaling FFAG. A key finding is that the Symon formula that describes the tunes in FFAGs does not hold in the general case where the orbits do not scale. In our inquiry, we established a new formula that explains the tune variations both qualitatively and quantitatively. This enabled us to generalize the concept of a scaling FFAG to a non-scaling one for which the tune variations obey a well defined law. This result allows to define a way to control the tune variations in FFAGs and therefore provides a correction scheme to scaling imperfections. Several schemes on how to remediate the problem are presented. Furthermore, one showed that achieving the scaling law is not a necessary condition in order to obtain a fixed tune FFAG. Instead, one can cancel the tune variations in FFAGs by alternating the sign of the difference of the scaling factors of the Focusing and Defocusing magnets (k F -k D ), which is the cornerstone of the novel FFAG concept that we present. The significance of these results is also for high power applications where the tune shift induced by space charge forces can be controlled.

Tune calculation and beam stability analysis

Technically, it is impossible to make a field which corresponds exactly to the designed one. Therefore, it is important to understand the effect of small imperfections of the field on the beam dynamics. The general equations of motion including non-linear terms and imperfections are defined by:

           d 2 x ds 2 = P (x, y, s) d 2 y ds 2 = Q(x, y, s) (4.1)
where x and y represent the deviation around the closed orbit, and P and Q are real analytic functions of x, y and s.

In order to better take into account the effect of scaling imperfections, an extension of the mean field index k of the magnets, as defined by Symon, is introduced: the mean field index of the focusing (F ), defocusing (D) magnet and the drift (drif t) are defined in the following way:

k i = R B i dB i dR ; i = F, D, drif t (4.2)
where B i is the vertical component of the magnetic field in the median plane of the FFAG, averaged over the width of the element.

In the following, we investigate the stability of the particle trajectories that can arise with different possible pairs of scaling factors (k F , k D ). We use three different approaches to investigate the beam stability due to field errors: the first approach is based on the hard edge model: we assume that the magnetic field drops abruptly to zero at the edges of the magnet which greatly simplifies the analysis and allows the development of a geometrical model of the lattice. The second model is the BKM average method which we use to seek an approximate solution of the particle equations of motion. The last and third model is the ZGOUBI tracking model of the magnet, which is the most accurate one. The ZGOUBI model solves the non-linear equation of motion using field maps or user-implemented analytical models. To conclude, we establish a comparison between the different results and comment on the outcomes.

Hard edge model

In this first approach, we restrict our analysis to a linear motion around the closed orbits. Therefore, Eqs. (4.1) simplify to:

           d 2 x ds 2 + K x (s)x = 0 ; K x = 1 -n ρ 2 d 2 y ds 2 + K y (s)y = 0 ; K y = n ρ 2 (4.3)
which is the set of linearized equations with periodic coefficients valid in a small neighborhood of the equilibrium orbit. There is a different linearized set of equations for each closed orbit, i.e different pairs (K x , K y ) so that a solution of the above equations may not be generalized to the whole system. For the following stability analysis, we prefer to use the term "stability in the vicinity of the closed orbit", to recall that a linear system is only a first order approach to the general problem. However, for small perturbations, the stability of a non-linear system may often be inferred from the stability of its associated linear system. To make the calculations tractable, we assume that the field inside each magnet is constant with respect to the azimuth and drops abruptly to zero at the edges as illustrated in Fig. 4.1. Also, we keep the radial dependence of the field (k is R-dependent) and allow the scaling law to be different between the F and D-magnet. We start from the field index n where we neglect the scalloping of the orbits,

n = - ρ B dB dx ≈ - ρ B dB dR = - ρ R k (4.4)
Now, we use the index i to make the distinction between the F and D magnet and authorize the field index to change with the radius. This yields,

n i (R i ) = - ρ i R i k i (R i ) ; i = F, D, drif t (4.5)
In order to solve the equation of motion, we use a piecewise method that consists in following the particle motion through a series of elements defined by their transfer matrices (Focusing/Defocusing magnets and drifts). This requires the development of the geometrical model to obtain the relationship between ρ i and R i , therefore the field index n i . For the focusing magnet, K = K x,y > 0 :

M foc =   cos( √ Ks) 1 √ K sin( √ Ks) - √ K sin( √ Ks) cos( √ Ks)   (4.6)
For the defocusing magnet, K = K x,y < 0 :

M defoc =   cosh( p |K|s) 1 p |K| sinh( p |K|s) p |K| sinh( p |K|s) cosh( p |K|s)   (4.7)
Both the horizontal and vertical transfer matrices for half the sector ( \ ( # » OE, # » OI) in Fig. 4.2) can then be computed:

M x sector/2 = M x F/2 × M x edge,F × M drif t1 × M x edge1,D × M x D × M x edge2,D × M drif t/2 M y sector/2 = M y F/2 × M y edge,F × M drif t1 × M y edge1,D × M y D × M y edge2,D × M drif t/2 (4.8)
The number of betatron oscillations per cell ν x,y can then be calculated using:

cos (2π.ν x,y ) = 1 2 T r (M x,y sector ) (4.9)
and a condition sine qua non for the stability of the transverse motion is:

|T r (M x,y sector ) | < 2 (4.10)
Now, we set the geometrical model of the DFD triplet (see Fig.

ρ F = O F A = O F I = O F B ; ρ D = O D C = O D D ; OC 6 = OD R F = OI ; R D = OC ; L B = BC ; L drif t = DF = 2DE (BC) ⊥ (O F B) ; (BC) ⊥ (O D C) ; (O D D) ⊥ (DE) α F 2 + α B + α D + α drif t 2 = π N ; θ F -2θ D = 2π N ; θ F , θ D > 0. 1 (4.11)
Applying the sinus law in the triangle OO F B yields:

ρ F R F = sin α F 2 sin α F 2 + sin θ F -α F 2 (4.12)
Applying the sinus law in the triangles OBI, O F BI and OBC yields:

L B ρ F = L B OB × OB ρ F = sin (α B ) cos θ F -α F 2 -α B × sin θ F 2 sin α F 2 (4.13)
Applying the sinus law in the triangle OBC yields:

R D L B = cos θ F -α F 2 sin(α B ) (4.14)
Applying the sinus law in the triangles OCD and O D CD yields:

ρ D R D = - sin (α D ) cos α D + α B + θ D -θ F + α F 2 × cos θ D 2 sin(θ D ) (4.15)
Applying the sinus law in the triangles ODE and OCD yields:

L drif t R D = L drif t OD × OD R D = 2 sin α drif t 2 × cos α B + θ D -θ F + α F 2 cos α D + α B + θ D -θ F + α F 2 (4.16)
Some remarks:

•The second cardinal condition of the similarity (see Eq (3.60)) implies that the ratio ρ i /R i is constant and shall be the same for different closed orbits. This necessarily implies that the bending angles θ F and θ D are constant and all the relations established above as well.

•Whenever the scaling factors k F and k D are different, the similarity of the orbits is no longer satisfied. To show this, let's proceed using a reductio ad absurdum: let's suppose that ρ i /R i = const. Then, from Eq (4.12) to (4.15), it results that ρ F /ρ D = const and R F /R D = const. However, if we write,

               Bρ = p q = B F ρ F = B D ρ D B F (R F ) = B F 0 R F R 0 k F B D (R D ) = B D0 R D R 0 k D (4.17) 1 \ ( # » OI, # » OE) = π/N = \ ( # » OF I, # » OF B) + \ ( # » ODC, # » ODD) = θF /2 -θD ; θF , θD > 0.
this yields,

ρ F ρ D = B D B F = B D0 B F 0 × R k D D R k F F × R k F -k D 0 ; R D ∝ R F ∝ R k D -k F F (4.18)
which is R-dependent if k F 6 = k D . This is in contradiction with ρ F /ρ D = const, quod erat demonstrandum. Now, we distinguish two cases:

1 st case: we let the scaling factors k F and k D be arbitrary functions of the radius. Yet, we assume that the similarity of the orbits remains, i.e ρ i /R i ≈ const. From Eq (4.5) it results that n

F (R F ) ∝ k F (R F ) and n D (R D ) ∝ k D (R D )
. This is generally a good approximation if k F and k D are slightly different.

Based on that, the transfer matrices of the cell can be obtained and the tunes computed using Eq (4.9).

In the original paper of Symon in 1956, it was shown that the number of betatron oscillations per turn obtained from the smooth approximation are given by [START_REF] Symon | Fixed-Field Alternating-Gradient Particle Accelerators[END_REF]:

ν 2 x ≈ k + 1 (4.19) ν 2 y ≈ -k + f 2 2 (4.20)
where k is the average scaling factor and f the flutter factor that are supposed to remain constant. However, this does not take into account any scaling imperfection. Therefore, one decided to extend the previous result by making the following conjecture:

ν 2 x = x 1 k F + x 2 k D + (x 3 k F + x 4 k D + x 5 ) 2 (4.21) ν 2 y = y 1 k F + y 2 k D + y 3 (4.22)
where (x 1 , x 2 , x 3 , x 4 , x 5 ) and (y 1 , y 2 , y 3 ) are constants that are strictly related to the geometry of the cell. If the conjecture is true, then solving a system of 5 non-linear equations in the horizontal plane, and a system of 3 linear equations in the vertical plane, will provide for a given cell the horizontal and vertical tunes, respectively. For instance, we solve the following system in the vertical plane,

     k F 1 y 1 + k D1 y 2 + y 3 = ν 2 y1 k F 2 y 1 + k D2 y 2 + y 3 = ν 2 y2 k F 3 y 1 + k D3 y 2 + y 3 = ν 2 y3
where (k F i , k Di ) i=1,3 are provided to the hard edge model implemented in Mathematica and (ν yi ) i=1,3 are calculated using Eq (4.9). This yields (y 1 , y 2 , y 3 ) = (-1.1777, 0.5859, 9.1117). For the horizontal plane, we proceed in the same way and obtain (x 1 , x 2 , x 3 , x 4 , x 5 ) = (2.4990, -0.4833, -0.5262, -0.0121, 3.1184). The results are finally summarized in Fig. In conclusion, for small perturbations, the number of betatron oscillations in the hard edge model of FFAG can be approximated by the formula 4.21 and 4.22 above.

2 nd case: we let the scaling factors be arbitrary functions of the radius. Yet, we drop the assumption of the similarity of the closed orbits. This implies that the bending angles θ F and θ D are functions of the radius. We have 6 variables in the system (ρ F , ρ D , R F , R D , θ F , θ D ) requiring a system of 5 equations in order to express the solutions as a function of R F . From the previous results, we set the system:

                                               ρ F R F = sin α F 2 sin α F 2 + sin θ F -α F 2 ρ D R D = - sin (α D ) cos α D + α B + θ D -θ F + α F 2 × cos θ D 2 sin(θ D ) R D R F = R D L B × L B ρ F × ρ F R F ρ F ρ D = B D0 B F 0 × R k D D R k F F × R k F -k D 0 θ F -2θ D = 2π N ; θ F,D � 0 (4.23)
The first 3 equations result from the hard edge model. The 4th equation expresses the change of the bending radii with the particle momentum while the last one expresses the periodicity of the equilibrium orbit, with 2π/N deflection in one sector. This system is highly non-linear and we do not solve it here. Instead, the Bogoliubov method of averages is invoked in the next section.

Method of averages

Using the Bogoliubov-Krilov-Mitropolsky (BKM)'s method of averages [START_REF] Bogoliubov | Asymptotic methods in the theory of nonlinear oscillations[END_REF], one can compute approximately the frequencies of betatron oscillations and their dependence on the average field index of the F and D magnets. This method is characterized by the expansion of the dependent variables in power series of a small parameter as detailed below.

We start from Eqs. (??) where the damping term is removed, and re-write the equation of motion for each closed orbit in cylindrical coordinates (R, θ): differentiating with respect to the azimuthal angle θ, one obtains:

dx dθ = dx ds ds dθ (4.24) d 2 x dθ 2 = d 2 x ds 2 ds dθ 2 + dx ds d 2 s dθ 2 (4.25)
Given that the longitudinal position s is quasi-proportional to the azimuthal angle, one can neglect the 2 nd term derivative (d 2 s/dθ 2 ≈ 0). Also, the arclength in cylindrical coordinates is given by:

ds 2 = dR 2 + R 2 dθ 2 + dz 2 (4.26)
Thus, assuming that the motion is restricted to the median plane, this yields:

ds dθ = r R 2 + ( dR dθ ) 2 ≈ R (4.27)
so that the equations of motion write:

     d 2 x dθ 2 + µ 2 (R, θ)[1 -n] x = 0 ; µ(R, θ) = R ρ d 2 y dθ 2 + µ 2 (R, θ)n y = 0 (4.28)
Now, Eqs. (4.28) can be written in cylindrical coordinates, in the standard form:

d 2 x dθ 2 + g(R, θ)x = 0 (4.29)
where R is a parameter that refers to the local radius of the closed orbit at a given azimuth θ and g(R, θ) is periodic in θ with period 2π/N . In order to apply the averaging method, we introduce the two dimensional vector X = x dx/dΘ as well as the independent variable Θ = N θ. Then, the above equation can be written in the matrix form:

dX dΘ = 1 N GX (4.30)
where

G = 0 1 -g(R, Θ) 0 (4.31)
The main idea of the averaging method is based on the assumption that the derivative dX/dΘ is small. Therefore, one can represent X as the sum of two terms, a slowly varying term ξ, and small rapidly oscillating terms. If N is large (N = 12 in our case), this is generally the case and a solution of (4.30) can be obtained in the form of power series of 1/N [START_REF] Kolomensky | THEORY OF CYCLIC ACCELERATORS[END_REF]:

X(R, Θ) = 1 + 1 N G + 1 N 2 G G -G�G� + ... ξ (4.32)
where ∼ is the tilde operator defined by:

G(R, Θ) = Z [G(R, Θ) -hG(R, Θ)i] dΘ (4.33)
and hi is the average operator defined by:

hG(R, Θ)i = lim T →∞ 1 T Z T 0 G(R, Θ)dΘ (4.34)
Now, a solution of Eq. (4.32) is described by the formula

x = Ae iνθ 1 - 1 N 2 g + ... (4.35)
where ν is an approximate solution of the number of betratron oscillations per turn and is given by:

ν 2 (R) = hgi + 1 N 2 g2 + ... ≈ hg(R, θ)i + * Z {g(R, θ) -hg(R, θ)i}dθ 2 + = g 1 (R) + g 2 (R) (4.36)
From the geometry of the equilibrium orbit and the approximation made in Eq. ( 4.27), it can be shown that:

ds dθ = R cos(φ) ≈ R ∂θ ∂x = sin(φ) R ≈ 0 (4.37) ∂R ∂x = cos(φ) ≈ 1
where φ = \ (x, R) ∼ 0 is the angle between the radial vector R and the Serret-Frenet vector x. Now, we calculate the field index n:

n = - ρ B ∂B ∂x = - ρ B ∂B ∂R ∂R ∂x + ∂B ∂θ ∂θ ∂x ≈ - ρ R k(R) (4.38)
Now, back to Eqs. (3.58). The preceding formalism applies by simply replacing g(R, θ) by2 :

       g x (R, θ) = R ρ 2 × (1 -n) ≈ R ρ 2 × h 1 + ρ R k(R) i g y (R, θ) = R ρ 2 × n ≈ - R ρ × k(R) (4.39)
We calculate each term in Eq. (4.36) separately by introducing the azimuthal variation of the scaling factor in the following way:

k(R, θ) =      k F (R), if θ ∈ θ F k D (R), if θ ∈ θ D k drif t (R), if θ ∈ θ drif t (4.40)
We start from the vertical plane and calculate the first order term which is the mean value of the forcing term in Eq. (4.29):

g 1y (E) = hg y (R, θ)i = -1 2π/N Z 2π/N 0 R ρ × k(R, θ)dθ = -1 2π/N " Z θ F µ(R, θ)dθ k F (R) + Z θ D µ(R, θ)dθ k D (R) + Z θ drif t µ(R, θ)dθ ! k drif t (R) # = X i α i (E) × k i (E) ; i = F, D, drif t (4.41)
where E is the kinetic energy of the particle and α i (E) is what we will refer to in the future as "the 1 st order index of similarity of the closed orbits":

α i (E) = -1 2π/N Z θ i µ(R, θ)dθ = -1 2π/N Z θ i R ρ dθ (4.42)
The assumption on the similarity of the closed orbits is withdrawn here given that α i is allowed to change with the radius. Also, note that for the F magnet, R/ρ>0, so that α F < 0. However, for the D magnet, α D > 0. Therefore, by means of Eq. (4.41), one can predict that, to the first order in the method of averages, the D magnet increases the vertical focusing, while the F magnet plays the opposite role. Also, the drift is included in this formula. Since the drift space between the magnets is susceptible to contain the fringe fields, it may be important in some cases to assign a scaling factor to it to determine its effect on the beam dynamics. However, in the ideal case, ρ drif t → ∞ so that α drif t ≈ 0.

The second order term can be calculated using:

g 2y (E) = * Z {g y (R, θ) -hg y (R, θ)i}dθ 2 + = * Z {g y (R, θ) -g 1y (E)}dθ 2 +
The azimuthal variation of the field yields an additional positive edge focusing effect (not calculated here), so that, to the first order, the vertical tune writes:

ν 2 y (E) = P i α i (E) × k i (E) + F 2 1 + 2 tan 2 (ξ) (4.43)
where ξ is the spiral angle (see Appendix 5.6) and F is the magnetic flutter (not the flutter function F ) defined by:

F 2 = B 2 -hBi 2 hBi 2 = F 2 -hF i 2 hF i 2 (4.44)
In general F is E-dependent.

We repeat the same steps in the horizontal plane by calculating the first order term

g 1x (E) = hg x (R, θ)i = 1 2π/N Z 2π/N 0 " R ρ 2 + R ρ × k(R, θ) # dθ = 1 2π/N Z 2π/N 0 µ(R, θ) 2 dθ -g 1y (E) = X i β i (E) - X i α i (E) × k i (E) ; i = F, D, drif t (4.45)
where β i (E) is what we will refer to in the future as "the 2 nd order index of similarity of the closed orbits":

β i (E) = 1 2π/N Z θ i µ(R, θ) 2 dθ = 1 2π/N Z θ i R ρ 2 dθ (4.46)
Thus, to the first order in the method of averages, the number of betatron oscillations in the horizontal plane can be approximated by:

ν 2 x (E) = P i β i (E) - P i α i (E) × k i (E) (4.47)
This result is general and applies to any type of FFRA because no assumption is made regarding the analytical form of the magnetic field except that it can be written in the general form of Eq.(3.81) . The idea is that any magnet can be decomposed into a sum of smaller magnets for which one can attribute an average field index that is constant inside (no azimuthal dependence) yet can change with the radius. In other words, if one knows precisely the closed orbits from injection to extraction and the average field index associated to each magnet, then the number of betatron oscillations can be determined by means of the previous results.

Application example: the KURRI 150 MeV FFAG ring:

According to the previous results, the number of betatron oscillations per turn in the KURRI DFD triplet is given, to the first order, by the following formula:

ν 2 x (E) = β F (E) + 2β D (E) + β drif t (E) -α F (E)k F (E) -2α D (E)k D (E) -α drif t (E)k drif t (E) (4.48) ν 2 y (E) = α F (E)k F (E) + 2α D (E)k D (E) + α drif t (E)k drif t (E) + F 2 (4.49)
where E refers to the particle kinetic energy. We calculate the number of betatron oscillations per turn using the tracking code Zgoubi and the TOSCA field map of the DFD triplet. As was shown earlier from benchmarking, the results are in agreement with the other codes (see Fig. 3.5). Next, the tracking results provide the closed orbits for various particle energies which were then exploited to calculate the horizontal and vertical tunes by applying the 1st order Bogoliubov method of averages (see Eqs. The Symon formula is less in agreement with the benchmarking results than the other calculations. There are two reasons for such a disagreement: first, the Symon model assumes that the scaling factor is constant and is azimuthally invariant. Such an assumption is withdrawn in our calculation given that we define an average field index for each magnet (azimuthal variation) that changes with the radius as well. This gives rise to an orbit distortion that is accounted for by the first and second order index of similarity of the closed orbits. These quantities exhibit strong dependence on the radius and provide an additional focusing in the horizontal plane, which is the second reason for the discrepancy. Fig. 4.6 shows how the second order index of similarity (β(E) = β F (E) + 2β D (E) + β drif t (E)) changes with the particle energy.

ZGOUBI model

For a scaling FFAG, an important quantity is the tune which is supposed to remain constant over the entire range of energies. From the hard edge model and the Bogoliubov method of averages, it was found that the tune is sensitive to the average field index k F and k D of the F and D magnet respectively. In other words, breaking the scaling law, although a major source of imperfection in scaling FFAG, can also be utilized in order to control the tune path in FFAG. In order to quantify the source of imperfection, we introduce two new quantities in the calculation: the average value of the tunes ν m x = hν x i and ν m y = hν y i over the range of energies to quantify the average focusing strength of the applied forces on the beam, and the rms value of the tunes ν rms x = σ νx and ν rms y = σ νy to quantify the scaling imperfections in terms of tune variations3 .

Given that the energy gain per turn is small (2 keV per turn), one can reasonably assume that the accelerated particle trajectory for a given energy is quasi the same as the closed orbit trajectory. This is illustrated in Annex 5.6 where a comparison of the calculation of the horizontal and vertical tunes using the closed orbits as well as the accelerated orbit shows good agreement. There are two criteria to ensure the stability of the lattice: the first is to find all closed orbits between injection energy and extraction (or final) energy. In our case, we limit this study to E inj = 11 MeV and E ext = 100 MeV, suitable for the KURRI FFAG. The second criteria is that the number of betatron oscillations of the particles, which is obtained from multi-turn tracking or from a first-order mapping around each of the closed orbits, be strictly positive. The first criteria is justified by the fact that a solution to the system of Eqs (4.23) does not always exist. The second criteria ensures that the motion is stable around each closed orbit.

In summary, the procedure is detailed below:

1. Generate a field map for a given (x, k F , k D ). The field can be written in cylindrical coordinates as:

B(R, θ) = B F 0 × R R 0 k F × F F (θ) + x × B D0 × R R 0 k D × F D (θ) (4.50)
where F F (θ) and F D (θ) are the flutter functions describing the azimuthal dependence of the field and equal to zero in the D and F magnet, respectively, and x is a scale factor that allows to vary the FD ratio of the magnetic field (x � 0). If k F = k D , the field writes in the standard form of Eq. (3.63). Note that the field is separable in radial and azimuthal coordinates since the flutter functions merge to zero in the drift space between the F and D magnets as can be seen in Fig. 2. Search for NCO closed orbits between injection and extraction using the built-in fitting routines in ZGOUBI. NCO was chosen to be 30 in order to have good statistics and ensure the convergence of the calculated quantities. A typical example of 4 closed orbits search is illustrated in Fig. (4.8).

θ (rad) B D (R) = B D0 (R/R 0 ) k D B F (R) = B F0 (R/R 0 ) k F

For each closed orbit, the number of betatron oscillations in both planes is computed using a Discrete

Fourier Transform (DFT), from which the average as well as the rms values of the tunes are deduced:

ν m x,y = hν x,y i = 1 N CO N CO X i=1 ν x,y,i (4.51)
ν rms x,y = σ νx,y = 1 N CO N CO X i=1 � ν x,y,i -ν m x,y 2 ! 1/2 (4.52)
The merit of the DFT transform is that the assumption of the linear motion made by using Eq (4.9) is now withdrawn. This method simply relies on the particle coordinates stored turn after turn from particle tracking in ZGOUBI.

The tracking code Zgoubi that we use for our analysis, solves the non-linear equation of motion4 using truncated Taylor expansions of the field and its derivatives up to the 5th order. Therefore it is more accurate than the linear approach.

Stability analysis based on the average field index k and the F D ratio of the magnets:

First, one investigates the stability of the lattice by varying the F D ratio (thus the scale factor x) and the average field index k of the magnets (k F = k D = k): a scan on k and x provides the stability diagrams below: Fig. 4.9 displays the results. As expected, the RMS tune variations are negligible and all lattices have a constant tune. Also, as can be observed, increasing the average field index of the magnets increases the horizontal cell tune and lowers the vertical one. It can also be seen that increasing the F D ratio allows to increase the tunes in both planes. The stability diagram showing both ν x and ν y dependence on the above parameters can be seen in Fig. 4.10: it can be observed that the F D ratio has a more profound infuence on the vertical cell tune than in the horizontal one (the tune changes are twice more important in the vertical than in the horizontal plane). However, the dependence on k is almost the same (although opposite). Furthermore, one can observe that, on the top and bottom right, the stability limits are set by the horizontal and vertical cell tunes, respectively. On the left side, the physical size of the magnets (here we choose a radial excursion limited to 10m) determines the boundary limits. A more complete stability diagram where the radial extent of the magnets, R max , is no longer a limitation was calculated and is shown in the appendix 5.6. Stability analysis based on the average field index of the F and D magnets:

� � � � � �� �� �� �� ��� ��� ��� ��� ��� ��� ��� -� �� �� �� * � ��� �� � �� � � � �� � � � �� � � ��� �� � �� � � ��� �� �� �� � �� �� � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ����
� � � � � �� �� �� �� ��� ��� ��� ��� ��� ��� ��� -� �� �� �� * � � � � � � � �� � � � �� � � � �� � � � �� � � � � � � � � � � � � � �� � � � �� � � � �� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ����������������� ���������� ��������
Now, we choose to focus our analysis on the average field index of the magnets. For that, we fix the FD ratio of the magnets. We choose x = 1 which corresponds to the green line in Fig. One can also observe that for high k values, the stability diagram shrinks, thus any design imperfection will make the orbits quickly unstable 5 . This is expected since the second order term in the field expansion becomes non-negligible so that the system is no longer slowly (linearly) responding to perturbations: This study revealed some interesting results: although the rms tune variations become quite sensitive to scaling imperfections near the boundary of stability, it is observed that in the vicinity of the central line, these quantities are mostly proportional to the difference of the scaling factors as shown in Fig. (4.15) and can be written in the form:

B(R) = B 0 × ( R R 0 ) k = B 0 × (1 + x R 0 ) k ; x ≪ R 0 = B 0 × 1 + k x R 0 + k(k -1)( x R 0 ) 2 + ... + k(k -1)...(k -n + 1) n! x R 0 n + O ( x R 0 ) n+1 � � � � � �� �� � � � � � �� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ����� � �� � � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � � �� � � � �� � � �� �� � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ����� � �� � � � �� � � � � � � � ����� � �� � � � �� � � ��� �� �� �� � � �� � � � �� � � � ��
� � � � � �� �� � � � � � �� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � �� �� � � �� � � � �� � � �� �� � � �� � � � �� � � �� �� � � �� � � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � � �� � � ���� � �� �� � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � ��� �� � �� � � � �� � � � �� � � �� �� � � �� � � ����� ���
� � � � �� �� � � � � � �� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ���� � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ����� � �� � � � ��
� � � � � �� �� � � � � � �� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � � �� � � �� �� � � �� � � � �� � � � �� � � �� �� � ��� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � �� �� � � �� � � � �� � � ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � ��
ν rms x,y ≈ a x,y |k F -k D | = a x,y |κ| ; κ = k D -k F (4.53)
This result can also be inferred from the contour plot (see Figs. ) . We will exploit this result later on when we establish a way to counteract the rms tune variations introduced by the space charge effects in scaling FFAG, the idea being to perturb k F or k D . Another important result that this study revealed is that the vertical tune exhibits the same behavior vis à vis k F and k D , however, for the horizontal tune, the situation is strongly related to whether k D > k F or the opposite (see Fig. (4.14)). This will be discussed in more detail in the next section. The implication of that is essentially, by implementing trim coils, as in cyclotrons, one can play with the tune in a flexible way: modifying k F or k D or k F and k D at the same time, would allow to change the tune path in both planes and avoid crossing any resonances.

� � � � � �� �� � � � � � �� �� ����� � � � � � � �� � � � �� � � � �� � � � � � � � � �� � � � � � � � � �� � � ����� � � � � � � �� � � � � � � � � �� � � ���� � � � � � � � �� � � �� �� � � � � � � � � � � � � �� � � � � � � � � � � � � � � � � � � �� � � � � � � � � � � � � � �� � � � � � � � � �� � � � �� � � � � � � � � �� � � � �� � � � � � � � � �� � � � � � � � � �� � � � ��
Finally, this study showed that the stability boundary on the left side, i.e when κ = k D -k F > 0, is limited by the radial π-mode stop-band resonance. This defines the termination to the growth of the horizontal cell tune which cannot exceed 0.5, since then, the condition (4.10) is violated. On the right side, i.e when κ < 0, the stability limit is set by the condition that ν 2 y is to remain positive (Floquet resonance). Recalling the above established Eq. ( 4 must exceed the mean field index for all energies, which is caused by the sharp increase of the field with the radius (B ∝ R k ). However, these boundaries can be overcome and larger κ-values explored in the novel FFAG concept that we present later on in this chapter.

Discussion

From the results of the hard edge model, it was shown that, under the assumption that the similarity of the orbits is maintained, the square of the number of betatron oscillations in both planes is, to the first order, a linear combination of the average field indices k F and k D of the focusing and defocusing magnets respectively. Later on, when this assumption was dropped, and, solving the linearized particle equations of motion by means of the Bogoliubov method of averages, it was found that this result can be generalized by defining the 1 st and 2 nd order indices of similarity of the closed orbits α i and β i respectively. According to this method, the effect of perturbing the field indices of the magnets may produce an orbit distortion that manifests through a radial dependence of the index of similarity of the orbits that adds to the radial dependence of the scaling factors. Furthermore, the magnetic flutter is susceptible to change when the scaling law is perturbed. This changes the tune path in both planes. Therefore, in this section we compare the results of solving the non-linear equation of motion using the tracking code Zgoubi with those predicted by the Bogoliubov analytical calculation based on the linearized particle equations of motion. In what follows, we exploit the results of the Zgoubi model based on the procedure that was detailed in the previous section (see Eq. (4.50)). The first test that was carried out to characterize the similarity of the orbits is to calculate the 1 st order index of similarity of the orbits averaged over the width of the DFD triplet. As can be seen in Fig. 4.16, if k F = k D , then the index of similarity is energy-independent. This is valid within the focusing and defocusing magnet as well (see Fig. 4.17) so that the second cardinal condition in Eq. (3.60) is satisfied. We will refer to the case where all parameters x scale, i.e k F = k D = k ref as x ref . Now, we introduce a perturbation of the average field indices of the magnets by perturbing k D and maintaining k F to a fixed value. Two regimes are observed:

•If κ < 0, i.e k D < k F (down-pointing triangle in the figures) then: α F < α ref F and α D > α ref D . Besides α F (E) is strictly decreasing while α D (E) is strictly increasing. •If κ > 0, i.e k D > k F (upward-pointing triangle in the figures) then: α F > α ref F and α D < α ref D .
Besides α F (E) is strictly increasing while α D (E) is strictly decreasing. Thus, both magnets act in opposition vis-à-vis scaling imperfections. This result is qualitatively in agreement with the hard edge model: solving the system of equations (4.23) in an iterative way yields: The 2 nd order index of similarity of the orbits as calculated from the closed orbits is shown in Fig. 4.19. The trend is similar to the α F term as expected: hR/ρi F ≫ hR/ρi D =⇒ β = hR/ρi 2 ≈ hR/ρi 2 F . The β term is particularly important as it represents the additional focusing provided by the horizontal restoring force. In order to better understand the importance of this term and the differences between the scaling and the non-scaling FFAG case6 , we compare the particle orbits for two samples (see Fig. 4.18): for the scaling FFAG case, the orbits are parallel to each other (their distance at various azimuthal positions does not change (δR = const)). However, for the non-scaling FFAG case (k F 6 = k D ), the distance increase between the orbits is essentially dominant at the centre of the F-magnet and grows with the radius. This justifies the importance of the indices α i and β i to characterize the particle orbit.

If α F (E) = R ρ
Next, the magnetic flutter (Fig. 4.21) as well as the average scaling factor (Fig. 4.20) of the DFD triplet were computed. Notice that the average scaling factor remains nearly constant (Δk/k < 1%) when scaling imperfections are introduced (k F 6 = k D ) so that the Symon formula does not apply to explain the origin of the RMS tune variations. This also shows that having an average scaling factor k constant is not a sufficient condition to obtain a fixed tune machine. Based on all the above calculation, we compute the tunes using the Bogoliubov method of averages and compare its outcomes with the Zgoubi calculation which is based on multi-turn tracking and DFT. The results in both planes are shown in Figs. [START_REF]Nuclear Fuel Cycle Information System[END_REF].22 and 4.23: both results exhibit the same behavior and the agreement is quite good except near the boundary of stability. Also it is observed that the tune variations in the vertical plane are mostly due to the variations of the magnetic flutter. However, in the horizontal plane, the index of similarity β carries most of the variations. All the above parameters characterize the transverse beam dynamics in an imperfect scaling FFAG. Example of tune calculation using the tracking code ZGOUBI (blue) and comparison with the 1st order analytical formula 4.47 (red). To the first order in the method of averages, good agreement (< 4%) is obtained in the cases (kF , kD) = (7.6, 7.0) and (kF , kD) = (7.6, 7.6). However, for the last case, i.e (kF , kD) = (7.6, 8.2) the large discrepancy observed (� 16%) is due to the fact that the beam is located nearby the boundary of stability. To the first order in the method of averages, the differences are less than 20%. The agreement could be better if the contribution of the second order term were taken into account (given that g2 > 0, see Eq. 4.36).

Proof:

Now, we will establish the following property: Property:"In presence of scaling imperfections, the number of betatron oscillations per turn increases (resp. decreases) with the energy if κ > 0 (resp. κ < 0). Besides, the RMS tune variations are proportional to |κ|, i.e ν rms

x,y ≈ a x,y |κ|."

Earlier on, we have established the above results through numerical simulations (see Figs. 4.22, 4.23 and 4.15). However, now, we will prove the above property: As shown earlier, the tune variations are imposed by the β term in the horizontal plane and the magnetic flutter in the vertical plane. Given that hR/ρi D ≪ hR/ρi F , one obtains,

β = * R ρ 2 + DF D ≈ * R ρ 2 + F = β F (4.54)
This means that, in order to explain the tune variations, one has to focus on explaining the changes of hR/ρi F with the energy of the particle: Let's consider the case of a scaling FFAG, for which k F = k D = k and the associated non-scaling FFAG for which k D is perturbed such that k D > k F , i.e. κ > 0. Now, writing the expression of the average magnetic field of the non-scaling FFAG, one obtains:

hBi ns = A 1 R R 0 k F + A 2 R R 0 k D ; A 1 > 0 & A 2 < 0 = A 0 R R 0 k " 1 + A 3 R R 0 k D -k F # ; A 0 > 0 & A 3 < 0 = hBi s A 4 1 + A 3 R R 0 κ ; A 4 > 0 (4.55)
where ns (resp. s) denotes the non-scaling (resp. scaling) case. Definining the average radius R E of a closed orbit of energy E such that hBi R E = p/q where p is the particle momentum, this yields:

R ns E R s E = hBi s hBi ns = A 4 1 + A 3 R R 0 κ ; A 4 > 0 & A 3 < 0 (4.56)
which is an increasing (resp. decreasing) function of the energy if κ > 0 (resp. κ < 0). Knowing that ρ F satisfies B F ρ F = p/q, one obtains ρ ns F = ρ s F so that

β ns F ≈ R ns E ρ ns F 2 =     A 4 β s F 1 + A 3 R R 0 κ     2 (4.57)
Given that β s F = const, this proves the first part of the property (the same approach can be used to prove the same results for the magnetic flutter). Formulated differently, the above results show that, if κ > 0, then the distance in radius between consecutive closed orbits increases with increasing energy compared to the scaling FFAG case. This is expected given that the average magnetic field hBi ns decreases. Now, the second part of the property: if we write R = R 0 + ΔR where R is the maximum radius of the closed orbits and R 0 is the injection radius, then a Taylor expansion of Eq.(4.57) yields (assuming κ > 0 and

ΔR ≪ R 0 ): ν rms x ≈ |β ns F (max) -β ns F (min)| = (A 4 β s F ) 2 1 + A 3 R R 0 κ 2 - (A 4 β s F ) 2 (1 + A 3 ) 2 (4.58) ≈ 2A 3 (1 + A 3 )(A 4 β s F ) 2 ΔR R 0 κ (1 + A 3 ) 2 1 + A 3 + A 3 κ ΔR R 0 2 ∝ |κ| (4.59)
(except near the boundary of stability) which proves the property above.

Verification test:

We start from Eq.(4.58) in order to compare its results with those obtained from tracking: the conditions of validity remain the same, i.e κ > 0. Also k F = 7.6 = const. Varying the average field index of the D-magnet yields the RMS tune variations per cell in the horizontal plane. Now, we compute the terms that intervene in Eq.(4.58):

R R 0 = R 0 + ΔR R 0 = 1.1363 ; A 3 ≈ -0.5 ; β s F =

Correction scheme and advanced FFAG concept

Correction scheme

In the following, we address the question of how to implement a correction scheme in order to remediate the tune variations of the 150 MeV FFAG. As shown earlier, based on the sign of κ, the tune variation with the energy exhibits antagonistic behavior. Therefore, the idea of the following scheme is to introduce a perturbation of the field, every two sectors in order to counteract the already existing imperfections: to be more explicit, the 12-fold symmetry of the FFAG is reduced to a 6-fold symmetry in the following way: let's note D i (resp F i ) the Defocusing (resp Focusing) magnet with scaling factor k Di (resp k F i ). The idea of the correction system is to replace the original design

D 0 F 0 D 0 -D 0 F 0 D 0 by D 0 F 0 D 0 -D 1 F 0 D 1
where k D1 is chosen in the following way: If k D0 > k F 0 then k D1 < k F 0 and vice-versa (see Fig. 4.25). In other words, achieving a fixed tune FFAG can be achieved by alternating the sign of κ (in consecutive cells). In this case, the phase advance of the combined two sectors (superperiod) will combine the two antagonistic effects of the scaling imperfections. The correction can be achieved by implementing trim coils along the radius of the D-magnet, the F -magnet or both, every two sectors. In order to show the validity of the previous Also, note that these results were not finely optimized. A step size of δk = 0.1 was employed to search for the minimum RMS tune variations, which can be further reduced.

Novel FFAG concept

This novel concept of a non-scaling FFAG allows to obtain a fixed tune machine and can be implemented in any alternating gradient accelerator type. Besides, it shows that the two cardinal conditions (??) and (??) are sufficient but non-necessary conditions to obtain a fixed tune FFAG. This concept is better illustrated in Fig. 4.30: when alternating the sign of κ, the phase advance per cell will alternate between the blue and red curves: this alternation will remain counterbalanced given that the phase advance per cell alternates its monotonic behavior. This yields the green curve.

The above results rely on the fieldmap tracking with an assumed fringe field form as shown in Fig. 4.7. There- fore, in order to verify the previous results for a variety of cases, one choose to vary the fringe field by using the analytical model "FFAG" in ZGOUBI [START_REF] Lemuet | Developements in the ray-tracing code Zgoubi for 6-D multiturn tracking in FFAG rings[END_REF]. As shown in Fig. 4.31, when generating an FFAG model, one paid attention not to overlap the fringe fields of the F and D magnets. In other words, the drift space between the magnets contains no fringe fields. This choice is made to simplify the analysis. However, as shown earlier one can account for such effects by adding k drif t . The approach is similar to what we discussed earlier. Fig. 4.32 shows the results of a scan of several lattices where one varies the average field index of the D-magnet, k D . This helps identify the optimum working point k D1 = 7.748 as summarized in Fig. 4.33: in blue, the tune per cell in a scaling FFAG is shown. The two non-scaling FFAGs with monotonic tune variations are shown in green and purple. Finally, the advanced fixed tune FFAG concept is shown in red. From now on, our analysis is restrained to this case.

First of all, one looked at the closed orbits inside the fixed tune FFAG lattice (Fig. 4.34): as can be seen, the orbits of the alternating κ FFAG do not scale: at any azimuthal position, comparison with the scaling FFAG (curves in black) show that the orbits do not scale, however the combined effect of the two cells seems to (k F ,k D )=(7.6,7.8) (k F ,k D )=(7.6,7.8), (7.6,7.4) (k F ,k D )=(7.6,7.4)

Figure 4.30: Tune variations as a function of the energy before and after correction: the corrected scheme is shown in green where (kF 0, kD0) = (7.6, 7.8) and (kF 1, kD1) = (7.6, 7.4). These results are obtained from multi-particle tracking assuming a Gaussian distribution of the particles in the transverse plane. The errorbars represent the overall tune variation from the distribution. Note that the space charge effects are not included so far in this study. cancel. To verify this, one calculated the 1 st order index of similarity of the orbits within each cell (Fig. 4.35): as expected, the orbits do not scale within each cell, but the scaling is achieved in average. In other words, R/ρ becomes energy-independent when (and only when) averaged over the width of the two DFD triplets.

To further pursue our analysis, one calculated the average field index within each cell, which is shown to increase with the energy (see Fig. 4.36). This, in particular, demonstrates that having an average field index k constant is not a necessary condition to obtain a fixed tune FFAG. In addition, one calculated the magnetic flutter as well as the 2 nd order index of similarity of the orbits (Figs. 4.37 and 4.38) which are in accordance with the expectations: scaling imperfections which manifest through an azimuthal variation of the average field index create a change in the magnetic flutter as well as the closed orbits which can be cancelled as shown here. Such an effect is negligible in synchrotrons for instance, given that the orbit is fixed so that the magnetic flutter may not change in a considerable way.

Conclusion:

Until now, a scaling FFAG concept was based on the second cardinal condition, i.e the assumption of the similarity of the orbits at any azimuthal position. However, according to the concept developed earlier, the similarity of the orbits can be achieved in an average sense by alternating the difference of the average field index of the focusing and defocusing magnets. This demonstrates that the second cardinal condition is a sufficient but not necessary condition in order to obtain a fixed tune FFAG. This allows more flexibility in the design of a fixed tune FFAG: since it is almost impossible to ensure that all orbits are similar at every azimuthal position, it is easier to ensure the similarity in an average sense by creating alternating scaling imperfections (by alternating κ) that maintain a constant magnetic flutter, similar orbits and therefore a constant tune. The condition that allows to achieve a fixed tune FFAG translates in the following way:

∂ ∂p * R ρ 2 + = 0 (4.61)
Another important outcome of this novel concept is that the stability limit (see Fig. 4.11) can be overcome: large κ values can be explored although we will not develop this further in the present dissertation. Closed orbits in a scaling FFAG (black) and in a fixed tune non-scaling FFAG with alternating κ (pink and light blue). For the sake of clarity, the distance between the closed orbits of the scaling and the non-scaling FFAG is amplified.

One particularly important parameter to study is the Dynamic Acceptance (DA) of the FFAG accelerator and how it relates to the field imperfections introduced above. In particular, one needs to evaluate how the concept of alternating the κ-values of the FFAG might impact the DA and how it compares to the DA of the scaling FFAG.

Dynamic Acceptance

To investigate the effects of the non-linear field imperfections on the DA of the FFAG, we construct a numerical model based on the procedure that allows to vary the average field index of the magnets. In our model, the acceptance corresponds to the maximum emittance that the beam can have without loss due to single particle dynamics effects. Particle tracking at fixed energy is used in our analysis: we launch several particles with different initial displacements around each closed orbit that we track over 1000 turns. Once the maximum stable orbit is calculated, one uses the linear approximation to determine the maximum emittance (The beam emittance damps adiabatically during the acceleration. Thus, we define the normalized emittance in our calculation):

ǫ x = βγ β x x 2 (4.62)
To simplify the analysis, we start by varying the average field index of the D-magnet while keeping the same configuration for the F-magnet. The results are shown in Fig. 4.39 where one can observe that increasing the field index of the D-magnet does increase the resulting DA. Next, we construct a heuristic model to explain the qualitative as well as the quantitative behavior of the DA calculation as a function of the average field indices (k F , k D ) of the magnets: We make the bold assumption that the DA can be calculated using a barycentric interpolation where the barycentric basis functions are the DAs defined with respect to k D1 and k D2 (in 1 Dimension) that are randomly chosen:

DA(k F , k D ) = k D2 -k D k D2 -k D1 × DA(k F , k D1 ) + k D -k D1 k D2 -k D1 × DA(k F , k D2 ) (4.63)
It results from the above model that, for a fixed k F -value, the DA of the FFAG has a linear dependence on k D , i.e:

DA(k D ) = DA(k D2 ) -DA(k D1 ) k D2 -k D1 × k D + k D2 k D2 -k D1 DA(k D1 ) - k D1 k D2 -k D1 DA(k D2 ) (4.64)
which shows that the DA is an increasing function of the average field index of the D-magnet. The above conjecture is tested against a variety of cases as illustrated in Fig. 4.40 which shows good agreement between the tracking results and the conjectured values.

Given the successful verification of the above conjecture in 1D, one decided to extend it to the 2D case by using a 2D triangular barycentric interpolation as illustrated in Fig. 4.41: the idea is to use 3 barycentric basis functions: The idea of the interpolation method is to use the areas of the triangles (A1,A2,A3) defined with respect to the working point (k F , k D ) where the DA is to be calculated (the working point need not belong to the triangle). The DA at the interpolation point is given by:

DA1=DA(k F 1 , k D1 ), DA2=DA(k F 2 , k D2 ) and DA3=DA(k F 3 , k D3 ).
DA(k F , k D ) = λ 1 DA(k F 1 , k D1 ) + λ 2 DA(k F 2 , k D2 ) + λ 3 DA(k F 3 , k D3 ) (4.65)
where λi = Ai/(A1 + A2 + A3). The results of the method are illustrated in Fig. 4.42 where one can observe a good agreement between the simulated and the conjectured values. Explicitly, the formulae for the barycentric coordinates at point (k F , k D ) are given by:

A1 DA1 DA2 DA3 A2 A3
                     λ 1 = (k D2 -k D3 )(k F -k F 3 ) + (k F 3 -k F 2 )(k D -k D3 ) (k D2 -k D3 )(k F 1 -k F 3 ) + (k F 3 -k F 2 )(k D1 -k D3 ) λ 2 = (k D3 -k D1 )(k F -k F 3 ) + (k F 1 -k F 3 )(k D -k D3 ) (k D2 -k D3 )(k F 1 -k F 3 ) + (k F 3 -k F 2 )(k D1 -k D3 ) λ 3 = 1 -λ 1 -λ 2 (4.66)
The method described above can also be extended to an N-sided polygon: knowing the discrete values of the DA on the boundary of stability, one can infer the continuous data over the entire domain by simply using the barycentric coordinates transformation.

Case of alternating-κ FFAG

We carry out the same simulations to calculate the DA of an alternating-κ fixed tune FFAG and compare its results with the scaling FFAG case. The results are illustrated in Fig. 4.43: increasing the k D -value increases the DA as shown earlier. However, when alternating κ, one can observe that the DA is lower than the scaling case (pink vs green curve). This is due to the cell with κ < 0.

As a general conclusion, this shows that in a realistic scaling FFAG machine, where the exact gradients of all the magnets are not known precisely, the DA is limited by the magnets for which κ < 0. A comparison of the horizontal phase space trajectories between the scaling FFAG and the alternating-κ fixed tune FFAG is shown in Fig. 4.44: for the latter, the phase space trajectories are symmetric with respect to X ′ = 0: this ensures that the phase advance on each cell (red and green) cancel each other so that the tunes become constant.

Finally, one looked at the effect of increasing κ on the DA of the alternating-κ FFAG (Fig. 4.45). Although large κ values can be explored, the main finding is that the DA decreases with increasing κ. A careful choice of this parameter is therefore needed. If we denote N the number of superperiods in the ring (N =12 for KURRI), then implementing the alternating κ concept will reduce N by a factor of 2 which will double the resonance population. An interesting problem to investigate in the future is therefore the crossing of the resonance lines in both cases and how the topology of the phase space trajectories is changed. (b) Case of alternating-κ fixed tune FFAG: (kF , kD) = (7.6, 7.4), (7.6, 7.8). In green, the phase space trajectories in the 1 st cell, and, in red, the phase space trajectories in the 2 nd cell. The phase space trajectories in both cells are symmetric with respect to the average closed orbits, i.e R = hRi and X ′ = hX ′ i = 0. The light blue lines illustrate the particle trajectories from one cell to the other.

Figure 4.44:

Horizontal phase space at 100 MeV including the separatrix. The limit of stable motion is obtained when the particle is lost in less than 1000 turns. An error-free lattice is assumed.

Conclusion

In this chapter, one analyzed the stability of the particle trajectories due to field errors. Several approaches to the problem were developed and analyzed. Comparison of the results showed that the linear approach is only valid in the vicinity of the closed orbits: all parameters of the hard edge model evolve due to the non-scaling of the orbits. Therefore, one relied on the non-linear approach based on tracking simulations. A crucial result was to establish a relationship between the number of betatron oscillations and the field defects. A key parameter to measure the amplitude of the defects is the κ-value defined as the difference of the average field index of the focusing and defocusing magnets. Based on these results, a new scheme to remediate the variation of the betatron oscillations with the energy was proposed. The main idea consists in alternating the κ-values of the magnets, every two (or more) sectors. This lead to the new concept of the fixed tune non-scaling FFAG that one developed in the last part of this chapter: in addition to the fact that this demonstrates the conditions of scaling are non necessary to obtain a fixed tune FFAG, the newly developed concept is easier to implement by means of trim coils that can be adjusted to find the condition of minimum tune excursion and avoid the (k F ,k D )=(7.6,7.6), (7.6,7.6) i.e κ=0.0 (k F ,k D )=(7.6,7.4), (7.6,7.8) i.e κ=0.2 (k F ,k D )=(7.6,7.2), (7.6,8.0) i.e κ=0.4 (k F ,k D )=(7.6,7.0), (7.6,8.2) i.e κ=0.6 (k F ,k D )=(7.6,6.8), (7.6,8.4) i.e κ=0.8 crossing of harmful resonances. Analysis of the Dynamic Acceptance showed that the lattice with alternating κ-values has lower values than the equivalent scaling FFAG case. This showed that the amplitude of κ has to be well controlled in order to maintain a large beam acceptance. However, this needs further investigation.

Une cause très petite, qui nous échappe, détermine un effet considérable que nous ne pouvons pas ne pas voir, et alors nous disons que cet effet est dû au hasard. Si nous connaissions exactement les lois de la nature et la situation de l'univers à l'instant initial, nous pourrions prédire exactement la situation de ce même univers à un instant ultérieur. Mais, lors même que les lois naturelles n'auraient plus de secret pour nous, nous ne pourrons connaître la situation initiale qu'approximativement. Si cela nous permet de prévoir la situation ultérieure avec la même approximation, c'est tout ce qu'il nous faut, nous disons que le phénomène a été prévu, qu'il est régi par des lois ; mais il n'en est pas toujours ainsi, il peut arriver que de petites différences dans les conditions initiales en engendrent de très grandes dans les phénomènes finaux ; une petite erreur sur les premières produirait une erreur énorme sur les derniers. La prédiction devient impossible et nous avons le phénomène fortuit.

-Henri Poincaré, Science et méthode (1908).

A very small cause which escapes our notice determines a considerable effect that we cannot fail to see, and then we say that that effect is due to chance. If we knew exactly the laws of nature and the situation of the universe at the initial moment, we could predict exactly the situation of that same universe at a succeeding moment. But, even if it were the case that the natural laws had no longer any secret for us, we could still only know the initial situation approximately. If that enabled us to predict the succeeding situation with the same approximation, that is all we require, and we should say that the phenomenon had been predicted, that it is governed by laws. But it is not always so; it may happen that small differences in the initial conditions produce very great one in the final phenomena. A small error in the former will produce an enormous error in the latter. Prediction becomes impossible, and we have the fortuitous phenomena.

-Translated by Francis Maitland (1914).

Chapter 5

Space charge effects in Fixed Field Accelerators

In this chapter, we present some aspects of the space charge effects in accelerators. We review the basics, from the non-relativistic Newton equations, to the Klimontovich as well as the Vlasov equations. Implementation of dedicated modelling in the tracking code ZGOUBI and benchmarking tests are shown. Application to the case of FFAGs is the objective and object of the second part of the chapter.

Equations of motion and mean field approximation

The basic model to describe the collective interaction of a bunch of N 0 particles consists in solving the system of Newton equations in R 3 :

m i d 2 x i dt 2 = X j6 =i F j→i (t) (5.1)
where m i is the mass of the particle i, x i (t) its position, d 2 x i dt 2 its acceleration, and F j→i is the force exerted by particle j on i. This model assumes a non-relativistic motion, which, as will be shown later, corresponds to the regime in which the space charge forces play a major role (e.g. at injection). Even though this problem is highly non-trivial to solve (the N-body problem is still an open problem in mathematics), it will be helpful to have an idea about its qualitative behavior. This problem fills textbooks in mathematics, and an enormous amount of literature can be found on this subject. An excellent paper that tells the intriguing story of the attempts, from Poincaré to Sundman and Wang, to obtain an explicit solution to this problem can be found in [START_REF] Diacu | The solution of the n-body problem[END_REF]. The paradox is, although a general solution to this problem was established in 1991 by a chinese student, Quidong Wang [START_REF] Wang | The global solution of the n-body problem[END_REF], this solution in the form of series is useless from a practical point of view since it converges extremely slowly. Now, let's be more specific to the problem we are tackling: the idea is to investigate the effect of the Coulomb interaction between the particles in an accelerator. For the moment, to simplify our analysis, we will assume a free transport (i.e. no applied external forces). Thus, the basic equations governing the position of particles in a free transport channel is:

m i d 2 x i dt 2 = q i 4πǫ 0 X j6 =i q j x i -x j |x i -x j | 3
(5.

2)

The major difficulty in solving this problem is that there are as many equations as there are particles (up to 10 12 -10 14 per bunch in an accelerator). Computer simulations are thus necessary to provide an exhaustive representation of these systems. We will seek approximate solutions to this problem in our analysis and discuss its mathematical correctness.

Mean Field limit

Taking the limit N 0 → ∞ allows to replace a very large number of differential equations, by a one complicated equation. Hence we define the mean field limit [START_REF] Villani | Landau damping[END_REF]. From a theoretical point of view, this approximation is crucial: not only does it simplify the analysis but also it shows that the behavior of the system does not depend much on the exact value of the number of particles. Now, we establish the correspondence between the Newton equations and the Klimontovich equation. This is a starting point to establish the mean field limit equation or Vlasov equation. This result is particularly important for the simulation model that we develop later on in this chapter as it shows that the particles interact with the bunch as a whole rather than with the individual particles. So the average self-field generated by the bunch can be calculated to account for the space charge forces.

Klimontovich equation

In the Klimontovich scheme [START_REF] Yu | The Statistical Theory of Non-equilibrium Processes in a Plasma[END_REF], the idea is to use the microscopic single particle distribution function that is a sum of δ-functions at phase-space coordinates. Thus, we define the Klimontovich function of a system of N 0 particles by:

N k (x, p, t) = N 0 X i=1 δ [x -x i (t)] δ [p -p i (t)] (5.3)
From the definition of the function N k , it follows that:

Z N k (x, p, t) d 3 x d 3 p = N 0 (5.4)
We will start by re-expressing the Newton equations in terms of the function N k . Then, by accommodating the Maxwell equations, we will derive the Klimontovich equation [START_REF] Nicholson | Introduction to Plasma Theory[END_REF]. The Newton relativistic equations of motion are:

p i = mγ i dx i dt = mγ i ẋi (5.5) ṗi = dp i dt = q [E m (x i , t) + v i × B m (x i , t)] (5.6)
where the superscript m denotes the microscopic quantity, i.e. not averaged locally over any volume and E and B represent the applied as well as the self fields of the beam. Deriving the evolution of N k by taking the time derivative, one obtains:

∂N k ∂t (x, p, t) = - N 0 X i=1 ẋi • ∇ x (δ [x -x i (t)] δ [p -p i (t)]) - N 0 X i=1 ṗi • ∇ p (δ [x -x i (t)] δ [p -p i (t)]) (5.7)
and by accommodating Eq. (5.6), one obtains:

∂N k ∂t (x, p, t) = - N 0 X i=1 v i • ∇ x (δ [x -x i (t)] δ [p -p i (t)]) - N 0 X i=1 q h E m (x i , t) + v i × B m (x i , t) i • ∇ p (δ [x -x i (t)] δ [p -p i (t)]) (5.8)
Also, one property of the Dirac Delta function is: aδ(a -b) = bδ(a -b). This allows to replace v i with v so that the Klimontovich equation writes:

∂N k ∂t (x, p, t) + v • ∇ x N k (x, p, t) + q h E m (x, t) + v × B m (x, t) i • ∇ p N k (x, p, t) = 0 (5.9)
Although Eq.( 5.9) is exact, it is intractable. In practice however, one is interested in the average quantities rather than the information on when and where each particle is located. Thus, taking the ensemble averages of Eq.(5.9) leads to the Vlasov equation that we recall and characterize in the following section.

Vlasov equation

In order to eliminate the singularity of the Klimontovich equation, one needs to use average quantities. 1 Thus, if we smoothen over the Klimontovich function N k , one can define the distribution function f by:

f (x, p, t) = 1 Δx 3 Δp 3 Z Δx 3 d 3 x Z Δp 3 d 3 p N k (x, p, t) = hN k (x, p, t)i (5.10)
We define the quantities δf , δE and δB by:

N k (x, p, t) = f (x, p, t) + δf (x, p, t) ; hδf i = 0 E m (x, t) = E(x, t) + δE(x, t) ; hδEi = 0 B m (x, t) = B(x, t) + δB(x, t) ; hδBi = 0
Inserting these definitions in Eq.(5.9) and taking the averages, one obtains:

∂f ∂t (x, p, t) + v • ∇ x f (x, p, t) + q h E(x, t) + v × B(x, t) i • ∇ p f (x, p, t) = -q Dh δE(x, t) + v × δB(x, t) i ∇ p δf (x, p, t) E (5.11)
The Left Hand Side term (LHS) describes quantities that vary smoothly while the Right Hand Side (RHS) term represents the contribution of the rapidly fluctuating quantities. The latter contains the information about the discrete particle effects and the collisional effects.

Neglecting the collisional effects yields the Vlasov equation:

∂f ∂t (x, p, t) + v • ∇ x f (x, p, t) + q h E(x, t) + v × B(x, t) i • ∇ p f (x, p, t) = 0 (5.12)
The system is thus described by means of a smoothed-out distribution. The Vlasov equation is sometimes also referred to as the collisionless Boltzmann equation.

As explained by Vlasov in his 1950 textbook, Many-Particle Theory and Its Application to Plasma, "calculating only the successive paired interactions is by no means an expression of the reality. An essential role must be played here by forces of interaction at distances greater than the main distance between particles (we shall call them distant forces, the action of which cannot be calculated by the ordinary pattern of the kinetic equation". Formulated differently, calculating the space charge effects on a given particle by taking into account only the near particles is inaccurate. Most of the effects are coming from the bunch as a whole. 2 A detailed analysis can be found in Appendix 5.6.

Properties and qualitative behavior of the equation

Liouville's theorem

In relation to accelerator physics, the theorem can be stated as follows [START_REF] Conte | An Introduction to the Physics of Particle Accelerators[END_REF]: "In the local region of a particle, the particle density in phase space is constant, provided that the particles move in a general field consisting of magnetic fields and of fields whose forces are independent of velocity."

Invariants and identities

There are four main identities associated with the Vlasov equation [START_REF] Villani | Landau damping[END_REF]:

1. The equation preserves the total energy of the system (kinetic+potential) which is time-independent.

1 A rigorous proof of the equivalence between the Newton equations of motion and the mean field limit equation is still an open problem in mathematics.

2 A typical example of this collective behavior is the Debye screening.

2. The nonlinear integrals of the density (sometimes called casimirs of the equation) are preserved. In particular, the entropy of the system is conserved.

3. The equation is time-reversible. This means that if we transport the distribution f (x, p, t) from time 0 to time T , then reverse the momentum and let it evolve again for a period of time T , the initial distribution is recovered. This property is particularly useful for the matching of the beam optics.

4. The last property is the virial theorem which states that the average kinetic energy of the system over large times equals half of the average potential energy.

Transverse equations of motion with space charge

Recalling the transverse equations of motion (3.13 and 3.17) established in chapter 2, adding the non-linearities arising from the external magnetic field of the element (multipoles of order 2 or higher, see [START_REF] Bazzani | A Normal Form Approach To The Theory Of Nonlinear Betatronic Motion[END_REF]) and neglecting the linear coupling (no longitudinal applied field), the damping term as well as any externally applied electric fields, one obtains:

           x ′′ + 1 -n ρ 2 x = - q mγ 3 b β 2 b c 2 ∂φ ∂x + Re P M n=2 k n (s) + ij n (s) n! (x + iy) n y ′′ + n ρ 2 y = - q mγ 3 b β 2 b c 2 ∂φ ∂y -Im P M n=2 k n (s) + ij n (s) n! (x + iy) n (5.13)
where k n and j n are respectively the normal and skew gradients, defined by:

k n (s) = 1 Bρ ∂ n B y ∂x n (0,0,s) ; j n (s) = 1 Bρ ∂ n B x ∂x n (0,0,s) (5.14)
Therefore, the non-linearities may result from both self and applied forces on the beam. In order to find explicit solutions to the above equation, one needs to calculate the self fields induced by the space charge forces. For that purpose, we make the distinction between the linear and the non-linear case.

Linear space charge

For a 2D uniform density beam with elliptical cross section, x 2 /a 2 +y 2 /b 2 ≤ 1, it was shown by Kapchinskij and Vladimirskij that space charge forces become linear so that the normalized emittances are constant in a lattice where the applied forces are linear. Their distribution, generally referred to as the KV beam distribution [START_REF] Kapchinskij | Proc. Int. Conf. On High Energy Accel[END_REF] is the only known self-consistent distribution formed from Courant Snyder invariants, valid for linear periodic focusing lattices. Self-consistency here refers to the fact that the KV distribution function satisfies the time independent Vlasov equation (i.e the shape function of the beam distribution does not change with time). Besides, all 2D projections of the distribution yield a uniform particle density. Although we will not show the derivation steps that lead Kapchinskij to those results, it is important to recall the expression of the self-induced electric field generated by a uniform density elliptical beam. In free space, the field writes [START_REF] Kapchinskij | Proc. Int. Conf. On High Energy Accel[END_REF]:

- ∂φ dx = λ πǫ 0 x r x (r x + r y ) (5.15) - ∂φ dy = λ πǫ 0 y r y (r x + r y ) (5.16)
where λ is the longitudinal linear charge density, which is constant for the KV beam and r x,y are the beam edge radii for a KV beam defined by r x = 2 x 2 1/2 = 2σ x . Note that the fields are linear functions of the transverse coordinates x and y, and that the space charge introduces coupling between the two planes via the beam radii r x and r y . However, these expressions of the field are only valid within the beam. A particle outside the beam core experiences a highly non-linear space charge field scaling as 1/r and consistent with the application of Gauss' law to a uniformly charged rod (see Fig. 5.1).

x E x core halo Substituting the expressions of the field into Eqs.(5.13), and considering only the linear applied forces, one obtains:

           x ′′ + K x - 2Q perv (r x + r y )r x x = 0 ; K x = 1 -n ρ 2 y ′′ + K y - 2Q perv (r x + r y )r y y = 0 ; K y = n ρ 2
(5.17

)
where Q perv is the generalized perveance (dimensionless), defined by:

Q perv = qλ 2πǫ 0 m 0 c 2 β 2 γ 3 = q A I I 0 2 β 3 γ 3
(5.18) ≈ 6.45 × 10 -8 q A I[Amp] (γ 2 -1) 3/2 for ions of mass number A and charge number q (q/A = 1 for protons). Q perv is a measure of the strength of the space charge effects: these effects are dominant at low energies and vanish at relativistic energies (γ → ∞). One can conclude from the above that space charge forces are defocusing in both planes. It is important to mention that charge neutralization effects can be incorporated into the linear model by introducing a charge neutralization factor f e into the perveance term (as introduced by Lawson):

Q perv = I I 0 q A 2 β 3 γ 3 � 1 -γ 2 f e (5.19)
However, in our analysis, we do not discuss such effects.

Non-linear space charge

In the preceding section, beams with linear space charge were considered. This is particularly useful to simplify the analysis and yields equations that are easy to solve and interpret. However, in reality, accelerator beams do not have a uniform density. And therefore non-linear studies are important to have more complete and thorough understanding of the beam dynamics. For instance, in presence of non-linearities, the tune shift is no longer a single quantity, but is a function of the amplitude: particles of smaller amplitude have in general a larger tune shift than large amplitude particles. This phenomenon can be further enhanced in presence of the space charge defocusing effects.

To illustrate this, let's assume a 2-Dimensional Gaussian distribution where the charge density function is given by:

ρ(x, y) = λ 2πσ x σ y exp - x 2 2σ 2 x - y 2 2σ 2 y (5.20)
where σ x and σ y are the horizontal and vertical rms beam radii. The associated scalar potential in free space, which satisfies Poisson's equation can be calculated [97][98]. This yields the transverse space charge potential [99]:

φ(x, y) = λ 4πǫ 0 Z ∞ 0 -1 + exp - x 2 2σ 2 x + q - y 2 2σ 2 y + q p 2σ 2 x + q + q 2σ 2 y + q dq ≈ - λ 4πǫ 0 x 2 σ x (σ x + σ y ) + y 2 σ y (σ x + σ y ) (5.21)
Thus, in free space, the field writes:

- ∂φ ∂x = λ 2πǫ 0 x σ x (σ x + σ y ) exp - x 2 + y 2 (σ x + σ y ) 2
(5.22)

- ∂φ ∂y = λ 2πǫ 0 y σ y (σ x + σ y ) exp - x 2 + y 2 (σ x + σ y ) 2
(5.23) Expanding Eq.(5.22) for y = 0 yields,

-4 -3 -2 -1 0 1 2 3 x (cm) -4 -3 -2 -1 0 1 2 3 y (cm) Potential �-∂φ/∂x�
- ∂φ ∂x = λ 2πǫ 0 σ x x σ x + σ y - x 3 (σ x + σ y ) 3 + x 5 2!(σ x + σ y ) 5 + O(x 7 ) (5.24)
Therefore, all odd multipoles are present in the expression above. Also, note that, for beams with the same linear charge density and the same rms beam radii, the linear part is twice more important than in the case of the uniform distribution. In general, it is well understood that uniform beams yield the minimum space charge forces. Fig. 5.2 shows the gaussian space charge electric field as a function of the transverse coordinates: the core of the distribution experiences the larger electric field. The latter vanishes for x, y ≫ σ x , σ y . This results in amplitude dependent effects.

The above system is highly non-linear. The study of such systems is almost inhibited due to the fact that linearly independent solutions of the differential equations cannot be used to construct general solutions, thus the superposition principle used earlier for the linearized equations is not valid anymore. However, as explained earlier, one is not interested in the motion of each particle within the accelerator. Some particles are unavoidably lost. However, the average behavior of the beam as a whole is more relevant. For that reason, the rms quantities are used as a way to characterize the beam evolution [START_REF] Sacherer | [END_REF]. In the next section, we recall the envelope equation that describes the beam evolution in non-linear systems and discuss some of its fundamental properties widely used in beam dynamics simulations.

Remark: Assuming that the beam has a gaussian distribution in the longitudinal plane with an rms bunch length σ s , then one can calculate the linear charge density:

λ = N ppb √ 2πσ s e (5.25)
where N ppb is the number of particles per bunch.

RMS quantities

When analyzing large amount of data, it is generally needed to have a statistical approach, one in which the evolution of the bunch of particles as a whole can be characterized. For instance, statistical averages over the beam transverse phase space are defined by:

h...i ⊥ = R d 2 x ⊥ R d 2 x ′ ⊥ ...f ⊥ R d 2 x ⊥ R d 2 x ′ ⊥ f ⊥ (5.26)
If we define r x,y as the envelope edge radii for the uniform density elliptical beam and ǫ x,y as the rms-edge emittances of the beam, then, for a KV beam, it can be shown that:

r x = 2 x 2 1/2 ⊥ ; r y = 2 y 2 1/2 ⊥ (5.27)
and

ǫ x = 4 h x 2 ⊥ x ′2 ⊥ - xx ′ 2 ⊥ i ; ǫ y = 4 h y 2 ⊥ y ′2 ⊥ - yy ′ 2 ⊥ i (5.28)
The above quantities are strongly dependent on the type of distribution assumed. A more detailed analysis of the above quantities can be found in [START_REF] Lund | Generation of initial kinetic distributions for simulation of long-pulse charged particle beams with high space-charge intensity[END_REF].

Envelope equation

It can be shown that the beam envelope radii (which are s-dependent) evolve according to the KV envelope equations:

r ′′ i + k i r i - 2Q r x + r y - ǫ 2 i r 3 i = 0 ; i = x, y (5.29)
where k i is the external focusing strength and is often defined as a piecewise function of the ordinate s.

Although these equations apply to the KV beam, they are also utilized to characterize the envelope evolution of a non KV beam via the concept of RMS-equivalent beam that was introduced by Lapostolle and Sacherer in 1971: According to this concept, two beams composed of the same particle species and having the same current and kinetic energy are equivalent in an approximate sense if the second moments of the distribution are the same [89][90].

Generation of RMS-equivalent initial kinetic distributions:

In order to test the previous concept, RMS-equivalent kinetic distributions were loaded both in ZGOUBI and OPAL [START_REF] Adelmann | The OPAL (Object Oriented Parallel Accelerator Library ) Framework[END_REF] and tracked along an injection line with 5 quadrupoles. The OPAL results are shown in Fig 5 .3 and 5.4 which demonstrates its validity. There is a general hidden assumption in this concept that the beam RMS emittances remain the same. However, this is generally not true as illustrated in

Implementation and benchmarking work

ZGOUBI is a single particle tracking code that solves the non-linear equation of motion using truncated Taylor expansions of the field and its derivatives up to the 5 th order: this is performed by integrating the Lorentz equation of motion, one particle at a time. Each particle is pushed till the exit of the element before the next particle is invoked. Using the concept of non-linear transfer maps developed in [START_REF] Bazzani | A Normal Form Approach To The Theory Of Nonlinear Betatronic Motion[END_REF], the second order non-linear differential equation (4.1) can be re-written as a first order differential equation, in the standard form:

dx ds = A(s)x + f (x, y, s) (5.30)
where x ≡ (x, p x , y, p y ), A(s) is the 4 × 4 matrix that is obtained from the linearized equations of motion (see Eq. (4.3)) and f is the non-linear content of the equation, defined by: Also, we note M (l) the transfer map of the element (l) where l is the length of the element. An elegant solution of the above equation is derived in [START_REF] Bazzani | A Normal Form Approach To The Theory Of Nonlinear Betatronic Motion[END_REF]. In what follows, we show the main steps of the latter in order to explain the cornerstone of the space charge calculation module that we implemented in the tracking code ZGOUBI, and to estimate the errors in view of validating the convergence of the numerical scheme: Let L(s, s l-1 ) be the fundamental matrix for the above linear system (i.e f = 0). Thus, L satisfies the equation:

A =     0 1 0 0 (n -1)/ρ 2 0 0 0 0 0 0 1 0 0 -n/ρ 2 0     f (x, y, s) =     0 f x (x, y, s) 0 f y (x, y, s)    
dL(s, s l-1 ) ds = A(s)L(s, s l-1 ) (5.31)
and the non-linear equation has a solution which writes in the integral form:

x(s) = L(s, s l-1 )x(s l-1 ) +

Z s s l-1 L(s, s ′ )f (x, y, s ′ )ds ′ = L(s, s l-1 ) x(s l-1 ) + Z s s l-1 L(s l-1 , s ′ )f (x, y, s ′ )ds ′ ! (5.32)
where x(s l-1 ) is the initial condition.

For the case of bunched beams, the longitudinal bunch length is calculated and the linear charge density updated by means of Eq. (5.25).

Thin lens approximation

In what follows, we use a method that has been developed in [START_REF] Bazzani | A Normal Form Approach To The Theory Of Nonlinear Betatronic Motion[END_REF] to construct approximate analytical solutions of Eq.(5.32): we consider that the nonlinearity f is concentrated in one place of the interval [s l-1 , s l ]. This method has the merit of providing approximate maps which are symplectic, an essential requirement for tracking codes. The idea is to solve the integral equation above by assuming that the nonlinear content f acts only at the beginning of the element. Thus, the following substitution is made:

f (x, y, s) → f (x, y, s)lδ (s -(s l-1 + 0)) (5.33)
where δ is the Dirac delta function.

The kick has no effect on the particle coordinate x (where applied). However, one has a discontinuity in the momentum:

Z s l-1 +ǫ s l-1 dp x ds ds = p x (s l-1 + ǫ) -p x (s l-1 ) = lf x (x, y, s l-1 ) + O(ǫ) (5.34)
In the limit where ǫ → 0, the kick map is given by:

    x p x y p y     s=s l-1 +0 =     x p x + lf x (x, y, s) y p y + lf y (x, y, s)     s=s l-1
Finally, the one-kick approximation will be given by:

x one-kick (s l ) = L(s l , s l-1 ) (x(s l-1 ) + lf (x, y, s)) (5.35) or formulated in a matrix form for an FFAG element and using the notation of the previous chapters:

    x x ′ y y ′     s l =     0 1 0 0 (n -1)/ρ 2 0 0 0 0 0 0 1 0 0 -n/ρ 2 0         x x ′ + lf x (x, y, s) y y ′ + lf y (x, y, s)     s l-1

Error estimate

It is possible to estimate the error e resulting from the thin lens approximation by using the general rules of integration (see [START_REF] Bazzani | A Normal Form Approach To The Theory Of Nonlinear Betatronic Motion[END_REF] for details). The approximate solution can read:

x(s l ) = x one-kick (s l ) + e ; e = O(l 2 )

(5.36)

Thus, the convergence is ensured if l → 0. And one can readily reduce the error of the method by splitting the magnetic element into smaller slices. If the element is split into m equal slices, then one still obtains a symplectic map and the error of the method is estimated to be O(l 2 /m). The convergence is thus guaranteed when m → ∞. The above result is only an estimate and shall depend on whether the orbits remain bounded or not. However, in the next section, we will show the dependence of the error of our numerical simulation on the number of kicks.

Implementation of the method

The idea of the implementation is to calculate the space charge kick. Once determined, the particles are advanced through the lattice by multiplying the transfer maps of the applied forces on the beam with those generated by the internal forces, i.e. the space charge kick. The ZGOUBI code solves the non-linear equation of motion without the internal forces of the beam, i.e. the left term on the RHS of Eqs.(5.13) is dropped. Given that the transfer map is obtained by the composition of the element maps, the thin lens approximation is invoked by incorporating the space charge terms. The latter is determined when all particles in the distribution reach the exit of the element by:

Δx ′ sc = - q mγ 3 b β 2 b c 2 ∂φ ∂x Δs = 2πǫ 0 Q perv λ - ∂φ ∂x Δs (5.37) Δy ′ sc = - q mγ 3 b β 2 b c 2 ∂φ ∂y Δs = 2πǫ 0 Q perv λ - ∂φ ∂y Δs (5.38)
where Δs is the length of the element. Thus, the problem reduces to calculating the transverse electric fields of the distribution at the entrance of each slice,

E x E y s=s l-1 =      - ∂φ ∂x - ∂φ ∂y      s=s l-1
(5.39)

and the space charge kick map for each element is given by:

    x x ′ + Δx ′ sc y y ′ + Δy ′ sc     s l-1 (5.40)
A frozen space charge model is employed for our numerical simulations, where one assumes that the analytic expression of the particle distribution remains the same, so that the analytical solution of the self-induced electric fields does not change. For instance, if we employ a KV beam model, then the problem reduces to calculating the edge radii r x and r y of the beam elliptical cross section. Besides, if we cut the magnet into m thin slices, one may assume that the beam radii do not change much within each slice and the error scales as 1/m. However, as one may expect, the error relates to how fast the change of r x and r y occurs, i.e. to the frequencies of the betatron oscillations.

Convergence of the calculation:

For better convergence, the step size has to be chosen carefully. It is generally found that a step size less than 1 cm ensures good convergence of the numerical calculation: given that the kick applied to each particle is concentrated at the entrance of each slice, and that the non-linearities inherent to the lattice may start to act on the beam, it is important to integrate the equation of motion as soon as the kick is applied, therefore with a fine step size. Fig 5.6 shows the tracking results in a drift where one varied the number of slices utilized for the calculation of the space charge effects: as expected, increasing the number of slices helps reducing the error (see Fig. 5.7). The latter scales as 1/m. A relaxation regime is observed where the error does not change much. This is expected since the space charge becomes weaker with increasing beam radii. Note that the element length must remain larger than the integration step size in order to avoid any issues. In order to test the results, one solved the non-linear differential equation in a drift where one assumed a symmetric KV beam distribution such as r x = r y = r b . The envelope equation writes: Note that the agreement is quite good for all cases. However, the results are quite sensitive to the difference in the initial conditions which arises from the noise of the generated kinetic distribution. Besides, the beam emittance remains constant in accordance with the KV beam model.

r ′′ b - Q r b - ǫ 2
Q/Q 0 =0 case1 Q/Q 0 =1 case1 Q/Q 0 =2 case2 Q/Q 0 =0 case2 Q/Q 0 =1

Space charge effects in FFAGs

When carrying out parametric studies of space charge effects, it is found that the depressed phase advance is the natural parametrization, simply because the phase advances are the most relevant parameters for analysis of resonance-like effects. The undepressed phase advance per lattice period provides a measure of the strength of the applied focusing describing the lattice. Namely, for FFAGs, we have shown that the undepressed betatron tune can be expressed by:

ν 2 x (R E ) = X i β i (R E ) - X i α i (R E ) × k i (R E ) (5.42) ν 2 y (R E ) = X i α i (R E ) × k i (R E ) + F 2 1 + 2 tan 2 (ξ) (5.43)
This shows that in order to vary the undepressed phase advance, one can vary the scaling factors of the F and D magnets (also the fringe field). These parameters are particularly useful for non-linear magnets where the local field index changes with the radius as well as the azimuthal angle of the particle orbit. As shown in Fig. 5.9, the phase space is no longer elliptic for large amplitude particles: the presence of space charge effects may therefore augment these non-linearities which are due to the presence of high order multipoles in the magnets. Therefore, the study of such systems requires appropriate tools which account for the non-linearities of the applied as well as the internal forces of the beam. For that reason, the tracking code ZGOUBI is employed with the newly implemented space charge module. Based on this, we will investigate the lattice sta- .9: Horizontal phase space at injection for the KURRI 150 MeV FFAG: as can be observed, when the initial amplitude increases, the motion invariant is deformed (here, by proximity of a 3rd integer tune). The limit of stable motion is obtained when the particle is lost in less than 10000 turns. This represents that largest phase-space area that will pas through the FFAG lattice, also referred to as the horizontal aperture (black curve).

bility from the point of view of space charge. Some of the results that we present here can also be found in [START_REF] Haj Tahar | Space Charge Effects In FFAG[END_REF].

Procedure:

1. Generate a median plane field map for a given (k F , k D ) (see Eq. (4.50)).

2. Search for NCO closed orbits between injection and extraction using the built-in fitting routines in ZGOUBI.

3. For each closed orbit, the matching condition ensuring a periodic motion around the closed orbit is obtained and the matched distribution is generated and tracked over several turns (120 turns). One assumes a symmetric beam distribution, i.e. ǫ x = ǫ y = ǫ. The damping of beam emittance with momentum is taken into account.

4. For each closed orbit, the number of betatron oscillations of the distribution is computed using the ZGOUBI Discrete Fourier Transform (DFT):

ν co x,y = hν x,y i = 1 N par N par X i=1 ν x,y,i (5.44) σ co νx,y = ν 2 x,y 1/2 = 1 N par N par X i=1 � ν x,y,i -ν co x,y 2 ! 1/2
(5.45)

From that, one then determine the average tune of the lattice ν m x,y as well as the rms tune over the closed orbits ν rms x,y . 5. The same steps are repeated in presence and in absence of the space charge effects: given that the emittance growth time is much shorter than the synchrotron period, our numerical simulation consists of a frozen longitudinal phase space. This provides the tune depression. We assume that the FFAG is operating in emittance-dominated regime. Therefore, the space charge is treated as a small perturbation.

Simple tests of convergence include varying the step size, the number of particles as well as the number of turns.

KV beam distribution

For a fast parametric study, a frozen space charge model is employed. One has to keep in mind though that the KV model does not include any non-linear effects that increase the emittance. This test is particularly useful to validate the built model. In order to compare our results with the analytical model, the initial emittances are chosen to ensure that the undepressed particles trace elliptical orbits. As described in [START_REF] Tahar | Space charge studies in FFAG using the tracking code ZGOUBI[END_REF], each particle in the distribution experiences a linear space charge kick given by:

Δx ′ = 2Q perv (r x + r y )r x xΔs (5.46
)

Δy ′ = 2Q perv (r x + r y )r y yΔs (5.47)
In order to calculate the space charge tune shift, we recall the Laslett tune shift formula applied to a KV beam:

Δν y = 1 4π I C 0 β y (s) 2Q perv r y (r x + r y ) ds (5.48)
After simplification, and recalling the hypothesis ǫ x = ǫ y = ǫ, this yields :

Δν y ν y 0 ≈ 1 2π Q perv ǫ 2 r 3 y r x + r y ≈ RQ perv ǫ ν -3/2 y 0 ν -1/2 x 0 + ν -1/2 y 0 (5.49)
where ν y 0 is the undepressed vertical betatron tune and R is the average radius of the particle orbit in the accelerator. Interchanging x and y gives the horiontal tune shift.

Introducing the normalized emittance ǫ n = βγǫ yields the scaling Δν x,y ∝ R/(βγ 2 ).

Scaling FFAGs

Assuming that the average field index of the focusing and defocusing magnets are the same, thus the undepressed tune is constant, one can calculate the effects of space charge forces on the transverse beam dynamics of scaling FFAG: the procedure described above is applied to a scaling FFAG case by varying the average field index k of the magnets as well as the linear charge density λ of the beam, thus the perveance Q perv . The results of tracking using the space charge module implemented in ZGOUBI are shown in Fig. 5.10. In order to interpret the results, we inject the approximated tune formula: ν 2 x 0 ≈ k + 1 and ν 2 y 0 ≈ -k + F 2 into Eq. (5.49). This yields: Given that the magnetic flutter as well as the radius of the closed orbits have a weak dependence on the scaling factor, it results that:

Δν x ν x 0 ≈ RQ perv ǫ (k + 1) -3/4 (k + 1) -1/4 + (-k + F 2 ) -1/4 ≈ f x (k, Q perv , ǫ) (5.50) Δν y ν y 0 ≈ RQ perv ǫ (-k + F 2 ) -3/4 (k + 1) -1/4 + (-k + F 2 ) -1/4 ≈ f y (k, Q perv , ǫ) (5
� �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � � � � � � � � � � � � � � � � � � � � � � � � ��
� ���� � ��� �� ��� �� � �� � � �� �� � � �� � � �� �� � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � � �� � � � � �� � � � � ��
� �� � � � �� � � � �� � � � � �� � � � � �� � � � � � � � � � � �� � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � � ��
k ր =⇒ ν x ր and ν y ց =⇒ Δν y ν y ր and Δν x ν x ց (5.52)
and vice versa. Besides, the space charge tune shift has a quasi linear dependence on the perveance term: increasing the linear charge density of the KV beam increases the tune shift. Also, increasing the beam emittance does reduce the tune shift as the physical beam size increases and thus the space charge forces reduce. A contour plot of the previous formula (Fig. 5.11) shows good agreement with the tracking results: increasing the k value increases the horizontal focusing and decreases the vertical one. Thus the tune shift exhibits the opposite behavior. Besides, the tune shift is quasi linear with respect to the perveance term, which is a good indication that the beam matching is well ensured. Some differences are observed between the analytical formula and the tracking results, mainly in the vertical plane: this is due to our assumption that the magnetic flutter F is independent of the k-value of the lattice, which is only an approximation. Now, we set the limit of the maximum beam current that can be transported through the FFAG lattice as the limit at which the tune is fully depressed so that the net focusing force is zero (conventional criterion). Given that the vertical focusing is in most cases weaker than the horizontal one in fixed field accelerators, this limit (for symmetric beams!) is reached when Δν y = ν y0 or f y (k, Q, ǫ) = 1. However, the effect of the dispersion due to the momentum spread in the particle distribution should also be taken into account and included in the calculation of the tune shift due to space charge. In order to allow this, we re-write the off-momentum linearized equation of motion in the vertical plane [START_REF] Ruggiero | Design Criteria of a Proton FFAG Accelerator[END_REF]:

y ′′ + n ρ 2 1 1 + δ y = 0 ; n ρ 2 ≈ k R 2 µ (5.53)
which transforms into:

y ′′ + k 1 + δ µ R 2 y = 0 (5.54)
Therefore, including the chromatic effects in the vertical plane, is equivalent to replacing the average field index k of the magnets by k/(1 + δ) (One has to pay attention though to the change of the off-momentum particle orbit and therefore to the change in the magnetic flutter F). Finally, the beam current limit is reached when:

f y (k, Q perv , ǫ, δ) = RQ perv ǫ - k 1 + δ + F 2 -3/4 k 1 + δ + 1 -1/4 + - k 1 + δ + F 2 -1/4 = 1 (5.55)
Recalling the expression of the perveance term for protons, this yields:

I max [A] = C ǫ R � γ 2 -1 3/2 k 1 + δ + 1 -1/4 + - k 1 + δ + F 2 -1/4 - k 1 + δ + F 2 -3/4 ; C = 1.55038 × 10 7 A (5.56)
As expected, the maximum beam current that can be transported through the FFAG lattice increases with increasing beam injection energy and is proportional to the geometrical beam emittance at injection. However, I max decreases when increasing the average field index k of the magnets as illustrated in Fig. 5.12. This is explained by the fact that the vertical focusing becomes weaker while the horizontal one increases when increasing the k value. The above analysis is valid for a continuous or unbunched beam. However, this can be extended to the case of bunched beams by replacing the beam current I by the peak current in the bunch Î. For that, as shown in [START_REF] Reiser | Theory and Design of Charged Particle Beams[END_REF], it is customary to define the bunching factor B f , the ratio of the average current to the peak current:

B f = I Î (5.57)
Thus, B f has the range 0 < B f � 1 and B f = 1 represents the unbunched beam treated above. The most noticeable space charge effect on the beam dynamics of a scaling FFAG is that the property of having a fixed tune machine is no longer valid. This is explained by the fact that the tune shift of a KV beam scales as Δν x,y ∝ RQ/ǫ ∝ R/(βγ 2 ). To show this, one calculates the ratio R σν of the tune variations of a scaling FFAG with and without space charge and normalized by the mean value of the tunes, which is a measure of the strength of the perturbation: lattice that maintains the fixed tune property in presence of space charge. For that reason, the non-scaling FFAG case is treated next: this allows more degrees of freedom to play with, such as the κ values of the DFD triplet which enable the tune variations in the cell.

Non scaling FFAGs

In the general case of a non-scaling FFAG, the similarity of the orbits is no longer valid and therefore the situation is more complicated than in the scaling case: the tunes become energy dependent and the magnetic flutter as well. The procedure described above is applied again: the perveance term is now fixed and a scan on (k F , k D ) allows to obtain the tune shift of a space charge perturbed beam. The results are shown in Figs. 5.14 and 5.15 and are in good agreement with Eq.(5.49). In an average sense, increasing the horizontal tune while maintaining the vertical one constant decreases the horizontal tune shift although the vertical one seems quite insensitive. Now, in an effort to answer the previous question, i.e. to find an FFAG lattice that does not suffer the tune shift induced by the space charge forces, we define the following parameter:

A ν = ln � R x σν R y σν (5.59)
A ν is a measure of the change of the tune excursion in presence of the space charge forces: if A ν < 0, then adding the space charge effects helps reduce the overall tune variation. As can be seen in Fig. 5.16, this condition can be satisfied if κ < 0. This result is expected: as shown earlier, if κ < 0, the bare tune is a decreasing function of the energy in both planes. Therefore it counteracts the space charge tune shift which is an increasing function of the energy in both planes. To be more quantitative, the problem writes in the following way: where all quantities are energy-dependent. In order to determine (δk F , δk D ) that needs to be applied in order to counteract the space charge effects, one equates the RMS tune excursion of the KV beam, ν rms x,y (ǫ, Q perv ), with the RMS tune excursion of the bare tunes. Recalling the expression of the latter (established in chapter 3) and solving for κ, this yields:

( ν x = ν x 0 + δν x (δk F , δk D ) + δν SC x (ǫ, Q perv ) ν y = ν y 0 + δν y (δk F , δk D ) + δν SC y (ǫ, Q perv ) ��� ��� ��� ��� ��� ��� ��� ��� ��� ��� �� � � � � � �� �� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � ����� � �� � � �� �� � � �� � � ��� �� �� �� � � �� � � ��� �� � � � � � �� �� � ��� ��� ��� ��� ��� ��� ��� ��� ��� ∆� � �� � ���� (a) (kF , kD) map. ��� ��� ��� ��� ��� ��� ����� ��� ��� ��� ��� ��� ��� ����� ����� � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � ��
� �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � �� � � � � � � � � �� � � � � � � � � �� � � �� �� � � �� � � �� �� � � �� � � �� �� � � �� � � �� �� � � �� � � �� �� � � � �� � � �� �� �� � � �� � � �� �� �� � � �� � � �� �� �� �� �� �� � � �� � � �� �� �� �� �� �� � � �� � � � � �� � � �� �� �� � � �� � � � � �� � � �� ��� � �� �� �� � � �� � � � � ��
|κ| ∈ ν rms x (ǫ, Q perv ) a x : ν rms y (ǫ, Q perv ) a y (5.60)
In order to test the previous scheme, we choose a lattice corresponding to k F = k D = 7.6. If we choose to maintain k F fixed, then Eq. (5.60) yields

k D ∈ k F - ν rms x (sc) a x : ν rms y (sc) a y ∈ [7.0 : 7.11] (5.61)
The results are shown in Fig. 5.17 which proves the validity of this method. The undepressed phase advance of the lattice is sensitive to the average field index of the magnets, which in turn provides a useful way to control the tune path in an FFAG. As explained earlier, this is mainly due to the non-scaling of the orbits which experience a continuously changing magnetic flutter. In what preceded, this effect was exploited in such a way as to counteract the continuous decrease of the space charge forces which is due to the cancellation of the electrostatic defocusing effect by the self-induced magnetic field. However, this result assumed that all forces were linear as guaranteed and verified by the uniform beam distribution. In reality, the RMS emittance can go up and down due to non-linearities. The tune spread does increase accordingly. One might think of different ways to correct the beam non-linearities: plugging multipoles for instance can help reduce it (although it rapidly reduces the dynamic aperture). However the non-linearities due to space charge forces are wave-like functions. From that point of view it might be difficult to correct (the correction is a point-like action while the space charge forces act continuously on the beam). This means that when aiming to remediate the problem of space charge effects, this should be considered in the preliminary design phase of the accelerator. For that reason, analysis of the space charge effects in non-linear FFAGs is crucial in order to assess their potential for high power applications.

Gaussian beam distribution

The previous study assumed a KV beam distribution. Therefore, the emittances were unchanged. However, for a Gaussian beam distribution, the situation is more complicated: each particle in the distribution experiences a non-linear space charge kick given by:

Δx ′ = Q perv (σ x + σ y )σ x exp - x 2 + y 2 (σ x + σ y ) 2 xΔs (5.62) Δy ′ = Q perv (σ x + σ y )σ y exp - x 2 + y 2 (σ x + σ y ) 2 yΔs (5.63)
The emittance increases with increasing turn number and with increasing beam current as shown in Fig. 5.18: a matched gaussian distribution is injected at 11 MeV and tracked over 100 turns for a scaling FFAG where k = 7.6. By varying the beam current or equivalently the number of particles per bunch, one can observe that the RMS emittance increases considerably. Such an increase is due to two contributions: the non-linear forces arising from the space charge fields as well as those of the externally applied fields of the FFAG. For this reason, computer simulations taking into account the non-linear terms of the equation of motion are crucial in order to study the beam dynamics of FFAGs. The tune spread is one important illustration of this: while the betatron tune shift is the largest for the core particles, the large amplitude particles experience a near zero tune shift (see Fig. 5.19). In order to better understand the space charge effects on the beam dynamics of a non-scaling FFAG, we consider the case where the lattice cells have κ > 0. As shown in Fig. 5.20, although the tune spread of the bunch reduces with the increasing beam energy, the tune excursion of the core remains quite similar to that of the bare tunes. Therefore, the previous approach of perturbing the field index of the magnets does not remediate the tune spread of the beam (which is absent in case of linear forces). However, such an approach remains valid for the core of the distribution which can be positioned in the tune diagram in a way to avoid overlapping of harmful resonances.

Another interesting investigation is to compare the tune spread of a scaling FFAG with that of an equivalent fixed tune alternating κ lattice. Tracking simulations in ZGOUBI based on field maps revealed that both lattices behave in the same way vis-à-vis space charge effects (see Fig.

Conclusion

In this chapter, simulations of the space charge effects were undertaken: to begin with, one focused on the mathematical correctness of the models by analysing the correspondence between the Newton equations and the Vlasov equations. After reviewing the basics of the space charge (or collective) regime, the implementation work based on the thin lens approximation was shown. Although the benchmarking work is still on-going to validate the simulation results, the results obtained so far are in good agreement with the analytical models. Application to the case of FFAGs where some general scaling laws were derived was also shown. In this analysis one made the distinction between the scaling and the non-scaling FFAG case. Although the fixed tune property of the scaling FFAG is no longer valid in presence of the space charge effects, it is found that the tune excursion can be controlled and even cancelled by introducing scaling imperfections, hence the study of the non-scaling FFAG. Study of the alternating-κ fixed tune FFAG and comparison with the equivalent scaling FFAG showed that the tune spread in both cases is very similar. In the future, it will be interesting to investigate the effect of crossing the resonance lines in both cases: given that the alternating-κ FFAG doubles the superperiod, an effect of which is to double the resonance population, one might think that the impact of the resonance crossing is more severe than the scaling FFAG case. This needs to be investigated, particularly in presence of space charge effects.

Conclusion

The nuclear waste is still an open problem in the nuclear industry. There is no intrinsically safe way to handle the spent fuel of the existing Light Water Reactors. The most viable options involve transmutation of the most radiotoxic elements in advanced reactor concepts and storage in underground repositories. Although there is no accurate estimate of the cost of deploying each technology, analysis of the concept of ADSR that we present in chapter 1 and further analyse in chapter 2, showed an increase in the levelized cost of energy by at least 20 %. Thus, besides the high risk of deploying an immature technology, ADSR cannot compete on a cost level with conventional technologies. The analysis that followed focused on defining the beam requirements of the accelerator required as well as compare the existing technologies: essentially, cyclotrons, linacs and FFAGs were considered here. Although cyclotrons and linacs can compete to produce a reliable high power beam to the target, the costs as well as the need for a multi-accelerator complex are in favour of the cyclotron technology. The FFAG, although a potential candidate for this application, has yet to demonstrate the high beam power capability. This requires several upgrades of the existing facilities

In KURRI, a demonstrator of the ADSR technology is built where the accelerator was based on the scaling FFAG technology. As shown in chapter 3, a simulation campaign established in 2014 to benchmark the different codes yielded excellent agreement. Furthermore, analysis of the field profile revealed some major sources of imperfections related to the radial and azimuthal variations of the average field index of the magnets.

In order to understand the effects of the latter on the beam dynamics, both analytical and simulation models were developed in chapter 4, which demonstrated that the tune variations obey a well defined law that involves the non-scaling of the orbits as well as the average field index of the magnet. As a result of its monotonic behavior with κ, the difference of the average field index of the focusing and defocusing magnets (property established in chapter 4), one was able to explain and quantify the impact of field defects on the transverse beam dynamics of the KURRI FFAG. As a remedy to these field defects, a correction scheme is proposed which is based on the idea of alternating κ. It is shown that the variation of the betatron oscillations with the energy is reduced by at least one order of magnitude when this scheme is employed. This is the cornerstone of the novel FFAG concept that we present and discuss in the last part of chapter 4. A key finding of this study is that the cardinal conditions of a scaling FFAG are non-necessary in order to obtain a fixed tune FFAG. In summary, this novel FFAG concept guarantees a strong focusing by alternating the gradient, but also ensures that the focusing is constant by alternating the azimuthal differences of the average field index of the magnets. Further analysis of this concept revealed that the dynamic acceptance is lower with increasing κ-value. This means that such values have to be well controlled in order to maintain a large beam acceptance necessary for high power applications.

Chapter 5 is dedicated for the latter: in order to undertake the beam dynamics simulations of high intensity beams, one used the thin lens approximation to implement the space charge kick in the tracking code ZGOUBI. Several tests were carried out to validate the simulation results. Although this still requires lots of dedicated work, the obtained results are in good agreement with the analytical models. Application to the FFAG was the object of the second part of this chapter where one showed that the tune excursion can be controlled and even cancelled. In the future, it will be important to investigate the effects of the resonance crossings on the beam dynamics: of particular interest is the comparison of the two cases: the scaling FFAG case and the equivalent alternating κ fixed tune FFAG.

Appendix A: Derivation of the self forces As shown in [START_REF] Lund | Transverse equilibrium distribution[END_REF] [START_REF] Barnard | Introduction[END_REF], the Lorentz force can be decomposed into applied and self-field forces where:

     F em = F a + F self = q(E + v × B) E = E a + E self B = B a + B self (64)
The self-field Maxwell equations, in steady state, are then invoked:

   E self = -∇φ B self = ∇ × A ; A = s β b c φ (65)
The main assumption is that the beam evolution is sufficiently slow to approximate the self interactions with electrostatic interactions in the frame moving with the beam. Inserting these expressions into F self yields:

F self = qE self + qv × B self = q - ∂φ ∂x + qβ b cs × ∂ ∂x ⊥ × s β b c φ = q - ∂φ ∂x + qβ 2 b ∂φ ∂x ⊥ (66)
where ∂φ ∂x ⊥ = ∂φ ∂x x + ∂φ ∂y y.

Finally, inserting ∂φ ∂x = ∂φ ∂x ⊥ + ∂φ ∂s into (66) yields:

F self = -q 1 γ 2 b ∂φ ∂x ⊥ -qs ∂φ ∂s (67)
where ξ(r) = ξ 0 + ξ 1 × r + ξ 2 × r 2 + ξ 3 × r 3 , θ is the azimuthal angle (origin θ=0) and ω is a parameter used to position the EFB with respect to the azimuthal position θ=0.

According to this model, the magnet gap is also allowed to vary: g is given as a function of the radius by: g(r) = g 0 + g 1 × r + g 2 × r 2 (71)

The field is then extrapolated off median plane by means of Taylor series: for that, the median plane antisymmetry is assumed and the Maxwell equations are accommodated. Note that the "CYCLOTRON" element allows the overlapping of 5 such dipole fields. This follows the method described in [START_REF] Lemuet | Developements in the ray-tracing code Zgoubi for 6-D multiturn tracking in FFAG rings[END_REF].

In the case of a cyclotron machine, the isochronicity is a crucial point: Because the revolution time has to be constant (f rev = qB 2πγm 0

), this implies that the radial dependence of the field must be proportional to γ, so that R(R) ∝ γ(R), where R is the average radius of the orbit. Since

f rev = v C
, where C is the path length of the particle for one closed orbit, one obtains, with a good approximation, that R ∝ β. Thus,

R(R) ≈ 1 r 1 -( R R 0 ) 2 (72)
A model of a 6-sector cyclotron is built where one varied the spiral angle as illustrated in Au cours des dernières années, le concept de réacteur sous-critique piloté par accélérateur de particules (ADSR) a gagné plus d'intérêt dans le monde entier en tant que candidat potentiel pour résoudre le problème des déchets nucléaires. Ce concept, objet de la présente dissertation, nécessite le couplage d'un accélérateur à protons de forte puissance avec un réacteur sous-critique. Parmi les nombreuses options qui peuvent produire le faisceau de forte puissance, on note les accélérateurs à champ fixe et gradient alterné (FFAG). Le travail présenté dans cette thèse est dédié à démontrer la capacité de haute puissance des FFAG. Dans notre étude, nous nous concentrerons sur une installation de recherche ADSR: l'Assemblée critique de l'Université de Kyoto (KUCA), une évaluation de faisabilité de l'ADSR.

F.1.1 ADSR à KURRI

À l'institut de réacteurs de recherche de l'université de Kyoto (KURRI), la première injection au monde de neutrons de spallation générés par des protons à haute énergie dans un coeur de réacteur a été réalisée le 4 mars 2009. Le KUCA est utilisé comme noyau sous-critique pour cette expérience. Un synchrotron FFAG accélère le faisceau de protons jusqu'à 100-150 MeV qui est délivré à une cible de tungstène (W).

Un FFAG est un accélérateur circulaire à champ fixe avec une forte focalisation, l'alternance du gradient du champ magnétique permettant la forte focalisation dans les deux plans.

F.2 Couplage accélérateur-réacteur

En guise de conclusion de ce chapitre, le principal inconvénient de la technologie ADSR est le coût élevé de l'énergie (LCOE) par rapport à d'autres concepts de réacteurs avancés qui n'emploient pas d'accélérateur. De nos jours, il s'agit d'un bouclier pour toute application industrielle visant à produire de l'énergie (sans traiter le problème des déchets). En outre, le réacteur n'est pas intrinsèquement plus sûr que les concepts critiques du réacteur, compte tenu de la complexité de la gestion de l'interface cible entre l'accélérateur et le noyau du réacteur (les transitoires de surtension sont un exemple). Cependant, plusieurs des principaux défis (haute puissance, sécurité, fiabilité) de l'installation peuvent être satisfaits en utilisant une configuration de cible multi-faisceaux où la puissance du faisceau est répartie entre plusieurs accélérateurs. Cela réduit les densités de puissance sur la cible, améliore les caractéristiques de transmutation qui devient plus uniforme et assure une meilleure sécurité de l'installation. Dans cette configuration, les cyclotrons ont le potentiel d'atteindre les exigences, car cette technologie fournit un coût moins cher que tout linac supraconducteur.

F.3 Étude de la dynamique des faisceaux dans les accélérateurs circulaires à champ fixe

Dans ce chapitre, un travail d'analyse comparative a été effectué permettant de valider les outils de simulation et d'établir le niveau de précision requis dans la simulation pour une bonne convergence des résultats calculés. On trouve en particulier que la solution approximative des équations linéaires du mouvement est en accord avec celle des équations non linéaires basées sur les développements en séries de Taylor du champ et ses dérivées jusqu'au 5ème ordre. En outre, l'analyse des variations du nombre d'onde montre que ces derniers ne sont pas négligeables, de sorte que la formule de Symon pour les morceaux ne s'applique pas. Ainsi, la mauvaise transmission globale peut résulter de la traversée des résonances Betatron. Afin d'expliquer l'origine de ces défauts, on a entrepris une analyse de l'indice du champ qui a révélé des imperfections radiales et azimutales du champ. Certaines définitions ont été introduites pour tenir compte de ces imperfections et une analyse détaillée de la stabilité du faisceau est présentée dans le chapitre suivant.

F.4

Étude de la stabilité du faisceau et concept nouveau de FFAG Dans ce chapitre, on a analysé l'effet des erreurs de champ sur la stabilité des trajectoires de particules. Plusieurs approches du problème ont été développées et discutées. La comparaison des résultats a montré que l'approche linéaire n'est valable qu'au voisinage des orbites fermées: tous les paramètres du modèle de bord évoluent en raison de la non-mise à l'échelle des orbites. Par conséquent, on s'est basé sur l'approche non linéaire basée sur des simulations de suivi. Un résultat crucial a été d'établir une relation entre le nombre d'oscillations betatron et les défauts de champs. Un paramètre clé pour mesurer l'amplitude des défauts est la valeur κ définie comme la différence de l'indice de champ moyen des aimants focalisats et défocalisats. Sur la base de ces résultats, un nouveau schéma pour corriger la variation des oscillations betatron avec l'énergie a été proposé. L'idée principale consiste À alterner les valeurs κ des aimants, tous les deux (ou plus) secteurs. Cela a conduit à un nouveau concept de FFAG à focalisation invariante sans mise à l'échelle que l'on a développé dans la dernière partie de ce chapitre: en plus du fait que cela démontre que les conditions de mise à l'échelle ne sont pas nécessaires pour obtenir un FFAG à focalisation invariante, le concept est plus facile à mettre en oeuvre au moyen de bobines de correction qui peuvent être ajustées pour minimiser la variation du nombre d'onde et éviter la traversée de résonances. L'analyse de l'acceptation dynamique a montré que celle-ci est plus faible pour le FFAG avec des valeurs κ en alternance. Cela a montré que l'amplitude de κ doit être bien contrôlée afin de maintenir une grande acceptation du faisceau. Cependant, cela nécessite une enquête plus approfondie.

F.5 Effets de la charge d'espace dans les accélérateurs à champ fixe Dans ce chapitre, des simulations des effets de la charge d'espace ont été entreprises: d'abord, l'accent sur la validité mathématique des modèles en analysant la correspondance entre les équations de Newton et les équations de Vlasov. Après avoir examiné les bases du concept de charge d'espace, le travail d'implémentation basé sur l'approximation de la lentille mince a été montré. Bien que le travail de benchmarking soit encore en cours pour valider les résultats de la simulation, les résultats obtenus jusqu'à présent sont en accord avec les modèles analytiques. L'application au cas des FFAG où certaines lois générales d'échelle ont été dérivées a également été démontrée. Dans cette analyse, on a fait la distinction entre la mise à l'échelle et la non-mise à l'échelle de FFAG. Bien que la proprià c tà c de focalisation invariante du FFAG ne soit plus valide en présence des effets de charge de l'espace, on constate que l'excursion du nombre d'omde peut être contrôlée et même supprimée en introduisant des imperfections de mise à l'échelle, d'oû l'étude du FFAG avec κ alternatif. À l'avenir, il sera intéressant d'étudier l'effet de la traversé des lignes de résonance dans les deux cas: étant donné que le FFAG à κ-alternatif double la population de rÉsonance, on pourrait penser que L'impact de la traversée des résonances est plus sévère que celui des FFAG où les indices de champ magnétique dans les aimants F et D est identique. Cela doit être étudié, en particulier en présence d'effets de charge d'espace.

F.6 Conclusion

Les déchets nucléaires constituent toujours un problème ouvert dans l'industrie nucléaire. Il n'existe pas de moyen intrinséquement sûr de manipuler le combustible usé des réacteurs d'eau légére existants. Les options les plus viables impliquent la transmutation des éléments les plus radiotoxiques dans les concepts de réacteurs avancés et le stockage dans les dépôts souterrains. Bien qu'il n'y ait pas d'estimation précise du coût du déploiement de chaque technologie, l'analyse du concept d'ADSR que nous présentons au chapitre 1 et qui est présenté plus en détail au chapitre 2, a montré une augmentation du coût d'énergie nivelé par au moins 20 %. Ainsi, outre le risque élevé de déploiement d'une technologie immature, ADSR ne peut pas concurrencer au niveau des coûts avec les technologies conventionnelles. L'analyse qui a suivi a porté sur la définition des exigences de faisceau de l'accélérateur requises ainsi que sur la comparaison des technologies existantes: essentiellement, cyclotrons, linacs et FFAG ont été considérés ici. Bien que les cyclotrons et les linacs puissent concurrencer pour produire un faisceau de haute puissance, les coûts ainsi que la nécessité d'un complexe multi-accélérateur sont en faveur de la technologie cyclotron. Le FFAG, bien qu'un candidat potentiel pour cette application, n'a pas encore démontré le faisceau de haute puissance. Cela nécessite plusieurs améliorations des installations existantes. Dans KURRI, un démonstrateur de la technologie ADSR est construit là où l'accélérateur était basé sur la technologie FFAG de mise à l'échelle. Comme indiqué dans textbf chapter 3, une campagne de simulation établie en 2014 pour comparer les différents codes a donné un excellent accord. En outre, l'analyse du profil de terrain a révélé certaines principales sources d'imperfections liées aux variations radiales et azimutales de l'indice de champ moyen des aimants.

Afin de comprendre les effets de ce dernier sur la dynamique du faisceau, des modèles analytiques et de simulation ont été développés au chapter 4, ce qui a démontré que les variations de l'accord obéissent à une loi bien définie qui implique la non-mise à l'échelle des orbites ainsi que l'indice de champ moyen de l'aimant. En raison de son comportement monotone avec κ, la différence de l'indice de champ moyen des aimants de focalisation et de défocalisation (propriété établie au chapitre 4), on a pu expliquer et quantifier l'impact des défauts de champ sur la dynamique transversale du faisceau du FFAG de KURRI. Un remède à ces défauts de champ a été proposé: ceci repose sur l'idée d'alterner κ. La variation des oscillations betatron avec l'énergie est réduite par au moins un ordre de grandeur lorsque ce schéma est utilis'e. C'est la pierre angulaire du nouveau concept FFAG que nous présentons et discuterons dans la dernière partie du chapitre 4. Une découverte clé de cette étude est que les conditions cardinales d'un FFAG à focalisation invariante ne sont pas nécessaires pour éliminer les variations du nombre d'onde. En résumé, ce nouveau concept FFAG garantit une forte focalisation en alternant le gradient, mais assure également que la focalisation est constante en alternant les différences azimutales de l'indice de champ moyen des aimants. Une analyse plus approfondie de ce concept a révélé que l'acceptance dynamique diminue avec κ. Cela signifie que de telles valeurs doivent être bien contrôlées afin de maintenir une large acceptance du faisceau nécessaire aux applications de haute puissance.

Le chapitre 5 est dédié pour ce dernier: afin d'entreprendre les simulations de dynamique du faisceau de haute intensité, on a utilisé l'approximation de la lentille mince pour implémenter la charge de l'espace dans le code de simulation ZGOUBI. Plusieurs tests ont été effectués pour valider les résultats de la simulation.

  fieldx ⊥i particle transverse coordinate in the Frenet-Serret frame (xu x + yu y ) same species considered along this paper)

12 CHAPTER

 12

Figure 1 . 1 :

 11 Figure 1.1: Inventory of the fresh and spent fuel in the Uranium dioxide (UOX) fuel cycle.

Figure 1 . 3 :

 13 Figure 1.3: Conceptual diagram of an ADSR.

Figure 1 . 4 :

 14 Figure 1.4: The PHOENIX concept [6].

Figure 1 . 5 :

 15 Figure 1.5: General layout of the Energy Amplifier complex [9].

Figure 1 . 6 :

 16 Figure 1.6: General layout of the ADSR complex in KURRI (2011-present)[15].

Figure 1 . 7 :

 17 Figure 1.7: FFAG accelerator complex connectd to KUCA [15].

Figure 1 . 8 :

 18 Figure 1.8: FFAG complex at KURRI [21]: the extraction line is shown in the background. The 150 MeV FFAG is the larger of the two rings.

Figure 2 . 1 :

 21 Figure 2.1: Core thermal power as a function of k ef f and k1. Here we assume that S = 1 1k ef f so that P th ∝ k1 1k ef f

10 )Fig. 2 .

 102 Fig. 2.2 illustrates the impact of the perturbation of the first neutron generations on the core power change

Figure 2 . 2 :

 22 Figure 2.2: Power variation as a function of k ef f by perturbing the first p terms of the series: p=3 to p=9 are shown here: δP th P th0 = δS S0 is shown by the solid lines. The dashed lines show the effect of the perturbation limited to the generation number p.

2 . 11 )Fig. 2 .Figure 2 . 3 :

 211223 Fig. 2.3 shows the convergence rate of the series for different multiplication factors. For instance, if k ef f = 0.95,

Figure 2 . 4 :

 24 Figure 2.4: Comparison of LCOE costs in four countries: France, Korea, United Kingdom and United States.System costs represent the total costs above plant-level costs to supply electricity at a given load and given level of security of supply[START_REF]The Economics of Nuclear Power[END_REF]. A discount rate of 7% is used.

Figure 2 . 5 :

 25 Figure2.5: LCOE for the Light Water Reactor burning fresh UOX fuel. A 7.6 % discount rate is assumed. The details of the model assumed can be found in[START_REF]The Future of the Nuclear Fuel Cycle[END_REF].

Figure 2 . 6 :

 26 Figure 2.6: Corrected ADS thermal efficiency as a function of the effective multiplication factor and the accelerator efficiency.

Figure 2 . 7 :

 27 Figure 2.7: Corrected ADS thermal efficiency as a function of the accelerator efficiency ηacc for various values of the effective multiplication factor k ef f .

  0.95 and η acc = 20 % k eff = 0.95 and η acc = 40 % Uncorrected thermal efficiency

Figure 2 . 8 :

 28 Figure 2.8: Corrected ADS thermal efficiency as a function of the 1st generation multiplication factor k1 for various values of ηacc.

Figure 2 . 9 :

 29 Figure 2.9: LCOE for ADSR.

Figure 2 . 10 :Figure 2 . 11 :

 210211 Figure 2.10: Monte Carlo simulations with MCNP6[START_REF] Goorley | MCNP6.1.1-Beta Release Notes[END_REF] of the neutron cost from a cylindrical lead target (L=60cm, d=20 cm). Note that below 100 MeV, energy losses due to ionization dominate and the spallation process is quasi-inexistent.

Figure 2 . 12 :

 212 Figure2.12: Spallation neutron spectrum compared to the neutron spectrum of a sodium cooled blanket[START_REF] Heidet | Accelerator-Reactor Coupling for Energy Production in Advanced Nuclear Fuel Cycles[END_REF].

Figure 2 . 13 :

 213 Figure 2.13: Average beam current required as a function of the thermal power of the core.

Figure 2 . 14 :

 214 Figure 2.14: Beam time structure.

Fig. 2 .Figure 2 . 15 :

 2215 Fig.2.[START_REF] Ishi | Status of KURRI Facility[END_REF] shows the two contributions to the rapid decay of the neutron population following an injection of a bunch of neutrons in the blanket: The first term in Eq. (2.32) represents the rapidly decaying transient, often

Figure 2 . 16 :

 216 Figure 2.16: Time variation of the neutron population for various values of the scale factor f .

40 )Figure 2 . 17 :

 40217 Figure 2.17: Domain of variation of the neutron population as a function of the scale factor f .

17

 17

Figure 2 . 18 :

 218 Figure 2.18: Monte Carlo simulations of the effect of the beam size on the target multiplicity for various target radii: a cylindrical lead target of 30 cm length is simulated. The 1 GeV proton beam has a uniform distribution with a square cross section.

Figure 2 . 19 :

 219 Figure 2.19: Beam trip frequency for operating high power proton accelerators [50].

Figure 2 . 20 :

 220 Figure 2.20: Pie chart of the PSI subsystems causing failure [51].

Figure 3 . 1 :

 31 Figure 3.1: Particle position in the Frenet-Serret coordinate system (x, s, y).

Figure 3 . 2 :

 32 Figure 3.2: Accelerating field as a function of time for a positive charge. Example of Stable (S) and Unstable phases (U) are shown as well as Early and Late particles with respect to the stable synchronous phase.

Figure 3 . 3 :

 33 Figure 3.3: Electric field lines in an accelerating gap [68].

Figure 3 . 4 :

 34 Figure 3.4: Injection orbit into the main ring with the matching point (yellow). Only one quadrant is shown here. The two tracking codes, EARLIETIMES and ZGOUBI[START_REF] Méot | Zgoubi Users Guide[END_REF] show good agreement (green and thick black curves superposed). The EARLIETIMES results were provided by Y. Ishi. The color density plot represents the magnetic field associated to the magnets composing the lattice.

Figure 3 . 5 :

 35 Figure 3.5: Betatron tunes from 11 to 139 MeV (left to right) calculated with several codes [64]. The Zgoubi model is in good agreement with the others. The solid lines, in black, show the resonance lines up to the 4th order.

Figure 3 . 6 :

 36 Figure 3.6: Betatron tunes from 11 to 100 MeV (left to right).The details of the measurement as well as the methods and tools that were used to characterize the 150 MeV KURRI FFAG are discussed in[START_REF] Sheehy | Characterization techniques for fixed field alternating gradient accelerators and beam studies using the 150 MeV proton FFAG[END_REF].

Figure 3 . 7 :

 37 Figure 3.7: Revolution frequency vs acceleration time.

  (a) Analytical solution of the longitudinal phase space trajectories including the separatrix (dashed line).The stable fixed point lies at φ = φ s ≡ -π/2[2π] (green point). Longitudinal phase space trajectories including the separatrix as obtained from the tracking code ZGOUBI.

Figure 3 . 8 :

 38 Figure 3.8: Longitudinal phase space trajectories including the separatrix from tracking and comparison with the analytical solution.

Figure 3 . 9 :Figure 3 . 10 :

 39310 Figure 3.9: Comparison of the single particle trajectory in the longitudinal phase space. The SCODE results were provided by S. Machida.

Figure 3 . 11 :

 311 Figure 3.11: Magnetic field for one sector of the ring obtained for different radii.

Figure 3 . 12 :

 312 Figure 3.12: Average scaling factor as a function of the radius (the reference radius is R0=445 cm).

  Local scaling factor as a function of the radius for different azimuthal positions. Field map for one sector in the median plane. X and Y are the Cartesian coordinates in the lab frame with origin the machine centre.

Figure 3 . 13 :

 313 Figure 3.13: Scaling factor as computed from the field map for different azimuthal positions.

Figure 3

 3 Figure 3.14:Average scaling factor for different reference radii: Eq. (3.74) above was applied to transform (R0 = 445cm, B0, k0(R)) into (R1 = 465cm, B1, k1(R)) which gave perfect agreement with the result obtained from the field map. The discontinuity point simply corresponds to the case where R = R1 for which the formula is not defined. (all the other possible transformations in this plot were also checked and gave perfect agreement as expected).

Figure 3 . 15 :

 315 Figure 3.15: Average scaling factor for different reference radii R0.

Figure 3 . 16 :Figure 3 . 17 :

 316317 Figure 3.16: Average scaling factor for the F and D magnets.

Figure 3 . 18 :

 318 Figure 3.18: Average magnetic field as a function of the radius: the field map results are shown in blue while the scaling laws with k constant and variable are shown in red and green respectively. The constant k model shows some discrepancy with the field map data (|ΔB/B| ≤ 4%), dotted red curve, right vertical scale.

Figure 3 . 19 :

 319 Figure 3.19: Flutter function F as a function of the azimuthal angle θ for different radii.

Figure 3 . 20 :

 320 Figure 3.20: Turn separation in the KURRI 150 MeV FFAG (fieldmap) as a function of the kinetic energy and comparison with the analytical formula (3.91).

  √ β b √ γ b ) 1/(√ β b √ γ b ) (a) Position of the single particle trajectory in the vertical plane and comparison with the damping law. √R √ β b √ γ b ) 1/(√ β b √ γ b ) (b) Angle of the single particle trajectory in the vertical plane and comparison with the damping law.

Figure 3 . 21 :

 321 Figure3.21: Single particle trajectory in the vertical plane as a function of the path length from solving the non-linear equation of motion using the tracking code ZGOUBI and the user-implemented scaling "FFAG" analytical model (red). Comparison with the damping laws (3.97 and 3.98) established earlier using the WKB approximation of the linearized equations showed good agreement in both cases (blue). The R-independent damping (that would be the case for instance in a pulsed synchrotron) is shown for comparison (in green)

Figure 3 . 22 :

 322 Figure 3.22: Beam edge radius as a function of the orbit radius from solving the non-linear equation of motion using the tracking code ZGOUBI and the user-implemented scaling FFAG analytical model (red). Comparison with the R-dependent damping law (Eq.(3.99)) established using the WKB approximation of the linearized equations showed good agreement in both cases (blue).

Figure 3 . 23 :

 323 Figure 3.23: Unnormalized emittance change as a function of the orbit radius from ZGOUBI tracking and comparison with the damping law (Eq.(3.101).

Figure 4 . 1 :

 41 Figure 4.1: Magnetic field in the hard edge model.

4 . 2 )Figure 4 . 2 :

 4242 Figure 4.2: Example of trajectory in the hard edge model of a KURRI-like DFD triplet.

4 . 3 .

 43 (a) Horizontal tune as a function of the radius. (b) Vertical tune as a function of the radius.

Figure 4 . 3 :

 43 Figure 4.3: Comparison of the tune calculation using the hard edge model with the formula 4.21 and 4.22.

( 4 .

 4 [START_REF] Joho | High intensity problems in cyclotrons[END_REF] and(4.49)). The hard edge model formula that was established earlier is also tested (see Eqs. (4.21) and (4.21)). The main difference between the Bogoliubov method of averages and the hard edge model formula is that the only parameters that may vary in the latter are the average field indices k F and k D . The results are also compared to the 1st order Symon formula (ν 2 x = k + 1 and ν 2 y = -k + F 2 ) in which the scaling factor k as well as the magnetic flutter F 2 are allowed to change with the radius and were obtained from the closed orbits calculation by Zgoubi. Comparison between the different results are summarized in Figs. 4.4 and 4.5 below.

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Tunes calculation in the horizontal plane of the KURRI 150 MeV FFAG ring and comparison with the analytical formula as well as the measurement.

Figure 4 . 6 :

 46 Figure 4.6: 2 nd order index of similarity of the closed orbits β.

Figure 4 . 7 :

 47 Figure 4.7: Magnetic field along several closed orbits.

Figure 4 . 8 :

 48 Figure 4.8: Example of several closed orbits for a scaling FFAG .

  Tune per cell of the DFD triplet in the horizontal plane. Tune per cell of the DFD triplet in the vertical plane.

Figure 4 . 9 :

 49 Figure 4.9: Contour plot of the tune variations as a function of the FD ratio and the average field index of the magnets.

Figure 4 . 10 :

 410 Figure 4.10: Stability diagram as a function of the FD ratio and the average field index k: the green line shows the value of the FD ratio that we choose for the study that follows. Note that the region on the left side is limited by the physical size of the magnets, which are limited to Rmax = 10 m.

  4.10. A scan on k F and k D provided the stability diagrams of the DFD triplet in the transverse plane (see Figs. 4.11 and 4.12 ). Qualitatively, it shows that, in the case where k F = k D = k, the average cell tune exhibits the right behavior predicted by the Symon formula (4.20): increasing k increases the horizontal tune and decreases the vertical one. Besides, the RMS tune exhibits the expected behavior in the vicinity of the line k F = k D where it becomes negligible.

  Average tune of the DFD triplet in the horizontal plane.

  RMS tune variations of the DFD triplet in the horizontal plane.

Figure 4 . 11 :

 411 Figure 4.11: Contour plot of the average and rms tune variations in the horizontal plane as a function of the scaling factor kF and kD.

�

  Average tune of the DFD triplet in the vertical plane.

  RMS tune variations of the DFD triplet in the vertical plane.

Figure 4 . 12 :

 412 Figure 4.12: Contour plot of the average and rms tune variations in the vertical plane as a function of the scaling factor kF and kD: as can be observed from this plot, whenever kF 6 = kD the rms tune variations increase considerably.

  4.11b and 4.12b). Modifying the colorbar range in Fig. (4.11b) enabled us to clearly observe the rms tune variations in the inner region of the stability diagram (see Fig. (4.13)

Figure 4 . 13 :

 413 Figure 4.13: RMS tune variations of the DFD triplet in the horizontal plane. (Same result as Fig. 4.11b, but with a modified colorbar range).

  (b) Average tune per cell in the vertical plane.

Figure 4 . 14 :

 414 Figure 4.14: Average tune variations per cell in the transverse plane as a function of kF or kD. A scan on kF or kD is performed by fixing kD or kF respectively, to the design value 7.6.

  (b) RMS tune variations of the DFD triplet in the vertical plane as a function of kF or kD..

Figure 4 . 15 :

 415 Figure 4.15: RMS tune variations per cell in the transverse plane as a function of kF or kD.

FD

  increases with the energy =⇒ θ F increases with the energy =⇒ θ D increases with the energy =⇒ α D (E) = R ρ of similarity <R/ρ> Kinetic Energy (MeV) (k F ,k D )=(7.6,7.0) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,8.2) (k F ,k D )=(8.2,7.6) (k F ,k D )=(8.2,8.2) (k F ,k D )=(8.2,8.8)

Figure 4 .

 4 Figure 4.16: 1 st order index of similarity of the closed orbits.

  ,k D )=(7.6,7.0) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,8.2) (k F ,k D )=(8.2,7.6) (k F ,k D )=(8.2,8.2) (k F ,k D )=(8.2,8.8) (a) Index of similarity of the F-magnet. ,k D )=(7.6,7.0) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,8.2) (k F ,k D )=(8.2,7.6) (k F ,k D )=(8.2,8.2) (k F ,k D )=(8.2,8.8)(b) Index of similarity of the D-magnet.

Figure 4 .

 4 Figure 4.17: 1 st order index of similarity of the focusing and defocusing magnets.

Figure 4 . 18 :Figure 4 . 19 :

 418419 Figure 4.18: Closed orbits patterns for two cases: scaling (red) vs non-scaling FFAG (blue) in which kF 6 = kD.

Figure 4 .

 4 Figure 4.20: 2 nd order index of similarity of the closed orbits β.

Figure 4 . 21 :

 421 Figure 4.21: Magnetic flutter F 2 as a function of the kinetic energy.

  ,k D )=(7.6,7.0) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,8.2) (k F ,k D )=(7.6,7.0) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,8.2)

Figure 4 .

 4 Figure 4.22: Example of tune calculation using the tracking code ZGOUBI (blue) and comparison with the 1st order analytical formula 4.47 (red). To the first order in the method of averages, good agreement (< 4%) is obtained in the cases (kF , kD) = (7.6, 7.0) and (kF , kD) = (7.6, 7.6). However, for the last case, i.e (kF , kD) = (7.6, 8.2) the large discrepancy observed (� 16%) is due to the fact that the beam is located nearby the boundary of stability.

Figure 4 . 23 :

 423 Figure 4.23: Example of tune calculation using the tracking code ZGOUBI (blue) and comparison with the 1st order analytical formula 4.43 (red). To the first order in the method of averages, the differences are less than 20%. The agreement could be better if the contribution of the second order term were taken into account (given that g2 > 0, see Eq. 4.36).

Figure 4 . 24 :

 424 Figure 4.24: Comparison of the analytical formula (4.58) with the ZGOUBI tracking results. The first order linear approximation, i.e ax,yκ is also shown in blue.

Figure 4 . 25 :

 425 Figure 4.25: Correction scheme based on alternating κ value. correction scheme, tracking in ZGOUBI was performed using built-in fieldmaps of DFD triplets where the average field indices of the F and D magnets, k F and k D , can be adjusted as illustrated in Fig. 4.7. The results for two different test cases are shown in Figs. 4.26 and 4.27. As summarized in Table4.1, implementing the correction scheme does reduce the tune variations in both cases by at least an order of magnitude: Case1: Before correction, we choose (k F 0 , k D0 ) = (7.6, 7.5). After correction, we choose k D1 = 7.7. The results are summarized in Figs.(4.26) and (4.27) below. Case2: Before correction, we choose (k F 0 , k D0 ) = (7.6, 7.8). After correction, we choose k D1 = 7.5. The results are summarized in Figs.(4.28) and (4.29) below. Also, note that these results were not finely optimized. A step size of δk = 0.1 was employed to search for the minimum RMS tune variations, which can be further reduced.

Figure 4 . 26 :

 426 Figure 4.26: Case 1: Tune variations as a function of the energy before and after correction: the corrected scheme is shown in red where (kF 0, kD0) = (7.6, 7.5) and (kF 1, kD1) = (7.6, 7.7).

Figure 4 . 27 :

 427 Figure 4.27: Case 1: Tune variations as a function of the energy before and after correction: the corrected scheme is shown in red where (kF 0, kD0) = (7.6, 7.5) and (kF 1, kD1) = (7.6, 7.7).

Figure 4 . 28 :

 428 Figure 4.28: Case 2: Tune variations as a function of the energy before and after correction: the corrected scheme is shown in pink where (kF 0, kD0) = (7.6, 7.8) and (kF 1, kD1) = (7.6, 7.5).

Figure 4 . 29 :

 429 Figure 4.29: Case 2: Tune variations as a function of the energy before and after correction: the corrected scheme is shown in pink where (kF 0, kD0) = (7.6, 7.8) and (kF 1, kD1) = (7.6, 7.5).

Figure 4 . 31 :

 431 Figure 4.31: Magnetic field of the analytical model "FFAG" along several closed orbits.

Figure 4 . 32 :Figure 4 . 33 :

 432433 Figure 4.32: Tune variations as a function of the energy for various lattices with alternating κ.

Figure 4 .

 4 Figure 4.34: Closed orbits in a scaling FFAG (black) and in a fixed tune non-scaling FFAG with alternating κ (pink and light blue). For the sake of clarity, the distance between the closed orbits of the scaling and the non-scaling FFAG is amplified.

  of similarity <R/ρ> Kinetic Energy (MeV) (k F ,k D )=(7.6,7.400) (k F ,k D )=(7.6,7.748) (k F ,k D )=(7.6,7.4),(7.6,7.748)

Figure 4 .

 4 Figure 4.35: Index of similarity α of the orbits in the advanced FFAG concept: down-pointing and upward-pointing triangle show α averaged over the width of the 1 st and 2 nd cell respectively, while the cross points represent α averaged over the width of the two cells combined.

Figure 4 . 36 :

 436 Figure 4.36: Average field index in the advanced FFAG concept.

  ,k D )=(7.6,7.400) (k F ,k D )=(7.6,7.748) (k F ,k D )=(7.6,7.4),(7.6,7.748)

Figure 4 . 37 :

 437 Figure 4.37: Magnetic flutter in the advanced FFAG concept.

  Figure 4.38: 2 nd order index of similarity in the advanced FFAG concept.

  ,k D )=(7.6,8.0) (k F ,k D )=(7.6,7.4) (k F ,k D )=(7.6,7.7) New conjecture

Figure 4 . 39 :

 439 Figure 4.39: DA of the FFAG lattice as a function of the energy for various values of kD.

Figure 4 . 40 :

 440 Figure 4.40: DA of the FFAG lattice as a function of the energy for various values of kD: the points represent the results of tracking simulations while the solid lines are the results of the above conjecture (Eq. (4.64)).

Figure 4 . 41 :

 441 Figure 4.41: Schematic of the concept of the 2D barycentric interpolation: the black point is the working point where one needs to estimate the value of the DA while the red, green and white points represent the FFAG lattices for which the values of the DA are known.

Figure 4 . 42 :

 442 Figure 4.42: DA of the three basis functions of the barycentric interpolation (orange, green and blue), the conjectured values (solid lines) as well as the corresponding calculated DA from tracking.

  ,k D )=(7.6,7.8) (k F ,k D )=(7.6,7.6) (k F ,k D )=(7.6,7.4) (k F ,k D )=(7.6,7.4),(7.6,7.8)

Figure 4 . 43 :

 443 Figure 4.43: DA of the alternating-κ FFAG and comparison with the scaling case.

  Case of a scaling FFAG with kF = kD = 7.6

Figure 4 . 45 :

 445 Figure 4.45: DA of the alternating-κ FFAG as a function of κ.

Figure 5 . 1 :

 51 Figure 5.1: Space charge electric field in a KV beam.

Figure 5 . 2 :

 52 Figure 5.2: Horizontal space charge electric field in a Gaussian beam. The black line shows the projection on the plane y = 0.

FigFigure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Horizontal beam size as a function of the path length for RMS equivalent beams.

Figure 5 . 5 :

 55 Figure 5.5: Beam RMS emittance as a function of the path length for RMS equivalent beams. Same notation as in the previous plot.

Figure 5 . 6 :

 56 Figure 5.6: Radial beam envelope in a drift.

  various initial conditions and perveance values yields the results shown in Fig 5.8 below.

Figure 5 . 7 :

 57 Figure 5.7: Convergence of the space charge calculation in a drift for different slicing length.

Figure 5 . 8 :

 58 Figure 5.8: Radial beam envelope in a drift: the ZGOUBI results are shown in solid lines while the points depict the analytical solution of the envelope equation.

Figure 5

 5 Figure 5.9: Horizontal phase space at injection for the KURRI 150 MeV FFAG: as can be observed, when the initial amplitude increases, the motion invariant is deformed (here, by proximity of a 3rd integer tune). The limit of stable motion is obtained when the particle is lost in less than 10000 turns. This represents that largest phase-space area that will pas through the FFAG lattice, also referred to as the horizontal aperture (black curve).

. 51 )

 51 ��� ��� ��� ��� ��� ��� ��� ��� ��� ���

Figure 5 . 10 :

 510 Figure 5.10: Contour plot of the space charge tune shift as a function of the average scaling factor k and the perveance Q (from tracking).

Figure 5 . 11 :

 511 Figure 5.11: Contour plot of the predicted space charge tune shift as a function of the average scaling factor k and the perveance Q (analytical formula 5.50 & 5.51).

Figure 5 . 12 :

 512 Figure 5.12: Maximum current that can be transported through the scaling FFAG lattice as obtained from Eq. 5.56.

Figure 5 . 13 :

 513 Figure 5.13: Tune variations R y σν as a function of the k factor and the perveance term.

Figure 5 . 14 :

 514 Figure 5.14: Contour plot of the space charge tune shift in the horizontal plane.

Figure 5 . 15 :Figure 5 . 16 :

 515516 Figure 5.15: Contour plot of the space charge tune shift in the vertical plane.

Figure 5 . 17 :

 517 Figure 5.17: Tune per cell before and after perturbation of the average field index of the D-magnet: before correction (kF , kD ) = (7.6, 7.6) while after correction (kF , kD ) =(7.6, 7.11).

  5.21).

Figure 5 . 18 :

 518 Figure 5.18: Emittance evolution of the first 100 turns (without acceleration) for a scaling FFAG. The labels correspond to the single bunch number of particles.

  (a) Density plot of the Gaussian distribution in the transverse plane. (b) Density plot of the tune spread.

Figure 5 . 19 :Figure 5 . 20 :

 519520 Figure5.19: Contour plot of the particle distribution and its equivalent tune spread: in red the density is maximum while in blue the particle density is minimum. As can be seen, particles at the centre of the distribution (red) experience the largest tune shit.

Figure 5 . 21 :

 521 Figure 5.21: Space charge tune spread for various lattices: the scaling case (in red) and the non-scaling case where the bare tune is fixed by alternating κ (in pink).

  Fig 25a: in red, the spiral angle is fixed, ξ = 30 deg, while, in blue, the spiral angle increases with the radius, ξ = 30 + 0.02 × R deg. The tune calculation is shown in Fig 25b, were one can see that increasing the spiral angle allows to increase the vertical tune. deg orbit for ξ=30+0.02*r deg (a) Plot of the orbits as well as the EFB of the entire ring. Tune of the cyclotron with and without increasing spiral angle.

Figure 25 :

 25 Figure 25: Comparison of the cyclotron analytical models: with fixed spiral angle and with varying spiral angle.

Figure 26 :

 26 Figure 26: Contour plot of the tune variations as a function of the FD ratio and the average field index of the magnets.

Figure 27 :

 27 Figure 27: Stability diagram as a function of the FD ratio and the average field index k.

Table 1 .

 1

		Charge	Discharge
	U(ton)	29.67	27.77
	Pu(ton)	0	0.35
	MA(ton)	0	0.03
	Total(ton)	29.67	28.15
	Fissile fraction (%)	4.2	1.64

1: Mass flow of PWR for a reactor power of 1000 MWe and a cycle length of 1.5 years.

Table 1 . 2 :

 12 Inventory of an irradiated LWR fuel 10 years after discharge. Nota Bene: these proportions vary depending on the cooling time after discharge. tHM refers to tons of Heavy Metal.

	Z
				163 d	29 y	18 y
	-� Uranium cycle
	-� Thorium cycle			432 y	16 h	7370 y	10 h
		88 y	fissile	6571 y	14 y	375 000 y	5 h
		2.14 �� � y	2.12 d	2.36 d
	fissile	fissile		fertile	23 min
	27 d
	fertile	22 min
	Cm		0.02 kg/tHM
	Am		0.6 kg/tHM
	Pu		8.5 kg/tHM
	Np		0.5 kg/tHM
	U		955.4 kg/tHM
	LLFP		2 kg/tHM

Figure 1.2: NZ diagram of the actinides formation in the Uranium and Thorium fuel cycle. N (resp. P) is the number of neutrons (resp. protons) in the nucleus. The most problematic elements that contribute to the long term radiotoxicity of the waste are framed in yellow.

Table 1 .

 1

	Specifications of the 150 MeV FFAG
	Focusing structure	DFD
	N cells , number of cells	12
	Injection energy	11 MeV
	Extraction energy	100 to 150 MeV
	Average radius hRi Field index k	4.57 → 5.4 m 7.6
	f RF , RF frequency V RF , RF peak voltage	1.6 → 5.2 MHz 4 kV
	Synchronous phase	-30 deg
	Repetition rate (max)	30 Hz
	Average current at injection (from the Linac)	600 nA
	Average current at extraction	1.5 nA
	Injection efficiency	40 %
	Transmission up to 1 ms in the acceleration cycle	2.5 %
	Transmission from 1 ms to extraction timing	25 %
	Extraction efficiency (assumed)	100 %
	Overall transmission	0.25 %
	Specifications of the sub-critical core
	Spectrum	Thermal
	k ef f	0.99
	P th,c , thermal power of the core	10 W
	P th,c max, maximum thermal power allowed	100 W

3: Table of parameters of the ADSR in KURRI.

Table 2 . 1 :

 21 Table of efficiencies for the PSI's high intensity proton cyclotron (P beam ≈ 1.3M W )

	ADSR dedicated

Table 2 . 2 :

 22 Capacity Factors for Utility Scale Generators in the USA [35]: record reached in 2015 for the nuclear industry.

	.2, the

Table 2 .

 2

3: Repetition rate for various reactor concepts.

Table 2 . 4 :

 24 List of parameters for ADSR concept.

  .[START_REF] Ludewig | Study of Multi-Beam Accelerator Driven Thorium Reactor[END_REF]) for the vertical cell tune, this simply means that the term involving the magnetic flutter

	RMS horizontal tune	0.002 0.004 0.006 0.008 0.01 0.012					kD fixed kF fixed		RMS vertical tune	0.001 0.002 0.003 0.004 0.005 0.006 0.007 0.008 0.01 0.009					kD fixed kF fixed
		0								0
		7	7.2	7.4	7.6	7.8	8	8.2	8.4	7	7.2	7.4	7.6	7.8	8	8.2	8.4
					kF, kD							kF, kD
	(a) RMS tune variations of the DFD triplet in the
	horizontal plane as a function of kF or kD.

  10.7 ; A 4 ≈ 0.0088(4.60) Comparison between the analytical formula (4.58) and the tracking results is shown in Fig 4.24. As can be observed, the agreement is quite good except near the boundary of stability, i.e. for large κ values.

		0.016				Analytical
		0.014		Analytical 1st order ZGOUBI Tracking results
	RMS horizontal tune per cell	0.004 0.006 0.008 0.01 0.012
		0.002
		0
		7.6	7.8	8	8.2	8.4	8.6	8.8
					kD

Table 4 . 1 :

 41 Comparison of the tune variations per cell before and after the implementation of the correction scheme.

  Compte tenu de la propriété de champ fixe, le FFAG combine les propriétés du cyclotron et du synchrotron. Il existe deux types différents de FFAG: le FFAG à focalisation invariante où le nombre d'oscillations betatron reste constant pendant l'accélération, et le FFAG à focalisation non-invariante où le nombre d'oscillations betatron varie avec l'énergie. Le premier type a été proposé indépendamment au milieu des années 1950 par T. Ohkawa au Japon, K. Symon aux États-Unis et A. Kolomensky en Russie. La découverte du gradient alternatif (AG) est survenue en 1952. Il s'en est suivi deux périodes de développement des FFAGs: la première, dans le Wisconsin, par la MURA (Midwestern Universities Research Association), où trois machines à électrons ont été construites et testées entre 1956 et 1967. Puis, après une période de silence de près de 30 ans avec un intérêt limité dans la technologie FFAG en raison de plusieurs défis. Les FFAG ont finalement été relancés à la fin des années 1990 au Japon dans le cadre de l'usine de neutrinos et des applications médicales: un prototype de proton a été construit et exploité en 1999 par le groupe de Mori chez KEK (machine de 500 keV, 2,5 m de diamètre). De plus, au cours des années suivantes, plus de 6 autres FFAG ont été construits au Japon. La caractérisation et les expériences sont toujours en cours.

The discovery of the Alternating Gradient (AG) occurred in 1952.

A rigorous approach to investigate the kinetics of subcritical multiplying systems relies on the adjoint flux calculation[START_REF] Lewins | Importance: The Adjoint Function: The Physical Basis of Variational and Perturbation Theory in Transport and Diffusion Problems[END_REF], where one uses a weighting function to evaluate the "importance" of producing fission neutrons.

The pulsed efficiency is defined as the ratio Beam pulse to AC power pulse length: η pulsed = τ beam τ modulator

If one were to start from the proton bunch in the calculation above, the result will be a double convolution.

From now on, unless otherwise specified, the magnetic field B refers to the vertical component of the field (that is also denoted Bz).

This result is a generalization to what can be found in the Symon paper[START_REF] Symon | Fixed-Field Alternating-Gradient Particle Accelerators[END_REF]. We introduce the radial dependence of the tune. In other words, no assumption is made regarding the constancy of the tune or the similarity of the orbits. Also, we distinguish between the F and D magnets, which is equivalent to introducing an azimuthal variation of the scaling factors.

One could instead use the |max -min| value of the tunes to account for the oscillations. However, the rms quantities have the merit to be average quantities, thus more appealing to use in order to obtain smooth variations of the described quantities.

Only externally applied forces are considered here, i.e no collective effects such as the space charge effects.

Another stability island exists for larger k values[START_REF] Machida | Scaling Fixed-Field Alternating Gradient Accelerators with a Small Orbit Excursion[END_REF], but we restrict our analysis to a KURRI type FFAG for which the design value of k ≈ 7.6

Definition: the concept of a scaling FFAG is based on the assumption of the similarity of the orbits at each azimuthal position i.e R/ρ = const. We refer to the case where this condition is violated as the non-scaling FFAG.

In our treatment, we assume only one particle species in a perfect vacuum.

Acknowledgements CHAPTER 5. SPACE CHARGE EFFECTS IN FIXED FIELD ACCELERATORS Bibliography

Appendices

Appendix B: Comparison of the tunes from the accelerated and fixed orbit

Appendix C: "Cyclotron" with varying spiral angle

The newly implemented [START_REF] Tahar | FFAG14 workshop tutorial[END_REF] element "CYCLOTRON" in zgoubi provides a model of a dipole field. The field along the particle's trajectory is computed as the particle motion proceeds, by using the magnet's geometrical boundaries: At any position P(R,θ) along the particle trajectory (see Fig 24), the value of the vertical component of the mid-plane field is calculated using:

(according to the Enge fringe model) where:

• B norm is a normalization coefficient,

• F(R, θ) is the fringe field coefficient, given by:

where The effective field boundaries are modelled by a logarithmic spiral for which the angle ξ is allowed to increase radially, namely:

Appendix E: Collisional/collective regime

In order to evaluate which regime governs the particle interactions, one can give some estimates of the RHS term in equation (5.12):

In a heuristic approach [START_REF] Barnard | Introduction[END_REF], one can write

where ν c is the collision frequency (the inverse of the time for a particle to suffer a collision): ν c ∼ σnv ∼ πr 2 c nv and σ is the collision cross section (cm 2 ), n the particle 3 density (cm -3 ).

where ω p is the angular frequency of plasma oscillations. This yields:

For a large scattering angle, the kinetic energy of a particle is of the order of the potential energy at closest approach. Thus the collision radius can be defined by [START_REF] Barnard | Introduction[END_REF]:

where

and T is the beam temperature in the laboratory frame. It results that

Collision effects

Collective effects

where

is the number of particles within a Debye sphere of radius λ D . Therefore, if N D ≫ 1, the collisional effects can be neglected compared to the collective effects and the smooth functions employed to calculate the space charge effects. The previous definitions only apply to a non-relativistic charged particle. However, these results can be generalized if we replace the beam temperature in the laboratory frame T by the beam temperature in the beam rest frame T b [START_REF] Reiser | Theory and Design of Charged Particle Beams[END_REF]. The relationship is T b = γT . Also replacing m by γm in the previous expressions yields:

where all quantities are defined in the lab frame (except T b ).

Based on all the above, the effect of the Coulomb interaction between N particles in an accelerator can be classified into two regimes:

(i) collisional regime: for instance intra-beam scattering and collisions with the gas molecules: in that case the particle sees its immediate neighbours and the forces are dominated by the binary, even tertiary collisions. Therefore, the fluctuations induced by these discrete effects, δf , δE and δB cannot be neglected. Such effects play a major role in charge neutralization. For all beams of interest for our study, the collisional forces are neglected.

(ii) collective regime or what we refer to in general as space charge regime: in that case, the interaction is smoothed-out and the bunch as a whole acts on each particle. This approximation is valid as long as the Debye length remains small compared to the bunch size.

Appendix F: Résumé en français F.1