Yu-Chao Jiang

Some issues in the physics and particle simulations of low temperature partially magnetized plasmas and ion sources

Chapter 1 Introduction

Our imagination is stretched to the utmost, not, as in fiction, to imagine things which are not really there, but just to comprehend those things which 'are' there.

Richard Feynman

This manuscript is devoted to the study of some specific issues related to the physics and simulation of low temperature partially magnetized plasmas. These plasmas are typical of ion sources such as those used in ion thrusters (gridded ion thrusters, Hall thrusters, Highly Efficient Multistage Plasma Thrusters -HEMPT, . . . ), in negative ion sources for neutral beam injection in fusion or in magnetron discharges for processing applications.

In partially magnetized plasmas electrons are strongly magnetized (electron Larmor radius much smaller than the plasma dimensions) while ions are not (ion Larmor radius larger than the plasma dimensions). This gives these plasmas very specific properties different from those of space plasmas or fusion plasmas. Since electrons are confined by the magnetic field, their residence time in the system is increased and they can undergo collisions with neutral atoms or molecules and ionization, allowing plasma maintenance even at very low pressure. Ions, on the other hand, are not sensitive to the magnetic field, can be practically collisionless if the pressure is low enough, and can be accelerated by electric fields that can be present in the cathode sheath, as in magnetron discharges, in the plasma volume, as in Hall thrusters, or between a plasma grid and an extracting grid as in gridded ion sources or in negative ion sources for fusion.

As said above, the magnetic field is used in these plasma sources to confine the electrons. Electrons can be efficiently confined if magnetic cusps are placed on the chamber walls and/or in so-called E × B configurations. E × B configurations correspond to conditions where an external magnetic field is placed perpendicular to the applied electric field or to the discharge current. For efficient confinement, such configurations must have a cylindrical symmetry with E × B in the azimuthal direction (i.e. E axial and B radial or E radial and B axial). In these conditions the E × B drift current is closed on itself and does not go to a wall (devices with E × B configurations are also called "closed drift" devices). Magnetic cusps are used to limit the electron and ion flow to the walls. They are generally created by magnets placed along the chamber walls so that the magnetic field is parallel to the walls over a large part of the wall and perpendicular to the wall over a smaller part. This is equivalent to reducing the effective wall area viewed by the electrons and therefore to reducing the charged particle losses to the walls. Classical theory of charged particle transport predicts that electron transport across a uniform magnetic field is possible only if the electrons undergo collisions. This is called "classical, collisional cross-field transport". For example, in an E × B configuration, the classical collisional transport theory predicts that while electron mobility or diffusion coefficient without magnetic field are inversely proportional to the electron collision frequency, these coefficients are proportional to the collision frequency in the direction parallel to the electric field and perpendicular to the magnetic field when a magnetic field is present and inversely proportional to the square of the magnetic field. In practice, when the electron collision frequency is much smaller than the electron cyclotron angular frequency (large Hall parameter), crossfield electron transport is larger than predicted by the classical collisional theory because of the formation of instabilities and turbulence. Instabilities and turbulence are ubiquitous in partially magnetized plasmas as well as in space plasmas or fusion plasmas, but they are very specific because electrons are strongly magnetized while ions are practically not magnetized. Due to these instabilities, electron transport across the magnetic field is larger than predicted by the classical collisional theory and this "anomalous electron transport" is difficult to quantify.

In this thesis we will illustrate some of the issues in the physics and modeling of partially magnetized plasmas with three specific examples that correspond to ongoing studies in the GREPHE group of the LAPLACE laboratory:

1. Electron extraction in negative ion sources for neutral beam injection in fusion 2. Instabilities in magnetron discharges and Hall thrusters CHAPTER 1. INTRODUCTION

Plasma confinement by magnetic cusps

In the study of the negative ion source for fusion, the aim of the GREPHE group is to better understand the physics of the negative ion source, and more specifically, the questions of plasma transport across the magnetic filter and of negative ion extraction from the plasma. One of the important issues in these negative ion sources is to minimize the current of electrons that are coextracted with the negative ions. In this thesis we focus on this aspect and we try to understand and quantify how electrons can be extracted through a grid aperture when a magnetic cusp is placed in front of the aperture. We will show with the help of 3D PIC simulations that E × B electron drift in the meniscus region plays an important role in this process.

Hall thrusters and magnetron discharges are E × B cylindrical devices with radial magnetic field and axial electric field. It has been known for a long time that instabilities are present in the discharges of Hall thrusters and magnetrons, leading to important anomalous electron transport. In this thesis we focus on one particular type of instability which is known to be present in Hall thrusters and magnetron discharges and is apparent in the experiments as a luminous non-uniformity rotating in the azimuthal direction. This macroscopic instability, called "rotating spoke" in the literature has been recently analyzed by the GREPHE group with 2D axial-azimuthal Particle-In-Cell Monte Carlo Collisions (PIC MCC) simulations. In this thesis we have used such a simulation to perform a parametric study of this instability.

Magnetic cusps have been used for more than 60 years to confine the plasma in a large variety of conditions. An important parameter characterizing plasma confinement by cusps is the effective loss area in the presence of magnetic cusps. Some semi-empirical theories have been proposed to quantify the effective loss area and their predictions have been compared with numerous experimental results. In spite of these efforts there is no fully reliable expression of the effective wall loss as a function of different parameters such as magnetic field, electron temperature, ion mass, gas pressure, etc. . . The availability of powerful simulation tools and computers should allow us to progress toward a more quantitative description of the confinement by magnetic cusps. We describe in this thesis an attempt at obtaining scaling laws for the effective loss width of magnetic cusps, based on 2D PIC MCC simulations.

The manuscript is organized as follows. In the remaining of this introduction we provide more details on the general context of this work and on the three topics that CHAPTER 1. INTRODUCTION we are addressing. Chapter 2 presents the basics principles of particle simulations in the context of magnetized plasmas and discusses some aspects of the accuracy of these simulations. Chapter 3 describes our attempts at obtaining scaling laws for magnetic cusps confinement in low β plasmas. The question of electron extraction across the grids of a negative ion source is addressed in chapter 4. Chapter 5 is devoted to the parametric study of an ionization instability predicted by 2D PIC MCC simulations and that seems to correspond to the observed rotating spokes in magnetron plasmas or Hall thrusters.

Negative ion source for neutral beam injection in fusion 1.Principles of magnetic fusion

Nowadays, energy consumption increases every year. With the increasing population, more energy is needed which makes finding new generation energy sources a primary task for humanity. As shown in figure 1.1, there is a significant difference between different energy producing methods. For coal and oil, it costs million tonnes to produce the same amount of energy which can be harvested by fusion only using 100 kg Deuteron and 150 kg Triton. Based on this tremendous advantage, sustainable fusion reaction has been studied since the early 20th century. We say: "Without fusion, there would be no life on earth", because the light and the warmth we feel right now is the results of fusion reaction in the core of the sun: hydrogen nuclei collide into heavier helium atoms and release huge amounts of energy during the process. For fusion reaction in the core of sun, the gravitation force confines high energy nuclei with temperature of about 15 million degrees Celsius. But on earth, such highly demanding confinement can not be conducted by gravitation. There are two ways designed by scientists to confine particles in a typical fusion device: inertial confinement fusion and magnetic confinement fusion.

Inertial confinement fusion(ICF) is a type of fusion energy research that attempts to initiate nuclear fusion reactions by heating and compressing a fuel target, typically in the form of a pellet that most often contains a mixture of deuterium and tritium. Generally ICF systems use a single laser, the driver, whose beam is split up into a number of beams which are subsequently individually amplified a trillion times or more. These are sent into the reaction chamber (called a target chamber) by a number of mirrors, positioned in order to illuminate the target evenly over its whole surface. The heat applied by the driver causes the outer layer of the target to explode, just as the outer layers of an H-bomb's fuel cylinder do when illuminated by the X-rays of the fission device.The material exploding off the surface causes the remaining material on the inside to be driven inwards with great force, eventually collapsing into a tiny near-spherical ball. In modern ICF devices the density of the resulting fuel mixture is as much as one-hundred times the density of lead, around 1000 g/cm 3 . This density is not high enough to create any useful rate of fusion on its own. However, during the collapse of the fuel, shock waves also form and travel into the center of the fuel at high speed. When they meet their counterparts moving in from the other sides of the fuel in the center, the density of that spot is raised much further. Given the correct conditions, the fusion rate in the region highly compressed CHAPTER 1. INTRODUCTION by the shock wave can give off significant amounts of highly energetic alpha particles. Due to the high density of the surrounding fuel, they move only a short distance before being "thermalized", losing their energy to the fuel as heat. This additional energy will cause additional fusion reactions in the heated fuel, giving off more highenergy particles. This process spreads outward from the center, leading to a kind of self-sustaining burn known as ignition. In magnetic confinement fusion, magnetic fields are used to confine the fusion fuel in a plasma to produce thermonuclear fusion power. Charged particle gyrates around the magnetic field line, so it is possible to constrain particles with magnetic field. In the early days of research, magnetic mirror was used to confine the plasma. A basic magnetic mirror is shown in figure 1.3. The confinement of a magnetic mirror is the result of conservation of magnetic moment µ:

µ = W ⊥ B = 1/2mv 2 ⊥ B = const (1.7)
where W ⊥ is the energy perpendicular to the magnetic field. Assuming we have a particle going to the left with initial velocity of v 0 ⊥ and v 0 (velocity perpendicular and parallel to the magnetic field) at B min , and assuming the particle energy is CHAPTER 1. INTRODUCTION conserved during the motion in magnetic field. So:

µ = W ⊥ B = 1/2mv 2 ⊥ B min = const (1.8) and W = 1/2m(v 2 ⊥ + v 2 ) = W 0 = const (1.9)
As particle moves to the left, B increases, so based on equation 1.8, W ⊥ = 1/2mv 2 ⊥ increases. If W ⊥ ever reaches W 0 where magnetic field is B max , there will be no motion parallel to the magnetic field. So particles will be reflected back to the right. This happens if:

µ = 1/2mv 0 ⊥ 2 B min > W 0 B max = 1/2m(v 0 ⊥ 2 + v 0 2 ) B max (1.10) or v 0 ⊥ 2 v 0 ⊥ 2 + v 0 2 > B min B max (1.11)
Defining the pitch angle θ = tan -1 (v 0 ⊥ /v 0 ), we have from equation 1.11:

sinθ > B min B max 1/2 = 1 R 1/2 (1.12)
the mirror ratio R = B max /B min . Particles with pitch angle satisfying equation 1.12 will be contained by the magnetic mirror. Later on, toroidal machines like z-pinch device, Stellarators, Tokamak were developed to study sustainable fusion reaction. The tokamak device was designed in USSR in 1951 by physicist Andre Sakharov and Igor Tamm [START_REF] Smirnov | Tokamak foundation in ussr/russia 1950-1990[END_REF]. A tokamak device is usually surrounded by two different magnetic field configurations: toroidal field and CHAPTER 1. INTRODUCTION poloidal field. The toroidal field is generated by sets of coils, and an inner poloidal field is produced due to the plasma current, an additional poloidal field is generated by coils in order to maintain the shape and stability of the plasma inside the device.

Figure 1.4: Sketch of tokamak magnetic field, including toroidal magnetic field and poloidal magnetic field [START_REF] Theiler | Basic investigation of turbulent structures and blobs of relevance for magnetic fusion plasmas[END_REF] Nowadays, ITER( International Thermonuclear Experimental Reactor shown in figure 1.5) is the largest fusion device on earth and the most ambitious energy project in the world. ITER is located in southern France, 35 nations are collaborating to build the largest tokamak device in the south of France. ITER will be the first fusion device to produce a net energy. The aim for ITER is to produce high amplification factor Q (reach 10) for a long sustainable time(> 300 s). The amplification factor Q represents the ratio of output energy to input heating energy.

Heating by neutral beam injection

In ITER, the temperatures must reach more than million Celsius in order for the gas to achieve the plasma state and ignite the fusion reaction. The Tokamak relies on three sources of external heating: the ohmic heating [START_REF] Haines | The joule heating of a stable pinched plasma[END_REF], the wave heating and the particle heating.

In a tokamak device, the toroidal current flowing in the plasma induces poloidal magnetic field, and also heats the plasma by Joule effect (ohmic heating). The heating process strongly depends on the resistivity of plasma and the applied current. As the plasma temperature increases, the plasma resistivity decreases. Because of this effect, ohmic heating becomes very inefficient when applied to fusion device at CHAPTER 1. INTRODUCTION very high plasma temperature. So additional heating methods are required.

Another heating method is wave heating at the cyclotron resonance. The principle is that if the wave frequency is the same as the particle cyclotron frequency, the charged particle rotation is synchronous with the electric field rotation, i.e. the charged particle "sees" a constant field and its energy increases as in a dc field. Cyclotron resonance heating can be used for electrons and ions at different frequencies (since the cyclotron frequency depends on the charged particle mass). For magnetic fields in fusion reactors, ion cyclotron heating occurs for frequencies in the MHz range and electron cyclotron heating for frequencies of a few hundred GHz. The third way of plasma heating in fusion reactors is heating by neutral beam injection which consists in injecting a large flux of energetic neutral atoms in the tokamak to heat the charged particles through collisions. Neutral beam injection will play an essential role in the heating of ITER. The ITER plasma will be heated by two 1 M eV negative ion beams, corresponding to a total current of 40 A. The negative ions are produced in a negative ion source (a low temperature partially magnetized plasma), are extracted by a system of biased grids, accelerated at 1 MeV, neutralized by collisions with neutral atoms, and injected in the tokamak. A schematic of the neutral beam injection system of ITER is shown in Fig 1 .7. A high energy of 1 M eV is needed to heat the plasma at the required temperature. The reason for using negative ions instead of positive ions (much easier to generate and extract) is the fact that the neutralization efficiency of 1 M eV positive ions is very small while the neutralization of 1 M eV negative ions by collisional detachment is reasonably large as can be seen in In this thesis we focus on one specific aspect of the negative ion source, i.e. on the question of co-extraction of electrons with negative ions. The negative ion source of ITER is described below.

Negative ion source:

The structure of negative [START_REF] Speth | Overview of the rf source development programme at ipp garching[END_REF] ion source showed below is a prototype used in IPP Garching. Basically negative ion source has 3 parts: source area(driver), expansion area and extraction region. In source area (driver), the plasma is induced either by heated filaments(arc discharge) or rf antennas(inductively coupled discharge). Based on the disadvantages of big caesium consumption and short lifetime in the filament-ed sources, the ITER consortium has decided to adopt the rf driven source as the main source of neutral beam injectors (see in Fig 1.9). In the expansion region between the driver and the extraction region a magnetic field (magnetic filter around 100 Gauss) is placed perpendicular to the expanding plasma flow. The driver generates high energy electrons that can produce vibrationallly excited H 2 (D 2 ) and atomic hydrogen (deuterium). Hydrogen (deuterium) atoms hit the surface of the plasma grid (which is coated with a cesium layer) and are converted into negative ions that can be extracted by the extracting voltage. Negative ions are also produced in the plasma volume from dissociative attachment of vibrationally excited states by low energy electrons (dissociative attachment of H 2 (ν > 3) , where ν is the vibrational level, is the largest for electron energy on the order of 1 eV ). The high energy electrons generated in the driver can destroy negative ions. The purpose of the magnetic filter in front of the plasma and extracting grids is to lower the electron energy in the extraction region and therefore to limit the destruction of negative ions by high energy electrons. Electrons are supposed to be trapped by CHAPTER 1. INTRODUCTION the magnetic field in the filter and their residence time in the expansion region is therefore increased. They can therefore undergo collisions and lose energy before reaching the extraction region.

Since a positive voltage is applied to the extracting grid, electrons can also be extracted with negative ions. An important issue in negative ion sources is to limit the current of co-extracted electrons. This is done in the ITER negative ion source by a magnetic field (called suppression magnetic field) generated by magnets in front of each grid aperture. Negative ions are not very sensitive to this magnetic field and can be freely extracted from the plasma meniscus in front the grid aperture. Electrons are strongly magnetized in front of the grid aperture and should be in principle deflected to the plasma grid and not extracted. However experiments on the ITER negative ion sources show that the co-extracted electron current is not negligible and is on the order of the negative ion current. The purpose of this thesis in this context is to understand and quantify the co-extracted electron current. Various experimental studies have been conducted in these years. Figure 1.10 shows a few prototypes of negative ion source: BATMAN (bavarian test machine for negative ions) and MANITU (multi ampere negative ion test unit) [START_REF] Franzen | Progress of the development of the ipp rf negative ion source for the iter neutral beam system[END_REF] are test facilities with a small prototype source (0.3 × 0.6 m 2 base area); ELISE and RADI [START_REF] Heinemann | The negative ion source test facility elise[END_REF] are larger facilities (0.9 × 1 m 2 ) which have half the size of the ITER source(0.9 × 2 m 2 ); SPIDER AND MITICA [START_REF] Toigo | Progress in the realization of the prima neutral beam test facility[END_REF] are full size test facility at RFX in Padua. Different diagnostic methods are applied to study the parameters both in expansion volume and extraction region. A preferable way to measure the plasma propertities is non-invasive optical emission spectroscopy(OES) [START_REF] Fantz | A comparison of hydrogen and deuterium plasmas in the ipp prototype ion source for fusion[END_REF]; Also langmuir probes are used in the boundary layer for detecting the plasma parameters; Cavity ring-down spectroscopy(CRDS) is used to study the negative ion density [START_REF] Berger | Cavity ring-down spectroscopy on a high power rf driven source for negative hydrogen ions[END_REF][START_REF] Wimmer | Investigation of the boundary layer during the transition from volume to surface dominated h-production at the batman test facility[END_REF]; Tunable diode laser absorption spectroscopy(TDLAS) is applied to determine the neutral Cs density [START_REF] Wimmer | Cesium dynamics and h-density in the extended boundary layer of negative hydrogen ion sources for fusion[END_REF]. For beam study: tunsten wire calorimeter is used to assess the shape of the extracted beam and get a rough estimation of the beam divergence [START_REF] Nocentini | Beam diagnostic tools for the negative hydrogen ion source test facility elise[END_REF]. Beam emission spectroscopy(BES) is used to diagnose the intensity and divergence of the beam [START_REF] Nocentini | Beam diagnostic tools for the negative hydrogen ion source test facility elise[END_REF].

Modeling of negative ion sources

Three types of simulation methods have been used to study negative ion sources:

1. The global model or 0D model is based on the assumption that the plasma discharge can be fully characterized by volume, surface area, absorbed power,mass flow and gas components and by solving a set of differential equation of particles and energy conservation to quantitatively study the plasma [START_REF] Monahan | Global models of electronegative discharges: critical evaluation and practical recommendations[END_REF]. A general particle(mass) conservation equation can be expressed as:

dn k dt = F k V - S p V n k + N k i=1 a i K i N r i j=1 n i,j (1.13)
where n k is the particle density of the kth plasma species, F k is the flow rate of the species k into the chamber, S p is the chamber pump speed, N k is the total number of chemical reactions involving species k, a i is an integer presenting the number of k-particle produced/lost per reaction in the ith process, K i is the reaction rate for the ith process, N r i is the number of reactant involved in the ith reaction and n i,j is the number density of the jth reactant in the ith chemical process. The first two terms of right-hand-side(rhs) of equation 1.13 represent the flow of particles into and out of the chamber. The last term CHAPTER 1. INTRODUCTION accounts for particle production and loss through chemical reactions in volume or on surface.

An energy balance equation is also solved for the electrons:

d dt ( 3 2 n e k B T e ) = P V - N + i=1 A V ε ei,i h i u B,i n i - Ne i=1 ([ε th,i - 3 2 δ e i (k B T e -k B T g )]a i K i N i j=1 n i,j ) (1.14)
where k B T e and k B T g are the electron and gas temperature respectively, ε ei is the energy loss per electron-ion pair cross the sheath, N + is the number of positive ion species in the model, h i is the ith positive ion transportation parameter, u B,i is the Bohm speed of the ith ion, N e is the number of reactions that have electrons as reactants, ε th,i is the threshold energy associated with process i, and δ e i is the fraction of electrons' kinetic energy loss through the ith processes. Equation.1.14 has three terms on the rhs. The first one correspond to the rate of change of electron energy density due to the power absorption, the second one is the electron energy loss to the wall and the last term is the energy loss through electron heavy-particle collisions, respectively. By solving the N k + 1(N k is the total number of particle species in simulation) sets of differential equations, the global model can be a powerful tool to study the complex plasma chemistry process in negative ion source chambers. Capitelli et al [START_REF] Capitelli | Open problems in the physics of volume h/sup-//d/sup-/sources[END_REF] had addressed the production of vibrationally excited states of H @ and D 2 by atomic recombination on surface through different mechanisms. Zorat and al [START_REF] Zorat | Global model of a radiofrequency h 2 plasma in denise[END_REF][START_REF] Zorat | Global model for an rf hydrogen inductive plasma discharge in the deuterium negative ion source experiment including negative ions[END_REF] used a global model to study the effect of cross-section data on H -density and calculated the density of different species in the deuterium negative ion source experiment(DENISE).

The challenge for the Global model is the assumption of uniform charged particles. In the ITER negative ion source, this assumption is a poor one, because of the presence of the magnetic filter and of the non-uniform plasma density in the expansion region and extraction region.

2. A fluid model can be used to study the plasma in terms of particles densities n, mean velocity u and mean energy ε. The model solves fluid equations such as the continuity, momentum, and energy equations which are partially differential equations in time and space, quasi-neutrality is often assumed in fluid models so Poisson's equation is not solved.

CHAPTER 1. INTRODUCTION

Continuity equation:

∂n ∂t + ∇ • (n u) = S (1.15)
the source term S is the net number of particles created per unit time per unit volume in all types of reactions:

S = i N i n 1,i n 2,i k i (1.16)
where N i is the numbers of particles produced in ith collision, n 1,i and n 2,i are the incident particle densities of the ith collision and k i is the rate coefficient.

Momentum equation:

∂n u ∂t +∇•(n u u)+ 1 m ∇• P - qn m ( E + u× B) = -n j m j m + m j n j k m,j ( u-u j )
(1.17

)
where

P = m ( v -u)( v -u)f d 3 v (1.18)
is the pressure tensor and k m,j is the rate coefficient for momentum transfer to particles of other species j. The complex momentum equation is usually simplified by different approximation under different simulation conditions.

Energy equation:

∂enε ∂t + ∇ • (en uε + P • u + Q) -qn u • E = Π (1.19)
where Π is the net power density gained in reactions and Q is the heat flux vector defined as:

Q = 1 2 m | v -u| 2 ( v -u)f d 3 v (1.20)
By solving the fluid equations, JP Boeuf and al [START_REF] Boeuf | Model of an inductively coupled negative ion source: Ii. application to an iter type source[END_REF] used 2D fluid model to study the plasma parameters as function of power(10-80 kW ) and pressure(0.2-0.8 P a). A saturation of the plasma density growth with power of 40kW and pressure around 0.3Pa was observed. H 2 and H temperatures and densities spatial distribution were studied. By using a PIC-MCC model, F Taccogna and al [START_REF] Taccogna | Particle modelling of the hybrid negative ion source[END_REF] had studied the electron induced collision for the production of vibrational excited H 2 molecules and the effect of plasma potential on the transport of surface produced H -in the extraction region. G Fubiani and al [START_REF] Fubiani | Role of positive ions on the surface production of negative ions in a fusion plasma reactor type negative ion source-insights from a three dimensional particle-in-cell monte carlo collisions model[END_REF] used 3D self-consistent PIC model with Monte Carlo Collisions to calculate the plasma parameters of the BATMAN prototype in Garching. The detail of the production of H -on the PG surface results from the impact of H and H + . A Hatayama and al [START_REF] Hatayama | Particlein-cell modeling of negative-ion transport and extraction processes in a hydrogen negative-ion source[END_REF] have discussed the effect of the magnetic field on the H -extraction and the spatial structure of potential near the PG was studied. D Wunderlich and al [START_REF] Wünderlich | Pic code for the plasma sheath in large caesiated rf sources for negative hydrogen ions[END_REF] studied the influence of H + , H -, Cs + on plasma sheath with a 1D PIC MCC simulation code.

Magnetron discharges and ion thrusters

Magnetron discharges are used in plasma processing, e.g. in sputtering-deposition reactors. HiPIMS (High Power Impulse Magnetron Sputtering) have been the subject of many research works and publications in the last ten years. In these devices, electrons are strongly magnetized by a transverse magnetic field while the practically collisionless ions are accelerated in the cathode sheath toward the cathode and sputter atoms from the cathode surface.

The magnetic configuration of Hall thrusters is very similar to that of magnetrons discharges but the purpose here is not to accelerate ions toward a surface but to eject ions at high speed away from the thruster. An electron emissive cathode is used in order to allow ion acceleration and extraction through a quasineutral plasma and to neutralize the ion beam space charge.

The large electron drift and negligible ion drift (ions are practically not magnetized) in the azimuthal direction of E × B devices tend to create a charge separation and trigger instabilities. Instabilities in a variety of wavelengths and frequencies ro-CHAPTER 1. INTRODUCTION tating in the azimuthal direction have been observed in these devices. In this thesis we focus on one type of instability, qualified as "rotating spoke" that is observed in both magnetron discharges and Hall thrusters. More details are given below. In Chapter 5 of this thesis, we use a PIC MCC model to better understand the physics of this instability.

Magnetron discharges

Magnetron discharges are cylindrical devices with radial magnetic field and axial electric field, or axial magnetic field and radial electric field. We consider here only planar magnetron discharges, with radial magnetic field and axial electric field (see Figure 1.11) In planar magnetrons, a magnetic barrier is formed in front of the cathode by ring magnets. The radial magnetic field is maximum on the cathode surface and decreases away from the cathode. Electrons emitted by the cathode under ion bombardment are trapped along the magnetic lines and can in principle cross the magnetic barrier only under the effect of collisions. In practice, when the Hall parameter (ratio of the electron cyclotron frequency and electron collision frequency) is sufficiently large, collisional transport across the magnetic field is no longer dominant, and instabilities can contribute to cross-field electron transport.

One of the instabilities that seem to play an important role in electron transport in magnetron discharges is called a Rotating spokes instability. Rotating Spokes are CHAPTER 1. INTRODUCTION rotating region of enhanced light emission that are commonly observed in planar magnetron discharges and a number of papers have been recently devoted to this phenomenon. Although there is no completely clear understanding of rotating spokes in magnetron discharges, some progress have been made in this direction thanks to the use of detailed space and time resolved optical and electrical diagnostics. The groups of A. Anders in California and of A. von Keudell in Germany have played an important role in the characterization of rotating spokes in high power magnetron discharges as well as in dc magnetron discharges. Reviews articles by these two groups have been recently published (A Hecimovic and A von Keudell [START_REF] Hecimovic | Spokes in high power impulse magnetron sputtering plasmas[END_REF]; A Anders [START_REF] Anders | Tutorial: Reactive high power impulse magnetron sputtering (rhipims)[END_REF]). Spoke rotation has been observed in the +E × B direction as well as in the -E×B direction. The -E×B rotation is generally seen in low power dc magnetron discharges or at the beginning of the current rise in HiPIMS. During the high power regime of HiPIMS, the direction of spoke rotation changes to +E × B. .13 shows the space distributions of electric potential and light emission at a given time, for a spoke rotating in the -E × B direction, in a low current, dc regime (270 V, argon, 0.27 Pa). We see that a rather large potential drop is present at the spoke's front. The authors mention that the ionization enhancement at the spoke's front is associated with a double layer and a potential "hump" but there is no clear understanding of this structure and of the fact that spokes can rotate either in the E × B direction or in the -E × B direction.

Hall thrusters and HEMPT

Hall thrusters (Fig 1 .14) are used for satellite propulsion. They are able to produce a thrust on the order of 70 mN/kW , with ion velocity around 20 km/s. Xenon is generally used because it is more efficient than other rare gases for propulsion (larger mass and easier to ionize). Hall thrusters are E × B devices with axial electric field and radial magnetic field. The magnetic field configuration of Hall thrusters is similar to that of planar magnetrons, as shown in Fig 1 .14, except that the cathode is not in the magnetic barrier. The magnetic field is maximum in the exhaust plane at the end of a channel between two concentric dielectric cylinders where the plasma is generated. Electrons generated by an external emissive cathode enter the channel across the magnetic field barrier under the effect of collisions and instabilities, and ionize the gas injected through the anode. A large axial electric field forms in the magnetic barrier because of the lowering of the electron conductivity by the radial magnetic field in this region. This electric field forms in the quasineutral plasma (not in a cathode sheath) and can accelerate ions out of the plasma without grid (Hall thrusters are gridless ion thrusters). The transport of electrons across the magnetic field seems to be different on each side of the magnetic field maximum. On the cathode side (increasing magnetic field from cathode to exhaust plane) it has been shown that an azimuthal instability called Electron Cyclotron Drift Instability, ECDI, contributes to electron transport across the magnetic field (J. C. Adam, A. Heron, and G. Laval [START_REF] Adam | Study of stationary plasma thrusters using two-dimensional fully kinetic simulations[END_REF]). In this region the gas density is very low, because a large part of the gas flow is ionized upstream, and electron transport seems to be dominated by this instability as well as by electron-wall interaction (electron scattering by the channel walls and secondary electron emission by electron impact on the walls can contribute to cross-field electron transport). Collective light scattering has been used to analyze the instability(S.Tsikata et al [START_REF] Tsikata | Collective light scattering for the study of fluctuations in magnetized plasmas: The hall thruster case study[END_REF]), and there is experimental evidence that the ECDI is present in Hall thrusters(J.Cavalier et al [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the e× b electron drift instability: Modeling and fitting on experimental data[END_REF]). The ECDI is a microinstability (wavelength on the order of electron Larmor radius) that seems to evolve into an ion acoustic instability (wavelength of a few Debye length) (T. Lafleur, S. D. Baalrud, and P. Chabert [START_REF] Lafleur | Theory for the anomalous electron transport in hall effect thrusters. ii. kinetic model[END_REF], J. P. Boeuf and L. Garrigues [START_REF] Boeuf | E× b electron drift instability in hall thrusters: Particle-in-cell simulations vs. theory[END_REF]). On the anode side CHAPTER 1. INTRODUCTION of the maximum magnetic field in the exhaust plane, another type of instability is present is the magnetic field in the anode region is sufficiently large (M. McDonald and A. D. Gallimore [START_REF] Mcdonald | Rotating spoke instabilities in hall thrusters[END_REF], C. L. Ellison, Y. Raitses, and N. J. Fisch [START_REF] Ellison | Cross-field electron transport induced by a rotating spoke in a cylindrical hall thruster[END_REF], M. J. Sekerak, B. W. Longmier, A. D. Gallimore, D. Brown, R. R. Hofer, J. E. Polk [START_REF] Sekerak | Azimuthal spoke propagation in hall effect thrusters[END_REF] and S.Mazouffre [START_REF] Mazouffre | Rotating spoke instabilities in a wall-less hall thruster: experiments[END_REF]). This instability is macroscopic and its wavelength is not small with respect to the circumference of the channel. This instability appears as a region (or several regions) of enhanced luminosity rotating in the azimuthal direction, and the term "rotating spoke" is generally used to define it. Indeed, the properties of this instability are similar to those of the rotating spokes observed in magnetron discharges. HEMPT (Highly Efficient Multistage Plasma Thrusters) or mutli-cusped field thrusters (see Fig 1.16) are another type of gridless ion thrusters. The plasma is formed in a cylinder with dielectric walls. As in Hall thrusters, an external emissive cathode is used to ionize the neutral flow coming from the anode at the end of the cylinder, and to neutralize the extracted ion beam. Instead of a magnetic barrier, a cusped magnetic field is used to lower the axial electron conductivity, leading to a potential drop around the exhaust plane of the cylinder (in principles, several cusp can be used, with a potential drop between two consecutive cusps).

CHAPTER 1. INTRODUCTION

Magnetic cusps are often used in plasma applications not to lower the electron conductivity parallel to the discharge axis, but to limit charged particle losses to the walls in plasmas sources as discussed in the next section and in Chapter 3 of this thesis.

Plasma confinement by magnetic cusps

Magnetic cusp are used in many plasma sources to confine the plasma by limiting the particles losses to the walls. Magnetic cusps were initially investigated as a way to confine particles in thermonuclear fusion plasma [START_REF] Haines | Plasma containment in cusp-shaped magnetic fields[END_REF][START_REF] Spalding | Advances in plasma physics[END_REF]. The basic idea of magnetic cusps is to have a magnetic field parallel to the walls over large regions of the vacuum chamber. This can be easily done for example by placing magnets along the wall surface. In that case, the magnetic field lines are parallel to the walls over of the length of the magnet and are directed toward the walls in the pole regions. The region where the magnetic field goes to the walls are called cusps. Since transport of charged particles across the magnetic field is strongly reduced, they can reach the walls only by following the magnetic field lies converging to the cusps. The charged particle loss area is therefore reduced when magnetic cusps are present, and one can define an effective loss area or loss width associated with the cusps. Figure 1.17 shows different possible configurations of magnetic cusps. Depending on the way the magnetic field is generated, the magnetic field lines can converge to a line ("line cusp"), to a ring ("ring cusp"), or to a point ("point cusp"). Plasma confinement by cusps has been studied since the 1970's when Limpaecher and MacKenzie [START_REF] Limpaecher | Magnetic multipole containment of large uniform collisionless quiescent plasmas[END_REF] showed that they could increase the plasma density in a low pressure discharge by two orders of magnitude by employing multipolar cusp confinement. A number of experimental as well as theoretical papers were published in the 1970-1980's on magnetic cusps. They were initially studied as possible confinement methods for thermonuclear fusion (M.G. Haines [START_REF] Haines | Plasma containment in cusp-shaped magnetic fields[END_REF], I. Spalding [START_REF] Spalding | Advances in plasma physics[END_REF]), but they were actually more used for basic plasma studies. Magnetic cusps are now commonly used in different types of ion sources or in plasma processing reactors, but, in spite of the important efforts to quantify the role of magnetic cusps on plasma confinement, there is still no clear scaling laws that give the effective charged particle loss area or loss width in the presence of cusps, as a function of plasma parameters such as electron temperature, gas pressure, magnetic field and ion mass. In their 1975 paper, Hershkowitz et al. (Hershkowitz et al [49]) mention (and this is still true today) that "Although the motion of charged particles in magnetic cusps has received considerable attention, a satisfactory description of how plasma leaks through cusps

CHAPTER 1. INTRODUCTION
has not yet emerged. In particular, there seems to be considerable disagreement between theory and experiments". We briefly summarize below the experimental and analytical or modeling efforts that have been devoted to the characterization of cusps. We are considering here only low-β plasmas, i.e. plasmas where the plasma pressure is much less than the magnetic pressure. In these plasmas the external magnetic field is not modified by the charged particles currents. In the following, the leak width w L is defined as the half-with of the profile of the current density to the wall. For a metallic surface the electron and ion current density profiles are not necessarily identical and one could define specific leak widths for electrons and ions. The electron and ion leak width should however become close to each other when the plasma density increases (because the plasma becomes quasineutral in a larger volume and close to the surface). For a dielectric surface one expects the electron and ion loss widths to be identical.

Most theories for high-β plasmas predicted a loss aperture or leak width w L between the electron gyroradius ρ e and the hybrid gyroradius defined as (ρ e ρ i ) 1/2 (geometric mean of the electron and ion gyroradii ρ e and ρ i ) while experiments found different results lying between the ion gyroradius, the hybrid gyroradius and the electron gyroradius (see references in Hershkowitz et al. [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]). The first measurements of the leak width in low-β plasmas where performed in 1975 by Hershowitz et al. [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]. They found that the electron leak width defined as the full width at half maximum (FWHM) was about four times the hybrid gyroradius:

w L,H = 4(ρ e ρ i ) 1/2 (1.21)
for low enough gas pressure and that the ion leak width was slightly larger. The measurements were performed in a plasma source divided in two regions, the driver and the target, separated by a "picket fence" structure consisting of an array of parallel conductor wires with currents in adjacent conductors flowing in alternating directions (see Fig 1.18). The picket fence magnetic cusps were generated by the currents flowing in this conductor array. The plasma was sustained by high energy electrons emitted from hot filaments in the driver. The plasma and high energy electrons were confined by the magnetic field of permanent magnets arranged in a full line-cusp geometry. The leak width was deduced from measurements of the electron and ion current through the picket fence (i.e. between the wires of the picket fence). The measurements were performed in helium, argon, and xenon, and confirmed the M 1/4 dependence of the leak width present in the expression of the hybrid gyroradius. The 1/B dependence of this expression was also checked. Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF] performed similar experiments in ring and point cusps over a large range of magnetic field strengths (between 10 and 260 G) and neutral gas pressures (between 10 -2 and 10 mtorr). The magnetic field was generated by two water-cooled coil of 17 cm inner diameter that produced a spindle cusp magnetic CHAPTER 1. INTRODUCTION field using currents of up to 1000 A. The maximum magnetic field in the center of the ring cusp was 160 G, while the maximum field in the center of the point cusp was 260 G. The current profiles in the cusps were measured by Langmuir probe. Figure 1.19 shows the measured profiles of ion saturation current, electron saturation current, and primary electron current in a 110 G ring cusp at 5.5 × 10 -5 torr in argon. We see on this figure that the current profiles present a sharp maximum in the cusp center, and that the ion leak width is larger than the electron leak width. The measurements of Bosch and Merlino also confirmed that the electron and ion leak widths become closer to each other when the plasma density is increased, as can be seen in Fig 1 . The measurements showed that the variations of the leaked widths with gas pressure was close to p 1/2 and that the leak widths varied as B -1 at relatively low magnetic field, and as B -1/2 at higher magnetic fields. This is illustrated in Fig 1 .21. The breakdown of the B -1 dependence of the leak width at high magnetic field and low pressure was attributed to the fact that collisional, ∝ B -2 diffusion was replaced by anomalous or Bohm, ∝ B -1 diffusion in these conditions. They proposed a semi-empirical theory where the leak width was proportional to the square-root of the cross-field electron diffusion, leading to the B -1 leak width at low magnetic field and B -1/2 at higher magnetic fields. The analytical expression of the leak width proposed by Bosch and Merlino was:

w L,B = (2 DR/C s ) 1/2 (1.22)
where D is an average effective ambipolar diffusion coefficient perpendicular to the magnetic field, R is related to the cusp length, and C s is the ion acoustic velocity.

It is interesting to note that although the expressions w L,H and w L,B of papers by Herskowitz et al. and by Bosch and Merlino are very different, they both scale as B -1 M 1/4 where B is the magnetic field and M the ion mass. The expression w L,B of Bosch and Merlino is certainly more realistic or more useful since it takes into account the effect of collisions which naturally appears in the diffusion coefficient. In that case, D is proportional to the electron-neutral collision frequency ν eN ( D ≈ KT e ν eN /(mω 2 ce )) i.e. the leak width scales as ν 1/2 eN or as P 1/2 where P is the gas

CHAPTER 1. INTRODUCTION
pressure. The scaling in electron temperature will be discussed in chapter 3.

Finally, we mention that Koch and Mathieussent [START_REF] Koch | Collisional diffusion of a plasma in multipolar and picket fence devices[END_REF] presented a theoretical derivation of the leak width that led to the following expression:

w L,K = 2R π ( ρ e ρ i λ e λ i ) 1/2 (1.23)
The scaling of this expression with magnetic field and ion mass is similar to that of Hershkowitz et al. and Bosch and Merlino, but w L,K is proportional to the pressure P and not to P 1/2 , which is in disagreement with the measurements of Bosch and Merlino.

The question of the quantitative description of the confinement by magnetic cusps is still an issue since, in a recent paper, Cooper et al [START_REF] Cooper | Direct measurement of the plasma loss width in an optimized, high ionization fraction, magnetic multi-dipole ring cusp[END_REF]. explored experimentally and with a simplified model (similar to the model of Bosch and Merlino), the loss width of their large multipole magnetic ring cusp WiPAL used in the Madison plasma dynamo experiment.

The available computing power is now becoming sufficient to explore quantitatively and more precisely the confinement by magnetic cusps in low-β plasmas. To verify the scaling of the leak width with the different parameters, we defined a simplified 2D Particle-In-Cell Monte Carlo (PIC MCC) model. This is described in Chapter 3. Another effort to solve this problem in the frame of fluid models is also on-going in our laboratory, in collaboration with the Institute of Mathematics of Toulouse.

Chapter 2

Particle simulations of magnetized plasmas

Physics is the belief that a simple and consistent description of nature is possible.

Niels Bohr

Numerical model

In this chapter, we review the basic principles of the Particle-In-Cell (PIC) method with Monte-Carlo Collisions (MCC). First we will describe the algorithm in details and calculate in the 1D limit the properties of a plasma in the ambipolar regime. We will compare the latter to an analytical description of the plasma. Next, we will discuss the preservation of the adiabatic invariants using the Boris [START_REF] Boris | Relativistic plasma simulation -optimisation of a hybrid code[END_REF] numerical solver to calculate the particle trajectories in the case of magnetised electrons trapped in a magnetic bottle with mirrors at both ends [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Nicholson | Introduction to plasma theory[END_REF][START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF]. This analysis is relevant to the work of this thesis as we study plasma confinement in a magnetic cusp and related instabilities.

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS

Decoupling the terms in the Boltzmann equation

The PIC-MCC technique solves the Boltzmann equation,

∂f i ∂t + v • ∂f i ∂x + a • ∂f i ∂v = ∂f i ∂t c , (2.1)
in two distinct steps [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF][START_REF] Nanbu | Direct simulation scheme derived from the boltzmann equation. i. monocomponent gases[END_REF], that is, (i) one calculates the particle trajectories in the force field and (ii) the collision processes between particles are evaluated.

∂f i ∂t c = t (f i f t -f i f t ) v r σ T t dΩdv t , (2.2)
is the collision operator, f i (f t ) is the distribution function for the incident (target) specie, respectively, m i the mass, F = m i a the force field, a(r, v, t) the acceleration, r (v) the position (velocity) vector, v r = |v i -v t | the relative velocity, σ T t (v r ) the total differential cross-section (summed over all the collision processes between the incident and the target particles) and, lastly, Ω the solid angle. r, v and t are independant variables. Primes denote the distribution function after the collision. One as follows, rewrites the eq. (2.1) as

∂f i ∂t = Df i + Jf i , (2.3) with D(f i ) = -v • ∂f i ∂x -a • ∂f i ∂v , (2.4)
and

J(f i ) = ∂f i ∂t c , (2.5)
respectively. Expanding ∂f i /∂t for small ∆t around t = t 0 [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF][START_REF] Nanbu | Direct simulation scheme derived from the boltzmann equation. i. monocomponent gases[END_REF], one gets

f i (x, t 0 + ∆t) = f i (x, t 0 ) + ∆t ∂f i ∂t t=t 0 , (2.6)
which is equivalent to

f i (x, t 0 + ∆t) = f i (x, t 0 ) + ∆tDf i (x, t 0 ) + ∆tJf i (x, t 0 ) , = (1 + ∆tD + ∆tJ) f i (x, t 0 ) , (2.7)
where eq. (2.3) evaluated at t = t 0 was inserted into eq. (2.6). For sufficiently small time steps, eq. (2.7) may be recasted in a form which is second order correct in ∆t,

f i (x, t 0 + ∆t) = (1 + ∆tJ) (1 + ∆tD) f i (x, t 0 ) . (2.8)
CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS One can hence calculate f i (x, t 0 + ∆t) in two steps [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF][START_REF] Nanbu | Direct simulation scheme derived from the boltzmann equation. i. monocomponent gases[END_REF]; first apply the operator (1 + ∆tD) to f i (x, t 0 )

f * i (x, t 0 + ∆t) = (1 + ∆tD) f i (x, t 0 ) , (2.9)
being equivalent to solving the Vlasov equation,

∂f * i ∂t + v • ∂f * i ∂x + a • ∂f * i ∂v = 0 , (2.10)
then the collisions processes are implemented via the operator (1 + ∆tJ) applied to the particle distribution function f * i (x, t 0 ) after the evaluation of the particle trajectories at t = t 0 ,

f i (x, t 0 + ∆t) = 1 + ∆t t f * i f * t -f * i f * t v r σ T t dΩdv t , (2.11) 2.1.2 Particle-in-Cell technique 2.1.2.

Equations of motion

The principles of the Particle-In-Cell (PIC) method are described in textbooks [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF]. The technique is equivalent to a characteristic solution of Vlasov equation, eq. (2.10). Characteristic solutions are found by interpreting the Vlasov equation as the total derivative of f i along a particle orbit defined as X(t) and V(t) where X and V are the Lagrangian (time dependant) representation of the particle motion. One finds

d dt f i [X(t), V(t), t] = ∂f i ∂t + 3 k=1 Ẋk ∂f i ∂X k + 3 k=1 Vk ∂f i ∂V k = 0 , (2.12)
where Ẋ(t) = V(t) and V(t) = a(t) = F i /m i taken along a particle orbit, respectively. F i is the force at the particle location. The PIC technique maps f i using a large number of so-called macroparticles which represents a large number of real particles of a given specie. The mass and charge of the macroparticle is hence larger than the one of a single particle, m i = N m M i and q i = N m Q i where N m is the macropaticle weight. The charge over mass ratio of the macroparticle is identical to the one of the real particle and so is the macroparticle trajectory.

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS

In magnetised plasmas, the force felt by the particles is the Lorentz force,

F i [X(t), t] = q i {E [X(t), t] + V(t) × B [X(t), t]} , (2.13)
with E [X(t), t] the electrostatic and B [X(t), t] the magnetic field both evaluated at the macroparticle location. In cold plasmas the self electromagnetic fields generated by particle currents can in general be neglected. This will be the case for the examples studied in this thesis. The magnetic field is prescribed while the electric field is induced solely by the charge separation. The Newton's equation of motion are solved numerically using the Boris scheme [START_REF] Boris | Relativistic plasma simulation -optimisation of a hybrid code[END_REF],

v n+1/2 -v n-1/2 ∆t = q i m i E n (x n ) + v n+1/2 + v n-1/2 2 × B (x n ) . (2.14)
The latter preserves the phase space volume but is not symplectic [START_REF] Qin | Why is boris algorithm so good?[END_REF][START_REF] Ellison | Comment on "symplectic integration of magnetic systems": A proof that the boris algorithm is not variational[END_REF]. In addition, this numerical technique does not need any intermediate storage of observables as required by non-symplectic integrators such as Runge-Kutta for instance [START_REF] Press | Metcalf[END_REF]. This greatly reduces the computer memory usage of the algorithm. Equation (2.14) is solved numerically as follows; first define intermediate velocity vectors v -and v + such as [START_REF] Birdsall | Plasma physics via computer simulation[END_REF][START_REF] Hockney | Computer Simulation Using Particles[END_REF],

v n-1/2 = v -- q i E n m i ∆t 2 , (2.15)
and,

v + = v n+1/2 - q i E n m i ∆t 2 , (2.16)
respectively. Inserting eqs. (2.15) and (2.16) into (2.14) yields,

v + -v - ∆t = q i m i v + + v - 2 × B , (2.17)
which is a rotation of v around the magnetic field vector during a time step ∆t. This rotation may be decomposed into two intermediate steps,

v = v -+ v -× t . (2.18)
and,

v + = v -+ v × s , (2.19)
where

t = tan θ 2 B |B| , (2.20)
s = 2t 1 + t 2 , (2.21)
CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS and θ = ω c ∆t which is the rotation angle during ∆t (ω c = q i |B|/m i ). Equations (2.15) through (2.21) characterise the standard Boris algorithm but in a form which preserves both the amplitude and phase of the gyromotion [START_REF] Zenitani | On the boris solver in particle-in-cell simulation[END_REF]. In practice one assumes that tan(θ/2) θ/2 which introduces a phase error 2 tan (θ/2) /θ 1 -θ 2 /12 + • • • . An angle θ < 0.35 provides an error less than 1%. The latter is cumulative over the time steps but we will see below that this has no consequences on the phase-space volume-preserving properties of the algorithm. In all cases the Larmor radius is correctly calculated. of an electron with a kinetic energy of 2 eV gyrating around a magnetic field vector of strength 100G. The x-axis is normalised to the period of rotation t/τ where ω c τ = 2π while the y-axis to the Larmor radius x/r L , respectively. The time step is ω c ∆t = 0.3 (i.e., ∆t/τ 5 × 10 -2 ). We compare the Boris algorithm with the exact phase (red solid line), eq. (2.20), to the approximated solution (blue solid line) with, tan(θ/2) ω c ∆t/2, and the analytical solution (orange dots),

x = r L cos[ω c (t -t 0 )],
for reference. One can see the phase shift induced by the approximation of t in eq. (2.20). The amplitudes are identical for the three cases. . Schematic representation of a magnetic bottle (taken from [START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF]). Trapped electrons will initially follow a given field line and rotate slowly around the magnetic bottle describing a surface where the field line number density remains constant.

Trajectory of an electron confined in a magnetic mirror calculated with the Boris algorithm

In this section, we show that the Boris algorithm calculates correctly and over a long integration time the trajectory of an electron confined by a magnetic bottle with mirrors at both ends [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Nicholson | Introduction to plasma theory[END_REF][START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF]. A schematic description of a magnetic mirror is shown in fig. 2.2. Such a field may be generated by a pair of coils with the current circulating in the same direction. The equation for the magnetic field vector components is

B x = -5xB 0 z L 2 , (2.22) B y = -5yB 0 z L 2 , (2.23) B z = B 0 1 + 5 z 2 L 2 , (2.24)
where B 0 = 100G and L = 60 cm in our case (i.e., a bottle of 1.2 m length in total). The magnetic field is divergence free,

∇ • B = 0 . (2.25)
This field is a second order approximation which is valid only close to the axis of the magnetic bottle, i.e., for x/L 1 and y/L 1. In this case, the curl of B is approximatively null, ∇ × B 0 .

(2.26)

The motion of a charged particle in a magnetic mirror may be represented by three adiabatic invariants [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF][START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF][START_REF] Nicholson | Introduction to plasma theory[END_REF]. The first one is the magnetic moment, µ = m i v 2 ⊥ /2B which may be combined with the total particle energy, being conserved, to deduce CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS the minimum angle for the particle velocity vector in order to trap the latter in the bottle. One writes,

v 2 ⊥0 B 0 = v 2 ⊥ B m , (2.27)
which can be combined with the conservation of energy

v 2 ⊥ = v 2 ⊥0 + v 2 0 = v 2 0 , i.e., sin (θ c ) = v ⊥0 v 0 = B B m . (2.28)
B m is the maximum field at the exit of the bottle where the particle can be reflected and hence v = 0 at that location. v ⊥0 and v 0 are the particle velocity components at a given location along the trajectory inside the bottle (associated with a magnetic field strength B). Any particle with an angle θ < θ c will hence not be trapped. The angle θ c does not depend on the charge and mass of the particle; ions and electrons will consequently be equally well confined. For the conditions of our simulations, we find for instance R m = B m /B 0 = 6 at the center of the bottle corresponding to θ c 24 • . Trapped electrons are bouncing between mirror points following the field lines (this periodic motion is the second adiabatic invariant). The off-axis field lines have a curvature which induces a centrifugal force. In addition the electron experience a magnetic field gradient. Both result in a drift in the same direction which lies in the plane perpendicular to the magnetic bottle axis (the z-axis in our case). The velocity vector of the combined drifts may be written as follows [START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF]

v GC = - m i q i B 4 v 2 + 1 2 v 2 ⊥ ∇ 1 2 B 2 × B . (2.29)
This drift drives the particle trajectory away from a given field line. The particle slowly rotates around the magnetic bottle (this is the third adiabatic invariant) following a surface which encompases the field line associated with the initial conditions of the particle trajectory. Such a surface is shown in fig. 2.3 where we have used the Boris algorithm to calculate the particle trajectory. The latter reproduces precisely the particle motion. Note that one of the interesting fact we can deduce from the second adiabatic invariant, which determines the length of a given field line between two mirror points, is that the particle always returns to the field line corresponding to the initial condition (this is also true for the earth magnetosphere where the magnetic field is not symmetric, but slowly varying spatially, due to the solar wind).

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS

Error introduced by the discretisation of the magnetic field map

In practice, the magnetic and electric field maps are discretised on a grid in a PIC algorithm. In 2D for instance, the grid typically consists of n x + 1 (n y + 1) points defined in our case from i = 1 to i = n x + 1 (j = 1 to j = n y + 1), which are the location of the boundary conditions, with a distance between grid points of ∆x (∆y), respectively. Particles move inside the simulation domain due to the interaction with the force field and may be scattered by collisions. The electromagnetic fields must be interpolated at the particle position in order to estimate the Lorentz force and hence calculate the particle trajectory. The interpolation is linear in our model. The evaluation of the y-component of the magnetic field for instance at the particle location is

B y (x p , y p ) = p x p y By (i, j) + (1 -p x )p y By (i + 1, j)+ = p x (1 -p y ) By (i, j + 1) + (1 -p x )(1 -p y ) By (i + 1, j + 1) ,(2.30)
where By is the magnetic field array corresponding to the discretised values on the mesh nodes, x p and y p are the particle position and p x , p y , the fraction of the distance CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS between the particle and the next mesh node,

p x = X i+1 -x p ∆x , (2.31)
p y = Y j+1 -y p ∆y , (2.32)
and

1 -p x = x p -X i ∆x , (2.33) 1 -p y = y p -Y j ∆y , (2.34)
respectively. X i , Y j are the physical position of the mesh nodes and the interpolation is hence bilinear. Generalization to a 3D grid is straightforward. Next we evaluate the error introduced by the discretisation of the magnetic field map for the calculation of the particle trajectory confined inside the magnetic bottle described in the previous section. Figure 2.4 shows the exact solution for the z-component of the magnetic field versus the discretised values on the mesh nodes. The width between nodes is ∆x = ∆y = ∆z = 5 mm, providing a precise piecewise fit of the actual field map. The particle trajectory, corresponding to the same initial conditions as in fig. 2.3, is shown in fig. 2.5. The two are nearly identical, confirming the fact that one can calculate accurately the characteristics of a magnetic mirror using a discretised field map as long as the size of the mesh is much smaller than

L i = (∂ i B i /B i ) -1
where i stands either of x, y or z. Lastly, fig. 2.6 plots the axis position of the electron, which is bouncing back and forth between mirror points, versus the time step in the Boris algorithm.

Hybrid OpenMP and MPI parallelisation of the particle pusher

The parallelization is performed in an hybrid manner using OpenMP [65] and "Message-Passing-Interface" MPI libraries. We use a particle-decomposition scheme for the particle pusher where each core (thread) have access to the whole simulation domain (as opposed to a domain-decomposition approach). The number of particles per core is nearly identical. We further implemented a sorting algorithm [START_REF] Bowers | Accelerating a particle-in-cell simulation using a hybrid counting sort[END_REF] in order to limit the access to the computer memory (RAM) and boost the execution time, ∆t push , of the pusher subroutine. In the estimates below (figure 2.7), the pusher includes electron heating (absorption of power from an external source), field interpolations, update of the velocities and positions together with the charge deposition on the grid nodes. Particles are sorted per grid cell. The field and density arrays are hence accessed sequentially. ∆t push is shown in figure 2.7 normalized to the number of particles in the simulation. The best performance is obtained by attaching one MPI thread per socket and a number of OpenMP threads identical to the number of cores per socket. For the simulations of Fig. 2.7, we set the number of OpenMP threads to 10. We sort particles every 10 time steps without any loss of performance.

The calculation is performed with the 3D version of our PIC-MCC model and the numerical resolution is either 96 × 64 × 128 grid nodes or eight times larger with 80 particles-per-cell (ppc). The time gained in the pusher with the particle sorting is a factor ∼ 4. The sorting algorithm remains efficient as long as there is on average at least one particle per cell per thread. Beyond this limit ∆t push converges toward the value without sorting as shown in Fig. 2.7. We define the efficiency of the pusher without sorting as,

β = ∆t (1) push ∆t push N core , (2.35)
where N core is the number of cores (threads) and ∆t

(1)
push the execution time of the pusher for N core = 1. β should be equal to 1 for a perfect parallelization of the pusher. We find β 78% for 20 cores, 70% for 320 cores and lastly, dropping to ∼ 60% for 640 cores (i.e., about 23% loss in efficiency with respect to 20 cores). 6.002 6.003 6.004 6.005 6.006 6.007 6.008 6.009 6.01

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS

GHz).

There is 2 sockets per CPU, 20 cores in total.

Parallelization of the Poisson solver

We have developed a multi-grid solver for Poisson's equation,

∇ 2 φ = - ρ ε 0 , (2.36)
which is parallelised in a hybrid manner using OpenMP and MPI libraries.

ρ = N i q i n i , (2.37)
where N is the total number of particle species in the simulation. φ is the plasma potentiel and ρ is Poisson's equation source term which is derived from the particle positions at every time steps. n i is the density and q i the charge. In 2D, the plasma density is calculated by bilinear extrapolation of the particle positions on the grid

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS nodes n i,j = Np k=1 p x p y δn , (2.38) n i+1,j = Np k=1 (1 -p x )p y δn , (2.39
)

n i,j+1 = Np k=1 p x (1 -p y )δn , (2.40
)

n i+1,j+1 = Np k=1 (1 -p x )(1 -p y )δn , (2.41)
where δn = N m /(∆x∆y) is the contribution to the density of a single macroparticle, N m is the macroparticle weight, p x , p y are defined by eqs. (2.31)-(2.34) and N p is the number of macroparticles of a given specie, respectively. i and j correspond to the nearest grid location of the macroparticle

i = x p /∆x + 1 , (2.42)
j = y p /∆y + 1; . (2.43)
We discretise Poisson's equation on the grid nodes using the finite difference method [START_REF] Rosenberg | Methods for the Numerical Solution of Partial Differential Equations[END_REF]; in 2D

∂ 2 ∂x 2 + ∂ 2 ∂y 2 φ = - ρ ε 0 , (2.44) becomes a e φ i-1,j + a w φ i+1,j + a n φ i,j+1 + a s φ i,j-1 + a c φ i,j = b i,j , (2.45)
where

a e = a w = ε 0 /∆x 2 , (2.46) a s = a n = ε 0 /∆y 2 , (2.47) a c = -(a n + a s + a w + a e ) , (2.48)
and b i,j = -ρ i,j , respectively. Equation (2.45) is second order correct in ∆x and ∆y and may be written in matrix form

A φ = b , (2.49)
where φ and b are now vectors. Defining the residual at iteration k as

r = b -A φ , (2.50)
CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS which is simply the amount by which the approximation φ differs from the exact solution φ (i.e., r = 0 when φ = φ) and defining in turn the error as ê = φφ , we find, Aê = r .

(2.51)

We solve Poisson's equation iteratively on the grid nodes with a multi-grid solver [START_REF] Briggs | A Multigrid Tutorial[END_REF].

A relaxation method (so-called Successive-Over-Relaxation, SOR, in our case) is applied successively on different grid levels (from fine to coarse grid levels and viceversa). The SOR algorithm has the higher frequency components of its error vector dampened typically over a few iterations on a given grid. These frequencies are bounded by the grid size and hence lower frequencies (whose wavelength λ > ∆x and ∆y) will be removed on a much slower scale. SOR algorithms hence require a large number of iterations to converge, typically O(N 3/2 ) when optimized. Note that for the case of Gauss-Seidel and Jacobi methods, convergence is reached in O(N 2 ) iterations. One efficient mean to accelerate the convergence of the error ê and hence approaching in a lower number of iterations the solution φ, is to iterate the SOR algorithm a couple times on a given grid and then to perform the same step on a new coarser grid where ∆x and ∆y are twice as large, and so on. The extremely effective multigrid idea is hence to change to a coarser grid, on which "smooth becomes rough" and low frequencies act like higher frequencies. A multigrid cycle is as follow: one goes from the initial grid toward coarser grids until the grid is just a few nodes in each directions and then extrapolate the solution back toward finer grids up to the initial grid. This is so-called the V cycle, which is just one of many possibilities [START_REF] Briggs | A Multigrid Tutorial[END_REF]. The technique converges in O(N ) iterations.

We parallelized our solver via a domain-decomposition approach. The latter is decomposed in rectangular parallelepipeds along (Oz) where each of them is assigned an MPI thread. The do-loops are parallelized using the OpenMP framework. The SOR algorithm is semi-implicit in nature and hence its parallelisation is performed via the so-called "red-black" ordering. The algorithm is solved twice per time step, the first pass is

φ k+1 i,j = (1 -ω) φ k i,j + ω b i,j -a w φ k i-1,j -a s φ k i,j-1 -a n φ k i,j+1 -a e φ k i+1,j /a c , (2.52)
for all i + j odd where 0 < ω ≤ 2 is the SOR relaxation factor. The second pass is

φ k+1 i,j = (1 -ω) φ k i,j + ω b i,j -a w φ k+1 i-1,j -a s φ k+1 i,j-1 -a n φ k+1 i,j+1 -a e φ k+1 i+1,j /a c , (2.53)
for all i+j even, respectively. In our 2D multi-grid algorithm, the transfer of data to a coercer grid (so-called "restriction") is performed by a 9-points bilinear interpolation CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS while to a finer grid (prolongation) via a 4-point bilinear extrapolation. Note that once there is less that one node per MPI thread in the direction where the physical domain is decomposed during a V-cycle then the numerical grid is merged between all the MPI thread. The parallelization for the coarsest grids in consequently only achieved by the OpenMP threads. This is a limiting factor and more work is needed to further improve the algorithm. As an example, using a mesh of 512

f * i (x, v, t) = n i N i N i 1 δ (v -V i ) , (2.54)
where n i is the density of specie i, N i is the associated number of macroparticles and V i (t) is the particle velocity. x and t are independent variables. After a lengthy calculation, which can be found in details in refs. [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF][START_REF] Nanbu | Direct simulation scheme derived from the boltzmann equation. i. monocomponent gases[END_REF], one can deduce the probability for an incident particle to undergo an elastic or inelastic collision with CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS a target particle during a time step ∆t

(P i ) max = ∆t Nc c=1 (n c σ c v r ) max , (2.55)
with N c corresponding to the total number of reactions for the incident specie, n c the density of the target specie associated with a given collision index and

v r = |v i -v c |.
The product n c σ c v r in eq. (2.55) is artificially set to its maximum value (n c σ c v r ) max and hence (P i ) max is greater than the real probability and is constant over the entire simulation domain. There is consequently a probability,

(P i ) null = 1 - Nc c=1 P c (P i ) max , (2.56)
that a particle undergoes a fake collision (dubbed "null" collision), which will be discarded.

P c = n c σ c (v r )v r ∆t.
The total number of incident particles which will hence collide during a time step ∆t (including a "null" collision) is,

N max = N i (P i ) max , (2.57)
where N i is the number of incident macroparticles in the simulation. N i must be replaced by (N i -1)/2 for collisions with another particle of the same specie [START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF]. (P i ) max is equiprobable for any pairs of incident-target particles and consequently the latter may be chosen randomly inside the simulation domain. In the model, one checks first if the incident macroparticle experienced a real collision,

r ≤ 1 -(P i ) null , (2.58)
where r is a random number between 0 and 1. The probabilities P c for each reactions (whose total number is N c for a given incident specie) are ordered from the smallest to the largest and a reaction k occurred if,

r ≤ k c=1 P c (P i ) max . (2.59)
Once a collision type is selected then the macroparticles are scattered away in the center-of-mass (CM) frame (see next section). In the model, neutrals are considered as a non-moving background specie with a given density profile. Collisions between charged particles and neutrals are hence performed by artificially extracting a neutral particle velocity from a Maxwellian distribution function. In the case of collisions between charged particles, the actual velocities of both species are used CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS in the calculation of the scattering process.

Elastic and inelastic collision processes

Collisions in the PIC-MCC algorithm (both elastic and inelastic) used in this thesis, are implemented assuming that particles (incident, target or newly created) are scattered isotropically in the center of mass (CM). Energy and momentum is conserved in the model and we posit for simplicity that each byproduct partner after the collision have identical momentum in the CM frame 1 . This implies that the lightest particles will equally share most of the available energy [START_REF] Fubiani | Modeling a high power fusion plasma reactor-type ion source: Applicability of particle methods[END_REF]. Cross-sections for light versus heavy or similarly heavy-heavy particle collisions are often solely function of the relative velocity (especially when originating from experimental measurements), i.e., information about the differential cross-section is lacking. It is the case for nearly all of the cross-sections associated with molecular hydrogen (or deuterium) gas chemistry. Consequently, we implemented a simple MC collision model derived from the isotropic character of a collision. This has the advantage of being versatile (easily adaptable to different types of collision processes both elastic and inelastic) and to conserve exactly energy and momentum. In the center of mass (CM) of the two interacting particles, one assume that each byproduct of the collision have identical momentum, that is,

|p 1 | = |p 2 | = • • • |p n | , (2.60)
where |p n | = m n v n , with m n the mass of the nth byproduct particle and v n its velocity, respectively. Note that the use of Eq. (2.60) impose a strict energy equipartition between particles of equal mass. We find after the collision,

E kr = µv 2 r /2 -E m , (2.61)
with E kr the relative kinetic energy in the CM frame, v r = |v 1 -v 2 | the relative velocity in the laboratory frame, µ = m i m t /(m i + m t ) the reduced mass of the system, m i (m t ) the incident (target) particle mass and E m the threshold energy of the reaction. The relative kinetic energy is shared between all byproduct particles, i.e.,

E kr = m 1 v 2 1 /2 + m 2 v 2 2 /2 + • • • + m n v 2 n /2 , (2.62)
and using the equality defined in Eq. (2.60) , one deduces the kinetic energy of each particle,

m k v 2 k 2 = E kr /m k 1/m 1 + 1/m 2 + • • • + 1/m n , (2.63)
where k is the index of the kth particle. For instance in a system with three particles after the collision, say two electrons and one ion, the electrons share the same and almost all the available energy, that is,

E ke E kr 2 1 - m e 2m i , (2.64)
while the ion takes the remaining part,

E ki E kr m e 2m i 1 .
(2.65)

Lastly, momentum conservation is preserved by assuming equal angle spread between momentum vectors in the CM frame, i.e., θ k = θ 1 + 2π(k -1)/n. The angle θ 1 = arccos(1-2r 1 ) is calculated using a random number r 1 between 0 and 1. The particle velocity in the laboratory frame is derived from,

v k = v CM + v k e k , (2.66)
where v CM = (m i v i +m t v t )/(m i +m t ) is the CM velocity. The vector e k is calculated as follow: first calculate e kx = cos θ k , e ky = sin θ k sin φ , e kz = sin θ k cos φ ,

(2.67)

with φ = 2πr 2 and r 2 is another random number. The angles θ k being generated by one random number, this generates correlations which may be cancelled out by performing an additional rotation of the frame defined by the vector e k . Another set of angles θ and φ must hence be generated to calculate e k . θ , φ and φ are identical for all byproduct particles.

1D analytical solution

In this section we derive a 1D analytical solution to characterise a plasma in the ambipolar regime. The analytical solution will be compared in the next section to a 2D PIC-MCC calculation with periodic boundary conditions in order to model a 1D-like configuration. The flux equation for the ions, neglecting pressure terms, is CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS written as follows [START_REF] Nicholson | Introduction to plasma theory[END_REF]70],

∂n i u i ∂t + ∂n i u 2 i ∂x = en i m i E -n i j ν m,j (u i -u j ) , (2.68)
where ν m is the momentum transfer frequency, E the ambipolar electric field, u i the mean velocity, m i (n i ) the ion mass (density), respectively. Positive ions generated by ionization processes are assumed at rest. The ionization frequency hence does not appear in Eq. (2.68). Adding the continuity equation and Boltzmann electrons,

∂n i ∂t + ∂n i u i ∂x = n i ν i , (2.69)
-eE = T e ∂ ln(n e /n 0 ) ∂x , (2.70)
we have a closed set of equations. ν i is the ionization frequency. Assuming quasineutrality (n i = n e = n), ion-neutral collisions exclusively (we further neglect neutral velocities) and lastly steady state conditions, we find,

u i dn dx + n du i dx = ν i n , (2.71)
u i du i dx + (ν i + ν m )u i = - 1 u 2 B d ln(n/n 0 ) dx .
(2.72)

Normalizing the latter with ũi = u i /u B , x = ν i x/u B , ñ = n/n 0 and z = ln ñ (note that ionization appears in the momentum equation as a friction), we have

ũ dz dx + dũ dx = 1, (2.73) ũdũ dx + (1 + k)ũ = - dz dx , (2.74)
where k = ν m /ν i . Combining the two equations, we get

dũ dx = 1 + (1 + k)ũ 2 1 -ũ2 . (2.75)
ũ varies from 0 to 1. The equation is diverging for ũ → 1 (i.e., u → u B ). Solving for x instead,

dx = (1 -ũ2 )dũ 1 + (1 + k)ũ 2 , (2.76)
and integrating, we find,

x = 2 + k (1 + k) 3/2 arctan ũ√ 1 + k - ũ 1 + k . (2.77)
At the plasma sheath presheath interface we have ũ = 1 and we can deduce an expression for the electron temperature versus k = ν m /ν i ,

h = Lν i 2u B = 2 + k (1 + k) 3/2 arctan √ 1 + k - 1 1 + k , (2.78)
where the sheath length was neglected and L 2x s was assumed (the electron Debye length is of micrometre size in typical cold plasma sources). ν i and u B are both function of T e . For k = 0 (i.e., without any ion-neutral collisions) we find h 0.57. Equation (2.78) can also be derived by integrating the continuity equation, eq. (2.73), over the plasma volume. The equilibrium electron temperature in the plasma can hence simply be deduced from the equality between plasma production inside the volume and losses at the walls in steady state and is a function of the dimensions of the plasma device and the ionisation cross-section. Furthermore, note that from eqs. (2.73) and (2.75), we may deduce the density as a function of ion velocity, that is,

z = ln ñ = - ũ 0 -(2 + k) ũdũ 1 + (1 + k)ũ 2 , (2.79) = - (2 + k) 2(1 + k) ln 1 + (1 + k)ũ 2 .
(2.80) For ũ = 1 (i.e., u = u B ) and k = 0, we find n s = n 0 /2, where n s is the plasma density at the sheath edge. The ambipolar potential φ s at this location is hence,

φ s = φ 0 -T e ln 2 , (2.81)
with φ 0 the potential at the center of the discharge. The height of the sheath potential may be deduced from the particle fluxes to the walls. If one assumes that the plasma sheath is collision-less then the electron and ion fluxes at the sheath interface is equal to the one impacting the walls. The ion flux is simply the Bohm flux

Γ is = n s u B , (2.82)
while for the electrons [START_REF] Golant | Fundamentals of plasma physics[END_REF] Γ

es = n s +∞ √ 2eφs/me v x f e (v x )dv x , (2.83)
where f e (v x ) is the electron distribution function which is assumed Maxwellian where T e is in electron-volts. Taking the ratio of the electron flux to the ion flux, we get

f e (v x ) = m e 2πT e 3
φ s = T e ln m i 2πm e Γ is Γ es . (2.86)
For a configuration where the plasma profiles are symmetric and the walls have the same potential, then Γ is = Γ es and

φ s = T e ln m i 2πm e .
(2.87)

Lastly one can deduce the ion kinetic energy at the wall from the knowledge that the ions have the Bohm velocity at the sheath presheath interface and gain an additional energy eφ s in the sheath

i = T e 1 2 + ln m i 2πm e .
(2.88)

The electron distribution function is assumed to be Maxwellian and the average kinetic energy at the wall is hence e = 2T e .

(2.89)

As an example, integrating eq. (2.78), 2u B h = Lν i , we deduced an electron temperature T e 4.9 eV considering an atomic neutral hydrogen background gas with a density n n = 10 20 m -3 , a plasma of electrons and hydrogen H + ions, collisions consisting only of ionisation (i.e., k = 0 and hence h = 0.57) with a cross-section

σ i =    10 -20 m 2 , if E k ≥ E m = 15.8 eV 0, otherwise (2.90)
and a device length L = 10 cm. E m is the threshold energy. The ionisation frequency in eq. (2.78) may hence be expressed as where φ 0 18V corresponds to the maximum potential at the center of the discharge.

ν i = n n σ i (2πv 2 th ) 3/2 ∞ vm exp -v 2 2v 2 th 4πv

2D PIC-MCC simulation of ambipolar diffusion in a plasma

Next we simulate the plasma described analytically in the previous section with a 2D PIC-MCC model using periodic boundary conditions. We implement a mesh size of 192 nodes in the axial direction and 96 nodes along (Oy). The simulation domain is 10 by 5 cm 2 ; the boundary conditions are of Dirichlet type on the left and right (x = 0 and 10 cm, respectively) with an applied potential V w = 0V. Particles impacting these surfaces are absorbed (removed from the simulation domain). The top and bottom boundary conditions are periodic. A particle leaving the domain at the top for instance will be re-injected at the location y = y -y m where y m is the size of the domain along (Oy). The velocity vector is not changed. For Poisson's equation, we enforce that φ(x, y m ) = φ(x, 0) at each iteration. This configuration allows us to model a 1D plasma with a 2D PIC-MCC algorithm. We use 40 particles per cell and assume an absorbed power of 5 W/m (i.e., 100 W/m 2 in 1D) which generates a plasma with an average density n i = 5 × 10 14 m -3 such that ∆x/λ De = ∆y/λ De = 0.7 in the simulation. Lastly the time step is ω p ∆t = 0.15 corresponding to a Courant-Friedrichs-Lewy (CFL) condition of v th ∆t/∆x = 0.3 where v th is the electron thermal velocity. External power is absorbed by macroparticles as follows: every time step, macroparticles which are found inside the region of power deposition (everywhere in this example) are heated according to some artificial heating collision frequency. Electrons, being the lightest particles, are assumed to absorb all of the external power. Redistribution of energy to the heavier ions and neutrals is done through collisions and the ambipolar potential. Electrons undergoing a heating collision have their velocities replaced by a new set sampled from a Maxwellian distribution with a temperature calculated from the average kinetic energy inside the power deposition region added to the absorbed energy per colliding particles, i.e., 3 2

T h = E k h + P abs eN eh ν h , (2.94)
where T h (eV) is the heating temperature in electron-Volts (eV), E k h is the average electron energy, P abs (W) is the absorbed power, ν h the heating frequency and N eh the number of electrons, respectively. For a given time step, N em ν h ∆t colliding macroelectrons are chosen randomly where N em is the total number of macroparticles inside the heating region. We chose ν h = 10 8 s -1 in this example. Figure 2.9 shows the electron and ion plasma density profiles over half the axial length of the domain.

CHAPTER 2. PARTICLE SIMULATIONS OF MAGNETIZED PLASMAS

One recovers the drop n s = n 0 /2 at the sheath entrance. The latter is found at the location where the plasma is no longer quasi-neutral and the positive ions reach the Bohm velocity. Using eq. (2.89) to estimate the electron temperature, we found T e 5.2 eV which is correlated with the value derived from 2/3 E k inside the plasma volume. The electron temperature is about 6% higher than the analytical estimate. This is due to the fact that in the simulation the sheath occupies ∼ 14% of the volume. This was neglected in eq. (2.78) where we made the assumption that L 2x s (x s is the location of the sheath entrance). Using T e 5.2 eV, we found a good agreement with the analytical theory as shown in table 2.1. The difference is partly explained by the fact that the electron distribution function remains truncated in the model (due to wall losses) beside the mean-free-path of about 1.3 cm for an electron to be heated (which replenishes the Maxwellian distribution). The plasma potential profile is displayed in figure 2.10.

Chapter 3 Plasma confinement by magnetic cusps in low-β plasmas

If a man never contradicts himself, the reason must be that he virtually never says anything at all.

Erwin Schrödinger

In this chapter we study, with the help of 2D PIC MCC simulations, plasma confinement by magnetic cusps. The main issues are described in the introduction chapter, section 1.3, with a brief bibliography about magnetic cusps. We have seen in this section that there is still no consensus on the scaling laws characterizing the confinement by magnetic cusps. In this chapter the goal is to derive clear scaling laws from numerical simulations of magnetic cusps. Although this should be now possible thanks to the availability of powerful computers and simulation tools, we note that there is no recent published work addressing the simulation of magnetic cusps and derivation of scaling laws for the leak width in low-beta plasmas. To our knowledge the only attempt at deducing or checking scaling laws of magnetic cusps with particle simulations was made by Marcus et al. in 1980[72]. At that time Particle-In-Cell simulations were extremely expensive in terms of computation time and only a few test cases were considered by Marcus et al.

More details on previous derivations of the scaling laws and on the objectives of this chapter are given in section 3.1. We discuss in section 3.2 issues related to the nature of the walls (boundary conditions) and to the method of plasma generation. Section 3.3 presents the numerical model. The simulation results are presented in CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS section 3.4. In this section, the leak width, defining the confinement by the magnetic cusps, is studied in argon as a function of magnetic field, pressure, electron and ion temperatures, and ion mass.

Leak width scaling laws and objectives of the chapter

We have seen in section 1.3 three expressions for the leak width of magnetic cusps derived theoretically or inferred from experiments by different authors. The leak width is defined by most authors as the half-width (full width at half maximum) of the profile of the current density to the wall.

Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF] deduced from their experiments that the leak width at a pressure of 0.1 mtorr could be well represented by the expression below, proportional to the hybrid gyroradius, for helium, argon and xenon and for a magnetic field (at the center of the cusp) in the range 50 -200 G:

w L,H = 4(ρ e ρ i ) 1/2 (3.1)
where ρ a = v a /ω ca , a = e, i , v a = (kT a /em a ) 1/2 is the charged particle velocity, and w ca = e/(m a B) is the corresponding cyclotron frequency. Since ρ i ρe = ( T i m i Teme ) 1/2 , the leak width can also be written as:

w L,H = 4ρ e ( T i m i T e m e ) 1/4 ∝ T 1/4 e T 1/4 i m 1/4 i B -1 (3.2)
Note that if the ion cyclotron radius ρ i is expressed for an ion velocity equal to the Bohm velocity C s = (kT e /em i ) 1/2 (i.e. T i is replaced by T e in the expression of the ion velocity above), we can write ρ i = ρ i,Bohm = C s /ωc a , and the scaling of becomes:

w L,H ∝ T 1/2 e m 1/4 i B -1 (3.3)
Considering ambipolar diffusion across the magnetic field and assuming that the plasma in the cusps is flowing outward at approximately the ion acoustic speed, CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF] derived the following expression for the leak width:

w L,B = (2 DR/C s ) 1/2 (3.4)
where D is the electron cross-field diffusion coefficient, and C s is the ion Bohm velocity. In the experiment of Bosch and Merlino the magnetic field was generated by a system of coils producing a spindle cusp magnetic field with a ring cusp and two point cusps. R was the radius of the ring cusp. We generalize this expression to a picket fence or a multicusp geometry with a distance d between cusps, we take d = 2R and the leak width becomes:

w L,B = ( Dd/C s ) 1/2 (3.5)
This expression is actually derived semi-empirically and the coefficient in front of in the expression above is not well defined. The collisional diffusion coefficient D can be written as:

D = kT e e e m ν eN ω 2 ce (3.6)
Introducing the electron mean free path λ eN = v e /ν eN , we get:

w L,B ≈ (d D/C s ) 1/2 = ( d λ eN ) 1/2 (ρ e ρ i,Bohm ) 1/2 (3.7)
The scaling of the leak width in that case is:

w L,B ∝ λ -1/2 eN T 1/2 e m 1/4 i B -1 (3.8)
and is identical to the Hershkowitz et al. scaling if λ eN does not depend on electron temperature (which is generally not the case).

The advantage of the Bosch and Merlino expression of the leak width is that it takes into account the gas pressure. It is clear that the leak weak should increase with gas pressure due to the increased collisional diffusion of electrons across the magnetic field. Since the electron mean free path is inversely proportional to the gas density or to the gas pressure p, the leak width of Bosch and Merlino scales as p 1/2 . Numerically, if the magnetic field in the electron and ion Larmor radii are taken at the maximum magnetic field in the expression of the leak width, the value of w L,B can be actually much smaller than in the experiments. For example, at 0.1 mtorr, the electron mean free path is several meters, much larger than d (10 cm) in the CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS experiment of Bosch and Merlino so that w L,B = ( d λ eN ) 1/2 w L,H w L,H . Bosch and Merlino and other authors (T. Morishita et al [START_REF] Morishita | Estimate of cusp loss width in multicusp negative ion source[END_REF], A. Kumar and V. Senecha [START_REF] Kumar | Cusp loss width in multicusp negative ion source: A rigorous mathematical treatment[END_REF], Cooper et al [START_REF] Cooper | Direct measurement of the plasma loss width in an optimized, high ionization fraction, magnetic multi-dipole ring cusp[END_REF]) argue that in the expression of w L,B , the diffusion coefficient (or the Larmor radii) should not be estimated at the point of maximum magnetic field. Bosch and Merlino estimate the diffusion coefficient at a point midway between the cusp and the filament. Cooper et al. (Cooper et al[52]), mention that the Bosch and Merlino model uses "ill-defined fitting parameters B and d which can always be selected to fit the data over a small scaling". These authors develop a 1D numerical model based on the 0D model of Bosch and Merlino (i.e. with perpendicular crossfield diffusion and outward plasma flow through the cusp at the Bohm velocity). This allows to get a more self-consistent dependence of the leak width on the magnetic field. They find a good scaling with their experiment on the WiPAL multicusp plasma source (see section 3.2).

Finally, a third expression of the leak width was derived by Koch and al. [START_REF] Koch | Collisional diffusion of a plasma in multipolar and picket fence devices[END_REF], and can be written as:

w L,K = 2R π ( ρ e ρ i λ e λ i ) 1/2 (3.9)
The scaling of this expression with ion mass and magnetic field is m

1/4 i B -1 , i.
e. is identical to that of Herskowitz et al. and Bosch and Merlino. However the leak width of Koch and Mathieussent is proportional to the gas pressure p instead of p 1/2 . Note also that the leak width of Koch and Mathieussent is proportional to the distance d between cusps while the expression of Bosch and Merlino varies as d 1/2 . Bosch and Gilgenbach [START_REF] Bosch | Leak widths resulting from plasma diffusion in a magnetic cusp[END_REF] give a simple and interesting discussion of the scaling laws deduced from the different models. For example, they show that the expression w L,K of Koch and Mathieussent can be considered as the continuation of that of Bosch and Merlino,W l,B , when the pressure is large enough so that particle losses parallel to the magnetic field, in the cusps, can no longer be described by ions at the Bohm velocity, but must be represented by parallel ambipolar diffusion. In that case the expression of Bosch and Merlino can be re-written by replacing C s by D /d , with D = kT e /(m i ν iN ) and w L,B becomes:

w L,B ≈ ( Dd/C s ) 1/2 → w L,K ≈ d( m i m e ν eN ν iN ω 2 ce ) 1/2 = d( ρ e ρ i λ e λ i ) 1/2 (3.10)
where we have used λ eN = v e /ν eN and λ iN = v i /ν iN . This expression is similar to that of Koch and Mathieussent. Again, the coefficients in front of these expressions are not accurate since these derivations are very approximative and only the scaling with the different parameters should be considered. The above derivation of Bosch and Gilgenbach provides a good justification of the P 1/2 scaling of the leak width of Bosch and Merlino, in comparison with the P scaling of Koch and Mathieussent. Note also that, for consistency, if the expression of Koch and Mathieussent is used, the velocity in the ion Larmor radius should be the ion thermal velocity and not the Bohm velocity. Therefore, the scaling of w L,K with temperatures is as (T e T i ) 1/4 instead of T 1/2 e for w L,B . It is interesting to look at the ranges of pressure where the leak width w L,B or w L,K should be used, for example in argon. The charged exchange cross-section σ of argon ions is on the order of 4 × 10 -19 m 2 . Therefore, for a pressure of 1 mtorr (gas density on the order of 3.5 × 10 19 m -3 ), the ion mean free path λ iN = 1/(N g σ) is about 7 cm. If the distance between magnets is on the order of 1 to 2 cm, the expression w L,B of Bosch and Merlino should provide a good scaling up to several mtorr, while the expression w L,K of Koch and Mathieussent would be probably better above 10 mtorr (i.e. for pressures where cusp confinement is actually not very efficient).

The objectives of this chapter is to use a 2D PIC MCC simulation to study the scaling laws of the leak width of a line cusp with magnetic field, ion mass, electron temperature, ion temperature and distance between cusps and to compare them with those of the different theories. It is clear that the results may depend on the method of plasma generation and on the charged particles losses to the walls. This question is discussed in the next section in the light of previous experiments and applications.

Plasma generation and particle losses

In most previous experimental studies, the plasma was generated by hot filaments. As said in section 1.3 of the Introduction, the magnetic field in the first experiments (Hershkowitz et al. [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF], K.N. Leung et al. [START_REF] Leung | Plasma confinement by localized cusps[END_REF], Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF]) was generated by parallel conductor wires with currents in adjacent conductors flowing in alternating directions, the so-call "picket-fence" magnetic field (see The picket-fence separates the plasma source from the "target". The plasma source is generated by hot filaments at a negative) potential ( -60 V ) with respect to the grounded chamber walls. Electrons and ions flow from the plasma source to the target through the magnetic cusps formed by the picket-fence. The leak-width, obtained by measuring the profiles of the electron and ion current densities (red line) along the midplane follows a hybrid radius scaling. Hershkowitz and Christensen [START_REF] Christensen | Mass scaling of permanent magnet line cusp plasma leaks[END_REF] showed that the profile of the electron current collected on an equipotential plate placed in the mid-plane along the picket-fence line followed the same hybrid radius scaling. This arrangement was useful for experiments, but practical confinement configurations using magnetic cusps are "multicusp plasma sources" where the magnetic cusps are generated by magnets placed on the chamber walls, as shown in Fig 3 .2. Figure 3.2a shows a schematic of the Kamaboko negative ion source for neutral beam injection in magnetic fusion. This source uses 48 hot cathodes filaments at a negative potential with respect to the grounded chamber walls. Multicusps are also used in positive ion sources for applications to ion thrusters for space propulsion, as shown in Fig 3 .2b.

In these examples the magnetic cusps are placed on the equipotential chamber walls, which play the role of the anode. The anode collects the electron current emitted by the cathode plus the current of plasma electrons generated by ionization.

CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

Some of the ion current can be collected by the anode but the largest part of the ion current is expected to be collected by the hot cathodes. [START_REF] Cooper | Direct measurement of the plasma loss width in an optimized, high ionization fraction, magnetic multi-dipole ring cusp[END_REF]). In this configuration the magnets are covered with an insulating layer. Therefore, the cusp surface cannot collect a net current. This is the reason why supplementary anodes are present in the chamber (see Fig 3 .3). LaB6 emissive hollow cathodes are used as electron sources in WiPAL.

The plasma generation in the multicusp plasma sources described above is made by hot filaments. RF or microwave (ECR) plasma generation is also possible as illustrated in Fig 3 .4 with an example of RF negative ion source.

In the plasma sources described above, except for WiPAL, the magnetic cusps are placed behind a metallic walls that play the role of an anode. The electron current collected by the metallic surface above the magnets is therefore in principle much larger than the ion current. In WiPAL the magnet poles are covered with a dielectric layer and the electron and ion currents collected at each position on this surface must be identical.

Another example where the walls above the magnetic cusps are made of a dielectric material is the High Efficiency Multistage Plasma Thruster -HEMPT) of Thales shown in Fig 3 .5. In this cylindrical thruster the plasma is generated by a discharge between an external emissive hot cathode outside a ceramic cylinder while the anode is inside the cylinder and at the other end. The surface charges on the ceramic walls adjust in such a way that the electron and ion fluxes at each position on the surface are identical. In this chapter we want to study the confinement by magnetic cusps of charged particles in the bulk plasma as much as possible independently of the method of plasma generation. We therefore simplify the problem by assuming uniform plasma generation in a volume sufficiently far from the cusps, i.e. in a region where the magnetic field is small. To simplify the interpretation of the results and to make possible meaningful scaling laws we also impose a Maxwellian electron velocity distribution and a uniform electron temperature. This is described in more detail in section 3.3.

Numerical model

We consider a 2D geometry with line cusps infinite in the direction perpendicular to the simulation domain and periodic in the direction parallel to the walls, as shown in Fig 3 .6. The magnetic field is in the simulation plane (its component perpendicular to this plane is zero). The walls are made of dielectric materials. We summarize below some aspects of the simulation that are specific in this magnetic cusp model. The plasma consists of electrons, positive Argon ion, neutral Ar. The set of reaction is presented in table 3.1. Null collision is performed based on the cross section(from Bolsig+ [START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF]).

Index

Reaction 1 e + Ar → e + Ar (collision) 2 e + Ar → e + Ar * (excitation) 3 e + Ar → e + Ar + +e (ionization) 4

Ar + + Ar → Ar + + Ar (collision) 5

Ar + + Ar → Ar + Ar + (charge exchange) And the cross-section of reaction related with electron is plotted in figure 3.7. The cross-section for elastic collision between electron and Argon is very interesting. It is very small from 0.1 eV to 1 eV , then increases dramaticly from 1 eV to 10 eV and after that decreases again. The threshold energy for excitation and ionization is 11.5 eV and 15.8 eV respectively.

Charging of the dielectric walls

The presence of dielectric walls is approximated in the simulation by capacitances placed at each segment of the grid along the dielectric surface. For example, if a particle reaches the left dielectric surface between the nodes (i = 1, j), and (i = 1, j + 1), the potential at (i = 1, j) is incremented as follows:

V 1,j = V 1,j ± p j αδt (3.11)
where α is a constant (in V /s), δt is the time step, and the ± sign corresponds to ions or electrons. Physically, αδt = δQ/C is proportional to the charge carried by the superparticle and inversely proportional to a capacitance. p j is an interpolating factor associated with the y position of the particle reaching the surface: p j = (y j+1 -y)/(y j+1 -y i ). At the same time the potential at V 1,j+1 is incremented as:

V 1,j+1 = V 1,j+1 ± (1 -p j )αδt (3.12)

magnetic field distribution

We use an analytical expression of the magnetic field derived by Lieberman and Lichtenberg [START_REF] Lieberman | Principles of plasma discharges and materials processing[END_REF] using a first order development in the assumed small parameterwhere ∆/d where ∆ is the width of the magnets and d is their distance (see Fig 3 .6).

CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

In this approximation, the x and y components of the magnetic field due to the magnets on the left side of the simulation domain can be written as:

plasma generation and electron heating

In these simulations the plasma is not generated according to a specific way of electron heating, and ionization is not treated self-consistently. We start the simulation with a given, uniform, plasma density. Each time an ion is lost to the walls, an electron-ion pair is randomly generated in the center of the plasma (the region defined by dashed lines in Fig 3 .6). Nothing is done when an electron is lost to the walls. Doing this, the total number of ions in the domain, i.e. the averaged ion density is kept constant. This procedure is much easier to handle that describing self-consistently the plasma generation since the plasma density is somewhat imposed. In order to impose a constant electron temperature, in the simulation domain, electron present in the "electron heating" region of Fig 3 .6 (identical to the plasma generation region) are heated at a constant rate. This means that their velocity is changed according to a Maxwellian distribution at a constant frequency. Of course, this is equivalent to imposing fictitious collisions in this region, which contributes to un-trapping the electrons from the magnetized CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS trajectories. If the electron heating region is sufficiently far from the walls, the magnetic field is low in that region and the electron heating only represents a small perturbation.

definition of the leak width

The leak width in the experiments is defined as the full width at half maximum of the profiles of electron and ion fluxes to the walls. It is also possible to use this definition in the simulations. Another, more accurate way of defining the leak width in the simulations is simply to calculate the total electron (or ion) current to the walls (per unit length in the perpendicular direction) for a given magnetic field, I B , and to compare it with the current obtained without magnetic field I 0 . The leak width can then be simply defined by:

w L = L y I B I 0 (3.14)
The ratio r L = w L /L y = I B /I 0 fully characterizes the confinement capability of the cusps.

Analysis of the simulation results

Space distribution of the plasma properties

In this sub-section we discuss the space distribution of the plasma properties for a particular but typical case: argon, 0.1 mtorr, B 0 = 400 G, T e = 2 eV , T i = 0.5 eV . Figure 3.9 shows a contour plot of the electric potential around one cusp, in the greyed region of Fig 3.6. We see that the potential along the left dielectric wall is about 3 V above the plasma potential between the cusps and is about 11 V below the plasma potential in the cusps. Therefore, positive ions tends to be pushed away from the wall by the electric field between the cusps and are guided by the potential toward the cusps. Physically, this means that in the regions between the cusps, where the magnetic field tend to be parallel to the walls, electron transport toward the wall is strongly impeded while ion transport is not. Therefore, the electric field tends to attract electrons toward the wall and to repel ions. In the cusps, where the magnetic field is perpendicular to the wall, the situation is opposite. Electrons are moving freely along the magnetic field lines and the field accelerates ions toward the

CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

The electron and ion fluxes to the dielectric wall are very small between the cusps and increases sharply around the cusps, as can be seen in Fig 3 .11 which displays the profiles of electron current density to the walls for different values of the magnetic field B 0 . The ion current density to the wall is identical with an opposite sign (the current to the wall must be zero at each position on the dielectric). As expected, the full width at half maximum of the current density to the wall decreases when the magnetic field increases (a detailed study of the magnetic field dependence of the leak width is done in section 3.4.2).

It is often mentioned, in the literature on magnetic cusps, that the leak width corresponding to high energy electrons tends to be smaller than the leak width of bulk electrons. This is important when the plasma is sustained by hot filaments, since the filaments are usually at a potential significantly lower than the plasma potential, e.g. -60 V with respect to the plasma potential. Therefore, at low pressure electrons of energy up to 60eV are likely to interact with the cusps. In our case, the electron velocity distribution is Maxwellian but it is possible to estimate the leak width as a function of electron energy. In the PIC MCC simulation, this can be done by calculating the electron flux or current density to the wall as a function of electron energy. Figure 3.12 shows the current density profiles on the dielectric surface, of electrons in different energy ranges. The simulations confirm that the leak width decreases with increasing electron energy. For electrons with energy higher than 20 eV , the average leak width is two to three times smaller than the leak width of electrons of energy below 5 eV . Figure 3.13 shows a contour plot of the ion density and figure 3.14 displays the axial profiles of the electron and ion densities at different y positions. The electron and ion density are very close together and very small in the region between cusps. In the cusp (y = 5 mm and y = 5.8 mm, Fig 3 .14) an ion sheath is present next to the dielectric surface. We see that plasma ions start to be accelerated (presheath) toward the cusp at the limit of plasma generation region.

CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

before the decay of the charged particle densities toward the wall. In this region the electron temperature is larger, on the order of 3 eV . This may be due to the fact that the potential increases from the plasma to the dielectric wall, between cusps (see Fig 3.9).

The axial profile of the ion mean energy is shown in Fig 3 .16 at two y locations, in the middle of the cusps (y = 5 mm), and between two cusps (y = 10 mm). In the cusp, we clearly see the ion acceleration toward the wall in the presheath and in the sheath. The sheath entrance is located at about 1 mm from the surface (see the electron and ion density profiles at y = 5 mm in Fig 3 .14). The ion mean energy at the sheath entrance is between 1 and 2 eV , which is consistent with an energy gain of 1 eV in the presheath (the potential drop in the presheath should be on the order of T e /2, i.e. close to 1 eV ). We also note that the ion temperature in the plasma and between cusps is slightly less than 0.4 eV , i.e. significantly lower than the mean energy (3/2T i = 1.5 eV ) corresponding to the temperature (T i = 1.5 eV ) at which ions are injected in the plasma. This is attributed to the fact that the high energy ions of the plasma are lost to the dielectric wall between cusps.

Scaling laws

We performed systematic PIC MCC calculations in order to compare the scaling laws that can be deduced from the simulations to those of the different (empirical) theories. The parameters that can be varied in the simulations are the magnetic field intensity at the surface, B 0 , the gas pressure p, the electron temperature, ion temperature and ion mass. All the calculations were performed in argon (when the ion mass was used as a parameter, the argon cross-section were kept the same, and the ion mass was varied).

We used the numerical model described in section 3.3. As in the previous sections the simulation domain was reduced to the grey region of Figure 3.6 to take into account the symmetry of the problem while reducing the computation time.

In the results below, unless indicated otherwise, the leak width is deduced, as described in section 3.3.4, from the ratio of the total collected electron or ion current , I B , at the dielectric wall with magnetic field to the current collected without magnetic field, I 0 .

w = L y I B I 0 (3.15)
The azimuthal length L y is noted d, in this section, and we plot the dimensionless ratio (leak width normalized to the distance between cusps) w/d. w/d should tend to 1 when there is no confinement (low values of the magnetic field, of high pressure). w should be comparable with the FWHM (Full Width at Half Maximum) leak width measure in the experiments. Figure 3.17 displays the variations of the normalized leak width with magnetic field, for different values of the pressure in argon, and for fixed electron and ion temperature (2 eV and 0.5 eV , respectively). The hybrid radius (ρ e ρ i ) 1/2 is also shown for comparisons.

As expected, w/d → 1 when the magnetic field goes to zero (no confinement). The decrease of leak width with magnetic field at constant pressure does not scale perfectly with 1/B , as seems to be the case in the experiments and "theory" reported by Hershkowitz et al. or Bosch and Merlino (see the discussion of section 3.1) at low pressures. We recall here that the measurements of Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF] were in agreement with the expression w L,H = 4(ρ e ρ i ) 1/2 of the leak width , while Bosch and Merlino proposed a leak width of the form w L,H = (2 D/C s ) 1/2 . Both expressions lead to a 1/B dependence of the leak width with the magnetic field. The results of Figure 3.17 show that the leak width variations at low pressure (0.1 mtorr) could be fitted by a B -1 curve at low magnetic fields, but that the variations with are slower than B -1 (and get closer to B -1/2 ) at high magnetic fields.

At higher pressure (see the 1.6mtorr case in Figure 3.17) the variations with B are much slower than 1/B at low magnetic fields because the plasma is more collisional and the confinement is less effective.

The results are plotted as a function of gas pressure p for different values of the magnetic field in Figure 3.18. We see that the normalized leak width scales well with p 1/2 at high enough pressure. This is consistent with the scaling of the expression of the leak width provided by Bosch and Merlino, w L,B = (2 D/C S ) 1/2 where D is the classical, collisional, electron diffusion across the magnetic field and is proportional to the electron collision frequency, and hence to the gas pressure. The p dependence if the ion velocity in the hybrid radius is taken as the Bohm velocity. In the expression of Bosch and Merlino, the leak width scales as T 1/2 e with electron temperature assuming a constant electron mean free path. It is therefore interesting to look at the predictions of the PIC MCC simulation for the dependence of the leak width with electron and ion temperature. Figure 3.19 displays the normalized leak width as a function of electron temperature for different combinations of magnetic field and gas pressure. We see in this figure that the variations of the leak width with electron temperature is slightly faster than T 1/2 e in some cases and is always much faster than T 1/4 e . As discussed in section 3.1, the expression w L,B = (d/λ eN ) 1/2 (ρ e ρ i,Bohm ) 1/2 of Bosch and Merlino varies as T 1/2 e due to the (ρ e ρ i,Bohm ) 1/2 term if the electron mean free path λ eN is constant. In argon λ eN decreases with electron temperature because of the Ramsauer minimum so the fact that the calculated leak width is faster than T 1/2 e is consistent with the theory of Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF].

We also performed simulations with a constant electron temperature of 2 eV , and for ion temperatures between 0.05 eV and 1 eV . The leak width was found to be practically independent of ion temperature. Therefore, the use of the hybrid radius gives a reasonable order of magnitude of the leak width (see Figure 3.17) if the ion velocity in the expression of the ion Larmor radius is taken as the Bohm velocity and not the thermal velocity. Note that Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF] used an expression of the leak width equal to 4 times the hybrid gyroradius. This is probably because they used the thermal ion velocity in the expression of the ion Larmor radius. Using the Bohm velocity gives a reasonable order of magnitude of the hybrid radius leak width but without the factor of 4 of the expression of Hershkowitz et al., as can be seen in Figure 3.17.

conclusion

In this chapter we have addressed the question of confinement by magnetic cusps in a low-beta plasma, using 2D PIC MCC simulations. The leak width w of a line cusp is an effective loss length, i.e. the average particle loss flux (or current, I B ) to the wall is equal to the average loss flux (or current, I 0 ) without magnetic field times the ratio of the leak width to the cusp length d. This can be written as w/d = I B /I 0 . The confinement is better for lower values of w/d.

CHAPTER 3. PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

A usual but empirical expression of w/d, proposed by Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]. and in plasma textbooks (see, e.g., M.A. Lieberman and A.J. Lichtenberg [START_REF] Lieberman | Principles of plasma discharges and materials processing[END_REF]) is proportional to the hybrid gyroradius 4(ρ e ρ i ) 1/2 . However this expression does not take into account the effect of collisions, which can significantly enhance the charged particle losses to the wall when the pressure increases, due to diffusion across the magnetic field. The goal of this chapter was therefore to provide a way to quantify the losses of charged particles to the wall in the presence of magnetic cusps.

We considered a 2D periodic simulation domain with dielectric walls and with a plasma source sufficiently far from the walls. Electrons and ions were generated with Maxwellian distributions at given temperatures in the source, and the electron temperature in the source was artificially maintained by performing a "Maxwellianization" of the electron velocity at a given frequency. The simulations were performed at constant averaged ion density, and ionization was not treated self-consistently: one electron and ion pair was generated in the source region each time an ion was lost to the wall. A parametric study of the calculated leak width was performed, the parameters being the magnetic field at the dielectric surface, the gas pressure, the electron temperature, ion temperature, and ion mass. The results are summarized below.

The hybrid gyroradius (ρ e ρ i ) 1/2 provides a reasonable order of magnitude of the leak width without the factor of 4 of the expression w L,H of Hershkowitz et al. and if the ion velocity in ρ i is taken as the Bohm velocity. The leak width is practically independent of the ion temperature T i and scales relatively closely to T 1/2 e with the electron temperature. The hybrid gyroradius only gives an order of magnitude, and the 1/B scaling of the hybrid gyroradius is not satisfied in the calculated leak width. The simulations predict a 1/B dependence of the leak width only at low pressure and low magnetic fields. The variations of the leak width with magnetic field are slower than 1/B at high magnetic fields.

The leak width becomes practically independent of pressure below 0.1 mtorr and increases as p 1/2 with pressure above 0.5 mtorr. This is consistent with a leak width proportional to the square-root of the cross-field diffusion coefficient, as in the theory of Bosch and Merlino, which predicts a leak width proportional to

w L,B = (2 D/C s ) 1/2 .
The expressions of the leak width of Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]. and of Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF] both predict a M 1/4 dependence with the ion mass. This comes from the square root of the inverse of the ion acoustic velocity in the expression of Bosch and Merlino, and from the square root of the ion Larmor radius in the hybrid gyroradius In the simulations with dielectric walls presented in this chapter, the electron and ion fluxes must be identical at each point of the walls. This was not the case in the experiments of Hershkowitz et al. and Bosch and Merlino where the charged particle fluxes were measured in a picket-fence configuration (the differences between the local electron and ion fluxes was shown to decrease when the plasma density was increased). In the simulations, the perpendicular electric field on the dielectric wall is attractive for ions in the cusps (as in a regular un-magnetized sheath) and becomes attractive for electrons between cusps. Therefore, between cusps, electrons are accelerated to the wall by the electric field while ions are repelled. The electron flux to the wall is therefore defined by the (very small) cross-field mobility and diffusion, while the ion transport to the wall is defined by the ion pressure gradient.

Although we tried to study the leak width over a large range of parameters, more investigations would be useful. More systematic simulations are needed to understand the dependence of the leak width on the ion mass, to better characterize the role of the specific Ramsauer shape of the electron momentum cross-section in argon (it would be useful to perform simulations with constant collision cross-section or constant collision frequency) and to study the dependence of the leak width the cusp length. Finally, we have checked the numerical convergence of the simulations by performing simulations with increasing number of particles per cell and with mesh refinement. In the conditions of the simulations with dielectric walls presented in this chapter convergence have been checked and confirmed on several points of the parametric study. We also performed systematic simulations and convergence tests for conditions where the walls are not dielectric but are metallic (and therefore equipotential). Surprisingly, we had much more difficulties to reach convergence with the number of particles and with the mesh size in these conditions. We do not understand, at this moment, the reasons for this convergence problem and we leave this question open for further studies.

Chapter 4 Electron transports across magnetic field in negative ion sources

Let the future tell the truth, and evaluate each one according to his work and accomplishments. The present is theirs; the future, for which I have really worked, is mine.

Nikola Tesla

Introduction

This work is relevant to the question of co-extraction of electrons from negative ion sources (and in particular of those used for fusion applications [START_REF] Speth | Overview of the rf source development programme at ipp garching[END_REF][START_REF] Mcneely | A langmuir probe system for high power rf-driven negative ion sources on high potential[END_REF][START_REF] Kraus | Deuterium results at the negative ion source test facility elise[END_REF][START_REF] Wünderlich | Formation of large negative deuterium ion beams at elise[END_REF][START_REF] Kashiwagi | Development of negative ion extractor in the high-power and long-pulse negative ion source for fusion applicationa)[END_REF][START_REF] Ikeda | Identification of the extraction structure of hions by h α imaging spectroscopy[END_REF]) and of micro-ECR neutralizers in space propulsion [START_REF] Takao | Electron extraction mechanisms of a micro-ECR neutralizer[END_REF][START_REF] Sato | Electron extraction enhancement via the magnetic field in a miniature microwave discharge neutralizer[END_REF]. In the former the magnetic cusp field (also called suppression field) is used as a barrier to filter out electrons in order to minimize the electron current co-extracted with negative ions toward an electrostatic accelerator (with a 1 MV acceleration voltage in the case of ITER). Electrons in this case can be damaging to the accelerator because they can induce a high parasitic power density deposition on the accelerator grids [START_REF] Fubiani | Analysis of the two accelerator concepts foreseen for the neutral beam injector of the international thermonuclear experimental reactor[END_REF].

We analyse the transport of electrons across an aperture surrounded by a magnetic cusp field profile in a simplified model without negative ions.. The electrons

CHAPTER 4. ELECTRON TRANSPORTS ACROSS MAGNETIC FIELD IN NEGATIVE ION SOURCES

are strongly magnetized by the cusp field while ions are somewhat magnetized but only over a short distance very close to the aperture surface. In this model, the plasma is generated on the left-hand-side (LHS) of the simulation domain which is separated from a biased electrode by a grid with apertures. The biased electrode is necessary to generate an electric field to extract the electrons. A meniscus forms around the aperture and separates the quasineutral plasma from the non-neutral region between the grid and the biased electrode. The cusp field amplitude decreases exponentially from the permanent magnet position, which is located in the model on the biased electrode inside the accelerator region (same as in ITER and DEMO accelerators).

In practice, the cusp field is generated by two permanent magnet bars of reversed polarity parallel to the slit apertures. The magnetic vector field hence lies approximatively in the plane perpendicular to the magnets when the latter are much longer than the width of the aperture. We find that the electron flux diffusing across the aperture is driven by an E × B drift where E is the electric field in the vicinity of the plasma meniscus. The E × B drift occurs in the plane perpendicular to the cusp magnetic field lines. As a consequence, in an hypothetical configuration where a slit aperture is infinite, this E × B drift does not contribute to electron extraction. The flux of electrons through the aperture is much smaller in that case and is only due to collisional cross-field drift. We show that a 2D Particle-In-Cell (PIC) algorithm with Monte-Carlo-Collisions (MCC) is sufficient to capture the underlying physics of an infinite slit, and is consistent with results from a 3D PIC-MCC model. In the case of a finite slit, only a 3D model can describe properly electron extraction due to E × B drift through the aperture In the next section, we describe in details the simulation parameters. In section 4.3, we analyse the mechanisms of electron extraction in the case of an infinite slit and for a finite slit aperture. In section 4.4, we describe some micro-instabilities occurring in the vicinity of the plasma meniscus in the case of a finite slit aperture. Lastly, section 4.5 and 4.6, provide a discussion and conclusions on the work presented in this thesis.

Overview of the 2D and 3D PIC-MCC models

We use the 2D and 3D PIC-MCC electrostatic explicit models described in details in the previous chapters. For better readability we summarise below the main characteristics of the algorithm. The electric field is calculated on the mesh nodes by solving Poisson's equation iteratively via a multi-grid algorithm [START_REF] Briggs | A Multigrid Tutorial[END_REF]. The magnetic cusp field profile is derived from analytical considerations [START_REF] Lieberman | Principles of plasma discharges and materials processing[END_REF],

B x = B 0 sin π (z -z 0 ) d exp - π (x 0 -x) d , (4.1)
B z = B 0 cos π (z -z 0 ) d exp - π (x 0 -x) d , (4.2)
where B 0 is the maximum amplitude of the field, d is the distance between the magnet bars, x 0 is their axial coordinate and z 0 the center position between the two magnets, respectively. The norm of the magnetic field vector, |B|, is decreasing exponentially from the location of the magnets. The particles trajectories are calculated by solving the Newton's equation of motion at the particle location (the electric and magnetic fields are interpolated),

m a dv a dt = q a (E + v a × B) , (4.3)
where a either stands for electrons or ions, E(x, y, z) is the electrostatic field, v a (t) the particle velocity and m a (q a ) its mass (charge), respectively. The force on the RHS is the Lorentz force. In a PIC-MCC algorithm, charged particles are actually macro-particles but with the same charge over mass ratio as real electrons and ions.

Once the particle positions have been updated, the density is calculated on the mesh nodes. This is used as a source term for Poisson's equations. Lastly, recall that collisions between charged particles and the neutral background gas are implemented using the Monte-Carlo technique [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF][START_REF] Nanbu | Probability theory of electron-molecule, ion-molecule, moleculemolecule, and coulomb collisions for particle modeling of materials processing plasmas and cases[END_REF][START_REF] Nanbu | Direct simulation scheme derived from the boltzmann equation. i. monocomponent gases[END_REF][START_REF] Vahedi | A monte carlo collision model for the particle-in-cell method: applications to argon and oxygen discharges[END_REF]. The PIC-MCC cycle described above is repeated iteratively. The collision processes are simplified in this problem. We considered only the elastic collisions between electrons and the positive ions or the neutrals (table 4.1). The neutral background gas profile is assumed constant with a density of n n = 5×10 19 m -3 . The average plasma density is n i = 5×10 15 m -3 and the plasma generation is modelled by mimicking ionisation through the reinjection of an electron-positive ion pair for each positive ion lost on the physical boundaries of the simulation domain (i.e, of Dirichlet type). Note that the electron Debye length is typically about eight times larger than in the extraction region of negative ion source used in ITER and DEMO [START_REF] Wimmer | Extraction of negative charges from an ion source: Transition from an electron repelling to an electron attracting plasma close to the extraction surface[END_REF][START_REF] Kraus | Deuterium results at the negative ion source test facility elise[END_REF]. The particles are injected into the gray area shown schematically in fig. 4.1 (dubbed "plasma source"). The electron

CHAPTER 4. ELECTRON TRANSPORTS ACROSS MAGNETIC FIELD IN NEGATIVE ION SOURCES

distribution function is artificially maintained in that area as a Maxwellian at a temperature of T e = 2 eV via the periodic replacement of the electron velocity (we use an artificial "heating" frequency). The positive ions are injected assuming also a Maxwellian distribution at a temperature T i = 2 eV (although we do not enforce it in the plasma volume, only at the time when the ions are generated). The latter is the typical value for the average kinetic energy deduced from 3D PIC-MCC models of a fusion-type negative ion source [START_REF] Fubiani | Role of positive ions on the surface production of negative ions in a fusion plasma reactor type negative ion source-insights from a three dimensional particle-in-cell monte carlo collisions model[END_REF]. We implemented a numerical resolution of 256 × 384 × 192 grid nodes associated with a simulation box of physical dimensions 3.2 × 4.8 × 2.4 cm 3 (length versus width and height). The number of particles per cell is 40 and the time step ω p ∆t 0.2. In addition, the size of the two slits is 4.3 × 0.8 cm 2 (length times width) and the amplitude of the cusp magnetic field on the plasma grid (PG) is 200G, which corresponds to B 0 = 875G, d = 1.2 cm, x 0 = 3.2 cm, z 0 = 0.6 cm in eqs. (4.1) and (4.2), respectively. The cusped magnetic field profile encompasses each of the slit apertures, as shown in fig. 4.1. The direction of the magnetic field being reversed between slits. In addition, the numerical model is parallelelized in an hybrid manner using the OpenMP and MPI libraries [START_REF] Fubiani | Modeling of plasma transport and negative ion extraction in a magnetized radiofrequency plasma source[END_REF]. For the simulations reported in this paper, we typically used 48 OpenMP threads together with 4 MPI threads (one per node), i.e, 192 cores in total (on Intel Xeon CPU E5-2697 v4 @ 2.30GHz). The top, bottom, left and right boundary conditions are periodic and we are hence modelling the equivalent of an infinite number of slit apertures. Lastly, the space between the PG and the extraction grid (EG) is free of plasma (the gap is 3.5 mm and the width of the PG is 2 mm). The LHS of the simulation domain, PG, and EG are of Dirichlet type. The EG is biased at a voltage of 400V in the simulations while the other absorbing boundaries (PG and LHS) are set at 0V. In this section, we compare the plasma characteristics derived from a 2D versus 3D PIC-MCC model of a slit aperture of infinite length. The simulation domain for the 2D calculation is shown in fig. 4.2. We model the plane which includes the cusp field lines. The particle drifts are out of the simulation domain and drift-wave type plasma instabilities are hence not described. The electrons which are extracted toward the accelerator diffuse across the field lines solely through collisions.

Next, we compare the 2D calculation to a 3D PIC-MCC model where the simulation domain is identical to the one of fig. 4.1 except that the slit aperture is of infinite length along (Oy); its width being 8 mm in the z-direction. Any occurence of plasma instabilities in the plane perpendicular to the magnetic field lines, which could potentially increase the electron current transmitted across the aperture, will be modelled in this case. We found instead that the electron current collected on the extraction grid (EG) is identical in 2D and 3D, as shown in fig. 4.3. The latter plots the ratio (in percent) of the electron current impacting the EG divided by the total current flowing toward the PG inside the plasma volume. The level of numerical noise is logically higher in the 2D calculation due to the larger weight of the macroparticles. To further support the observation that the plasma properties are identical in 2D and 3D in the case of slit apertures of infinite length, we show in fig. 4.4 the plasma potential profile along (Ox) for a location on the PG in the middle of the interstice between slit apertures. We compare the results from the 2D calculation versus estimates from the 3D PIC-MCC algorithm for either the case with two apertures as in fig. 4.1 or using a reduced domain in 3D where we modelled only one cusp, i.e., half the domain of fig. 4.1 along (Oz).

3D PIC-MCC model of a slit aperture of finite length

In this section, we analyse the electron kinetics for the case of a slit aperture of finite length. The simulation domain is described in details in section 4.2 and figure 4.1. occurs (XY plane at Z = 0.6 cm, crossing the middle of the slit aperture), (iii) the front of the PG and, lastly, (iv) the EG, respectively. In (i), one can see the telltale sign of the electrons oscillating along the cusp field lines while in (ii), the electrons drift in the direction perpendicular to the magnetic field (E × B drift, see figure 4.6) and they can exit toward the accelerator only through one side of the slit aperture. Figure 4.6(a) shows one electron trajectory in the Z = 0.6 cm plane assuming that the Z component of its velocity is zero. When the electron reaches the edge of the plasma grid aperture, the electric field increases rapidly due to the large voltage of the extraction grid. This large non-uniformity of the electric field generates an electron drift across the magnetic field, i.e. the electron guiding center is no longer confined along an equipotential line as can be seen in Fig. 6. This is a finite Larmor radius effect due to the electric field non uniformity (see for instance FF. Chen [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF] ).

In the example of figure 4.6(a), the electron trajectory does not reach the extraction grid and ends up between the plasma grid and the extraction grid. In the examples of figure 4.6(b), the electrons reach the extraction grid. The trajectories in this figure are calculated assuming that the electrons move both along the magnetic field toward the upper edge of the plasma grid of figure 4.5 and drift in the E × B direction in a XY plane toward the right edge of the plasma grid. We assume, for simplicity and to illustrate our point, that the potential does not depend on Z, i.e. the electrons see the same XY distribution of the potential as in the plane Z = 0.6 cm. B z decreases when the electron moves along the magnetic field in the direction of increasing Z away from the plane Z = 0.6 cm meanwhile the X component of the magnetic field, which is perpendicular to the EG, increases. The trajectories of figure 4.6(b) are calculated with the B z component of eq (4.2) only. Therefore, while the electrons drift in the E × B direction toward the right aperture edge, the magnetic field that confines their trajectories in the XY plane decreases and the large and increasing electric field on the edge of the grid aperture extracts them to the EG. The reason why the maximum electron current on the extraction grid is seen on the two right corners of the aperture (see figure 4.5) is therefore due the combined effect of the E × B drift and to the electron motion along the magnetic field lines which spreads the drift current along (Oz) toward the corners 5) while drifting in the E × B direction toward the right edge of the PG (the B z component of the magnetic field, perpendicular to the XY plane, decreases during this motion). For simplicity, these trajectories are calculated assuming that the XY distribution of the potential does not vary while the electrons moves along the magnetic field line (i.e. the electric potential does not depend on Z).

CHAPTER 4. ELECTRON TRANSPORTS ACROSS MAGNETIC FIELD IN NEGATIVE ION SOURCES !"# !"$ %"# # & ! % ' ()*+,- .)*+,- !"# !"$ %"# ()*+,- Φ !"# $%% &%% % ' ( ) * &% $% !+# !,#
of the aperture where the electrons are extracted to the EG.

Comparing the case of a slit aperture of infinite (described in section 4.3.1) versus finite length, one observe than in the latter the electron current extracted toward the accelerator vessel is greatly enhanced by the E × B drift motion in front of the plasma meniscus. The ratio of the electron current impacting the EG divided by the total current flowing into the PG area (inside the plasma volume) is shown in figure 4.7. We find that the current on the EG strongly varies versus time within a range for the transmission ratio between 50% and 100% (i.e., all the current). This is to be compared to the ∼ 22% in the case of a slit aperture of infinite length. Such a large extracted current is induced by a plasma potentiel much higher than the average electron kinetic energy in the model as shown in figure 4.4 (the temperature is T e = 2 eV).

Unstable plasma behaviour in the vicinity of the PG

The large time variations of the current transmitted toward the EG in figure 4.7 are due to strong fluctuations in the plasma quantities in the vicinity of the plasma meniscus. Figures 4.8 and 4.9 plots the electron flux and plasma potential profiles versus several time shots which are also displayed in figure 4.7 for reference. The time window corresponds to approximatively two periods of oscillations. One can see that the flux profile is strongly modified, which is correlated with stiff changes of the plasma potential. We observe that (i) the electron flux follows approximately the potentiel lines in the aperture mid-plane plane as discussed previously in section 4.3.2 (see also figure 4.6), (ii) case A and I (D et G), which correspond to a maximum (minimum) transmission rate, share similar features, respectively. In D et G, the isopotential lines are twisted and we observe the appearance of potential wells (which attract positively charged ions). (∼ λ De ), the space charge (n e /n i -1) switches from globally positive (excess of electrons) to negative. The direction of the electron flux is upward in figure 4.8 and electrons experience a strong E × B drift in the vicinity of the aperture top edge in cases D et G (the 10V and 20V isopotential lines are very close in that area, hence the amplitude of the electric field is large, of order 20 kV/m for a magnetic field strength of 200G and a drift velocity v D ∼ 10 6 m/s, respectively). Such a situation is accompanied by a surge in transmission rate, as shown in figure 4.7. In cases A and I, the electric field strength (and hence the drift velocity) is much lower in front of the aperture, which explains the drop in the electron current extracted toward the EG. We can conclude that the current oscillations observed in the simulations are due to a drift instability induced by the charge separation associated with the electron loss to the extraction grid.

Discussion

We modelled in this chapter a slit aperture of either finite or infinite length. We also performed simulations for cylindrical apertures and reached similar conclusions. In addition, for the case of negative ion sources for fusion applications, there is typically a large negative ion to electron ratio in the vicinity of the aperture. In Christ-Koch et al. for instance [START_REF] Christ-Koch | Laser photodetachment on a high power, low pressure rf-driven negative hydrogen ion source[END_REF] the measured ratio n -/n e is about 0.4 in hydrogen and as high as 1.5 in deuterium 2 cm from the PG, inside the plasma volume. In these sources, the negative ions are produced on the surface of the PG as a byproduct of the impact of neutral hydrogen or deuterium atoms on the grid. The latter is covered by a thin layer of caesium to lower the work function of the metal increasing, as a consequence, the negative ions production yield. The adjunction of negative ions in our model did not modify the conclusions that an E×B drift near the plasma meniscus drives the electron current across the aperture toward the EG.

Lastly, in section 4.3.2, we found that a high fraction of the electron current flowing from the plasma volume toward the PG may be extracted (figure 4.7). The extracted current can be significantly lowered in the model by biasing the PG with a voltage positive with respect to the ion source walls. In that case, the potential difference between the biased PG and the sheath potential is reduced and hence the electron current collected on the grid may be significantly increased. This technique is commonly employed in fusion-type ion sources used for Neutral Beam Injection (NBI) [START_REF] Wimmer | Extraction of negative charges from an ion source: Transition from an electron repelling to an electron attracting plasma close to the extraction surface[END_REF] .

Conclusion

We analysed the transport of electrons across a grid with a slit aperture surrounded by a magnetic cusp field profile. An electrode biased positively with respect to the grid is placed behind in order to produce an electric field which can extract negatively charged particles. In the case of negative ion sources, the role of the cusp field is to act as a barrier, preventing electrons from being extracted. The slit was either of finite or infinite length in the model. In the former, one found identical results between 3D and 2D PIC-MCC calculations indicating that the transport across the aperture is solely collisional. The plane which includes the cusp field lines was simulated in the 2D PIC-MCC model, hence any transport due to the magnetised drift dynamics was not included. In the 3D case, we found that the electron current transmitted onto the bias electrode (EG) was greatly enhanced CHAPTER 4. ELECTRON TRANSPORTS ACROSS MAGNETIC FIELD IN NEGATIVE ION SOURCES compared to the 2D simulations. This is due to an E × B drift motion for the electrons where E is the electric field of the plasma meniscus. The electrons escape parallel to one of the two edges of the slit aperture, being the only location where the drift velocity is directed toward the EG (the drift is generated by the B z and E y components in our configuration). Furthermore, we observed that the plasma parameters (potential, density, electron flux) are strongly fluctuating in front of the plasma meniscus resulting in a transmitted electron current ratio which varies within a factor of two over a time scale of 1 µs. The maximum corresponds to an extraction of most of the electrons flowing onto the PG from the plasma volume in the configuration that we modelled. These oscillations are microinsatbilities due to the charge separation induced by the loss of electrons to the extraction grid.

Chapter 5

Ionization instability in magnetron discharges

We have seen in the Introduction that instabilities are ubiquitous in magnetron discharges and Hall thrusters. In this chapter we focus on one of the important types of instabilities that are observed in these devices, called "Rotating spoke" instability. The rotating spoke is a macroscopic (i.e. long wavelength) instability that can be easily detected by CCD imaging. This instability appears for sufficiently large magnetic fields and rotates in the azimuthal direction of these E × B devices at velocities typically less than 10 km/s. As mentioned in the introduction, an interesting and intriguing feature of this instability is that is can rotate either in the +E × B direction or in the -E × B direction.

To study this instability we have chosen the conditions of an experiment performed in the group of M. Cappelli at Standford University (T. Ito, C. V. Young, and M. A. Cappelli [START_REF] Ito | Self-organization in planar magnetron microdischarge plasmas[END_REF].) This experiment is described in section 5.1. In section 5.2.1 we present the Particle-In-Cell Monte Carlo Collision model that is used to simulate this experiment. The results are analyzed in section 5.2.2. A parametric study of the influence of the cathode magnetic field is performed and analyzed in section 5.3 under conditions of dimensions, pressure and magnetic field scaled by a factor of 10 (d × 10, p/10, B/10) with respect to the conditions of the Stanford experiment.

Rotating instability in magnetron discharges

The experiment of Stanford University is relatively simple and very nice self-organized rotating structures were observed in the conditions of that experiment. This is the reason why we chose to model and simulate this experiment. Figure 5.1 shows the discharge geometry, magnetic field, and optical diagnostics used in this experiment. The dimensions of the magnetrons are small: 2 mm gap and discharge diameter on the order of 5 mm. The operating pressure is large (20 P a, i.e. about 0.15 torr) and the magnetic field decreases from about 1 T at the cathode surface to 0.1 T at the anode. Note that according to the basic discharge similarity laws and assuming that non-linear effects (such as electron-ion recombination or stepwise ionization) are negligible, a discharge operating at a pressure 10 times lower, with dimensions 10 times larger and magnetic field intensity 10 times lower would have similar properties as this miniature magnetron. These scaled conditions (i.e. 10 mtorr, 2 cm gap, 5 cm diameter, magnetic field intensity of 0.1 T at the cathode and 10 -2 G at the anode) would be closer to those of usual magnetron discharges.

In this miniature magnetron, Ito et al [START_REF] Ito | Self-organization in planar magnetron microdischarge plasmas[END_REF]. observed the formation of well-defined regions of enhanced luminous intensity rotating in the -E × B direction as shown in Figure 5.2. In these conditions of applied dc voltage (261 V ), the rotating structures exhibits 5 spokes rotating in the -E × B direction (mode m = 5). The number of modes decreased from 5 to 3 when the applied voltage was increased from 261 V to Simulating the whole length of the azimuthal direction (about 15 mm) would be computationally expensive so we restrict the azimuthal dimension to 4 mm and we assume periodic boundary conditions on the azimuthal domain boundaries. Of course, this imposes a maximum wavelength to the simulated instabilities, i.e. this model cannot reproduce rotating spokes with wavelength larger than 4 mm. The axial azimuthal domain is therefore 2 mm x 4 mm. The magnetic field is perpendicular to the simulation domain. Electrons are emitted from the cathode (secondary electron emission) due to ion bombardment with a secondary emission coefficient γ(for one ion reaching the surface, on the average γ electrons are emitted). Since the detailed mechanisms of electron emission in the conditions of large magnetic field at the cathode surface are not well understood we use a net secondary electron emission coefficient. Electrons emitted by the cathode return to the cathode surface due to their cyclotron trajectory and can be reflected or absorbed by the cathode. Practically, we assume that each electron coming back to the cathode is re-emitted ("net" emission coefficient). This net emission coefficient has been adjusted to 0.005 to get results that are consistent with experiments. The electrons are emitted by the cathode according to a Maxwellian distribution at 2 eV . The simulations were started with a uniform initial plasma density in the range 2 -5 × 10 16 m -3 .

In these conditions steady state was reached in less than 20 µs and the maximum plasma density at steady state was around 5 × 10 17 m -3 . The simulations in the conditions of Figure 5.3 were performed on a 128 × 256 grids. The number of particles per cell was on the order of 200.

CHAPTER 5. IONIZATION INSTABILITY IN MAGNETRON DISCHARGES

Simulations in these conditions were first performed by JP Boeuf and M Takahashi (JP Boeuf and M Takahashi [START_REF] Boeuf | Rotating spokes, ionization instability, and electron vortices in partially magnetized e × b plasmas[END_REF]). We have reproduced these simulations in section 5.3.1. In section 5.3.2 we present a parametric study of the influence of the magnetic field profile (these simulations are performed with a longer azimuthal dimension and on a 128 × 512 grids).

Description of the simulated rotating spoke

In the conditions described above, non-uniformities in the azimuthal direction appear in the simulation in a few 100 ns. After a few µs, the results exhibit a relatively stable azimuthal non-uniformity that propagates in the -E × B direction. We call this instability "rotating spoke" in the following since its properties are similar to those of the experimental rotating spoke described by Ito et al. and by other authors (see e.g. the paper by M. Panjan and A. Anders [START_REF] Panjan | Plasma potential of a moving ionization zone in dc magnetron sputtering[END_REF], and the reviews by A. Anders [START_REF] Anders | Tutorial: Reactive high power impulse magnetron sputtering (rhipims)[END_REF] and A. Hecimovic and A. von Keudell [START_REF] Hecimovic | Spokes in high power impulse magnetron sputtering plasmas[END_REF] ).

The simulated instability is characterized by a region of enhanced ionization moving in the retrograde E × B direction at a velocity of about 10 km/s. The rotating spoke is due to an ionization instability. It is remarkable that most of the ionization takes place in the spoke. In these conditions of large magnetic field in the cathode region, the electrons do not gain sufficient energy in this region to significantly ionize the gas (the energy gain term n e v e E in the electron energy equation is proportional to ν eN /B 2 ). In the cathode sheath, the electron temperature is large because of the larger electric field there but the electron density is small so ionization is also small in the sheath. The reason for electron heating in the spoke is analyzed in section 5.3.4.

The equipotential contours in Figure 5.4c show that the plasma is divided into two distinct regions. A region of quasi-constant potential close to the anode potential (the red region if Figure 5.4c) and a region where the electric field is relatively large, between the cathode and the quasi-equipotential region. The low field, quasi-potential region, and the larger field region are clearly apparent on the contour plot of the ion mean energy in Figure 5.4f. Some local maxima of the potential can be seen in the quasi-equipotential region (contour lines in ) . These maxima are associated with "holes" in the electron density. These electron holes and the associated potential maxima move along the equipotential line separating the two regions. Note that since the electric field is perpendicular to the equipotential lines, the E × B electron drift follows the equipotential lines.

There is a deficit of electrons along the line separating the two regions mentioned above. This is because this line is connected to the anode and electrons are therefore lost at the points where this line intersects the anode. The deficit of electrons can better seen on Figure 5.5. The ion density along this line is significantly larger than the electron density and we can call this line a double layer. The electron holes moving along the double layer are actually electron vortices that are formed because of the velocity shear in the double layer. The vortices move along the double layer at a velocity E/B where E is an average value of the electric field in the double layer. The physics of these electron vortices is very similar to that of the electron vortices observed in the diocotron instability of pure electron plasmas. Electron vortices in pure electron plasmas generally correspond to maxima of the electron density and minima of the potential. Such electron vortices are called "vortex clumps" in the paper by A. Kabansev et al [START_REF] Kabantsev | Experiments and theory on 2d electron vortex dynamics in sheared flows[END_REF]. These authors show that it is also possible to generate "vortex holes" in pure electron plasmas by removing electrons locally.

The motion of electron vortices and local plasma maxima along the double layer is the result of electron-wave interaction. This electron-wave interaction clearly contributes to electron heating in the double layer since both the electron mean energy and the ionization rate are maximum at the location of the electron vortices. However, as seen in Figure 5.4d and Figure 5.4e electron heating takes place in a region significantly larger than the double layer or the electron vortex size. JP Boeuf and M. Takahashi [START_REF] Boeuf | Rotating spokes, ionization instability, and electron vortices in partially magnetized e × b plasmas[END_REF] have shown that a large part of the electron heating in these conditions is due to cross-field electron transport due to gradB drift. We discuss this in more detail in section 5.3.4.

Parametric study -Influence of the magnetic field profile

The simulations above were performed in the conditions of Ito et al [START_REF] Ito | Self-organization in planar magnetron microdischarge plasmas[END_REF], i.e. for a small magnetron of 2 mm gap length operating in argon at 0.15 torr. These small dimensions were probably convenient for the diagnostics used by Ito et al. but are not usual for practical applications of magnetrons.

Conditions of the simulations -Similarity laws

According to the classical similarity laws of gas discharges, similar results can be obtained for a pressure divided by a scaling parameter s if the discharge dimensions are multiplied by s and the magnetic field is divided by s (same applied voltage). Provided that non-linear reactions such as electron-ion recombination or step-wise ionization are not important, these scaling laws are valid and it is easy to show that the discharge invariants are sx, st, E/s, B/s, n e,i /s 2 , j e,i /s 2 , etc. . . , where x represent a dimension, t the time, n e,i and j e,i the electron and ion densities and current densities (for the same applied voltage). We checked that these similarity laws were satisfied in the simulations.

The results presented in this sub-section correspond to a pressure ten times lower and dimensions ten times larger than in the example of section 5.2.1. The gap length is therefore 2 cm and the argon pressure is 15 mtorr. These conditions are closer to typical conditions of magnetron discharges in applications. To be able to look at situations with a larger number of azimuthal modes, we performed simulations with an azimuthal dimension twice longer than in the example above (i.e. 8 cm at 15 mtorr). We study in this section the effect of the magnetic field, with an axial profile of the magnetic field similar to the simulations above but with different values of the magnetic field on the cathode and anode. The magnetic field distribution is still given by: B(x) = ae - In section 5.2.1, electrons were emitted at the cathode due to ion bombardment, with a given emission coefficient by ion impact noted γ. In this section, instead of using a secondary electron emission by ion impact, we impose a given electron current density, J e0 , uniformly distributed on the cathode surface. The results are very similar to those obtained with secondary emission by ion impact. The advantage of using a given electron current density on the cathode surface is that it is easier to control the plasma density. Electron emission due to ion bombardment introduces a supplementary non-linearity that makes it more difficult to predict the resulting plasma density (we must limit the plasma density in order to avoid having to increase the number of grid points). In other words, the plasma density is very sensitive to the value of the secondary electron emission coefficient. As in the previous section, electrons coming back to the cathode are re-emitted, i.e. the imposed electron current density is a net current density.

Results for moderate magnetic fields

We present below simulations results obtained with increasing magnetic field intensity B(0), from 0.01 T to 0.04 T , with B(d) = 0. The electrons emitted by the cathode are accelerated by the large sheath electric field. In the case B(0) = 10 -2 T the magnetic field is sufficiently small that some high energy electrons are still present in the plasma. The ionization rate is nonnegligible in the plasma (see Figure 5.7) even though the electric field is small and changes sign in the negative glow (Figure 5.8). In other words, ionization is non-local in the negative glow, i.e. is not due to electron heating by the local electric field but results from the presence of high energy electrons accelerated in the cathode sheath. The calculated electron current per unit length at the anode, I ea , is 0.2 A/m, which corresponds to an electron multiplication M = I ea /I e0 of 3.3 between cathode and anode, and to a net secondary electron emission coefficient γ = 1/(M -1) equal to 0.4. This value of the net secondary electron emission coefficient is very large and is unrealistic, meaning that a discharge could not be sustained in these conditions (the actual secondary emission coefficient for argon ions on a metal is much smaller than 0.4) unless the voltage were significantly increased. This example is however useful to illustrate a case where the magnetic field is not large enough to lead to the development of instabilities.

For a magnetic field at the cathode surface equal or larger than B(0) = 0.02 T , azimuthal instabilities develop in the discharge gap. Time averaged electron and ion density profiles, ionization rate, and electric field are also shown on Figure 5.7 and Figure 5.8 for B(0) = 0.02, 0.03, and 0.04 T . The imposed net electron rate for B(0) = 0.03 T . In these conditions, the plasma properties exhibit very regular self-organized structures with a mode m = 3, i.e. an azimuthal periodicity of 2.7 cm. The plasma non-uniformities (or rotating spokes) move in the +E × B azimuthal direction (i.e. downward in Figure 5.9) at a velocity of 20 km/s. For these larger values of the magnetic field, the axial transport of electrons in the plasma is strongly reduced and electrons cannot reach the anode as in a usual, unmagnetized negative glow. A negative (directed toward the cathode) axial electric field must be present in the plasma to draw the electrons to the anode. On the other hand, some ions must also reach the anode to ensure plasma quasineutrality. This situation leads to instabilities. These conditions, where the plasma density gradient (Figure 5.7b, c, d) and electric potential gradient (Figure 5.8) are in opposite directions, are typical to those of the so-called Simon-Hoh instability (A. Simon [START_REF] Simon | Instability of a partially ionized plasma in crossed electric and magnetic fields[END_REF], F. C. Hoh [START_REF] Hoh | Instability of penning-type discharges[END_REF], A. I. Smolyakov et al. [START_REF] Smolyakov | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of discharges[END_REF], J.P. Boeuf [START_REF] Boeuf | Micro instabilities and rotating spokes in the near-anode region of partially magnetized plasmas[END_REF]). We see in Figure 5.8 that the averaged axial electric field in the plasma is directed toward the cathode (in the same direction as the plasma density gradient) and increases when the magnetic field increases. The enhancement of electron transport from cathode to anode due to the instability can be understood by looking at the space distribution of the plasma potential in Figure 5.9c and Figure 5.9d. The electric field is perpendicular to the equipotential contours so that the E × B direction follows the equipotential contours. Since electrons tend to drift in the E × B direction, we see on Figure 5.9b that electrons can move back and forth between cathode and anode along the equipotential lines. When they are close to the anode, the magnetic field is much lower, and some electrons can reach the anode. This is also illustrated on Figure 5.10 where the intensity and lines of electron current density at a given time at steady state are shown for three values of the magnetic field. We see on this figure how the structure in the plasma density and plasma potential allows the electrons to cross the region of large magnetic field between cathode and anode.

As seen above, cross-field electron transport is clearly enhanced by the azimuthal oscillations of the plasma potential and electric field. Collisions also contribute to cross-field electron transport especially for low values of B(0) or in regions of low magnetic field. An important parameter that characterizes the "collisionality" of the plasma, is the Hall parameter, defined as the ratio of the electron gyrofrequency w ce to the electron collision frequency ν eN . This parameter is also equal to the ratio of the E × B drift velocity to the "collisional drift velocity" parallel to the electric field (provided that w ce be sufficiently larger than ν eN : h = wce ν eN : Note that the electron temperature and ionization frequency are smaller in the cathode sheath and larger in the plasma in the B(0) = 0.04 T case compared with the B(0) = 0.01 T case. In the sheath, the larger magnetic field for the B(0) = 0.04 T case limits the electron energy gain. The collisional electron energy gain per unit volume per unit time is J e E where the electron current density J e can be written, neglecting diffusion, as J e = e 2 m ν eN w 2 ce E. For the same electric field, the electron energy gain therefore decreases quickly with increasing magnetic field, which explains the smaller electron mean energy in the cathode sheath for B(0) = 0.04 T . In the plasma, the time averaged electric field is larger for larger magnetic fields (see Figure 5.8) therefore the electron energy gain and electron mean energy increases leading to an increase of the collision frequency, which in turn enhances the collisional electron energy gain. One can therefore say that the electron energy deposition and ionization in the plasma in the 0.04 T case are "more local", i.e. due to the local field, than in the 0.01 T case. It is interesting to note that for much larger magnetic fields (e.g. B(0) = 0.1 T at 15 mtorr, as described below, or B = 1 T at 0.15 torr, as in the previous subsection), ionization in the cathode sheath becomes un-significant while ionization in the plasma is dominant. We can compare for example the ionization rates of Figure 5.4e and Figure 5.9d. Clearly most of the ionization takes place in the plasma in Figure 5.4e, while it is not the case for Figure 5.9d. We can say that the rotating spoke in the conditions of Figure 5.4 is an ionization instability triggered by a Simon-Hoh instability (initial phase of the instability). In the conditions of Figure 5.9 the rotating spoke is the consequence of a plasma non-uniformity resulting from a Simon-Hoh instability but where ionization does not play the dominant role (i.e. the instability could exist even with an azimuthally uniform ionization frequency).

Figure 5.12: Electrons currents I eA1 and I eA2 calculated on two segments of the anode of lengths 80% and 20% respectively of the total anode length in the azimuthal direction, for B(0) = 0.02, 0.03, 0.04 T . Same conditions as above. The total electron anode current I eA = I eA1 +I eA2 is also shown for the B(0) = 0.03 T and B(0) = 0.04 T cases. In all cases the instability moves in the +E×B direction. The number of modes m is indicated (and is not clearly defined for B(0) = 0.04 T )

We can also note that the instability in the conditions of this sub-section, i.e. for magnetic fields B(0) below 0.04 T at 15 mtorr would be difficult to observe experimentally since the ionization rate and therefore the excitation rates are not strongly modulated in space and time. The spokes would be therefore more difficult to detect optically than in the example of section 5.2.2. Another way of detecting plasma rotation is to measure the current on a segmented anode. This has been done in several experiments on Hall thrusters (see, e.g. C. L. Ellison, Y. Raitses, and N. J. Fisch [START_REF] Ellison | Cross-field electron transport induced by a rotating spoke in a cylindrical hall thruster[END_REF].) and in simulations of magnetized plasmas (J. P. Boeuf [START_REF] Boeuf | Rotating structures in low temperature magnetized plasmas-insight from particle simulations[END_REF]). In the simulations, one can for example divide the anode in two segments and calculate separately the current on each. This has been done for the cases reported above (B(0) = 0.02, 0.03, and 0.04 T ) where the anode is divided in two segments of 20% and 80% of the total length of the anode in the azimuthal direction. The results are shown in Figure 5.12.

One can see in Figure 5.12 that the plasma non-uniformity cannot be detected on the segmented anode current for B(0) = 0.02 T . For a magnetic field B(0) = 0.03 T , where the plasma non-uniformity is clearly apparent on the plasma density (see Figure 5.9), the total anode electron current is constant and the current on each segment exhibits small amplitude fluctuations. For a magnetic field B(0) = 0.04 T the instability is not as regular as in the B(0) = 0.03 T case. The electron currents on each of the two segments of the segmented anode present irregular fluctuations. Note that the total current also exhibits some fluctuations in that case.

The spoke velocity can be directly deduced from the motion of the instability in the azimuthal direction. For example, in the case B(0) = 0.03 T of Figure 5.9, the spoke velocity was found to be around 20 km/s. The spoke velocity can also be deduced from the frequency of the current oscillations on the anode segments, provided that the mode number m (i.e. the number of plasma density maxima along the azimuthal direction) is known. For Figure 5.9, we have m = 3, and the periodicity of the spokes is L y /(m -1) = 8/2 cm. From Figure 5.12, the period of the current oscillations is about T = 10/5 ms (5 current maxima over 10 ms). The spoke velocity is therefore on the order of L y /[T (m -1)] = 20 km/s, in agreement with the direct estimation of the spoke velocity.

We present below some results for larger magnetic fields and we will see that the spokes become more and more apparent as the magnetic field increases, and that the direction of rotation changes for B(0) between 0.04 and 0.05 T .

Results for larger magnetic fields

For magnetic fields B(0) equal or above 0.05 T , the motion of the non-uniformity changes direction with respect to the previous cases, and the spokes move in the -E × B direction. The calculated currents on the segmented anode for B(0) = 0.05, 0.07, and 0.1 T (with B(d) = 0) are shown in Figure 5.13. For B(0) = 0.05 T , the plasma non-uniformity is not regular and its motion in the -E × B direction is relatively slow. This can be seen on the oscillation frequency of the current in Figure 5.13a. For B(0) = 0.07 T , the spoke velocity, deduced from the spoke motion, increases to about 20 km/s, and for B(0) = 0.1 T , the spoke velocity is 10 km/s. The anode electron current for the B(0) = 0.07 T and B(0) = 0.1 T cases of Figure 5.13 on each part of the segmented anode present very well defined and regular large amplitude oscillations which seems to be typical of the spoke ionization instability. One long spoke is present for B(0) = 0.07 T , while three smaller size spokes are apparent and well defined for B(0) = 0.1 T . The latter case is close to the case of section 5.2 except for the scaling factor of 10 and for the value of the magnetic field at the anode surface, which is 0 in the present section and was ten times smaller than the cathode magnetic field in section 5.2. As in section 5.2, we note that ionization in the plasma is large in region where electrons flow along the equipotential lines (E×B drift) from the anode region to the cathode region. We will show in section 5.3.4 that the grad B drift, combined with the large azimuthal electric field is responsible for electron heating in that region, and for electron cooling in the region where the E × B drift along the equipotential lines is directed from cathode to anode.

As in the previous sub-section, the electron current density intensity and lines displayed in Figure 5.14d, Figure 5.15d, and Figure 5.16d show how the distortion of the electric potential contributes to cross-field electron transport by creating a kind of "short-circuit" across the magnetic field. To end this section, we show in Table 5.2 the electron current imposed at the cathode, the calculated electron current at the anode, the electron multiplication and the equivalent secondary emission coefficient that can be deduced from the simulations in the conditions of

Electron heating

In this section we discuss the question of electron heating in a magnetron configuration. As we have seen above, when the magnetic field in the cathode region is small enough, electrons are heated in the cathode sheath and can transport their energy relatively deeply in the plasma, where they release their energy through inelastic collisions, excitation and ionization. This is typical of non-magnetized or moderately magnetized conditions such as those of section 5.3.2.

When the cathode magnetic field increases (section 5.3.3), for the same applied voltage, electrons gain less energy in the sheath and do not transport this energy in the plasma. Ionization is "local", i.e. the electron energy gain is balanced by the losses in the cathode sheath. If electron heating in the sheath were the only heating mechanism, the plasma could not be sustained. On the other hand, a very large axial electric field, sufficient to allow electron heating and ionization in the plasma cannot form because such negative electric field combined with a negative plasma density gradient leads to a Simon Hoh instability. It turns out that the azimuthally non-uniform re-distribution of the electric potential and field in the plasma provides a new mean of electron heating that is sufficient to sustain the discharge. Such conditions were met in the examples of section 5.3.3, i.e. for cathode magnetic fields larger than 0.04 T. We noticed in the results above, (see the ionization rate plots of Figure 5.14 -Figure 5.16) that the ionization rate was large in regions (spokes) where the E×B drift along the distorted equipotential lines was directed from anode to cathode. Clearly, this large ionization rate indicates that electrons are heated in this region.

To understand this electron heating, it is useful to look at a specific electron trajectory along an equipotential line. Figure 5.17a shows an example of electron trajectory in the conditions of Figure 5.16. The potential and the electron trajectory are displayed on Figure 5.17a while the electron energy along the trajectory is plotted in Figure 5.17b together with the mean energy along the trajectory assuming that the electron energy gain is due to the ∇B drift. The grad B drift is a consequence of the spatial variation of the magnetic field. It is possible to show (see, e.g. F. F. Chen [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF]) that in the presence of a gradient ∇B in the magnetic field, the electrons experience a drift in the direction perpendicular to B and ∇B. This "∇B drift velocity", is given by:

v ∇B = 1/2ρ e v ⊥ B × ∇B/B 2 (5.3)
where v ⊥ is the electron velocity in the axial-azimuthal plane, ρ e is the electron gyroradius. In our conditions the magnetic field gradient is axial and v ∇B is directed downward in the azimuthal direction (Figure 5.17c). In the presence of a non-zero electric field E , this drift leads to electron heating or cooling depending on the sign of v ∇B • E. This can be written as:

∂ε ⊥ = -v ∇B • E = -ε ⊥ E y /(BL) (5.4)
where ε ⊥ is the electron kinetic energy perpendicular to the magnetic field and L = |B/∂ x B|. Therefore, the electron energy gain per unit time is positive (electron heating) when the electron drifts along the equipotential line from anode to cathode and is negative (electron cooling) when the electron drifts along the equipotential line from cathode to anode (see the arrows in Figure 5.17a and, in Figure 5.17c, the directions of the electric field, E × B drift velocity and ∇B drift velocity in the heating region ).

The electron heating frequency α = |E y /(BL)| can be as large as 10 8 s -1 in the heating region of Figure 5.17 (E y /B = 5 × 10 5 m/s and L=0.5cm) leading to significant heating (the electron transit time in this region is in the 10 ns range). This simple estimation is confirmed numerically on Figure 5.17b, where the red line corresponds to the electron mean energy calculated by integrating the expression above along the electron trajectory: ε ⊥ = ε ⊥0 -v ∇B • Edt . This clearly shows that the ∇B drift plays an important role in electron heating in the conditions of section 5.3.3.

conclusion

We have described in this chapter some instabilities that appear in 2D axial-azimuthal simulations of a magnetron discharge. The properties of the instabilities ("rotating spokes") predicted by the simulations in section 5.2 in the conditions of the experiments of Standford (small, 2 mm gap magnetron operating in argon at 0.15 torr with a magnetic field of 1 T at the cathode surface) are in excellent qualitative agreement with those observed in the experiments. The PIC MCC models satisfies (as expected) the classical discharge scaling laws, as confirmed by the simulations reported in sections 5.3 with a scaling factor of 10 with respect to the Stanford experiments (dimensions 10 times larger, pressure and magnetic field 10 times smaller, i.e. gap length 2 cm, pressure 15 mtorr, and cathode magnetic field 0.1 T ). For relatively low values of the magnetic field, the spoke moves in the +E × B direction. The spoke rotation changes directions for larger values of the magnetic field (rotation in the -E × B direction). In our 15 mtorr argon conditions, this change of direction of motion occurs when the magnetic field at the cathode surface is above 0.05 T .

During the discharge formation a cathode sheath and negative glow form, as in a usual glow discharge. If the magnetic field is large enough, electrons become less mobile than ions in the plasma and an axial electric field tends to form in the plasma to pull the electrons to the anode and the plasma density decreases from the cathode sheath-plasma boundary to the anode. This situation (electric field and plasma density gradient in the same direction) is typical of the Simon-Hoh instability. The development of this instability leads to long wavelength azimuthal oscillations of the plasma density and plasma potential which enhance cross-field electron transport. Electrons flow along the equipotential lines (E × B drift), which oscillate back and forth between the cathode and anode regions.

If the magnetic field at the cathode surface is large enough, electron heating in the cathode sheath region decreases while electron heating in the plasma increases. The combination of large azimuthal electric field in the plasma with the magnetic field gradient leads to a significant electron heating or cooling due to ∇B drift. The spoke forms in the region of electron heating associated with the ∇B drift. The simulations clearly show that the rotating spoke in these conditions is associated with an ionization instability (even though, in the linear stage, the plasma nonuniformity is triggered by a Simon-Hoh instability).

Conditions where rotating spokes are very easily observed in experiments seem to correspond to situations where the magnetic field is large enough so that electron It is interesting to note that the direction of rotation of the spoke can change depending on the conditions. In the conditions of our simulations, the spokes were rotating in the +E×B direction for lower values of the magnetic field at the cathode surface, and in the -E × B direction for higher magnetic field (the transition was around 0.04 -0.05 T in argon at 15 mtorr and for a decay of the magnetic field over 2 cm). In the experiments the transition from -E × B to +E × B rotation is observed when the discharge power or plasma density is increased, for a fixed magnetic field distribution. The -E × B rotation occurs when the ionization region is shifted toward the -E × B side of the maximum plasma density, which occurs, in our conditions, for large enough magnetic field, i.e. the plasma seems to "move" in the direction of the maximum ionization.

Chapter 6 Conclusion

You cannot teach a man anything; you can only help him discover it in himself.

Galileo

In partially magnetized plasmas, electrons are strongly magnetized while ions are not. The plasma pressure in these plasmas is smaller than the magnetic pressure (low-β plasmas), so that the externally applied magnetic field is not modified by the presence of the charged particles. This situation is typical to a variety of plasma sources or ion sources used in practical applications. Often, in these plasmas, a current is driven through the magnetic field by applying an external voltage perpendicular to the magnetic field. This is the case for example in magnetron discharges, Hall thrusters, and negative ion sources for neutral beam injection in fusion. The difference in the magnetization of electrons and ions gives these E × B devices specific properties (for example with respect to instabilities) that are different from those of fusion plasmas.

In this thesis we have studied three specific issues of the physics of partially magnetized plasmas that are relevant to some ongoing work in the GREPHE group of the LAPLACE laboratory. These three topics are:

Plasma confinement by magnetic cusps

Co-extraction of electrons in negative ion sources for fusion Rotating spoke instabilities in magnetron discharges and Hall thrusters CHAPTER 6. CONCLUSION These issues have been studied with the help of particle simulations with Monte Carlo collisions (Particle-In-Cell Monte Carlo Collisions, PIC-MCC simulations). To become familiar with this method, we have studied in detail the PIC MCC algorithm, developed and tested a 1D PIC MCC model and used a 2D PIC MCC model to study ambipolar diffusion. We also studied the accuracy of the methods of integration of the charged particle trajectory in electric fields and magnetic fields and checked the conservation of magnetic moment. We reviewed the question of code parallelization using OpenMP and MPI techniques. In the rest of the thesis we have used 2D and 3D PIC-MCC simulation codes developed in the group to study the three topics of this thesis. We summarized below the main results.

Plasma confinement by magnetic cusps:

Our work on plasma confinement by magnetic cusps is the first attempt at quantifying in a systematic way, with PIC-MCC simulations, the confinement efficiency of magnetic cusps as a function of various parameters such as magnetic field, gas pressure and charged particle temperatures. The confinement efficiency was characterized by the effective leak width of one magnetic cusp deduced from the simulations . We found that the hybrid gyroradius (ρ e ρ i ) 1/2 provides a reasonable order of magnitude if the ion gryroadius is calculated for an ion velocity equal to the Bohm velocity. The early experiments of Hershkowitz et al [START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]. at low pressures were consistent with a leak width of 4 times the hybrid gyroradius (but with a velocity in the ion gyroradius equal to the ion thermal velocity). We found that the leak width is practically independent of the ion temperature T i and scales relatively closely to T 1/2 e with the electron temperature. The hybrid gyroradius only gives an order of magnitude, and the 1/B scaling of the hybrid gyroradius was not satisfied in the calculated leak width. The simulations predict a 1/B dependence of the leak width only at low pressure and low magnetic fields. The variations of the leak width with magnetic field are slower than 1/B at high magnetic fields. Note that the presence of the ion gyroradius in the expression of the leak width is a pure coincidence and has no particular physical meaning since the ions are practically not magnetized in the conditions of most experiments in low-beta plasmas and in our simulations. This is clearly seen in the semi-empirical expression of Bosch and Merlino [START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF], (2 DR/C S ) 1/2 , where the leak width is proportional to the square root of the electron diffusion coefficient perpendicular to the magnetic field divided by the ion Bohm velocity. This expression provides the same scaling as the hybrid gyroradius with magnetic field, ion mass, and electron temperature, but takes in to account the effect of collisions (in the diffusion term). This expression therefore scales as the square root of the gas pressure. In the simulations, the leak width is practically independent of pressure below 0.1 mtorr and increases as p 1/2 with pressure above 0.5 mtorr.

Although we tried to study the leak width over a large range of parameters, more investigations would be useful. For example, more systematic simulations are needed to understand the dependence of the leak width on the ion mass: the simulations predict a dependence of the leak width as M 0.1 with the ion mass, while the semiempirical theories give a M 1/4 dependence. One also needs to better characterize the role of the specific Ramsauer shape of the electron momentum cross-section in argon (it would be useful to perform simulations with constant collision cross-section or constant collision frequency) and to study the dependence of the leak width on the cusp length.

Electron extraction from negative ion sources:

In negative ion sources, negative ions are extracted by applying a positive voltage bias to the extraction grid, placed behind the plasma grid in contact with the plasma. Therefore, electrons tend to be co-extracted with negative ions. It is essential to limit the current of co-extracted electrons, since energetic electrons could damage the grids and should absolutely not reach the accelerator. In negative ion sources for fusion, a "suppression magnetic field" is generated by magnets at each grid aperture to prevent electrons from being extracted from the plasma. The magnetic field around each grid aperture forms a magnetic cusp that is supposed to block electron transport though the grid hole and deviates the electron trajectories to the plasma grid. In practice, the co-extracted electron current is of the same order as the extracted negative ion current. We had noticed in earlier works of the group that the co-extracted electron current in 2D PICC simulations was much smaller than in the experiments. 2D PIC MCC simulations where the cusp magnetic field is in the simulation plane can only describe extraction through an infinite slit and not through a circular or rectangular grid aperture as in a real grid. We checked that 3D PIC MCC simulations give the same results for an infinite slit. In a more realistic, 3D PIC MCC simulation of electron extraction through a rectangular aperture we found that the co-extracted electron current is much larger than in the case of an infinite slit. We proved that E × B electron drift through the grid aperture is responsible for the large co-extracted electron current in a rectangular aperture. For an infinite slit aperture, the E × B drift is parallel to the slit and does not transport electrons to the aperture. In a rectangular (or circular) aperture, the E × B drift is tangent to the equipotential meniscus and pull the electron towards one side of the aperture. We showed that since the cusped magnetic field around the aperture becomes more and more perpendicular to the grid in the corners of the rectangular aperture, electrons are less and less confined along this direction and can be extracted by the large electric field when they approach these corners. Future work in this direction should benefit from the better understanding and CHAPTER 6. CONCLUSION quantification of electron extraction brought by the PIC MCC simulations. For example, the simulations could provide guidance for optimization of the shape and aspect ratio of the grid apertures. Another aspect that should be studied is the role of the presence of negative ions. To simplify the problem we have considered electron extraction without negative ions. Although it seems clear that the mechanisms of electron extraction will be the same in the presence of negative ions, it is necessary to check in future work how the presence of negative ions modifies quantitatively the results.

Rotating spokes in magnetron discharges and Hall thrusters:

Rotating spokes are regions of enhanced light emission rotating in the azimuthal direction of E × B devices. Rotating spokes are ubiquitous in planar magnetron discharges and are often observed in Hall thrusters.

Following a study recently performed in the GREPHE group, we have used a 2D PIC MCC model to explore the physics of rotating spokes in the conditions of the experiments of Standford (small, 2 mm gap magnetron operating in argon at 0.15 torr with a magnetic field of 1 T at the cathode surface) and for different magnetic field distributions. The results predict the formation of a rotating spoke, characterized by the rotation of a region of enhanced ionization. The properties of the simulated rotating spoke are consistent with the experimental observations. The spoke can be considered as an ionization instability that rotates due to the E × B configuration. A very interesting feature of the simulated spoke is that the electron heating in the spoke is the result of ∇B drift of electrons in the distorted electric potential. Electrons undergoing E × B drift along the distorted equipotential lines are heated due to ∇B drift when moving from anode to cathode and are cooled when moving from cathode to anode. Ionization is therefore enhanced on one side of the distorted potential. This result and the role of ∇B drift in electron heating in rotating spokes brings a new insight in the question of "electron energization" in the quasineutral plasma of magnetron discharges, that has been the subject of many speculations in recent papers.

The triggering of the instability observed in the simulations is of the Simon-Hoh type, but this instability evolves into an ionization instability sustained by electron heating due to ∇B drift. Another interesting feature is the fact that the simulated spoke rotates in the -E × B direction in the conditions of the Stanford experiment (as in the experiments), but the direction of rotation can change depending on the discharge parameters. In the experiments on magnetron discharges the direction of rotation changes from -E × B to +E × B when the discharge power is increased. 123 CHAPTER 6. CONCLUSION

In the simulations, we showed that the spoke can rotate in the +E × B direction if the magnetic field is decreased in the cathode region.

These results open up a new direction for diagnostics, simulation and theory of E × B plasmas, rich in potential consequences in applications.

Some issues in the physics and particle simulations of low temperature partially magnetized plasmas and ion sources Résumé

In this thesis we will illustrate some of the issues in the physics and modeling of partially magnetized plasmas with three specific examples that correspond to ongoing studies in the GREPHE group of the LAPLACE laboratory:

1) Electron extraction in negative ion sources for neutral beam injection in fusion 2) Instabilities in magnetron discharges and Hall thrusters 3) Plasma confinement by magnetic cusps -In the study of negative ion sources for fusion, the aim of the GREPHE group is to better understand the physics of the negative ion source, and more specifically, the questions of plasma transport across the magnetic filter and of negative ion extraction from the plasma. One of the important issues in these negative ion sources is to minimize the current of electrons that are co-extracted with the negative ions. In this thesis we focus on this aspect and we try to understand and quantify how electrons can be extracted through a grid aperture when a magnetic cusp is placed in front of the aperture. We discuss, with the help of 3D PIC MCC (Particle-In-Cell Monte Carlo Collisions) simulations, the contributions of different electron drifts (ExB drift, Grad B drift and curvature drift) and instabilities to electron extraction through a grid aperture. -Hall thrusters and magnetron discharges are ExB cylindrical devices with radial magnetic field and axial electric field. It has been known for a long time that instabilities are present in these discharges, leading to important anomalous electron transport. In this thesis we focus on one particular type of instability, called "rotating Spoke", which is known to be present in Hall thrusters and magnetron discharges and is apparent in the experiments as a luminous non-uniformity rotating in the azimuthal direction. In this work we use a 2D PIC MCC simulation to perform a parametric study of this instability. We show that, in some conditions where rotating spokes have been observed in the experiments, Grad B electron drift plays a major role in electron heating and in the formation and maintenance of the rotating spokes. -Magnetic cusps have been used for more than 60 years to confine the plasma in a large variety of conditions. An important parameter characterizing plasma confinement by cusps is the effective loss area in the presence of magnetic cusps. Some semi-empirical theories have been proposed to quantify the effective loss area and their predictions have been compared with numerous experimental results. In spite of these efforts there is no fully reliable expression of the effective wall loss as a function of different parameters such as magnetic field, electron temperature, ion mass, gas pressure, etc… We describe in this thesis an attempt at obtaining scaling laws for the effective loss width of magnetic cusps, based on 2D PIC MCC simulations.

Keyword: low-temperature plasma, partially magnetized, negative ion source, magnetron discharge, magnetic cusp

3

  Plasma confinement by magnetic cusps in low-β plasmas 3.1 Leak width scaling laws and objectives of the chapter . . . . . . . . . 3.2 Plasma generation and particle losses . . . . . . . . . . . . . . . . . . 3.3 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Charging of the dielectric walls . . . . . . . . . . . . . . . . . 3.3.2 magnetic field distribution . . . . . . . . . . . . . . . . . . . . 3.3.3 plasma generation and electron heating . . . . . . . . . . . . . 3.3.4 definition of the leak width . . . . . . . . . . . . . . . . . . . . 3.4 Analysis of the simulation results . . . . . . . . . . . . . . . . . . . . 3.4.1 Space distribution of the plasma properties . . . . . . . . . . . 3.4.2 Scaling laws . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Electron transports across magnetic field in negative ion sources 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1 . 1 :. 3 )D + He 3 =. 4 ). 5 )

 113345 Figure 1.1: Fuel consumption in different energy production methods[1]

Figure 1 . 2 :

 12 Figure 1.2: Schematic of inertial confinement fusion device[2]

Figure 1 . 3 :

 13 Figure 1.3: Basic magnetic mirror configuration[3]

Figure 1 . 5 :

 15 Figure 1.5: ITER tokamak device[6]

Figure 1 . 7 :

 17 Figure 1.7: ITER Neutral beam injector configuration[9]

Fig 1. 8 CHAPTERFigure 1 . 8 :

 818 Figure 1.8: Neutralization efficiency vs energy[START_REF] Barnett | Atomic data for fusion. volume 1: Collisions of h, h2, he and CHAPTER 7. BIBLIOGRAPHY li atoms and ions with atoms and molecules[END_REF]

Figure 1 . 9 :

 19 Figure 1.9: Schematic view of the IPP RF source with one driver

Figure 1 . 11 :

 111 Figure 1.11: Schematic of a planar magnetron discharge[30].

Figure 1 .

 1 Figure 1.12 shows an example of measurements of light emission during a HiPIMS current pulse. The spokes rotated in the +E × B direction during the high current phase while the rotation alternates between the +E × B and -E × B directions during the low current phase.

Figure 1 . 12 :

 112 Figure 1.12: Images of spokes in argon at 2.7 P a illustrating the change of direction in the spoke rotation during a current pulse. (a) and (b) correspond to streak camera images with a 50 ms sweep time. (c)-(e) are images taken at different times with a ICCD camera and a 20 ns exposure time. The spokes rotates in the +E × B direction at a velocities around 5 × 10 3 m/s during the early stage of the pulse, when the current is above approximately 6 A.Later in the pulse, when the current has fallen below 6 A, the spoke velocity is not so well defined and the direction of rotation alternates many times between the +E × B and -E × B directions.[START_REF] Yang | Propagation direction reversal of ionization zones in the transition between high and low current magnetron sputtering[END_REF].

Figure 1 . 13 :

 113 Figure 1.13: Plasma potential distribution (left) 5 mm from the cathode, measured at a given time with a Langmuir probe, and light emission of a rotating spoke in argon at 0.27 P a, dc applied voltage of 270 V , discharge current 100 mA. The spoke is rotating in the -E × B direction. [33].

Figure 1 . 14 :

 114 Figure 1.14: Schematic of a Hall thruster. [34].

Figure 1 . 15 :

 115 Figure 1.15: Color image of a rotating spoke in a wall-less Hall thruster [43].

Figure 1 . 16 :

 116 Figure 1.16: Schematic of a multi-cusped field thruster [44].

Figure 1 . 17 :

 117 Figure 1.17: Magnetic field lines(smooth curves) in some cusped geometries. The symbols x and o represent current into and out of the drawing respectively[47].

Figure 1 . 18 :

 118 Figure 1.18: Schematic diagram of the multipole device[START_REF] Hershkowitz | Plasma leakage through a low-β line cusp[END_REF]

Figure 1 . 19 :

 119 Figure 1.19: Measured electron and ion current profiles in a ring cusp (110 G, argon, 5.5 × 10 -5 torr)[START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF]

  20.

Figure 1 . 20 :

 120 Figure1.20: Leak widths of electron and ions in the ring cusp, as a function of ω pe /ω ce as the plasma density was varied by increasing the filament temperature with a neutral pressure of 7.1 × 10 -5 torr.[START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF]

Figure 1 . 21 :

 121 Figure 1.21: Leak widths as a function of magnetic field at three different pressure in argon (left), and as a function of gas pressure (right) in a ring cusp for a magnetic field of 110 G in the center of the cusp. The hybrid gyrodiameter and the primary electron gyrodiameter are also shown for comparison on the left plot[START_REF] Bosch | Confinement properties of a low-beta discharge in a spindle cusp magnetic field[END_REF]

FigureFigure 2 . 1 :

 21 Figure 2.1: (Color). X-axis coordinate of the gyration of a magnetised electron versus time in normalized units, x/r L and t/τ where r L is the Larmor radius and ω c τ = 2π the orbital period. Magnetic field strength B = 100G, time step ω c ∆t = 0.3 and kinetic energy of 2 eV. The analytical solution (orange dots) is compared to the Boris algorithm with either the exact phase calculation (red solid line), eq. (2.20), or the approximated one (blue solid line) t (ω c ∆t/2)(B/|B|), respectively.

Figure 2

 2 Figure2.2: (Color). Schematic representation of a magnetic bottle (taken from[START_REF] Bittencourt | Fundamentals of Plasma Physics[END_REF]). Trapped electrons will initially follow a given field line and rotate slowly around the magnetic bottle describing a surface where the field line number density remains constant.

Figure 2

 2 Figure 2.3: (Color) Trajectory of a trapped electron inside a magnetic bottle of length 2L = 120 cm. The magnetic field map corresponds to eqs. (2.22)-(2.24). The initial conditions for the trajectory are x 0 = 0, y 0 = 5 cm, z 0 = 15 cm, a velocity angle θ 0 = 28 • and a kinetic energy E k = 2 eV, respectively. We implemented the Boris algorithm, eqs. (2.15)-(2.21), with the approximation t (θ/2) B/|B| and a time step θ = ω c ∆t = 0.3.

Figure 2 . 4 :

 24 Figure 2.4: (Color) Axial component of the magnetic field B z . Profile from eq. (2.24) is shown versus the piecewise fit provided by the array Bz of length n z = 240 and grid size along (Oz) ∆z = 5 mm. L = 60 cm and B 0 = 100G, respectively.

Figure 2

 2 Figure 2.5: (Color) Trajectory of a trapped electron inside a the same magnetic bottle as fig. 2.3 but discretised on a mesh of size ∆x = ∆y = ∆z = 5 mm. The length of the bottle is 2L = 120 cm and B 0 = 100G on axis in its center.

Figure 2

 2 Figure2.6: (Color) z-coordinate of the electron bouncing between mirror points inside the magnetic bottle versus the time step in the Boris algorithm, ω c ∆t = 0.3 and 0.03 respectively. The two calculations are identical. The trajectory is shown for a time corresponding to approximatively 6 × 10 5 gyroperiods τ . The latter is calculated using the maximum magnetic field strength at the exit of the bottle, i.e., 600G.

Figure 2 . 9 :Figure 2

 292 Figure2.9: (Color) Electron and positive ion density profiles over half the length of the simulation domain. The sheath entrance is at the location x s 6.8 mm where the plasma is no longer quasi-neutral (corresponding to n s /n 0 = 0.5). In the configuration of the simulation ν m = 0 and hence k = 0 and h 0.57, respectively.

CHAPTER 3 .

 3 PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMAS

Figure 3 . 1 :

 31 Figure 3.1: Magnetic field distribution in a picket-fence.The picket-fence separates the plasma source from the "target". The plasma source is generated by hot filaments at a negative) potential ( -60 V ) with respect to the grounded chamber walls. Electrons and ions flow from the plasma source to the target through the magnetic cusps formed by the picket-fence. The leak-width, obtained by measuring the profiles of the electron and ion current densities (red line) along the midplane follows a hybrid radius scaling. Hershkowitz and Christensen[START_REF] Christensen | Mass scaling of permanent magnet line cusp plasma leaks[END_REF] showed that the profile of the electron current collected on an equipotential plate placed in the mid-plane along the picket-fence line followed the same hybrid radius scaling.

Figure 3 .

 3 Figure 3.2: (a) Kamaboko negative ion source. Ichikawa et al[78] ; (b) example of gridded ion thruster with magnetic cusps. D. M. Goebel and I. Katz[79].

Figure 3 . 3 :

 33 Figure 3.3: WiPAL multi-dipole magnetic ring cusp plasma source. Cooper et al[52]. The magnet poles are covers with a 3 mm thick insulating, alumina layer.

Figure 3 . 4 :

 34 Figure 3.4: Schematic of a radiofrequeny (RF) negative ion source with multicusp confinement. D. Faircloth and S. Lawrie[80].

Figure 3 . 5 :

 35 Figure 3.5: High Efficiency Multistage Plasma Thruster (HEMPT) of Thales, with magnetic cusps. The magnets are place behind a ceramic (dielectric) wall. N Koch, HP Harmann, and GK Cornfeld[81].

Figure 3 .

 3 Figure 3.6: Simulation domain for the PIC MCC simulation of magnetic cusps. These are line cusps, infinite in the direction perpendicular to the simulation domain, and periodic in the direction parallel to the walls. The walls are made of dielectric materials, i.e. the electron and ion fluxes to the wall must be identical at steady state. The dimensions of the domain in most simulations presented here are L x = 4 cm and L y = 2 cm. The distance between cusps is therefore d = 1 cm. The magnetic field intensity (top) decreases exponentially from the walls. Electrons are heated in the region indicated by the dashed lines. Note that because of the obvious symmetries of the problem, the simulation domain can be divided in two in each direction. This can significantly speed up the simulations.

CHAPTER 3 .Figure 3 . 7 :

 337 Figure 3.7: Cross-section of reaction vs electron temperature. C1: elastic collision (index 1), C2: excitation (index 2), C3: ionization (index 3).

  magnetic field intensity B = B 2 x + B 2 y = B 0 e -πx/d does not depend on y. This magnetic field distribution is shonw in Fig 3.

8

 8

Figure 3 . 8 :

 38 Figure 3.8: Distribution of the magnetic field(intensity and field line) used in the simulations

Figure 3 . 10 :

 310 Figure 3.10: Electric potential profiles as a function of transverse position y for different axial positions x. Same conditions as Fig 3.9: B 0 = 400 G, Argon, 0.1 mtorr, T e = 2 eV , T i = 0.5 eV

Figure 3 . 12 :Figure 3 . 13 :

 312313 Figure 3.12: Profiles of the electron current density to the wall for electrons in different energy ranges. B 0 = 1200 G, 0.1 mtorr, T e = 2 eV , T i = 0.5 eV .

Figure 3 . 16 :

 316 Figure 3.16: Axial profiles of the ion mean energy at two y locations; in the center of a cusp, y = 5 mm and at a mid-position between two cusps, y = 10 mm. B 0 = 400 G, Argon, 0.1 mtorr, T e = 2 eV , T i = 0.5 eV .

Figure 3 . 17 :

 317 Figure 3.17: Leak width normalized to the distance between cusps as a function of magnetic field at the dielectric surface, for different values of the gas pressure. Argon, T e = 2 eV , T i = 0.5 eV . The variations of the normalized hybrid radius (1/B dependence) and a 1/B 2 curve are also shown comparisons. The hybrid radius is calculated for an ion velocity equal to the Bohm velocity.

CHAPTER 3 .

 3 PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMASof the leak width predicted by the expression of Koch and Mathieussent does not appear in the range of pressure considered in Figure3.18. At low pressures, the normalized leak no longer depends on pressure and reaches a constant value, which is the collisionless leak width.

Figure 3 . 18 :

 318 Figure 3.18: Leak width normalized to the distance between cusps as a function of gas pressure for different values of the magnetic field at the dielectric surface. Argon, T e = 2 eV , T i = 0.5 eV . The full black lines correspond to p 1/2 variations.

CHAPTER 3 .Figure 3 . 19 :

 3319 Figure 3.19: Leak width normalized to the distance between cusps as a function of electron temperature, for different combinations of magnetic field at the dielectric surface and gas pressure. Argon, T i = 0.5 eV . The full black line corresponds to T 1/2 e variations.

CHAPTER 3 .

 3 PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMASFinally, we looked at the influence of the ion mass on the leak width. Figure3.20 shows the variations of the normalized leak width with ion mass for different combinations of magnetic field and pressure. It is interesting (and surprising) to note that the leak width does not scale as M 1/4 as suggested by the results of Herskowitz et al. and Bosch and Merlino. The variations of the leak with ion mass is slower than M 1/4 and close to a M 0.1 law in the examples of Figure3.20. We do not have, at the moment a clear explanation of this result.

Figure 3 . 20 :

 320 Figure 3.20: Leak width normalized to the distance between cusps as a function of ion mass, for different combinations of magnetic field at the dielectric surface and gas pressure. T e = 2 eV , T i = 0.5 eV . The full black lines indicate M 1/4 variations (suggested by the results of Hershkowitz et al. and Bosch and Merlino), and M 0.1 variations (best fit to the PIC MCC results in these conditions).

CHAPTER 3 .

 3 PLASMA CONFINEMENT BY MAGNETIC CUSPS IN LOW-β PLASMASof Hershkowitz et al. Surprisingly, the simulations do not reproduce this dependence, and the calculated leak width varies only like M 0.1 with the ion mass.

CHAPTER

Figure 4 . 1 : 1 23%45%Figure 4 . 2 :

 41142 Figure 4.1: Schematic of the simulation domain. On the left-hand-side (LHS) is the plasma volume, which is separated from a biased electrode (EG) by a grid (PG) with two slit apertures. The volume between the grid and the biased electrode is free of plasma. The latter is generated numerically by the injection of electron-positive ion pairs in the gray area labelled "plasma source". We fixed the plasma density in the model and hence each positive ion lost on a physical boundary (Dirichlet) induce the re-injection of an electron-ion pair. The magnetic cusp field profile (in yellow) is shown of the right-hand-side (RHS) of the figure. In practice, this field would be generated by permanent magnet bars of reversed polarities (blue and red colours, N stands for "north" and S for "south", respectively). The boundary condition on the LHS together with the plasma grid (PG) and extraction grid (EG) are of Dirichlet type while the top/bottom and left/right boundaries are periodic in 3D.

4. 3

 3 Transport of electrons across a magnetized aperture 4.3.1 2D versus 3D PIC-MCC calculation of a slit aperture of infinite length

Figure 4 .Figure 4 . 3 :

 443 Figure 4.3: Electron current collected on the extraction grid (EG). 2D versus 3D PIC-MCC calculation for a configuration with a slit aperture of infinite length.

Figure 4 . 4 :

 44 Figure 4.4: Plasma potential profile along (Ox) for a position on the plasma grid (PG) corresponding to the middle of the interstice between two slit apertures. Estimates from the 2D PIC-MCC model (black solid line) are compared to the 3D PIC-MCC calculation for either a configuration with two slit apertures (blue solid line), as in fig. 4.1, or for half the simulation domain along (Oz), i.e., only one cusp (red solid line), respectively.

Figure 4 . 5 :

 45 Figure 4.5: 2D Electron flux profile. Four planes are displayed, that is, (i) the place encompassing the cusp magnetic field line (XZ at Y = 2.6 cm), (ii) the electron E × B drift plane (XY, at Z = 0.6 cm, i.e., in the middle of the slit aperture), (iii) the front of the PG and, lastly, (iv) the EG, respectively. Γ max = 3 × 10 21 m -2 • s -1 .

Figure 4 . 6 :

 46 Figure 4.6: Examples of electron trajectories, (a) in the Z=0.6 cm plane, (b) for electrons moving along a magnetic field line toward the upper edge of the PG (see Fig.5) while drifting in the E × B direction toward the right edge of the PG (the B z component of the magnetic field, perpendicular to the XY plane, decreases during this motion). For simplicity, these trajectories are calculated assuming that the XY distribution of the potential does not vary while the electrons moves along the magnetic field line (i.e. the electric potential does not depend on Z).

Figure 4 . 7 :

 47 Figure 4.7: Ratio of the electron current collected on the EG divided by the total current flowing toward the PG, I EG / (I EG + I P G ) (in percent) versus time. I EG (I P G ) is the current impacting the EG (PG), respectively.

CHAPTER 4 .Figure 4 . 8 :

 448 Figure 4.8: Electron flux profile in the plane Z = 0.6 cm of fig. 4.5 for different time shots labeled from A to I (which are also indicated in fig. 4.7 for reference). The isopotential lines for φ = 10V and 20V inside the plasma volume are also shown. Γ max = 2.3 × 10 21 m -2 • s -1 .

Figure 4 . 9 :

 49 Figure 4.9: Plasma potential profile in the slit aperture mid-plane (Z = 0.6 cm) for different time shots.

Figure 5 . 1 :

 51 Figure 5.1: Schematic of the Stanford experiment (Ito et al[98])

Figure 5 . 2 :

 52 Figure 5.2: Sequences of images taken with a high-speed camera in a 150 mtorr argon plasma through the anode (transparent ITO anode) of the magnetron of Fig 5.1

  The axial distribution of the magnetic field is shown in Figure5.3b and corresponds to the expression: B(x) = ae -x 2 2a 2 + b (5.1) with σ = 0.35d where d = 2 mm is the gap length. The constants a and b are calculated so that B(0) = 1 T , B(d) = 0.1 T .

Figure 5 . 3 :

 53 Figure 5.3: (a) Simulation domain for the model of the experiment of Ito et al.. The simulations are performed in argon at a pressure of 20 P a (0.15 torr) and the applied voltage is 260 V . Periodic boundary conditions are used in the azimuthal direction. The net secondary electron emission due to ion impact is γ = 0.005 ("net" means that electrons coming back to the cathode are re-emitted). The magnetic field is perpendicular to the simulation domain and its direction is indicated on the figure; (b) axial distribution of the magnetic field.

Figure 5 .

 5 Figure 5.4 shows the space distribution at a given time, of the ion density, electron density, electric potential, electron mean energy, ionization rate, and ion mean energy. The results of Figure 5.4 can be summarized as follows:

Figure 5 . 4 :

 54 Figure 5.4: Colour contours at a given time, of (a) ion density, (b) electron density, (c) electric potential, (d) electron mean energy, (e) ionization rate, (f ) mean ion energy, in the conditions of the experiments of Ito et al. (see text). The maximum values are: 4 × 10 17 m -3 for the electron and ion densities, 269 V for the electric potential, 15 eV for the electron mean energy, 1.3 × 10 25 m -3 s -1 for the ionization rate, and 30eV for the ion mean energy. Contour lines of ionization rate are superimposed on the ion density plot (c), extra contour lines at 250, 260, 262, 264, 268 V are shown on the electric potential (c) and on the electron mean energy (d), two equipotential contours at 255 and 265 V are superimposed on the ionization rate (e) and ion mean energy (f ) plots (these two lines approximately define the double layer -see text).

Figure 5 .

 5 Figure 5.4c and Figure 5.4d) . These maxima are associated with "holes" in the electron density. These electron holes and the associated potential maxima move along the equipotential line separating the two regions. Note that since the electric field is perpendicular to the equipotential lines, the E × B electron drift follows the equipotential lines.

Figure 5 . 5 :

 55 Figure 5.5: Axial profiles of the electron density, ion density, electric potential, and electron mean energy in the conditions of Fig 5.4, at the azimuthal position y = 2 mm. The units for the densities, potential and mean energy are respectively: 5 × 10 17 m -3 , 300 V , 18 eV .

  σ = 0.35d where d = 2 mm is the gap length. The constants a and b are calculated so that B(0) = 1 T , B(d) = 0.1 T (same as above, with a scaling of 10) to 10 -2 T , with B(d) = 0, 1 × 10 -2 T , 2 × 10 -2 T , etc. . . . The simulation domain and some of the magnetic field profiles used in this section are shown in Figure5.6.

Figure 5 . 6 :

 56 Figure 5.6: Simulation domain (a), and magnetic field profiles (b) with B(d) = 0 and B(0) = 0.1, 0.05, 0.02, and 0.01 T , for the simulations of section 5.3.3. The simulations were performed on a 128x512 grid with about 200 particles per cell.

Figure 5 . 7 :

 57 Figure 5.7: Axial profiles of the time averaged electron density, ion density, and ionization rate, for B(0) = 10 -2 , 2×10 -2 , 3×10 -2 , and 4×10 -2 T and B(d) = 0. Argon, 15 mtorr, 260 V . The imposed net electron current at the cathode is 0.06, 0.02, 0.02, and 0.01 A/m for B(0) = 10 -2 , 2 × 10 -2 , 3 × 10 -2 , and 4 × 10 -2 T respectively.

Figure 5 .

 5 Figure 5.7 displays the axial profiles of the time averaged electron density, ion density, and ionization rates for four values of the magnetic field at the cathode surface. For B(0) = 0.01 T and a net electron current density at the cathode 7.5 × 10 -3 A/m2 (i.e. electron current per unit length equal to 6 × 10 -2 A/m) instabilities do not form and the plasma properties are uniform in the azimuthal direction. The discharge in these conditions looks very much like a regular glow discharge, with a cathode sheath and negative glow. In these conditions of relatively low magnetic field, electrons are more mobile than ions in the negative glow, and the plasma diffuses ambipolarly to the anode. The axial profiles of the axial electric field are shown in Figure5.8.

Figure 5 . 8 :

 58 Figure 5.8: Axial profile of the electric field in the same conditions as Fig 5.7, for B(0) = 10 -2 , 2 × 10 -2 , 3 × 10 -2 , and 4 × 10 -2 T and B(d) = 0. Argon, 15 mtorr, 260 V .

Figure 5 . 10 :

 510 Figure 5.10: Lines (cyan) and intensity (grey scale) of the electron current density for B(0) = 0.02, 0.03, and 0.04 T . Lighter grey corresponds to Argon, 15 mtorr, 260 V .

Figure 5 .

 5 Figure 5.11a, b, c show respectively the axial profiles of the time averaged Hall parameter, electron temperature and electron-neutral collision frequency for B(0) = 10 -2 T , and B(0) = 4 × 10 -2 T . It is interesting to note that, contrary to what could be expected, the Hall parameter is not maximum in the cathode region, where the magnetic field is the largest. This is because of the important Ramsauer minimum in the electron-neutral cross-section in argon. The electron-neutral momentum cross-section increases by almost two orders of magnitudes when the electron energy increases from 0.3 eV to 10 eV . Because of the large electric field in the cathode sheath, the electron mean energy is large in this region (the maximum electron mean energy is more than 100 eV for B(0) = 0.01 T , and around 50 eV for B(0) = 0.04 T , see Figure 5.11b), and the electron-neutral collision frequency is much larger than in the plasma where the electron mean energy drops to much lower values. Because of the Ramsauer minimum, the collision frequency strongly decreases in the plasma (Figure 5.11c). Therefore the Hall parameter is relatively small in the cathode region (on the order of 10 for B(0) = 0.01 T , and of 50 for B(0) = 0.04 T ), increases sharply at the plasma entrance because of the electron

Figure 5 . 11 :

 511 Figure 5.11: Time averaged axial profiles of the, (a), Hall parameter, (b), electron mean energy, and (c), electron-neutral collision frequency and ionization frequency for B(0) = 10 -2 T and B(0) = 4 × 10 -2 T . Same conditions as Figure 5.7, Figure 5.8.

Figure 5 . 13 :

 513 Figure 5.13: Electrons currents IeA1 and IeA2 calculated on two segments of the anode of lengths 80% and 20% respectively of the total anode length in the azimuthal direction, for B(0) = 0.05, 0.07, 0.1 T . Same conditions as above. The total electron anode current I eA = I eA1 +I eA2 is also shown. In all cases the instability moves in the -E × B direction. The number of modes m is indicated (and is not clearly defined for B(0) = 0.05 T .

Figure 5 . 14 :

 514 Figure 5.14: Space distribution of (a), ion density, (b), electric potential, (c), ionization rate, (d), electron current density (lighter grey corresponds to larger current density, current lines in cyan). The min-max values are: [0 -3.6 × 10 15 ] m -3 for n i , [0 -262] V for φ, [0 -1 × 10 22 ) m -3 s -1 for n e ν i . Argon, 15 mtorr, B(0) = 0.05 T , B(d) = 0.

Figure 5 . 15 :

 515 Figure 5.15: Space distribution of (a), ion density, (b), electric potential, (c), ionization rate, (d), electron current density (lighter grey corresponds to larger current density, current lines in cyan). The min-max values are: [0 -5.8 × 10 15 ] m -3 for n i , [0 -262] V for φ, [0 -1.5 × 10 22 ) m -3 s -1 for n e ν i . Argon, 15 mtorr, B(0) = 0.07 T , B(d) = 0.

Figure 5 . 14 ,

 514 Figure 5.14, Figure 5.15, and Figure 5.[START_REF] Fantz | A comparison of hydrogen and deuterium plasmas in the ipp prototype ion source for fusion[END_REF] show 2D distributions of the plasma properties for cases B(0) = 0.05 T , 0.07, and 0.1 T respectively. As the magnetic field at the cathode surface increases, ionization in the cathode region or at the sheath-plasma boundary decreases while ionization in the bulk plasma increases. One long spoke is present for B(0) = 0.07 T , while three smaller size spokes are apparent and well defined for B(0) = 0.1 T . The latter case is close to the case of section 5.2 except for the scaling factor of 10 and for the value of the magnetic field at the anode surface, which is 0 in the present section and was ten times smaller than the cathode magnetic field in section 5.2. As in section 5.2, we note that ionization in the plasma is large in region where electrons flow along the equipotential lines (E×B drift) from the anode region to the cathode region. We will show in section 5.3.4 that the grad B drift, combined with the large azimuthal electric field is responsible for electron heating in that region, and for electron cooling in the region where the E × B drift along the equipotential lines is directed from cathode to anode.

Figure 5 . 16 :

 516 Figure 5.16: Space distribution of (a), ion density, (b), electric potential, (c), ionization rate, (d), electron current density (lighter grey corresponds to larger current density, current lines in cyan). The min-max values are: [0 -5.8 × 10 15 ] m -3 for n i , [0 -262] V for φ, [0 -2 × 10 22 ) m -3 s -1 for n e ν i . Argon, 15 mtorr, B(0) = 0.1 T , B(d) = 0.

Figure 5 . 14 ,Table 5 . 2 :

 51452 Figure 5.15, and Figure 5.16. CHAPTER 5. IONIZATION INSTABILITY IN MAGNETRON DISCHARGES B(0) I e0 (A/m) I ea (A/m) M = I ea /I e0 γ = 1/(M -Imposed net electron current emitted at the cathode, I e0 , calculated anode electron currentI ea , electron multiplication M , and equivalent net secondary electron coefficient, γ for different values of the cathode magnetic field B(0) (conditions of Figure 5.14, Figure 5.15, and Figure 5.16).

Figure 5 .

 5 Figure 5.17: (a) Equipotential contours (grey color and black lines) and example of electron trajectory (blue line); (b) Electron energy along the trajectory of (a) (blue line), and averaged electron energy deduced from the grad B heating (red line, see text); (c) diagram showing the directions of the electric field E, the E × B drift velocity and the ∇B drift velocity on a given point of the electron trajectory. Conditions of Figure 5.16: Argon, 15 mtorr, B(0) = 0.1 T , B(d) = 0.

Figure 5 . 18 :

 518 Figure 5.18: Distribution of (a) electron mean energy, and (b) electric potential at a given time at steady state in the conditions of Figure 5.16 and Figure 5.17: Argon, 15 mtorr, B(0) = 0.1 T , B(d) = 0. The min, max are [0 -15] eV for the mean energy, and [0 -265] V for the potential.

Figure 5 .

 5 Figure 5.18 displays the space distribution of the electron mean energy together with the electric potential at a given time of a steady state discharge with B(0) = 0.1 T . This figure confirms the arguments above concerning the electron heating and cooling regions. The spokes are located in the electron heating regions, as seen in Figure 5.16.

CHAPTER 5 .

 5 IONIZATION INSTABILITY IN MAGNETRON DISCHARGESheating in the cathode region is small and most of the electron heating takes place in the plasma due to the ∇B drift.

  CONTENTS 2.1.2.3 Error introduced by the discretisation of the magnetic field map . . . . . . . . . . . . . . . . . . . . . 2.1.2.4 Hybrid OpenMP and MPI parallelisation of the particle pusher . . . . . . . . . . . . . . . . . . . . . . . 2.1.2.5 Parallelization of the Poisson solver . . . . . . . . . . 2.1.2.6 Implementation of collisions in a particle model -MC and DSMC methods . . . . . . . . . . . . . . . . . . 2.1.3 Elastic and inelastic collision processes . . . . . . . . . . . . . 2.2 1D analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3 2D PIC-MCC simulation of ambipolar diffusion in a plasma . . . . .

  INTRODUCTIONthe total numbers of particles in plasma is enormous, so one particle in the simulation (test particle, super particle, or macro particles) represents a large number of real particles. And the model solves Newton's equation coupled with Maxwell's equation(or just Poisson's equation in electrostatic PIC simulation) to study the detail of particle behavior(trajectory). The detail of the PIC model will be discussed in the following chapter. A Monte Carlo collision module is added to PIC models to describe the effect of collisions (PIC-MCC simulations)

3. Particle model or Particle in cell model(PIC) describe the trajectories of individual plasma particles(ions electrons and neutrals) interacting with other plasma particles, neutral atoms or molecules and boundary materials. Usually CHAPTER 1.

  Figure 2.7: (Color) Execution time of the particle pusher (per time step) normalized to the number of macroparticles in the simulation versus the number of cores. The time is shown either with (red and grey lines) or without implementing a sorting algorithm (black-line).We use 80 particles-per-cell (ppc), a numerical resolution of 96×64×128 grid nodes (black and red lines) and 192 × 128 × 256 (grey line). The calculation is performed with a 3D PIC-MCC model on a 10 cores Intel Xeon processor E5-2680 v2 (25M cache, 2.

	w/o sorting, 7.9•10 8 nodes sorting, 7.9•10 8 nodes 6.3•10 9 nodes

  Execution time of the geometric multigrid Poisson solver (per time step) normalized to the number of grid nodes in the simulation versus the number of cores. the numerical resolution is 512 3 (black line), 1024 3 (red) and 2048 3 (grey) grid nodes, respectively. The calculation is performed on a 10 cores Intel Xeon processor E5-2680 v2 (25M cache, 2.80 GHz). We set the number of OpenMP threads to 10.

	10 10 ∆t (s/nodes)	-7 -6		512 3 nodes 1024 3 nodes 2048 3 nodes
	10	-8
		1	10	100	1000
			# of cores
	Figure 2.8: (Color)

Table 2 .

 2 1: 2D PIC-MCC calculation versus analytical theory. We used the electron temperature derived from the simulation, T e 5.2 eV, for the comparison.

		Model	Analytical
	φ s (V)	14.6 (x s 6.8 mm) 14.75 [eq. (2.87)]
	(φ 0 -φ s )(V)	3.4	3.6 [eq. (2.81)]
	i (eV)	17.2	17.3 [eq. (2.88)]
	n s /n 0	0.5	0.5 [eq. (2.80)]

Table 3 .

 3

1: Reaction used in simulation

Table 4 .

 4

1: Particle processes.

A. Later in the pulse, when the current has fallen below 6 A, the spoke velocity is not so well defined and the direction of rotation alternates many times between the +E × B and -E × B directions.[START_REF] Yang | Propagation direction reversal of ionization zones in the transition between high and low current magnetron sputtering[END_REF].

This assumption and the following collision model was derived by G. Hagelaar.

PLASMAS

wall and repel electrons, as in a usual, un-magnetized sheath. The space distribution of the electron temperature is displayed in Fig 3 .15. We see that the electron temperature is well maintained at 2 eV in a large part of the plasma region up to the cusp region, except between cusps, in the region just current emitted at the cathode, anode electron current, electron multiplication and equivalent net secondary emission coefficient are shown in Table 5.1. As expected, the electron multiplication in the gap significantly increases when the magnetic field increases due to a better confinement of the electrons.