The history of EP starts with the same pioneers than for chemical propulsion (Tsiolkovsky, Goddard, Oberth), over one century ago. But at that moment, the relevant scientific fields (plasma physics, electrical engineering, etc.) were not mature enough and their brilliant ideas have not been applied to build thrusters . It was many decades later, in 1964, that the first electric thruster was fired in space, aboard the SERT-1 spacecraft. This Gridded Ion Thruster (GIT) was later followed by a Hall thruster (HT) in 1972, aboard the Meteor 18 satellite. While these two types of electric thrusters were both studied independently in the US and in USSR, the former decided to use mainly the GIT while the latter focused on the HT development, lead by Morozov . The situation changed in 1991 when the Cold War ended and the HT Russian design plans were made available. It is interesting to note that the French spatial agency (CNES) decided to focus its EP research work on this particular thruster technology and, in collaboration with Safran Aircraft Engines (SAE), they successfully equipped the SMART-I satellite (that reached the moon orbit in 2003) with Hall Thrusters (PPS®1350 ). Since then, almost 600 spacecraft have been launched, with various types of electric thrusters. For a more detailed history of EP, the reader may refer to the review of Choueiri  for the first 50 years, or the more recent review of Lev et al. .
As mentioned ealier, what distinguishes an electric thruster from a chemical one is the decoupling between energization (a plasma is produced by adding en-ergy to a gas) and acceleration (thrust is produced by accelerating the plasma charged particles). Hence, an electric thruster is usually decomposed in three sep-arate zones: the region where the plasma is created (plasma coupling region), the one where it is accelerated (ion extraction/acceleration region) and the thruster plume (ion beam neutralization/detachment region), as seen in figure 1.3. The thruster plume needs to be neutral because charged particles could interact with other subsystems of the satellite and perturb their operation. Moreover, plume neutralization is necessary to avoid space charge region that can be a barrier for the following ions.
Electric thrusters are usually distinguished by the acceleration method used to produce the thrust. Hence, we can list three main categories:
— Electrothermal thrusters: the propellant is heated electrically (by a resistive heater for resistojets or with a high current arc for arcjets) and this thermal energy is then converted in kinetic energy by a nozzle. The propellant can also be heated by a radio-frequency source .
— Electrostatic thrusters: the propellant is accelerated by a static electric field, created by a potential diﬀerence applied between two close parallel grids (for GIT) or between an anode and a cathode at each end of the thruster (for HT).
— Electromagnetic thrusters: the propellant is accelerated by a combi-nation of magnetic and electric fields (i.e. by the Lorentz force). One can mention the Helicon Double Layer Thruster , along with the very promising VASIMR technology .
A broad range of needs
From the diversity of space missions stems the need of widening the thruster power capabilities. The specific impulse and power ranges of the main types of electric thrusters are given in figure 1.4. In the last decade, mainly for con-stellation purposes (like OneWeb), the number of small satellites (micro and nano-satellites) has grown exponentially, and with them, the need of low-power thrusters (10-500 W). Many companies have already developed these smaller thrusters, such as ThrustMe (iodine gridded thruster) and ExoTrail (HT) in France. Moreover, high-power thrusters (5-25 kW) are needed for deep-space missions  and manned flight, to the Moon, Mars, or even beyond.
However, this scaling is often done empirically (trial and error method), which limits the eﬃciency of the newly developed thrusters and increases the development time and cost. Innovative spacecraft technologies are also under-development  to meet these new needs and, as mentioned before, alternative propellants must be found to overcome the xenon problematic [15, 16]. The reader can refer to the review of Levchenko et al.  for more insights on the future of plasma propulsion.
The working principles of the two most commonly used thrusters (GIT and HT) are described in the following. For more general informations on electric propulsion devices, the reader may refer to detailed textbooks [17, 18] or review articles [19, 8, 3].
Working principles of the two main electric thrusters Gridded Ion Thrusters (GIT)
Whereas GIT can diﬀer by the way the plasma is created (Direct-Current (DC) heating or radio-frequency (RF) heating), the thrust is always produced by ions extracted and accelerated by a multi-aperture grid assembly, commonly two parallel grids, as seen in figure 1.5. Each grid is biased, and the potential diﬀerence induces an electric field which accelerates the ions produced upstream. A source of electrons is then needed to neutralize the ion beam of the plume.
The main characteristic of GIT is that the thrust is directly proportional to the ion current density. However, this extracted ion current density Ji is limited by space-charge saturation. Its maximum value is a function of the discharge voltage V and the distance between the grids d, given by the Child-Langmuir law Jmax = 9 0 r d2 ; (1.8) with 0 the vacuum permittivity, e the elementary charge and M the ion mass. Hence, while GIT have the advantage of a well-focused ion beam (no divergence), this fundamental limitation is a drawback for the maximum achievable thrust.
No more details on GIT will be given here, as these thrusters are out of scope of this manuscript, and the reader can refer to  for the history of GIT or  for an example of GIT global model.
Figure 1.5 – From : Schematic diagram of a DC electron bombardment Gridded Ion-Thruster.
Hall thrusters (HT)
Since their original development in the mid-1960s, mainly in the Soviet Union at the instigation of Morozov, many experimental, numerical and theoretical work have been undergone to better understand the Stationary Plasma Thrusters (SPT), or Hall thrusters (HT), that have now been propelling many satellites.
The Hall thruster working principle will be detailed more extensively com-pared to the GIT. However, for additional insights, the reader can refer to the initial work of Morozov et al. [23, 6] or various topical reviews [24, 25, 22], including one focused on the French research . Other closed-drift thrusters are often associated with HT, such as the Thruster Anode Layer (TAL) with a shorter acceleration region and conducting walls instead of dielectric ones . Even though most of the results that will be presented in this work can be applied to them, the study will be focused on typical SPT thrusters such as the SPT-100 or the PPS®1350.
A side and front view of a Hall thruster are displayed in figure 1.7. Their working principle is rather simple:
1. Neutral gas is injected through an anode at one end of the discharge channel. At the other end, electrons are emitted from a cathode. However, the gas pressure is very low and without an imposed magnetic field, the ionization mean free path is orders of magnitude longer than the thruster size and hence, the ionization eﬃciency (mass flow rate of output ions over the one of input neutrals) is low.
2. For this reason, coils are added to the thruster (often 5 coils, with one central and four on the corners) in order to create a radial magnetic field which acts to confine the electrons: when they enter the channel from the cathode, they are trapped by this magnetic field and due to the E B configuration, they drift in the azimuthal direction. Hence, the ionization eﬃciency will be greatly increased (they will encounter more neutrals to ionize), up to 95%.
3. Due to the potential diﬀerence imposed between the anode and the cathode, an axial electric field will form and accelerate the ions created inside the channel, to produce the thrust. Also, a part of the electrons emitted by the cathode are meant to neutralize this ion plume.
To better understand this working principle, the averaged axial profiles of the main parameters of interest in a HT have been displayed in figure 1.8. First, it highlights a fundamental feature of HT: the presence of two diﬀerent zones, an upstream one (near x = 1.5 cm) where the ions are created, called the ionization zone, and a downstream one (near x = 2.5 cm) characterized by a strong axial electric field, called the acceleration zone . The HT performance depends on how these two zones will overlap: if they are too far, the ions will not experience any electric field and hence will stay immobile after being created. If the overlap is too important, some ions will only experience a small electric field and hence, reach lower axial velocities. We can relate the position of the axial electric field peak to the one of the maximum of radial magnetic field (with Bmax often set to 150-200 G): the electric field is forming where the electron mobility is lower due to the magnetic trapping. In most of the HT, this position corresponds to the exit plane of the discharge channel. Finally, we can see that the ions reach axial velocities of the order of 15 km • s 1.
Despite the apparent simplicity of these thrusters, many questions can be raised during their development phase:
— Which neutral gas should we choose? Xenon is almost always used because of its low ionization energy Ei=12.1 eV and its high molar mass M=131u, but its high and fluctuating cost will probably be a major issue in the next years .
— Which material do we choose for the channel walls? Initally made of alumine Al2O3, dielectric walls of BNSiO2 are now commonly used because they exhibit an excellent combination of fabricability, performance [28, 29] and lifetime. The erosion of the walls should also be taken into account as it aﬀects the thruster performance and lifetime.
— Which magnetic topology should we use? Recently, Magnetically-Shielded (MS) HT have been used to shift the acceleration zone and limit the wall erosion. [30, 31].
— Where to position the cathode? The inherent asymmetry of having an exterior cathode might be an issue and work have been done to use instead a cathode at the center of the discharge channel [32, 33].
All these features aﬀect the performance of HT and should be optimized . Keeping in mind the never-ending human dream of conquering the stars, these many challenges needs to be overcome by studying these complex systems.
Hall thrusters – General concepts
At the heart of this work lies the concept of plasmas. Often referred as the fourth state of matter, plasma is constituted of ions and electrons, and while the bulk plasma is electrically neutral (quasi-neutrality ), a charge diﬀerence can appear at the walls (sheaths) due to the lower mass of electrons (higher velocity) compared to ions. This mass diﬀerence between these two species is the main driver of the collective phenomena in plasmas and is hence responsible for their complexity. Plasmas can be found in diverse systems, from solar wind to tokamaks, by way of Earth’s ionosphere and magnetosphere or neon lamps. They are often classified using their density and temperature: HT plasmas are belonging to the category of « low-temperature plasmas » (LTP). In these weakly ionized plasmas, there is no thermal equilibrium between the electrons and the ions/gas i.e. the electron temperature is often way higher than the ion temperature (which is close to the gas temperature): Te Ti Tg.
To understand the HT plasma physics, some fundamental concepts will be introduced in this section. The reader can refer to various textbooks for extensive details on general plasma physics [35, 36] or plasma instabilities .
Basic concepts of Hall thruster plasma physics
The two important parameters characterizing a plasma are the Debye length d and the electron plasma frequency !pe. The first one corresponds to the typical screening distance between charged particles, i.e. if two particles are separated by a distance lower than d, they can act on each other via Coulomb interaction and a static polarisation can appear. For distances higher than d, this force is screened and the plasma can be considered quasi-neutral.
Hall thruster instabilities
Despite their apparent simplicity, Hall thrusters are very complex systems because of the presence of gradients in many directions. Multiple instabilities can grow and interact together, in a broad range of frequencies. Here, the three main instabilities that will be observed and studied in this work are described: the low-frequency (kHz) Breathing Mode (BM), the intermediate frequency (hundreds of kHz) Ion Transit-Time instability (ITTI) and the high-frequency (MHz) Electron Drift Instability (EDI). A brief description of some other instabilities that will be encountered is also given, but the reader can refer to the review of Choueiri  and the more recent work of Smolyakov et al.  for extensive details.
This perturbation analysis can be applied to the relevant set of equations to obtain a fundamental relation which defines the eigenmodes of the system, the so-called dispersion relation (DR) linking the frequency to the wave vector:
! = !(k). ! is a complex number: its real part corresponds to the frequency !R, while its imaginary part corresponds to the growth rate .
The first hypothesis that have been suggested by Fife  to explain this low-frequency oscillation is a predator-prey model (Lotka-Volterra model ) between the neutrals (prey) and the plasma (predator) , i.e. a competition between the neutrals that are ionized and the ions that cannot exist without the neutrals. More specifically, a BM oscillation can be decomposed in several steps: 1. First, the discharge channel fills up with neutrals from the anode. These neutrals are strongly ionized in the high magnetic field region which depletes their density profile.
2. Then, because there are not enough neutrals in the channel any more, the atom front moves towards the anode where the ionization is less eﬃcient due to a higher electron mobility (lower magnetic field).
3. The discharge current hence decreases, the plasma is extracted by the axial electric field and the neutrals can replenish the channel. The ionization will take place again once the neutral density will be high enough.
The cycle repeats with a period depending on the time needed for the neutral atoms to fill the ionization region. This qualitative description was confirmed by a 1D model of Boeuf and Garrigues  who named this phenomenon, « breathing mode ». Later, Barral and Ahedo  generalized the 0D predator-prey model in 1D.
While this simple explanation gives good order of magnitudes for the oscilla-tion frequency f = pvivn (undamped harmonic oscillator, with vi and vn the ion 2 L and neutral velocity, respectively, and L the length of the acceleration region), two main assumptions are incorrect [50, 51]: the neutral mass flow rate is not taken into account and in 0D, the discretization length for the spatial convection terms is not well-defined (is it the channel length or the ionization zone length?). Hence, this mode is still studied, for example by Hara et al. [52, 51] who high-lighted the importance of electron energy equation to get a positive growth rate through a linear stability analysis; or by Romadanov et al.  who tried to control these oscillations by modulating externally the DC voltage applied on the HT. They also managed to reproduce their experimental results with a simplified 1D model of resistive-ionization mode in quasineutral plasma. This resistive be-haviour was also found to be the origin of this instability by Chable and Rogier : they related it to a Buneman instability, driven by a coupling between the electric field and the ion current.
All in all, the origin of this instability is still debated, with an onset criteria still not well-defined and hence, further investigations are needed.
These oscillations are characterized by a frequency in the range 100-500 kHz which corresponds to the transit time of ions in the acceleration region (around 1 cm): as they commonly travel with a velocity of around 1-5 km/s, their transit-time is on the order of 2-10 s. They were first measured and characterized experimentally by Esipchuk et al. . They described it as a « quarter-wave standing mode » with an antinode near the anode and a node near the exit plane, along with plasma potential variations of the order of 10 to 30 % of Ud. More recently, Vaudolon et al. [56, 57] have also observed this instability through Laser Induced Fluorescence (LIF) measurements of the local Ion Velocity Distribution Function (IVDF). They noticed that the ion exhaust velocity was higher than the maximum theoretical velocity defined by the imposed potential diﬀerence. This so-called wave-riding mechanism can be explained as the following. In 1D, when a particle of mass m and charge q is created with zero velocity at an electric potential 0, it will gain the following energy along its trajectory
Numerically, Fife  was the first to observe this instability in a 2D axial-radial hybrid code, appearing only when the cross-field electron mobility near the anode was low, compared to the one at the exit plane. He described it as 3 successive steps:
1. The axial electric field Ex rises near the anode which implies a local reduc-tion of the density.
2. The lower density further increases Ex by Ohm’s law, which increases Te through Ohmic heating.
3. The disturbance travels downstream led by the sharp electric field, up to the exit where it is mixed with other ions and exhausted.
Later, Hagelaar et al.  and Bareilles et al.  have also observed this ITT instabilities with a 2D axial-radial hybrid code. During a transit-time oscillation, they found that the position of the acceleration zone oscillates, as seen in figure 1.10:
1. The acceleration zone moves upstream: ions that are just accelerating see the electric field « slipping through their fingers » and do not experience the full voltage drop. A low-energy ions group appears.
2. Then, when it moves downstream, a « wave-riding » phenomenon makes the ions accelerate more than the total Ud.
3. It results in an ion beam with two populations: because the acceleration zone moves usually faster upstream, the low-energy ions are a smaller group of energy as low as 100 eV. First, we see that the slow group creates a maximum of ion density, which induces a maximum of potential. Second, when the fast group joins the slow group, the ion current can locally exceed the total current so that the electric field must change sign to maintain current continuity. The ion velocities are reduced by this electric field and the slow group slows down further to only 10 km • s 1 i.e 70 eV.
They state that the perturbation of the plasma density and ion flux cause in turn a perturbation of via the current conservation, providing the back coupling that maintains the transit-time oscillations. They also find that these transit-time instabilities can be reduced by increasing the anomalous electron conductivity.
More details on this instability will be given in chapter 6 in which its presence in our 2D axial-azimuthal simulations will be studied.
Electron Drift Instability (EDI)
The Electron Drift Instability (EDI), sometimes called Electron Cyclotron Drift Instability or E B Electron Drift Instability, was first observed and exten-sively studied theoretically in other fields (such as shockwave propagation across a magnetic field) in the 1970s by Forslund et al. [60, 61], Gary and Sanderson , Lampe et al. [63, 64] and Biskamp and Chodura . They found that this instability is kinetic in nature, with a dispersion relation that can be obtained by linearizing the electron Vlasov equation, coupled with Poisson’s equation and ion cold fluid equations. It results from a coupling between electron Bernstein modes and ion acoustic modes.
In the HT context, the EDI was first found in the axial-azimuthal PIC sim-ulations of Adam et al. , characterized by a short wavelength (in the mm range), a high frequency (in the MHz range) and a dominant propagation in the azimuthal direction. Later, the collective Thomson scattering measurements performed by Tsikata et al. [67, 68, 69] confirmed the presence of EDI in Hall thrusters by measuring the electron density fluctuations in the plume. Since then, EDI has been found in many 1D [70, 71, 72, 73] and 2D [74, 75, 76, 70, 77, 78, 79] PIC simulations and additional experimental diagnostics [80, 81].
Following the work of Ducrocq et al. , Cavalier et al.  have derived the EDI 3D dispersion relation and compared the numerical solutions to the measurements of  with a good agreement. Later, Lafleur et al. [70, 84, 85] have also extensively analysed these instabilities.
If we consider the 2D dispersion relation for this instability, neglecting the radial wave vector kz = 0 (parallel to B), we find that the azimuthal wave vector ky is discrete, with large resonances near the cyclotron frequency (cyclotron resonances) 
Table of contents :
1.1 A new age for spacecraft propulsion
1.1.1 Electric propulsion
1.1.2 A broad range of needs
1.1.3 Working principles of the two main electric thrusters
1.2 Hall thrusters – General concepts
1.2.1 Basic concepts of Hall thruster plasma physics
1.2.2 Hall thruster instabilities
1.2.3 Why do we keep studying Hall thrusters?
1.3 Hall thrusters – Modeling
1.3.1 Multi-fluid models
1.3.2 Particle-In-Cell models
1.3.3 Hybrid and Direct-Kinetic models
1.4 Scope and outline of the thesis
2 Code development: from radial-azimuthal to axial-azimuthal geometry
2.1 A quick overview of LPPic, a 2D PIC-MCC code
2.1.1 Code main structure
2.1.2 A powerful numerical tool
2.2 Preliminary modifications towards axial-azimuthal simulations
2.2.1 Magnetic field
2.2.2 Ionization process
2.2.3 Code optimization
2.3 Neutral dynamics
2.3.2 Model description
2.3.3 Solver verification and validation
2.4 Cathode model
2.4.1 Current Equality condition
2.4.2 Quasi-neutrality condition
2.4.3 Setting positions and energies
2.4.4 Which model should we choose?
2.5 Further improvements of the model
2.5.1 Fake radial dimension
2.5.2 Doubly-charged xenon ions
3 Code verification: 2D axial-azimuthal benchmark
3.1 Why do we need this benchmark?
3.2 Description of the model
3.2.1 Simulation domain
3.2.2 Imposed axial profiles
3.2.3 Boundary conditions
3.3 Code specificities
3.4.1 Main plasma parameters
3.4.2 Azimuthal instabilities
3.5.1 Numerical convergence
3.5.2 Case sensitivity
3.6 A tool to challenge numerical models – Example of the cathode injection
3.6.1 Studied cathode models
3.6.2 Main plasma parameters
3.6.3 Azimuthal instabilities
4 Instability-enhanced electron transport
4.2 Theoretical models for enhanced electron transport
4.2.1 Electron-ion friction force
4.2.2 Theoretical models
4.2.3 Comparison with the nominal simulation case
4.3 2D PIC simulations
4.3.1 Model & Parametric studies
4.3.2 Comparison with parametric PIC results of the friction force derivation accounting for non-Maxwellian electrons
4.3.3 Ion-electron friction force
4.4 Influence of ion-neutral collisions
4.5 Discussion and Conclusion
5 Towards self-consistent simulations of Hall thrusters
5.2 Numerical model
5.2.1 Case description
5.2.2 Importance of neutral dynamics and anode ion recombination
5.2.3 Cathode model
5.3 Axial electron transport during a breathing mode oscillation
5.3.1 Case closer to reality ( = 4)
5.3.2 Influence of vacuum permittivity scaling factor
5.3.3 VDF and plasma fluctuations
5.3.4 Azimuthal instabilities
6 Instabilities in Hall thrusters
6.1 Ion transit-time instabilities
6.1.1 Current understanding
6.1.2 Analysis of ITT: « 1D axial » simulations
6.1.3 Origin of the ITTI?
6.1.4 Evidence of ITTI in the self-consistent axial-azimuthal simulations
6.2 Influence of the azimuthal domain length Ly
6.2.1 Axial electron transport
6.2.2 Azimuthal instabilities
6.3 Interaction between ITTI & EDI
6.3.1 Preliminary analysis on the self-consistent simulations
6.3.2 Insights from the simplified benchmark case
7.1 Summary of the thesis
7.1.1 Importance of code verification
7.1.2 Origin of electron anomalous transport
7.2 Recommendations and prospects
7.2.1 Improvements of LPPic
7.2.2 Code validation with experiments
7.2.3 Towards an engineering tool?
A Sod test case