Keywords: fluid modeling, plasma, instabilities, ExB drift, anomalous transport, Hall thruster

It is well known from experiments that magnetized low-temperature plasmas in devices such as Hall thrusters and ion sources often show the emergence of instabilities that can cause anomalous transport phenomena and strongly affect the device operation. In this thesis we investigate the possibilities to simulate these instabilities self-consistently by fluid modeling. This is of great potential interest for engineering. We used a quasineutral fluid code developed at the LAPLACE laboratory, called MAGNIS (MAGnetized Ion Source), solving a set of fluid equations for electrons and ions in a 2D domain perpendicular to the magnetic field lines. It was found that in many cases of practical interest, MAGNIS simulations show plasma instabilities and fluctuations. A first goal of this thesis is to understand the origin of the instabilities observed in MAGNIS and make sure that they are a physical result and not numerical artifacts. For this purpose, we carried out a detailed linear stability analysis based on dispersion relations, from which analytical growth rates and frequencies were successfully compared with those measured in MAGNIS simulations for simple configurations forced to remain in a linear regime. We then identified these linear unstable modes and their responsible mechanisms (involving parameters such as the density gradient, electric and magnetic fields and inertia), known from the literature, that are likely to occur in these fluid simulations. Subsequently, we simulated the nonlinear evolution and saturation of the instabilities and quantified the anomalous transport generated in different cases relevant to ion sources, depending on various key parameters of the system (electric and magnetic fields and electron temperature). Finally, we highlighted several limitations of MAGNIS, and more generally of fluid models, due to the physical approximations made (quasineutrality, absence of kinetic effects). We showed that the fluid modes are sometimes most unstable at infinitely small scales for which the theory is no longer valid and which cannot be resolved numerically. We proposed, and tested in MAGNIS, ways to overcome this problem by introducing effective diffusion terms representing small scale processes (non-neutrality, Larmor radius).

Résumé

Il est bien connu que les plasmas froids magnétisés dans des dispositifs tels que les propulseurs Hall et les sources d'ions montrent souvent l'émergence d'instabilités qui peuvent provoquer des phénomènes de transport anormaux et affecter fortement le fonctionnement du dispositif. Dans cette thèse, nous étudions les possibilités de simuler ces instabilités de manière auto-cohérente par la modélisation fluide. Cela n'a jamais été fait auparavant pour ces conditions de plasma froid, mais cela présente un grand intérêt potentiel pour l'ingénierie. Nous avons utilisé un code fluide quasi-neutre développé au laboratoire LAPLACE, appelé MAGNIS (MAGnetized Ion Source), qui résout un ensemble d'équations fluides pour les électrons et les ions dans un domaine 2D perpendiculaire au champ magnétique. On a constaté que dans de nombreux cas d'intérêt pratique, les simulations MAGNIS produisent des instabilités et des fluctuations du plasma. Un premier objectif de cette thèse est de comprendre l'origine de ces instabilités observées dans MAGNIS et de s'assurer qu'elles sont un résultat physique et non un artefact numérique. Pour ce faire, nous avons effectué une analyse de stabilité linéaire basée sur des relations de dispersion, dont les taux de croissance et les fréquences qui en sont issus analyse ont été comparés avec succès à ceux mesurés dans les simulations de MAGNIS pour des configurations simples et forcés à rester dans un régime linéaire. Nous avons ensuite identifié les principaux modes et mécanismes de ces instabilités (induits par les champs électrique et magnétique, le gradient de densité et l'inertie), connus de la littérature, susceptibles de se produire dans ces simulations de fluides. Par la suite, nous avons simulé l'évolution non-linéaire et la saturation des instabilités et quantifié le transport anormal généré dans différents cas relatifs aux sources d'ions en fonction de divers paramètres clés du système (champs électriques et magnétiques et température des électrons). Enfin, nous avons mis en évidence plusieurs limitations de MAGNIS, et plus généralement de modèles fluides, dues aux approximations physiques (quasi-neutralité, absence d'effets cinétiques). Nous avons montré que les modes fluides sont parfois les plus instables à des échelles infiniment petites où la théorie n'est plus valable et ne peuvent donc être résolues numériquement. Nous avons proposé différentes manières de remédier à ce problème par l'introduction de termes diffusifs inspirés de la physique à petite échelle (non-neutralité, rayon de Larmor), que nous avons ensuite testés dans MAGNIS. mots clés : modélisation fluide, plasma, instabilités, dérive ExB, transport anormal, propulseur de Hall

Introduction Plasma devices and applications

Plasma, a quasi-neutral ionized gas at high temperature, first described by Langmuir in 1927, is defined as one of the four fundamental states of matter. Plasma exists in a natural state in space (stars, the sun, solar wind, Earth magnetosphere and ionosphere) and is believed to represent most of the visible matter in the universe. On Earth, plasmas are artificially generated, typically by application of strong electric current and/or electromagnetic radiation, for multiple laboratory and industrial use, such as for thermonuclear fusion, material processing and manufacturing, medical and environmental applications, and many others.

Many configurations of the plasma devices with different attached parameters exist depending on the application field. Depending on emphasis, one may classify different particular type of plasma according to specific applications and the plasma characteristics; for example, plasma regimes can be defined as high or low-temperature plasmas, magnetized or non-magnetized, fully or partially ionized, low or high (atmospheric) pressure, and so on. The existing applications attached to these regimes range from fusion reactors to plasma etchers, along with plasma arc jets for cutting and welding, plasma sources and reactors for sterilization, water and exhaust cleaning and many others. The fully ionized plasma is confined in a magnetic "cage" formed by poloidal and toroidal magnetic field to keep it away from the walls, and thus, to prevent it from cooling and damaging the chamber. Some of the largest plasma devices, tokamaks and stellarators, use strong magnetic field (1-3 T) to confine high temperature plasmas (10-20 KeV) for thermonuclear fusion applications (we show an example of a tokamak device figure 4). On smaller scales, there are various plasma sources extensively used in industry for material processing and film depositions such as capacitive and inductive discharges. Some of these plasma devices also resort to apply magnetic fields to provide and/or improve the confinement, or for other purposes; this results in the magnetic field being an additional control parameter of the system. Magnetically enhanced discharges typically, often use crossed electric and moderate magnetic fields configuration (E × B), to confine electrons and accelerate ions to produce thrust, bombard, material deposit on the surface, and other applications of ion and plasma beams. Such magnetized low-temperature plasma devices, shown in figures (1), ( 2) and (3), share much common physics. In this thesis, we focus on this specific configuration as used in magnetic filters and closed-drift plasma accelerator devices, e.g. ion sources for the neutral injection and the Hall Thruster for space propulsion.

Plasma confinement, instabilities and transport

The role of the magnetic field in plasma devices is to confine the plasma (both electrons and ions, or only electrons), thus insulating it from the walls and reduce energy losses so that the plasma can be heated to achieve high temperature, e.g. such as that required for thermonuclear reaction for fusion applications. For low temperature plasma systems, the magnetic field can be used to confine electrons and maintain sufficient electron temperature required to ionize a neutral gas (note that low temperature still means around 10 eV, corresponding to 10 5 K).

Already in early experiments, it was observed that the plasma displays a very "noisy" behaviour and exhibits a wide range of fluctuations of density and electric field ( [START_REF] Bohm | The characteristics of electrical discharges in magnetic fields[END_REF][START_REF] Hoh | Instability of pening-type discharge[END_REF], also see the references in [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF]). It was also realized that the plasma transport under these conditions strongly exceeds the classical diffusion values that would be expected from inter-particle collisions. This enhancement was attributed to instabilities in the plasma (plasma turbulence) and was called turbulent diffusion, also "Bohm" diffusion, named after D. Bohm who first highlighted and studied this phenomena [START_REF] Bohm | The characteristics of electrical discharges in magnetic fields[END_REF].

In general, instabilities may be viewed as a tendency of the system to transit towards a lower energy state; in this specific context, an important phenomenon occurring, formulated by the le Chatelier principle, is that the system reacts such as to counter the imposed perturbation in order to return back to its equilibrium. This action together with the system inertia, provides the restoring effect resulting in oscillations, and thus periodic motion and waves. As a result, the instabilities are closely related to the wave phenomenon, i.e. sound waves that exist in neutral gas. As any confined plasma and plasmas with flows are, by definition, away from the most equilibrium state, these systems are then trying via the instabilities to move towards the state of the thermodynamic equilibrium, which would be the state with the lowest possible energy.

In magnetized plasmas, there are many wave eigen-modes that facilitate the appearance of instabilities. For instance, one of the most violent and dramatic plasma instabilities were observed by fusion physicists in late 1950's and 1960's, while leading experiments aiming to trap a thermonuclear plasma with a magnetic field. Typical examples of instabilities found in such configurations are the Kink instability and the Sausage instability [START_REF] De | Mhd instabilities in tokamaks[END_REF] (seen in figure 5). These instabilities were called MagnetoHydroDynamic (MHD) modes since they could be explained on a basis of magnetohydrodynamics, a theory developed by H. Alfvén in 1942 to describe the behaviour of electrically conducting fluids under electromagnetic field. The specific MHD waves, which exist only due to the magnetic field, propagate with Alfven velocity, which is typically relatively fast compared to the sound velocity. As a result, an essential feature of MHD modes is their fast time scale which also means that they have large length scales. Further studies of magnetically confined plasmas in tokamaks (fig. 4), other confinement systems, and space plasma physics have revealed numerous instabilities [START_REF] Mikhaïlovskii | Instabilities in a confined plasma[END_REF] in high temperature magnetized plasmas. Since 1960's, instabilities in fusion grade plasmas has been a topic of active theoretical, experimental, and computational studies culminating in pretty mature understanding of the physics of large scale fusion devices such as ITER [START_REF] Ikeda | Progress in the ITER physics basis[END_REF].

While anomalous electron transport across the magnetic barrier was discovered and pointed out first in experiments by D. Bohm [START_REF] Bohm | The characteristics of electrical discharges in magnetic fields[END_REF] and confirmed in later experimental and theoretical studies of low temperature magnetized plasma devices [START_REF] Spitzer | Particle diffusion across a magnetic field[END_REF][START_REF] Hoh | Instability of pening-type discharge[END_REF][START_REF] Simon | Instability of a partially ionized plasma in crossed electric and magnetic field[END_REF][START_REF] Tilinin | Plasma oscillations in closed-drift accelerators with an extended acceleration zone[END_REF][START_REF] Sakawa | Excitation of the modified simon-hoh instability in an electron-beam produced plasma[END_REF][START_REF] Sakawa | Nonlinear evolution of the modified simon-hoh instability via a cascade of side-band instabilities in a weak beam-plasma system[END_REF], plasma turbulence in low temperature plasmas for industrial applications are not well understood. The properties of a magnetized low-temperature plasma, such as present in Hall Thrusters and similar closed-drift plasma accelerators technologies (introduced later in the late 1960's and 1970's), were found to be different in many ways than those of fusion or space plasma due to the presence and close proximity of walls, low ionization and only partial magnetization of ions. These specific features bring in new complexities such as the sheath theory, losses of particles, ionization and recombination, and other effects. Very early studies of Hall thrusters in the Soviet Union have already revealed ubiquitous presence of waves and instabilities and their likely important role in the performance of these devices [START_REF] Esipchuk | Drift instability in a hall-current plasma accelerator[END_REF][START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF][START_REF] Morozov | Reviews of plasma physics[END_REF][START_REF] Arefev | Instability of a current-carrying homogeneous plasma[END_REF]. Subsequent introduction of this technology in the West have further stimulated interest in fluctuations and transport in Hall thrusters [START_REF] Bouchoule | Physical investigations and developments of hall plasma thrusters[END_REF][START_REF] Choueiri | Anomalous resistivity and heating in current-driven plasma thrusters[END_REF][START_REF] Cappelli | The hall effect and rocket flight[END_REF][START_REF] Chesta | Stability of a magnetized hall plasma discharge[END_REF][START_REF] Boeuf | Low frequency oscillations in a stationary plasma thruster[END_REF][START_REF] Escobar | Low frequency azimuthal stability of the ionization region of the hall thruster discharge. i. local analysis[END_REF][START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF][START_REF] Hagelaar | Role of anomalous electron transport in a stationary plasma thruster simulation[END_REF].

The interest in these problems has recently renewed again due to the need for re-scaling of Hall thrusters to larger and lower powers, along with the necessity to understand new emerging devices with different plasma conditions such as neutral beams sources, magnetic filters [START_REF] Boeuf | Physics of a magnetic filter for negative ion sources. ii. e x b drift through the filter in a real geometry[END_REF], and other variations of electric propulsion system, e.g. employing magnetic nozzle [START_REF] Merino | Fully magnetized plasma flow in a magnetic nozzle[END_REF].

Recent advances in diagnostic capabilities [START_REF] Tsikata | Three-dimensional structure of electron density fluctuations in the hall thruster plasma: ExB mode[END_REF][START_REF] Tsikata | Collective light scattering for the study of fluctuations in magnetized plasmas: The hall thruster case study[END_REF][START_REF] Tsikata | Modulated electron cyclotron drift instability in a high-power pulsed magnetron discharge[END_REF], have allowed non-invasive detection of small scale fluctuations thought to be one of the mechanism of the anomalous transport. The kinetic simulations, in particular Particle-In-Cell (PIC) model, have revealed the presence of kinetic instabilities [START_REF] Garrigues | Exb electron drift instability in hall thrusters : Particle-in-cell simulations vs. theory[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in hall effect thrusters. ii. kinetic model[END_REF][START_REF] Lafleur | Anomalous electron transport in hall-effect thrusters: Comparison between quasi-linear kinetic theory and particle-in-cell simulations[END_REF], fig. 6, for the conditions of Hall thruster devices. Analytical theory [START_REF] Rc Davidson | Anomalous transport in high-temperature plasmas with applications to solenoidal fusion systems[END_REF] and nonlinear fluid simulations have demonstrated the existence of an anomalous transport due to fluid instabilities related to plasma gradients, drift and collisions effects [START_REF] Smolyakov | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of E x B discharges[END_REF]; large scale structures are often observed experimentally [START_REF] Ellison | Cross-field electron transport induced by a rotating spoke in a cylindrical Hall thruster[END_REF] as well as in numerical simulations [START_REF] Koshkarov | Self-organization, structures, and anomalous transport in turbulent partially magnetized plasmas with crossed electric and magnetic fields[END_REF]. It is widely thought that these instabilities and structures are responsible for an anomalous transport, and thus, are strongly affecting the performance of plasma sources and devices.

The scope of the thesis

As detailed above, the development and advancement of plasma technologies require better modeling of plasma processes and dynamics. Kinetic theory can provide the most complete and accurate description of plasmas [START_REF] Langdon | Plasma physics via computer simulations[END_REF][START_REF] Boeuf | Tutorial: Physics and modeling of hall thrusters[END_REF][START_REF] Lafleur | Anomalous electron transport in hall-effect thrusters: Comparison between quasi-linear kinetic theory and particle-in-cell simulations[END_REF][START_REF] Tavant | The effects of secondary electron emission on plasma sheath characteristics and electron transport in an e x b discharge via kinetic simulations[END_REF][START_REF] Croes | The effect of alternative propellants on the electron drift instability in hall-effect thrusters: Insight from 2d particle-in-cell simulations[END_REF][START_REF] Croes | 2d particle-in-cell simulations of the electron drift instability and associated anomalous electron transport in hall-effect thrusters[END_REF][START_REF] Carlsson | Multi-dimensional kinetic simulations of instabilities and transport in exb devices[END_REF], but they are also very costly numerically; even the simplest 2D cases with many physics phenomena omitted may take weeks and months to calculate. A full size Hall thruster, or other similar devices in a 3D geometry remains out of reach even for modern computers. An alternative for Ey and 0 m -3 for n, and the maximum values are respectively 5 × 10 4 V.m -1 and 5 × 10 17 m -3 . The domain simulation is represented in the left.from Ref. [START_REF] Garrigues | Exb electron drift instability in hall thrusters : Particle-in-cell simulations vs. theory[END_REF].

approach is to use fluid modeling, which is numerically less cumbersome with the benefit to capture and single out the major physical mechanisms in terms of macroscopic quantities. Fluid models importantly contributed to the understanding and the development of the tokamak physics, and still are indispensable for this purpose. Multi-fluid models, in which each species present in the plasma (ions, electrons and neutrals) are modeled as separate fluids interacting with the self-consistent electric field and the magnetic field, and with each other via collisions and ionization, have shown their effectiveness in many studies of instabilities and transport in partially magnetized low-temperature plasmas [START_REF] Smolyakov | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of E x B discharges[END_REF][START_REF] Hagelaar | Role of anomalous electron transport in a stationary plasma thruster simulation[END_REF][START_REF] Chesta | Stability of a magnetized hall plasma discharge[END_REF][START_REF] Koshkarov | Self-organization, structures, and anomalous transport in turbulent partially magnetized plasmas with crossed electric and magnetic fields[END_REF].

Recently, in the LAPLACE laboratory, a quasi-neutral multi-fluid code meant to describe magnetized low-temperature plasmas with E × B configuration has been developed. This code, called MAGNIS ("MAGNetized Ion Source") has already been used for the description and characterization of plasma sources, such as ion source concepts for the neutral beam injection system of ITER and DEMO [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF][START_REF] Simonin | R & d around a photoneutralizer-based nbi system (siphore) in view of a demo tokamak steady state fusion reactor[END_REF], shown fig. 3. The simulations results from MAGNIS revealed for some of these plasma configurations the appearance of instabilities, which, in some cases, even tend to dominate the plasma dynamics. [START_REF] Furno | A novel helicon plasma source for negative ion beams for fusion[END_REF][START_REF] Furno | Helicon wave-generated plasmas for negative ion beams for fusion[END_REF] plasma source for a magnetic field of 200G [START_REF] Fubiani | [END_REF].

This phenomenon was also highlighted in PIC models for these same configurations, and even though the physical model of MAGNIS is mostly limited to large scale dynamics, there are reasons to believe that these instabilities observed in MAGNIS results are physical, since many phenomena occurring under the magnetized low-temperature plasma conditions are of fluid nature.

This represents the first goal of this thesis; through comparisons between numerical simulation results and linear analysis of a simplified magnetized low-temperature configuration, we aim to confirm the physical nature of these observed instabilities, meaning we want to make sure they are a solution of the physical model of MAGNIS and not numerical artifacts. We proceed in two steps; first, we try to understand what mechanisms are behind these instabilities and relate them to unstable modes known from the literature, that are likely to develop under specified conditions of our simplified configuration. For this, an extensive linear analysis via a general linear dispersion relation, attached to our defined configuration, enables to bring out these unstable mode and to study their behaviour. After that, thorough comparisons are made between analytical quantities, deduced by the dispersion relation, and numerical ones given by the code, which allows to achieve this first goal.

The second part of the thesis is dedicated to the study of the evolution these instabilities in the non-linear regime; we propose a qualitative description of the formation of non-linear structures along with a characterization of their properties, and explanations of non-linear mechanisms and their effects. We also proceed to quantify the anomalous transport, a typical non-linear mechanism.

To finish, we will discuss the limitations and problems that we have encountered in our simulations, implying more generally the validity of fluid models. We will present and discuss some further formulations and possible enhancements that allow to extend the validity of fluid models, partially improving them to take into account some kinetic effects and deviations from quasineutrality, which are all important for small scale instabilities.

Plasma type

Plasma density Gas density Temperature Magnetic field

(m -3 ) (m -3 ) T e (

Fusion plasmas

Tokamaks 10 20 10 20 10 3 1 -10

Table 1 -Some reference of typical plasma values; throughout this thesis, we focus on magnetized low-temperature plasmas.

Chapter 1

Plasma basics and modeling

General introduction

Plasma modeling is the numerical description of the state of the plasma obtained by solving a set of physical equations. Contrary to a neutral gas, a plasma consists of different species of neutral and charged particles (neutrals, electrons, ions) which interact with each other not only through collisions but also via electromagnetic fields; furthermore the charged particles respond to external electromagnetic fields applied to the plasma by means of electrodes or magnets. Accordingly, a typical plasma model is based on separate equations describing the different particles species, coupled with equations for the electromagnetic fields. There exist different types of plasma models, using different types of physical equations and approximations. The main two model types are based either on the kinetic or fluid theory. Hybrid models may combine kinetic and fluid descriptions, e.g. kinetic equations are used for one species and fluid for another one.

The choice of a specific method to model a plasma depends on its conditions and, more precisely, on the importance and ordering of different length and time scales involved in the plasma dynamics. In this thesis, we focus on magnetized low-temperature plasma sources. Figure 1.1 introduces the different scales, in frequency and length, that are typically found in these plasma sources. These are:

• the plasma frequency ω p,e and the Debye length λ D , characterizing the electrostatic coupling between electrons and ions;

• the electron cyclotron frequency ω c,e and the electron Larmor radius ρ e , characterizing the electron cyclotron motion due to the magnetic field;

• similarly, the ion cyclotron frequency ω c,i and the ion Larmor radius ρ i ;

• the electron-neutral and ion-neutral collision frequencies ν e,col and ν i,col , respectively, and the corresponding mean free path λ col (in this figure, the range of λ col is considered to be similar for both species);

• a macroscopic variation frequency that we call τ -1 , inverse of a variation or transit time, for the whole considered plasma structure, associated with L n , a macroscopic structure size or gradient length. Let us consider, for example, the Debye length and plasma frequency. These parameters are important for the modeling of the electrostatic coupling between electrons and ions that is characteristic for a plasma. In fact, on length scales larger than the Debye length and time scales larger than the inverse plasma frequency, a plasma is a quasi-neutral medium, meaning that there are as many negative as positive charges. For these larger scales, it can be a good model approximation to assume the electron density directly equal to the ion density, at every point in space and time, which is known as a quasineutral plasma model. As one can see in the figure, for the plasmas of interest in this thesis, all the other scales are larger than the Debye length, so it seems reasonable to use a quasi-neutral model. However, quasi-neutrality becomes invalid at scales below the Debye length, where significant charge separation can occur. This happens for example in a boundary layer near the wall, called plasma sheath, whose size is of the order of a few λ D . In order to model these latter phenomena, it is necessary to describe the coupling between the charged particle dynamics and electric field via the Maxwell equations.

The mean free path is an important parameter for the choice of physical equations used to describe the particle dynamics. In standard fluid dynamics, it is generally considered that the continuum fluid equations are valid only if the mean free path is much smaller than the macroscopic structure size (i.e., the so-called Knudsen number, K n = λ col /L n , must be a small parameter, K n 1); this is necessary for the validity of all kinds of inherent approximations of fluid models (fluid closures). If the mean free path is not small, a kinetic description of the particles can be necessary, tracking the full evolution of the particles in phase space. However, in plasma models, thanks to the other interaction mechanisms and other scales involved, it is sometimes possible to obtain a meaningful description from fluid equations even if the electron or ion mean free paths are larger than the plasma size. The plasmas of interest in this thesis are in this long mean free path regime.

In magnetized plasmas, due to the magnetic Lorenz force, the charged particles are confined (magnetized) and follow complex cyclotron orbits, if the Larmor radius is smaller than both the mean free path and the plasma size. According to Fig. 1.1 this is not the case for the ions in the plasmas of our interest here, which can have quite a large Larmor radius. This implies that these ions are not very much affected by the magnetic field so that it is often a reasonable approximation to neglect to magnetic force in the ion equations. For this reason we also call these plasmas "partially magnetized plasmas" (only the electrons are magnetized).

In this chapter, we will first (section 1.2) review the elementary physical processes in plasmas that are underlying the different scales mentioned above, and then (sections 1.3-4) present the principal general plasma modeling approaches. We introduce basic plasma-physical quantities and equations that will come back throughout this thesis.

Elementary plasma processes 1.2.1 Motion of particles in uniform static fields

In this section we describe the motion of a single charged particle in a constant and uniform electric field E and magnetic field B, implying that these fields are not affected by the plasma dynamics. At the scale of one particle, its behavior can be described by Newton's equations of motion, relating the time variation of its momentum to the applied force :

m dw dt = F = q(E + w × B), (1.1)
dx dt = w, (1.2)
where we introduced the particle mass m, charge q, velocity w and position x, as well as the force F. In the case where we only apply an electric field (B = 0), the force is constant. As a result, solving (1.1) gives the following expression for the particle velocity:

w = q m Et + w(0). (1.3)
Integrating (1.3) leads to an equation for the particle trajectory in space :

x(t) = q m Et 2 + w(0)t + x(0). (1.4)
The particle moves with a constant acceleration in the direction of E if the particle charge is positive (q > 0) or the opposite direction if it is negative (q < 0). If we consider a purely magnetic case with no electric field (E = 0), then the force is proportional to the particle velocity and directed perpendicular to it, which results in a pure gyration motion. In order to show this, it is convenient to decompose the particle velocity into two parts parallel and perpendicular to the magnetic field :

w = w ⊥ + w (1.5)
where w = (w • b)b, with b = B/B a unit vector along the magnetic field. The equation of motion then writes:

m dw ⊥ dt + m dw dt = q(w ⊥ × B), (1.6)
so that:

m dw dt = 0 (1.7) dw ⊥ dt = ±ω c (w ⊥ × b) (1.8)
where ± indicates the sign of the particle charge and

ω c = |q|B m = eB m (1.9)
is the cyclotron frequency. Considering Cartesian position coordinates (x, y, z) with the z axis along the magnetic field direction b, the solution of these equations leads to

x(t) = ρ L sin(±ω c t -φ 0 ) + x 0 (1.10) y(t) = ρ L cos(±ω c t -φ 0 ) + y 0 (1.11)
which correspond to a circular orbit centered at (x 0 ,y 0 ) with a radius

ρ L = w ⊥ ω c , (1.12)
called the Larmor radius, and ω 0 is a phase determined by the initial conditions. Note that the constant particle motion term w(0)t of equation (1.4) is no longer present in these x and y directions: the particle motion is confined. A schematic of the motion is given in Fig. 1.2 along with some explanations.

Figure 1.2 -The gyration motion of a particle in a pure magnetic case, a constant acceleration perpendicular to both the particle velocity and the magnetic field. This does not affect the particle's motion parallel to the magnetic field, but results in circular motion at constant speed in the plane perpendicular to the magnetic field.

Let us now consider the general case where both a magnetic and an electric field are present. Decomposing also the electric field into parts parallel and perpendicular to the magnetic field,

E = E ⊥ + E , (1.13)
we get for the equations of motion:

m dw dt = qE (1.14) m dw ⊥ dt = q(E ⊥ + w ⊥ × B). (1.15)
Equation (1.14) describes a motion of constant acceleration along the magnetic field lines, with a solution for the parallel velocity and position similar to the non-magnetized expressions (1.3) and (1.4). To solve the perpendicular equation (1.15), it is convenient to decompose the perpendicular velocity in two parts :

w ⊥ = ŵ⊥ + w E (1.16)
where

w E = E ⊥ × B B 2 (1.17)
is a constant velocity called the E × B drift velocity, perpendicular to the electric and magnetic field. Substituting this into ((1.15)), we get

d ŵ⊥ dt = ±ω c ( ŵ⊥ × b) (1.18)
which is similar to (1.8) without electric field. Therefore, the solution for the particle trajectory in the perpendicular plane is a superposition of E × B drift and gyration motion (cyclotron motion). In Cartesian coordinates, like before:

x(t) = ρ L sin(±ω c t -φ 0 ) + x 0 + E y B t (1.19) y(t) = ρ L cos(±ω c t -φ 0 ) + y 0 - E x B t. (1.20)
A picture of these trajectories for electrons and ions is shown and explained in Fig. 1.3. Note that the direction of the gyration depends on the sign of the particle charge, but the drift velocity is the same for electrons and ions. Note also that the Larmor radius of ions is much larger than that of electrons because of the much larger ion particle mass.

Collisions

In the low-temperature plasmas of interest in this thesis, the electrons and ions undergo collisions mainly with neutral gas particles, because these plasmas are weakly ionized, meaning that only a small fraction of particles is charged. There are different types of electron-neutral and ion-neutral collisions, which can be classified in two categories: elastic and inelastic. An elastic collision implies that the total kinetic energy of both particles colliding is conserved. In the case where that energy is not conserved during the process, the collision is considered inelastic. An elastic collision of an electron or ion with a neutral gas particle causes elastic scattering: an abrupt random change in the direction of the velocity of the charged particle, leading to a deviation of its trajectory, due to the impact with the neutral (Fig. 1.4). This elastic scattering also involves a transfer of momentum and kinetic energy between the colliding particles, depending on the ratio of their masses. For ions, a common collision process that can be viewed as an elastic collision is resonant charge transfer. Resonant charge transfer is a process that happens between an ion and an atom of the same species, where an electron is transferred from the internal structure of the atom to a fast ion that is passing by. The ion then becomes a rapid atom whereas the atom becomes an ion at rest.

Inelastic collisions occur mainly for electrons and lead to changes in the internal structure of the neutral particle, such as excitation and ionization of that particle. Excitation is a process in which internal configuration and energy of the neutral particle is changed to a different quantum state, while ionization is the process in which an neutral loses an electron and becomes an ion. The colliding electron then usually loses a fixed large amount of kinetic energy (typically around 10-20 eV for ionization). A fundamental quantity to characterize the probability of collisions (on the level of individual particles) is the cross section [START_REF] Lichtenberg | Principles of plasma discharges and materials processing[END_REF] that we denote by σ. This quantity represents the area within which the two particles must meet in order to have a given type of interaction with each other (elastic scattering, ionization, excitation, etc...). If two particles interact upon contact (hard spheres), the cross section is constant and determined by their geometric features, but if the interaction implies a distant action then the related cross section is generally bigger and dependent on the relative velocity of the particles. The cross sections of electron-neutral collisions are usually considered as a function of the electron impact energy, and for inelastic processes they have a threshold, meaning that the cross section is zero below a certain value (because the electron does not have enough energy to enable the inelastic process).

From the cross section, it is possible to define another important quantity representative of an interaction which is the mean free path

λ col = 1 n g σ (1.21)
where n g is the neutral gas density. This is the average distance traveled by a particle between successive collisions. Along with the mean free path, one may also define the collision frequency, corresponding to the inverse average time between collisions and also the collision probability per unit time:

ν micro = n g σw (1.22)
This collision frequency is defined in a microscopic point of view at the particle scale, and must not be mistaken with the macroscopic collision frequency, which is an averaged value of the microscopic collision frequencies over the particle energy distribution functions:

ν α = ν micro .
(1.23)

Quasi-neutrality and plasma sheath

We already mentioned in the introduction that the fundamental characteristic of a plasma its quasineutrality due to the electrostatic coupling between electrons and ions via the Maxwell equations, or more specifically, the Poisson equation:

0 ∇ • E = -0 ∇ 2 φ = e(n i -n e ), (1.24)
where φ is the electrostatic potential, with E = -∇φ, and n i and n e are the ion en electron particle number densities. This equation describes how deviations from neutrality (n e = n i ) generate an electric field, which will then accelerate the charged particles such as to restore quasi-neutrality, i.e. n i ≈ n e or more precisely: .25) This restoring process takes place over a time and length scale given by the plasma frequency and Debye length, respectively, related to electron inertia and electron thermal motions:

|n i -n e | n i + n e . ( 1
ω p = e 2 n e 0 m e , (1.26)
λ D = 0 T e en e , (1.27)
where T e is the electron temperature (in units of Volt). Note that these scales depend on the plasma density n e and become shorter when this increases. When the plasma is in contact with a wall where charged particles are lost by wall recombination processes, an interesting phenomenon happens in the vicinity of said wall: the formation of a sheath. The sheath is a non-neutral layer between the wall and the plasma and maintains the quasi-neutrality condition within the plasma by equilibrating the electron and ion wall losses, repelling electrons from the wall and accelerating ions towards it (this is necessary because electrons have much faster thermal motions than ions due to their much smaller mass). This implies a greater density of ions in the layer, and thus an excess in positive charge, causing a drop of the electric potential (via Poisson's equation).

m i n i v i ∂v i ∂x = -en i ∂φ ∂x , (1.29
)

eT e ∂n e ∂x = en e ∂φ ∂x , (1.30)
Equation (1.28) is the ion continuity equation, (1.29) and (1.30) the momentum equations for respectively ions and electrons (we present these equations later in the following sections). If we consider that the coordinate x = 0 is at the sheath edge, we can then set the boundary conditions as below :

n e (0) = n i (0) = n s , (1.31)
φ(0) = 0, (1.32) v i (0) = v i,s , (1.33)
where n s and v i,s are the density and the ion velocity at the sheath edge, respectively. A simple integration of (1.30) leads to the Boltzmann relation for the electron density : [START_REF] Lafleur | Anomalous electron transport in hall-effect thrusters: Comparison between quasi-linear kinetic theory and particle-in-cell simulations[END_REF] density for electrons at equilibrium. The integration of the equation (1.29) gives the energy conservation of ions :

n e = n s exp -φ(x) T e (1.
1 2 m i v i (x) 2 = 1 2 m i v 2 i,s -eφ(x) (1.35)
and the equation (1.28) illustrates the fact that no ions are created in the sheath area, which implies that the flux is the same everywhere:

n s v i,s = n i (x)v i (x). (1.36)
Solving v i (x) from (1.35) and injecting it in (1.36), one gets the following expression for the ion density:

n i = n s 1 - 2eφ(x) m i v 2 i,s -1/2
.

(1.37)

A necessary condition for the existence of the sheath can be obtained by considering that when moving into the sheath from the sheath edge (in the -x direction in Fig. 1.5), the electron density must drop faster than the ion density (otherwise the space charge becomes negative):

∂(n i -n e ) ∂(-x) s = ∂(n i -n e ) ∂φ s ∂φ ∂(-x) s ≥ 0. (1.38)
Hence, using equations (1.37) and (1.34): .39) This implies that at the sheath edge, where the ions enter into the sheath, their velocity must reach or exceed the ion accoustic speed, also known as the Bohm velocity; this condition is essential for the sheath to appear and is known as the Bohm criterion:

∂(n i -n e ) ∂φ s = 2en s m i v 2 i,s - n s T e ≤ 0. ( 1
|v i,s | ≥ c s = eT e m i . (1.40)
In agreement with this criterion, significant ion acceleration and potential drop already take place within the quasi-neutral plasma region upstream from the sheath, called pre-sheath (see Fig. 1.5).

If we assume that the ions enter the sheath with exactly the Bohm velocity, and that the ion flux is equal to the thermal electron flux at the wall (x = x w ), we get:

n s c s = n e (x w ) 1 4 v th,e = n e (x w ) 1 4 8eT e πm e , (1.41)
where v th,e is the electron thermal speed. With the help (1.40) and (1.34), this then yields the total potential drop that occurs inside the sheath: .42) Typically this is 4 -5 times the electron temperature, depending on the ion mass. Finally, knowing the potential on both sides of the sheath and the ion and electron densities inside it, one can obtain the length of the sheath (along with the profile of Φ(x)) by integrating the Poisson equation. This is not a simple analytical calculation so we will not show it here, but from a dimensional analysis of the Poisson equation it is easy to see that the final result is of the order of a few times λ D , the Debye length.

φ s -φ(x w ) = 1 2 T e ln m i 2πm e . ( 1

Kinetic models

Kinetic models provide a complete statistical description of the plasma in which the particles of a particular species are represented by a distribution function f (x, w, t), particle density in phase space, where the independent variables x, w and t are position, velocity, and time respectively. The fundamental equation used in kinetic theory to describe the evolution of this distribution function is the Boltzmann equation and/or more generally, the Fokker-Planck type equation, in which the collisions effects are included (in the right-hand side):

∂f ∂t + w • ∇f + a • ∂f ∂w = δf δt col (1.43) where a = F m = q m (E + w × B) (1.44)
is the acceleration due to the electromagnetic force F. When the collisions are neglected, this becomes Vlasov equation. A kinetic description of the plasma is achieved by coupling these kinetic equations selfconsistently with Maxwell's equations for the electromagnetic fields, e.g. the Boltzmann-Poisson system or the Vlasov-Maxwell system. Kinetic theory is the most complete comprehensive model for plasma description.

Among the many existing kinetic models, the most used one in low-temperature plasma physics is undoubtedly the Particle-In-Cell (PIC) method [START_REF] Langdon | Plasma physics via computer simulations[END_REF]. This method does not solve the Boltzmann equation explicitly, but is based on tracking of individual electron and ion trajectories in continuous phase space while moments of their distribution functions such as densities and currents are simultaneously computed and coupled with the Maxwell equations. More precisely, if we take the example of an electrostatic case with a static magnetic field (which applies well for our magnetized low-temperature plasmas), the method involves integrating the Newton equations of motion (1.1) and (1.2) for a group of simulation particles, taking into account the self-consistent electric field, solved simultaneously from the Poisson equation (1.24). Collisions of the simulation particles with the neutral gas are included via Monte-Carlo collision sampling (PIC-MCC).

A PIC code follows a computational cycle in which, first started by given appropriate initial conditions for the particle positions x and velocities w, determines for each time step, the velocity and position at the particle frame, while the fields are solved on a discrete spatial grid. From there, the link between the particle quantities and the fields is made by calculating the charge and current densities on the grid (for that, the particle charges are weighted to the grid points surrounding each particle position). From these densities, we obtain the electric and magnetic fields, still on the grid. As the fields are known on the grid points and particles are scattered around those, it is required to interpolate the fields from the grid to the particle frame to apply the forces on the particle using again a weighting method (fig. 1.6). This cycle is illustrated in fig. 1.6. Usually, the numerical method used to integrate the motion equations is a simple second order Leap-Frog method, a good compromise between accuracy, stability and efficiency. In a pure electrostatic case, the scheme is stable when ω p ∆t ≤ 0.2 with ω p = (e 2 n e / 0 m e ) 1/2 the plasma frequency, however, in the electrostatic case with a static magnetic field, the time step must be small enough to resolve the electron gyration, namely the electron cyclotron frequency, such that ω c,e ∆t < 0.15. Despite its utmost physical accuracy, the main drawback in PIC is the heavy numerical description of the multiscale system; indeed, one has to resolve all scales present in the plasma, the smallest one being the Debye length λ D = ( 0 T e /qn 0 ) 1/2 in space, and the inverse of the plasma frequency ω p in time. As one can easily see in the Debye length expression, the higher is the density, the smaller λ D will be, which implies that the grid cells have to be smaller, and the same goes for the time step in relation to the plasma frequency, all also coupled with the constraint of a high enough number of particle per cell.

It is important to note that although the Boltzmann equation is not directly solved, the electrostatic PIC method is strictly equivalent to a solving a system of Boltzmann and Poisson equations:

∂f e ∂t + w • ∇f e - e m e (E + w × B) • ∂f e ∂w = δf e δt col (1.45) ∂f i ∂t + w • ∇f i + e m i (E + w × B) • ∂f i ∂w = δf i δt col (1.46) 0 ∇ • E = (f i -f e )d 3 w. (1.47)
Some other kinetic models, such as Lattice Boltzmann for instance, base their formalism on solving directly the Boltzmann equation. Kinetic models are the most fundamental way to describe a plasma, but as seen with the PIC method, they require cumbersome numerical computations. For for many realistic configurations and parameters, kinetic models are too expensive and remain out of reach even for modern computers. Alternatively, fluid models can provide a reduced description and offer the possibility to cut small scales to avoid the computational constraints met on kinetic models, even though the effects of these scales are then no longer captured. Kinetic models remain indispensable for some cases where the distribution function becomes very different from Maxwellian. In the following section, we will describe and detail thoroughly the general fluid theory for a partially magnetized plasma; particularly, we show how the fluid equations can be obtained from the Boltzmann equation.

Fluid models 1.4.1 Fluid quantities and equations

The fluid model meets the need of overcoming difficulties in the kinetic modeling by describing the plasma based on its macroscopic properties [START_REF] Bittencourt | Fundamentals of plasma physics[END_REF]: density, mean velocity, and mean energy. These quantities are defined as velocity moments of the distribution function mentioned above, meaning we integrate the function distribution, multiplied by a certain power of the microscopic velocity, over velocity space. We first define non-centered moments such as:

n(x, t) = f (x, w, t)d 3 w, (1.48) v(x, t) = w = 1 n(x, t) wf (x, w, t)d 3 w, (1.49) Σ(x, t) = mn(x, t) ww = m wwf (x, w, t)d 3 w, (1.50) Θ(x, t) = mn(x, t) www = m wwwf (x, w, t)d 3 w, (1.51)
respectively the particle number density n, mean velocity v, energy density tensor Σ, and total energy flux tensor Θ. More precisely, (1.48) represents the number of particles per unit volume for a given species, (1.49) is the corresponding average velocity of the particles, while (1.50) is the total kinetic energy of the particles per unit volume and (1.51) is the flux of kinetic energy crossing an unit area per unit time. Note that we defined the macroscopic average (for an arbitrary quantity X) as

X = 1 n(x, t) Xf (x, w, t)d 3 w. (1.52)
Making use of the mean velocity in (1.49), it is then possible to define the centered moments. For that purpose, we set a "centered velocity" which represents the deviation with respect to the mean velocity:

u = w -v. (1.53)
Thanks to this velocity, we can express centered moments of second and third order, representing respectively the pressure tensor and the heat flux tensor:

P(x, t) = mn(x, t) uu = m uuf (x, w, t)d 3 u (1.54) Q(x, t) = mn(x, t) uuu = m uuuf (x, w, t)d 3 u. (1.55)
These quantities are related to the non-centered moments of second and third order given above, substituting (1.53) in their expressions, expanding and using the definitions above:

P = Σ -mnvv (1.56) Q = Θ -mnvvv -(v, P) (1.57)
where (v, P)

= v i P jk + v j P ki + v k P ij .
All these macroscopic variables must be seen as average values of physical quantities involving the collective behavior of a large number of particles.

The pressure tensor

Usually, the kinetic pressure tensor is decomposed in two different parts, both with different contributions, as follows [START_REF] Braginskii | Transport processes in a plasma[END_REF]:

P = pI + π = pδ ij + π ij (1. 58
)
where I is the identity tensor and

p = 1 3 T rP (1.59)
is the trace of the total P tensor i.e. the sum of its diagonal elements, known as the normal scalar pressure, while π represents the viscous stress tensor, part of P ij resulting from a symmetry deviation and defined as

π ij = nm u i u j - 1 3 u 2 δ ij .
(1.60)

Both P and π are symmetrical tensors (M ij = M ji ), and in an anisotropic case, the diagonal elements of the viscous stress tensor are not zero (but their sum is).

The heat flux vector

In most fluid models, it is common and more convenient to reduce the 10-component heat flux tensor Q, defined above in (1.55), to a 3-component vector we name q, the heat flux vector. It is easy to obtain this vector by tensor contraction:

Q = mn uuu → q = 1 2 mn u 2 u (1.61)
In the same way, it is also possible to define a total energy flux vector thanks to the total energy flux tensor, so that:

Θ = mn www → e = 1 2 mn w 2 w . (1.62)
The relation established in (1.57) is conserved and valid for the reduced forms of the heat and energy flux tensors. We show that with some calculation, this relation becomes:

e = q + Pv + 3 2 p + 1 2 nmv 2 v. (1.63)

Transport equations

In the previous paragraphs, we defined the main macroscopic quantities controlling our system. They are related to each other by are the basic fluid equations, namely the continuity, momentum, and energy equations, which can be obtained by taking velocity moments of the Boltzmann equation. We proceed to demonstrate this procedure below; we first write the integral of the Boltzmann equation (1.43) in the velocity space, multiplied a particular order of the velocity w which we represent here as function L(w).

We thus write:

L(w) ∂f ∂t + w • ∇f + a • ∂f ∂w d 3 w = L(w) δf δt col d 3 w, (1.64)
Using definition (1.52) of the macroscopic average, one gets for the first two terms in the left hand side of (1.64):

L(w) ∂f ∂t + w • ∇f d 3 w = ∂ nL(w) ∂t + ∇ • ( nL(w)w ). (1.65)
If we rewrite the third term of the equation (1.64), it becomes:

L(w)a • ∂f ∂w d 3 w = L(w) ∂(af ) ∂w -f ∂a ∂w d 3 w (1.66)
where ∂a ∂w = 0, (1.67) due to the fact that the divergence of the electromagnetic acceleration (1.44) is non-existent. Thus, thanks to an integration by parts, we have:

L(w) ∂(af ) ∂w d 3 w = - f a ∂L(w) ∂w d 3 w = -n a ∂L(w) ∂w . (1.68)
We then define the last term on the right side of the equation such that:

L(w) δf δt col d 3 w = δ nL(w) δt col , (1.69)
and finally, equation (1.64) becomes:

∂ nL(w) ∂t + ∇ • ( nL(w)w ) = n a ∂L(w) ∂w + δ nL(w) δt col (1.70)
Equation (1.70) represents a general transport equation where the two terms on the right side can be considered as a source term, taking into account the external forces and the collisions. In order to obtain the macroscopic conservation equations, we use equation (1.70) where we replace the function L(w) by a particular order of w for each conservation law we aim to recover, meaning the first three equations in our case, being the continuity (mass conservation, here L(w) = m), momentum equation (L(w) = mw), and the energy equation (L(w) = mww). Hence, the macroscopic equations of transport can be written as follows:

∂(nm) ∂t + ∇ • (nm w ) = δ(nm) δt col (1.71) ∂(nm w ) ∂t + ∇ • (nm ww ) = n F + δ(nm w ) δt col (1.72) ∂(nm ww ) ∂t + ∇ • (nm www ) = n wF + Fw + δ(nm ww ) δt col . (1.73)
This set of equations is then rewritten thanks to the non-centered moment quantities defined previously (1.48, 1.49, 1.50, 1.51), the most common forms used in literature:

∂n ∂t + ∇ • (nv) = δn δt col (1.74) ∂(mnv) ∂t + ∇ • Σ = n F + δ(mnv) δt col (1.75) ∂Σ ∂t + ∇ • Θ = n wF + Fw + δΣ δt col , (1.76)
where the expression of the averaged Lorentz force F is given by

F = q(E + w × B) = q(E + v × B). (1.77)
In plasma physics, it is convenient to express the transport equations with centered moments rather than non-centered moments as we did above. We proceed to rewrite our equations with the centered moment quantities by substituting the relations (1.56) and (1.63):

∂n ∂t + ∇ • (nv) = S (1.78) ∂(mnv) ∂t + ∇ • (mnvv) + ∇ • P = n F + R (1.79) ∂( 3 2 p + 1 2 mnv 2 ) ∂t + ∇ • 3 2 p + 1 2 mnv 2 v + ∇ • (P • v) + ∇ • q = n F • w + C, (1.80)
where we reduced the energy equation to scalar form, thanks to the tensor contraction of Σ and Θ (equation (1.63)), and we introduced the following collisional sour terms:

S = δn δt col (1.81) R = δ(mnv) δt col (1.82) C = δ( 1 2 mn w 2 ) δt col .
(1.83)

The particle source term S represents the rate per unit volume at which particles of a considered species are produced or lost as a result of collisions, and R and C are the rates of change of their momentum density and the energy density, respectively, resulting from collisions. These terms gather contributions from different collision processes such as ionization, excitation, and so on. By combining the above transport equations with each other, one can rewrite them in different forms that are also often used. For example, taking the scalar product of the momentum equation (1.79) with the mean velocity v, and combining with the continuity equation (1.78), we get the following equation for the transport of directed energy 1 2 nmv 2 :

∂( 1 2 nmv 2 ) ∂t + ∇ • 1 2 nmv 2 v + v • (∇ • P) = n F • v + R • v - 1 2 mSv 2 , (1.84)
where we used

v • ∇ • (mnvv) = ∇ • 1 2 nmv 2 v + m∇ • (nv) v. (1.85)
We can then subtract the directed energy equation (1.84) from the total energy equation (1.80) in order to find an equation for the internal energy only. Noticing that

∇ • (P • v) -v • (∇ • P) = (P • ∇) • v (1.86) and n F • w -n F • v = 0, (1.87)
we then obtain the internal energy equation 3 2

∂p ∂t + 3 2 ∇ • (pv) + (P • ∇) • v + ∇ • q = C -R • v + 1 2 mSv 2 .
(1.88)

The momentum equation (1.79) can be expressed in the non-conservative form thanks to the continuity equation (1.78) by developing the second term as follows:

∇ • (nmvv) = mn(v • ∇)v + m∇ • (nv)v = mn(v • ∇)v -mv ∂n ∂t + mvS. (1.89)
Once the first term of the momentum equation (1.79) is developed and (1.89) is injected, we get the following non-conservative form for the momentum equation:

mn ∂v ∂t + nm(v • ∇)v + ∇ • P = n F + R -mSv. (1.90)
Similar non-conservative forms can be obtained for the other fluid equations, leading to the non-conservative system of respectively continuity, momentum and energy equations:

Dn Dt = S -n∇ • v (1.91) Dv Dt = F m - ∇ • P mn + R mn - S n v.
(1.92)

3 2 Dp Dt = -R • v -∇ • q + C + 1 2 mSv 2 - 3 2 p∇ • v -(P • ∇) • v, (1.93)
where we introduced the material derivative

D Dt = ∂ ∂t + v • ∇ (1.94)

Closures in fluid models

The fluid equations contain several quantities that need to be determined before the system is closed and can be solved such as the pressure tensor P, the heat flux vector q, and the collision terms S, R and C. Typically these parameters will be linked to the main fluid variables solved from the system (density, mean velocity, temperature) by so-called closure relations. To determine the closure relations, assumptions must be made about the behavior of the distribution function in velocity space. A standard assumption is that this velocity distrubtion function is very close to an isotropic, Maxwellian distribution function:

f (u) = f (0) (u) + f (1) (u) (1.95)
where f (1) f (0) and

f (0) (u) = n m 2πeT 3/2 exp - mu 2 2eT (1.96)
with T the kinetic temperature (which we express in units of Volt). This corresponds to the equilbrium distribution function that particles acquire under the influence of random interactions between each other (e.g. collisions within the same species). In classical fluid dynamics, this assumption is well justified as long as the mean free path is short (λ c L n ) but in plasmas this is more complicated due to the other interactions involved (as we discussed in the introduction section). Significant deviations f (1) from the Maxwellian distribution function can arise due to non-local kinetic effects or due to collisions with particles of other species that have a different distribution function, for example collisions with a cold background gas. For specific configurations, the behavior of these deviations can sometimes be predicted from a kinetic analysis, leading to special closures for that configuration. In such analysis, the distribution function is typically developed in an orthogonal polynomials basis:

f (x, w, t) = f (0) (x, w, t) β(x, t)P (x, w, t) (1.97)
where f (0) (x, w, t) is the chosen distribution function (possibly non-Maxwellian) around which is made the development, β(x, t) the coefficients resulting from the development and P (x, w, t) the polynomials. For example, this kind of approach is commonly used for electrons in weakly ionized gas discharges in order to obtain electron transport coefficients from a local Boltzmann analysis [52]. However, most often basic Maxwellian-type closures are used, involving the kinetic temperature T , often based mainly on phenomenological grounds. Below we briefly describe some of the most common closures used in plasma models.

Pressure tensor

A common approximation for weakly ionized plasmas, dominated by collisions with the neutral gas, is to neglect the viscous stress part π of the pressure tensor (1.58):

π = 0 ⇔ P = pI.
(1.98)

Injecting the Maxwellian distribution function (1.96) into the expression for the scalar pressure yields

p = 1 3 nm u 2 = neT. (1.99)
which is none other than the ideal gas law. The dynamics of this scalar pressure and temperature is then calculated from the energy equation or simply deduced from an assumed thermodynamic law (isothermal, adiabetic). In magnetized plasmas, significant anisotropy of the velocity distribution function may arise but the pressure tensor can still be diagonal when expressed in a coordinate system aligned with the magnetic field, with the particularity that the components perpendicular to the magnetic field of the tensor are different than the parallel component, so that the pressure tensor becomes [START_REF] Ramos | General expression of the gyroviscous force[END_REF]

P = p ⊥ I + (p -p ⊥ )bb ⇔ P =    p ⊥ 0 0 0 p ⊥ 0 0 0 p    (1.100)
where b = B/B a unit vector in the direction of the magnetic field and

p ⊥ = 1 2 mn u 2 x = 1 2 mn u 2 y = neT ⊥ (1.101) p = mn u 2 z = neT . (1.102)
Taking into account this kind of anisotropic pressure requires the use of separate energy equations for the perpendicular and parallel directions, not covered by our derivations of the previous section. (In this thesis we will neglect this anisotropy and assume that there is only one temperature.) Furthermore, in magnetized plasmas it can be important to take into account the viscous stress tensor π for the description of certain effects due to the finite Larmor radius of the particle trajectories (FLR effects), via a so-called gyro-viscosity closure [START_REF] Smolyakov | Gyroviscous forces in a collisionless plasma with temperature gradients[END_REF][START_REF] Ramos | General expression of the gyroviscous force[END_REF]. This results in a gyro-viscous force that is of the same order as the convective inertia term of the momentum equation and partially cancels out with this term, which is called gyro-viscous cancellation [START_REF] Smolyakov | Gyroviscous forces in a collisionless plasma with temperature gradients[END_REF][START_REF] Belova | Nonlinear gyroviscous force in a collisionless plasma[END_REF].

Heat flux

If an energy equation is included in the model, a closure is required for the heat flux. Usually, the closure for this quantity is either to put it to zero or to determine it with Fourier's law, which we write in its most simple form as follows:

q = -κ∇T (1.103)
where κ is a thermal conductivity coefficient. Equation (1.103) is a very basic form which does not take into account non-stationary effects and can apply for magnetized plasma if the magnetic field effects are included, such that

q = -κ ∇ T -κ ⊥ ∇ ⊥ T -κ × ∇ ⊥ T × b (1.104)
The thermal conductivity coefficients depend on collisions and involve averages of the microscopic collision frequencies over the distribution function. The dependence of these collision frequencies on particle energy can give rise to additional terms in the above heat flux expressions [START_REF] Golant | Fundamentals of Plasma Physics[END_REF].

Collision terms

The collision terms of the fluid equations generally involve averages of microscopic collision probabilities (1.22) over the distribution functions of the colliding particle species. For electron collisions, these averages depend mainly on the electron distribution function, and when assuming Maxwellian electrons, on the electron temperature. In addition, the collision terms are generally proportional to the particle densities of the colliding particle species. In weakly ionized plasmas dominated by collisions with a nonflowing neutral gas at rest, the momentum transfer term R (friction force) is often assumed to be of the following form:

R = -ν m nv, (1.105)
where ν m is a macroscopic momentum transfer frequency proportional to the gas density. In case the microscopic collision frequency depends strongly on the particle energy, there may be an additional contribution to R that is proportional to the temperature gradient, the so-called thermal force [START_REF] Golant | Fundamentals of Plasma Physics[END_REF].

Simplified fluid models

Drift-diffusion models

Drift-diffusion models are fluid models that are very used in the low-temperature plasma modeling, and assume that collisions between charged particles and neutrals are dominant in the transport phenomena. These models can either take into account a magnetic field or not. The equations of these models come as follows. We first define the velocities from the momentum equation (1.79, 1.90) in which the inertial terms are neglected due to their weak contribution compared to the collisions (the main assumption of the model): The same logic applies to a magnetized drift-diffusion model, except that now some mathematical manipulation is needed to solve the mean velocity from the momentum equation

∇ • P = n F + R (1.
v = ±µ(E + v × B) -D ∇n n , (1.110)
which, with some mathematical manipulation, becomes:

v = ±µE -D ∇ n n + 1 1 + h 2 ± µE ⊥ -D ∇ ⊥ n n + h 1 + h 2 µE ⊥ -±D ∇ ⊥ n n × b (1.111)
where h = µB and the parallel and perpendicular vector components are defined by E = (E • b)b and E ⊥ = E -E , with b = B/B a unit vector along the magnetic field, and similar for ∇. This can be

written as v = ±µ • E -D • ∇n n (1.112)
µ and D, being second order tensors, respectively mobility and diffusion, defined such that

µ =    µ ⊥ µ × 0 µ × µ ⊥ 0 0 0 µ    and D =    D ⊥ D × 0 D × D ⊥ 0 0 0 D    (1.113) where µ = µ µ ⊥ = µ 1 + h 2 µ × = ±µh 1 + h 2 (1.114)
are the mobility components for the directions parallel and perpendicular to the magnetic field lines, and for drift in the cross-field direction; the diffusion component are defined in a similar way. These components can be very different from each other depending on the so-called Hall parameter:

h = µB = ω c ν m . (1.115)
When this parameter is zero, the tensors defined above become isotropic and one recovers the unmagnetized drift-diffusion model defined before. The larger the value of this parameter, the smaller are the mobility and diffusion coefficients perpendicular to the magnetic field, i.e. the stronger the magnetic confinement of the particles. In the case of magnetized low-temperature plasma sources, the usual values for the Hall parameter for electrons can go up to h 10 2 or more, meaning that the perpendicular electron transport is very strongly reduced. Note that the parallel transport is not affected by the magnetic field.

Low-frequency approximation models

Although the drift-diffusion model described in the previous section can be adapted for magnetized plasmas by means of anisotropic mobility and diffusion tensors (as we showed), it is not appropriate for the description of magnetized plasma instabilities, because these often depend on inertia effects. Therefore, magnetized plasma models usually include the low-frequency, large scale effects of particle inertia via an approximation method that we will outline below. It is assumed that the particles are strongly magnetized, such that their Larmor radius and cyclotron period are much smaller than the macroscopic transport scales and the collisional scales (h 1). The equilibrium solution of the perpendicular momentum equation is then dominated by the electric and magnetic forces and the pressure force:

v (0) ⊥ = E × b B - ±∇(nT ) × b nB ≡ v E + v p (1.116)
where the two terms are the E × B drift velocity and the diamagnetic drift velocity, and ± is the again the sign of the particle charge, which only affects the diamagnetic term. Note that v E is the same as the drift velocity (1.17) of individual particle trajectories but that the diamagnetic drift is a macroscopic effect not seen on the trajectories. Note also that (1.116) corresponds to the limit of the drift-diffusion expression (1.111) for h → ∞.

The zeroth order solution (1.116) is injected into the inertia and collision terms terms of the momentum equation, in order to obtain a first order correction:

∂ ∂t + v (0) ⊥ • ∇ + ν m v (0) ⊥ + ∇ • π = q m v (1)
⊥ × B.

(1.117)

Here usually also the viscous stress force is taken into account because this cancels out with part of the inertia terms, the gyroviscous cancellation [START_REF] Smolyakov | Gyroviscous forces in a collisionless plasma with temperature gradients[END_REF]:

(v p • ∇)v (0) ⊥ + ∇ • π = 0. (1.118)
This yields:

v (1) ⊥ = m qB 2 ∂ ∂t + v E • ∇ E ⊥ - ±∇ ⊥ (nT ) n + mν m qB 2 E ⊥ - ±∇ ⊥ (nT ) n ≡ v pol + v c (1.119)
The first of these term is called polarisation drift or inertial drift, the second is collision drift and corresponds to the second (⊥) term of the magnetized drift-diffusion equation (1.111). The total mean velocity can then be written as a sum of drift velocities:

v ⊥ = v (0) ⊥ + v (1)
⊥ = v E + v p + v pol + v c . (1.120) Chapter 2
The fluid code MAGNIS

Introduction

In this chapter we will present the numerical fluid code that we used in this thesis. This code, called MAGNIS (MAGNetized Ion Source), is a quasi-neutral multi-species fluid code specially designed to describe magnetized low-temperature plasmas in the conditions of ion sources for neutral beam injection, Hall thrusters, Penning traps, magnetrons and similar magnetized low-temperature plasma devices. [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF] It was developed at the LAPLACE laboratory in Toulouse, initially to study magnetized plasma transport in the ion source for the ITER neutral beam injection system, and is currently being used in other studies and research projects. MAGNIS is implemented as a serial Fortran code and runs on a desktop computer. One of the main purposes of this thesis is to asses the capabilities of this code to describe magnetized plasma instabilities and to understand its limitations. The physical model of MAGNIS comprises different sets of fluid equations modeling each species (electrons, ions) as a separate fluid. These equations are kept in rather brute form as presented in Section 1.4.1 without further approximations such as the drift-diffusion approximation or low-frequency approximation (see Section 1.4.3), so that they can apply for different conditions and orderings of scales (collisional or non-collisional, magnetized or non-magnetized). These equations are all coupled together and solved to describe the plasma dynamics in the 2D plane perpendicular to the magnetic field lines, while the plasma behavior and losses in the third (parallel) direction are averaged along the field lines.

In the following sections, we describe the physical and numerical principles of the MAGNIS and then present a few examples of simulation results. For the presentation in this chapter, several points are important to note:

• We will present the equations solved by MAGNIS for the case of a simple electropostive plasma with only one positive ion species and minimal plasma chemistry, for the sake of simplicity and also because this is what was used during the thesis. However, MAGNIS can in principle be extended to take into account multi-component plasmas with more ion species and more complex plasma chemistry, and even neutral gas dynamics.

• Throughout this chapter, unless otherwise indicated, all spatial derivatives and differential operators and vectors refer to the 2D simulation plane perpendicular to the magnetic field lines, i.e. they are perpendicular to the magnetic field, but for simplicity we omit the '⊥' subscripts: ∇ → ∇ ⊥ , v → v ⊥ , etcetera. A separate section will explain how some effects of the parallel direction are included in the model.

Physical model 2.2.1 Main equations

The first equations of the model are continuity equations (1.78) for electrons and ions:

∂n e ∂t + ∇ • (n e v e ) = S e -L e , (2.1
)

∂n i ∂t + ∇ • (n i v i ) = S i -L i . (2.2)
Here S is the particle source term which represents ionization, while L represents losses of particles along the magnetic field lines (this term is explained later in Section 2.2.3). In our case with only one ion species, the electron and ion source terms are equal:

S e = S i = S = n e ν iz = n e n g k iz (T e ) (2.3)
where n g is the neutral gas density, ν iz is the ionization frequency (or electron/ion creation frequency), and k iz (T e ) an ionization rate coefficient that is a given function of the electron temperature. The neutral gas density is fixed. The mean velocities are described by momentum equations in the non-conservative form (1.90),

mn ∂v ∂t + nm(v • ∇)v + ∇ • P = qn(E + v × B) + R -mSv. (2.4)
It is assumed that the velocity distribution function is Maxwellian so that the pressure tensor is diagonal and isotropic (1.98, 1.99) so that

P = enT I ⇒ ∇ • P = e∇(nT ), (2.5)
and that the collision term is of the simple form (1.105), i.e. R = -ν m nv with ν m the momentum transfer frequency due to collisions with the gas at rest. The electron and ion momentum equations used in MAGNIS are then written as:

α 1 ∂v e ∂t + α 2 v e • ∇v e + ν e v e = - e m e (E + v e × B) - e∇(n e T e ) m e n e , (2.6
)

β 1 ∂v i ∂t + β 2 v i • ∇v i + ν i v i = e m i (E + v i × β 3 B) - e∇(n i T i ) m i n i , (2.7)
where ν is an effective frequency defined as:

ν = ν m + S n = ν m + ν iz = n g (k m + k iz ). (2.8)
This frequency includes a contribution from ionization which takes into account that new particles are created from the neutral gas and added to the fluid with zero momentum; for ions this can be dominant over the contribution ν m from collisions. In equations (2.6) and (2.7), the parameters α 1,2 and β 1,2 are model switch parameters allowing to retain (1) or remove (0) the inertial terms; setting all these parameters to zero corresponds to applying the drift-diffusion approximation (see Section 1.4.3). Similarly, the parameter β 3 enables to switch on/off the magnetic force on the ions. The above equations are coupled together by the assumption that the plasma is quasi-neutral, which means in our one-ion case that the electron and ion densities are considered to be equal:

n e = n i = n.
(2.9)

This implies that the self-consistent electric field is not calculated from the Poisson equation but is deduced from the current conservation equation, obtained by subtracting the ion and electron continuity equations:

∇ • j = e∇ • (n i v i -n e v e ) = e(L e -L i ), (2.10)
where j = e(n i v i -n e v e ) is the plasma current density. The right hand side of this equation comes from currents entering or leaving the plasma along the magnetic field lines (remember that the left hand side only represents the perpendicular plane). Notice also there is no displacement current term, corresponding to the the sum of the non-stationary terms of the continuity equations, as the model is quasi-neutral. The electric field E = -∇φ and potential φ are calculated such that the electron and ion mean velocities from (2.6, 2.7) satisfy the current conservation equation (2.10), while the plasma density is calculated from the ion continuity equation (2.5). The magnetic field is constant with a given spatial profile or uniform in space.

The electron temperature is solved from the electron internal energy equation (1.88), expressed as 3 2

∂en e T e ∂t + 3 2 ∇ • (en e T e v e ) + ∇ • q e = 2n e ν e -S e + L e 1 2 m e v 2 e -en e T e ∇ • v e + H + C -L T , (2.11)
where H is an external heating power density, C is the power density due to collisions with the gas and L T represents a loss of energy along the field lines (see Section 2.2.3). The external heating H represents electron heating mechanisms such as RF heating that are not self-consistently described in MAGNIS but just assumed as input, with a given total power P tot and spatial profile P such that:

H = P tot n e P V n e P dV , (2.12)
where the integral is over the total volume. The collisional term C is found from

C = -n e n g κ(T e ) = - 2m e m g n e ν m,e 3 2 e(T e -T g ) -n e n g j ε j k j (T e ), (2.13)
where the first term on the left-hand side represents elastic energy losses and the last term, inelastic losses from excitation and ionization collisions with threshold energy ε j and rate coefficient k j ; all these are grouped into one energy loss coefficient κ, again a function of T e . Finally, the heat flux vector q e is determined from [START_REF] Golant | Fundamentals of Plasma Physics[END_REF] ∂q e ∂t + ν e q e + e m e q e × B = -5 2

en e T e m e ∇T e . (2.14)
The stationary solution of this equation is of the same form as the magnetized version of the Fourier law closure (1.104). The electron energy equation can play an important role because the electron temperature controlls not only the electron pressure but also very strongly the ionization source term. The ion temperature is fixed. (MAGNIS also contains energy and heat flux equations for ions but these are not used because they have little effect on the results.)

In the above equations, the collision terms and frequencies are determined from rate coefficients that are input data for the model. The rate coefficients related to electrons are functions of the electron temperature, read in look-up tables that are pre-calculated with the BOLSIG+ software [52,[START_REF] Hagelaar | Solving the boltzmann equation to obtain electron transport coefficients and rate coefficients for fluid models[END_REF] from cross section data [START_REF]Siglo database[END_REF], for a Maxwellian electron energy distribution function (EEDF) such that:

k(T e ) = σ(ε) 2eε m e = 2e m e ∞ -∞ f ε (ε, T e )σ(ε) √ εdε. (2.15)
Here σ is the cross section of the particular process, as a function of the electron impact energy ε, and the EEDF is given by

f ε (ε, T ) = 2 √ π √ ε T 3/2 exp - ε T (2.

16)

. This formula applies for the ionization rate coefficient k iz and the rate coefficients in the energy source term C. However, it is worthwhile to note that the expression for the momentum transfer frequency is a bit different (e.g. see [START_REF] Golant | Fundamentals of Plasma Physics[END_REF]):

k m (T ) = 1 ε εσ m (ε) 2eε m = 2e m 2 3T ∞ -∞ f ε (ε, T )σ(ε)ε 3/2 dε, (2.17)
with σ m the momentum transfer cross section. Figure 2.1 shows these rate coefficients for hydrogen and argon.

Boundary conditions

The boundary conditions of the model, describing the losses of particles and energy to walls or electrodes in contact with the plasma, are obtained from theoretical expressions for the particle and energy fluxes normal to the wall, depending on the nature of the wall (the choice can be made between a dielectric wall or an electrode). The boundary conditions for the continuity equations are thus expressed as [START_REF] Hagelaar | Model of an inductively coupled negative ion source: I. general model description[END_REF]:

(nv) • n = nW, (2.18)
in which n is the unit vector normal to the wall and pointing towards it, and W is an effective wall loss speed decuded from theory. Since the model is quasi-neutral, the effect of the plasma sheath is taken into account in the effective speed W on the basis of the classical sheath theory, described in Section 1.2.3.

For electrons and ions this yields: [START_REF] Choueiri | Anomalous resistivity and heating in current-driven plasma thrusters[END_REF])

W e = eT e 2πm e 1/2 exp - (φ -φ w ) T e (2.
W i ≥ eT e m i 1/2 ⇔ W i = max eT e m i 1/2 , v i • n (2.20)
where φ is the potential in the plasma at the sheath edge (a variable of the model) and φ w is the potential at the wall surface, and (2.20) is the Bohm criterion. Subtracting the wall fluxes of ions and electrons, we find the boundary condition for the current conservation equation (2.10):

j • n = en i W i -en e W e (φ, φ w ), (2.21)
from which the boundary value of φ is determined self-consistently, depending on the chosen type of wall surface. For the dielectric wall (an electrical insulator), the current density must be zero whereas the wall potential φ w is free to adjust itself to this condition:

j • n = 0 ⇒ W e = W i ⇒ φ w = φ - 1 2 T e ln T e 2πm e W 2 i . (2.22)
In the case of an electrode (an electrical conductor) it is the opposite: the wall potential φ w is fixed as a model input parameter (e.g. φ w = 0 for grounded wall) while the current density adjust itself. Note that if the ions reach exactly the Bohm velocity (W i = c s ), then the dielectric wall boundary condition (2.22) reproduces the classical sheath potential drop (1.42), but in the magnetized plasmas of interest it often happens that the ions become supersonic (W i = c s ). The boundary conditions for the electron energy equation is obtained from the next energy flux carried by the electrons that are lost to the wall across the sheath. According to classical sheath theory:

5 2 T e en e v e • n + q • n = en e W e (2T e + φ -φ w ) ⇒ (2.23) q • n = en e eT e 2πm e 1 2 exp - (φ -φ w ) T e φ -φ w - 1 2 T e .
(2.24)

Transport losses along the magnetic field lines: 2 1 2 D approach

The physical model defined above, consisting of the set of fluid equations (continuity, momentum and energy equations) and their boundary conditions, applies for the 2D Cartesian plane perpendicular to the magnetic field lines. However, the effects of particle and energy losses by diffusive transport in the third, parallel, direction are taken into account in the 2D model via the parallel loss terms L that are included in the continuity and energy equations. These losses occur whenever the magnetic field lines are intercepted by walls, as is usually the case in the plasma sources of interest. They are estimated by assuming that the parallel direction is governed by classical equilibrium relations for low-pressure plasma (Maxwell-Boltzmann equilibrium for electrons, Bohm criterion for ions) which are then averaged over the length of the magnetic field lines. In this view, all the variables of the model actually represent averages over the field line length, for example in case of the continuity equation:

∂n ∂t + ∇ ⊥ • (nv ⊥ ) ≈ ∂n ∂t + ∇ ⊥ • (nv ⊥ ) = S -L (2.25)
where the overlines indicate averinging and

L = ∇ • (nv ) = 1 L z Lz/2 -Lz/2 ∇ • (nv )dz = 2 L z n w W = 2a L z nW . (2.26)
Here L z is the plasma chamber size along the magnetic field lines, W and n w are the wall loss speed and the plasma edge density at the end of the lines. The parameter a = n w /n takes into account the plasma density drop between the averaged value and the edge value. This parameter is usually known as the h factor (but we do not call it that to avoid confusion with the Hall parameter) and for normal low pressure plasmas it has a constant value of 0.5 -0.7, depending on the plasma geometry [START_REF] Lichtenberg | Principles of plasma discharges and materials processing[END_REF]. The wall loss speeds W at the end of the field lines are given by the same expressions as those given in the previous section for W in the boundary conditions, with the additional simplification that the ions cannot be supersonic in the parallel direction, so W i = c s . Hence, omitting the overlines as in the other sections, the expressions for the parallel loss terms become:

L e = 2 L z n e eT e 2πm e 1/2 exp - (φ -φ w ) T e (2.

27)

L i = 2a L z n i eT i m i 1/2 (2.28)
where φ w is the wall surface potential at the end of the field lines, depending on whether this wall is a dielectric or an electrode. For a dielectric wall one gets:

L e = L i ⇒ φ w = φ - 1 2 T e ln m i 2πm e . ( 2

.29)

Note that the (2.27) does not contain the factor a (because the density drop is already included in the Boltzmann factor). The expression for the parallel electron energy loss term (at the end of equation (2.11)) is:

L T = L e 2T e + φ -φ w . (2.30)
MAGNIS also proposes the choice to set all these parallel loss terms to zero, L e = L i = L T = 0, in which case it is a pure 2D model; this is what we considered for most of the results in this thesis.

Numerical aspects and procedure

The system of model equations described in the previous sections is solved numerically by an original FORTRAN 90 code which runs in serial mode on a normal desktop computer. All the input and output is handled via text files. In this section we will briefly describe the numerical methods used to solve the equations. Since most of the development and implementation of these methods was done outside the scope of this thesis, we will only sketch here the overall solution strategy and highlight the most essential aspects. For more details we refer to [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF][START_REF] Hagelaar | Model of an inductively coupled negative ion source: I. general model description[END_REF]. The code solves the set of equations on 2D rectangular Cartesian computational domain, discretized by a uniform grid, with the magnetic field B perpendicular to this domain. The equations are solved one by one in a computational time stepping cycle starting from some simple initial spatial profiles for all the variables. The methods used to solve the equations are essentially the finite volume method [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] for the conservation equations (ion continuity, current conservation and electron energy equations) and an original finite difference scheme [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF] for the momentum and heat flux equations. These methods are based on a representation of the scalar variables (densities, temperatures, source terms) and different components of vector viarables (velocities, fluxes) on a different grids displaced by half a cell, which means that the scalars are defined at the grid points (nodes) while vector components are midway between the grid points, as illustrated in figure 2.2. The finite volume method makes use of control volumes centered around the grid points in order to ensure that the fluxes are strictly conserved, also in the numerical approximation. In the following subsections we describe the computational cycle and the methods for each equation.

Time integration cycle

In the main computational cycle, all the model variables are advanced in time over small time steps ∆t by the following sequence:

1. The code first calculates all frequencies ν, source terms S and C, and wall loss speeds W .

2. Then, the new velocities v are predicted from the momentum equations (2.6) and (2.7).

3. The electric field E and potential φ are solved from the current conservation equation (2.10), taking into account the response of the particle velocities to changes in E (such as to correct the velocities to satisfy current continuity).

4. The particle velocities are corrected with the changed E 5. The code then solves the ion density n i from the ion continuity equation (2.2), using the new corrected ion velocity.

6. The electron density n e is updated from the quasineutrality condition.

7. After that, the new electron heat flux q is predicted thanks to the equation (2.14) (in analogy with the particle velocities).

8. The new electron temperature T e is solved from the electron energy equation (2.11), taking into account the response of the heat flux to changes in T e . 9. The heat flux is corrected for the changed T e .

This cycle is repeated continuously until the end of the simulation.

Momentum equation

Because of the magnetic cross term, the momentum equation shows to be the hardest to handle numerically. We write it as follows:

α 1 ∂v ∂t + νv + ωb × v = q m F ≡ q m E - e q ∇(nT ) n -α 2 m q v • ∇v (2.31)
where ω = qB/m, b = B/B, F is an effective field variable (not the force here) and α 1,2 are the model switch parameters included in equation (2.6). This equation is then discretized in time as:

δ(v k+1 -v k ) + νv k+1 + ωb × v k+1 = q m F, (2.32)
with superscripts indicating discrite time (t k+1 = t k + ∆t) and

δ = max α 1 ∆t , |ω| , (2.33)
where the lower limit |ω| = ω c is included to under-relax the velocity on the cyclotron timescale for larger time steps, so that the scheme remains stable for ∆t >|ω| -1 , and this even without the inertial terms (meaning α 1 = 0), although it is of course necessary to respect ∆t <|ω| -1 if one wishes to resolve plasma dynamics related to cyclotron motion. We can then obtain the new velocity:

v k+1 = A(F + δv k ) + B(b × F + δb × v k ) (2.34)
with coefficients A and B are given by:

A = q m δ + ν (δ + ν) 2 + ω 2 B = - q m ω (δ + ν) 2 + ω 2 (2.35)
However, after many time steps, the grid interpolation needed for the calculation of the cross term b × v causes an accumulation of numerical errors leading to convergence problems. To limit these errors, both the velocity v and crossed velocity v × = b×v are tracked at each point of the grid (at the control volume edges, see figure 2.2) at the same time:

v k+1 = A(F + δv k ) + B(b × F + δ((1 -ξ)v k × + βb × v k )) (2.36) v k+1 × = -B(F + δv k ) + A(b × F + δ((1 -ξ)v k × + βb × v k )). (2.37)
Here, ξ is a parameter giving control over the cumulative velocity interpolation errors. When ξ = 0, no velocity interpolation is necessary but v × and v tend to get desynchronized and the scheme doesn't work. Some velocity interpolation (ξ > 0) is in fact necessary to keep v × consistent with v, but too much interpolation (ξ = 1) causes the interpolation errors to explode at small time steps. One can show that a good choice for ξ is [START_REF] Futtersack | Modélisation fluide du transport magnétisé dans les plasmas froids[END_REF]:

ξ = ν 2 + ω 2 (δ + ν) 2 + ω 2 1/2 . ( 2.38)
With this value for ξ the scheme converges nicely for all conditions and orderings of ω, ν and (∆t) -1 . Alternatively, for the case of non-or weakly magnetized ions, MAGNIS can also solve the ion momentum equation via a more standard discretization method [START_REF] Hagelaar | Model of an inductively coupled negative ion source: I. general model description[END_REF], better adapted to handle the (v • ∇)v term in the presence of shocks.

Current conservation equation

Here we explain how the electric field is calculated from the current conservation equation (2.10). A first prediction of the electron and ion velocities, v k+1 pre , is done from the momentum equations using the old electric field E k . We then wish to update the field in such a way that these velocities become consistent with current conservation. From equation (2.34) it is easy to see that if the electric field changes by ∆E = -∇ φ, the velocities change by

∆v = A∆E + Bb × ∆E = -A∇ φ + Bb × ∇ φ (2.39)
where φ = φ k+1 -φ k is the potential change. Therefore, the predicted velocities can be corrected by writing

v k+1 cor = v k+1 pre -A∇ φ -Bb × ∇ φ, (2.40)
and substituting this into the current conservation equation (2.10):

∇ • n i v k+1 i,pre -n e v k+1 e,pre -(n i A i -n e A e )∇ φ -(n i B i -n e B e )b × ∇ φ = L e -L i . (2.41)
This equation is solved for φ, after which the potential and the field are updated from φ k+1 = φ k + φ, and

E k+1 = E k -∇ φ

Continuity equation for ions

Once the corrected velocities are obtained, we get new the ion density through the continuity equation (2.5). The discretization is based on an explicit second-order MUSCL (Monotonic Upwind Scheme for Conservation Laws) scheme with an appropriate flux limiter:

n k+1 i = n k i + ∆t S i -∇ • (n k i v k+1 i ) (2.42)
The particularity of this scheme is the use of several nodes upstream from the flux to give a good estimation of the density at the edges of the control volume (see figure 2.2) needed to calculate the transport term (last term), while preventing numerical oscillations to happen. It can provide highly accurate numerical solutions, even in cases where the said solutions exhibit shocks, discontinuities, or large gradients. However, the time step must be limited by a CFL (Courant Friedrichs Levy) condition:

(usually not severe) ∆t < 1 2

|v i,x | ∆x + |v i,y | ∆y -1
.

(2.43)

Electron heat flux and energy

The heat flux vector is calculated from equation (2.14) in analogy with the calculation of the velocities explained in Section 2.3.2, where this time F = -∇T e and the A and B coefficients are A T = (5/2)en e T e A e and B T = (5/2)en e T e B e . First a prediction of the new heat flux, q k+1 pre , is made from the old electron temperature T k e , then this is corrected by inserting

q k+1 cor = q k+1 pre -A T ∇(T k+1 e -T k e ) -B T b × ∇(T k+1 e -T k e ) (2.44)
in the electron energy equation (2.11) and solving for T k+1 e .

Examples of MAGNIS results

In this section we illustrate the results obtained from the MAGNIS by some examples.

Magnetic filter

Our first example is the simulation of the plasma Hall-effect in ion sources using a magnetic filter, such as the the negative ion source for the ITER neutral beam injection system [START_REF] Speth | Development of powerful rf plasma sources for present and future nbi systems[END_REF][START_REF] Stephen | Positive and negative ion sources for magnetic fusion[END_REF] or the PEGASES thruster at Ecole Polytechnique [START_REF] Aanesland | Electron energy distribution function and plasma parameters across magnetic filters[END_REF]. Studying these sources was one of the initial motivations for the development of MAGNIS. A schematic picture of these ion sources is shown in Figure 2.3. They are charaterized by recangular plasma chamber that is divided into two parts by means of a magnetic field barrier called "magnetic filter": (1) a "driver" region (left side in Figure 2.3, orange dashes) where plasma is created by an inductive radio-frequency (RF) gas discharge; and (2) an extraction region (right side in Figure 2.3) where negative ions are produced and extracted from the plasma by a system of biased grids (blue dashes on the right). The electrons in the driver region being very energetic due to the RF field, the role of the magnetic filter is to reduce the electron temperature in the extraction region, which is necessary to production of negative ions by plasma chemistry and the survival of the negative ions created, and also to limit the electron current that is co-extracted when these ions are exacted by the grids. The plasma transport from the driver region across the magnetic filter plays a central role in the operation of the source. The main dynamics of this transport occurs in the 2D plane perpendicular to the magnetic field lines (Figure 2.3 left view). Figure 2.4 shows MAGNIS simulation results of this 2D plane, for a source geometry corresponding approximately to that of PEGASES. For the sake of simplicity, the discharge gas in these simulations is argon and there are no negative ions (as in early experiments on this source [START_REF] Aanesland | Electron energy distribution function and plasma parameters across magnetic filters[END_REF]). The 2D domain size is 0.12 m by 0.12 m and losses along the magnetic field lines are included for a length L z of also 0.12 m in the third direction. The RF discharge in the driver region is modeled in via the external heating power term H in the electron energy equation (2.11,2.12) with a total injected power of 100 W and a narrow spatial profile P near the left boundary. The magnetic field has a Gaussian profile in the horizontal direction with a maximum at a position indicated by the vertical dashed lines in the figure. All the walls are considered to be grounded, meaning that the wall potentials φ w and φ w are fixed to 0 V in the electron boundary conditions and parallel loss terms, except at the right domain boundary, where a wall potential φ w = 15 V is applied, representing the bias voltage of the first plasma grid in contact with the plasma. The results in figure 2.4 illustrate several interesting points:

• The electrons are heated on the left side of the source, where they reach a temperature of around 6 eV, necessary to sustain the plasma by ionization. The electron temperature on the right side is much lower (about 2 eV), simply because of the absence of local heating there and the poor heat transport across the magnetic filter.

• As a result of this, the electron pressure drops strongly from the right to the left side of the source. This pressure gradient, in combination with the applied potential of 15 V at the right boundary, drives the electrons in the horizontal direction from the left to the right. Naturally, this generates combination of diamagnetic drift and E × B drift in the vertical direction, from the bottom to the top.

• However, since this magnetic drift is blocked by the bottom and top walls, the plasma develops an assymmetric gradient in the vertical direction, which tilts the whole density and potential profiles such that the drift can cross the magnetic filter in an oblique way, as one can see in the right-most figure. This effect is similar to the Hall effect in solids and is sometimes called the plasma Hall effect [START_REF] Hagelaar | Plasma transport across magnetic field lines in low-temperature plasma sources[END_REF][START_REF] St | Physics of a magnetic barrier in low-temperature bounded plasmas: insight from particle-in-cell simulations[END_REF][START_REF] Boeuf | Physics of a magnetic filter of negative ion sources. ii. exb drift through the filter in a real geometry[END_REF]. The Hall effect also involves the electron heat flux and causes an inclination of the isothermal lines in the magnetic filter. (These plasma asymmetries can be a problem for the uniformity of the negative ion current over the extraction grids under real negative ion source operation.)

MAGNIS simulations of this kind of magnetic filter configuration are generally in good agreement with results from fully kinetic particle-in-cell (PIC) simulations [START_REF] St | Physics of a magnetic barrier in low-temperature bounded plasmas: insight from particle-in-cell simulations[END_REF][START_REF] Boeuf | Physics of a magnetic filter of negative ion sources. ii. exb drift through the filter in a real geometry[END_REF] (see Section 1.3). To illustrate this, Figure 2.5 shows some PIC simulation results from [START_REF] St | Physics of a magnetic barrier in low-temperature bounded plasmas: insight from particle-in-cell simulations[END_REF] that strongly resemble our fluid results in Figure 2.4 (although the plasma conditions are slightly different and the discharge gas is hydrogen instead of argon). In addition, various MAGNIS results were successfully compared with measurements on a magnetic filter experiment at the LAPLACE laboratory [START_REF] Baude | Étude du transport de particules chargées dans un plasma froid basse pression magnétisé[END_REF] (profiles of electron temperature, plasma density and collected electron and ion currents).

Periodic magnetic filter with closed drift

In the second example we consider the same geometry as in the previous section, except that we remove the top and bottom walls and replace them by periodic boundary conditions. This is shown in Figure 2.6. Naturally, the Hall effect disappears because there are no more walls to induce it. The magnetic drift is now closed along the vertical direction, meaning that it goes out at the top and comes back in at the bottom, so it does not cross the magnetic field barrier anymore. As a result, the total electron current collected at the right domain wall is much smaller than in the previous case and the electron temperature downstream from the filter drops much more (¡ 1 eV, not shown here). [START_REF] Hagelaar | Model of an inductively coupled negative ion source: I. general model description[END_REF] Of course this periodic simulation case does not correspond to a real magnetic filter ion source, but it is similar to the design of many other magnetized plasma devices such as Hall thrusters and magnetrons, in the sense that these devices have a cylindrical geometry where the magnetic drift is closed along the azimuthal direction, in order to avoid the Hall effect and optimize the magnetic confinement of the electrons.

However, both in these cylindrical plasma devices and in the periodic MAGNIS simulation case of figure 2.6, the closed drift is not stable but develops small scale fluctuations of the electric field and plasma density in the vertical direction, which propagate in this direction with a wave velocity much smaller than the electron drift velocity (from top to bottom). It turns out that such plasma instabilities occur very generally in MAGNIS simulations of closed drift configurations. On first sight, these instabilities sometimes resemble those found in PIC simulations. For example, the results in Figure 2.6 look very similar to those in PIC simulations of an idealized 1D magnetic filter (see Figure 2.7 reproduced from [START_REF] Boeuf | Physics of a magnetic filter of negative ion sources. i. collisional transport across the filter in an ideal, 1d filter[END_REF]) and even to the electron-cyclotron drift instabilities causing anomalous transport in PIC simulations of Hall thrusters [START_REF] Adam | Physics, simulation, and diagnostics of hall effect thrusters[END_REF] (see also Figure 6 in the introduction of this thesis). However, it is questionable whether these instabilities are really the same. Many magnetized plasma instabilities are known or presumed to involve kinetic phenomena that are not captured by standard fluid equations and fluid model closures. It is thus questionable whether fluid models like MAGNIS are capable of capturing properly the behavior and effects of instabilities in low-temperature plasma sources, and if so, which types of instabilities and under which conditions. It is the goal of this thesis to answer these questions, or to make progress towards answering them.

Chapter 3

Waves and instabilities in partially magnetized plasma

Preliminaries and basic eigen-modes

Wave phenomenon is a typical macroscopic behaviour happening in plasmas. A wave can be defined as a disturbance in the system state that propagates through space and time. For large amplitudes, plasmas can support non-linear waves, such as shock waves, and other large amplitude phenomena, which require nonlinear theory and are unpredictable in linear models. For small amplitudes of perturbations, a linear analysis is appropriate. In such cases, the underlying equations are linearized with respect to the perturbations, Ψ (x,t) = Ψ 0 (x) + Ψ (x,t) , Ψ Ψ 0 . Here Ψ 0 (x) and Ψ are respectively equilibrium and perturbed quantities. In our case, Ψ(x, t) can be any variable of the system defining the plasma (density n, particle velocity v, electric field E, ...), representative of the wave motion [START_REF] Krall | Principle of plasma physics[END_REF][START_REF] Pecseli | Waves and oscillations in plasmas[END_REF].

In a linear situation, the amplitude of the oscillation is small enough and we can look for the general solution in form of the Fourier integral

Ψ(x, t) = Ψ k,ω exp(ik • x -iωt)d 3 kdω. (3.1)
But the most common is to only use a single harmonic, such that :

Ψ(x, t) = Ψ 1 exp(ik • x -iωt) + c.c. (3.2)
where Ψ 1 represents the amplitude and the "c.c." term indicates the complex conjugate. However, in the linear analysis it suffices to work with a complex quantity, the first term in (3.2).

The wave can be characterized by its phase velocity, which defines the wave-front propagation, and can be written like :

v ϕ = ω k 2 k (3.3)
This velocity is expressed with two quantities, ω, the wave's angular frequency and k, the wave vector, in the direction of the wave propagation, defined as a spatial gradient of the wave phase. The wave number k, norm of the wave vector, can be related to λ, the wavelength, spatial period of the wave such that :

λ = 2π k . (3.4)
Considering now a wave packet, namely the superposition of several harmonic waves, which frequency and wave vector respectively are between ω + dω and k + dk, then we define the group velocity such that:

v g = ∂ω ∂k . (3.5)
This is the velocity with which the overall envelope shape of the wave propagates through space. One can also show that the group velocity describes the velocity of the energy propagation.

It is important to note that ω and k are functionally linked through a dispersion relation, ω = ω (k) . The dispersion relation is obtained as a solubility condition for the underlying linearized system. In other words, dispersion relation defines the eigenmodes of the system. Normally, we will consider the boundary value problem, that is for a given real value of the wavevector k, that is determined by boundary conditions, eg. by the size of the simulation domain, we will look for the corresponding value of ω from the dispersion relation ω = ω (k) . Thus the dispersion relation identifies the eigen-modes that are to be found in a specific system. A dispersion relation can give either purely real or complex value of the frequency,

ω(k) = ω r (k) + iγ(k) (3.6)
Here, ω r (k) is the real part, frequency of the wave, and γ(k) is the growth or damping rate, nothing else than the imaginary part. If γ(k) = 0 or negative, it means that the medium is stable or damped respectively. More precisely, an oscillation is considered heavily damped if -γ |ω r | while the term "wave" is more suitable for a weakly damped oscillation satisfying the condition -|ω r | ≥ γ. In a case where γ > 0, the medium becomes unstable and the growing solution resulting from the dispersion relation is called an instability, a case that we detail further in the following sections. Fluid theory usually lead to the dispersion relations in form of the various orders (depending on the model) algebraic polynomials of ω.

In system without dissipation, all coefficients of the polynomials are real, so the complex roots of ω may appear in complex conjugate pairs, this is so called reactive instabilities. When the dissipation is included, some coefficients can be complex, and then the instabilities are called dissipative. Both examples will occur in partially magnetized plasmas as is discussed below. In kinetic theory, the dispersion relations can be transcendental with different types of complex roots.

The values for the phase and the group velocities are defined above in (3.3) and (3.5) from the real part of ω. Plasmas are often strongly dispersive, so the real part of frequency is a complex function of the wavevector. For the linear dispersion laws ω kv , typical for sound waves or electromagnetic waves (light) in vacuum or weakly dispersive transparent material, the phase velocity does not dependent on k and both velocities v ϕ and v g are equal and constant. In dispersive case, the phase velocity is a function of k and as a consequence, leads to v ϕ = v g . Plasmas often show strong dispersion cases, e.g. such as

ω 2 = ω 2 pe + k 2 c 2 for electromagnetic, or electrostatic ω 2 = ω 2 pe + 3k 2 v 2
T e /2 waves, or ion sound waves

ω 2 = k 2 c 2 s / 1 + k 2 λ 2
D in unmagnetized plasmas. For our work, the ion sound waves are most relevant and will be discussed further below.

Many different types of waves in plasmas can exist and can be classified depending on the nature of the perturbations (whether electrostatic or electromagnetic), the direction of the propagation with respect to the magnetic field, k B, or perpendicular k ⊥ B, or general case, range of frequencies with respect to some characteristic frequencies, eg. cyclotron or collisions frequencies, range of the wavelengths, and in many other different ways.

In this work we consider only the electrostatic case, the propagation in the plane perpendicular to the magnetic field, the plasma being partially magnetized and quasi-neutral. The types of waves that are relevant to our situation are the ion-sound, lower hybrid, and drift (gradient) waves.

Ion sound waves. These waves are oscillations of electron and ions in the plasma, analogous to classical acoustic waves traveling at the sound speed in a neutral gas, but achieved by a mechanism in which the electric field provides the coupling between the ion inertia and the electron pressure.

Electrons are considered at dynamic equilibrium, massless and quickly establishing the equilibrium in the time dependent electric field. Thus they can be described by the Boltzmann relation and expanded such that :

n e = n = n 0 exp eφ T e = n 0 1 + eφ T e + ... (3.7)
So that the perturbation for the density is

ñ = e φ T e n 0 (3.8)
The ions are described by the continuity and momentum equations

∂n ∂t + ∇ • (nv i ) = 0 (3.9) m i ∂v i ∂t = eE (3.10) 57
Then, we proceed to linearize the system 3.9 and 3.10 leading to :

∂ ñ ∂t + ∇ • (nṽ i ) = 0 (3.11) ∂ṽ i ∂t = - e m i ∇ φ (3.12)
Assuming (3.2) for the perturbed quantities ψ and inserting it in (3.11) and (3.12), one gets the following system :

-

eφ 1 T e iω + ik • v i,1 = 0 (3.13) v i,1 = e m i ω φ 1 k (3.14)
And finally, inserting (3.14) in (3.13), the dispersion relation for the ion sound wave :

ω 2 = k 2 c 2 s (3.15)
with c s = T e /m i , the ion sound speed, velocity with which this mode is propagating.

It is important to note that the above derivation is for the plasma without magnetic field. We consider the plasma with unmagnetized ions, so the ions are not affected by the magnetic field and all equations for ions remain valid. The electron dynamic is more complex with magnetic field. It will be considered below. For now, we only note that for relatively short wavelength k 2 ρ 2 e 1, the electrons become unmagnetized and their density is again described by the equation (3.8) as it will be shown below, so one gets the so called short wavelength ion sound waves propagating perpendicular to the magnetic field.

The lower hybrid waves. These are oscillations of ions and electrons across the magnetic field.

Both species are considered cold. The physics of this mode is based on the balance of the unmagnetized ion inertia against the magnetized electron transverse inertia. We write the full equations of the system which consist of continuity and momentum equations for ions and electrons.

∂n ∂t + ∇ • (nv i ) = 0 (3.16) ∂n ∂t + ∇ • (nv e ) = 0 (3.17) m i ∂v i ∂t = eE (3.18)
As it is said above, the electron transverse inertia plays a role in this mode; thus the momentum equation for cold electrons is m e ∂v e ∂t = -e(E + v e × B).

(3.19) 58
We linearize the set of equations giving the following perturbed system :

∂ ñ ∂t + ∇ • (n 0 ṽi ) = 0 (3.20) ∂ ñ ∂t + ∇ • (n 0 ṽe ) = 0 (3.21) m i ∂ṽ i ∂t = e Ẽ (3.22) m e ∂ṽ e ∂t = -e( Ẽ + ṽe × B) (3.23)
and, considering perturbations in the form of (3.2), one gets :

n 1 n 0 ω -k • v i,1 = 0 (3.24) n 1 n 0 ω -k • v e,1 = 0 (3.25) v i,1 = eφ 1 k m i ω (3.26) v e,1 = em e B 2 φ 1 ωk - ieφ 1 B 2 k × B. (3.27)
Inserting (3.26) and (3.27) respectively in (3.24) and (3.25), and combining them, one finally obtains the following dispersion relation for the lower-hybrid waves :

ω 2 = ω c,i ω c,e (3.28)
where ω c,i = eB/m i and ω c,e = eB/m e , are the ion and electron cyclotron frequencies. Contrary to the ion sound mode, the mode in this simple model does not depend on k . Also, it is interesting to notice the presence of ω c,i in this expression even if we considered the ions unmagnetized; in fact, the presence of this term must be considered just as a way to write this frequency's expression and should not be literally interpreted as a magnetization of ions.

The gradient-drift waves. These wave result from a density gradient appearing in the electron dynamics described by continuity and momentum equations, in which we consider cold electrons. Thus, the representative equations of this mode are for the electrons just like (3.17) and (3.19), except with no inertia, such that:

∂n e ∂t + ∇(n e v e ) = 0 (3.29) v e = E × B B 2 , (3.30)
and for ions, the same system as in the ion sound wave :

∂n i ∂t + ∇(n i v i ) = 0 (3.31) m i ∂v i ∂t = eE. (3.32)
Linearizing the above system and considering a Boussinesq approximation (k x >> g) leading to the removal of the gradient term in the ions continuity, one may get the following equations for perturbations:

∂ ñe ∂t + ṽe ∇n 0 + n 0 ∇ṽ e = 0 (3.33) ṽe = Ẽ × B B 2 = -∇ φ × B B 2 (3.34) ∂ ñ ∂t + n 0 ∇ • (ṽ i ) = 0 (3.35) ∂ṽ i ∂t = - e m i ∇ φ (3.36)
and using the Fourier form of fluctuations in the linearized system leads to :

-iω n e,1 n 0 + v e,1 • ∇n 0 n 0 = 0 (3.37) v e,1 = -iφ 1 k × B B 2 (3.38) - n i,1 n 0 iω -k • v i,1 = 0 (3.39) v i,1 = - e m i ω iφ 1 k. ( 3
= g • k × B eωB 2 T e eφ 1 T e (3.41) n i,1 n 0 = k 2 T e m i ω 2 eφ 1 T e = k 2 c 2 s ω 2 eφ 1 T e (3.42)
and equalizing (3.41) and (3.42) while considering quasineutrality, one finally find the expression of the gradient-drift frequency :

ω = k 2 c 2 s ω * , (3.43)
such that :

ω * = eT e k × B B 2 • g (3.44)
is the drift frequency. This mode shows a linear dependency on the wavector for the frequency, thus the phase velocity is constant. Due to the inverse dependence on the drift frequency it is also called the anti-drift mode. The expression of the frequency displays an electron temperature (the ion sound mode velocity due to the electron pressure), however it is cancelled by the temperature in the drift frequency. For g =g x, and k =k y y one has for (3.43)

ω = ω ci k y g , (3.45)
These eigen-modes, typical waves developing in the sources we chose to study, can be destabilized due to some mechanisms also present in the plasma leading to appearance of instabilities. In the following sections, we explain what are these mechanisms and how they impact on the eigen-modes defined above and their dispersion relation.

Mechanisms for instabilities

As we identified the main eigen modes that are likely to develop in our case, we focus on the mechanisms that may destabilize them, and give rise to linear instabilities. An instability can be defined as a disturbance or variation in the quasi-equilibrium state of the system that reduces the free energy. Therefore, one can consider a classification [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF] according to the type of free energy source triggering and driving the instability. This is presented below with some additional remarks :

Rayleigh-Taylor type instabilities. In this case, there is a non-uniformity of the plasma (density gradient) present and an applied external force such as gravity balances the pressure gradient force. Classical Rayleigh-Taylor instability was considered in fluid mechanics, occurring for the heavy fluid maintained by the gravity on top of the light one. In laser plasmas, the role of gravity is played by the acceleration, when the plasma heated by laser beam expands outward. In other situation, in configuration similar to ours, the applied potential acts like gravity so some instabilities in partially ionized ionosphere also called of the Rayleigh-Taylor type.

Streaming instabilities. This instability occurs when a driven current or a beam of energetic particles travels through the plasma, causing drifts between the different species or within one same species. Concerning the latter, these instabilities can be considered as caused by gradients in the phase space (velocity space), rather than gradients in real space, e.g. via inverse Landau damping, when ∂f ∂w > 0, which implies a deformation of the distribution function from Maxwellian form (c.f. kinetic instabilities).

Drift instabilities (Universal instabilities). A finite plasma which is magnetically confined has gradients in plasma density, temperature and other parameters. The plasma gradients in magnetic field lead to plasma drifts of various species. Thus relative drifts of ions and electrons is responsible for drift instabilities. Since gradients were always present, it was thought that plasma will be always unstable, hence the name of universal instabilities. This was how they were called first [START_REF] Chen | Introduction to plasma physics and controlled fusion[END_REF][START_REF] Mikhaïlovskii | Instabilities in a confined plasma[END_REF][START_REF] Kadomtsev | Plasma turbulence[END_REF][START_REF] Kadomtsev | Instability of the positive column in a magnetic field and the 'anomalous' diffusion effect[END_REF][START_REF] Galeev | Universal instability of an inhomogeneous plasma in a magnetic field[END_REF]. Actually, it was later realized that magnetic shear stabilizes it, so they are not considered universal anymore [START_REF] Mikhaïlovskii | Instabilities in a confined plasma[END_REF]. Nowadays, it is clearly shown that these drift instabilities due to gradients in density, temperature and magnetic field, considered the most important for fully magnetized fusion plasmas, in particular, the ion and electron temperature gradient instabilities [START_REF] Mikhaïlovskii | Instabilities in a confined plasma[END_REF][START_REF] Kadomtsev | Plasma turbulence[END_REF]. Kinetic instabilities. They appear when the processes in phase space are involved, e.g. deviation from thermodynamic equilibrium due to a distortion of the Maxwellian distribution function, and kinetic resonances are involved, such as Landau damping, electron and/or ion-cyclotron resonances. Many instabilities may exist in fluid and/or kinetic regimes, e.g. ion/electron temperature gradient instability when the Landau or toroidal resonance are involved [START_REF] Mikhaïlovskii | Instabilities in a confined plasma[END_REF].

One can notice the tricky character of instabilities nomenclature. Many definitions overlap and the classification is not so rigid and depends on the particular aspect, emphasis, or application, and different names can be used to emphasize various aspects of the same mode. For example, universal (drift) instabilities and streaming instabilities, can be also classified as kinetic part as they involve Landau damping and/or a distortion of the Maxwellian distribution function. Thus the same instability can be classified as kinetic, streaming and universal (drift), depending on what aspect one may want to emphasize.

However, even if not rigid in absolute terms, this classification is useful to provide an insight on the main mechanisms that may trigger unstable modes when we refer to our own configuration. Indeed, if we take a closer look on the partially magnetized plasma sources, we can identify the different types of free energy sources leading to instabilities. For instance, the applied external force (here the electric field E) along with gradients (density gradients) may lead to the Rayleigh-Taylor type instability. The inertia of ions (in both collisional and collisionless regimes), and electron flow cause a relative drift between electrons and ions results in the streaming instability. For our applications, the drifts between different species are most important as ions are not sensitive to the magnetic field.

It is thus possible to summarize the destabilizing mechanisms : the electron E × B drift flow, the ion flow across the magnetic field, and the space gradients present in the system such as the density, electron temperature and magnetic field gradients. We will consider these mechanisms to the eigen-modes identified in partially magnetized plasmas in the previous section, and then identify the instabilities that will appear under these conditions. The following section will be dealing with the study of such unstable modes.

Linear instabilities in partially magnetized plasmas

This section deals with linear instabilities that are likely to be found in our plasma source of interest. We have discussed previously the main eigen modes possible in this configuration and the mechanisms that may destabilize them. Here we present more details characterising the linear instabilities one may observe in our plasma source and present relevant dispersion relations.

Farley-Buneman instability

The Farley-Buneman instability, studied by D. T. Farley [START_REF] Farley | Two-stream plasma instability as a source of irregularities in the ionosphere[END_REF] and O. Buneman [START_REF] Buneman | Excitation of field aligned sound waves by electron streams[END_REF] in the ionosphere context, is a low-frequency plasma instability driven by a strong electric field E 0 (considered quasistationary), perpendicular to the geomagnetic field B 0 . The instability occurs in a particular region of the ionosphere, the weakly ionized E-region, where the electrons are magnetized, while ions are unmagnetized due to frequent collisions with neutral particles, ν i ω ci. The instability does not require any density gradients. The electrons are moving relative to ions with the drift velocity v e,0 ≈ E 0 × B 0 /B 2 0 , so the mechanism could be seen as a variant of two-stream instability. The instability occurs when when the electric field E 0 exceeds some threshold value (usually around 10-20 × 10 -3 V m -1 ), corresponding to the condition that the electron velocity with respect to the ions exceeds the ion acoustic speed. Even if the collisions dominate ion dynamics, their inertia is important for the instability [START_REF] Smirnov | Modeling of the farley-buneman instability in the e-region ionosphere: a new hybrid approach[END_REF]. In general, the Farley-Bunemann instability can be seen as the destabilization of the ion sound wave by the E × B flow in presence of dissipation.

The dispersion relation is obtained from the quasi-neutrality, considering unmagnetized ions (ν i,m >> ω c,i ) and massless electrons (m e << m i ). The continuity equation for both species is used

∂n α ∂t + ∇ • (n α v α ) = 0, (3.46)
α = (i, e), where the velocity found from the momentum equations :

-e(E + v e × B) - ∇(nT e ) n -m e nν e,m v e = 0, (3.47)
m i ∂v i ∂t = eE - ∇(nT i ) n -m i nν m,i v i . (3.48)
The temperatures and the magnetic field are assumed constant. We linearize the above equations considering fluctuations in density, velocities and electric field in the form ψ ≈ exp(ik Similarly, we can obtain the ion velocity from (3.48):

• x -iω).
v i,1 = k iν m,i + ω e m i φ 1 + v 2 th,i n 1 n 0 (3.51)
where v th,e and v th,i are the electron and ion thermal velocities respectively. Inserting (3.50) and (3.51) into the linearized continuity equations for both species leads to :

[(ω -k • v e,0 ) + ik 2 ν m,e ω 2 c,e v 2 th,e ] n 1 n 0 = i e m e ν m,e ω 2
c,e

k 2 φ 1 (3.52) n 1 n 0 (ω(ν m,i -iω) + ik 2 v 2 th,i ) = -i e m i k 2 φ 1 (3.53)
which, by combination, leads finally to the dispersion relation for the Farley Buneman instability :

ω -k • v e,0 ω(ν m,i -iω) -ik 2 c 2 s = ψ ν m,i (3.54)
with c 2 s = (T e + T i ) /m i the ion acoustic speed, and :

ψ = ν m,e ν m,i ω c,e ω c,i (3.55)
a parameter representing an effect of collisions, typically ψ 1. In the limit γ < ω r , it is easy to obtain the real and the imaginary parts of ω in the form:

ω r = k • v e,0 1 + ψ , (3.56)
γ = ψ (1 + ψ)ν m,i (ω 2 r -k 2 c 2 s ) (3.57)
where (3.56) and (3.57) are respectively the frequency the growth rate of this unstable mode. One may easily see from (3.57), the necessary condition for the instability is :

ω 2 r -k 2 c 2 s > 0, (3.58)
hence the following stability criteria :

(k • v e,0 ) 2 k 2 c 2 s > (1 + ψ) 2 . (3.59)
Expression (3.59) underlines the condition for the drift velocity to exceed the ion sound speed for this instability to occur [START_REF] Hassan | Ion thermal and dispersion effects in farleybuneman instabilities[END_REF]. This instability is important in the ionosphere [START_REF] Oppenheim | Saturation of the farley-buneman instability via nonlinear electron exb drifts[END_REF][START_REF] Hamza | A fully self-consistent fluid theory of anomalous transport in farley-buneman turbulence[END_REF], however it may also occurs in our plasma of interest, such as in magnetized low-temperature plasma devices. As we will see later, the condition (3.59) also occurs in somewhat different but related conditions of the plasmas of interest.

Simon-Hoh, or gradient-drift instability.

The gradient drift instability, often called the Simon-Hoh instability [START_REF] Simon | Instability of a partially ionized plasma in crossed electric and magnetic field[END_REF][START_REF] Hoh | Instability of pening-type discharge[END_REF], occurs when the electric field is collinear with the density gradient (fig. 3.1):

E • ∇n > 0.
(3.60) This instability is intrinsically related to the gradient-drift (anti-drift) eigen-mode (3.43). The most transparent insight on the destabilization mechanism can be obtained by looking at the ion and electron perturbations that were obtained above for the anti-drift mode as follows.

n e n 0 = ω * ω e φ T e (3.61)

n i n 0 = k 2 c 2 s ω 2 e φ T e (3.62)
These two expressions , together with qusineutrality, give the equation (3.43). Now, we add the equilibrium electric field to this configuration. The electron response is modified due to the Doppler shift induced by the E × B drift [START_REF] Frias | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of e x b discharges[END_REF], ω → ω -ω 0 , while the unmagnetized ions do not feel the magnetic field and therefore no effect of ω 0 . The modification of the electron response in (3.61), with (3.76) give the dispersion relation (3.69), discussed in more details below. Addition of the equilibrium ion flow results in the so called modified Simon-Hoh instability [START_REF] Sakawa | Excitation of the modified simon-hoh instability in an electron-beam produced plasma[END_REF][START_REF] Sakawa | Nonlinear evolution of the modified simon-hoh instability via a cascade of side-band instabilities in a weak beam-plasma system[END_REF]. The modified ion response is obtained from ion dynamics, considered cold and unmagnetized ions, as below :

∂n ∂t + ∇ • (nv i ) = 0 (3. 63
)
∂v i ∂t + (v i • ∇)v i = - e m i ∇φ (3.64)
This gives the expression for the perturbed density ion density in the form:

ñ n 0 = k 2 (ω -k • v i,0 ) 2 e φ m i . (3.65)
The E × B drift mainly determines the density evolution for the electrons. Considering the simplest case where there is no inertia and a response of the electron only due to the E × B flow, one may get the following linearized continuity equation :

∂ ñ ∂t + v E,0 • ∇ñ + ṽE • ∇n 0 + n 0 ∇ • ṽE = 0, (3.66)
where v E = E × B/B 2 . Considering the perturbations as in (3.2), (3.66) leads to :

ñ n 0 = ω * ω -ω 0 e φ T e , (3.67)
with ω 0 = k • v E,0 . Using equations (3.65) and (3.67), one can obtain the following dispersion relation :

ω * ω -ω 0 = k 2 c 2 s (ω -ω i,0 ) 2 (3.68) in which ω i,0 = k • v i,0
. In a simpler case, where we do not take into account an ion flow, (3.68) becomes simply :

ω * ω -ω 0 = k 2 c 2 s ω 2 , (3.69)
which was eluded to with simple arguments from (3.61,3.62).

One can see from (3.69), that for ω ω 0 , the necessary condition for the instability is

ω * ω 0 > 0. (3.70)
More generally from the expression of ω * and ω 0 , the Simon criterion, can be written as follows :

(k × b • g)(E × b • k) > 0. (3.71)
In the case where we consider a Cartesian plane in which b in the z direction, and E and g in the x direction then (3.71) simply becomes (3.60).

More accurate criterion can be derived 3.69, solving it as a second order polynomial. The the exact condition for the instability then is :

ω * ω 0 > k 2 c 2 s ω 2 * (3.72)
The expression (3.69) is known in the literature as the collisionless Simon-Hoh instability [START_REF] Smolyakov | Fluid theory and simulations of instabilities, turbulent transport and coherent structures in partially-magnetized plasmas of E x B discharges[END_REF], while (3.68) is known as the modified Simon-Hoh regime [START_REF] Sakawa | Excitation of the modified simon-hoh instability in an electron-beam produced plasma[END_REF]. For our applications, a Doppler shift in the ion response induced by the stationary ion velocity is present. This contribution modifies the real part of the frequency and for large values of the equilibrium ion velocity may lead to the emergence of a new unstable mode. The ions can be accelerated by an applied electric field which is balanced by friction due to ion collisions with neutrals.

The Simon-Hoh instability also exists in the collisional regimes as in the Simon's original paper [START_REF] Simon | Instability of a partially ionized plasma in crossed electric and magnetic field[END_REF]. Again, the instability is triggered by the relative drift between the ions and electrons due to the differences in the electron and ion mobility and diffusion. In this case, the ion and electron motion can be considered in drift-diffusion approximation , so the ion and electron currents in the equilibrium can be written as

Γ α,0 = -D α,⊥ ∇n 0 (x) ± n 0 µ α,⊥ E 0 (x), (3.73)
α = (i, e).
The perturbed currents are found in the next order ∼ O (ν/ω cα ), that require the addition of the Hall conductivity and mobility perpendicular to the magnetic field and to the direction of the perturbations,

Γ = -D α,⊥ ∇ n ± µ α,⊥ nE + b× -D α,H ∇ n ± µ α,H nE . (3.74)
The quasineutrality condition is written in terms of the current conservation equation

∇ • J α = 0 (3.75)
with J α = e α Γ α being the total current for a given species α. This approach leads to a complex form of the dispersion relation, different from (3.69); yet, the same criterion for the instability as in (3.60) appears for this case too.

General case: transition to the lower-hybrid and ion sound instabilities

The derivation for gradient-drift (or Simon-Hoh) instability neglected the effect of the inertia. One of the important consequences of this approximation is that the growth rate obtained from (3.69) increases indefinitely with the wave vector, γ ∼ k y . This is clearly nonphysical and such a model would not be suitable for numerical simulations as the maximal growth rate would be determined by the smallest resolution of the spatial grid. It turns out that the account of the electron inertia suppress the growth at short wavelengths and thus introduces physics based cutoff at short scales. We now give a brief outline of the appropriate dispersion relation. We also include electron-ion collisions connecting it to the collisional regimes of Simon-Hoh instability.

The ion response remains the same as in (3.62)

n i n 0 = k 2 c 2 s ω 2 e φ T e (3.76
which now, contrary to (3.47), includes the electron inertia effects on the left hand side. In the linear approximation we have approximated the total fluid (substantive) derivative

d dt = ∂ ∂t + v e • ∇ ∂ ∂t . (3.78
)

∂n α ∂t + ∇ • (n α v α ) = 0, (3.79)
The electron momentum equation (3.77) is solved perturbatively using the small parameters (ν e /ω ce , ω/ω ce ) < 1. Then one has for the electron velocity

v = cb B 0 × ∇ ⊥ φ - cT e en e b B 0 × ∇ ⊥ n + -iω + ν e,m ω 2
n n 0 = ω * + k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) ω -ω 0 + k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) eφ T e . (3.81)
Finally inserting this expression with the quasineutrality condition and (3.76) one obtains the dispersion relation

ω * + k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) ω -ω 0 + k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) = k 2 c 2 s ω 2 (3.82)
It is easy to see various interesting limits from the expression. In the long wavelength limit, k 2 ⊥ ρ 2 e → 0, one recovers the Simon-Hoh dispersion relation (3.69). In the opposite limit, short wavelength limit

k 2 ⊥ ρ 2
e → ∞, one has the stable ion sound mode

1 = k 2 c 2 s ω 2 (3.83)
In fact, ion sound appears as the finite temperature generalization of the lower -hybrid mode , namely

ω 2 = ω ce ω ci 1 + k 2 ⊥ ρ 2 e (3.84) For k 2 ⊥ ρ 2 e 1, ω 2 ω ce ω ci k 2 ⊥ ρ 2 e k 2 ⊥ c 2
s . This is the limit when the electrons become demagnetized and their response is simply Boltzmann from (3.81)

n n 0 eφ T e (3.85)
Lower-hybrid mode is recovered from (3.82) in neglect of the drift term, ω * → 0, then one has

k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) ω -ω 0 + k 2 ⊥ ρ 2 e (ω -ω 0 + iν e,m ) = k 2 c 2 s ω 2 (3.86)
With k = k ⊥ and ω 0 = 0 and ν e,m = 0 one obtains the stable lower-hybrid mode for k 2 ⊥ ρ 2 e < 1 and the ion sound for k 2 ⊥ ρ 2 e 1. The expression (3.86) shows that with ω 0 = 0 and ν e,m = 0 the lower

Conclusion

This chapter aimed to set the basis of the linear study of the plasma stability in partially magnetized plasmas. We gave the most basic definitions for waves and instabilities, overviewed the different eigenmodes relevant to plasmas of our interest, discussed mechanisms of destabilization, and finally considered the resulting linear instabilities in more details. The following chapter will deal with the linear analysis as applied to our fluid code MAGNIS. We will derive general dispersion relation appropriate to the code, aiming to check the linear simulation results. We will show that we can recover the unstable modes described in this chapter and study further their behaviour. Our analytical results from the dispersion relation will be used to verify the MAGNIS results when it is run under the conditions of the linear approximations.

Chapter 4

Linear analysis in MAGNIS

In the previous chapters, we introduced some basics of waves and instabilities in plasmas and specifics of our fluid code MAGNIS, a self-consistent multi-species and quasi-neutral fluid model. Our goal is to use MAGNIS to identify and characterize relevant instabilities and determine various consequences such as anomalous transport, plasma confinement and formation of nonlinear structures in conditions of some plasma sources. The first step of such analysis is to identify and verify linear instabilities in MAGNIS.

The instabilities and dispersion relations in Chapter 3 were derived by considering linear perturbations of a given equilibrium state of the plasma, defined by an equilibrium electric field, plasma density, plasma density gradient, ion flow velocity, and so on. This equilibrium state is simply assumed to exist and to be fixed and constant, so that the linearized equations have constant coefficients. However, in fully selfconsistent plasma models like MAGNIS, the equilibrium state is not fixed but is self-consistently calculated from the system of model equations. Therefore, in order to make a proper linear analysis of MAGNIS, it is essential to consider an equilibrium state that corresponds to a proper equilibrium solution of the MAGNIS equations. This leads to a number of complications. First of all, the self-consistent equilibrium solution evolves as a result of sources (ionization) and sinks (losses) and transport processes. For example, a plasma density gradient does not just exist but is caused by ionization (e.g. in the center) and losses (e.g. at the walls). Yet, all the analysis in Chapter 3 was done neglecting the source terms. Furthermore, the self-consistent equilibrium solution of MAGNIS tends to have inhomogeneous profiles for many variables such that the linearized perturbation equations do not have constant coefficients. This implies that the linear dispersion relation analysis of Chapter 3 can only be applied locally in the so-called Boussinesq approximation, neglecting the inhomogeneity of the equilibrium state on the scale of the wave length, so that only very small wave lengths can be described ( gradient length scale) and compared with the simulation results. Another complication is that the equilibrium solution in MAGNIS is usually found by simulating the evolution of the plasma from some initial state until it reaches equilibrium. During this evolution also the linear instabilities start to grow and it may be difficult to distinguish them from the overall plasma evolution. Moreover, the instabilities may quickly enter into the nonlinear regime and generate anomalous transport, thereby changing the equilibrium dynamics.

Thus, there are two questions: (a) how to describe and characterize linear instabilities in case of inhomogeneous profiles generated by ionization and sink sources; (b) how to set up MAGNIS in such a way that it is easy to observe and analyze the simulated instabilities in the linear regime, so that they can be compared with linear theory.

In this chapter, we describe a novel approach to generalize the theory of linear instabilities to a simple yet representative case of inhomogeneous equilibrium. This simple inhomogeneous equilibrium case can be simulated self-consistently with MAGNIS and at the same time, its linear perturbations can be described by a dispersion relation without the need for the local Boussinesq approximation. Thus, we derive here a modified linear dispersion relation that takes into account the inhomogeneities of the equilibrium. Secondly, we describe a technique that we use to force MAGNIS to run in the linear regime so we can compare it with linear theory. We analyze the modified dispersion relation and show that we recover the instabilities referenced earlier in Chapter 3. Then, we analyze the numerical results run in the linear regime and show that numerical diagnostics of frequencies and growth rates give results consistent with the results obtained from our dispersion relation.

Linear dispersion relation modified with the account of inhomogeneous equilibrium profiles

Here we present a derivation of the modified dispersion relation that include plasma sources (sinks), equilibrium electric field and density gradients. A simplified geometry of a plasma source will be presented, along with the boundary conditions, used the linear analytical and numerical studies. The linear analysis will be based on the continuity and momentum equations, the energy equation is not taken into account (T e and T i are constant and uniform), the magnetic field is also uniform. MAGNIS will be also runs under the same conditions.

Geometry

As part of our linear study, we set a simple model geometry in the 2D plane perpendicular to the magnetic field in Cartesian coordinates (x, y), where an electric field is applied in the x direction and the E × B drift is along the y direction, which is periodic. By analogy with cylindrical plasma devices such as Hall thrusters, we will sometimes refer to the x direction as the axial direction and to the (periodic) y direction as the azimuthal direction, but our coordinates are really Cartesian. The following boundary conditions are used for the axial direction of a length L x . A fixed potential difference applied is in the axial direction φ(0, y) = 0 (4.1)

φ(L x , y) = V (4.2)
and different densities are fixed at the ends

n e (0, y) = n i (0, y) = n L (4.3) n e (L x , y) = n i (L x , y) = n R . (4.4)
The E × B drift is closed in the azimuthal direction and periodic boundary conditions are applied for all variables and derivatives, u = (n, φ), u(x, 0) = u(x, L y ) (4.5)

∂ m u(x, 0) ∂y m = ∂ m u(x, L y ) ∂y m , ∀m = (0, ..., M ). ( 4.6)
This simple model geometry will be used for the set up of our dispersion relation and the simulations in MAGNIS.

Plasma equilibrium

As it was noted above the equilibrium in MAGNIS is set self-consistently by ionization and sinks. Here we describe the base equilibrium state in the model geometry defined above, using the equations that constitute the MAGNIS model. The full continuity equations for ions and electrons are

∂n ∂t + v α • ∇n + n∇ • v α = S α -L α , (4.7)
where α = (e, i) is a general index referring to both electrons and ions. Quasineutrality is assumed so n e = n i = n. Following the MAGNIS equations in Chapter 2, the source terms S α are due to ionization such that S e = S i = ν iz n (4.8)

with ν iz the ionization frequency. The sink terms L α represents losses by parallel transport which are proportional to n but may be different for each species:

L e = ν e n = L i = ν i n, (4.9)
where we introduced parallel loss frequencies ν α . Thus, we may have two different net source terms, and we will see that this is necessary for the stationary profiles to exist.

The stationary continuity equations are

v e,0 • ∇n 0 + n 0 ∇ • v e,0 = (ν iz -ν e )n 0 (4.10) v i,0 • ∇n 0 + n 0 ∇ • v i,0 = (ν iz -ν ,i )n 0 (4.11)
The key simpflication is the assumption that the equilibrium electric field is uniform; we will show below that this leads to the the conditions that equilibrium ion and electron velocity, v i,0 and v e,0 , are uniform as well, so ∇ • v e,0 = ∇ • v i,0 = 0. Thus the exponential density and a linear potential (Fig. 4.1) profiles are the equilibrium solutions:

n 0 = N 0 exp(g • x) (4.
12)

φ 0 = -E 0 • x. (4.13)
It is useful then to define a constant density gradient vector

g = ∇n 0 n 0 (4.14)
The velocities of the ions and electrons are related to the sources from (4.10) and (4.11) : The stationary velocities are from the total momentum equations

v e,0 • g = ν iz -ν e (4.15) v i,0 • g = ν iz -ν i (4.16) hence (v i,0 -v e,0 ) • g = ν e -ν i (4.17
m a ∂v α dt + (v α • ∇)v α + ν α v α = q α (E + v α × B) -eT α ∇n n (4.18)
In steady state (and with ∇v α,0 = 0) one has

m α ν α v α,0 = q α (E 0 + v α,0 × B) -eT α ∇n 0 n 0 (4.19)
Multiplying (4. [START_REF] Choueiri | Anomalous resistivity and heating in current-driven plasma thrusters[END_REF]) by B gives the following expression :

m α ν α v α,0 × B = q α (E 0 × B -B 2 v α,0 ) -eT α ∇n 0 n 0 × B (4.20)
Combining (4. [START_REF] Choueiri | Anomalous resistivity and heating in current-driven plasma thrusters[END_REF]) and (4.20) to get rid of the v × B term, one may get the following expressions

     v e,0 = -µ ⊥e (E 0 + T e g) + µ ⊥e (E 0 + T e g) × b v i,0 = µ ⊥i (E 0 -T i g) + µ ⊥i (E 0 -T i g) × b (4.21)
where

µ ⊥α = µ α 1 + h 2 α (4.22)
and we remind that g = ∇n 0 /n 0 is the density gradient, h α = eB/mν α is the Hall parameter, µ α = e/m α ν α the classical mobility, and b = B/B is a unit vector along the direction of B. In the simple model geometry defined at the beginning of this chapter, E 0 and g are both directed along the x axis, E 0 = E 0 x and g = gx, and the magetic field along the positive z axis, b = ẑ, in the Cartesian coordinate system of MAGNIS, so that

     v e,0 = -µ ⊥e (E 0 + T e g) x -µ ⊥e (E 0 + T e g) ŷ v i,0 = µ ⊥i (E 0 -T i g) x -µ ⊥i (E 0 -T i g) ŷ (4.23)
Since MAGNIS does not separate the equilibrium and perturbations, we have had to build some sort of equilibrium for which the relatively simple analytical solutions are possible in which the various instabilities can be analyzed and treated separately. The exponential density profile allows us to do just that and then compare with numerical results. Thus, expressions (4.21) describe non-ambipolar (v e,0 = v i,0 ) equilibrium state of plasma with non-uniform plasma density. Such solutions, with exponential density profile, constant velocities and electric field as given by (4. 21, 4.15, 4.16), are possible due to different source terms for the ions and electron as shown in equation (4.17). Although these conditions may not be valid for many practical situations (in particular, the assumption of the constant electric field), it allows us to consider a non-ambipolar cases with current flow and develop analytical solutions for the linear dispersion relation which then can be directly compared with MAGNIS solutions. Thus we are able to verify the capability of MAGNIS to recover the linear stage of the instabilities.

Alternatively, including the full profiles effects (density and temperature gradients, non constant fields and velocities, etc), leads to the full global non-local problem, [START_REF] Raitses | Structure of nonlocal gradient-drift instabilities in hall exb discharges[END_REF]. Such analysis, in the full non-local case, is difficult. Most of the time, it is solved numerically, the picture of the instabilities is complex and difficult for the comparison with the MAGNIS numerical results.

Perturbed equations and final general dispersion equation

In this part, we present the linearized equations and show how the inhomogeneous equilibrium density can be accounted in the linear dispersion relation. The temperatures and the magnetic field are assumed uniform and constant in time. The perturbed continuity equations 4.7 become

∂ ñ ∂t + ∇ • (n 0 ṽα ) + ∇ • (ñv α,0 ) = Sα -L α = (ν iz -ν α )ñ (4.24)
where the perturbed source and loss terms are Sα = ν iz ñ and L α = ν α ñ, while the frequencies ν iz and ν α for ions and electrons related to the equilibrium flows by the expressions (4.15) and (4.16).

The linearization of the total momentum equations (4.18) give for the electrons and ions

           α 1 ∂ṽ e ∂t + α 2 v e,0 • ∇ṽ e + ν e ṽe = - e m e (-∇ φ + ṽe × B) - eT e m e ∇(n 0 + ñ) (n 0 + ñ) - ∇n 0 n 0 β 1 ∂ṽ i ∂t + β 2 v i,0 • ∇ṽ i + ν i ṽi = e m i (-∇ φ + β 3 ṽi × B) - eT i m i ∇(n 0 + ñ) (n 0 + ñ) - ∇n 0 n 0 (4.25)
where ν α is a total momentum frequency given by:

ν α = ν iz + ν m,α (4.26)
with ν m,αn being the momentum transfer frequency for collisions with the neutral gas.

The last term of both equations is linearized by a Taylor expansion such that it gives

eT m ∇ñ n 0 - ñ n 0 ∇n 0 n 0 (4.27)
The final linearized momentum equations are:

           α 1 ∂ṽ e ∂t + α 2 v e,0 • ∇ṽ e + ν e ṽe = - e m e (-∇ φ + ṽe × B) - eT e m e ∇ñ n 0 - ñ n 0 ∇n 0 n 0 β 1 ∂ṽ i ∂t + β 2 v i,0 • ∇ṽ i + ν i ṽi = e m i (-∇ φ + β 3 ṽi × B) - eT i m i ∇ñ n 0 - ñ n 0 ∇n 0 n 0 (4.28)
We have added the controlling factors before inertial terms (α 1 , α 2 , β 1 and β 2 ) for both equations and before the magnetization term for ions (β 3 ) so that we can investigate the roles of these effects separately. We now note an important property of the linearized equations (4.28). As usual in the linear analysis we consider the Fourier perturbations near the stationary state:

η = η 0 + η, (4.29)
and η = η 1 exp(-iωt + ik • x). (4.30)
in which η 1 is constant. From the structure of equations (4.28) it is obvious that the proper form of the potential, density and velocity perturbations are the coefficient defining the amplitude of the perturbation.

In our case, referring to the explanations given for the choice of different source terms, if we want this coefficient to remain constant for all fluctuating variables of the system all while being consistent with the linearized system above, one may see that it is essential to consider an exponential stationary density, along with a fluctuating density under the following form:

φ, ṽα , ñ n 0 (x) ∼ e -iωt+ik•x (4.31)
Thus we write ṽα = v α,1 e -iωt+ik•x (4.32)

ñ n 0 (x) = z 1 e -iωt+ik•x (4.33) φ = φ 1 e -iωt+ik•x (4.34)
These give the expressions for the fluctuating velocity components

             v e,1 = - e m e [ e me ik × B -ikν e,ef f ] ν 2 e,ef f + ω 2 c,e
(φ 1 -T e z 1 )

v i,1 = - e m i [β 3 e mi ik × B + ikν i,ef f ] ν 2 i,ef f + β 2 3 ω 2 c,i (φ 1 + T i z 1 ) (4.35)
where ν α,ef f is an effective frequency proper to each species

ν e,ef f = ν e -α 1 iω + iα 2 v e,0 • k ν i,ef f = ν i -β 1 iω + iβ 2 v i,0 • k (4.36)
After inserting the velocities (4.35) into the linearized continuity equations, one may obtain the following equations:

-z 1 iω + v e,1 • (g + ik) + z 1 v e,0 • (g + ik) = (ν iz -ν e )z 1 -z 1 iω + v i,1 • (g + ik) + z 1 v i,0 • (g + ik) = (ν iz -ν i )z 1 (4.37)
and using the equilibrium conditions (4.15) and (4.16) , one obtains:

                 z 1 iω + [ e 2 m 2 e ik × B -ik e me ν e,ef f ](φ 1 -T e z 1 ) ν 2 e,ef f + ω 2 c,e • (g + ik) -z 1 v e,0 • ik = 0 z 1 iω + [β 3 e 2 m 2 i ik × B + ik e mi ν i,ef f ](φ 1 + T i z 1 ) ν 2 i,ef f + β 2 3 ω 2 c,i • (g + ik) -z 1 v i,0 • ik = 0 (4.38)
Finally, eliminating the perturbation amplitudes, we find the dispersion relation:

(iω -v e,0 • ik)(ν 2 e,ef f + ω 2 c,e ) e me (ω c,e g • ik × b -ν e,ef f (g • ik -k 2 )) - (iω -v i,0 • ik)(ν 2 i,ef f + β 3 ω 2 c,i ) e mi (β 3 ω c,i g • ik × b + ν i,ef f (g • ik -k 2 )) = T i + T e , (4.39)
which is in fact a fourth order polynomial equation in ω (to see this, realize that ω is included in ν α,ef f ).

As it was discussed, this general dispersion relation takes into account the inhomogenous equilibrium that can be generated and mantained within the full fluid equations employed in MAGNIS. It also has most of the physics of unstable modes discussed in the previous chapter and expected to be active in MAGNIS as it will be outlined in the following section.

Linear instabilities from the generalized dispersion equation

In Chapter 3, we have discussed several unstable modes that can be relevant to plasma sources of interest in conditions of the local approximation. In the previous section, we built a more general linear dispersion relation that takes into account the effects inhomogeneous equilibrium state that can be reproduced in MAGNIS. Now in this section, we show that the typical instabilities are defined in Chapter 3 are all included in our dispersion relation. We investigate these unstable modes independently and then as a combination of both gradient and ExB drift modes, as functions of the key parameters. Thus, we propose to investigate our dispersion relation following the scheme in Figure 4.3 in order to find the unstable modes defined in Chapter 3.

Case 1 : a gradient type instability, the classical Simon-Hoh

We defined the Simon-Hoh instability in section 3.3.2, characteristics, equations and dispersion relation as commonly seen in literature. However, it is possible to recover this instability with our general dispersion relation as defined in fig. 4.3, meaning without the ion mass but with collisions as in Simon's initial publication [START_REF] Simon | Instability of a partially ionized plasma in crossed electric and magnetic field[END_REF], with the advantage to have the axial contribution included, so that we focus only on the gradient driven mode to see its behaviour and effects.

In this configuration, if we refer to our dispersion relation, we do not take into account the inertial terms, and we remind that ions are unmagnetized. Then, once these terms neglected, the dispersion relation related to this instability becomes:

(iω -v e,0 • ik)(ν 2 e + ω 2 c,e ) e me (ω c,e g • ik × b -ν e (g • ik -k 2 )) - (iω -v i,0 • ik)ν i e mi (g • ik -k 2 ) = T i + T e . (4.40)
Injecting the expression of the velocities defined in (4.23) and reorganizing the variables, one obtains:

(iω -(-µ ⊥e (E + T e g)(x + h e ŷ) • ik) µ ⊥e (h e g • ik × b -(g • ik -k 2 )) - (iω -µ i (E -T i g)x • ik) µ i (g • ik -k 2 ) = T i + T e (4.41)
which, once renormalized, becomes:

Ω(1 + h 2 e ) + (1 + Σ)(K x + h e K y ) h e K y -(K x + iK 2 ) - Ω -ψK x (Σ -T ) ψ(K x + iK 2 ) = 1 + T. (4.42)
Here we used the following dimensionless parameters:

Ω = ω µ e T e g 2 , (4.43)
the dimensionless frequency, in which we normalized the frequency ω by an effective frequency that could be seen as an unmagnetized electron-diffusion frequency (≈ D e g 2 ),

T = T i T e , (4.44)
the temperature ratio,

h e = ω c,e ν e , (4.45)
the Hall parameter,

Σ = E 0 gT e , (4.46)
ratio between the electric field and the pressure,

ψ = µ i µ e , (4.47)
the ion-over-electron unmagnetized mobilities ratio, and finally

K = k g , (4.48)
the wave vector normalized by the density gradient g.

An interesting thing to study is the behaviour of wave frequency and growth rate of this instability as a function of k, the wave vector, first for some fixed values of the four other dimensionless parameters listed above. Figure 4.4 shows the dimensionless growth rate and real frequency in the (K x , K y ) plane, for the set of the dimensionless parameter values given in Table 4.1, in which we allowed negative values for the axial component of the wave vector but kept only the positive part for the azimuthal component given the symmetry of the system. The dimensionless growth rate enables to characterize the instability's growth, while the frequency gives information about the wave propagation characteristics; the velocity and its direction. One can see that the growth rate is maximum for intermediate values of K, and is stabilizing and decreasing progressively as K increases; moreover, it is interesting to see that the instability velocity strictly depends on the azimuthal wave vector K y .

ψ 3.44×10 -4 Σ 0.17
Equation (4.42) being a first order equation in Ω, it is straightforward to solve it for the total frequency, which leads to :

Ω = (Σ+1)(Kx+heKy) heKy+Kx+iK 2 -(Σ-T )Kx Kx+iK 2 -(1 + T ) (h 2 e +1) heKy+Kx+iK 2 + 1 ψ(Kx+iK 2 )
.

(4.49)

From (4.49), one can obtain the growth rate and the real frequency by isolating the imaginary and real parts, respectively. Once we obtain an expression for the growth rate, we can study it analytically to understand the observed trend in Figure 4.4; a first thing we are interested in is the behaviour of the instability as a function of K. If we investigate Ω considering the limit K → ∞, then one may get the following expression:

Ω(K → ∞) → -i ψh e (T + 1) ψ(1 + h 2 e ) + 1 K 2 . (4.50)
Equation (4.50) indicates that the growth rate becomes negative at high K and decreases as a function of K 2 , which explains the observed trend fig. 4.4 and confirms it in general for all possible values of the input parameters. Also, a stability criterion can be analytically deduced from this growth rate; indeed, as we explained in chapter 3, one has just to find the growth rate's threshold such that Ω i = (Ω) > 0, condition for the instability to appear. Once calculated, this stability criterion shows for this mode a link between the parameters Σ, ψ and h e , and is rigorously written as follows :

[ψ(-Σ + T )K x -(Σ + 1)(h e K y + K x )]h e K y + (T + 1)[(ψ + 1)(K 4 + K 2 x ) + h e K y (h e K y -2K x )] < 0. (4.51)
A representation of this stability criterion depending on the different dimensionless parameters listed above is given in Figure 4.5. We choose to represent the variation of the growth rate over the (h e , Σ) plane for a fixed value of the ψ parameter. An interesting fact to notice is that the growth rate remains positive even in the case where Simon's instability criterion E 0 g > 0 (see equation (3.60)) is not verified ( Σ < 0 in Figure 4.5). This discrepancy can be simply explained: in our approach, we allowed for all existing k x along with all k y whereas Simon in fact established this criterion in a configuration in which the axial direction contained profiles for the velocities, density and electric field and thus, required to be averaged in order to obtain a dispersion relation. Consequently, this criterion applies only for waves in the transverse direction (y direction) and we reproduced it in Figure 4.6, which is similar to Figure 4.5 except that did not allow for axial wave numbers, i.e. we imposed K x = 0. Here, one may see that this time, Simon's criterion is verified. Similarly, we also show that Simon's initial instability criterion can be recovered with our own criterion given in (4.51); indeed, for that purpose, we just have to remove the axial contribution of our expression, which leads to: Expression (4.54) is equivalent to Simon's criterion for the weak potential field limit; this criterion is satisfied when k y is small, meaning large wavelength, and T i g low compared to E.

-Σh 2 e K 2 y + [(T + 1)(ψ + 1)K 4 y + T h 2 e K 2 y ] < 0. ( 4

Case 2 : a drift type instability, Farley-Bunemann

In this case, we remind that we do not take into account any plasma inhomogeneity (g = 0) and the electron inertial terms are considered negligible (α 1 = α 2 = 0). The ion-neutral collision frequency is considered greater than the ion cyclotron frequency (ν i ω c,i ), which makes the latter negligible and as a result, ions are unmagnetized (β 3 = 0), and finally we consider them inertial (β 1 = β 2 = 1). Considering all these conditions, then, our general dispersion relation becomes a simple second order polynomial equation in ω:

m e m i (iω -v e,0 • ik) (ν 2 e + ω 2 c,e ) ν e + (iω -v i,0 • ik)(-iω + v i,0 • ik + ν i ) = c 2 s k 2 (4.55)
where we choose to write c s = e(T i + T e )/m i . We rewrite equation (4.55) as follows:

m e m i (i ω -v ei • ik) (ν 2 e + ω 2 c,e ) ν e + i ω(-i ω + ν i ) = c 2 s k 2 (4.56)
with ω = ω -v i,0 • k the frequency in the ion frame, and v ei = v e,0 -v i,0 , the relative velocity. If we divide (4.56) by ν 2 i , we obtain:

1 ν i (i ω -v ei • ik) m e /m i ν i ν e /(ν 2 e + ω 2 c,e ) + i ω ν 2 i (-i ω + ν i ) - c 2 s ν 2 i k 2 = 0, (4.57)
so that after some rewriting, the dimensionless dispersion relation becomes:

ψ Ω 2 + (ψ + 1)i Ω -iV ei • K -ψK 2 = 0 (4.58) where ψ = ν e ν i (ν 2 e + ω 2 c,e ) me mi (4.59)
is the ratio between magnetized electrons and ions mobilities,

V ei = v ei c s (4.60)
is the normalized relative velocity,

K = c s k ν i (4.61)
is the normalized wave vector, in which c s /ν i can be seen as an averaged ion mean free path, and finally,

Ω = Ω + V i • K = ω ν i + v i,0 • k ν i (4.62)
defines the normalized frequencies in the laboratory and ion frames. Figure 4.7 displays a first result for this unstable mode in the normalized (K x ,K y ) plane. As the value for the input parameters, we fixed ψ, while V ei is calculated self-consistently such that:

V ei = - E 0 c s (T e + T i )ν i ((1 + ψ)x + h e ψŷ). (4.63)
Also,

V i,0 = E 0 c s (T e + T i )ν i x. (4.64)
We show the values we chose for the input parameters in Table 4 In Section 3.3.1 of Chapter 3, we presented the Farley-Buneman instability and showed that in literature, the growth rate is given such that (see equation (3.57)):

Ω i = ψ 1 + ψ (K 2 -Ω 2 r ). ( 4

.65)

This growth rate, as already mentioned in the said section, based on the approximation that Ω i (Ω r , ν i ), implies that the instability grows indefinitely with K 2 ; however, one may see in Fig. 4.7 that expression (4.65) is not valid for high K. Indeed, if the assumption made to facilitate the calculation of the growth rate in the case of the ionosphere is correct most of the time, it true is in our case only for a very small part of the K domain of interest. If we solve the polynomial (4.58) analytically without approximations, we obtain two roots as follows:

Ω = (ψ + 1)i 2ψ ± 1 2ψ 4ψ(ψK 2 + iV ei • K) -(ψ + 1) 2 , (4.66)
and isolating the imaginary part of the most unstable root, we get:

Ω i = ψ + 1 2ψ + 1 2 (4ψV ei • K) 2 + [4ψ 2 K 2 -(ψ + 1) 2 ] 2 - 1 2 [4ψ 2 K 2 -(ψ + 1) 2 ]. (4.67)
Expression (4.67) is the general growth rate for the Farley-Bunemann instability without any approximation made. If we now consider the high K limit of (4.66), we get

Ω(K → ∞) = ±K + i 2ψ ± V ei • K K -(ψ + 1) . (4.68)
The most unstable root is then given by:

Ω i (K → ∞) = i 2ψ |V ei • K| K -(ψ + 1) (4.69) Ω r (K → ∞) = ±K + V i0 • K, (4.70)
where the ± sign in the expression for Ω r corresponds to the sign of (V ei • K). This analytical solution shows that the growth rate Ω i of this mode tends asymptotically towards a finite value at high K. In Figure 4.7 one can see the growth rate is indeed increasing to reach a constant value on the (K x , K y ) plane. In fact, this case presents a stabilization issue at high K, as the growth rate is always positive and does not decrease when K increases so that instabilities keep arising; if K is too high, this means the wavelength attached to this mode will be too small for the validity of the model, which is problematic in fluid models as we are limited to certain scales (in general for our plasma of interest, it should not be under the electron Larmor radius) as mentioned in Chapter 3. Usually, it is common to damp the growth rate and thus, to stabilize this mode by implementing a diffusive term in the ion momentum equation representative of an ion viscosity (we discuss this matter later, in Chapter 6). Equation (4.70) explains variation of the real frequency Ω r in Figure 4.7; this is almost independent of K y and proportional to K x because V i0 1 in this case, meaning that the ions are supersonic. According to the values of Table 4.2, the normalized ions speed (Mach number) is 4.09, indeed much greater than 1. In fact, this equation predicts (4.70) the following wave propagation velocity:

v wave = ±c s + v i0 • k k (4.71)
in the direction of the wave vector. Note that for the shown plane, the minus sign should apply in the first term.

The same stability criterion found in Chapter 3 can be recovered from expression (4.66), one finds that ω i > 0 if:

|V ei • K| K > 1 + ψ ⇔ |v ei • k| k > (1 + ψ)c s , (4.72)
equivalent to (3.59). We show in Figure 4.8 the normalized growth rate, real frequency and wave vector in the (V ei , ψ) plane; the figure clearly illustrates the boundary between the unstable and stable regions described by (4.72). We also note that the above expressions (4.67, 4.69, 4.72) show that the most unstable wave vector is always directed along the relative velocity:

K = ±K V ei V ei ⇔ |V ei • K| = V ei K. (4.73)
Using the expression of the relative velocity, (4.63), one finds that for this optimal wave direction:

K x K y = 1 + ψ h e ψ . (4.74)
This tilted direction can clearly be seen in Figure 4.7.

Case 3 : a gradient-drift type instability, a combination of the first two cases

As mentioned in fig. 4.3, we call the gradient-drift instability the coupling of the Farley-Bunemann instability (Case 2) with a density gradient. This mode could also be seen as our collisional Simon-Hoh instability (Case 1) in which we add the ion inertial terms. Referring to our dispersion relation (4.39), one may only remove the electron inertial terms (α 1 = α 2 = 0, which leads to:

(iω -v e,0 • ik)(ν 2 e + ω 2 c,e ) e me (ω c,e g • ik × b -ν e (g • ik -k 2 )) - (iω -v i,0 • ik)ν i,eff e mi (g • ik -k 2 ) = (T i + T e ). ( 4

.75)

If we develop ν i,eff and pose ω = ω -v i,0 • k, then one can rewrite (4.75) as follows:

(i ω -v e,i • ik)(ν 2 e + ω 2 c,e ) e me (ω c,e g • ik × b -ν e (g • ik -k 2 )) + (i ω + ν i )i ω e mi (g • ik -k 2 ) = (T i + T e ), (4.76)
and reorganizing the different variables in (4.76) enables to rewrite it again:

i ω -v e,i • ik miνe me(ν 2 e +ω 2 c,e ) (h e g • ik × b -(g • ik -k 2 )) + (i ω + ν i )i ω (g • ik -k 2 ) = e(T i + T e ) m i . (4.77)
From (4.77), it is then easy to obtain a normalized dispersion relation:

Ω -V ei • K ψ h e G•K×b G•K+iK 2 -1 + (i Ω + 1) Ω -(G • K + iK 2 ) = 0 (4.78)
where the dimensionless parameters are:

ψ = m i m e ν i ν e ω 2 c,e + ν 2 e , (4.79)
which is the same electron-over-ion ratio mobility as in (4.59),

V ei = v ei c s , (4.80)
the relative velocity normalized as in (4.60),

G = c s ν i g, (4.81)
the normalized density gradient, and finally, the normalized wave vector and growth rates as given in (4.61) and (4.62):

K = c s ν i k, (4.82
)

Ω = Ω + V i • K = ω ν i + v i,0 • k ν i . (4.83)
The normalized dispersion relation (4.78) is again second order in Ω and can be written as:

a 2 Ω 2 -a 1 i Ω + a 0 = 0 (4.84)
with the coeficients

a 2 = ψ, (4.85
)

a 1 = ψ - 1 h e ξ -1 , (4.86)
a 0 = V ei • iK h e ξ -1 + ψ(iK • G -K 2 ). (4.87) where ξ = K × b • G K • G + iK 2 . (4.88)
Equation ( 4.84) has the benefit to be analytically solved just like the Farley-Bunemann case. If we proceed to find the roots, we get:

Ω = i 1 2ψ ψ - 1 ξh e -1 ± 1 2ψ ψ - 1 ξh e -1 2 + 4ψ V ei • iK ξh e -1 + ψ(iK • G -K 2 ) . (4.89)
Seeking the limit for high K (1, G) leads to the following two roots:

Ω(K → ∞) = ±K + i 2ψ ± (V ei -ψG) • K K -(ψ + 1) , (4.90)
where the most unstable root is given by:

Ω i (K → ∞) = i 2ψ |(V ei -ψG) • K| K -(ψ + 1) (4.91) Ω r (K → ∞) = ±K + V i0 • K. (4.92)
These high K expressions are similar to those in the previous section, the only difference being that the relative velocity V ei is replaced by (V ei -ψG). Therefore the same trend observed in the previous section is reproduced for this case; the growth rate tends asymptotically towards a constant value as K increases, which again, leads to high-k issues, since it shows no signs of damping. This trend can be observed in Figure 4.9, where we show the normalized growth rate and real frequency in the (K x , K y ) plane. One can also recognize in this figure the behavior of the real frequency Ω r as predicted in equation (4.92). Contrary to the case in the previous section, this is now dominated by the first term that is proportional to the norm of K, with a negative sign in our case (because V ei •K < 0 for the K plane shown in the figure).

In the previous chapter, we discussed in detail the case of the Simon-Hoh instability, and how this regime transits towards a lower-hybrid mode when the electron inertial terms are added to the system, which stabilizes the instability at high k and thus, cuts small scales. We now wish to investigate the effect of these inertial terms for our current case, which differs from the Simon-Hoh case by the presence of the collisional and equilibrium drift contributions for ions.

In this case, we show that the addition of the electron inertial terms (α 1 = α 2 = 1) in our initial dispersion relation (4.75) for this case results in:

(iω -v e,0 • ik)(ν 2 e,eff + ω 2 c,e ) e me (ω c,e g • ik × b -ν e,eff (g • ik -k 2 )) - (iω -v i,0 • ik)ν i,eff e mi (g • ik -k 2 ) = (T i + T e ) (4.93) ⇔ (iω -v e,0 • ik)((ν e -iω + v e,0 • ik) 2 + ω 2 c,e ) e me (ω c,e g • ik × b -(ν e -iω + v e,0 • ik)(g • ik -k 2 )) + (iω -v i,0 • ik)(iω -v i,0 • ik -ν i ) e mi (g • ik -k 2 ) = (T i + T e ).
(4.94) In order to analytically lighten expression (4.94), we will consider our case only in a low-frequency regime where all the frequencies in the system are negligible compared to the electron cyclotron frequency (|ω| ω c,e ) (see Section 1.4.3); thus, we can remove some terms and simplify (4.94) as:

(iω -v e,0 • ik)ω 2 c,e e me (ω c,e g • ik × b -(ν e -iω + v e,0 • ik)(g • ik -k 2 )) + (iω -v i,0 • ik)(iω -v i,0 • ik -ν i ) e mi (g • ik -k 2 )
= (T i + T e ) (4.95) which can be written as:

(i ω -v ei • ik) e meω 2 c,e (ω c,e g • ik × b -(ν e -i ω + v ei • ik)(g • ik -k 2 )) - i ω(i ω -ν i ) e mi (g • ik -k 2 ) = (T i + T e ) (4.96)
Finally, rewriting this equation with normalized variables, as before, one obtains the normalized dispersion relation for the gradient-drift instability with electron inertia included:

( Ω -V ei • K) ψ ξh e -1 + ν N (i Ω -V ei • iK) + (i Ω -1) Ω -(G • K + iK 2 ) = 0 (4.97)
where Ω, K, V ei , ξ, h e and G are the same quantities defined above in (4. Comparing equation (4.97) with the dispersion relation (4.78) without electron inertia, we see that the only difference lies in the term multiplied by ν N in the denominator on the left, consistening of two parts, i Ω and -V ei • iK, or in the laboratory frame, iΩ and -V e0 • iK, coming from the temporal and spatial inertia terms of the electron momentum equation. We now wish to investigate the role of these terms and see if they are sufficient to stabilize the large wave numbers, thereby making the model suitable for numerical simulations. Note that in a numerical code like MAGNIS, these terms are handled very differently (see Chapter 2) and the ∂v e /∂t term is much more easy to include than the (v e • ∇)v e term. Furthermore, according to the gyroviscous cancellation, the (v e • ∇) used above and in MAGNIS should actually be (v E • ∇) without diamagnetic drift (see 1.118).

Therefore, we first look at the effect of the time derive term alone. Figure 4.10 shows the results for the normalized growth rate and real frequency in the (K x , K y ) plane, taking into account the effects of ∂v e /∂t alone; one may easily see that for the given conditions, the time derivative term works very well to stabilize the growth rate at high K. Moreover, a close look at the figure enables to see that the general dynamics of the system is preserved; indeed, in comparison with Figure 4.9, the maximum and minimum values for the real frequency and growth rate remain the same along with the real frequency trend in K.

However, the stabilizing effects of the time derivative term alone is limited to a narrow range of parameters. Indeed, for a configuration identical to that in Figure 4.10 but a larger value of the electric field (E 0 = 200 V/m), Figure 4.11 shows a growth rate that never stabilizes at high K. We then proceed to add the electron equilibrium velocity contribution via the (v e • ∇)v e term, meaning that we now include the full electron inertial terms. It turns out that in our current case of interest, this stabilizes the growth rate at high K all the time, as illustrated by Figure 4.12. We may therefore conclude that all the contributions in the electron inertia are required to avoid stabilization issues at high K if one wants to consider a larger range of values of the different input parameters.

Comparisons with MAGNIS 4.3.1 Setting up MAGNIS for the simple model conditions

In the previous section, we investigated all the possible unstable modes present in our dispersion relation, and thus, instabilities that can appear in our simple model. In this section, we proceed to make a verification of MAGNIS by running simulations of these cases to make sure that the code is numerically capable to return the solutions of the physical model. For this purpose, MAGNIS is first set in the simple model configuration defined in Section 4.1.1 in which:

• The considered simulation domain is exactly as in Figure 4.1, along with the exact same boundary conditions listed (4.1 -4.6), meaning a periodic y direction with a constant electric field imposed in the x direction and resulting in the closed E × B drift in the y direction.

• The electron energy equation is switched off.

• Given an electric field E 0 , density gradient g, electron temperature T e , ion temperature T i and plasma density value at the left boundary n L = N 0 , we fix the applied voltage V = E 0 /L x , the density at the right boundary n R = N 0 exp(gL x ), as well as the following ionization and loss frequencies:

ν iz = µ i (E 0 -T i g)g, ( 4
.99)

ν i = 0 ⇒ ν e = ν iz + µ ⊥e (E 0 + T e g)g, ( 4
.100) such that the equilibrium equations (4.15, 4.16, 4.17) are satisfied. All these parameters are kept fixed during the simulations, and together they generate and maintain the desired stationary solution.

• As an initial condition, we use the equilibrium profiles: n = N 0 exp(gx), φ = E 0 x, and the constant electron and ion velocities from equation (4.23).

A next and very important adaptation of MAGNIS for this study, is the implementation of a simple method to force the solution to stay in a linear instability regime at all time. This method works as follows. Thanks to the periodic boundary conditions in the y direction, it is possible to seperate the solution into a fluctuating and equilibrium part. The equilbrium part is simply the average of the solution over the y direction:

X 0 (x) = X(x, y) y , (4.101)

where X(x, y) represents every variable of the model, and so the fluctuating part can be calculated from X(x, y) = X(x, y) -X(x, y) y . (4.102)

A global amplitude of this fluctuating part is monitored continously during the simulation. Whenever this amplitude becomes too large, the whole fluctuating part of the solution, for all the variables, is rescaled by the same factor, and the full solution is reconstructed by adding the equilibrium part:

X(x, y) → X y (x) + α X(x, y) -X y (x) . (4.103)
where α < 1 is a rescaling factor for the fluctuations which must small enough to compensate for the overall fluctuation growth per time step, α < (γ∆t) -1 , but large enough that the fluctuations remain well resolved by the machine precision; we typically use α = 10 -3 . The monitoring of the fluctuation amplitude is done from the electric potential fluctuations, which we calculate every time step and compare with the electron temperature; the above fluctuation rescaling is done if

max |φ(x, y) -φ y |/T e > β (4.104)
where the maximum is taken over the full domain and β is a parameter defining the rescaling threshold, where the value must be small enough to keep the system in a linear regime, typically β = 10 -3 . This simple technique allows us to simulate the linear instabilities indefinitely, so that we have time to perform accurate diagnostics on them that can be checked against the dispersion relation. Note the method does not work for purely axial modes, some fluctuation component along y is necessary.

Implementation of the diagnostics tools

Lastly, the final modification consists in implementing the necessary tools in order to realize diagnostics of the simulation; in particular, we wish to deduce from the simulation results the growth rate, frequency and wave vector of the instabilities, in order to compare them with those given by the dispersion relation.

To obtain these quantities, we perform a Fourier transform (with a standard Fast Fourier Transform FORTRAN routine) on the spatial profile in the periodic y direction, of the azimuthal component of the electric field E y , characteristic variable of the emerge and development of instabilities. For a given point x and time t, the y profile of E y can be written as

E y (x, y, t) = Ny/2-1 n=1 A n (x, t) exp(ik y y) (4.105)
where A n (x, t) are complex coefficients given by the Fourier transform for all discrete wave numbers

k y = 2πn L y , (4.106)
and N y is the number of numerical grid points in y. We consider that

A n (x, t) = C n exp(ik x,n x -iω n t) (4.107)
where k x,n and ω n are complex quantities, respectively the axial wave vector and the total frequency for a given mode n. From the Fourier coefficients at two different time steps t and t + ∆t and in the same x 99 coordinate, it is then possible to extract a growth rate and a frequency:

A n (x, t + ∆t) = exp(-iω n ∆t)A n (x, t) ⇒ (4.108) ω n = i ∆t ln A n (x, t + ∆t) A n (x, t) (4.109)
where ω n is the total complex frequency gathering the frequency and growth rate of a particular unstable mode n as follows:

ω n = ω num r,n + iγ num i,n . (4.110)
where we added the "num" superscripts to indicate that these are values measured in the numerical simulations.

In order to obtain the corresponding axial wave number k x , the same logic is applied. Considering the Fourier transforms at two successive axial points x and x + ∆x for the same time t, we have

A n (x + ∆x, t) = exp(ik x ∆x)A n (x, t) ⇒ (4.111) k x,n = - i ∆x ln A n (x + ∆x, t ref ) A n (x, t ref ) . ( 4

.112)

Note that this axial wave number can in principle be complex, depending on the boundary conditions for the x direction. However, in the dispersion relation it is considered real (corresponding to a plane wave without boundaries), so in most of the comparisons below we just neglected the imaginary part. After a running the simulation for a sufficiently long time, there is usually a single mode n that dominates the profile of E y , the mode with the largest growth rate. It is the wave number, frequency and growth rate of this mode that we compare with the dispersion relation.

We described how MAGNIS is set in the specific simple model configuration, forced in a linear regime, and finally how are the numerical frequencies and growth rates are determined in this linear regime. In the following subsection, we describe how simulations are carried out in MAGNIS; we explain the cases we choose to run, and how we proceed to select for each case a correct set of input parameters. After that, we compare the frequencies and growth rates numerically measured in the simulations for each case with the ones predicted by the dispersion relation.

Comparisons : results

Section 4.2 presents the behaviour of the different linear instabilities that are expected to appear in our configuration. As a result of this study, we choose to select the case 1 and the case 3 split in two parts; first, with partial electron inertia that we choose to call case 3bis (a), in which we take into account only the time derivative term, and total electron inertia, that we name (b) respectively. We justify this choice thanks to our analytical study since we demonstrated that they do not cause issues at high k for a reasonable range of the input parameters. Note that we considered case 3bis (b) given the sensitivity of case 3bis (a) to the variation of the input parameters.

Then, we proceed to find thanks to our dispersion relation a reasonable set of input parameters for each case to enter in MAGNIS; for fixed values of the density gradient and temperatures, we browse a wide range of values for the electric and magnetic fields to find for which couple of (E, B) the most unstable roots can be properly resolved by MAGNIS, basing on the fact that k y L y must not be too large or too small, meaning we want a reasonable number of periods in the domain for fixed values of its dimensions L x and L y (note that the condition is meaningful only for the azimuthal direction since it is the only one which is periodic; however k x has an effect on the system that we detail later). An attentive observation of the resulting figures (an example is shown for case 1 in fig. 4.13) enables to identify a good pair of (E, B) values for the checking of MAGNIS; we choose for each studied case to limit the number of periods to five in the domain fixed by the input dimensions. For that purpose, our approach consists in looking for an unstable growth rate in the (E, B) plane for which the attached wave number k y does not give more than four periods in the domain for a given L y . The values we chose as input parameters (see figures in fig. 4.13) are reported in tables 4.4, 4.5 and 4.6, and we present for each of these cases the expected values given by the dispersion relation of the frequency, growth rate and their respective axial and azimuthal wave numbers for the most unstable root in table 4.7.

E (V/m) 125
B (G) 100
T e (eV) 1
T i (eV) 0.2 g (m -1 ) 50
L x (cm) 4 L y (cm) 12
N g (m -3 ) 10 20
E (V/m) 20
B (G) 40
T e (eV) 1
T i (eV) 0.2 g (m -1 ) 25
L x (cm) 4 L y (cm) 4 N g (m -3 ) 10 20
Table 4.5 -Input values for case 3bis (a).

E (V/m) 50
B (G) 20
T e (eV) 1
T i (eV) 0.2 g (m -1 ) 40
L x (cm) 4
L y (cm) 4

N g (m -3 ) 10 The expected number of periods in the domain for the azimuthal directions is given such that n p = L y k y,max /2π, with n p the number of periods and k y,max , the azimuthal wave number attached to the most unstable growth rate γ max and frequency ω r,max .

Cases ω r,max (s -1 ) γ max (s -1 ) k x,max (m -1 ) k y,max (m -1 ) n p case 1 -2.425×10 6 5.033×10

Mesh convergence and error assessment

We show the simulation results of MAGNIS for the azimuthal component of the electric field for each case fig. Each case displays the expected number of periods in the domain predicted by the dispersion relation for all meshes, except for case 3bis (a) where the right number of period appears only for the 512 mesh. Tables 4.8, 4.9 and 4.10 show all the values on each mesh for each case with an error assessment for each value of each mesh. The error percentage is calculated such that

δ σ = |σ th -σ num | |σ th | (4.113)
where σ th is the analytical quantity deduced from our dispersion relation and listed Globally speaking, all the cases show good mesh convergence with values very close to those predicted by the dispersion relation for each considered mesh grid. However, some interesting discrepancies can be pointed out, along with some remarks :

• In general, for the three cases, the more we refine the mesh, the more the error tends to decrease while the results accuracy increases; this is particularly true for case 3bis (a), in which one can clearly see that the results obtained on mesh 512 are closer to the analytical ones than the results given on the two other mesh. In fact, in the simplest case, it may happen that the mesh is just not fine enough to get the correct solution; thus, the correct expected mode will not appear in favor of another mode close to the expected one. However, in another hand, what might also occur is that basing on a fixed domain size, the code will hesitate between two very close values of the growth rate with two different attached wave numbers; thus, this could lead to the observation of another mode close to the one expected.

• Another interesting thing to notice is the significant error percentage δ kx due to the considerable difference between the analytical and numerical k x . The most probable explanation for this difference is the axial boundary conditions put in MAGNIS are such that the perturbations are set to zero (due to the fact that the boundary conditions imposed for the axial profiles apply for the total variable, meaning the equilibrium plus the fluctuating part), which is not the case in our dispersion relation in which we presumed the instability to be a single plane wave. The axial boundary conditions in MAGNIS however can only be rigorously satisfied by standing wave in x, i.e. a superposition of two waves with the same amplitudes and opposite real k x values:

η(x, y, t) = iη 1 2 exp (ik x+ x + ik y y -iωt) - iη 1 2 exp (ik x-x + ik y y -iωt) (4.114)
where

k x± = ± πn x L x + ik x,i (4.115) such that η(x, y, t) = η 1 exp(-k x,i x) sin(πn x x/L x ) exp (ik y y -iωt) (4.116)
is always zero at x = 0 and x = L x , with n x = 1, 2, 3, .... Moreover, the axial wave number in MAGNIS can have an imaginary part, as is illustrated by the above equations, but which we did not consider in the dispersion relation analysis. Indeed, a close look to the observed oscillations of the azimuthal component of the electric field let appear some signs of an imaginary part for k x ; the instabilities tend to develop more on one side of the domain (in the left for all cases).

• Lastly, one may notice in some cases that either the real frequency is better solved than the growth or the opposite (noticeable in case 1 via the error percentage, for instance). A possible reason for this may also be the boundary conditions in the axial direction; indeed, the axial dependency could be more present in one quantity (real frequency or growth rate) rather than the other, which ends up being more impacted.

Extension of the comparisons for different E and B values

In this last part, we aim to show the good agreement of analytical and numerical results in general for the three cases as a function of different values of E and B. From the observations made in the previous part for the first results, it appears that the comparisons involve in fact two different questions:

1. For a given domain size, does MAGNIS reproduce the exact same mode number that is expected from the dispersion relation? We saw that for this question there can be complications such as two successive modes having almost the same growth rate, axial boundary conditions... A general observation that stands out from these figures for the three cases is how in most of them, the axial boundary conditions was indeed responsible of the issues; this is particularly conspicuous in figures 4.17, 4.19, 4.22 and 4.21, where the analytical variables calculated with the numerical wave numbers are the closest to the numerical ones compared to the values expected by theory. However, the axial boundary conditions do not seem to be solely responsible of the small disparities between the numerical and analytical results, as one may see in figure 4.16 for instance.

Overall, the results show a good agreement between the theoretical and numerical values for all cases. This allows to be confident about the code's ability to properly solve the physical equations of the system, and thus, to validate the fact that these instabilities are a result of this same system, meaning they are indeed physical, and not numerical artifacts.

Conclusion

This chapter showed the linear analysis of a system based on a simple geometry consisting of a cartesian plane perpendicular to the magnetic field with periodic boundary conditions in the azimuthal direction and imposed equilibrium density and potential profiles in the axial direction; we called this system the simple model. We then linearized the same fluid equations present in MAGNIS (but assuming constant uniform magnetic field and electron and ion temperatures) to build a general linear dispersion relation attached to this configuration. First, we showed that the linear unstable modes introduced in chapter 3 and observed in the plasma sources of our interest can be recovered with our dispersion relation, and proceeded to study their behaviour in function of their respective key parameters. Thanks to this theoretical study, we highlighted the possible problematic cases for the code, and isolated the favorable ones.

After some required modifications in MAGNIS to force it in a linear regime and to implement the tools necessary for the diagnostics, we simulated the most favorable cases and compared their attached numerical growth rates, real frequencies and wave numbers with the ones predicted by the theory. We proposed an error assessment and discussed the possible reasons of the slight differences observed between the numerical and theoretical results, but overall, we observed a very good agreement between them. This process allowed us to confirm the physical origin of these instabilities since they clearly are a result of the system of equations and not numerical artifacts.

Confident about the abilities of our code, we will, in the following chapter, focus on more realistic configurations close to real plasma sources.

Chapter 5

Non-linear regime : Analysis of non-linear effects and anomalous transport

Introduction

The previous chapter dealt with the linear analysis of a simple model defined by a Cartesian geometry and specific conditions for the characteristic variables representative of the model. Based on this simple model configuration, we built a dispersion relation with which we recovered the main unstable modes expected in this configuration and listed in Chapter 3. We studied their behaviour according to the key parameters of each mode and highlighted ideal cases along with problematic ones due to a non-stabilizing growth rate at high wave numbers k for the considered instability. We then compared MAGNIS to the ideal cases found with our dispersion relation and were able to confirm that the instabilities observed in the simulation results are a result of the physical equations and not numerical errors.

In this chapter, we aim to understand the evolution of these linear instabilities in a more realistic framework; a self-consistent configuration that we let evolve in a non-linear regime. In comparison with the linear regime which is easily analyzed with a linear theory via dispersion relations based on the assumption of one dominant unstable mode in a weakly perturbed system, the non-linear regime constitutes a much more complex subject to apprehend due to many unpredictable effects induced by the behaviour and coupling between numerous unstable modes with comparable amplitudes, with no existing theory that allows to make an analysis. Moreover, this regime has the particularity to see the different variables of the system strongly deviate from their equilibrium profiles; this non-linear phenomenon is typical of an anomalous diffusion of the plasma across the magnetic field induced by an anomalous transport, a notion we briefly mentioned in the introduction, unpredictable with the classical theory. This chapter will be the opportunity to characterize this non-linear regime and its effects; for this purpose, we still remain in a Cartesian domain with a periodic azimuthal direction in which we perform diagnostics to characterize frequencies, velocities and wave lengths and quantify an anomalous transport, from the simplest cases to the most complex ones. After that, our interest focuses on real plasma sources, the magnetized plasma column and the Hall thruster, which configuration are close to the cases studied beforehand.

Evolution of instabilities in a non-linear regime : qualitative description

In this section we will present MAGNIS simulations of non-linear plasma instabilities for a number of simple test cases. In these cases, we consider a Cartesian domain with periodic boundary conditions in the y direction, as in the previous chapter. Unlike the previous chapter, we now apply realistic sheath boundary conditions in the x direction and let the plasma density profile evolve self-consistently as a result of sources and transport losses, such a to satisfy the equilibrium continuity equations, e.g. for electrons:

∂ n e v e,x y ∂x = S y , (

where the angle brackets indicate averaging over the y direction, and parallel losses are not taken into account (L = 0). Without instabilities, all variables in this equation are independent of y (i.e. we have a pure 1D system) and the axial equilibrium electron flux obeys the drift-diffusion equation:

n e v e,x = -µ ⊥e n e E x + ∂(n e T e ) ∂x . (5.2)
However, in the presence of non-linear instabilities and fluctuations in the y direction, one typically finds that the average axial electron flux is larger than expected from this expression:

n e v e,x y > -µ ⊥e n e E x + ∂(n e T e ) ∂x y , (5.3)
which is known as anomalous transport. It is then common to define an effective mobility µ ef f that takes into account the anomalous transport, such that the classical equation holds on average:

n e v e,x y = -µ ef f n e E x + ∂(n e T e ) ∂x y . (5.4)
As a result of this anomalous transport, the axial equilibrium profiles of n e y (x) and all other variables may change. We will now study these phenomena in MAGNIS.

Definition of test cases

We study different cases, the most basic one to begin with, in which we choose to look at the effects produced by each key parameter of the system independently. Then we proceed to study more complex cases, closer to real plasma sources, where all parameters are coupled and determined self-consistently, so that we can bring out some tendencies that are likely to happen in some real sources. For each case, as said above, we consider a Cartesian plane with periodic azimuthal direction (y) and classical sheath boundary conditions in the axial direction (x), and parallel losses (z) are not taken into account (L = 0). The ion mass and rate coefficients are that of molecular hydrogen. We define the cases as follows:

• Case A : Isothermal case with a uniform magnetic field -In this case, we fix constant values for the magnetic field and the electron temperature. With this fixed electron temperature, it is not possible to describe the ionization source term self-consistently as a function of T e (via the ionization rate coefficient k iz (T e ) defined in a look-up table, see Chapter 2). Therefore we apply the ionization aritficially, assuming that the ionization frequency is uniform with a value

ν iz = I iz /e V n e dV (5.5)
such that the total volume integrated ionization source term (total number of electrons and ions created in the domain per unit time) has a fixed value I iz /e, where I iz is this total source term value expressed like a current in units Ampere.

• Case B : Isothermal case with a magnetic field profile -Unlike Case A, this case presents a Gaussian axial profile for the magnetic field that we fix externally such that:

B(x) = B max exp - (x -x 0,b ) 2 2∆x 2 b . (5.6)
where x 0,b , ∆x b and B max are input values. To make this case similar to a magnetic filter plasma source (see Chapter 2), we also define a profile for the ionization frequency such that the plasma is created only on one side of the magnetic field barrier. As in Case A, the magnitude of the ionization frequency is determined from a fixed value I iz for the total source term:

ν iz (x) = P iz (x) I iz /e V P iz n e dV (5.7)
where P iz is the profile of the ionization frequency, which we again take to be Gaussian:

P iz (x) = exp - (x -x 0,iz ) 2 2∆x 2 iz . (5.8)
with x 0,iz and ∆x iz input values.

• Case C : Non-isothermal case with a constant magnetic field -Here, the magnetic field is constant as in Case A, along with the total injected heating power, but the electron temperature is determined self-consistently by solving the energy equation; this implies that the ionization is also determined self-consistently, since the ionization rate is function of the electron temperature in our model (see Chapter 2).

• Case D : Non-isothermal case with a magnetic field profile -Finally, Case D gathers the electron energy and the magnetic field inhomogeneity effects. The magnetic field profile is, as in Case B, fixed manually, and the electron temperature, along with the ionization, are determined by the energy equation. An injected heating power density profile is set such that (see equation (2.12)):

H(x) = P tot n e P H (x)
V n e P H dV ,

where P tot is the total injected power and P H (x) is again a Gaussian profile, similar to (5.8) above but with input parameters x 0,h and ∆x h . For each of these different cases, the neutral gas density is fixed as n g = 10 20 m -3 and a constant ion temperature is assumed of T i = 0.1 eV. The domain dimensions are 5 cm by 5 cm, and for the isothermal cases, we choose to fix the electron temperature at T e = 6 eV and a the total ionization source at I iz = 40 A. For non-isothermal cases, the total injected power is P tot = 100 W. Moreover, these cases present no applied voltage, and we run simulations with five different values of the magnetic field for every four case, such that B =[25 G, 50 G, 75 G, 100 G, 125 G]. Figures 5.1 and 5.2 display the input parameters and profiles for each case.

Diagnostics

Once have we introduced each different case we chose to study, we proceed to show how we implemented the necessary diagnostics to investigate and characterize the nonlinear modes resulting from these cases. First, we aim to characterize the frequency and velocity of the nonlinear instabilities that we observe, after which we compare the obtained values of these variables with the characteristic plasma scales and velocities in these configurations. Secondly, we quantify the anomalous transport these instabilities generate for each case. The following parts of this section explain how these diagnostics are set up.

Diagnostics for the measurement of frequencies and wave velocities

Since the plasma fluctuations in the non-linear regime contain a wide range of Fourier components and higher harmonics, we found that in this regime it was very difficult to extract useful information from the Fourier transform diagnostic we used in the previous Chapter. We therefore implemented another diagnostic method that provides a good estimate of the dominant wave propagation velocity. This method involves recording the time evolution of one of the model variables, for which we use the plasma density, at successive points in space (along the direction for which we wish to measure the wave velocity) If we consider the plasma density (or another variable) to be of the following form:

n(x, t) = F (k • x -ωt) (5.10)
where F (k • x -ωt) is an arbitrary function, then the derivatives of the density can be written as follows:

∂n ∂t = -ωF (5.11)
∂n ∂y = k y F. (5.12)
A reliable estimate of the wave velocity in the y direction can be found from the above expressions by approximating them by finite differences and then averaging them over time as follows: (5.13) density and potential profiles, and for non-isothermal cases, the electron temperature profile. We propose to analyze the obtained profiles and their tendencies, and show the diagnostics results.

v prop = ω k y = -

Case A

In this case, we remind that all input parameters are constants and thus uniform. We first run this case linearly for the magnetic field value B=25 G; figure 5.3 shows the results for the azimuthal component of the electric field on 128x128, 256x256 and 512x512 mesh grids. One can see the instabilities are developing at the edges of the domain with a total number of nine periods. We let this case evolve non-linearly, along with the other values for the magnetic field. Figure 5.4 shows the density and potential profiles in a 256x256 and 512x512 mesh grid respectively, for the five different values of the magnetic field given above. A first thing to notice is the good mesh convergence for this case, as one may see that the number of periods and the minimum and maximum values are conserved in both mesh for both quantities.

The instabilities appear for all values of the magnetic field; the density profiles display long structures propagating along the azimuthal direction, all disposed around the density peak situated at the center of the domain. One may notice that the more the magnetic field value increases, the more periods appear in the domain. However, if we compare with the linear regime (see chapter 4) for this uniform configuration, the corresponding profiles present more periods in the domain than those in the nonlinear regime; such a behaviour could be explained by an inverse energy cascade phenomenon, in which the smaller scales of motion communicate their energy to the larger scales of the system [START_REF] Kadomtsev | Plasma turbulence[END_REF][START_REF] Sudan | Theory of strongly turbulent two-dimensional convection of lowpressure plasma[END_REF][START_REF] Diamond | Physical Kinetics of Turbulent Plasmas[END_REF] (in opposition to the classical energy cascade, where the large scales transfer their energy to the small ones). As a result, the larger wave lengths tend to dominate the system over the smaller ones in the nonlinear regime. We propose a characterization of the nonlinear instability in table 5.1, which lists the measured frequency, wave velocity and wave length, along with the different characteristic frequencies, velocities and lengths of the plasma dynamics.

Case A B=25 G B=50 G B=75 G B=100 G B=125 G v prop (m.s -1 ) 1.479×10 4 1.367×10 4 1.381×10 4 1.452×10 4 1.430×10 4 λ instab (m) 8.33×10 -3 7.14×10 -3 6.3×10 -3 4.55×10 -3 4.16×10 -3 ω instab (s -1 ) 1.116×10 7 1.203×10 7 1.377×10 7 2.005×10 7 2.160×10 7 v p (m.s -1 ) 1.377×10 4 4.956×10 4 3.794×10 4 3.319 ×10 4 3.022×10 4 v E (m.s -1 ) 1.752×10 4 6.947×10 3 3.642×10 3 2.083×10 3 1.262×10 3 c s (m.s -1 ) 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 ρ e (m) 2.336×10 -3 1.168×10 -3 7.788×10 -4 5.840×10 -4 4.674×10 -4 L n (m) 3.887×10 -2 2.529×10 -2 2.096×10 -2 1.721×10 -2 1.485×10 -2 ω LH (s -1 ) 7.282×10 6 1.456×10 7 2.184×10 7 2.913×10 7 3.641×10 7 ω c,e (s -1 ) 4.397×10 8 8.794×10 8 1.319×10 9 1.759×10 9 2.198×10 9 ω cs (s -1 ) 1.283×10 7 1.497×10 7 1.696×10 7 2.349×10 7 2.569×10 7 Table 5.1 -Case A : Table of characteristic values.
The characteristic plasma quantities in this table are determined only in a very approximative way; first, the ion sound speed in this isothermal case is simply c s = eT e /m i , where we remind that T e = 6 eV, and the ion mass is the hydrogen mass. The diamagnetic and the drift velocities, v p and v E respectively, are absolute averaged values of the ratio between the axial variables (pressure term and electric field) and the magnetic field, defined as such :

v p = - 1 n e B ∂n e T e ∂x
x,y,t

)

v E = - E x B x,y,t . (5.22)
We also determine the absolute averaged gradient length L n such that :

L n = ∂ ln n e ∂x -1
x,y,t .

(

As for the listed frequencies, lower hybrid, electron cyclotron and ion sound respectively, we determined them with their basic expressions (see Chapter 3). However, a detail must be mentioned for the calculation of the ion sound mode ω cs , which is the determination of the attached wave number of this frequency;

we estimated a k y via the estimated wave length, such that :

k y = 2π λ instab . (5.24)
A first observation is that these instabilities seem to propagate at a velocity close to the ion sound speed for the major part of these five cases, in the same direction as v p . As a result, their averaged measured frequency remains equally close to the ion sound frequency. This suggest that case A in the nonlinear regime seems mainly driven by the ion inertia dynamics. Note that this trend is in agreement with the linear dispersion relation, which also predicts a wave velocity close to the ion sound speed (see equation (4.92)). Note furthermore that the instability frequency in table 5.1 is always far below the electron cyclotron frequency, and the characteristic wavelength of the structures is smaller than the gradient length L n and lager than the electron Larmor radius ρ e , always between these two characteristic length.

The figures displayed in 5.5 show two types of axial density and potential profiles:

1. Azimuthally and time averaged 2D profiles (solid lines): profiles n e y,t (x) and φ y,t (x) from the full 2D MAGNIS simulations shown in Figure 5.4, averaged over the y direction and in time 2. Pure 1D stationary profiles (dashed lines): profiles n e (x) and φ(x) obtained from purely 1D simulations, for which we adapted MAGNIS to solve only in the axial direction By comparing these two profiles for each case, one can get an idea about the anomalous transport and how it affects the plasma equilbrium. Considering this case is ambipolar and uniform, the axial equilibrium equation (5.1), once coupled with the ion transport, takes the form of a simple diffusion equation:

-D amb ∂ 2 n e ∂x 2 y = S y = ν iz n e y , (5.25)
where

D amb = µ i µ ef f µ i + µ ef f (T e + T i ) ≈ µ ef f T e (5.26)
is the ambipolar diffusion coefficient, which in magnetized plasma case is approximately the electron diffusion coefficient. The solution of this equation is given by

n e y ≈ N 0 sin (πx/L x ) , (5.27)
and thus one obtains for the average plasma density: The tables below (tab.5.2 and 5.3) display the average values of these velocity profiles for both region and each value of B max , along with other characteristic plasma quantities. Since this case has the specificity to work with a magnetic field profile, some regions of the domain, particularly the edges and their vicinity, have very low magnetic field values (near zero); which leads to undefined v p and v E velocities for these regions and thus unreliable intel for comparisons. Instead, we propose the total electron velocity v e,y , sum of these velocities along with the collisional contribution as a reference velocity of the system. The quantities dependent on the magnetic field B were calculated with an averaged value of B in both regions. One may see that the measured velocities, along with the attached calculated frequencies, for each magnetic field value, are very close to the ion sound speed and frequency, as in Case A, in both regions; which again, illustrates that the ion dynamics dominates the system. As for the propagation direction, these instabilities all propagate in the same direction with v e,y . Note again that this is in agreement with the linear dispersion analysis of equation (4.92).

n e x,y ≈ L 2 x π 2 S x,y D amb ≈ L 2 x π 2 S x,y µ ef f T e . ( 5
Case B B max =25 G B max =50 G B max =75 G B max =100 G B max =125 G Region 1 v prop (m.s -1 ) 1.422×10 4 1.507×10 4 1.502×10 4 1.540×10 4 1.425×10 4 λ instab (m) 2.5×10 -2 1.67×10 -2 1.67×10 -2 1.67×10 -2 1.25×10 -2 ω instab (s -1 ) 3.574×10 6 5.670×10 6 5.651×10 6 5.794×10 6 7.163×10 7 v e,y (m.s -1 ) 1.285×10 4 1.710×10 4 1.828×10 4 2.160×10 4 2.391×10 4 c s (m.s -1 ) 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 ρ e (m) 4.5×10 -3 1.95×10 -3 1.43×10 -3 1.06 ×10 -3 8.02×10 -4 L n (m) 5.71×10 -2 8.33×10 -2 1×10 -1 1.1×10 -1 1.2×10 -1 ω LH (s -1 ) 3.78×10 6 8.73×10 6 1.19×10 7 1.60×10 7 2.12×10 7 ω c,e (s -

Case B

B max =25 G B max =50 G B max =75 G B max =100 G B max =125 G Region 2 v prop (m.s -1 ) 1.495×10 4 1.631×10 4 1.765×10 4 1.810×10 4 1.725×10 4 λ instab (m) 6.25×10 -3 4.17×10 -3 3.13×10 -3 2.77×10 -3 2.5×10 -3 ω instab (s -1 ) 1.503×10 7 2.458×10 7 3.543×10 7 4.106×10 7 4.335×10 7 v e,y (m.s -1 ) 4.554×10 4 4.157×10 4 4.361×10 4 4.596×10 4 4.623×10 3 c s (m.s -1 ) 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 1.701×10 4 ρ e (m) 4.5×10 -3 1.95×10 -3 1.43×10 -3 1.06 ×10 -3 8.02×10 -4 L n (m) 3.57×10 -2 3.33×10 -2 2.86×10 -2 2.22×10 -2 1.89×10 -2 ω LH (s -1 ) 3.78×10 6 8.73×10 6 1.19×10 7 1.60×10 7 2.12×10 7 ω c,e (s -
A quick look at fig. 5.10, representation of the stationary 1D and averaged 1D axial density and potential profiles for the different values of the magnetic field, enables to see the sign of an anomalous transport; indeed the time-averaged profiles differ considerably from the stationary ones. As we did previously in case A, we quantify an anomalous transport by measuring an effective mobility in the given results; fig. 5.11 displays the classical and anomalous mobilities for all magnetic fields obtained for two different mesh numbers. As one may notice, there is an important difference between the effective mobilities and the classical ones, which confirms the presence of anomalous transport in this case too. Interestingly, the anomalous transport is much higher on the right side of the barrier than on left side of the barrier, possibly related to the absence of ionization on the right side. The effective mobilities appear to follow the classical theory until the magnetic field reaches high values, after which they become higher than the mobilities predicted by the classical theory until the other end of the domain where they finally follow again the classical mobilities, since the magnetic field is getting weaker. Moreover, the mesh convergence issues encountered above appear also in the diagnostics of the measured mobilities to a small extent.

Case C

In this case, we remind that the conditions are almost like Case A, except we activated the energy equation meaning that the creation, along with the electron temperature, are self consistent. Figures 5.12 and 5.13 show the density, potential and electron temperature profiles respectively attached to this case for two different meshes. If we compare these results with the ones obtained for an isothermal case, we notice that the same tendency is recovered for the density and potential profiles. However, some differences can be pointed out, such as the lower values for this case, due to the different input conditions; indeed, we fixed high values for the electron temperature and the ionization rate for the isothermal case, which as a consequence, leads to higher values for the displayed variables. Another difference here is the slight discrepancy between the results on the two different meshes probably due to a mesh convergence issue, especially for the highest values of the magnetic field.

Tab.5.4 presents the characterization of the instabilities properties along with the other expected variables of the system; in this case where the electron temperature is self-consistently determined, the temperature-dependent quantities in the table below were calculated with an averaged value of T e . A look at table 5.4 enables for one to see that his case does not differ from Case A, since here also, the instabilities properties (velocities and frequencies) are mainly driven by the ion dynamics and the propagation direction is again the same than v p . We proceed to quantify this anomalous transport just like the two previous cases, that we present also under the form of an effective mobility in Figure (5.15); the classical and effective mobility profiles in this case also show no difference compared to the ones given in Case A; there again, the same trend is recovered. And as we did for Case A, we also show a comparison between the average effective mobility and the one predicted by Bohm, in Figure 5.16; with no surprise, one can observe the same trend, with however a little gap between the effective mobility and Bohm values.

Case

C B=25 G B=50 G B=75 G B=100 G B=125 G v prop (m.s -

Case D

In this part, we aim to investigate the last case described in the introduction of this chapter, which is defined with a magnetic barrier similar to case B, but where the ionization and the electron temperature are coupled and determined self-consistently. In figures 5.18 and 5.17, we show the simulation results for this case, but this time only for the first three different values of B max ; indeed, a close look to the results enables to realize this case is problematic, as significant mesh convergence issues already appear for B max = 50G, which only get worse as we increase the magnetic field values. As one may see in figures 5.18 and 5.17, two distinct regions are forming, where the characteristic wavelength of the observed instabilities changes, as in Case B; in analogy with this case, we also propose a characterization of the instabilities for both regions in tables 5.5 since the diagnostics results are similar on both mesh grids. In this case, the characteristic values dependent on T e were calculated with an averaged value of the latter in both regions. This table underlines that just like the three previous cases, every characteristic property of the system follows the ion dynamics, since, again, the values for the measured velocities, along with the resulting frequencies, are comparable to the ion sound speed and frequency.

Case D

B max =25 G B max =50 G B max =75 G Region 1 v prop (m.s -
Although we were able to characterize the main properties for this case, it was not possible to determine the equilibrium profiles and an anomalous transport; indeed, it appears that the 1D solutions are unstable and oscillate strongly. In other words, our 1D MAGNIS simulations could not find any stationary equilbrium solutions for this case, so perhaps they do not exist. In contrast, the 2D simulations in figures 5.18 and 5.17 do achieve an axial equilibrium state. This suggests that anomalous transport plays an important role here: the axial equilibrium can be achieved only thanks to anomalous transport. However, when we tried to measure this anomalous transport, it turned out that the diagnostic tools we implemented were unable to measure it in a meaningful way, e.g. the effective mobility profile showed strange oscillations and peaks and was negative in some places. This was the case both with our diagnostic based on equation (5.16) and with that based on (5.18) via the Hall parameter.

A possible explanation for this failure could be that the characterization of an anomalous transport via an effective mobility cannot be systematically applied; indeed, the concept of effective mobility implies that the axial electron flux is proportional to the local axial electric field and pressure gradient, but it could be that such proportionality is not garanteed for anomalous transport. For example, the net force in the denominator of equation (5.16) could be zero but not the flux in the numerator. We did not further investigate this explanation. Yet, it is remarkable that this effective mobility problem occurs for precisely this case for which there also appear to be no stationary equilibrium solutions.

The whole study of these four different systems enabled to highlight the main trends that are expected for each specific case. We now aim to simulate two systems with more realistic conditions, namely the magnetized plasma column and the Hall thruster.

Magnetized plasma column

In this section we present MAGNIS simulation results of a configuration that resembles our Case C above but is closer to that of certain real plasma sources of interest for a novel neutral beam injection system concept for fusion [START_REF] Simonin | R & d around a photoneutralizer-based nbi system (siphore) in view of a demo tokamak steady state fusion reactor[END_REF][START_REF] Hagelaar | Modeling of plasma transport and negative ion extraction in a magnetized radio-frequency plasma source[END_REF][START_REF] Jacquier | First b-dot measurements in the raid device, an alternative negative ion source for demo neutral beams[END_REF][START_REF] Furno | A novel helicon plasma source for negative ion beams for fusion[END_REF][START_REF] Furno | Helicon wave-generated plasmas for negative ion beams for fusion[END_REF][START_REF] Agnello | Negative ion characterization in a helicon plasma source for fusion neutral beams by cavity ring-down spectroscopy and langmuir probe laser photodetachment[END_REF][START_REF] Simonin | Negative ion source development for a photoneutralization based neutral beam system for future fusion reactors[END_REF] (CYBELE or RAID devices in figures 3 and 7 in the Introduction of this thesis): a magnetized plasma column. These sources have a very long rectangular chamber with a uniform axial magnetic field (all along their length) in which the plasma is created and sustained by RF heating. Here we consider an idealized configuration in which the RF heating is uniform and all effects of parallel transport along the column are completely neglected. The plasma is thus described in a square 2D domain perpendicular to the magnetic field, with grounded walls on all boundaries, as shown in Figure 5.19. This geometry is very similar to Case C above, except that the periodic boundary conditions on the top and bottom boundary (y direction) are now replaced by the classical sheath boundary conditions. As a consequence, the ∇p × B that is closed in the y direction in our Case C, is now closed in a loop around the center. So this configuration still has closed drift, even though it is not explicitly periodic. We simulated this configuration for a large range of magnetic field values from 0 to 400 Gauss and on different numerical meshes. Other input parameters we used in the simulations are given in table 5.6. We found that when increasing the magnetic field, there appear to be different regimes where the instabilities and plasma structure look different. This is illustrated by Figure 5.20, which shows typical profiles of the main plasma variables (plasma density, electron and ion flux vectors, potential and electron temperature) for each of these regimes:

1. At low magnetic field (B =15 G case in first column of Fig. 5.20) the plasma is totally stable and has a classical diffusive profile.

2. From B =20 G, instabilities appear in the form of thin structures stretching out from the center to the walls and rotating around the center (Fig. 5.20, second column). These structures get thinner and thinner as we increase the magnetic field, and until B =150 G they look very similar to what we get for Case C.

3. Beyond that, these structures start to gradually disappear and a more stable star-like plasma configuration appears, typically with four arms (Fig. 5.20, third column).

4. Finally, for strong magnetic field above B =200 G (Fig. 5.20, last column), another kind of instability occurs with larger and different looking structures, probably related to the fact the that the ions become magnetized (this regime is not there if we remove the magnetic force from the ion momentum equation in MAGNIS).

We then characterized these instabilities by different diagnostics. However, since there is no periodic direction in these simulations, this characterization is more difficult than for the periodic cases above. We adapted our diagnostic for the wave propagation speed so that it measures the both the x and y components of the wave velocity. From this, we found that once again, the plasma structures propagate with a speed close to the ion sound speed and in the direction of the ∇p × B drift, for all instabilitiy regimes. This is illustrated in Figure 5.21, which shows the absolute wave velocity and the ion sound speed as a function of magnetic field, both for the plasma column and for the periodic Case C above. Determining an effective mobility for these plasma column simulations turned out to be very difficult, as the methods used above rely on averaging over the direction of the magnetic drift, which is simply the y direction in the periodic cases, but which is not so easy to identify here. We therefore tried to characterize the anomalous transport in a different way, by means of a global plasma transport time (or confinement time), defined as:

τ = V ndV t V SdV t = n x,y,t S x,y,t , (5.30)
the ratio of the average plasma density over the average ionization source, or in other words, the total number of electrons present in the domain over the total number of electrons created per unit time. This is shown in Figure 5.22 as a function of magnetic field and for different meshes. One can see that the global transport time increases when the magnetic field increases, as expected, but that the rate of this increase changes in different regimes. There are also some significant differences between the meshes, in particular for higher magnetic fields, although the 256 and 512 points mesh simulations are in good agreement up to B = 200 G.

In order to get an idea about the anomalous transport, we have to compare these global transport times with those corresponding to classical stationary transport. For this, we used stationary solutions from pure 1D simulations with a domain length of 7.07 cm, representing a 1D approximation of the 2D domain of 10 cm by 10 cm for the case of diffusive transport. (Note that stationary classical solutions cannot be obtained directly for the 2D domain.) The global transport times for classical transport increase approximately as B 2 (because µ -1 ⊥e ∝ B 2 ). In view of our previous findings that the effective mobility was close to the Bohm mobility, we also calculated the global transport times for 1D solutions including the Bohm mobility via an anomalous frequency such that [START_REF] Adam | Physics, simulation, and diagnostics of hall effect thrusters[END_REF]: The figure confirms the presence of an anomalous diffusion, given the noticeable difference between the 1D and 2D results in particular for intermediate magnetic field strength. It is interesting to make the link between these transport times and the different instability regimes observed on the profiles of Figure 5.20 above. For low values of the magnetic field the transport time follows the classical trend. Then, in the second regime with strong thin instability structures, it becomes very close to the Bohm transport values, but after that, in the third regime, the transport time rapidly increases to get back to the classical values, until the point where the larger scale structures appear and it stops increasing. So, strong anomalous transport is linked with strong instabilities visible on the profiles. The trend at high magnetic field is not completely clear due to mesh convergence issues.

ν ef f = ν Bohm + ν m = ω c,e 16
The right hand panel of the Figure 5.22 shows the global transport times for Case C, corresponding to the mobilities shown earlier for this case. These curves are very similar to the first part of those for the column simulations, but the magnetic field range is much smaller so that the high-field regimes do not appear. This comparison along with the wave velocity measured for both cases, confirms Case C and the present plasma column simulations despite the different boundary conditions.

Finally, we remark that the idealized plasma column configuration is currently being considered as a benchmark case between MAGNIS and PIC simulations developed at the LAPLACE laboratory. Figures 5.23 and 5.24 show some preliminary results of this benchmark that seem to confirm the capabilities of MAGNIS to correctly describe this plasma configuration.

Magnetic barrier : Hall thruster alike case

Another representative source for the periodic cases studied previously, specifically with a non-uniform magnetic field (cases B and D), can be the Hall thruster, which shares some similarities with case D in the configuration. In this section, we try to reproduce the Hall thruster conditions in MAGNIS and make simulations of this configuration. We first define the case and then show the results and discuss them. In the case of the Hall thruster, we focus on the channel in which we choose as a simulation plane for MAGNIS the axial direction and a portion of the azimuthal direction as displayed in fig. 5.25, perpendicular to the magnetic field. As for the input parameters, the Hall thruster configuration is very similar to case D, but here we choose to include the gas density dynamics and to apply a strong electric field along with a cathode current.

For this, some modifications in the code were required; first, to take into account the neutrals dynamics in the system. A way to proceed is to either fix a profile, or, like we chose to do, to model it thanks to a continuity equation such that :

∂n g ∂t + ∇ • (n g v g ) = -S (5.32)
where v g is a fixed gas flow velocity of the order of the Maxwellian thermal speed :

v g ≈ 1 2
8eT g πm g x (5.33)

T g and m g being respectively the neutrals temperature and mass. After that, we remove the basic axial sheath boundary conditions at the end of the channel and replace it with a cathode current, while the azimuthal boundary conditions remain unchanged. We show the input profiles fig. We run the simulations on different meshes for this case in an attempt to see how it behaves when every variable in the system is coupled to one another and solved self-consistently. We present the results in the following subsection, fig.5.27. Figure 5.27 displays the results on four different mesh for the density, potential and electron temperature profiles. As one may see, the profiles seem to be correctly solved on a coarse mesh (64x32), yet, if we progressively refine it, some serious mesh convergence issues appear, already for a mesh little finer than the first one, in a similar way than met in case D. However, this Hall thruster case demonstrates more dramatic issues than said case D, since the numerical results are diverged in the finest mesh (256x128).

These issues are most likely due to some limitations in our model which provoke the appearance of small scales that get always smaller than the mesh cell's size and cause the simulation to crash for the particular case of the Hall thruster. Another attempt to model this configuration is to make simulations with uncoupled variables; this time, we consider an isothermal case where the electron temperature is given as an input and the ionization is determined via a profile, just like in case B, which was a successful case. As for the input values, the electron temperature is fixed at T e = 3 eV, and the ionization rate at I iz = 0.5 A. The simulation results for this isothermal case is shown fig. 5.28. Unfortunately, one can clearly observe that the same issue is met even for this simplified case, in an even worse way; indeed, if we compare with the previous self-consistent case, the numerical results diverged abruptly for a less fine mesh. While the first attempt showed signs of mesh convergence issues in the first three grids before numerically crashing in the last mesh grid, this isothermal case didn't display any signs of possible dramatic issues for the first mesh grids but show sudden serious issues for the last one. A general observation for the Hall thruster case is the corrupted simulations resulting from small scales problems. As explained above, these issues underlines the weakness of our code and fluid models in general, to model this kind of plasma sources, eventually due to numerous approximations made to achieve a macroscopic description when the multi-scale nature of the system plays an important role.

Conclusion

Throughout this chapter, we studied and characterized the non-linear regime through basic and simplified cases in which we progressively add complexities (self-consistent ionization and electron temperature, magnetic field inhomogeneity...), so that we can understand each mechanism effects separately, and this, for different values of the magnetic field. We implemented the necessary tools to perform the diagnostics in order to characterize the main properties (wavelength, propagation direction, and wave frequency and velocity) of each case.

The whole study of these four different systems enabled to highlight the main trends that are expected for each specific case. So far, we saw that for all four cases, the nonlinear instabilities propagate with a velocity close to the ion sound speed and in the direction of the magnetic electron drift, as is also predicted by linear theory. Also, we saw that MAGNIS seems to perform best for cases with a uniform magnetic field, regardless of the presence of the energy equation. However, the combination of an inhomogeneous magnetic field with self-consistent ionization and electron temperature seems to be more challenging as some issues are arising visible via mesh convergence problems due to the instabilities small structures.

Moreover, we were able to quantify an anomalous transport via an effective mobility in the first three cases, and showed for the magnetically uniform cases that they follow Bohm's predictions. As for the last case, we encountered difficulties in determining a meaningful effective mobility, calling into question the concept of anomalous mobility for this case.

After that, we considered more realistic conditions as we chose to model real plasma sources close to the previous studied cases, namely the magnetized plasma column and the Hall thruster. Concerning the magnetized plasma column, we showed that the MAGNIS results very are similar to those for the periodic uniform case and also in good agreement with preliminary PIC simulations, and that MAGNIS predicts different instability and transport regimes as a function of magnetic field. However, we showed that the Hall thruster case, close to the problematic last Case D, gives really dramatic results as it causes simulation breakdowns when we progressively refine the mesh. We attributed this issue to the limitations in fluid models for this specific case, and discuss further in Chapter 6.

Chapter 6

Exploratory developments and extensions of fluid models

Discussion on approximations in fluid models

As expected and already explained in chapter 1, fluid models present some limitations. In previous chapters, we highlighted the issues met in our own fluid model, namely high-k issues in the linear regime, and small scales effects causing a mesh convergence problem in the nonlinear regime due to the non saturation of the instabilities for a lot of realistic plasma sources. Several reasons can be mentioned, all of them linked to the approximations made on fluid models (see chapter 1), which questions their reliability for their capability of modeling some of our plasma sources of interest. In chapter 1, we have seen that the most common fluid models applied to these sources are :

• drift-diffusion models where electrons and ions are considered massless so that their inertial terms are neglected with respect to the collisional terms. This model shows no particular issues in terms of stabilization or saturation, but is rigorously justified only for the strongly collisional limit, an assumption that may not be applicable to many cases of plasma sources. Some phenomena, such as ion sound waves or the plasma presheath are neglected in this approximation.

• Electron drift-diffusion models, similar to a classical drift-diffusion model except the ion inertial terms are taken into account. These models can present issues since the ion inertia dynamics, when coupled with the electron pressure by via quasi-neutrality assumption, generate unstable ion sound modes that never reach saturation. Usually, some diffusive terms, meant to mimic kinetic effects, are implemented in the equations to overcome this problem (we detail it in the following sections). Another solution is to complete the model with electron inertia, or to couple the system with the Poisson equation.

• Low-frequency approximation models, more complete than the electron drift-diffusion model since electron inertia is added. In this configuration, the issues met in the electron drift-diffusion model are no longer present; however, an assumption on the electron velocity is made so that only the contributions where ω ω c,e are conserved, which, as a result, cuts high frequencies that are equivalent or larger than the electron cyclotron frequency. To sum up, this model is a way to include every contribution in the system of fluid equations, but in which we deliberately cut high frequencies to focus only on the low ones, where the limit ω ω c,e , and thus electron cyclotron resonance effects, are not treated.

In all of these approximated models, one may make an additional assumption based on the quasi-neutral state of the plasma, which is simply the quasi-neutral assumption n e = n i , instead of solving Poisson equation, and thus, neglecting some possible non-neutrality effects. Finally, we could end the list with our own model MAGNIS, which gathers all contributions in the system of fluid equations, but also contains some approximations.

In our model, the main approximations concern the quasi-neutrality assumption and the assumed Maxwellian distribution with the resulting chosen closures for the electron and ion pressure tensors and the electron heat flux. The quasi-neutrality assumption may not take into account some stabilizing mechanisms induced by non-neutrality effects, while the Maxwellian distribution attached to the chosen closures for some of the variables of the system underlines the necessity to neglect kinetic effects due to microscopic behaviour and non-Maxwellian distribution to achieve a macroscopic description. These approximations result, in some cases of partially magnetized plasma sources (including the cases seen in chapter 5), in a failure of the numerical simulations due to skipped physical effects from kinetic contributions.

However, it is possible to overcome these issues and capture the real physics of a particular plasma source if one is aware of the present limitations in the model. Most of the time, these problems are fixed by completing the physics of the model or studying the possible stabilizing mechanisms, most of the time under the form of diffusive terms that will smooth the profiles and disperse the small scales. In the following sections, we propose three improvements to fix these issues; two diffusive terms, each of them implemented either in the ions dynamics (ion viscosity) or in the electron dynamics (effective potential), and the coupling of our system with Poisson equation. We propose a linear analysis of these three possibilities in high-k issues context, and see how they behave in MAGNIS once in the non-linear regime.

However, one must have in mind that these extensions described in this chapter may not solve all problems at high-k values, in particular, in nonlinear regimes, and, at this stage, should be viewed as some exploratory and preliminary developments.

Modeling of kinetic ion effects via an ion viscosity 6.2.1 Definition

In the linear analysis we made chapter 4, we highlighted the issues encountered with the infinite growth rates in (k x , k y ) plane of some unstable modes, causing small scales problems in simulations. In this subsection, we aim to investigate effects of a diffusive term added on the momentum equation for ions that may be seen as an ion viscosity, in analogy with Navier-Stokes equations in fluid mechanics. Just like in the Navier-Stokes equations, this term is mathematically written as a diffusivity multiplying the Laplacian of the ion velocity :

η diff ∇ 2 v i (6.1)
where η diff , the diffusivity, is usually expressed as the ratio between the ion thermal speed and their collision frequency :

η diff = v 2 th,i ν i,n . (6.2)
The physical origin of this term issues from the ion viscosity stress tensor, assuming a present anisotropy in the total pressure tensor for ions due to shear forces induced by thermal action from ions, and thus the presence of non-diagonal elements in the total pressure tensor, such that :

∇ • P i = T i ∇n + ∇ • π i (6.3)
with π i being the ion stress tensor, containing the non-diagonal part of the total pressure tensor and written as follows [START_REF] Smolyakov | Multiscale equatorial electrojet turbulence: Baseline 2-d model[END_REF]:

π i,xy = - nT i ν i,n ∇v i + ∇ t v i - 2 3 δ xy ∇ • v i . (6.4)
Thorough calculation of the components in (6.4) leads to the ion viscosity expression in (6.1).

Dispersion relation and linear analysis

Dispersion relation

The system and the attached equations remain unchanged as in Chapter 4, except that we take into account (6.1) in the momentum equation for ions, which leads to :

∂v i ∂t + v i • ∇v i + ν i v i = e m i (E + v i × B) - eT m ∇n n -η diff ∇ 2 v i (6.5)
The initial general dispersion remains as seen in (4.39), except the ion viscosity has the effect to add a new term in the effective frequency ν i,eff we defined (4.36) such that :

ν i,eff = iω -v i,0 • ik + ν i + η diff k 2 (6.6)
in which η diff is an input parameter of our system.

Linear effects of the ion viscosity on Farley-Bunemann and gradient-drift instabilities

We aim to investigate the effects of the ion viscosity defined above for the problematic cases 2 and 3 described in chapter 4. We remind that the case 2, commonly called Farley-Buneman instability in the literature, never stabilizes at high-k while case 3, that we called the gradient-drift instability, stabilizes when the electron inertia is added to the dynamics of this unstable mode. First, we show figures 6.1 and 6. Looking now at the effects of this diffusive term on the gradient-drift instability as introduced in chapter 4, eq.( 4.78), we present figures 6.3 and 6.4, in which fig. 6.3 is the pure case without any diffusive term, fig. 6.4 the gradient-drift instability with the ion viscosity. Here also, the ion viscosity works very well to stabilize the growth rate. Note that, again in this case, the dynamics isn't conserved since the values change when ion viscosity is added.

The ion viscosity showed its ability for the stabilization at high k. In addition to their effectiveness, this diffusive term has the benefit to not increase the order of the polynomial dispersion relation, which enables to investigate their effects analytically. After calculation, we show that at high K, the growth rate for the Farley-Bunemann and the gradient-drift instabilities with ion viscosity is :

Ω i (K → ∞) = - c 2 s η diff ν i , (6.7)
where the growth rates become negative; both instabilities are damped equally at high K when the ion viscosity is added.

Effects of the ion viscosity in non-linear regime : implementation in MAGNIS

In the beginning of the chapter, we also introduced an issue other than the high-k problems met in our model in a nonlinear regime, being the problem of the instabilities saturation. Usually, this issue manifests via a mesh convergence problem or sometimes even via simulation breakdowns; we met typically these specific cases in the previous chapter, with the case D and the Hall thruster examples. The proposed stabilizing term in this section constitute a possible solution [START_REF] Chaturvedi | Non-linear stabilization of the exb gradient drift instability in ionospheric plasma clouds[END_REF] for these non linear issues. However, if we are sure that this mechanism works linearly, it is hard to predict if it will work successfully for non linear cases, since a non linear regime implies the appearance of new unstable modes coupled to each other, leading to unpredictable phenomena.

The ion viscosity in MAGNIS

As already mentioned in the linear section, the ion viscosity consists in the addition of a Laplacian term in the ion momentum equation (see eq.(6.5)). The diffusive coefficient is determined as the ratio between the ion sound speed and a diffusion length that we give as an input parameter, such that :

η diff = c s L diff , (6.8)
L diff being the diffusion length.

Results

In chapter 5, we came across some unsuccessful cases, particularly the Case D, a fully self-consistent configuration with an inhomogeneous magnetic field. We proceed to show some results for this case in which was taken into account the ion viscosity, where we chose a diffusion length of L diff = 1.2 × 10 -4 m, for Case D with a B max = 100G. Until now, we assumed the quasi-neutrality condition n i = n e = n, which simplifies the model since we avoid solving the Debye length by coupling our system of equations with the Poisson equation, and thus enables to lighten the formalism. However, the quasi-neutrality assumption implies that neutrality is also assured locally, in every point of the system; even if most of the time, this simple condition is a correct assumption and allows to obtain good simulation results, it is physically not really relevant for our considered plasma sources, and in some cases it may happen that this approximation is neglecting local non-neutrality effects that are important in the stabilization mechanisms of the concerned instabilities. In this case, there is little to no other choice left for one to solve this issue than to add the Poisson equation to the system, that we remind below :

∇ 2 φ = - ρ 0 (6.9)
where : ρ = e(n i -n e ) (6.10)

is the volume charge density defined by the difference between the ion and electron densities.

Dispersion relation and linear analysis

Dispersion relation

This sections shows the calculation of a dispersion relation taking into account non-neutrality effects by adding Poisson equation. We remind the configuration of reference is still the simple model geometry (fig. 4.1). Temperatures and the magnetic field are still considered constants and continuity and momentum equations for both species are linearized the same way as we did to obtain (4.39), with the same perturbed forms and the same stationary solutions; indeed, quasi-neutrality is conserved in the stationary regime : n i,0 (x) = n e,0 (x) = n 0 (x) (6.11) with n 0 (x) as defined in (4.12); as a result, the stationary potential is also the same as defined in (4.13).

As for the perturbed part of the density, we consider ñi = ñe and define each as follows :

ñe = ñ - 1 2 n ρ (6.12) ñi = ñ + 1 2 n ρ (6.13) in which : ñ = 1 2 (ñ i + ñe ) = n 0 (x) z0 exp(-iωt + ik • x) (6.14)
can be seen as an average perturbed density, and :

n ρ = ñi -ñe = n 0 (x)z 0 exp(-iωt + ik • x) (6.15)
a general perturbed term being the difference between the two perturbed densities. Linearizing equation (6.9) leads to:

∇ 2 φ = - e (ñ i -ñe ), (6.16)
which, once we inject expressions (6.12) and (6.13) in it, becomes :

∇ 2 φ = - e n ρ . (6.17)
For the very same reasons we explained in chapter 4 where we assumed different source terms for both species to get constant velocities and electric field along with a constant amplitude for the perturbations, in order to include a density profile and thus, gradient effects, we assumed a permittivity profile under the following form : = 0 exp(g • x) (6.18) so that no approximations in our equations are needed to keep our stationary solutions and perturbed quantities compatible with the system. We now proceed to linearize the system of equations. Taking into account the required modifications detailed above to add Poisson equation, we express first the continuity equations for both species as follows :

∂ ∂t (ñ - 1 2 n ρ ) + ∇ • (n 0 ṽe ) + ∇ • ((ñ - 1 2 n ρ )v e,0 ) = Se (6.19)
for electrons, and for ions, we get :

∂ ∂t (ñ + 1 2 n ρ ) + ∇ • (n 0 ṽi ) + ∇ • ((ñ + 1 2 n ρ )v i,0 ) = Si . (6.20)
As for the momentum equations, the linearization gives the following expressions for electrons and ions respectively :

α 1 ∂ṽ e ∂t + α 2 v e,0 ∇ • ṽe + ν m,e ṽe = - e m e (-∇ φ + ṽe × B) - eT e m e ∇(ñ -1 2 n ρ ) n 0 (x) - (ñ -1 2 n ρ ) n 0 (x) ∇n 0 (x) n 0 (x) (6.21) β 1 ∂ṽ i ∂t + β 2 v i,0 ∇ • ṽi + ν m,i ṽi = e m i (-∇ φ + β 3 ṽi × B) - eT i m i ∇(ñ -1 2 n ρ ) n 0 (x) - (ñ -1 2 n ρ ) n 0 (x) ∇n 0 (x) n 0 (x) . (6.22)
Basing on (6.17) in which we inject the Fourier form given in (4.30) as in chapter 4, it is easy to obtain an expression for n ρ such that: ), then, one may obtain the following set of equations, respectively electron and ion continuity equations followed by their momentum equations:

n ρ = - e ∇ 2 φ = e k 2 φ 1 exp(-iωt + ik • x). ( 6
-(z 0 n 0 - 1 2 0 e k 2 φ 1 )iω + n 0 v e,1 • (g + ik) + v e,0 • ik(z 0 n 0 - 1 2 0 e k 2 φ 1 ) = 0 (6.24) -(z 0 n 0 + 1 2 0 e k 2 φ 1 )iω + n 0 v i,1 • (g + ik) + v i,0 • ik(z 0 n 0 + 1 2 0 e k 2 φ 1 ) = 0 (6.25) (-α 1 iω + α 2 v e,0 • ik + ν m,e )v e,1 = - e m e (-ikφ 1 + v e,1 × B) - eT e m e ik(z 0 - 1 2 0 e k 2 φ 1 ) (6.26) (-β 1 iω + β 2 v i,0 • ik + ν m,i )v i,1 = e m i (-ikφ 1 + β 3 v i,1 × B) - eT i m i ik(z 0 + 1 2 0 e k 2 φ 1 ). (6.27)
Determining the expressions of the fluctuating velocity components for both species with (6.26) and (6.27), injecting them in their respective continuity equations (6.24) and (6.25), and finally combining them, we obtain the following dispersion relation in which are taken into account non-neutrality effects :

-iω + v e,0 • ik -eTe me 1 ν 2 e,eff +ω 2 c,e (-ω c,e ik × b • g -ν e,eff (-ig • k + k 2 )) - k 2 λ 2 D 2 (-iω + v e,0 • ik) + eTe me 1+ k 2 λ 2 D 2 ν 2 e,eff +ω 2 c,e (-ω c,e ik × b • g -ν e,eff (-ig • k + k 2 )) (6.28) = -iω + v i,0 • ik + Ti Te c 2 s 1 ν 2 i,eff +ω 2 c,i (-β 3 ω c,i ik × b • g + ν i,eff (-ig • k + k 2 )) k 2 λ 2 D 2 (-iω + v i,0 • ik) + c 2 s 1- T i Te k 2 λ 2 D 2 ν 2 i,eff +ω 2 c,i (-β 3 ω c,i ik × b • g + ν i,eff (-ig • k + k 2 ))
with λ D , the Debye length and ν e,eff and ν i,eff as defined in (4.36).

Study of a Bunemann instability and transition to the linear electron cyclotron drift instability (ECDI) via non-neutrality effects

The electron cyclotron drift instability is a kinetic unstable mode resulting from the electron cyclotron wave growing unstable due to the coupling of the electron dynamics (mainly their inertial terms) with the ion acoustic mode [START_REF] Lampe | Anomalous resistance due to crossfield electron-ion streaming instabilities[END_REF][START_REF] Lampe | Nonlinear development of beam-cyclotron instability[END_REF]. This instability's conditions to trigger presents mechanisms that can be found in some of our plasma sources of interest, the Hall thruster and the magnetron precisely [START_REF] Adam | Study of stationary plasma thrusters using two-dimensional fully kinetic simulations[END_REF][START_REF] Adam | Physics, simulation, and diagnostics of hall effect thrusters[END_REF][START_REF] Tsikata | Three-dimensional structure of electron density fluctuations in the hall thruster plasma: ExB mode[END_REF][START_REF] Tsikata | Modulated electron cyclotron drift instability in a high-power pulsed magnetron discharge[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in hall effect thrusters. ii. kinetic model[END_REF].

Even if this unstable mode is perceived as kinetic (some describe this instability as the result of a coupling between kinetic modes called the Bernstein mode and the ion acoustic mode [START_REF] Gary | Instabilities in perpendicular collisionless shock waves[END_REF][START_REF] Gary | Longitudinal waves in a perpendicular collisionless plasma shock . 1. cold ions[END_REF]) ) and its study remained essentially kinetic [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the e x b electron drift instability: Modeling and fitting on experimental data[END_REF][START_REF] Janhunen | Nonlinear structures and anomalous transport in partially magnetized exb plasmas[END_REF][START_REF] Janhunen | Evolution of the electron cyclotron drift instability in two-dimensions[END_REF], it is actually possible to recover some of this instability characteristic features, namely the effect of the electron cyclotron resonances, from fluid equations coupled to the Poisson equation.

Basing on the dispersion relation in (6.28), we first consider a simplified quasi-neutral case in which we remove the magnetization for ions, the density gradient, the collisions and the temperatures for both species; in this case, the dispersion relation then becomes :

ω 2 c,e -(ω -v E,0 • k) 2 - m i m e (ω -v i,0 • k) 2 = 0 (6.29)
A 2D plot of (6.29) in the (k x , k y ) plane shows this mode is constant in k x ; the instability is thus purely azimuthal, which leads us to focus only on what happens in k y . Then, we can rewrite (6.29) as follows :

ω 2 c,e -(ω -v E,0 k y ) 2 - m i m e ω 2 = 0, (6.30)
a very simple second-order polynomial on ω. Calculating the polynomial's most unstable root and its limit for k y → ∞, one may obtain the following expression :

lim |ky|→∞ γ → ±v E,0 k y 1 + mi me m i m e ≈ ±v
-iω + v e,0 • ik -eTe me 1 ν 2 e,eff +ω 2 c,e (-ν e,eff k 2 ) - k 2 λ 2 D 2 (-iω + v e,0 • ik) + eTe me 1+ k 2 λ 2 D 2 ν 2 e,eff +ω 2 c,e (-ν e,eff k 2 ) = -iω + v i,0 • ik + Ti Te c 2 s 1 ν 2 i,eff (ν i,eff k 2 ) k 2 λ 2 D 2 (-iω + v i,0 • ik) + c 2 s 1- T i Te k 2 λ 2 D 2 ν 2 i,eff (ν i,eff k 2 ) (6.32)
which, if we carefully put T e and T i to zero, gives :

(-iω + v e,0 • ik)(ν 2 e,eff + ω 2 c,e ) - k 2 2 0 2en0 (-iω + v e,0 • ik) -e me ν e,eff k 2 = (-iω + v i,0 • ik)ν i,eff k 2 2 0 2en0 (-iω + v i,0 • ik)ν i,eff + e mi k 2 . (6.33)
Inverting expression (6.33) and developing the expressions for ν e,eff and ν i,eff , one can get :

k 2 0 en 0 = -ek 2 me (iω -v e,0 • ik) 2 + ω 2 c,e - ek 2 mi (iω -v i,0 • ik) 2 , (6.34)
which multiplying by en 0 / 0 and considering only the azimuthal direction, leads to the final expression :

1 - ω 2 pi ω 2 - ω 2 pe (ω -kv e,0 ) 2 -ω 2 c,e = 0 (6.35)
which is none other than the dispersion relation for the Buneman instability [START_REF] Buneman | Instability of electrons streaming through ions across a magnetic field[END_REF][START_REF] Buneman | Excitation of field aligned sound waves by electron streams[END_REF], a low-temperature plasma limit for the kinetic ECDI. This time, as shown in figure 6.6, one can see that the non-neutrality effects contribute to stabilize the growth rate in this specific case. Note that in general for this case, the growth rate is maximized around the resonance point [START_REF] Janhunen | Evolution of the electron cyclotron drift instability in two-dimensions[END_REF]: In fact, Poisson equation is the best solution to solve all small scales issues in our configuration, and at the same time, add some real physical effects that correctly adjust the current physics of the model and get rid of the quasi-neutrality assumption. However, it is really tough to numerically implement Poisson equation in a formalism such as MAGNIS; thus we propose in the following section a way to capture some non-neutrality effects which can be handled numerically. We named it the effective potential.

kv e,0 = ω c,e . (6.36)

Effective potential

Definition and dispersion relation

Here, we propose a method to get rid of small scales effects in our fluid model which is to add a diffusive term that comes under the form of an effective potential in the electron dynamics.

For this, the idea is to take into account the effect of non-neutrality only in the pressure term of the electron momentum equation, allowing the electron density to be different from the ion density via Poisson's equation: n = n i n e = n i + 0 e ∇ 2 φ, (6.37) which, when injected into the electron pressure term, leads to the following momentum equation : where we write the effective potential φ eff as follows :

φ eff = φ - 0 T e en ∇ 2 φ = φ -λ 2 D ∇ 2 φ. (6.39)
This allows to capture some effects induced by Poisson equation in a way that ions react to the basic potential φ while electrons see the effective potential φ eff . If we refer to our general dispersion relation with non-neutrality effects in (6.28), we can recover a dispersion relation with the effective potential if we remove the non-neutral terms in the electron continuity (6.24) and in the ion continuity and momentum equations (6.25),(6.27); as a result, we obtain the following dispersion relation : This dispersion relation contains the effective potential effects via the diffusive term λ 2 D k 2 , with λ D none other than the Debye length. From its expression, it is easy to see that we will observe effects from this term in a linear analysis when :

Effects of the effective potential in linear issues

In this section, we show how this effective potential described previously enables to solve some problematic linear cases. Here, we show a result for the Farley-Buneman instability (the classical form without the effective potential is shown above figure 6.1). This diffusive term happens to be very effective to damp the growth rate at high K but seems to have an impact on the dynamics of the system since the values are not preserved for the observed quantities.

However, it is harder to find a case where the effective potential is effective in the case of gradient-drift instability; this particular diffusive term may work for a narrow range of cases, but not in general.

Effective potential effects in MAGNIS

In this subsection, we aim to detail briefly how we added the effective potential in our code and see some results of its effects in the non-linear regime for the same problematic case than treated in the ion viscosity section.

Referring to equation (6.38), it is easy to see how the momentum equation for electrons is modified in MAGNIS to take into account the effective potential. From there, MAGNIS will proceed to solve φ eff with the electron system of equations, and then the real potential seen by the ions φ thanks to equation (6.39). The density n is solved with the ions equations. Here, λ D is an input parameter. We show the results figure 6.8. Also, we show how ions and electrons do not perceive the same potential as mentioned before, since φ (seen by the ions) corresponds to a smoothed version of φ eff (seen by the electrons), as one can see figure 6.9.

Conclusions and prospects

Before giving the general conclusions from this thesis, it is useful to remind its objectives. This thesis aims to explore the capabilities of fluid models to simulate instabilities in magnetized low-temperature plasma devices such as negative ion sources, Hall thrusters, magnetrons, and so on. The most accurate modeling method to describe these devices is by kinetic particle-in-cell simulation but this is very costly computationally and often out of reach even for modern computers. Alternatively, fluid models are more efficient computationally and have shown their effectiveness in many studies of these magnetized plasma devices, but they are not accuracte and fail to capture certain physical phenomena; in particular, it is very questionnable if they can capture the different plasma instabilities that play an important role in many of these devices. In this thesis we considered one fluid simulation code in particular, the MAGNIS code developed at the LAPLACE laboratory. This code already served to characterize the operation of several ion sources (e.g. for neutral beam injection) and in many cases, plasma instabilities were observed in the simulations results, but the physical revelance of these instabilities was not clear. The questions that this thesis set out to answer are then:

• Are the instabilities observed in MAGNIS proper solutions of the physical model and free from numerical artifacts?

• What are the physical mechanisms of these instabilities and how do they relate to the basic instability modes known from the plasma literature?

• What is the effect of different model approximations and fluid closures on the behavior of these instabilities, and is it possible to improve the model so that instabilities are better described?

• What is the effect of these instabilities on the overall dynamics of the plasma, and in particular, do they generate anomalous transport, and how much?

In order to answer these questions, we revisited the theory of plasma models and magnetized plasma instabilities, carried out a detailed linear analysis of the MAGNIS equations, implemented diagnostics to characterize the simulated instabilities and anomalous transport, performed and analyzed numerous numerical simulations, and investigated several modifications of the physical model. In the next subsection, we give a detailed summary of all this work in the chronological order of the thesis, after which we highlight the main results and conclusions in a separate subsection.

Chronological summary of the thesis

In the introduction of this thesis, we gave a brief overview of magnetized low-temperature plasma sources, which are main topic of our interest, along with the history of instabilities in plasma physics, that were studied in many applications of fusion and space plasmas. We then defined the scope of this thesis. In Chapter 1, we reviewed different elements of basic plasma physics that are essential to this thesis, presented the main general plasma modeling approaches, and introduced various physical quantities and basic equations used in the rest of the thesis.

This general introduction to plasma modeling was necessary to be able to properly present, in Chapter 2, the physical principles of the numerical fluid model MAGNIS used in the thesis. For a 2D Cartesian domain perpendicular to the magnetic field, MAGNIS solves self-consistently the continuity and momentum equations for ions and electrons with an energy equation for electrons, all coupled by the quasi-neutrality assumption. We specified that MAGNIS includes the full inertia terms for electrons and ions and avoids the common drift-diffusion or low-frequency approximations so that it can be applied over a wide range of conditions. We also outlined the main numerical principles of MAGNIS and presented a few examples of simulation results to illustrate its usefulness, for example to help understand the Hall effect in magnetic filter sources. We concluded this chapter by illustrating the magnetized plasma instabilities previously observed in MAGNIS simulations, whose origin was not clear then.

Chapter 3 was intended to set the basis for waves and instabilities that are likely to develop in the plasmas our of interest; and introduced some some generalities about waves and their characteristics and dispersion relations. We then outlined the basic eigen-modes from the general equations of our system, as well as the mechanisms that destabilize these eigen-modes, and finally identified and described the main instabilities found in our configuration.

In Chapter 4, we presented the linear analysis based on a simplified version of MAGNIS geometry of a Cartesian plane perpendicular to the magnetic field with periodic boundary conditions in the azimuthal direction and imposed equilibrium density and potential profiles in the axial direction.

Based on this simple model, we linearized the fluid equations in MAGNIS assuming a constant and uniform magnetic field and species temperatures; from this, we were able to build a dispersion relation related to our simple model system, in which the linear unstable modes introduced in Chapter 3 and observed in the plasma sources of our interest can be recovered. We then proceeded to study their behaviour as a function of their respective key parameters; this theoretical study permitted to identify an appropriate set of benchmark cases for MAGNIS.

After that, we made some required modifications in MAGNIS to force it in a linear regime and to implement the tools necessary for the diagnostics. We simulated the benchmark cases and compared their measured and numerical growth rates, real frequencies and wave numbers with the ones predicted by our dispersion relation. Overall, we observed a very good agreement between them. We speculated that the occasional slight differences between them could be due to a difference in the axial boundary conditions (non-periodic in the simulations). This linear study reassured us on the physical origin of these instabilities since we showed they clearly are a result of the system of equations and not numerical artifacts.

Chapter 5 then followed with the study of the evolution of linear instabilities in the non-linear regime; we focused on four cases, a first very simple one, where we uncoupled every input parameter of the system and kept them uniform and constant, and gradually coupled them and even considered a profile for some (the magnetic field especially), so that we can understand each mechanism effects separately.

After that, we studied the non-linear regime through the characterization of important properties such as the wavelength, propagation direction, and wave frequency and velocity, for which we implemented the required diagnostics into MAGNIS.

The whole study of these four different systems enabled to highlight the main trends in every case, in particular that the wave velocity is close the that of the ion sound mode. Also, we highlighted that the cases where MAGNIS performs best are generally the ones with a uniform magnetic field, with or without the energy equation. On the other hand, we saw that the combination of an inhomogeneous magnetic field with self-consistent ionization and electron temperature seems to be more problematic since we saw the appearance of instabilities with structures that are getting smaller as we refine the mesh, in other words, mesh convergence problems. Moreover, we proceeded to quantify an anomalous transport via an effective mobility; we showed for the magnetically uniform cases (case A and C) that they follow Bohm's predictions. As for the last case, we questioned the concept of an anomalous transport being expressed via a mobility since we were not able to correctly determine an anomalous transport for this configuration.

After that, we considered more realistic conditions such as the magnetized plasma column and the Hall thruster. Concerning the magnetized plasma column, we showed that MAGNIS predicts different instability regimes as a function of magnetic field, probably related to magnetization of the ions, and consistent with preliminary results from PIC simulations. Also in this configuration the wave propagate with appromately the ion sound speed, in the direction of the diamagnetic electron drift. We charactized the anomalous transport in this configuration by means of an effictive global transport time and found that for a certain range of magnetic field values, this was again close to predictions from the anomalous Bohm diffusion theory. However, we showed that the Hall thruster case, close to the problematic last Case D, leads to even more dramatic results than the latter since the simulations diverge as we gradually refine the mesh; this issue clearly highlighted the limitations in fluid models for this specific case, that we discussed in Chapter 6.

Chapter 6 was meant to be an opening to further developments of fluid models in order to cope with the limitations induced by the many approximations and assumptions made to build such a formalism; we listed these assumptions and explained the drawbacks of the fluid model in general, and then focused on our own fluid model and proceeded to identify its main weaknesses, namely the absence of kinetic effects that play a fundamental role on the instabilities saturation, and the quasi-neutrality assumption, erasing some local non-neutrality effects that also are important for the stabilization of unstable cases with small scales effects. For that purpose, we proposed some terms we may implement to capture these missing effects, such as the ion viscosity, that can be perceived as some added ion kinetic effects via a Landau damping, which shows to be effective in both linear and non-linear regimes. We also proposed to study linearly the non-neutrality effects via a general dispersion relation in which we added Poisson equation; in the case of a problematic Bunemann instability when quasi-neutrality is first assumed, we showed that when Poisson equation effects are added, this unstable mode transits towards a stabilizing linear electron cyclotron drift instability. However, the coupling with Poisson's equation being particularly tough to handle numerically, we proposed to capture some of its effects via the electron pressure tensor, which leads to a diffusive term added to the electron momentum equation that we call an effective potential due to the impact on the basic potential; so far, this term however is hardly effective for every regime, but can be an interesting alternative to work further.

Main findings and conclusions

The most important and original findings from this thesis can be summarized as follows:

• In order to make a detailed comparison between self-consistent fluid simulations and linear dispersion relation analysis, one needs to consider a carefully prepared model system such that perturbations of the equilibrium solution take the exact form of harmonic waves. We showed that this is possible for the MAGNIS equations if one considers a system with an exponential equilibrium plasma density profile and different source terms for electrons and ions. This approach allowed us to actually benchmark MAGNIS against a general dispersion relation for this system, and thus, verify its numerical capabilities.

• The instabilities observed in MAGNIS in the linear regime are proper solutions of the physical model system; their behavior is in good agreement with the general dispersion relation that we constructed, provided that the numerical mesh is sufficiently refined to resolve the most unstable wave length (at least 20 grid points per period).

• Depending on the configuration, it can happen that no mesh convergence is possible; the size of the instabilities is then determined by the numerical grid even if this is refined indefinitely. The origin of this problem is not numerical but mathematical: the system of model equations is linearly most unstable at infinite k. This unphysical property of the model is due to the fact that some of its physical approximations are not valid at small scales, in other words, the small scale physics is not properly included. We identified several approximations that can lead to high-k issues and should therefore better be avoided, such as the quasi-neutrality assumption (used in MAGNIS for numerical reasons) and the neglect of electron inertia (optional in MAGNIS). We also demonstrated that the high-k problem can be prevented by adding a viscosity term in the ion momentum equation, but more work is needed to justify this physically.

• The main instability mechanisms are resulting from the coupling between the electron E × B drift, the ion flow across the magnetic field, and the space gradients present in the system such as the density, electron temperature and magnetic field gradients, along with the collisional processes between the charged species and the neutrals.

• When the instabilities evolve into a nonlinear regime and saturate, their dominant wave length can be quite different from that observed when they first appear in the linear regime; typically it is larger. However, for all configurations we investigated, their wave velocity was close to the ion sound speed and directed in the direction of the electron magnetic drift, exactly as predicted by the linear analysis.

• In plasma configurations with a uniform magnetic field and diffusive transport driven by the pressure gradient, the nonlinear fluid instabilities generate an anomalous transport that is well fitted by the Bohm diffusion formula, when averaged over the plasma volume. However, this holds only up to some value of magnetic field (≈ 150 G for the hydrogen plasma column we studied), above which the instabilities and anomalous transport first decrease and then give way to larger nonlinear structures, related to magnetization of the ions.

• In plasma configurations with a magnetic barrier, it was sometimes impossible to characterize the anomalous transport by an effective mobility in a meaningful way.

• MAGNIS in its present form performs best for configurations with diffusive transport in a uniform magnetic field (plasma column), for which there are relatively few numerical converges issues and the results seem consistent with PIC simulations. However, configurations with a magnetic field barrier can be problematic, in particular when a strong voltage is applied (Hall thruster).

Prospects

Although this thesis made a good start with the study of fluid simulation of instabilities in magnetized plasma sources, there are still many open questions and issues that were beyond the scope of this thesis and should be addressed in future work. We foresee the following future developments for this study:

• In the linear analysis of this thesis, we assumed a constant uniform electron temperature and magnetic field. It would be interesting to extend the linear analysis to include the electron energy equation and a non-uniform magnetic field profile (or at least a magnetic field gradient). For these extensions it will most probably not be possible to keep fully self-consistent equilibrium profiles as we did in this thesis, i.e. it will be necessary to make the Boussinesq approximation k L -1 T,B , with L T,B the gradient length of the equilibrium electron temperature or the magnetic field. Hence, it will be difficult to use the analytical results for detailed benchmarking of MAGNIS, but the analysis will still be interesting to get better understanding of the instabilities. From the MAGNIS simulation results in this thesis, we expect that the magnetic field gradient in particular leads to significant modifications of the unstable modes. The electron energy equation on the other hand does not seem to cause essential differences, but all this needs to be confirmed.

• In order to solve the mesh convergence issues and instability saturation issues that we observed in the MAGNIS simulations, it is essential to modify the system of model equations. In Chapter 6 we explored several of such modifications (addition of ion viscosity, Poisson equation, and effective potential) but this work is far from finished. Preferably one would like to make modifications on physical grounds, so that not only they ensure good mesh convergence and nonlinear saturation, but also they represent real physical mechanisms in a self-consistent way. One obvious modification that should be made to MAGNIS is to remove the quasi-neutrality assumption and add the Poisson equation to the system. As we showed, this will stabilize the length scales below the Debye length. However, adding the Poisson equation appears to be a big numerical challenge, which probably requires to redesign the numerical schemes and time integration cycle of MAGNIS. In order to investigate this further, a first test with adding the Poisson equation was made during this thesis for a simplified 1D version of MAGNIS, only solving the model in the azimuthal direction. This test was indeed successful, meaning that we managed to perform 1D simulations with Poisson that were in agreement with the linear dispersion relation, but we had to use a very small numerical time step and also make use of the 1D simplifications of the numerical system. Much more work is needed to extend this approach to the full 2 1 2 D MAGNIS code. Furthermore, in these 1D simulations we observed that although the Poisson equation stabilizes the small scales, it does not prevent nonlinear saturation problems (e.g. for the conditions of the Hall thruster). Other modifications of MAGNIS that should be further explored concern the closure relations, in particular regarding the pressure tensor, which could probably be extended to include specific physical mechanisms that play a role at small scales (e.g. Larmor radius) or for the saturation of the instabilities.

• To check the physical validity of the MAGNIS simulation results, it is necessary to compare them with data obtained by other methods. A first important step is to compare them with results from PIC simulations. In fact, systematic comparisons between MAGNIS and PIC simulations at the LAPLACE laboratory are already underway for a magnetized plasma column configuration (see figures 5.23 and 5.24 for preliminary results from this comparison). A next step will be to make systematic comparisons with experimental results (e.g. results obtained on the RAID device [START_REF] Furno | A novel helicon plasma source for negative ion beams for fusion[END_REF]).

Figure 1 -

 1 Figure1-Hall thruster used for space propulsion. In this device, the neutral gas (Xenon) is ionized and accelerated by an applied electric field and the electron current trapped in a magnetic field.

Figure 2 -

 2 Figure 2 -Magnetron devices are widely used for material processing and surface modification. Crossed electric and magnetic field configuration is used to confine electrons creating dense plasma. The configuration and many physics processes are similar to those in the Hall thruster.

Figure 3 -

 3 Figure 3 -Magnetized plasma column in the CYBELE device [1, 2] (left) and a ITER negative ion source [3] (right) use magnetic filter configuration to extract ions.

Figure 4 -

 4 Figure 4 -Plasma discharge in European Joint Tokamak JET (left), and schematic of the tokamak device (right).The fully ionized plasma is confined in a magnetic "cage" formed by poloidal and toroidal magnetic field to keep it away from the walls, and thus, to prevent it from cooling and damaging the chamber.

Figure 5 -

 5 Figure 5 -Representation of the plasma in equilibrium state, Sausage instability (m=0) and Kink instability (m=1) in a torus portion.

Figure 6 -

 6 Figure 6 -2D PIC simulations of a channel region in Hall thruster (right): Azimuthal instabilities in the azimuthal electric field (right-top) and density (right-bottom) profiles. the minimum values are -5 × 10 4 V.m -1 for Ey and 0 m -3 for n, and the maximum values are respectively 5 × 10 4 V.m -1 and 5 × 10 17 m -3 . The domain simulation is represented in the left.from Ref. [31].

Figure 7 -

 7 Figure 7 -PIC (top) and fluid (bottom) simulations of a RAID (Resonant Antenna Ion Device)[START_REF] Furno | A novel helicon plasma source for negative ion beams for fusion[END_REF][START_REF] Furno | Helicon wave-generated plasmas for negative ion beams for fusion[END_REF] plasma source for a magnetic field of 200G[START_REF] Fubiani | [END_REF].

Figure 1 . 1 -

 11 Figure 1.1 -Typical frequency and length scales present in the magnetized plasma sources of our interest.

Figure 1 . 3 -

 13 Figure1.3 -The motion of positive ions and electrons in uniform the electric and magnetic fields case. The particle drifts in the direction perpendicular to both the electric field and the magnetic field with the drift velocity (1.17) and this phenomenon is called E × B drift.

Figure 1 . 4 -

 14 Figure1.4 -The motion of electrons in uniform electric and magnetic fields for a collisional case. The particle drifts in the E × B direction and is deviated from its initial trajectory when colliding with another particle in the case of an elastic scattering.

Figure 1 . 5 -

 15 Figure 1.5 -Image of the physics near the wall. At the top, the density drops in the sheath (and pre-sheath to a lesser extent) until the medium becomes non-neutral due to an excess of positive charges from ions (and the deceleration of electrons which cannot get through the sheath because of the potential drop), while the profile of the increasing ion velocity in the sheath illustrates the Bohm criteria. The bottom image shows also the potential drop, reaching the value V f called the floating potential.

Figure 1 . 6 -

 16 Figure 1.6 -Numerical cycle of the PIC code.

Figure 2 . 1 -

 21 Figure 2.1 -Rate coefficients km, kiz and κ (see equation (2.13)) for electron processes in molecular hydrogen (left) and in argon (right). Note that kiz has been multiplied by 10.

Figure 2 . 2 -

 22 Figure 2.2 -Image of the grid used on the 2D simulation plane, where appears the control volume. As explained above, scalars are defined at the nodes (dots), and vector components in the middle between two nodes, at edges of the control volume (arrows).

  and the velocities are corrected from(2.40). Equation (2.41) is an anisotropic elliptic equation for φ which is discretized in a straightforward way by the control volume method with central differences for ∇ φ and b × ∇ φ; the resulting nine-point linear system is solved by a self-made iterative multigrid solver that can handle a large number of grid points. The elliptic character of equation (2.41), essential for the functioning of the solver, requires that the ∇ φ terms are dominant over the b × ∇ φ terms, which is always true here because from equation (2.35) we have |A| > |B| thanks to δ ≥ |ω|.

Figure 2 . 3 -

 23 Figure 2.3 -Two schematic views of the same magnetic filter type ion source.

Figure 2 . 4 -

 24 Figure 2.4 -MAGNIS simulation results illustrating the Hall effect in a magnetic filter source geometry.

Figure 2 . 5 -

 25 Figure 2.5 -PIC simulation of in a magnetic filter source operated in hydrogen: electron density (left), electron temperature (middle) and electron flux (right). From [63].

Figure 2 . 6 -

 26 Figure 2.6 -MAGNIS simulation results illustrating the plasma instabilities that arise when the magnetic drift is closed via periodic boundary conditions.

Figure 2 . 7 -

 27 Figure 2.7 -Ion density (left) and potential (right) profiles in PIC simulations of a hydrogen plasma in a magnetic filter with periodic boundary conditions. Figures from [66].

  The electron momentum in(3.47) gives an expression for the fluctuating velocity :

Figure 3 . 1 -

 31 Figure 3.1 -Instability mechanism : example of potential and density profiles in a Simon-Hoh case.

)

  The crucial modification is the electron momentum equation in the form m e n e ∂v e ∂t = en e (-∇φ + 1 c v e × B) -∇p e -m e n e νv e ,

80 )

 80 Note that the second term in this equation generalizes the electron diffusion and mobility in (3.50) with time dependent term, which is the effect of the inertia. Using (3.80) in the electron continuity equation (3.46) one obtains the expression for the perturbed electron density in the form

Figure 3 . 2 - 2 e 1 , 2 e 1 2 e 1 ,

 32212121 Figure 3.2 -The gradient drift (Simon-Hoh) instability, transition to the destabilized lower hybrid mode, and to the ion sound mode , note that ωr ∼ ky . (a) The ion-sound mode Destabilized by the density gradient and collisions, v * = 1.3 10 4 , νe = 10 6 , ω0 = 0, (b) Destabilization by density gradient, E × B drift and collisions, v * = 1.3 10 4 , vE0 = 5 10 5 νe = 10 6 . Reproduced with permission from A Smolyakov et al, PPCF 2016

Figure 4 . 1 -

 41 Figure 4.1 -Model 2D geometry

  )

Figure 4 . 2 -

 42 Figure 4.2 -Stationary density and potential profiles in the axial direction, as in the domain (Fig. 4.1)

Figure 4 . 3 -

 43 Figure 4.3 -A tree illustrating the instability mechanisms included in our dispersion relation.

Figure 4 . 4 -

 44 Figure 4.4 -Representation of the normalized growth rate and real frequency (defined in (4.43)) in the normalized (Kx, Ky) (defined in (4.48)) plane for the particular case of a pure gradient unstable mode. The unstable region of the growth rate is delimited by the red dots.

Figure 4 . 5 -

 45 Figure 4.5 -Representation of the stability criterion for normalized growth rate, frequency (defined in (4.43)) and wave numbers (defined in (4.48)) in the dimensionless (he, Σ) parameters plane, defined in (4.45) and (4.46) respectively.

Figure 4 . 6 -

 46 Figure 4.6 -Representation of the stability criterion for normalized growth rate, frequency (defined in equation (4.43)) and wave numbers (defined in (4.48)) in the dimensionless (he, Σ) parameter plane (defined in (4.45) and (4.46)) when imposing Kx = 0 in order to reproduce the Simon-Hoh configuration.

Figure 4 . 7 -

 47 Figure 4.7 -The Farley-Bunemann instability normalized growth rate (top) and wave frequency (bottom), defined in (4.62), in the normalized (Kx, Ky) plane, defined in (4.61), in the general case.

Figure 4 . 8 -

 48 Figure 4.8 -The Farley Buneman normalized growth rate, frequency and wave numbers on the normalized (ψ, Vnorm = vei/cs) plane, as defined in (4.59) and (4.60) respectively, representative of the stability criteria.

Figure 4 . 9 -

 49 Figure 4.9 -Case 3 : the normalized growth rate and real frequency (defined in (4.83)) in the normalized (Kx, Ky) plane, defined in (4.82).

Figure 4 . 10 -

 410 Figure 4.10 -Case 3 : the normalized growth rate and real frequency (defined in (4.83)) in the normalized (Kx, Ky) plane with the electron time derivative inertial term.

Figure 4 . 11 -

 411 Figure 4.11 -Case 3: the normalized growth rate and real frequency (defined in (4.83)) in the normalized (Kx, Ky) plane with the electron time derivative inertia term and in a large E limit.

Figure 4 . 12 -

 412 Figure 4.12 -Case 3: the normalized growth rate and real frequency (defined in (4.83)) in the normalized (Kx, Ky) plane with all electron inertial terms.

Figure 4 . 13 -

 413 Figure 4.13 -An example of how we select our (E,B) values, here for the Case 1.

  4.14 and fig.4.15 in three different mesh grids.

Figure 4 . 14 -

 414 Figure 4.14 -Azimuthal electric field Ey.

Figure 4 . 15 -

 415 Figure 4.15 -Azimuthal electric field Ey.

Figure 4 . 16 -Figure 4 . 18 -

 416418 Figure 4.16 -Growth rate values for different values of E Figure 4.17 -Real frequency values for different values of E

Figure 4 . 19 -Figure 4 . 20 -

 419420 Figure 4.18 -Growth rate values for different values of B.Figure 4.19 -Real frequency values for different values of B

Figure 4 . 21 -

 421 Figure 4.21 -Case 3bis (b) : Growth rate values in function of E.

Figure 4 . 22 -

 422 Figure 4.22 -Case 3bis (a) : Real frequency values in function of B.

Figure 4 . 23 -

 423 Figure 4.23 -Case 3bis (b) : real frequency values in function of E.

Figure 5 . 1 -

 51 Figure 5.1 -Case B input profiles : the artificial ionization and magnetic field profiles for an isothermal case, where Te = 6 eV.

Figure 5 . 2 -

 52 Figure 5.2 -Case D input profiles : the injected heating power and magnetic field profiles for selfconsistent electron temperature and ionization.

Figure 5 . 3 -

 53 Figure 5.3 -Simulation results of the azimuthal component of the electric field Ey for Case A forced in a linear regime at B=25 G.

Figure 5 . 4 -

 54 Figure 5.4 -Case A : density and potential profiles on 256x256 and 512x512 mesh.

. 28 )Fig. 5 .

 285 Fig.5.7 shows that our average values of each different measured mobility profile closely match the Bohm mobility. The agreement is remarkably good.

Figure 5 . 6 -

 56 Figure 5.6 -Case A : Effective and classical mobility profiles for each magnetic field value.

Figure 5 . 7 -

 57 Figure 5.7 -Case A : Comparison of the domain-averaged effective mobility and the Bohm mobility.

Figure 5 . 9 -

 59 Figure 5.9 -Case B : Density and potential profiles on 256x256 and 512x512 mesh.

Figure 5 .

 5 Figure 5.10 -Stationary 1D and averaged 2D axial profiles for the density (left) and the potential (right).

Figure 5 . 11 -

 511 Figure 5.11 -Classical and anomalous mobilities calculated from MAGNIS simulations on a 256x256 (left) and 512x512 (right) numerical mesh.

Figure 5 .

 5 Figure 5.12 -Case C : density, potential and electron temperature profiles in function of the magnetic field values in a 512x512 mesh grid.

Figure 5 .

 5 Figure 5.13 -Case C : density, potential and electron temperature profiles in function of the magnetic field values in a 512x512 mesh grid.

Figure 5 .

 5 Figure 5.14 show the stationary 1D axial profiles and the averaged 2D ones for the density, potential and electron temperature. No differences with Case A stand out, except for the presence of a profile for the electron temperature and its effects on the values (as mentioned just before). Here again, the differences between the averaged and stationary quantities underline the presence of an anomalous transport responsible of a change in the profiles of the observed variables.

Figure 5 . 14 -

 514 Figure 5.14 -Stationary and averaged 1D profiles for the density, the potential and the electron temperature.

Figure 5 .

 5 Figure 5.15 -Case C : classical and effective mobility profiles in 256x256 and 512x512 mesh grids.

Figure 5 .

 5 Figure 5.16 -Case C : averaged values of the effective mobility (black) and values of the mobility deduced from Bohm's theory (red) for the different magnetic field values.

Figure 5 .

 5 Figure 5.17 -Case D : density, potential and electron temperature profiles for the first three magnetic field values in a 256x256 mesh grid.

Figure 5 .

 5 Figure 5.18 -Case D : density, potential and electron temperature profiles for the first three magnetic field values in a 512x512 mesh grid.

Figure 5 . 19 -

 519 Figure 5.19 -The magnetized plasma column geometry used in MAGNIS.

Figure 5 . 20 -

 520 Figure 5.20 -Instantaneous profiles of the electron density, electron and ion flux vectors, electric potential and electron temperature, respectively (rows of the figure), for different magnetic field values (columns), simulated on a 256x256 mesh.

Figure 5 . 21 -

 521 Figure 5.21 -The measured wave velocity for different magnetic field values, along with ion sound velocity, for the column (black and red) and the periodic case C (blue and cyan).

  + ν m , (5.31) where ν ef f is an effective collision frequency gathering the classical collisions (ν m ) and the anomalous Bohm collisions. Naturally, the global transport time including Bohm transport increases approximately like B. The 1D stationary classical and Bohm results are shown by the dashed curves in Figure 5.22.

Figure 5 . 22 -

 522 Figure 5.22 -Global plasma transport times as defined in equation (5.30) for different MAGNIS simulations as a function of magnetic field and for different numerical meshes. Left panel: plasma column. Right panel: periodic Case C.

Figure 5 . 23 -

 523 Figure 5.23 -Plasma density in PIC simulations of the plasma column case with B =400 G [48].

Figure 5 . 24 -

 524 Figure 5.24 -Plasma density in MAGNIS fluid simulations of the plasma column case with B =400 G.

Figure 5 . 25 -

 525 Figure 5.25 -Basic and simplified representation of a Hall thruster's considered geometry : the channel.

Figure 5 . 26 -

 526 Figure 5.26 -Input profiles in MAGNIS for the Hall thruster case.

Figure 5 .

 5 Figure 5.27 -Density, potential and electron temperature profiles for a self-consistent case of an ideal version of the Hall thruster.

Figure 5 . 28 -

 528 Figure 5.28 -Density and potential profiles for an isothermal case of an ideal version of the Hall thruster.

  2 the Farley-Buneman case with and without the diffusive term; fig.6.1 shows the standard result of Farley-Buneman dispersion relation and fig.6.2 the same unstable mode with the addition of the ion viscosity.

Figure 6 . 1 -

 61 Figure 6.1 -Normalized real frequency and growth rate (Ω = ω/νi) in the normalized (Kx, Ky) plane (K = csk/νi) for the classical Farley-Bunemann instability.

Figure 6 . 2 -

 62 Figure 6.2 -Normalized real frequency and growth rate (Ω = ω/νi) in the normalized (Kx, Ky) plane (K = csk/νi) for the classical Farley-Bunemann instability with the ion viscosity.

Figure 6 . 3 -

 63 Figure 6.3 -Normalized real frequency and growth rate (Ω = ω/νi) in the normalized (Kx, Ky) plane (K = csk/νi) for the classical gradient-drift instability.

Figure 6 . 4 -

 64 Figure 6.4 -Normalized real frequency and growth rate (Ω = ω/νi) in the normalized (Kx, Ky) plane (K = csk/νi) for the gradient-drift instability with the ion viscosity.

Figure 6 . 5 -

 65 Figure 6.5 -Case D for Bmax = 100G in a 256x256 mesh grid, without ion viscosity (a) and with ion viscosity (b).

Figure 6 .

 6 Figure 6.5 shows positive effects of the ion viscosity in this case; indeed the profiles appear cleaner, with their initial values and trends conserved.

6. 3

 3 Deviation from quasi-neutrality : Poisson equation and effective potential 6.3.1 Poisson equation

Figure 6 . 6 -

 66 Figure 6.6 -The normalized growth rate (Ω = ω/ωLH ) in function of the normalized azimuthal wave vector (K = kρe) in the quasi-neutral case (dashed line) and with Poisson effects (red line).

α 1 n ∂v e ∂t + α 2 m e ∇n e = en m e ∇ φ - 0 T e en ∇ 2 φ

 22 nv e ∇ • v e + nν m,e v e = en m e (∇φ -v e × B) -eT e v e × B -eT e m e ∇n = en m e (∇φ eff -v e × B) -eT e m e ∇n (6.38)

(

  -iω + v e,0 • ik)(ν 2 e,eff + ω 2 c,e ) -T e e me (-ω c,e g • ik × b -ν e,eff (-g • ik + k 2 )) (1 + λ 2 D k 2 ) e me (-ω c,e g • ik × b -ν e,eff (-g • ik + k 2 )) = (-iω + v i,0 • ik)(ν 2 i,eff + β 3 ω 2 c,i ) e mi (-β 3 ω c,i g • ik × b + ν i,eff (-g • ik + k 2 ))+ T i . (6.40)

λ 2 D k 2

 22 >> 1. (6.41) Otherwise, if the opposite happens, then this diffusive term becomes negligible and (6.40) is equivalent to the dispersion relation (4.39) we established chapter 4. Hence, it is important to fix a consequent value of λ D in the input parameters in order to see the effects.

Figure 6 . 7 -

 67 Figure 6.7 -Normalized real frequency and growth rate (Ω = ω/νm,i) in the normalized (Kx, Ky) plane (K = csk/νm,i) for the classical Farley-Bunemann instability with the effective potential.

Figure 6 . 8 -

 68 Figure 6.8 -Case D for B = 50 G with the effective potential in a 256x256 mesh grid. Here, we set λD = 3.125 × 10 -3 m.

Figure 6 . 9 -

 69 Figure 6.9 -The normal potential (left) seen by ions, and the effective potential (right) seen by the electrons.

Table 4 .

 4

1 -Input values for the dimensionless parameters, case 1.

Table 4 .

 4 .2.

	ψ	7.31×10 -4
	h e	4.4×10 2

E 0 c s /(T e + T i )ν i

4.09 2 -Input values for the dimensionless parameters, case 2.

  83), (4.82), (4.80),(4.88), (4.45) and (4.81) respectively, and we introduced a new parameter:

	ν N =	ν i ν e	,	(4.98)
	the ion-over-electron momentum frequency ratio.
	ψ		0.35
	h e		1.76
	E 0 c s /(T e + T i )ν i 0.82
	G		2.45

Table 4 .

 4 3 -Input values for the dimensionless parameters, case 3.

Table 4 . 4

 44

-Input values for case 1.

Table 4 .

 4

	6 -Input values for case
	3bis (b).

Table 4 . 7

 47

	5	157.1	209.4	4

-table of the expected values.

Table 4 .

 4 table 4.7, while σ num is the numerical value obtained in MAGNIS. Case 1 ω r,num (s -1 ) γ num (s -1 ) k x,num (m -1 ) k y,num (m -1 )

	Mesh 128	-2.396×10 6	4.439×10 5	142.1	209.4
	Mesh 256	-2.398×10 6	4.592×10 5	142.1	209.4
	Mesh 512	-2.398×10 6	4.623×10 5	142.1	209.4
	Error percentage	δ ωr	δ γ	δ kx	δ ky
	Mesh 128	1.2 %	11.8 %	9.5 %	0 %
	Mesh 256	1.1 %	8.8 %	9.5 %	0 %
	Mesh 512	1.1 %	8.1 %	9.5 %	0 %
	Mesh 128	-3.411×10 6	1.288×10 6	287.6	628.3
	Mesh 256	-3.355×10 6	1.498×10 6	301.6	628.3
	Mesh 512	-4.278×10 6	1.719×10 6	382.9	785.4
	Error percentage	δ ωr	δ γ	δ kx	δ ky
	Mesh 128	24 %	11.6 %	39 %	20 %
	Mesh 256	25.4 %	2.8 %	36 %	20 %
	Mesh 512	4.9 %	18 %	18.7 %	0%

8 -table of the obtained values.

Case 3bis (a)

ω r,num (s -1 ) γ num (s -1 ) k x,num (m -1 ) k y,num (m -1 )

Table 4 .

 4 9 -table of the obtained values.

	Mesh 128	-3.401×10 6	3.313×10 6	168.3	471.2
	Mesh 256	-3.434×10 6	3.507×10 6	175.6	471.2
	Mesh 512	-3.361×10 6	3.513×10 6	183.1	471.2
	Error percentage	δ ωr	δ γ	δ kx	δ ky
	Mesh 128	1.5 %	7.6 %	7.2 %	0%
	Mesh 256	2.5 %	2.2 %	11.8 %	0%
	Mesh 512	0.3 %	2.1 %	16.5 %	0%

Case 3bis (b)

ω r,num (s -1 ) γ num (s -1 ) k x,num (m -1 ) k y,num (m -1 )

Table 4 .

 4 10 -table of the obtained values.

Table 5 .

 5

	1 )	2.29×10 8	5.27×10 8	7.21×10 8	9.67×10 8	1.28×10 9
	ω cs (s -1 )	4.27×10 6	6.39×10 6	6.39×10 6	6.39×10 6	8.55×10 6

2 -Case B : Table of characteristic values for region 1.

Table 5 .

 5 3 -Case B : Table of characteristic values for region 2.

	1 )	2.29×10 8	5.27×10 8	7.21×10 8	9.67×10 8	1.28×10 9
	ω cs (s -1 )	1.71×10 7	2.56×10 7	3.41×10 7	3.85×10 7	4.27×10 7

Table 5 .

 5 4 -Case C : Table of characteristic values.

	1 ) 1.335×10 4	1.264×10 4 1.2434×10 4 1.182×10 4	1.121×10 4
	λ instab (m)	8.3×10 -3	6.25×10 -3	5×10 -3	4.5×10 -3	4.2×10 -3
	ω instab (s -1 )	1.011×10 7	1.271×10 7	1.562×10 7	1.650×10 7	1.677×10 7
	v p (m.s -1 )	1.09×10 5	8.46×10 4	3.86×10 4	3.08×10 4	2.56×10 4
	v E (m.s -1 )	2.55×10 4	8.05×10 3	4.15×10 3	2.44×10 3	1.36×10 3
	c s (m.s -1 )	1.7815×10 4 1.6672×10 4 1.6058×10 4 1.5339×10 4 1.4660×10 4
	ρ e (m)	2.45×10 -3	1.15×10 -3	7.35×10 -4	5.27×10 -4	4.03×10 -4
	L n (m)	2.14×10 -2	1.68×10 -2	1.43×10 -2	1.23×10 -2	1.06×10 -2
	ω LH (s -1 )	7.282×10 6	1.456×10 7	2.184×10 7	2.913×10 7	3.641×10 7
	ω c,e (s -1 )	4.397×10 8	8.794×10 8	1.319×10 9	1.759×10 9	2.198×10 9
	ω cs (s -1 )	1.35×10 7	1.68×10 7	2.02×10 7	2.14×10 7	2.19×10 7

Table 5 .

 5 5 -Case D : tables of values for regions 1 and 2.

					Case D	B max =25 G B max =50 G B max =75 G
					Region 2
	1 )	2.13×10 4	2.25×10 4	2.11 ×10 4	v prop (m.s -1 )	1.73×10 4	1.50×10 4	1.20 ×10 4
	λ instab (m)	6.25×10 -3	5×10 -3	4.17×10 -3	λ instab (m)	3.13×10 -3	2.5×10 -3	2.1×10 -3
	ω instab (s -1 )	2.14×10 7	2.83×10 7	3.18×10 7	ω instab (s -1 )	3.47×10 7	3.77×10 7	3.59×10 7
	v e,y (m.s -1 )	1.04×10 5	9.37×10 4	8.62×10 4	v e,y (m.s -1 )	1.01×10 5	6.61×10 4	4.25×10 4
	c s (m.s -1 )	2.30×10 4	2.41×10 4	2.50×10 4	c s (m.s -1 )	1.70×10 4	1.43×10 4	1.20×10 4
	ρ e (m)	4.15×10 -3	2.02×10 -3	1.68×10 -3	ρ e (m)	3.07×10 -3	1.20×10 -3	8.10×10 -4
	L n (m)	4.8×10 -2	4.3×10 -2	3.6×10 -2	L n (m)	4×10 -2	4.2×10 -2	4×10 -2
	ω LH (s -1 )	5.53×10 6	1.19×10 7	1.48×10 7	ω LH (s -1 )	5.53×10 6	1.19×10 7	1.48×10 7
	ω c,e (s -1 )	3.34×10 8	7.21×10 8	8.97×10 8	ω c,e (s -1 )	3.34×10 8	7.21×10 8	8.97×10 8
	ω cs (s -1 )	2.03×10 7	2.44×10 7	2.85×10 7	ω cs (s -1 )	3.41×10 7	3.59×10 7	3.59×10 7

Table 5 .

 5 6 -Table of the input parameters for the column system.

Table 5 .

 5 [START_REF] De | Mhd instabilities in tokamaks[END_REF] -Input values for the Hall thruster case.

	5.26, and the input

  in a quasi-neutral case, the growth rate of this unstable mode is endlessly growing in function of k y without stabilizing; this represents clearly a high-k issue. An illustration of this demonstration is shown fig.6.6. Now, taking into account Poisson equation's effects on the previous dispersion relation, one gets the following new dispersion relation :

	E,0 k y	m e m i	(6.31)
	which shows that

Once the wave velocity is known, we determined the associated frequency simply from ω instab = 2πv prop λ instab (5.14) where λ instab is the dominant wave length which we deduced by eye from:

n y being the number of periods counted in the domain and L y , its size in y.

Diagnostics for the measurement of an anomalous transport

In order to quantify the anomalous transport in each given configuration, we calculate an effective mobility from averaged measured quantities, namely the axial electron flux and velocity, all averaged in the azimuthal direction and in time, as follows (see equation (5.4) above):

(5.16)

An alternative way to way obtain an effective mobility from the simulation results is to first deduce the effective Hall parameter from the ratio of the azimuthal and axial electron fluxes:

(5.18)

We found that both these effective mobility diagnostics gave almost identical results for all the cases. For comparison purposes, we also calculate the classical Hall parameter and classical mobility such that:

(5.20)

Analysis of non-linear modes and anomalous transport

This section shows and discusses the results obtained for each case. We show the simulation results of each case for different numerical mesh grids and for the different magnetic field values, namely the According to these simple equations, the plasma density has a sinusoidal profile, and for a given average source term S x,y (as we have in our simulations), it increases as the inverse of the (effective) magnetized electron mobility. These trends are indeed well reproduced by the purely 1D results in Figure 5.5, with an inverse mobility that in creases as

For the averaged 2D results, the plasma density increases much less with B, a clear sign of anomalous transport (µ -1 ef f µ -1 ⊥e ). Note also the averaged 2D densities do not exactly have a sinusiodal profile but are much more peaked in the center, which indicates that the anomalous transport is not uniform but varies along x. We thus proceed to quantify the anomalous transport responsible of the observed profile trends. For that purpose, we show fig. 5.6 the classical (pure 1D) and effective (averaged 2D) mobilities measured for each case and for the two different numerical meshes. One can see that effective mobility is much higher than the classical mobility and it has a strong spatial profile: it is minimum in the center (but still larger than the classical mobility) and then increases towards the walls. Note that there is a slight difference in the effective mobility between the two meshes, for the higher magnetic field values it tends to be slightly lower on the finer mesh. An idea is to compare our measured mobilities with the anomalous mobilities predicted by the well-known Bohm formula [START_REF] Bohm | The characteristics of electrical discharges in magnetic fields[END_REF]; for this purpose, we calculated a simple average value for each mobility profiles displayed in fig. 5.6 along with the expected Bohm mobility for each different magnetic field value such that :

.

(5.29)

Case B

Unlike Case A, Case B has the particularity to take into account a profile for the magnetic field. Here again, if we take a look at the simulation results fig. 5.9, showing the density and the potential profiles, we obtained a correct mesh convergence. However, for the highest values of the magnetic field, thinner structures appear on the finer mesh; indeed, a careful look to the cases for B max =100 and 125 G on both meshes enables to see these slight differences between the obtained results.

A general observation to be made for this case is that the instabilities do not behave in the same away everywhere in the domain, but that there are two regions with different behavior, before and in the magnetic barrier; first, before the magnetic barrier, which is the creation area, the density is at its maximum and presents less periods in the domain than inside and beyond the barrier, where the density tends to have lower values. The same behavior can be observed for the potential. Thus, we can divide the domain in two different regions where we characterize separately the two observed tendencies; we call region 1 the area which is approximately the first half of the domain and region 2 the second half. We proceed to the characterization below for these two different regions. Fig. 5.8 shows the obtained time-averaged measured wave velocities for each values of the maximum of the magnetic field profile; these profiles highlight too the fact that this case presents two different regions where the wave velocities display different profiles depending on the region we study, contrary to the uniform case A where the absolute wave velocity is constant.