J Porto
email: jcportoh@gmail.com

P Q Elias

A Ciardi

Anisotropic Electron Heating in an Electron Cyclotron Resonance Thruster with Magnetic Nozzle

In a grid-less Electron Cyclotron Resonance (ECR) plasma thruster with a diverging magnetic nozzle, the magnitude of the ambipolar field accelerating the positive ions depends of the perpendicular energy gained by the electrons. This work investigates the heating of the electrons by electromagnetic waves, taking their bouncing motion into account in a confining well formed by the magnetic mirror force and the electrostatic potential of the thruster. An electromagnetic Particle-In-Cell (PIC) code is used to simulate the plasma in a magnetic field tube. The code's Maxwell solver is based on a semi-Lagrangian scheme known as the Constrained Interpolation Profile (CIP) which enables larger time steps. The results show that anisotropic plasma heating takes place exclusively inside the coaxial chamber, along a Dopplerbroadened zone. It is also shown that a trapped population of electrons with a larger perpendicular energy exists in the plume.

I. INTRODUCTION

Electric thrusters play a fundamental role in the field of space propulsion. Their main advantage lies in an efficient use of the propellant mass, and therefore a reduced consumption of propellant. Hall Effect Thrusters or Gridded Ion Engines are examples of the most well-known and flight-proven technologies in the current propulsion market nowadays. Both technologies eject an ion beam which is subsequently neutralized to prevent the spacecraft from charging. Several components of these technologies, such as the acceleration grid or the neutralizer, are subject to erosion and wear and for this reason, meeting the challenging lifetime targets requires careful optimization and demanding testing [START_REF] Holste | Ion thrusters for electric propulsion: Scientific issues developing a niche technology into a game changer[END_REF] . The complexity of some of the components has driven the investigation of alternative concepts of propulsion devices that require neither a grid nor a neutralizer. The Electron Cyclotron Resonance (ECR) plasma thruster [START_REF] Jarrige | Characterization of a coaxial ECR plasma thruster[END_REF][START_REF] Vialis | Geometry optimization and effect of gas propellant in an electron cyclotron resonance plasma thruster[END_REF] is one of these concepts and it is the subject of the present study.

The ECR plasma thruster consists of a semi-open chamber where a quasi-neutral plasma is heated by electron cyclotron resonant microwaves at 2.45 GHz, and accelerated by a magnetic nozzle. This concept was first proposed in the 1960s in the works of Miller et al. [START_REF] Miller | Cyclotron resonance propulsion system[END_REF] and Nagatomo [START_REF] Nagatomo | Plasma Acceleration by High Frequency Electromagnetic Wave in Static Magnetic Field Gradient AlAA Electric Pronulsion and Plasmadynamics Conference[END_REF] , then further developed by Sercel [START_REF] Sercel | An experimental and theoretical study of the ECR plasma engine[END_REF] . These studies used a prototype with a wave-guide structure to couple the microwaves to the plasma. Their results showed that it was possible to achieve specific impulses and thrust values high enough to be of interest for space propulsion applications [START_REF] Sercel | An experimental and theoretical study of the ECR plasma engine[END_REF] . Nonetheless, the inefficiency, size and weight of the micro-wave sources and electromagnets at that time led to a stagnation of the research on ECR thrusters for several years. Interest for this technology arose again recently with experimental works [START_REF] Jarrige | Characterization of a coaxial ECR plasma thruster[END_REF][START_REF] Jarrige | Performance Comparison of an ECR Plasma Thruster using Argon and Xenon as Propellant Gas[END_REF] . In particular, the use of coaxial microwave coupling structures and compact rare-earth permanent magnets were instrumental in designing compact sources (a schematic of the design is shown in Fig. 1a).

More experimental and theoretical efforts has since been made in order to get a deeper understanding of the physical phenomena governing the plasma heating and acceleration in the thruster. Experimental characterizations of the plasma properties have been carried out using different measurement techniques such as Langmuir and Faraday probes, Laser Induced Fluorescence diagnostics, diamagnetic loops and thrust balances [START_REF] Jarrige | Characterization of a coaxial ECR plasma thruster[END_REF][START_REF] Hepner | Low Frequency Wave Detection in the Plume of a Low Temperature Magnetic Nozzle[END_REF][START_REF] Correyero | Plasma beam characterization along the magnetic nozzle of an ECR thruster[END_REF][START_REF] Correyero | Characterization of diamagnetism inside an ECR thruster with a diamagnetic loop[END_REF][START_REF] Vialis | Direct Thrust Measurement of an Electron Cyclotron Resonance Plasma Thruster[END_REF] . Unfortunately, most of the experimental studies so far have been limited to survey the plasma outside the thruster coaxial chamber. Recently, a resonant probe was developed to measure an electron density of about 1 × 10 11 cm -3 at the source exit plane, close to the coaxial chamber [START_REF] Boni | The curling probe: A numerical and experimental study. Application to the electron density measurements in an ECR plasma thruster[END_REF] . In the source, it is likely that the plasma density is higher (∼ 1 × 10 12 cm -3 ) with electron temperatures of a few tens of eV .

From a theoretical point of view, as a first step, global models describing the energy balance in the plasma source were proposed as a means to obtain the key parameters of the thruster [START_REF] Cannat | Optimization of a coaxial electron cyclotron resonance plasma thruster with an analytical model[END_REF][START_REF] Verma | Thrust evaluation of compact ECR plasma source using 2-zone global model and plasma measurements[END_REF] . While this approach yielded good agreement with measured electron temperature at high mass flow rate or high pressure, they failed at the lower mass-flow rate where the thruster achieves its best performance. Indeed, the assumptions of uniform electron temperature and isotropic Maxwellian electron distribution are too crude approximations when collisionality decreases and the electron mean free path becomes much larger than the source length: in that range non-local effects become prevalent, as electrons undergo a bouncing motion along the magnetic field line. Those electrons which cross the ECR surface can gain energy depending on their phase in the gyromotion [START_REF] Lieberman | Theory of electron cyclotron resonance heating. II. Long time and stochastic effects[END_REF] , which leads to a strong anisotropy of the distribution function. An attempt to account for this stochastic heating in the plasma was made by considering the electron heating as a random walk in phase space [START_REF] Peterschmitt | A diffusion model in velocity space to describe the electron dynamics in an ECR plasma thruster with magnetic nozzle[END_REF] .

2 While this model provided a qualitative agreement with the measured ion energies, it could not account for the plasma feedback on the waves (assumed constant) and the collisions along the bouncing motion. Recently Sánchez-Villar et al. [START_REF] Sánchez-Villar | Coupled plasma transport and electromagnetic wave simulation of an ecr thruster[END_REF] performed 2D axisymmetric simulations of the thruster with a hybrid model consisting of particle-in-cell (PIC) ions and a fluid model for the mass-less electrons. One of the main findings of this study was the identification of different regions in the source where the waves are either propagating or evanescent, with most of the power absorption taking place close to the inner conductor, near the ECR surface. By acting as a sink for the plasma, the inner conductor induces a decrease of the plasma density in its vicinity, enabling the propagation of electromagnetic waves downstream of the ECR surface. These features lead to the formation of a hot electron beam close to the inner conductor, with a colder plasma in the bulk of the source. While these 2D results provided important insights on the operation of the thruster, some assumptions of the fluid model limit the validity of the results obtained from these simulations. The most important one being the assumption of isotropic electron temperature which excludes anisotropic heating in the directions parallel and perpendicular to the magnetostatic field. This latter point is still an open question for this technology. Indeed, ECR heating is expected to lead to anisotropic heating of the electron translation modes. This difference affects the power losses near the source walls and the potential drop in the magnetic nozzle [START_REF] Hooper | Plasma flow resulting from electron cyclotron resonance heating on a magnetic hill[END_REF] . However, most of the electron temperature measurements performed in the thruster plume did not differentiate perpendicular and parallel electron temperature (with respect to the local magnetic field direction). A way to measure the electron temperature anisotropy is, for example, incoherent Thomson scattering [START_REF] Vincent | A compact new incoherent Thomson scattering diagnostic for low-temperature plasma studies[END_REF] , but this type of measurement is not presently available in the ECR source. At any rate, this heating is intimately linked to the absorption of the electromagnetic waves in the coaxial source.

Another issue is the non-local transport due to the bouncing magnetized electrons in the nozzle (the electron mean free path is greater than the source radius). In particular, the production and the heating of the electrons are not necessarily at the same location.

While gaining a better understanding of these issues should firstly rely on experimental measurements, the challenges associated with such an investigation are a strong incentive to use numerical models, even if simplified, to investigate the main physical processes at play in ECR thrusters. In particular, such a model should be able to account for the selfconsistent wave absorption and the anistropic heating, as well as the non-local effects and bouncing motion of the particles. Electromagnetic kinetic models, such as Particle-In-Cell (PIC) or Vlasov methods, are natural candidates for this task.

There are currently a few works using kinetic simulations of propulsion devices exploiting the ECR phenomenon, however the majority of these developments are concerned with gridded ion thrusters with ECR heating [START_REF] Nishiyama | 20-cm ECR plasma generator for xenon ion propulsion[END_REF][START_REF] Shen | Development of a Micro ECR Ion Thruster for Space Propulsion[END_REF][START_REF] Takao | Threedimensional particle-in-cell simulation of a miniature plasma source for a microwave discharge ion thruster[END_REF][START_REF] Yamashita | Numerical investigation of plasma properties for the microwave discharge ion thruster µ10 using PIC-MCC simulation[END_REF] , where the plasma acceleration is achieved by a grid-imposed electric field and not the plasma expansion in the magnetic nozzle, as in our design. Takao et al. [START_REF] Takao | Threedimensional particle-in-cell simulation of a miniature plasma source for a microwave discharge ion thruster[END_REF] successfully modelled a gridded ion thruster where the ions are produced in an ECR source at 4.2 GHz. The authors used a Particle-In-Cell (PIC) code considering the microwave electric field as a temporal modulation of its initial amplitude obtained by simulating the microwave propagation without plasma. Therefore, in this approach the plasma feedback on the wave was considered negligible.

The main purpose of the our study is to perform full-PIC electromagnetic simulations of the plasma in the thruster, taking into account the plasma feedback on the wave propagation, and to investigate the heating and confinement of the electrons. For this purpose it is necessary to simulate the microwave propagation and its interaction with the charged particles in the source and the nozzle region. However, due to the complexity and computational cost of simulating a full 3D configuration (which should include the nozzle region), we restrict our investigation to the simplified case of an isolated magnetic flux tube. This approximation effectively restricts the phase space to 4 dimensions (1 dimension in space, 3 dimensions in velocity space), and a 1D3V electromagnetic Particle-In-Cell can be used to model the ECR thruster.

We show that the electron heating takes place over a broad region in the thruster source and leads to a significant anisotropy (ratio T e⊥ /T e ∈ [2.5, 7.5]). The perpendicular electron temperature reaches a first maximum in the source and, surprisingly, has a second maximum in the downstream region. The explanation for these features lies in the confinement of electrons in the potential well formed by the combination of the diverging magnetic field and the electrostatic potential.

II. NUMERICAL MODEL

A. The Quasi-One-Dimensional Approach

In the coaxial ECR thruster, an axially magnetized cylindrical permanent magnet creates a diverging static magnetic field B MS in the source and in the plume region [START_REF] Vialis | Direct Thrust Measurement of an Electron Cyclotron Resonance Plasma Thruster[END_REF][START_REF] Cannat | Optimization of a coaxial electron cyclotron resonance plasma thruster with an analytical model[END_REF] . This shape for the magnetostatic field was chosen to ensure a magnetic confinement at the close end of the coaxial chamber (called backplate in Fig. 1) while allowing the electrons to get accelerated in the plume thanks to the divergence of the magnetic field lines.

In fact, the ECR uses a diverging magnetic field whose magnitude decreases from approximately 100 mT at the back of the source to around 5 mT 10 cm downstream of the thruster exit plane. Under these conditions, assuming an electron temperature around T e ≃ 10 eV, the Larmor radius of the electrons is between r L ≃ 0.07 mm -1.4 mm. Thus electrons are strongly magnetized in the source and in the near-field plume region, while ions remain mostly unmagnetized. As a consequence, before the onset of plasma detachment, electrons and ions are bound to the magnetic field tube. Several mechanisms may account for the plasma detachment [START_REF] Ahedo | On plasma detachment in propulsive magnetic nozzles[END_REF] : collisions, stretching of the magnetic field lines, electron demagnetization and plume instabilities. While it is out of the scope of this work to study the dominant mechanisms, several recent works have investigated some of these effect in 2D PIC simulations [START_REF] Andrews | Fully kinetic model of plasma expansion in a magnetic nozzle[END_REF][START_REF] Hu | Kinetic insights into thrust generation and electron transport in a magnetic nozzle[END_REF][START_REF] Emoto | Numerical investigation of internal plasma currents in a magnetic nozzle[END_REF] .In This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

this work, we decided to rely on experimental evidence to define the section of the nozzle where the plasma remains bound to the field lines. Recently, Little and Choueiri [START_REF] Little | Electron Demagnetization in a Magnetically Expanding Plasma[END_REF] have mapped the plasma potential in a magnetic nozzle to show that a good criterion for detachment is χ p = r L /L ∇B ≃ 0.1, where r L is the electron Larmor radius and L ∇B = (∇B/B) -1 is the characteristic length scale of the magnetic field gradient. In the region of the nozzle where χ p = r L /L ∇B < 0.1, the plasma remains attached to the magnetic field. In our case, we considered a magnetic field with L ∇B ≃ 5 -10 cm. Under this condition, we have χ p < 0.1 up to L = 10 cm downstream of the nozzle, and it is a reasonable assumption to consider that electrons do not detach from the magnetic field tube over this distance.

As a consequence of this assumption, we decided to consider the creation and formation of the plasma enclosed in a magnetic field tube of length L = 10 cm. More precisely, a portion of the thruster chamber and plume was represented by a quasi-1D model of a magnetic field tube with a varying cross-sectional area, as seen in Fig. 1b. There are several examples of the use of quasi-1D models in the space propulsion field. Niewood and Martinez-Sanchez 30 used it to model a Magnetoplasmadynamic thruster, while De Giorgi and Fontanarosa [START_REF] De Giorgi | A novel quasi-one-dimensional model for performance estimation of a Vaporizing Liquid Microthruster[END_REF] studied a Vaporizing Liquid Microthruster with this approach. Recently, Saini and Ganesh 32 also used this approach to model plasma expansion in a Radio-Frequency thruster. The moderate computational cost of a 1D3V model of the thruster facilitates the analysis of the plasma behavior in both the coaxial chamber and in the magnetic nozzle, and importantly, taking into account the nozzle is critical to resolve the bouncing motion of the electrons.

The quasi-1D model assumes that the electrons and the ions are confined within a diverging magnetic flux tube, whose area is related to the axial magnetic field intensity through the conservation of the magnetic flux:

A(x)B x (x) = A 0 B 0 (1)
The model further assumes that the electromagnetic fields and all the plasma properties are constant across the section of the flux tube. For the ECR thruster under consideration [START_REF] Vialis | Direct Thrust Measurement of an Electron Cyclotron Resonance Plasma Thruster[END_REF][START_REF] Cannat | Optimization of a coaxial electron cyclotron resonance plasma thruster with an analytical model[END_REF] , the shape of the magnetic field lines close to the antenna is well approximated by an exponential function. For the sake of simplicity we approximated the magnetic field as:

B x (x) = B 0 exp - x L B (2)
In addition, we considered cylindrical symmetry for the static magnetic field around the field tube centerline. These assumptions mean that the particles guiding centers remain on the centerline. Since the plasma is assumed uniform in the cross section, this approach does not allow the formation of a diamagnetic current and E × B drifts.

From now on, the term parallel and the subscript will refer to the direction parallel to the magnetostatic field lines. Similarly, perpendicular and the subscript ⊥ refer to the direction perpendicular to the magnetostatic field lines. The source region, which corresponds to the coaxial cavity in Fig. 1b, was defined by 0 ≤ x ≤ L S , where L S is the coaxial source length. The plume region, which corresponds to the plasma expansion in vacuum, was defined by x ≥ L S .

B. Particle-In-Cell Code overview

The simulations were carried out with the Particle-In-Cell (PIC) code Rhei, which was developed to simulate low pressure cold plasmas and is adapted to parallel architectures. It can be run with either a pure MPI or a hybrid MPI/OpenMp parallelization. The code integrates a Monte-Carlo Collision (MCC) module to simulate the collisions between the charged particles and a prescribed neutral background. At each iteration, once the electrostatic and the electromagnetic fields were computed, the position of each macro particle labeled "p" was updated using dx p /dt = v p , and the velocity using Eq. 3. Each macro-particle represents W physical particles. The value of W used in the simulation is given in table I.

m s dv p dt = q s E ESp + E EMp + v p × B MSp + B EMp (3)
In Eq. 3, q s is the charge of the particle, m s the mass, x p the position, and v p the velocity. Regarding the fields, they were computed at the location of the particle p using linear interpolation function, where E ESp is electrostatic field from the charge distribution, B MSp is magnetostatic field from the permanent magnets and E EMp and B EMp : are electromagnetic fields produced by the microwave source and by the plasma itself.

The equations of motion were integrated using the leap-frog method and the Boris scheme to get the v × B rotation from the Lorentz force 33 . Details of the integration in the context of the quasi-1D model are provided in appendix A. Particle quantities were projected on a uniform grid using linear shape functions.

The Rhei code development follows a test-driven approach to ensure the robustness and the maintainability of the code over time. Additionally, several test cases were run as a validation of the code. The first elementary test was the simulation of a magnetic bottle. The simulation domain, with a converging-diverging parabolic magnetic field, was uniformly loaded with a Maxwellian electron population. At the end of the simulation the electron distribution in velocity space v , v ⊥ was plotted to verify that the loss cone angle is coherent with the expected theoretical value arcsin B 0 /B Max . The second elementary test concerned the electromagnetic modes in a one-dimensional magnetized plasma. The simulation domain was initialized with a uniform Maxwellian distribution of electrons and cold ions. The random fluctuations excited the modes of the plasma. The resulting dispersion curves were obtained by computing the discrete 2D Fourier transform of the electric fields during the simulation. This was compared to the expected theoretical description of the extraordinary and the ordinary wave. Finally, the third test case was the classical capacitively coupled discharge in Helium, which verified in particular the collision module [START_REF] Turner | Simulation benchmarks for low-pressure plasmas: Capacitive discharges[END_REF] .

Collisions

The Monte-Carlo Collision module used the Null Collision technique [START_REF] Skullerud | The stochastic computer simulation of ion motion in a gas subjected to a constant electric field[END_REF] to speed-up the computation of the collisions by removing the velocity dependency of the total collision crosssection. Assuming N p collision processes defined by their respective cross sections σ i (v), i = 1..N p , a null collision crosssection is defined as σ 0 (v) such that:

σ 0 (v) = max v≥0 Np ∑ i=1 σ i (v) - Np ∑ i=1 σ i (v) (4)
A first test over all the particles of species s found the fraction of particles which undergo a collision with the background. In that case the total cross section σ T = ∑ i=0 N p σ i (v) (including the null collision process) did not depend on the velocity (thus avoiding a costly interpolation to get the cross section for all the particles). Then a second test among those selected particles computed all the collision cross sections for their given relative velocity and determined which cross section to use (including the null collision). When this test pointed to the null-collision cross section, then the particle did not experience an actual collision and was left unchanged. When the test pointed to another cross section, the the particles experienced a collision.

For the collisions of electrons with Xenon neutrals, we considered a simplified set of three processes: elastic, ionization, and excitation. Excitation processes were lumped into a single process. Electron impact ionization and excitation were taken from the Morgan (Kinema Research & Software) database, while the total elastic scattering is from Ref. 36. For all electronic processes, we assumed an isotropic scattering of the relative velocity vector between the electron and the target during the collision. For the ionization collisions, the kinetic energy of the projectile electron was equally split (after subtracting the threshold energy) between the secondaries. For the collisions of Xenon ions with Xenon neutrals, we considered isotropic scattering and backscattering [START_REF] Turner | Simulation benchmarks for low-pressure plasmas: Capacitive discharges[END_REF] . Ion cross sections comes from Ref. 37. All electronic and ionic processes conserved momentum and total energy (kinetic plus internal). In order to start with a simplified description simulating weakly ionized plasmas, in which the collisions with the neutral particles are the dominant process, Coulomb collisions were not considered in the code. Indeed, the electron-ion collision frequency ν ei , for Maxwellian electrons, is given by:

ν ei = ω p Λ ei ln Λ ei (5)
Here, ω p is the plasma frequency and ln Λ is the Coulomb logarithm. For the typical simulation conditions in the thruster source, as it will be shown below, the maximum plasma density was n e ∼ 1 × 10 11 cm -3 , the electron temperature was T e ∼ 10 eV and the electron-neutral elastic collision frequency was ν en ∼ 1 × 10 7 s -1 . This gave ln Λ ∼ 12 -15, ω p ∼ 1.8 × 10 10 rad s -1 . Consequently, the maximum electron-ion collision frequency was ν ei ∼ 1 × 10 5 s -1 , much less than the the electron-neutral collision frequency ν en .

The neutral gas in the thruster is injected at the backplate (see Fig. 1a) and expands in the source resulting in a decreasing density. Since there is no measurement of the neutral gas density profile in the thruster, and to avoid a costly particle simulation of the neutral particles, we modelled this expansion heuristically by assuming that the neutral background density followed an exponential profile:

n n (x) = n n 0 exp - x L n (6)
where n n 0 is the maximum density of neutrals found at the close end of the source, and L n is the neutrals density characteristic length. The assumption of a time-independent neutral gas density profile means that the simulation did not conserve the total mass, momentum and energy. In addition, it means that the neutral gas depletion due to ionizing collisions was not considered. However, both of those limitations are acceptable in the frame of this work which does not seek to compute the total thrust and energy balance but is concerned with the particle heating and trapping. To estimate the neutral depletion we note that the ion removal is driven by their velocity (at most 10 km s -1 ), while the neutral removal is driven by their thermal speed, ∼ 200 m s -1 . From mass balance, the neutral inflow is balanced by the ion flux and the neutral outflow. Using the characteristic speeds and the typical parameters for the gas density n g ∼ 1 × 10 14 cm -3 , and plasma density n max e ∼ 3 × 10 11 cm -3 , gives a neutral depletion of at most 10%, indicating that the assumption of a static background remained consistent with the assumed density profile.

Fields solvers

The Rhei code solved the Poisson equation to compute the electrostatic space potential Φ and the electric field (E ES = E x x) using Eq. 7 where the charge density is ρ s . The solver implements a second order finite difference discretization and the resulting linear system is inverted using an iterative method (GMRES) [START_REF] Stoer | Iterative Methods for the Solution of Large Systems of Linear Equations. Additional Methods[END_REF] .

∇ 2 Φ(x,t) = - ρ s (x,t) ε (7)
In addition, an electromagnetic solver computed the fields produced by the microwave source and by the plasma itself: E EM = E y y + E z z and B EM = B y y + B z z. This solver was based on the Constrained Interpolation Profile (CIP) method explained in detail in Ref. 39. This method considers not only the electromagnetic fields but also their spatial derivatives, therefore suppressing instabilities and providing lower numerical dispersion even when using coarse grids and large time steps [START_REF] Kajita | Computation of lightning electromagnetic pulses using the constrained interpolation profile method[END_REF] . The use of this method is a novel solution for a PIC code since most of the electromagnetic solvers are based on conventional approaches like the finite-difference timedomain method (FDTD). It was shown that it provides higher accuracy than the latter under the condition of identical cell size [START_REF] Okubo | Analysis of an electromagnetic field created by line current using constrained interpolation profile method[END_REF] .

The CIP method is a semi-Lagrangian scheme that circumvents the Courant-Friedrichs-Lewy (CFL) stability This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834 condition [START_REF] Yabe | The Constrained Interpolation Profile Method for Multiphase Analysis[END_REF][START_REF] Smolarkiewicz | A Class of Semi-Lagrangian Approximations for Fluids[END_REF] , i.e., (u∆t/∆x) < 1 where u is the magnitude of the velocity, ∆t is the time-step, and ∆x the length interval. This feature allows computations with CFL values ≥ 1.0, as can be seen in Ref. 44 and 45 where the authors performed simulations using a CFL value of 2.6 in a Cartesian coordinate system. The gain in computational time, that is afforded by using high CFL values, is a key factor that enables the self-consistent kinetic simulations presented here to reach steady-state. In this paper, CFL values close to 3 were used for the simulations. As a check, simulations were also run with CFL=0.6 and compared to the results obtained with larger time steps. The results were identical to the one at larger time-steps, within small variations due to the noise inherent to the statistical nature of the PIC simulations. Finally, the CIP scheme does not necessarily maintain the divergence-free condition for the dynamic field B EM . However, B EM is smaller than the magnetostatic field (which has divergence equal to zero by construction, see appendix A) by several orders of magnitude, over the whole computational domain. Therefore, the resulting error on the total divergence was considered to be negligible.

Boundary conditions

As it was shown in Fig. 1b when describing the model, the domain goes from x = 0 at the left side which corresponds to the backplate and the microwave input, to the right-end at x = L, as discussed in II A.

Electrostatic: At the right end of the computational domain x = L, we imposed a Dirichlet boundary condition, with Φ(L) = 0, to simulate the presence of a grounded vacuum chamber wall. The dielectric backplate, at x = 0, is in contact with the plasma and therefore its surface voltage Φ BP is changed by the collection of charged particles. This can be modeled as a capacitor. Hence, the evolution of Φ BP is given by ∆Φ BP = ∆Q/(C∆t), where ∆Q is the charge deposited at the backplate at each time step, and C is an equivalent capacitance under the assumption that the backplate is in contact with a grounded conductor. This capacitance is computed by considering that the backplate is a plane capacitor, its value is of a few picoFarads. Changing its magnitude modifies the charging rate of the backplate and thus the transient phase of the computation. However, its does not affect the steady-state voltage of the backplate. This approach guarantees that at steady-state, the ion flux equates the electron flux on the backplate. In principle the steady-state value of the backplate potential is also affected by other processes such as secondary electron emission or charge migration. However, for this study, these processes were neglected.

Electromagnetic: In the coaxial ECR thruster, the microwaves are injected as Transverse Electro-Magnetic (TEM) mode. For this 1D simulation, the TEM mode can be seen as a linearly polarized wave, where the radial component of the electric field is along the transverse y axis, the azimuthal magnetic field defines the z axis and the wavevector direction is along the longitudinal x axis. Therefore, the microwaves were injected at the backplate as a propagating wave with a linear polarization along the y-axis. The incident wave was parametrized by its power per unit area P in and its frequency f EM = ω/2π. The electric fields from the injected linearly polarized wave were computed as E y = √ µcP in sin (ωt) and

E z = 0.
The injected microwave input power per unit area P in could be fixed, or it could be adapted to keep a roughly constant predefined number of particles N target during the transient phase. This feature was intended to speed up the simulations by reproducing a faster plasma response to a given variation in the simulation's parameters. The value of P in can be regulated with an attenuation factor α ≤ 1 varying with the number of particles in the domain: α = exp(-N particles /N target ). A run performed without this regulation confirmed that it did not have an effect on the final steady state but only on the duration of the transient phase.

Particles: We imposed a loss condition at both ends of the domain, for both ions and electrons. Particles crossing these boundaries are suppressed from the simulation. As a simplifying assumption, secondary emission processes on the backplate were not considered in this first approach.

Cross field diffusion loss model

Electron cross-field diffusion is an important mechanism to model to get a more accurate representation of the discharge loss mechanisms. Previous works using PIC codes for electric thrusters took it into account as wall losses by artificially increasing the collision rate or by using a profile of the cross field diffusion based on empirical evidence [START_REF] Blateau | PIC simulation of a ceramic-lined Hall-effect thruster[END_REF][START_REF] Fox | Advances in Fully-Kinetic PIC Simulations of a Near-Vacuum Hall Thruster and Other Plasma Systems[END_REF] . The electron balance equation is:

∂ n e (r,t) ∂t + ∇ ⊥ • n e u ⊥ + ∇ • n e u = k ion n e (r,t) (8)
Where u ⊥ and u are the electron macroscopic velocity perpendicular and parallel to the local magnetic field, respectively, and k ion is the ionization rate.

For our 1D3V simulations, the transport along the magnetic field is taken into account by the kinetic model. However, the perpendicular transport cannot be modeled with a 1D model. Therefore we simulated the particle losses into the coaxial chamber walls using a phenomenological, Monte Carlo loss model, as shown in Fig. 1b. The probability of an electron impacting the walls of the coaxial chamber was calculated from the diffusion equation of electrons across the magnetic field based on the assumption that their number density profile in the radial direction was independent of time and axial position. In a cylindrical coordinate system it can be expressed as the product n e (x, r,t) = n e 0 (x,t)g(r). The balance equation (for a constant diffusion coefficient D) integrated over the radius of the flux tube r max was then given by:

∂ n e 0 (x,t) ∂t + ∂ n e 0 (x,t)u x (x,t) ∂ x = -ν L n e 0 (x,t) + k ion n e 0 (x,t) (9)
With the loss frequency given by:

ν L = -r max g ′ (r max ) S D (10)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

Where r max is the radius of the flux tube, and the weighted cross section S is given by:

S = rmax 0 rg(r)dr and (11)
A first choice for the diffusion coefficient D would be a coefficient based on classical diffusion obtained from theories on standard electron-neutral collisions. It can be seen in Eq. 12 where τ = 1/ν is the collision period with the neutral background. However, the electron mobility tends to be higher than the value predicted by this classical diffusion approach [START_REF] Chen | Introduction to Plasma Physics and Controlled Fusion[END_REF] . The cause of this discrepancy is an active area of research in the electric propulsion field [START_REF] Garrigues | Empirical electron cross-field mobility in a Hall effect thruster[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] . As a consequence, we decided to use the Bohm coefficient, which is a phenomenological coefficient accounting for the anomalous cross-field diffusion.

D Bohm = 1 16 k B T e eB or D classical = ω c τ 1 + (ω c τ) 2 k B T e eB (12)
The probability for a given particle to be lost between t and t + ∆t is given by p L = ν L ∆t. In our quasi-1D model, the flux conservation relates the magnetostatic field to the crosssectional area of the magnetic field tube as shown in Eq. 1.

In this work we assumed g(r) = J 0 (k 0 r/r max ), with k 0 the first zero of the Bessel function. Then Eq. 10 has the form ν L ∝ k 2 0 /r 2 max D. Since D ∝ B -1 and in the model BS = Bπr 2 max is a constant, the loss probability does not depend on the position along the flux tube and is given by:

p L (x) = 2 3 π 16 k 2 0 1 2 m v(x) 2 dt eA 0 B 0 ( 13
)
Where ∆t is the time step, m v(x) 2 /2 is the electron's mean kinetic energy, and A 0 = A(0) is the cross-section of the magnetic field tube at x = 0. The losses are computed at each time-step. For all electrons in the source (such as x ≤ L S , L S being the length of the coaxial chamber as shown in Fig. 1b), the probability p L is computed using equation 13. A random number x is drawn from a uniform distribution. If x ≤ p L , the electron and a neighboring ion are removed from the simulation.

C. Simulation Setup

The electron dynamics and the electromagnetic solver were updated every iteration. For these one-dimensional calculations the real Xenon mass for the ions was used and to speed up the calculations a subcycling was used so the ion's position and velocity that were updated every 10 time steps as given by ∆t ions in Table I. The collisions were also computed every 10 time steps as given by ∆t coll . The charged particle's population was seeded using a uniform density distribution (N ∼ 10 3 ). The electron's initial energy along each of the x-y-z axis was set to T e = 20 eV, while ions were assumed cold T i = 0.03 eV. These values were intended to reproduce a non-equilibrium plasma at low density. The choice of the initial electron temperature T e = 20 eV is somewhat arbitrary. Checks run with several energy values between 10 eV and 30 eV showed no impact of the initial electron energy on the final characteristic of the steady state. To sustain the plasma at the beginning, a plasma source located at 2 mm from the backplate injected electrons at 3 × 10 5 m s -1 and ions at 3 × 10 2 m s -1 during the first 150 ns of the simulation. These velocities were specified along each of the x-y-z axis. Here the idea was to sustain the initial plasma long enough for the ionization to pick up and the plasma density to grow. The conditions for the simulation presented below are shown in table I. With this choice of magnetic field profile, the resonance condition f EM = eB/2πm was met at x = 6.7 mm. The simulation was run until it reached a steady state, usually after around 30 µs which represents between 5 to 8 ion transit times. The definition of this steady state was done by following up the variation of the total number of particles in the domain, its mean kinetic energy, and the particle flux at the backplate and the plume since an equal number of ions and electrons must be impacting both surfaces, as shown in Fig. 2. At the end of the simulation, when the steady state was reached, the plasma properties were obtained by calculating the time average for each parameter over several time steps. Overall, the wall time of the simulation was 44 hours, with 12 OpenMP threads.

III. RESULTS

Figure 3 shows the steady-state plasma potential distribution over the whole computational domain. Except for slight random fluctuations on the instantaneous potential, no large scale fluctuations were observed. Time-averaging improved the signal to noise ratio but did not blur the shape of the profile. The backplate reached a positive steady state potential of around 70 V. The peak of the plasma potential was 105 V, and it was reached at around 3 mm, interestingly not at the ECR surface (indicated with a vertical dashed line). Indeed, the shapes of the plasma density and potential are driven by the ionization rate. For this simulation, the ionization rate This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

was monotonically decreasing, because the background density decrease was faster than the ionization rate increase due to the plasma heating. As a consequence, the maximum ionization was upstream of the ECR surface. This peak in the plasma potential formed a barrier. As a result, ions collected on the backplate were necessarily created in a region where x ≤ 3 mm, while ions collected downstream were created in a region where x ≥ 3 mm and accelerated into the nozzle by the ambipolar electric field. Sheaths were formed at the backplate and the vacuum chamber wall. The plasma sheath width on the backplate was ∼ 0.1 mm. At the downstream end x = L of the domain, as shown by Fig. 4, the plasma sheath began at around 90 mm. This size was consistent with a Debye length λ D ∼ 1 -5 mm for a plasma density around 1 × 10 8 cm -3 . The electron and ion peak number density was 1.12 × 10 11 cm -3 at x = 1.5 mm. Recall that the ECR condition is met at x = 6.7 mm. In Fig. 5 we plotted the electron's mean kinetic energy in both the axial (e ) and the perpendicular (e ⊥ ) direction as a function of the axial position on the domain. First, we observed that the mean parallel kinetic energy remained nearly constant, around 4 -5 eV, over the whole simulation domain. The perpendicular energy was higher than the parallel component, which underlined the anisotropic heating of the electrons in this thruster. More precisely, the mean perpendicular kinetic energy e ⊥ reached a first peak at around x = 9 mm and then decreased before reaching a global maximum of 25 eV at x = 45 mm. After this point, e ⊥ decreased until the end of the simulation domain. Over the whole simulation domain, the anisotropy ratio T e,⊥ /T e, was found to vary between 2.5 and 7.5. Given that the ECR heating increases the perpendicular energy of the electrons, it was expected to see an anisotropic behavior depending on the direction parallel or perpendicular to the magnetic field lines. However, the second broad energy peak in the downstream part of the magnetic nozzle was puzzling. To get a better understanding of these feature, it was necessary to evaluate the energy deposition by the electromagnetic field.

A. Electromagnetic Energy deposition in the source

To understand how the field energy was transferred to the particles, we considered the energy balance equation, including the Poynting flux (its derivation is provided in appendix B).

∂ ε EM + ε ∂t + ∇ • (Q + Π) = S coll (14)
In this equation, ε and ε EM stands for the electron kinetic energy density and the electromagnetic energy density, respectively. Q and Π are the kinetic energy flux and electromagnetic energy flux; S coll , whose expression is given in Eq. B4, is the volume power loss term due to the collisions and the diffusion. This latter term account for the energy lost by elastic and inelastic collisions with the neutral background and by the particles removed by the loss model. Since we were interested in the steady state regime, and considering that the field quantities depend on x only, this was further simplified to:

1 A ∂ ∂ x A(Q e + Q i + Π) = S e,coll + S i,coll (15)
where we separated the time-averaged total energy flux into a kinetic contribution due to the electrons Q e , the ions Q i and the electromagnetic contribution Π. The kinetic energy flux of the electron was further separated into a flux of parallel energy Q e, and perpendicular energy Q e,⊥ (see appendix B).

To quantify the magnitude and direction of the energy exchanges between the electromagnetic field and the particles, the different terms of Eq. 15 were evaluated. To do so, the particles and field quantities were sampled in the steady state phase (after t = 30 µs, see fig. 2). First, particles were sorted in 120 spatial bins equally spaced along the axial direction. In each bin, the moments of particle distribution provided the total energy flux Q e and Q i , as detailed in appendix B. Second, the cross product of the electric and magnetic field provided the axial component of the Poynting vector. This vector was time-averaged over a period corresponding to an integer number of wave periods.

The results, plotted in Fig. 6, showed first that the energy source terms are negligible in the plume region (x ≥ 20 mm) compared to the source region (x < 20 mm). In the plume region, the magnitude of the source terms is below 1 × 10 4 W m -3 . For that magnitude, the signal is dominated by the statistical noise of the PIC simulation. This noise drowns the finer features of the source terms, especially for the electrons which are more prone to statistical noise. Nevertheless, this underlines that most of the energy exchanges take place in the source region and cannot explain the secondary peak for mean perpendicular kinetic energy e ⊥ observed in the plume region (Fig. 5).

Second, the collision source term is negative over the whole source region. Its magnitude is maximum near the backplate, where the neutral and plasma densities are higher, and decreases along the source axis. This behaviour is not unexpected, since this term lumps together the contribution of inelastic collisions and the diffusion model: these two processes are loss mechanisms for the plasma. Given that the plasma density and the neutral gas density decrease as we move away from the backplate, the collision frequency drops and the magnitude of the energy loss decreases.

Third, the sign and magnitude of the source terms reveal different phenomena. For x ≤ 3 mm, the terms linked to the Poynting vector are positive, while the source term due to the electron energy flux is negative. This indicates an energy conversion from the electron kinetic energy to the field energy. In parallel, the ion source term is positive, which points to a gain of energy in the sheath. For the 3 mm < x ≤ 10 mm range, the Poynting source term shows a negative peak, while the source term due to the electron energy flux is positive, with a peak centered on the ECR surface location. In that case, there is a transfer of energy from the field to the electrons. This latter feature corresponds to the ECR heating of the perpendicular energy mode of the electrons. In fact, this appears when considering the parallel and perpendicular contributions to the source term in Fig. 7. The perpendicular source term This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

dominates over the parallel source term, with a peak centered on the ECR surface. The axial extent of this peak shows that the perpendicular mode of the electrons is heated in a zone of about ∆x ECR ≈ 6 mm, i.e., from x = 3 mm to x = 9 mm. Considering that the ECR condition is only met on a specific surface, the presence of an extended region may seem surprising. However, as will be discussed later, most of the electrons in the magnetic field tube are confined and undergo a bouncing motion in the potential well formed by the electrostatic field and the magnetic mirror force. These bouncing electrons can cross the resonance surface with a significant parallel velocity, thus one may expect a shift of the resonance condition due to the Doppler effect. The width ∆x ECR can be compared with the expected value for a Doppler broadened resonance ∆x D in Eq. 16 [START_REF] Williamson | Selfconsistent electron cyclotron resonance absorption in a plasma with varying parameters[END_REF] .

∆x D = 2πv ωc B MS ∂ B MS ∂ x (16)
Using the electrons' mean axial velocity v is 1.1 × 10 6 m s -1 we obtain ∆x D = 4.7 mm. However, as it will be shown later in Fig. 8, there is a high dispersion for the values of v . Therefore, we can expect a much larger Doppler broadening for the fastest electrons. The electrons' axial velocity can reach values up to 3.0 × 10 6 m s -1 around the ECR zone. With this velocity, we can compute a maximum Doppler broadening of 7.8 mm, which means that 4.8 mm < ∆x D < 7.8 mm. The ECR heating zone obtained in the simulations is consistent with the one expected analytically, indicating that the Doppler effect is a good candidate to explain the width of the heating zone observed in Fig. 6.

Because of Doppler-broadening, ECR heating of the electrons can even occur when the ECR surface is outside the plasma source. Indeed, the fact that the plasma in the thruster can be sustained even with an ECR zone outside the coaxial chamber was demonstrated experimentally by Vialis [START_REF] Vialis | Développement d'un propulseur plasma à résonance cyclotron électronique pour les satellites[END_REF] , where the location of the resonance surface was placed at x = -0.17 mm and x = -0.77 mm. Doppler broadening is a possible explanation to this finding and this hypothesis was tested in simulations using the same parameters described in Table I but with an input microwave frequency of f EM = 2.9 GHz. It was possible to sustain a discharge for this frequency even though the ECR condition was meet at x = -1.78 mm (i.e., upstream of the plasma source).

In summary, the analysis of the power deposition shows that the energy transfer occurs mainly in the source region. The electrons absorb the wave energy over a Doppler-broadened volume. This energy goes preferentially to the perpendicular energy mode and leads to an anisotropy ratio T e,⊥ /T e, ∈ [2.5, 7.5]. This energy deposition from the field to the perpendicular energy mode can explain the first peak in perpendicular energy seen in Fig. 5. However, the source terms are negligible in the plume region and thus cannot explain the broad peak observed in this region.

B. Electron confinement in the magnetic nozzle

To determine the factors driving the evolution of the mean electron energy, we considered the electron distribution in the nozzle region. In Fig. 8 we plotted the normalized electron distribution in the velocity space v , v ⊥ plane. The distribution was plotted at different locations in the computational domain. The results in Fig. 8 show that as we move downstream into the nozzle, we see an increased electron population with high energies in the perpendicular direction. To understand this phenomenon, we must first get a more detailed description of the electron confinement, i.e., how they get trapped and under which conditions they can leave the thruster.

Three loss pathways are identified for the electrons produced in the source.

1. Cross field losses : electrons can diffuse across the magnetic field, due to anomalous transport, collision, etc. This is modeled by the phenomenological cross-field diffusion model presented in section II B 4.

2. Losses at the downstream of the nozzle electrons which have a kinetic energy sufficient to overcome the confining potential well are lost, alongside ions accelerated by this same potential drop.

3. Electrons than can overcome both the repelling potential of the plasma sheath at the close-end of the source and the mirror-force are collected on the dielectric plate and will contribute to its surface charge.

While the first loss mechanism does not depend on the kinetic energy of a single electron but rather on the mean kinetic energy at a given location, the two other mechanisms depend on the electron kinetic energy. An electron moving along the magnetic field will see both an electrostatic potential Φ (Fig. 3) and the magnetostatic field B. Depending on its initial kinetic and potential energies, its trajectory might have turning points (where v = 0) within the domain, or out of the domain. In the first case, this electron is confined in the potential well. In the latter case, the electron is lost, either downstream (loss pathway 2) or at the backplate (loss pathway 3). Now the limiting case between confined / unconfined electrons is when the turning points are located at the boundaries. This will define the necessary conditions for the electrons confinement. Let us note Φ BP the potential of the backplate. Neglecting the plasma-wave interaction, and noting µ the magnetic moment, the equation for the total energy of an electron is:

E total = 1 2 mv 2 + µB -eΦ (17)
From Eq. 17 we can say that the electron is oscillating in an effective potential given by U e f f = µB -eΦ, where µB represents the magnetic confinement as the electron moves towards the backplate while -eΦ is the electrostatic confinement given by the plasma potential. It can be seen in Fig. 9 for different arbitrary values of the magnetic moment taken from the results of the simulation. Its concave shape explains the confinement of the electrons in the ECR thruster inside this potential well. This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

Let us now consider an electron moving from an arbitrary initial point to a turning point at a position x 0 along the longitudinal direction, such that v (x 0 ) = 0. The energy conservation between any initial location and the turning point gives:

v 2 + v 2 ⊥ 1 - B x 0 B = - 2e m ∆Φ (18)
Where ∆Φ = Φ x 0 -Φ. Now let us consider the loss pathways 2 and 3 identified above.

For an electron lost to the downstream boundary (pathway 2), the confinement condition is obtained by setting x 0 = L. In that case, given the divergence of the magnetic field, we have 0 < 1 -B(L)/B < 1 and ∆Φ < 0. Thus equation 18 describes an ellipse in the vv ⊥ plane. Electrons in the ellipse have there turning points x 0 ≤ L and remain confined by the electrostatic well. Electrons out of the ellipse can overcome this electrostatic confinement and are lost downstream.

For an electron lost to the backplate (pathway 3), the confinement condition is obtained by setting x 0 = 0. Because the magnetic field is monotonically decreasing, 1 -B(0)/B < 0. In addition, Φ(0) = Φ BP . We define the loss cone angle as:

sin(θ ) = B B(0) (19)
Equation 18 can be recast as:

v 2 ⊥ = tan 2 (θ ) v 2 + 2e m ∆Φ (20)
Depending on the sign of ∆Φ, three cases are possible:

• ∆Φ = Φ BP -Φ = 0: The confinement of the electron at the backplate is given exclusively for the topology of the magnetostatic field according to the loss cone angle θ . Those electrons with values for v , v ⊥ such as they are located in the loss cone and will therefore be lost at the backplate (Fig. 10a).

• ∆Φ = Φ BP -Φ < 0: The backplate potential repels the negative charges and thus confines the electron. Consequently, a fraction of the electron in the loss cone will be reflected back and stay confined (Fig. 10b).

• ∆Φ = Φ BP -Φ > 0: The backplate potential attracts the electrons. Thus, even electrons out of the loss cone are collected. Therefore, the confined electrons are those that meet two conditions: they are not on the loss cone for the magnetic field, and they are energetic enough in the perpendicular direction to avoid being lost at the backplate thanks to the electrostatic acceleration towards it (Fig. 10c).

Fig. 3, displays the sign of ∆Φ in the nozzle. The combination of the loss conditions at the backplate (pathway 3) and downstream (pathway 2) delimits a confinement volume in phase space where electrons are confined, as shown in Fig. 10. Pitch-angle scattering, either caused by collisions with the neutral background or by the electromagnetic field, enables electrons to cross the confinement volume boundaries. Depending on which boundary is crossed, electrons are lost at the backplate or at the downstream side of the nozzle. Indeed, it is important to recall that the electron deconfinement is mainly driven by these two phenomena in the source region. Since the neutral background density decreases exponentially, most of the collisions occur in the source. In addition, as shown above, wave absorption happens over a few millimeters in the source. As a consequence, the electron deconfinement rate is driven by the wave interaction and the collisions:

• Interaction with the electromagnetic wave is akin to a scattering of the electron momentum [START_REF] Kaganovich | Electron Boltzmann kinetic equation averaged over fast electron bouncing and pitchangle scattering for fast modeling of electron cyclotron resonance discharge[END_REF] . After several passages through the ECR heating zone, the electron may gain enough energy to overcome the electrostatic barrier and escape into the plume.

• If the electron undergoes an elastic collision, it will randomly scatter its velocity vector. If the electron scattered momentum falls in the loss region defined by Eq. 20 and shown in Fig. 10, the particle is lost at the backplate.

If we now go back to the results in Fig. 8 for the electron distribution in the velocity space v , v ⊥ plane, we notice that as we move downstream into the plume, the confinement boundaries change. There is a transition from a confinement boundary as the one in Fig. 10b (∆Φ < 0) to the one in Fig. 10c (∆Φ > 0). This transition is a consequence of the fact that, as shown in Fig. 3 the plasma potential is greater than the backplate potential in the source Φ > Φ BP and lesser than Φ BP downstream. Thus, inside the coaxial chamber, the plasma potential is such that ∆Φ = Φ BP -Φ < 0 (Fig. 10b), and in the plume section it is such that ∆Φ = Φ BP -Φ > 0 (Fig. 10c). As a consequence, as we move downstream into the magnetic nozzle, the mean perpendicular kinetic energy can increase. However, this is not given by an additional heating phase but as a result of confining only a highly energetic electron population in the perpendicular direction. Those electrons with a low perpendicular kinetic energy (i.e., below the dashed line) are lost at the backplate as previously described. It can be seen as a filtering process where only the hot electrons are confined, and this is what we see when plotting the electron perpendicular kinetic energy in Fig. 5. Further downstream of the magnetic nozzle, the magnitude of the potential well to the end of the nozzle decreases, while the magnitude of the attracting potential drop to the backplate increases. This results in a narrower distribution for the confined population and finally a decrease in the mean perpendicular kinetic energy.

IV. CONCLUSIONS

We have performed electromagnetic full-PIC simulations of the ECR thruster using a 1D3V model that allowed us to shed light onto some of its working principles. The results confirmed the expected anisotropic behavior for the electrons' energies in the direction perpendicular and parallel to This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

the magnetic field lines and a peak for the mean perpendicular energy near the resonance zone. The microwave energy injected at the backplate of the thruster propagates through the coaxial chamber while being absorbed by the electrons increasing their kinetic energy perpendicular to the magnetic field lines. The absorption takes place exclusively inside the coaxial chamber on a zone of 6 mm around the resonance condition. This zone is coherent with the predicted value from Doppler broadening. The width of this heating zone may explain why the thruster works even with a configuration in which the resonance condition is met outside the coaxial chamber [START_REF] Vialis | Développement d'un propulseur plasma à résonance cyclotron électronique pour les satellites[END_REF] . From a practical point of view, this feature improves the reliability of the thruster, since it means that the thruster can still operate even when the magnetostatic field decreases, for example due to excessive heating of the magnets. The results also show, unexpectedly, that the electrons' mean perpendicular energy has a second peak in the plume due to the confinement of highly energetic electrons. The confinement is determined by the backplate's potential, the magnetostatic field, and the potential drop on the plume. As a consequence, there is a population of trapped electrons with significant perpendicular kinetic energy in the downstream region of the magnetic nozzle. The existence of doubly trapped electron population has been investigated using a kinetic model with a paraxial approximation similar to this work [START_REF] Ahedo | Macroscopic and parametric study of a kinetic plasma expansion in a paraxial magnetic nozzle[END_REF] . While this work was assuming the shape of the initial distribution function, it has also been seen in the case of an anistropic distribution that the perpendicular electron temperature could increase in the diverging part of the nozzle [START_REF] Plaza | Physics of plasma plumes accelerated by magnetic nozzles: An experimental and theoretical research[END_REF] . We can speculate that these high temperature electrons trapped in the plume could be important to drive some instabilities observed in the plume that are thought to enhance cross-field transport of electrons and thus play a role in the detachment. In particular, Lower Hybrid drift instabilities have been recently observed in diverging magnetic nozzle. In these experiments, diamagnetic drift v D = ∇p e⊥ × B/enB 2 was identified for the primary energy source for the instability [START_REF] Hepner | Wave-driven non-classical electron transport in a low temperature magnetically expanding plasma[END_REF] . The trapped electrons could enhance the radial gradient in perpendicular pressure and thus enhance the diamagnetic drift.

∂ ∂t

m e v 2 2 f e d 3 v + ∂ ∂ x • m e v 2 2 v f e d 3 v (B1) +q e v 2 2 (E + v × B) • ∂ f e ∂ v d 3 v = m e v 2 2 ∂ f e ∂t col d 3 v
Where the right hand side lumps the contribution of collisions and the loss model detailed in section II B 1 and II B 4, respectively. Eq. B1 can be rewritten as:

∂ ε e ∂t + ∇ • Q e = -j e • E + S coll Q e = m e v 2 2 v f e d 3 v (B2) ε e = mev 2 2 f e d 3 vr (B3) S e,coll = mev 2 2 ∂ fe ∂t coll d 3 v (B4)
The same procedure can be applied to the ions:

∂ ε i ∂t + ∇ • Q i = -j i • E + S coll Q i = M i v 2 2 v f i d 3 v (B5) ε i = M i v 2 2 f i d 3 v (B6) S i,coll = M i v 2 2 ∂ f i ∂t coll d 3 v (B7)
This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

Considering the electron population, the total heat flux can be written as:

Q e = q e + P e • u e + n e u e (e K + E K ) (B8)
Where the density is given by n e = f e d 3 v and the macroscopic velocity by n e u e = f e vd 3 v. The random part of the velocity is c = vu e The different terms are then:

q e = m e 2
A similar expression can be written for the ions. if we make use of the Poynting theorem, we can relate the time averaged joule term to the divergence of the pointing flux Π = E×B µ 0 :

A(x) (j e + j i )

• E + ∂ ∂ x Π = 0 (B14)
If we drop the < ... > symbol for simplicity, recalling that all quantities are time-averaged, one obtains equation (15).

Images/fig_ECRA-eps-converted-to.pdf This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

c 2 cd 3

 3 the one-dimensional approximation, we are considering only the axial (parallel) part of the total heat flux. Thus, after averaging over a period of the incoming wave : < ... >= 1 T ...dt, we obtain: ∂ ∂ x A(x) q e + A P e u + P e⊥ u ⊥ +A(x) u e + n e u e; (e K + E K ) = -A(x) j e • E + A(x) S e,coll

FIG. 2 :FIG. 3 :FIG. 5 :FIG. 6 :FIG. 7 :FIG. 8 :FIG. 10 :

 23567810 FIG. 1: ECR thruster: (a) Schematic view of the coaxial source. The magnetic field lines are shown in red. The dashed surface corresponds to the flux tube (b) Schematic view flux tube used for the quasi-1D model. The exit plane of the coaxial source of length L S is represented by the dashed line. The end of the computational domain is reached at x = L. The axial magnetic profile and tube cross section along the axis are shown in red and blue, respectively.

TABLE I :

 I Simulation parameters for the electromagnetic full PIC simulations using the quasi-one-dimensional model.

	Parameter	Description	Value
	∆t	Time step	1.6 ps
	∆x	Mesh spacing	167 µm
	C	CFL condition	2.87
	f EM	Microwave frequency	2.45 GHz
	L S	Coaxial chamber length	20 mm
	x ECR	ECR surface position	6.7 mm
	W L D	Weight for the charged particles Computational domain length	2 × 10 5 100 mm
	n n0 L n	Maximum number density of neutrals 8 × 10 19 m -3 Neutral density characteristic length 1.0 cm
	A L	Cross-sectional area for the loss module	1 cm 2
	∆t ions	Time step to push the ions	10∆t
	∆t coll	Time step for collisions	10∆t

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0124834

Appendix A: Electron motion integration in the Quasi-1D model On one hand, the model assumes that the axial static magnetic field in the flux tube depends on x only

In the quasi-1D model, the field remains constant across the section of the tube. On the other hand, in the ECR thruster the electromagnetic part of the magnetic field (the part due to the propagation of the electromagnetic power injected in the source) remains negligible compared to the static part. Indeed, assuming a locally plane wave, the Poynting vector S = E × B/µ 0 . Since E ≃ cB and S ≃ 1 W cm -2 , the order of magnitude for the electromagnetic component of the magnetic field is B ≃ 1 × 10 -2 mT. The static part of the field is between 10 mT and 10 mT, much greater than the electromagnetic component. Therefore, using the divergence equation for the static magnetic field and neglecting the electromagnetic component, it is possible to obtain the radial component of the magnetic field :

Assuming this value for the radial part of the magnetic field ensures that the divergence condition is automatically enforced for the static part of the field. The electron guiding center is on the flux tube centerline. The Larmor radius of the electrons is given by:

Where

The gyromotion of the electron is shown in Fig. 11. Thus, knowing the particle velocities it is possible to obtain the phase angle θ in its gyromotion, as given in Eq. A5.

Knowing the phase angle, sine and cosine, the B y an B z components can be deduced from eqs. A2 and A5.

Appendix B: Energy equation

For the energy equation for the particles we consider the second order moment of the Vlasov equation. Multiplying the Vlasov equation for the electrons by m e v 2 /2, we obtain:


S Bernard
email: sbernard@mnhn.fr

O Beyssac

A Ollila

G Lopez-Reyes

J Manrique

S Le

Perseverance S Bernard

S Le Mouélic

P Beck

O Forni

P Pilleri

A Cousin

O Gasnault

P Y Meslin

G Travis

E Clavé

C Royer

R C Wiens

S Maurice

IRRADIATION OF ORGANICS ON MARS: EVOLUTION OF THE RAMAN SIGNAL OF THE ERTALYTE TARGET ABOARD PERSEVERANCE

whether they are published or not. The documents may come

Introduction:

The Perseverance NASA rover is exploring the Jezero crater on Mars using an arsenal of spectroscopic tools [1,2] to interrogate and select the samples that will be returned to Earth. Among the main targets are rocks or sediments possibly containing ancient organic traces of life [3]. In fact, although the surface of Mars is globally no longer habitable, life may have existed on early Mars and traces of it may be preserved within the rocks still lying at the surface. Yet, because the surface of Mars is bombarded by ionizing radiation, it remains hostile for organic molecules. The thin CO 2 atmosphere of Mars absorbs most X-rays and far-UV radiation, but it lets most UV photons reach the surface [4]. Even though UV photons do not penetrate that much below the surface [5], a number of experimental studies have shown that the exposure to UV radiation leads to the rather fast degradation of organic compounds [6][7][8]. But these experimental studies may not completely mimic what is really happening on Mars. In fact, the surface of Mars is not only bombarded by UV, but also by gamma rays, solar energetic protons, and galactic cosmic rays [9]. Documenting the degradation of organic materials under true Martian conditions, i.e. at the surface of Mars, thus appears necessary.

Methods: Among the SuperCam calibration targets [10][START_REF] Cousin | th Lunar and Planetary Science Conference 2023[END_REF], Perseverance is carrying an organic target made of polyethylene terephthalate (the Ertalyte® target), i.e. an organic polymer relatively resistant to UV radiation and exhibiting a strong Raman signature (Fig. 1). The Ertalyte target has been measured using the SuperCam Raman every 40 to 50 sols since the landing. Each spectrum corresponds to the accumulation of 100 shots collected using a gate of 100 ns, with the laser tuned at 110 A and an intesifier gain of 3200.

Results: White at landing, the Ertalyte target has turned yellow to brown with time as revealed by RMIs (Fig. 2). The Raman spectrum of the Ertalyte has also changed as a function of time (Fig. 3), with a modification of the shape of the background and a decrease of the Raman signal relative to the background. Of note the background has strongly increased with time but this increase remains difficult to quantify, the total amount of counts collected on the spectrometer varying from one sol to another. The Raman spectrum of the irradiated Ertalyte target still exhibits all the features observed in the spectrum of pristine Ertalyte, even after 627 sols on Mars. Still, exposure at the surface of Mars has strongly impacted the molecular structure of the Ertalyte target by triggering the creation of electronic defects and/or colored centers, which is consistent with experimental reports [6][7][8]. Of note, such good agreement with experimental studies suggests that using only UV photons to simulate in the lab the radiative environment of the surface of Mars can be considered as a good firstorder approximation.

The evolution of the area of the peak at 1615 cm -1 (i.e. the peak attributed to the presence of aromatic cycles) roughly follows a first-order logarithmic law with time (Fig. 4) and reveals that the degradation of the Ertalyte is still ongoing even after almost a year on Mars. Whether or not the degradation of the Ertalyte target will reach a plateau remains on open question.

Conclusion:

The present results confirm that the conditions existing at the surface of Mars are detrimental to the preservation of organic materials. In addition to demonstrating that the SuperCam Raman can detect organic materials on Mars, the Ertalyte target aboard Perseverance allows us to document the impact of irradiation on such materials, providing crucial information for the search for traces of life on Mars. Good constraints on the degradation of the Ertalyte target will help determining what could/should be expected to be found or not in Martian rocks.

Figure 1 :

 1 Figure 1: Raman spectrum the pristine Ertalyte Target obtained with a time-resolved setup at IMPMC.

Figure 2 :Figure 3 :

 23 Figure 2: RMIs of the Ertalyte Target on different sols

Figure 4 :

 4 Figure 4: Evolution of the area of the main Raman peak at 1615 cm -1 of the Ertalyte target with time

Acknowledgments:

The authors acknowledge the M2020 engineering and scientific teams.