W Villafana
email: wvillafa@pppl.gov

B Cuenot

O Vermorel

3D Particle-In-Cell studyof the Electron Drift Instability in a Hall Thruster using unstructured grids 3D Particle-In-Cell study of the Electron Drift Instability in a Hall Thruster using unstructured grids

    3D particle-in-cell study of the electron drift instability in a Hall Thruster using unstructured grids W. Villafana, B. Cuenot, O.

3D Particle-In-Cell study of the Electron Drift Instability in a Hall Thruster using unstructured grids 2 in this approach, Minelli, Taccogna et al. [START_REF] Minelli | How to build pic-mcc models for Hall microthrusters[END_REF][START_REF] Taccogna | Three-dimensional particle-in-cell model of Hall thruster: The discharge channel[END_REF] proposed another strategy to build a 3D PIC model for a miniaturized version of an SPT100. Their methodology consists in reducing the three space dimensions by a factor f , set to 10 in their work. In an effort to properly "scale down" the problem, some physical parameters had to be adjusted. Results showed that rescaled values of main parameters fairly agree with data that Taccogna, Minelli et al. could gather in the literature. Besides, they found that secondary electron emission has a modest effect on the anomalous transport but acts as an additional saturation mechanism for the EDI.

Overall, Taccogna and Minelli paved the way for insightful and more accessible 3D PIC simulations, but they also recognized limitations to their approach. First of all, scaling down the geometry by a factor f unavoidably increases the area to volume ratio, which may enhance the effects of the walls on the plasma dynamics. Moreover, other important physical parameters such as the ion convection time or the cyclotron frequency Ω ce are respectively multiplied by a factor f -1 and f 1 . Increasing Ω ce leads to an overestimation of the growth rate of the ECDI [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] while reducing the convection time is likely to modify the saturation mechanism of the instability. Finally, the study focused on the sole channel region and therefore plume effects could not be studied.

In this paper, based on previous work [START_REF]Numerical Particle-In-Cell studies of Hall thrusters using unstructured grids[END_REF] , 3D simulations of HT are presented with the objective to overcome these different limitations. First, the computational domain considers both the channel and the near plume region, with reasonable axial and radial lengths to avoid additional assumptions at the channel exit and preserve a reasonable area to volume ratio. Second, the true values of ion to electron mass ratio m i /m e and vacuum permittivity are retained to preserve plasma dynamics. Simulations were performed with the in-house PIC code, AVIP PIC, which was successfully verified in three international benchmarks [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF][START_REF] Turner | Simulation benchmarks for low-pressure plasmas: Capacitive discharges[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . A specificity of AVIP PIC is its capability to deal with unstructured grids. To our knowledge, such a functionality is unparalleled in the HT community and is a feature that interests the industry by its capability to model any geometry, whatever its complexity.

A detailed description of the present 3D model, along with its limitations, is given in Section II. In Section III the development of the EDI and other oscillatory phenomena during the transient state of the simulation is studied. Section IV gives an overview of the steady state. Based on theoretical results, a careful spectral analysis of the EDI and other fluctuations is also performed. Section V compares the results obtained in the current 3D configuration with results extracted from an analogous 2D configuration in order to highlight 3D effects. Finally, in Section VI, the statistical convergence of the 3D PIC simulation is explored.

II. NUMERICAL MODELING, ASSUMPTIONS, AND LIMITATIONS

All simulations were performed with the in-house code AVIP PIC, an explicit PIC-MCC code developed at CERFACS. AVIP PIC was developed from ABVP, a CFD code dedicated for reactive flows [START_REF] Schonfeld | Steady and unsteady flow simulations using the hybrid flow solver avbp[END_REF][START_REF] Esclapez | A statistical model to predict ignition probability[END_REF][START_REF] Male | Direct numerical simulations and models for hot burnt gases jet ignition[END_REF] and aerodynamics problems [START_REF] Queguineur | Large eddy simulations and global stability analyses of an annular and cylindrical rotor/stator cavity limit cycles[END_REF] in industrial geometries. AVIP PIC inherits AVBP features and in particular the capability of working with unstructured grids for highperformance applications [START_REF] Gicquel | High performance parallel computing of flows in complex geometries[END_REF] . The PIC-MCC method has been implemented with a domain decomposition, which is achieved with the external PARMETIS library [START_REF] Karypis | MeTis: Unstructured Graph Partitioning and Sparse Matrix Ordering System, Version 4.0[END_REF] .

A. Computational domain

The computational domain has a "T shape" which is extruded in the azimuthal direction. Therefore, no curvature effects is accounted for in this configuration. This choice simplifies the understanding of the ongoing physics and facilitates the comparison with previous work. Indeed, the dimensions of the domain were chosen based on those used in two now well-established 2D PIC benchmarks [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . The axial length L x is set to 2.5 cm while the azimuthal length L y is 1 cm. In the channel, the radial gap L z is 1 cm and in the plume region the total height is 2L z . The channel has a length L ch = 0.85 cm. At its exit, rounded edges avoid an unrealistically high electric field that would exit at sharp tips. There, the local curvature has a radius L R = 0.1 cm. A schematic view of the domain in the radial-axial plane is given in Figure 1 (a) and the whole 3D domain is shown in Figure 1 (b).

Because of these rounded edges, a classical Cartesian mesh is not used here. Instead, a fully unstructured mesh made of tetrahedral cells is considered and was generated with the Centaur software [START_REF] Borras | Centaur: The system[END_REF] . The typical cell size ∆l was chosen such that,

∆l = V 1 3 c < λ D (1)
where V c is the cell volume and λ D is the Debye length:

λ D = ε 0 k B T e /(n e e 2 ) (2)
with ε 0 the vacuum permittivity, k B the Boltzmann constant, n e the electron density and e the elementary charge. According to Langdon and Okuda [START_REF] Langdon | Theory of plasma simulation using finite-size particles[END_REF][START_REF] Okuda | Verification of theory for plasma of finite-size particles[END_REF] , Equation ( 1) is sufficient to obtain a good approximation of the behavior of a real plasma on a Cartesian mesh. It is assumed that this still holds for the present unstructured grid.

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. B. Charged particles

PLEASE CITE THIS ARTICLE AS

Initialization and pusher

The plasma is collisionless and made of singly charged ions X + e and electrons e -only. The simulation starts with a uniform plasma density n 0 at thermal equilibrium with electron and ion temperatures T e,0 = 10 eV and T i,0 = 0.5 eV. In each tetrahedral cell N ppc,ini = 120 macroparticles are initialized for each type of particle. As the cells have slightly different shapes and volumes, the statistical weight q f is automatically adjusted in each cell to ensure a density n 0 .

A Boris algorithm [START_REF] Boris | Relativistic plasma simulation-optimization of a hybrid code[END_REF] is used for the displacement of the electrons, which are sensitive to the magnetic field. For the ions that are only affected by the electric field, a Leap-Frog scheme [START_REF] Birdsall | Particle-in-cell charged-particle simulations, plus monte carlo collisions with neutral atoms, pic-mcc[END_REF] is preferred. The ion position and speed are updated every

f sub = 5∆t
, where ∆t is the fixed time step. The last computed electric field was used to update the velocity of ions. From the stability criterion from Birsdsall et al.,

ω pe ∆t < 0.2, (3)
where ω pe = n e e 2 /(m e ε 0 ) is the plasma frequency, the time step is set at ∆t = 5 × 10 -12 s with ω pe = 1.79 rad s -1 at initialization.

Ionization

Ionization is indirectly accounted for by using a source term similar to previous 2D configurations [START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . The 2D radial-axial ionization profile is given by:

     if (x, z) ∈ [x 1 , x 2 ] × [z 1 , z 2 ], S(x, z) = S 0 cos x-xm x 2 -x 1 cos z-zm z 2 -z 1 else, S(x, z) = 0, (4)
where S 0 = 9.55 × 10 23 m -3 s -1 is the strength of the source term, (x 1 , x 2 ) and (z 1 , z 2 ) are respectively the axial and radial limits of the ionization zone. Their values are taken from the 2D configurations: x 2x 1 = 0.75 cm like in Charoy et al. [START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] while z 2z 1 = 0.86L z = 0.86 cm as in Villafana et al. [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] . Finally, x m = (x 2 + x 1 )/2 and z m = (z 2 + z 1 )/2 are respectively the mean axial and radial values. An overview of the ionization zone is provided in Figure 2. The total extracted current can be computed with:

I tot = x 2 x 1 z 2 z 1 Ly 0 S(x, z) dx dy dz ⇐⇒ I tot = e 4 π 2 S 0 (x 2 -x 1 )(z 2 -z 1 )L y , (5)
Hence, the corresponding total current density J m through the surface channel L y L z can be calculated as:

J m = I tot L y L z = e 4 π 2 S 0 (x 2 -x 1 ) (z 2 -z 1 ) L z . (6)
The value J m = 400 A m -2 used in this work is significantly lower than the typical values encountered in HTs 1 (around 1000 A m -2 ) but is sufficient to capture oscillatory phenomena and the expected working conditions of a HT, as it will be shown in Section IV. Pairs of X + e /e -are injected in the ionization zone delimited by (x 1 , x 2 ) and (z 1 , z 2 ) following the profile S(x, z). From a practical implementation point of view, the (x P , y P , z P ) coordinates of each pair are given by three random numbers R 1 , R 2 and R 3 :

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.0133963

     x P = arcsin (2R 1 -1) x 2 -x 1 π + x m y P = R 2 L y z P = arcsin (2R 3 -1) z 2 -z 1 π + z m . (7)
Furthermore, the velocity components are sampled from a Maxwellian distribution with injection temperatures T e,0 and T i,0 for the electrons and ions, respectively. Finally, the weight of injected particles is constant and defined as:

q f = n 0 ∆l 3 N ppc,ini , (8)
with n 0 the initial plasma density, N ppc,ini , the initial number of particles per cell and ∆l, the mean cell size.This gives q f ≈ 256.

Using this ionization model inevitably modifies the actual HT physics but makes it easier to compare with analogous 2D PIC simulations, using a similar injection mechanism as in references [START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . This also allows to avoid simulating the actual ignition of the HT which poses various numerical challenges that are beyond the scope of this study.

Boundary conditions for particles

All particles leaving the domain are removed from the simulation except at the periodic boundary conditions at y = 0 and y = L y . In particular, this means that no dielectric layers are taken into account in this work. At the top, bottom and right boundaries (see blue edges in Figure 2) a quasi-neutral (QN) cathode model is used, following previous recommendations [START_REF] Charoy | Numerical study of electron transport in Hall Thrusters[END_REF][START_REF] Szabo | Fully Kinetic Numerical Modeling of a Plasma Thruster[END_REF][START_REF] Cho | The effects of cathode boundary condition on particle simulation of a spt-100-like Hall thruster[END_REF] .

The number of electrons to be injected N e,in j is calculated by summing the net number of electrons in each cell as,

if ∑ N cell,BC j N e -N i > 0, then N e,in j = ∑ N cell j N e -N i else N e,in j = 0 (9)
where N cell,BC is the number of cells in contact with the cathode boundary and N e and N i are respectively the number of electrons and ions in each cell. Electrons are uniformly injected along the cathode boundary at x = L x and their orthogonal velocities are sampled from a half-Maxwellian velocity at a temperature T e,0 . This process ensures a neutralized plasma in the plume, which is the expected role of the cathode.

C. Electromagnetic fields

Electric field

The Poisson equation is self-consistently solved with the external library MAPHYS [START_REF] Agullo | Robust preconditioners via generalized eigenproblems for hybrid sparse linear solvers[END_REF] , designed for unstructured meshes and combining both direct and iterative methods for fast and accurate results. Dirichlet conditions are imposed on all boundaries of the domain for the potential φ . The imposed constant values are specified in Figure 2. Note that in addition to the right side of the domain, the cathode also includes the top and bottom boundaries where the potential is zero. Neumann conditions were also tested at the top and bottom boundaries as they seemed more realistic but they led to the creation of an important potential bulk (several thousands of Volts) in the plume that was nonphysical. Furthermore, the absence of a dielectric layer at the walls implies to manually ensure the continuity of the potential. The downstream section of the walls, in yellow in Figure 2, is in contact with the cathode, and thus is set to φ W,2 = 0 V. Between the anode and the wall, the potential must be continuous to avoid electric arcs. This is achieved by assuming that the potential decreases linearly from φ a = 200 V to 0 V over a distance L W .

In this region, the potential is expressed as

φ W,1 = - φ a L W x + φ a , (10)
where L W is set to 0.375 cm. Other more realistic values for L W and shapes for the φ W,1 function would be possible. Nevertheless, imposing a linear decrease as described in Equation ( 10) turns out to be sufficient to reproduce HT operating conditions as shown by Joncquières et al. [START_REF] Joncquieres | A fluid formalism for low-temperature plasma flows dedicated to space propulsion in an unstructured high performance computing solver[END_REF] .

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.0133963

x z Wall y Anode Cathode (QN) φ W,2 = 0 φ W,1 (x) = -φa LW x + φ a φa = 200 φ c = 0 Figure 2.
Boundary conditions and ionization zone for the 3D PIC simulation (radial-axial view). The ionization source term shape is given by Equation ( 4) and represented with a blue gradient. The anode is in red and set at the potential φ a = 200 V. The cathode is in blue and uses a quasi-neutral (QN) model. Its potential is φ x = 0 V. The walls are absorbent and the potential is split into two parts. In the upstream channel, in orange, the potential φ W,1 is a linear function which decreases from the anode potential φ a down to 0 V. The other part of the walls, in yellow, has a zero potential φ W,2 .

Magnetic field

The choice of the magnetic field is less trivial than in the 2D cases as in reality, magnetic poles around the channel generate a 3D magnetic field topology. In order to calculate such a magnetic field, the ∇ • B = 0 equation has to be solved with appropriate boundary conditions. However, even for apparently basic configurations such as permanent ring cusp magnets [START_REF] Ravaud | Analytical calculation of the magnetic field created by permanent-magnet rings[END_REF] , analytical formulas are incredibly complicated. One solution is to use experimental measurements as Taccogna et al. did [START_REF] Taccogna | Three-dimensional particle-in-cell model of Hall thruster: The discharge channel[END_REF] , but in this case, the numerical setup must be as close as possible to the real thruster. Besides, experimental data can be difficult to obtain as they can be classified by manufacturers. If no data are available, one can turn to open source solutions. As an example, the Finite Element Method Magnetics (FEMM) software [START_REF] Meecker | Finite Element Method Magnetics[END_REF] can generate a 2D magnetic field topology for a variety of configurations, but requires knowledge of the electronic circuit of the magnetic coils. This approach was used by Garrigues et al. [START_REF] Garrigues | Model study of the influence of the magnetic field configuration on the performance and lifetime of a Hall thruster[END_REF] in their HT simulations. Finally, another interesting approach consists of a recently developed open source Python package capable of calculating the magnetic field of simple permanent magnets [START_REF] Ortner | Magpylib: A free python package for magnetic field computation[END_REF] . The total magnetic field stemming from several magnets is obtained by the superposition principle.

All the above strategies could have been used in the present work. As of today, they seemed however premature as it is important to understand and verify the behavior of the 3D plasma dynamics in a simplified configuration first. Therefore, it was chosen to remain as close as possible to the 2D configuration of Charoy et al. [START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] assuming a purely radial magnetic field, varying in the axial direction with a Gaussian shape:

B(x) = a k exp - (x -x Bmax ) 2 2σ 2 k + b k u z , (11)
where k = 1 for x ≤ x Bmax and k = 2 for x > x Bmax , and σ 1 = σ 2 = 0.625 cm. The four coefficients a k and b k are chosen such that

B(x = 0) = 6 mT, B(x = L x ) = 1 mT, B(x = x Bmax ) = 10
mT and by satisfying the continuity at x = x Bmax .

Such a choice means that iso-potential surfaces at the top (y = 2L z ) and bottom (y = 0) boundaries do not follow magnetic field lines, which should be the case in a real HT [START_REF] Garrigues | Model study of the influence of the magnetic field configuration on the performance and lifetime of a Hall thruster[END_REF] . However, a future work to be published will show that the main oscillatory features remain intact even with a more realistic B field [START_REF] Villafana | 3D Particle-In-Cell modeling of anomalous electron transport driven by the Electron Drift Instability in Hall thrusters[END_REF] .

The 3D PIC simulation was performed on Irene cluster at the Atomic Energy Commission (CEA) using AMD Rome nodes made of 2 × 64 cores each. For best performances, 1920 cores were used for a total cost of 1.4 million CPU hours to simulate 10 µs.

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.0133963

III. TRANSIENT STATE

In this section, the temporal evolution of the simulation is first described with emphasis on the transient phase. The development of the EDI and oscillatory phenomena in the computational domain is also highlighted.

The temporal evolution of the space-averaged ion density (over the whole domain) is shown in During the second stage of the transient, phase II, from 1.2 to ∼ 5 µs, axial losses at the cathode overtake the creation of particles, as ions leave the computational domain. At first, up to 2.1 µs, particle losses remain moderate before increasing significantly.

Finally, in phase III, from 5 µs, most of the wall and axial losses at the cathode start to be balanced out by the constant ionization source term and the ion profile stabilizes. More details on the transient state can be provided by looking at the main planes of interest of the domain, starting with the axial-azimuthal view xy. Figure 4 presents the centered xy plane for the axial electron current J e,x , and the azimuthal electric field E y at selected times of interest. During the linear-like phase, at t = 0.20 µs, the 2D map of the azimuthal electric field E y exhibits fluctuations in the azimuthal direction that start to develop upstream, near the anode. These fluctuations have initially a short wavelength of the order of 1 mm and propagate in the E × B direction, similarly to what was observed in previous 2D axial-azimuthal PIC studies [START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . Because of this early similarity, these propagating fluctuations are identified as the EDI. This will be confirmed by theoretical arguments provided in Section IV B. Shortly after initial growth, the EDI is convected toward the plume by the transient ion wave. Finally, from around 5 µs, the EDI seems to be established in almost the entire domain. In the channel, the wave vectors seem to have both axial and radial components while it is mostly axial in the plume. This was also observed in a pure 2D axial-azimuthal configuration [START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . In the present 3D configuration the azimuthal wavelength is measured to be λ y,ch ≈ 628 µm in the channel while it is around λ y,pl ≈ 913 µm in the plume. The axial electron current J e,x is also very interesting. Similarly to E y , J e,x presents short-wavelength azimuthal oscillations at the beginning of the simulation. When the EDI takes over the whole domain, from 5 µs, its associated short-wavelength fluctuations are only visible inside the channel, where the electron density is highest. Besides, another subtle pattern arises inside the channel in the azimuthal direction. At 5 µs and 7.5 µs, a large fluctuation is observed with an azimuthal wavelength of the order of L y = 1 mm. 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.0133963 To go further in the analysis the radial-azimuthal zy plane is examined in Figure 5 where the axial electron current J e,x , the in Figure 4, a pure azimuthal wave due to the EDI, develops in the first moments of the simulation (t = 0.20 µs). From t = 1.20 µs, the short-wavelength (λ EDI ∼ L y /10) is still clearly visible for the azimuthal electric field E y but a radial pattern of wavelength 2L z seems to emerge in the electron radial current J e,z . This larger oscillation also has an azimuthal component that distinctly appears at t = 7.50 µs. The corresponding azimuthal wavelength is approximately L y . Moreover, the pure azimuthal wave in the E y 2D map sees its wavelength roughly divided by two when approaching the steady state. Overall, this large radial-azimuthal wave resembles the MTSI encountered in previous works [START_REF] Janhunen | Evolution of the electron cyclotron drift instability in twodimensions[END_REF][START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] . This point will be further detailed in Section IV B.

0 2 4 6 8 10 t [µs] 0.8 1.0 1.2 1.4 n i [m -3 ] ×10 17 (I) (II) (III)

IV. STEADY STATE A. Main plasma variables

The steady state is reached after about 10 µs. The proper functioning of the thruster is verified by analyzing Figure 6. In Figure 6 (a), the great ions acceleration can be observed at the channel exit where the plasma flow becomes supersonic. The classical definition of the Mach number M is adopted here:

M = u i c s , (12)
where u i denotes the magnitude of the ion hydrodynamic velocity and c s , the ion speed of sound defined as:

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.0133963

c s = k B T e m i , (13)
where m i is the ion mass. Additionally, M ≥ 1 at the inner and outer walls of the channel and also at the vicinity of the anode at left, which marks the beginning of their respective sheath, at a wall distance of roughly 1 mm at the inner/outer walls, of the order of the local Debye length. The M = 1 threshold corresponds to the Bohm sheath criterion condition near these walls. As detailed in Figure 7 (a), ions accelerate in the axial direction from ∼ 0 to 5 × 10 3 m s -1 as the plasma travel toward the exit before expanding into the plume. In the channel, the velocity vector is mostly dominated by the radial v i,z component. Besides, Figure 7 (a) shows that, in the plume, the Mach number continues to increase in the radial direction. Some supersonic ions might be attracted back to the exterior wall, but most of them keep accelerating to the cathode boundary of the simulation domain. In Figure 7 (b), from the rounded edge (channel exit) to the top cathode boundary, the potential drop between the exterior wall and the plasma bulk decreases and becomes negligible at z ≈ 1.8 cm = 0.9 × 2L z (∆φ ≤ 5 V). This is primarily due to the depletion of the plasma as the latter expands into the plume. At the junction between the exterior wall and the top cathode boundary, the plasma density has been reduced by a factor of 100 with respect to the channel (where n i ≈ 5 × 10 17 m -3 ). Besides, the plasma has become quasi-neutral as imposed by the quasi-neutral (QN) cathode condition described in Section II B 3.

The ion acceleration in the axial direction is caused by an important potential drop seen in Figure 6 (b). The subsequent electric field streamlines are mostly directed toward the plume and have both radial and axial components, suggesting that the divergence beam is important. This is coherent with the choice of using Dirichlet boundary conditions for the Poisson equation at the top and bottom boundaries in the plume. An important radial gradient of potential is also present in the downstream section of the channel, also highlighted by the iso-φ lines in Figure 6 (e), leading to a strong radial electric field. In Figure 7 (b), providing a zoomed-in picture of the area, shows the potential drops from 100 V to 0 V over 0.25 mm, leading to a local radial electric field E z ∼ 400 kV m -1 . This is a direct consequence of the choice of imposing the φ W,1 and φ W,2 profiles along the channel walls. In the upstream section of the channel, the φ W,1 profile along the wall is also visible as shown by Figure 7 (b), as the profile falls from 200 (anode) to 0 V (at nearly half of the channel), passing by 100 V at approximately the quarter of the channel. Overall, the selected boundary conditions φ W,1 and φ W,2 artificially bend iso-potential lines in the channel in comparison to a standard Hall Thruster [START_REF] Garrigues | Operation of a low-power Hall thruster: comparison between magnetically unshielded and shielded configuration[END_REF] , which probably enhances losses at the walls. Besides, in the channel, a strong electron drift current is observed in the azimuthal direction under the effect of the E and B fields (Figure 6 (c)).

The global electron transport is shown in Figure 6 (d) with the help of the axial electron current J e,x = -en e u e,x . With this definition, J e,x > 0 means that u e,x < 0. From the streamlines, we can see that electrons effectively enter the channel despite the presence of recirculation zones near the walls. In the channel, the presence of these recirculation zones might be related to the radial azimuthal wave seen in Figure 5. In the plume, another recirculation zone is also visible. The situation is unclear here:

the streamlines seem to follow the curvature of the azimuthal instability carried by E y (see Figure 6 (f)). However, at this point, it is difficult to provide a satisfying explanation for these recirculation areas.

In Figure 6 (d), another observation can be made regarding the cathode. First some streamlines issuing from the right boundary reach the top and bottom plume boundaries, which means that electrons are lost and will be re-injected by the cathode. This behavior is even more evident when looking closely at the right plume boundary where we note the current being alternatively positive and negative. Actually, a significant portion of newly injected cathode electrons leaves the domain after a few iterations, which is not desirable as it artificially disturbs the local charge distribution. Cho et al. [START_REF] Cho | The effects of cathode boundary condition on particle simulation of a spt-100-like Hall thruster[END_REF] noted a similar problem with their cathode model. A possible solution would be to shift the injection surface slightly into the domain to help electrons stay inside the domain as it was done by Charoy [START_REF] Charoy | Numerical study of electron transport in Hall Thrusters[END_REF] . Another option would be to replace the random injection of electrons with a targeted injection as suggested by Szabo [START_REF] Szabo | Fully Kinetic Numerical Modeling of a Plasma Thruster[END_REF] . In this approach, electrons are preferentially injected where the charge difference between ions and electrons is the highest.

Figure 6 (e) shows the space distribution of the electron temperature T e . With a maximum of approximately 50 eV, it has a bell shape that is prescribed by the potential field topology. Where the iso-lines are the closest, the electric field is the highest, and Joule heating J • E is the most intense. The plume is mostly cold with a weak electric field.

In Figure 6 (f), the azimuthal electric field E y exhibits two distinct behaviors. In the channel, the azimuthal electric shows fluctuations in the axial direction with a wavelength of the order of ∼ 800 µm, which was also observed in the last row of Figure 4. In the plume, the azimuthal electric field also fluctuates in the radial direction, with a pattern resembling spherical-like waves, which seem to originate right downstream of the channel exit. Further studies, focused on the change of the radial gap, which goes from L z to 2L z in the plume, the decrease of the radial magnetic field B z , would provide a greater understanding of this pattern since the EDI is sensitive to these parameters as it will be shown in Section IV B. This seems to be a robust feature occurring in Hall Thrusters as recent experimental investigations found similar waves in the radial direction [START_REF] Tsikata | Coexistence of small-scale instabilities in Hall thrusters and hollow cathodes[END_REF] .

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.0133963 B. Spectral analysis

Dynamic Mode Decomposition

In Section IV A, the discussion focused on the radial and axial directions of the HT but the azimuthal direction is crucial as well and is best explored through spectral analysis. During the transient phase, detailed in Section III, several instabilities arise with different wavelengths in different directions. In order to highlight the modes present in the system, the Python library Antares 60 is used to perform a Dynamic Mode Decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF] . The strength of such a technique is its ability to single out the modes present in the system and to reconstruct them individually. In Figure 8 the main reconstructed modes are presented for the azimuthal electric field E y , the radial electron current density J e,z and the axial electron current density J e,x .

For the azimuthal electric field E y , in Figure 8 (a), a main short-wavelength azimuthal instability propagating in the whole domain is observed. The 3D view allows to show this azimuthal instability in the axial-azimuthal and radial-azimuthal planes, 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.0133963 the latter being very similar to the views of Figures 4 and5. Furthermore, it can be seen that this main mode is also responsible

for the E y pattern, resembling spherical-like waves, observed in the axial-radial plane in the plume region (Figure 6 (f)). This dominant mode can be interpreted as the presence of the EDI-ion acoustic instability. A more theoretical argument for this will be provided in Section IV B 2.

For the electron currents J e,x and J e,z , the corresponding main mode exhibits a radial-azimuthal pattern as shown in Figure 8 (b) and (c). For J e,x , the observed instability is mainly azimuthal, with a wavenumber k y = 2π/L y , although a small radial component is also present in the radial-azimuthal plane. In Figure 8 (b), J e,x exhibits a radial pattern with a wavenumber k z = π/L z . Again, the features observed on the raw data in Figure 5 are retrieved and a clear radial-azimuthal fluctuation is captured. Further discussion on this matter is conducted in Section IV B 3.

Comparison with theoretical results: EDI-ion acoustic wave

In order to consolidate the present analysis, the numerical results are now confronted to theory.

Nonlinear effects play an important role in the evolution of plasma instabilities [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF][START_REF] Lampe | Nonlinear development of the beam-cyclotron instability[END_REF][START_REF] Lampe | Theory and simulation of the beam cyclotron instability[END_REF] and dramatically complicate their analysis. Therefore the analysis first focuses on the beginning of the simulation when fluctuations have a small amplitudes and for which the linear theory of instabilities still applies. In this condition, the 3D dispersion relation (DR) derived by Ducrocq et al. [START_REF] Ducrocq | High-frequency electron drift instability in the cross-field configuration of Hall thrusters[END_REF] maybe apploied. It is recalled in Equation ( 14):

1+k 2 λ 2 D + g ω -k y V d Ω ce , k 2 x + k 2 y r 2 L , k 2 z ρ 2 - k 2 λ 2 D ω 2 pi (ω -k x v p ) 2 = 0, (14)
where ω pi is the ion plasma frequency, Ω ce the cyclotron frequency, V d = E x /B y the electron drift velocity, r L the Larmor radius of electrons, v p the velocity of the ion beam and k x , k y and k z the axial, azimuthal and radial components of the wave number k.

The Gordeev function 65 g(Ω, X,Y ) can be expressed in two ways:

g(Ω, X,Y ) = iΩ +∞ 0 e -X[1-cos(ϕ)]-1 2 Y ϕ 2 +iΩϕ dϕ, = Ω √ 2Y e -X +∞ ∑ m=-∞ Z Ω -m √ 2Y I m (X), (15)
where I m is the modified Bessel function of the first kind and Z the plasma dispersion function.

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. In contrast to earlier 2D PIC studies [START_REF] Janhunen | Nonlinear structures and anomalous transport in partially magnetized ExB plasmas[END_REF][START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF][START_REF] Petronio | Conditions of appearance and dynamics of the modified two-stream instability in ExB discharges[END_REF][START_REF] Hara | Cross-field electron diffusion due to the coupling of drift-driven microinstabilities[END_REF] , the DR does not need to be approximated and can fully take leverage of the present 3D PIC simulation. There is no analytical solution of Equation ( 14), but it can be numerically solved using a fixed point scheme described by Cavalier et al. [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] . A DR solver was then developed using Python [START_REF] Oliphant | Python for scientific computing[END_REF] and the built-in Fadeeva function from

PLEASE CITE THIS ARTICLE AS

Scipy. The Plasmpy package [START_REF] Community | The development of PlasmaPy has been partially supported by the U[END_REF] directly provides the plasma function Z as well. For given plasma parameters and wave numbers k x , k y and k z , the scheme converges to a complex number Ω = ω + iγ where ω and γ are the theoretical angular frequency and growth rate of the plasma wave. Input parameters are normalized with the ion plasma frequency ω pi , the Debye length λ D and the ion sound speed c s = (k B T e /m i ) -1/2 depending on whether they express a temporal frequency, spatial length or speed.

The temporal evolution of the spectrum of the dimensionless azimuthal wave number k y λ D is first presented in Figure 9 (a).

Here λ D is calculated at each instant as an average over the radial azimuthal zy plane at x = 0.12L x without the sheath that was previously measured to be around 1 mm. It appears that k y λ D remains fairly constant over time except at the beginning of the simulation where k y λ D ≈ 0.5, which corresponds to a wave number k y ≈ 6283 rad m -1 .

It can be observed in Figure 9 (b) how this single mode k y ≈ 6283 rad m -1 develops in the first moments of the simulation (the dark green solid line corresponds to the plane x = 0.12L x ). At this frequency the plasma instability quickly undergoes an exponential growth, with ∂ t E 2 y ≈ 2γE 2 y , where γ is the growth rate of the instability driven by the azimuthal electric field. With a least-square interpolation, the exponential factor 2γ is evaluated and compared with the expected theoretical growth rate in this radial-azimuthal plane, obtained in plasma conditions corresponding to the beginning of the exponential growth and reported in Table I.

The plasma parameters were averaged in the x = 0.12L x radial-azimuthal plane without the sheath. The azimuthal wave number spans the range 0 -5/λ D , where λ D is computed from the aforementioned plasma parameters. k x is set to 0, assuming that the instability does not propagate in the axial direction, in the linear stage at least. This is supported by Figure 4 (bottom row, first column), where it seems that the instability is indeed purely azimuthal. Furthermore, the effect of k x turns out to be neutral on the growth rate γ because k x is only involved in the k x v p term in Equation ( 14). The same holds true for v p . As shown by Cavalier et al. [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] , the final complex angular frequency ω can be written as ω = k x v p + κ where κ is some complex number with a non-negative imaginary part whereas k x v p ∈ R is a pure real number. Therefore assuming k x = 0 or even v p = 0 does not matter for the growth rate of the instability. Finally, the radial wave number k z is set to π/L z , because plasma sheath effects make possible the existence of shorter wave numbers below the geometrical constraint k z = 2π/L z 69,70 . This choice is consistent with the current simulation, as the radial-azimuthal structure shown in Figure 8 (c) presents a half wavelength fit between the walls in the channel. It should be noted that for this particular case, other values for k z (k z = 2πL z for example) are possible, but it was observed that the resulting growth rate was almost insensitive to this choice.

The resulting theoretical curve is shown in Figure 9 (c) together with the numerical growth rate γ, labeled as "PIC", previously computed in Figure 9 (b). An excellent agreement is observed with the linear theory, with a maximum growth rate reached for k y λ D = 0.5. Note that the DR curve presents oscillations around a mean value. This effect was predicted by Cavalier et al.'s [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] parametric study as the DR can transition from a discrete set of resonances, the growth rate being close to zero between each cyclotron resonance, to a continuous function that can be treated as an ion-acoustic DR. The corresponding continuous ionacoustic DR for the growth rate, derived by Lafleur et al. [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] , is also shown in Figure 9 (c) and its analytical expression, for the 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. given conditions, is reminded in Equation ( 16):

PLEASE CITE THIS ARTICLE AS

348 γ = πm e 8m i k y V d (1 + (k y λ D ) 2 ) 3/2 ( 16
)
The difference between the ion-acoustic DR (orange curve) and the actual DR (blue curve) indicates that the transition to an ion-acoustic mode is not fully achieved yet. The actual DR curve exhibiting a single resonance with small remaining oscillations, rather suggests that the simulation is already halfway between a pure discrete nature of the DR and a fully continuous ionacoustic DR. Besides, a close look at early times in Figure 9 (a) shows that the main azimuthal wave number quickly increases from k y λ D = 0.5 to higher values, close to the value of maximum growth rate for Equation ( 16), k y,max λ D = √ 2 -1 , highlighted in Figure 9 (a) with the red line. This is in agreement with a complete transition toward the ion acoustic DR after the first microseconds.

The previous analysis can be extended to other radial-azimuthal planes in the channel to verify how well the numerical results agree with linear theory. To do so, the most dominant mode in the first instants of the simulation is selected for each considered axial position and the corresponding growth rate is computed with a least-square interpolation as in Figure 9 (b). At the same time, the theoretical DR is also computed using local plasma parameters and its maximum growth rate is extracted. In all cases, the maximum growth rate is obtained approximately at the same dimensionless k y λ D wave number. The final axial profile of the growth rate is shown in Figure 9 (d). Overall, the numerical results are in excellent agreement with the linear theory, confirming that the dominant azimuthal mode is the EDI. However, the first growth rate obtained by the PIC simulation at x = 0.04L x is about two times lower than expected. This might be due to a too low temporal resolution, which makes the linear regression difficult to conduct as only a handful of points are available at this axial position. Besides, the position x = 0.04L x = 1 mm is very close to the anode sheath, which is likely to create plasma conditions violating the assumptions used to derive the DR from linear theory. The other axial positions show that the growth rate γ decreases with x, which suggests that the EDI grows first in the upstream section of the channel. This result is consistent with what was observed in the 2D maps of azimuthal electric field E y during the transient phase (see Figure 4).

Radial azimuthal mode

In addition to the EDI evidenced in the previous section, Figure 9 (c) suggests that only a single Bernstein mode is excited in the channel and thus no other modes should be visible. This is not contradictory with the radial azimuthal mode observed in Figure 8 (b) and (c), as the latter remains quite modest when considering the corresponding azimuthal electric field E y . Indeed, measured amplitude for this mode (not shown here) is around 200 V m -1 vs 30 kV m -1 for the main EDI mode. At this point, it is unclear why the radial azimuthal pattern is more visible in the current density fields J e,x and J e,z . As a reminder previous purely 2D radial-azimuthal cases, performed by Janhunen et al. [START_REF] Janhunen | Evolution of the electron cyclotron drift instability in twodimensions[END_REF] and Villafana et al. [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] , showed two modes: the MTSI and the EDI. Although it is premature to claim that the observed radial-azimuthal mode in the present 3D simulation is the MTSI, Figure 10 suggests that the similarity with the MTSI is nevertheless important. In this figure, the radial azimuthal mode from the 3D simulation is reported, as well as the identified MTSI from Villafana et al. [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] in the dimensionless (k y λ D , k z λ D ) plane. In 2D, the MTSI maximum growth rate γ max in the (k y λ D , k z λ D ) plane can be approximated by a simple parabola equation derived by [START_REF] Petronio | Conditions of appearance and dynamics of the modified two-stream instability in ExB discharges[END_REF] :

k z λ D = m e E x eB 2 z λ D (k y λ D ) 2 . ( 17
)
For both the 2D and 3D simulations, the corresponding parabolas were computed based on local plasma conditions. For the 2D simulation, the couple (k y λ D , k z λ D ) is coherently found to be close to the maximum growth rate of the MTSI. For the 3D simulation, the radial-azimuthal mode is also very close to the 2D maximum growth rate, which may be interpreted as a footprint of the MTSI. Obviously, since the parabola equation was derived for a 2D configuration, no formal proof can be inferred at this point. However, this finding encourages to explore other initial plasma conditions in order to assess whether or not it is possible to trigger a MTSI in 3D.

V. COMPARISON WITH 2D SIMULATION

In this section, the results of an analogous 2D PIC simulation are reported and compared to the present 3D results in order to identify 3D effects.

A. Definition of the analogous 2D model

A rectangular domain is considered, similar to the centered 2D plane obtained from a cut at z = L z of the 3D case. The analogous 2D computational domain is shown in Figure 11. Most parameters defining the original 3D domain can be directly used for the 2D simulations and are reported in Table II. In particular, the axial and azimuthal lengths, L x and L y respectively, are the same. The boundary conditions in the axial and azimuthal directions are also identical: periodicity at y = 0 with y = L y and 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. [START_REF] Petronio | Conditions of appearance and dynamics of the modified two-stream instability in ExB discharges[END_REF] . The 2D result is taken from previously published work [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] . The 3D result is obtained from the radial azimuthal mode, i.e., k y = 2π/L y and "exit" conditions for the particles at x = 0 and x = L x . The potential is set to 200 V at the anode and 0 V at the cathode, where the quasi-neutral (QN) model is used. The magnetic field strictly follows Equation (11) as in the 3D case. The 2D mesh also has a similar cell size with ∆x = ∆y = 50 µm, whereas it was ∆l = 67 µm in the 3D case. Computational parameters and initialization conditions are the same. Finally, the 2D ionization source term profile S 2D is given by:

PLEASE CITE THIS ARTICLE AS

k z = π/L z , at x = 0.12L x . B x y x Bmax Cathode (QN) Γ ai Γ ae Γ cd E Γ be Γ bi                      L y L x Anode
if x ∈ [x 1 , x 2 ] S 2D (x) = S 0 cos x-xm x 2 -x 1 else S 2D (x) = 0. (18)
Note that this profile exactly corresponds to the 3D profile S 3D in the centered axial-azimuthal plane. Mathematically this results in writing: S 2D (x) = S 3D (x, z = z M ). One consequence of this choice is that the current densities J m,2D and J m,3D cannot be equal for both simulations and from Equation ( 6) it can inferred that:

J m,2D = e 2 π S 0 (x 2 -x 1 ) J m,3D = e 4 π 2 S 0 (x 2 -x 1 ) (z 2 -z 1 ) Lz ⇒ J m,2D J m,3D = π 2 L z z 2 -z 1 > 1. ( 19
)
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.0133963

Therefore, this geometrical effect results in actually setting J m,2D ≈ 730.5 A m -2 vs J m,3D ≈ 400 A m -2 .

Parameters

B. Transient and steady state

In Figure 12 are shown the mean ion density profiles n i extracted from the 3D and 2D configurations (blue curves). In the first instants, the density linearly increases with the same slope, which is linked to the ionization source term creating similar plasma conditions for both the 2D and 3D simulations. Nevertheless, from around 1 µs, n i continues to rise for the 2D case whereas it starts to saturate in the 3D simulation. This difference may be due to the presence of the walls in 3D which probably limits the maximum possible density while in 2D, the losses are mitigated. In the 2D case, the density reaches a maximum value n i = 4 × 10 17 m -3 at 2.5 µs. It then drops until 5 µs before slowly increasing again to stabilize at around 3 × 10 17 m -3 at 15 µs.

These strong variations are due to the initial plasma wave leaving the domain at the cathode. Therefore, this phenomenon is similar to what was described in the 3D case in Figures 3 and4. However, n i does not vary with the same magnitude compared to the 3D simulation. In this regard, wall losses probably avoid in 3D the strong initial overshoot, which helps the simulation to reach a steady state more quickly. As a result, the 3D case only needs 10 µs to stabilize vs 15 µs in 2D.

A similar trend regarding the total electron temperature T e can be laid out as well. Although similar heating is noted for both configurations initially, T e remains significantly higher in 2D than in 3D at steady state. This is again due to the absence of losses at the walls in the radial direction, which overall makes the plasma hotter and denser in 2D than in 3D.

At steady state, Figure 13 In the plume, the wave is perpendicular to the y direction for both cases. However, in 3D the wavelength is very similar to that of the channel, whereas in 2D the wavelength is twice the wavelength in the channel. Overall, note that the magnitude of oscillations is three times bigger in 2D than in 3D. This difference in magnitude could be related to local plasma conditions (n i , T e ) that reach greater levels in 2D as seen in Figure 12.

In contrast to the top row of Figure 4, no large azimuthal wavelength fluctuations can be observed in 2D (not shown here).

In the absence of the radial direction, this second mode, which was observed in 3D, does not appear. This aligns well with the analysis of Section IV B 3, where it was assumed that this mode might be the signature of the MTSI. Indeed, this latter requires to include the direction parallel to the magnetic field (purely radial in this case) to develop [START_REF] Janhunen | Evolution of the electron cyclotron drift instability in twodimensions[END_REF][START_REF] Mcbride | Theory and simulation of turbulent heating by the modified two-stream instability[END_REF] . 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.0133963

C. Discussion on the EDI and the subsequent effect on anomalous electron transport

In order to understand why the EDI first appears in the near plume region in 2D, the linear theory of instabilities may be used as done in Section IV B 2. Unfortunately, Equation ( 14) cannot be used because the Gordeev function detailed in Equation (15) does not converge for k z = 0, which is a condition imposed in 2D. However, the azimuthal instabilities at steady state can still be analyzed. Ducrocq et al. [START_REF] Ducrocq | High-frequency electron drift instability in the cross-field configuration of Hall thrusters[END_REF] showed that in 2D the unstable, i.e., growing modes were located at multiples of the cyclotron frequency:

k y ≈ m Ω ce V d , with m = 1, 2, ... (20)
2D axial-azimuthal simulations agree with this theoretical result but surprisingly show a transition toward an ion-acoustic wave in some cases [START_REF] Boeuf | Exb electron drift instability in Hall thrusters: Particle-in-cell simulations vs. theory[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] . Indeed, a rigorous study of the full 3D dispersion relation [START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] showed that this transition can occur only for a non-zero wave-vector parallel to the magnetic field (k = 0), which is not the case in 1D azimuthal and 2D azimuthal-axial configurations. The mechanism responsible for this transition in 2D is not clear but [START_REF] Lampe | Nonlinear development of the beam-cyclotron instability[END_REF][START_REF] Lampe | Theory and simulation of the beam cyclotron instability[END_REF] suggests that an anomalous waveparticle interaction called resonance broadening could "smear out" discrete frequencies from Equation (20). There is however no consensus on this matter yet. When the transition toward an ion acoustic wave occurs, Lafleur et al. [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] derived the corresponding dispersion relation as:

ω R ≈ k • u i ± kc s 1 + k 2 λ 2 d , (21)
where v i is the ion beam velocity, c s is the ion speed of sound (Equation ( 13)) and k the wave vector.

For the present work, the numerical dispersion relation ω(k) has been obtained inside the channel at x = 0.12L x using a normalized 2D FFT computed over the last microsecond of each simulation. The axial and radial wave numbers are assumed negligible with respect to the azimuthal one. In Figure 14, the normalized DR in its dimensionless form (λ D k y , ω/ω pi ) is presented. The analytical form ω R is also plotted, presuming k • v i ≈ 0, i.e., the instability essentially propagates in a perpendicular direction to the ion beam. The main mode in both the 2D and 3D simulations is well predicted by the modified ion-acoustic DR. In 2D, the results seem to be more continuous as the dominant modes closely follow the theoretical DR whereas in 3D a single-mode appears to prevail. Note that in both simulations, the transition to the ion acoustic wave leads to different main modes: (k y,max λ D = 0.4, ω max /ω pi = 0.4) for the 3D case vs (k y,max λ D = 1.02, ω max /ω pi = 0.71) in 2D.

In recent studies [START_REF] Lafleur | Anomalous electron transport in hall-effect thrusters: Comparison between quasi-linear kinetic theory and particle-in-cell simulations[END_REF][START_REF] Charoy | Numerical study of electron transport in Hall Thrusters[END_REF][START_REF] Lafleur | The role of instability-enhanced friction on 'anomalous' electron and ion transport in hall-effect thrusters[END_REF] , the mobility of electrons in the axial direction was shown to be enhanced by the EDI. Thus, as the EDI differs in 2D and in 3D, one can expect to observe differences in electron mobility. In the present collisionless plasma, the mobility µ e f f is predicted to be solely due to the EDI. An analytical formula was given by Lafleur et al. [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF] :

µ e f f = - δ n e δ E y n e E x B z , (22)
where δ n e δ E y represents coherent fluctuations of the electron density and azimuthal electric field. As a first approximation, one can express these fluctuations around a certain time average operator denoted by . . . t :

δ n e δ E y = [ n e t -n e (t)] [ E t -E(t)] t (23)
The effective mobility stemming from the EDI can be compared to actual PIC measurements computed as:

µ PIC = - v e,x E x (24)
In Figure 15, both the expected µ e f f and measured µ PIC mobilities are compared. The computation was performed by using physical quantities averaged azimuthally along the axial direction. Results suffer from statistical noise as the electric field E x can reach very low values, especially in the plume and near the anode. Nevertheless it is observed that the mobility is ∼ 100 times higher in 2D than in 3D in the plume area. Near the anode the ratio µ PIC,2D /µ PIC,3D falls at 10 -1. Besides, the measured mobility µ PIC is minimum at the maximum E × B current, which occurs at distinct axial locations: at x ≈ 0.3L x in 2D versus

x ≈ 0.5L x in the 3D case. In both cases, the PIC mobility is well predicted by the expected anomalous transport due to the EDI. After examination of the µ e f f formula the coherent azimuthal fluctuations δ n e δ E y are found to be the primary factor explaining the overall greater mobility in 2D than in 3D. According to Lafleur et al. [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. i. insights from particle-in-cell simulations[END_REF] , an important saturation mechanism for the EDI is ion trapping. In a companion paper, Lafleur et al. [START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] provide an estimate of the amplitude of electric field fluctuations δ E y that are proportional to √ n e T e . From Figure 12, both the plasma density and electron temperature are approximately 2.5 times higher in 2D than in 3D, which is the correct order of magnitude measured in Figure 13. Besides, it also means that lower losses probably lead the EDI to saturate at a higher energy level in 2D than in 3D, which eventually is responsible for an enhanced mobility.

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

VI. STATISTICAL CONVERGENCE

PIC simulations are inherently noisy and a sufficient number of macroparticles is needed to avoid numerical artifacts [START_REF] Okuda | Verification of theory for plasma of finite-size particles[END_REF]73 .

In 2D, a few hundreds of particles per cell at steady state are needed to guarantee a statistical convergence [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF][START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] and N ppc,ini = 75 is enough in the considered 2D axial azimuthal case.

In 3D we tested three initial numbers of macroparticles of X + e /e -per cell N ppc,ini . Averaged macroscopic quantities such as the current or plasma density did not have important variations, indicating convergence was already sufficient. However, for spectral quantities, we note important differences. In Figure 16 the DMD spectrum for N ppc,ini = 30, 60 and 120 retrieves the same main EDI-ion acoustic mode. However for N ppc,ini = 30 and N ppc,ini = 60 we observe two higher frequency modes that disappear when we increase N ppc,ini up to 120. This suggests these two higher frequency modes are definitely not physical but rather numerical artifacts. These results confirm that at least N ppc,ini = 120 is needed for this 3D simulation.

VII. CONCLUSION

In this paper, a 3D PIC simulation of the channel and near plume region of a Hall Thruster was presented and analyzed. The setup is inspired by previous successful 2D PIC simulations by Charoy et al. [START_REF] Charoy | 2d axial-azimuthal particle-in-cell benchmark for low-temperature partially magnetized plasmas[END_REF] and Villafana et al. [START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] . In this 3D simulation, 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. 22)) and measured µ PIC (from Equation ( 24)) mobilities in absolute value. the plasma remains collisionless, curvature effects are neglected, dielectric effects are discarded and the magnetic field is purely radial.

PLEASE CITE THIS ARTICLE AS

During the transient state, the growth and development of the 3D structure of the EDI are in agreement with the linear theory of plasma instabilities. The EDI appears near the anode and the channel region before propagating into the plume. At steady state, the EDI presents non-zero wave numbers in all directions. Besides, a radial-azimuthal mode, characterized by small radial and azimuthal wave numbers k z = π/L z and k y = 2π/L y , is also noted. A preliminary analysis suggests that this mode might be the signature of the MTSI, previously observed in 2D [START_REF] Janhunen | Evolution of the electron cyclotron drift instability in twodimensions[END_REF][START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] and 3D [START_REF] Taccogna | Three-dimensional particle-in-cell model of Hall thruster: The discharge channel[END_REF] .

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.0133963

To explore 3D effects, a 2D analogous case was considered in the axial-azimuthal plane, to be compared with the centerline slice of the 3D case at z = L z . Same initial conditions, boundary conditions, and imposed electromagnetic fields as in the 3D case were used. In addition the imposed constant ionization source term helped to easily ensure the identical creation of charged particles in the domain. In both 2D and 3D cases, the EDI develops and eventually transitions to an ion acoustic wave, with amplitude fluctuations however three times higher in 2D than in 3D. The signature of the radial-azimuthal mode observed in 3D is absent in 2D

Still, the structure of the EDI is different in 2D and 3D as a higher axial wave number in the channel is found in 2D, while the azimuthal fluctuations in the plume have a larger wavelength. Larger amplitudes are obtained because the plasma is overall 50% hotter and denser in 2D due to lower plasma losses. Indeed, the walls and the plume expansion in the radial direction are absent in the 2D setup. Furthermore, the axial mobility in 2D is overall higher than in 3D. Following Lafleur et al.'s model, the axial electron mobility in a collisionless plasma is proportional to coherent azimuthal fluctuations of the plasma density and azimuthal electric field δ n e δ E y . A simple approximation of δ n e δ E y leads to a reasonable prediction of the axial mobility. This result is another evidence, in 2D and 3D, of the role of the EDI in the electron anomalous transport in a E × B discharge.

Beyond these findings on the physics, this work also demonstrates the capability of the PIC method to achieve highperformance simulations in 3D using an unstructured grid, a feature that is crucial for more realistic geometries. The subsequent collected 3D PIC data can serve as databases for hybrid or fluid models to improve and implement existing 3D effects. n e T e

Figure 1 .

 1 Figure 1. (a): Schematic view of the 3D computational domain in the radial-axial plane. The azimuthal direction has a length L y and is periodic. (b): Full 3D domain.

Figure 3 .

 3 Overall, the simulation reaches a steady state in approximately 10 µs. Besides, the transient stage can be divided into three main phases. In the first moments of the simulation, called phase I, the global density increases linearly thanks to the ionization source term because particle losses at the walls and in the plume remain limited. This situation goes on until ions traveling from the ionization zone reach the cathode boundary. The corresponding ion convective time t c can be roughly estimated as t c ∼ d c 2m i /(eE acc ) with E acc the typical electric field in the axial direction, i.e., E acc ∼ 50 kV m -1 , and d c the distance from the ionization zone to the cathode, i.e., d c ∼ 1 cm. This gives t c ∼ 0.7 µs, which approximately corresponds to the time for reaching the maximum density, at around 1.2 µs, as shown in Figure3. At this point, the total number of electron and ion macroparticles is around 4.26 billions.

Figure 3 .

 3 Figure 3. Temporal evolution of the space-averaged ion density n i in the whole domain. The three different phases (I,II,III) are delimited with black dashed lines.

Figure 4 .

 4 Figure 4. Timeline of variables of interest in the centered axial-azimuthal xy plane at z = L z . Top row: axial electron current J e,x ; bottom row: azimuthal electric field E y . Each column corresponds to a specific time: t = 0.20, 1.20, 5.00 and 7.50 µs. The dashed black line represents the channel exit at x = 0.34L x = 0.85 cm. The proportions between the x and y directions are conserved.

Figure 5 .

 5 Figure 5. Timeline of variables of interest in the radial-azimuthal z-y plane at x = 0.12L x = 3 mm. From the first to the last row are successively presented the axial electron current J e,x , the radial electron current J e,z and the azimuthal electric field E y . Each column corresponds to a specific time: t = 0.20, 1.20, 5.00 and 7.50 µs.The proportions between the y and z directions are conserved.

220

  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 ASDOI: 10.1063/5.0133963 3D Particle-In-Cell study of the Electron Drift Instability in a Hall Thruster using unstructured grids 9 radial electron current J e,z and the azimuthal electric field E y are presented in the channel plane located at x = 0.12L x = 3 mm. As

Figure 6 .

 6 Figure 6. Radial-axial view at steady state, at t = 10 µs. (a): axial ion speed v i,x with iso-contour of Mach number M. (b): potential φ with streamlines of electric field E in the plane. (c): electron azimuthal current density J e,y . (d): axial electron current density J e,x with streamlines of electron velocity v e in the plane. (e): total electron temperature T e with iso-contour of potential φ expressed in Volt. (f): azimuthal electric field E y .

Figure 7 .

 7 Figure 7. Sheath structure at the channel and exterior wall in the radial-axial plane at steady state, at t = 10 µs. This corresponds to the top left corner of the simulation domain, the shaded grey areas representing the walls of the thruster. (a): axial ion speed v i,x with iso-contour of Mach number M and streamlines of ion velocity v i . (b): Potential φ with iso-contours.

Figure 8 .

 8 Figure 8. 3D clip view of the main DMD reconstructed modes for (a) the azimuthal electric field E y , (b) the axial electron current density J e,x and (c) the radial electron current density J e,z . Shaded grey areas correspond to the part of the domain that is not plotted here.

xFigure 9 .

 9 Figure 9. Growth and temporal evolution of the EDI-ion acoustic instability. (a): temporal evolution of the dimensionless azimuthal wavenumber k y λ D in the radial-azimuthal plane at x = 0.12L x . The solid red line corresponds to the k y λ D for which the maximum growth rate is achieved for the ion-acoustic wave[START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] . (b): temporal evolution of the azimuthal electric field E y corresponding to the dominant azimuthal wavenumber k y = 6283 rad m -1 at different axial locations expressed as a fraction of the axial length L x (from x = 0.04L x to x = 0.32L x ). The exit channel is located at x = 0.34L x . The dotted lines indicate the linear regression obtained with a least square method during the exponential growth of the instability. (c): comparison of the dimensionless growth rate γ/ω pi (labeled as "PIC") measured in the x = 0.12L x plane with linear theory from Cavalier et al.[START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] (blue line) and with asymptotic ion acoustic DR from Lafleur et al.[START_REF] Lafleur | Theory for the anomalous electron transport in Hall effect thrusters. ii. kinetic model[END_REF] (orange line). (d) Axial profile of the growth rate in its dimensionless form γ/ω pi (left axis, blue) and dimensional form γ (right axis, green). Results are compared with predicted dimensionless values stemming from linear theory by Cavalier et al.[START_REF] Cavalier | Hall thruster plasma fluctuations identified as the exb electron drift instability: Modeling and fitting on experimental data[END_REF] (dashed line).

Figure 10 .

 10 Figure 10. Comparison of 2D and 3D PIC results with the theoretical maximum growth rate for the MTSI in 2D, derived by Petronio et al.[START_REF] Petronio | Conditions of appearance and dynamics of the modified two-stream instability in ExB discharges[END_REF] . The 2D result is taken from previously published work[START_REF] Villafana | 2d radial-azimuthal particle-in-cell benchmark for ExB discharges[END_REF] . The 3D result is obtained from the radial azimuthal mode, i.e., k y = 2π/L y and k z = π/L z , at x = 0.12L x .

Figure 11 .

 11 Figure 11. 2D axial azimuthal computational domain.

Figure 12 .

 12 Figure 12. Temporal evolution of space averaged ion density n i (in blue) and space averaged total electron temperature T e (in orange). Results of the 2D simulation are depicted with dashed lines whereas 3D results at z = L z are represented by solid lines.

Figure 13 .

 13 Figure 13. Comparison of the azimuthal electric field E y at steady state in (a) 3D and (b) 2D. In (a), E y values have been multiplied by 3 to fit on the same scale as (b).

Figure 14 .

 14 Figure 14. Numerical dispersion relation in 3D (top row) and 2D (bottom row) in the channel at x = 0.12L x . The red dots correspond to the dominant resonance. The modified ion acoustic dispersion relation ω R (k) is indicated by the white solid line.

4 µ [m 2 • V - 1 • s - 1 ]Figure 15 .

 421115 Figure 15. Comparison of expected µe f f (from Equation (22)) and measured µ PIC (from Equation (24)) mobilities in absolute value.

Figure 16 .

 16 Figure 16. DMD spectrum of azimuthal electric field E y obtained at steady state for different number of macroparticles N ppc,ini = 30, 60, 120 at initialization. The red rectangle highlights higher frequency resonances that disappear as we increase the number of particles per cell.

Furtherx{

  work based on the present 3D results may include consideration of collisions with neutral particles, similarly to what was done in 2D by Lafleur and Chabert 19 . A more realistic 2D magnetic field topology should be considered as the direction and strength of the E × B drift plays a crucial role in the EDI growth and development. 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= 0.12L x (|E y | 2 )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

DOI: 10.1063/5.0133963

	Parameters	Symbol Value	Unit
	Plasma parameters
	Ion density Electron temperature	n i T e	1.05 × 10 17 12.5	m -3 eV
	Ion beam velocity	v p	0	m s -1
	Axial electric field Radial magnetic field	E x B z	25.6 × 10 3 82.2	V m -1 G
		Wave numbers
	Axial wave number	k x	0	rad m -1
	Azimuthal wave number k y Radial wave number k z	rad m -1 π/L z = 3.14 × 10 2 rad m -1 0 -5/λ D

Table I .

 I Input parameters for the DR calculation measured in the radial-azimuthal plane at x = 0.12L x at 0.05 µs without sheaths.

Table II .

 II Simulations parameters for the 2D axial-azimuthal configuration.

		Symbol Value	Unit
	Simulation domain
	Cell size	∆x = ∆y 50	µm
	Domain size Number of triangular cells Number of nodes	L x × L y 2.5 × 1 N cell 500 × 200 × 2 N nodes 100,701	cm 3
		Initial conditions
	Plasma density Ion temperature	n 0 T i,0	1 × 10 17 0.5	m -3 eV
	Electron temperature	T e,0	10	eV
	Number of particles/cell	N ppc,ini 75
	Ionization module
	Axial limits	(x 1 , x 2 ) (0.25,1.00)	cm
	Electromagnetic fields
	Anode potential	φ a	200	V
	Cathode potential	φ c	0	V
	Maximum magnetic field	B max	100	G
	Position of maximum magnetic field x Bmax	0.75	cm
	Magnetic field at anode	B a	60	G
	Magnetic field at cathode	B c	10	G
	Computational parameters
	Time step	∆t	5	ps
	Average time range	N a	5000∆t	s
	Subcylcing frequency	f sub	5	iterations

A. Smolyakov, T. Zintel, L. Couedel, D. Sydorenko, A. Umnov, E. Sorokina, and N. Marusov, "Anomalous electron transport in one-dimensional electron cyclotron drift turbulence," Plasma Physics Reports 46, 496-505 (2020).

M. M. Turner, "Kinetic properties of particle-in-cell simulations compromised by monte carlo collisions," Physics of Plasmas 13, 033506 (2006), https://doi.org/10.1063/1.2169752.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.0133963

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.0133963

ACKNOWLEDGMENTS

The authors would like to warmly thank the Laboratoire de Physique des Plasmas (LPP), and in particular Thomas Charoy, Alejandro Alvarez Laguna, Federico Petronio, Pascal Chabert and Anne Bourdon for their advice and many useful discussions. W. Villafana acknowledges financial support from Safran Aircraft Engines doctoral research award as well as from the Association Nationale de la Recherche et de la Technologie (ANRT) as part of a CIFRE convention. This work has been partially funded by ANR (No. ANR-16-CHIN-003-01) and Safran Aircraft Engines within the project POSEIDON. This work used the HPC resources of IDRIS/TGCC/CINES under the allocation A0102A06074 made by GENCI.

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

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