Etienne Bertin

Elliot Brendel

Bruno Hérissé

Alexandre Chapoutot

Julien Alexandre Dit Sandretto

Prospects on the application of necessary optimality conditions on the resolution of the Goddard problem with unknown bounded parameters using interval arithmetics

Keywords: Intervals, Optimal control, Command, Ordinary differential equations, Validated simulation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Prospects on the application of necessary optimality conditions on the resolution of the Goddard problem with unknown bounded parameters using interval arithmetics

Etienne Bertin, Elliot Brendel, Bruno Hérissé, Alexandre Chapoutot, Julien

Alexandre Dit Sandretto

To cite this version:

Etienne Bertin, Elliot Brendel, Bruno Hérissé, Alexandre Chapoutot, Julien Alexandre Dit Sandretto. Prospects on the application of necessary optimality conditions on the resolution of the Goddard problem with unknown bounded parameters using interval arithmetics. Summer Workshop on Interval methods, Jul 2019, Palaiseau, France. hal-02372308

Introduction

Optimal control of aerospace systems is performed by modelling the considered system by dynamics depending on multiple uncertain parameters (for example, aerodynamic coefficients and maximal thrust). Usually, the optimal control problem is solved for the nominal values of these parameters and the robustness of the solution is demonstrated by dispersing the parameters around nominal values with Monte Carlo simulations. In addition to parameter uncertainties, the problem-solving method often introduces numerical approximation (for example the numerical solver of the ordinary differential equation representing the dynamics of the system or the optimization algorithm solving the optimal control problem). Interval Arithmetics has shown its ability to address several control problems, providing validated solutions while dealing with method uncertainties (numerical approximations) as well as with model uncertainties (unknown parameters). The Pontryagin Maximum Principle [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF] (PMP) provides necessary optimality conditions for the resolution of optimal control problems by transforming an optimal * Corresponding author.

control problem into a zero-finding problem : the dynamics of the system is extended with a co-state and necessary conditions are given by the PMP on that co-state, and the initial value of the co-state vector is the unknown to be found by the zero-finding algorithm. This method has proven its efficiency and its precision compared to direct methods [START_REF] Trélat | Optimal control and applications to aerospace: some results and challenges[END_REF], but its convergence depends strongly on its initialization and a prior knowledge of the solution structure is needed.

Our goal is to address the return version of the Goddard problem, which consists in performing the landing of the first stage of a rocket while minimizing its fuel consumption, combining interval arithmetics and the necessary optimality conditions given by the application of the PMP. Although this goal has not been reached yet, this paper presents preliminaries results on simplified problems, exposes the challenges encountered and suggests further developments. The optimal reentry trajectory for the Goddard problem is presented in Figure 1 with the ballistic phases in blue, and in red the first boost for the inversion of the speed vector, the intermediate boost for the dynamic pressure constraint and the landing boost. The evolution of the dispersions along the complete trajectory (launch and landing) of a rocket are described in [START_REF] Blackmore | Autonomous precision landing of space rockets[END_REF].

As a first step, a very simplified version of the Goddard problem is studied, namely a double integrator where [k] = k, k is an interval parameter and u is the control in- put. This double integrator fits in with a Goddard problem without gravity and aerodynamic forces, and with a constant mass. The interval parameter [k] represents the uncertainty on the maximal thrust force. Hence, the optimal control problem is

min T 0 |u|dt,        ṙ = v, v = [k]u, r(0) = 0, v(0) = 0, r(T ) = r T , v(T ) is free, T is fixed.
Two ways of combining interval arithmetics with the PMP are investigated : an open loop approach providing an enclosure on the system trajectory which can be used to assess robustness (regarding the interval parameter [k]), and a closed loop approach providing a closer enclosure of the optimal trajectories and can be used to initialize a non-interval algorithm.

Once the double integrator optimal control problem is solved, multiple approaches such as validated continuation methods are considered to solve the Goddard problem using this simplified version.

Open-loop approach

In this approach, the goal is to find the smallest initial co-state interval that contains every admissible co-state by combining numerical integration tools such as DynIbex1 [START_REF] Dit Sandretto | Validated explicit and implicit Runge-Kutta methods[END_REF] and an algorithm to solve the zero-finding problem. Many algorithms are considered, for example Krawczyk method, forward-backward operators and branch algorithms. The enclosure of the solution of the zero-finding problem provides a validated initialization for the co-state vector. Issues like discontinuous control input, control saturation, pure state constraints and mixed constraints are to be studied in order to solve a practical optimal control problem like the Goddard problem.

Closed-loop approach

In this approach, dynamic programming is used to find a finer enclosure of the optimal trajectories : if the system measures its interval state vector at a certain time, a new optimal control problem is solved from this interval, providing a better enclosure of the solution. Due to its algorithmic cost, this approach is irrelevant for solving practical cases online, but it can improve the enclosure of the solution when it is applied offline and therefore provides useful information for the online guidance algorithm.

Figure 1 :

 1 Figure 1: Optimal trajectory for the re-entry Goddard problem

http://perso.ensta-paristech.fr/∼chapoutot/dynibex/


J P Boeuf

Micro instabilities and rotating spokes in the near-anode region of partially magnetized plasmas

Micro instabilities and rotating spokes in the near-anode region of partially magnetized plasmas

J. Boeuf

I. INTRODUCTION

Partially magnetized plasmas are low pressure plasmas with moderate magnetic fields, i.e. where electrons are strongly magnetized (electron Larmor radius much smaller than plasma dimensions) while ions are not. These plasmas have specific properties, different from fusion plasmas, but are also subject to instabilities and turbulence phenomena. Because of these instabilities, electron transport across the magnetic field is "anomalous", and, hence, cannot be described by the classical, collisional cross-field mobility. Thus, the development of predictive models of these discharges is very challenging. In this paper we focus on the development of instabilities in the near-anode region of magnetized plasmas.

In partially magnetized discharges such as Hall thrusters [START_REF] Goebel | Fundamentals of Electric Propulsion: Ion and Hall Thrusters[END_REF][START_REF] Boeuf | [END_REF] , magnetron discharges 3 , ion sources for neutral beam injection 4 , Penning discharges 5 , etc…, the region of strong ionization is separated from the anode by a region of relatively low magnetic field and low electron temperature that we call "near-anode region" in this paper. In a Hall thruster, the length of this region in on the order of 1 to 2 cm, the gas pressure is in the mtorr range and the average magnetic field is on the order of or less than 50 G. Electrons are magnetized and collisional but the Hall parameter (ratio of the electron cyclotron frequency to the electron collision frequency) is large and their residence time is considerably increased by the magnetic field in the absence of instabilities.

The effect of instabilities is to enhance the electron momentum exchange rate, and this can be described by a a Electronic mail: jpb@laplace.univ-tlse.fr larger effective electron collision frequency, or a lower effective Hall parameter. Ions are practically not magnetized and will be considered as collisionless in this paper.

For low magnetic fields in the near-anode region, the electrons are sufficiently mobile that the sheath is repelling. That is, the anode potential is below the plasma potential, as in an unmagnetized plasma. This situation corresponds to a negative anode sheath and electron transport across the magnetic field can be considered as classical, collisional. When the magnetic field in this region is increased, the electron mobility across the magnetic field decreases and the anode sheath must become positive to pull electrons to the anode. In this work we study the transition between negative and positive anode sheaths and its relation to the development of instabilities and coherent structures of different wavelengths in the near-anode region.

The applied electric field and the discharge current in magnetized discharge devices are perpendicular to the external magnetic field, leading to E B and P B (diamagnetic) electron drift (E is the electric field parallel to the discharge current, B is the external magnetic field, and P is the electron pressure). The design of Hall thrusters or magnetrons is such that the E B and P B electron drift currents are in the azimuthal direction of a cylindrical geometry (these plasma devices are sometimes called "closed drift" devices 6 ). Since ions are essentially not magnetized, the azimuthal velocity of electrons is much larger than that of ions, and this difference in velocity can generate charge separation and trigger numerous types of instabilities [7][8][9][10] .

Instabilities in a large range of frequency and wavelength, from large scale "rotating spokes" to submillimeter structures have been reported experimentally or by simulations in the partially magnetized plasmas of Hall thrusters 7,8,[11][12][13][14][15][16][17][18][19][20][21][22][START_REF] Frias | 34th International Electric Propulsion Conference, IEPC-2015-370[END_REF][START_REF] Romadanov | [END_REF] , magnetrons 8,[25][26][27][28][29] , ion sources for neutral beam injection 30 , cylindrical magnetized plasma column [31][32][33][34][35][36] and Penning ion sources [37][38][39] . These instabilities are very dependent on the specific conditions or regions of a particular device. For example, in a Hall thruster, the large electric field that is present in the acceleration region (quasineutral, collisionless plasma) is responsible for the development of the Electron Cyclotron Drift Instability (ECDI) or E B Electron Drift Instability which has been the subject of many recent studies 7,8,18,19,[40][41][42][43] . This instability leads to anomalous transport but the effective Hall parameter associated with it is relatively large (larger than 100 according to the PIC-MCC simulations 43 ). Thus the effective cross-field electron mobility is still relatively small in this region. In this paper we consider a different situation, where the electric field is small but the combination of this low field with plasma density gradients and electron collisions can trigger instabilities. An important feature of the near-anode region when the sheath tends to be positive is that the plasma density gradient and the electric field are parallel and in the same direction. This situation is characteristic of the Simon-Hoh instability 44, 45 and we will see in this paper that the effective Hall parameter associated with this instability can be much smaller than the effective Hall parameter due to the ECDI in the acceleration region of Hall thrusters (i.e. larger effective collision frequency).

In this first approach we consider the simplifying assumption of a constant magnetic field in the near-anode region. Magnetic field gradients can modify the development of these instabilities 10,20,[46][47][48][49] and we leave the study of their effect for future work.

In section II we define the near-anode region and briefly summarize previous works related to this region and to the transition between positive and negative anode sheath. A simple relation between averaged magnetic field, electron temperature, length of the near-anode region and ion mass, that characterizes this transition is proposed in section III. The dispersion relation for the development of instabilities is discussed ion section IV. Finally, the formation of instabilities and coherent structures in the near-anode region, as predicted by PIC-MCC simulations, is described in section V and analyzed with the help of the dispersion relation. A general discussion of the results is provided in section VI.

II. THE NEAR-ANODE REGION

In this paper, we focus on the "near-anode region" of magnetized discharges. We define this region as the region between the anode and the region of large ionization. The near-anode region, shown in Figure 1 in grey color, is characterized by a low electric field, a relatively low magnetic field and a plasma density gradient. We will assume here that ionization is negligible in this region. In a Hall thruster (Figure 1a) this large field region, called the acceleration region, is located in the quasineutral plasma. The ion beam is extracted on the right of the acceleration region, and an emissive cathode allows neutralization of the ion beam as well as providing a source to sustain the discharge. In a magnetron discharge (Figure 1b) the large field region is the ion sheath next to the cathode. In both discharges, electrons accelerated in the large field region release their energy and ionize the gas in the ionization region. The thickness of the ionization region is relatively small because of the large magnetic field and is on the order of 1 cm in the conditions of Hall thrusters [START_REF] Goebel | Fundamentals of Electric Propulsion: Ion and Hall Thrusters[END_REF][START_REF] Boeuf | [END_REF] .

In "modern" Hall thrusters the radial magnetic field is maximum around the channel exhaust plane (Figure 1a) and the near-anode region is relatively short (less than 2 cm). Proper operation of a Hall thruster is achieved when most of the ionization takes place immediately upstream of the acceleration region and we will assume here, as said above, that ionization is negligible in the near-anode region. This supposes a significant decrease of the electron temperature over a short distance due to collisions with neutral atoms and electron-wall interaction. In the early experiments of Janes and Lowder 11 , the thruster channel was long (10 cm), with a large radial magnetic field along the entire channel length. In that case, ionization must be present in the plasma up to the anode to ensure quasineutrality and can be the result of electron heating associated with the formation of instabilities. Janes and Lowder 11 showed the presence of rotating instabilities "rotating spokes") in their experiments and suggested that the mechanism of Critical Ionization Velocity (CIV) was responsible for this plasma rotation and for ionization in the bulk plasma. The CIV concept, introduced by Alfven for space plasmas [START_REF] Alfven | On the Origin of the Solar System[END_REF] has been thoroughly discussed in the context of magnetron discharges 8,17,[START_REF] Brenning | [END_REF][52][53][54] .

In the conditions we are considering here, electrons must flow from the ionization region to the anode through a small but non-negligible magnetic field. A small ion flux (ion backflow) must also cross this region and reach the anode to ensure quasineutrality. Depending on the value of the magnetic field the plasma potential in the near-anode region may be above the anode potential (negative anode sheath, as in Figure 1) or below the anode potential (positive sheath). Ionization may or may not be negligible in this region in a real device. A positive anode sheath will tend to form for a large enough transverse magnetic field due to the small crossfield electron mobility. More generally, a positive anode sheath can be due to any mechanism that impedes electron flow to the anode (geometry, small anode or anode partially covered with dielectric as reported by Dorf et al. 55,56 ).

Note that in the case of Hall thrusters, the discharge current can exhibit strong oscillations in the 10-30 kHz range 57 (breathing oscillations). These oscillations are associated with a periodic depletion of the neutral atom density. Because of ionization, the front of the neutral flow moves upstream to a region where the electron mean energy is lower, and ionization stops until the neutral flow fills again the depleted region. Sekerak et al. 58 showed that this regime, termed as "global oscillation mode" is present when the magnetic field is below a given transition point. For magnetic fields above this transition point the discharge regime evolves into a "local oscillation mode" where the current oscillations due to the periodic axial neutral depletion disappear, while some rotating plasma structures appear in the azimuthal direction. The influence of ionization in the near-anode region on the formation of azimuthal instabilities in the local oscillation mode is not clear and we choose, in order to simplify the interpretation of the results in this first approach, to neglect ionization in that region. The influence of ionization in the near-anode region will be considered in future work.

In this paper we consider situations represented in Figure 1b, with a given, non-self-consistent ionization source. The electric field and ionization profiles on this figure are typical of magnetrons since the potential drop on the cathode side takes place in an ion sheath and not in the quasineutral plasma as in the case of Figure 1a. Since we are only interested here in plasma transport in the near-anode region, the general results and conclusions of this work can be applied to Hall thrusters as well as to magnetron discharges.

III. TRANSITION FROM NEGATIVE TO POSITIVE ANODE SHEATH

In this section we derive a simple criterion for the transition from a negative to a positive anode sheath, and we check the validity of this criterion with a 1D (onedimensional) PIC-MCC simulation.

We consider a simple, non-self-consistent 1D model of a magnetized discharge with the following assumptions (see Figure 2): -A given voltage 𝑈 𝐷 is applied between cathode and anode with a given, uniform magnetic field B perpendicular to the discharge current. -Ionization is not treated self-consistently, and the ionization source S is located near the cathode region, with a given profile (and associated injection current) as shown in Figure 2a. A constant electron temperature, 𝑇 𝑒 , is imposed in the ionization region. -Only elastic collisions between electrons and neutral atoms are considered, with a given, constant momentum exchange frequency  en . Ions are collisionless. This model will be used in the 1D PIC-MCC simulations described in the second part of this section and in the 2D PIC-MCC simulations of section V. We first derive a simple analytical criterion from the transition from a negative to a positive anode in these conditions.

A. Analytical model of the transition

In the absence of instabilities, the electron flux in the nearanode region at steady-state can be written in a drift diffusion form, accounting for the presence of the perpendicular magnetic field: Γ 𝑒 = 𝑛𝜇 𝑒,⊥ ∇ϕ -𝐷 𝑒,⊥ ∇𝑛 = 𝑛𝜇 𝑒,⊥ (∇ϕ -T e ∇𝑛 𝑛 ⁄ )

where ϕ is the electric potential, T e the electron temperature (T e is in units of volts; we note 𝑇 𝑒 the electron temperature in eV), 𝑛 is the plasma density, and 𝜇 𝑒,⊥ = 𝑒 𝑚𝜈 𝑒𝑛 1 1+ℎ 2 is the classical, collisional electron mobility along the discharge axis, i.e. perpendicular to the magnetic field. ℎ is the Hall parameter defined by: ℎ = 𝜔 𝑐𝑒 𝜈 𝑒𝑛 = 𝑒𝐵 𝑚𝜈 𝑒𝑛 𝐷 𝑒,⊥ = 𝜇 𝑒,⊥ 𝑇 𝑒 is the electron diffusion coefficient perpendicular to B. In this 1D model, ∇≡ 𝜕 𝜕 𝑥 ⁄ . We assume that the electron cyclotron frequency 𝜔 𝑐𝑒 is much larger than the electron-neutral collision frequency 𝜈 𝑒𝑛 , 𝜔 𝑐𝑒 ≫ 𝜈 𝑒𝑛 , so that the perpendicular electron mobility and diffusion coefficient can be written as:

𝜇 𝑒,⊥ = 1 𝐵ℎ , 𝐷 𝑒,⊥ = 𝜇 𝑒,⊥ T e =
T e 𝐵ℎ We define the transition from negative sheath to positive sheath as the "no-sheath" regime, corresponding to a situation where the drift term in the expression of the electron flux becomes small with respect to the diffusion term. Setting the axial electric field in the plasma to zero at this transition point, the electron flux to the anode writes: Γ 𝑒𝑎,𝑇 = -𝐷 𝑒,⊥ ∇𝑛 = -T e 𝐵ℎ ∇𝑛 ≈ -T e 𝐵ℎ 𝑛 0 𝐿 𝑛 where 𝑛 0 is the plasma density at the edge of the near-anode region (𝑥 = 𝑥 1 ) and 𝐿 𝑛 = 𝑛 0 ∇𝑛 ⁄ characterizes the plasma density gradient length (𝐿 𝑛 ≡ 𝑥 1 ). 𝑛 0 is close to the maximum plasma density in the conditions of this paper. The ion flux to the anode is small with respect to the electron flux so that the electron flux to the anode must be close to the ion flux to the cathode. In the case of the cathode sheath of a magnetron discharge the ion flux at the sheath edge can be approximated by: Γ 𝑖𝑐 = 𝛼𝑛 0 𝑐 𝑆 where 𝑐 𝑆 = √𝑒T e 𝑀 ⁄ is the ion acoustic velocity, 𝑀 is the ion mass, and 𝛼 is a parameter on the order of ½. Writing that these two fluxes are equal, |Γ 𝑒𝑎,𝑇 | ≈ |Γ 𝑖𝑐 | leads to the relation:

T e 𝐵ℎ𝐿 𝑛 = 𝛼𝑐 𝑆
This gives the following expression for the magnetic field at the transition point (i.e. the magnetic field above which the anode sheath becomes positive), as a function of electron temperature, ion mass, electron collision frequency and length of the density gradient. 3: Magnetic field at the negative-positive sheath transition point, as a function of ion mass for an electron temperature 𝑇 𝑒 =5eV, an electron-neutral collision frequency of 5 × 10 6 s -1 , and for two values of the density gradient length, with 𝛼 = 0.5. The colored full lines are from the analytical expression (see text) and the circle symbols are obtained from the PIC-MCC simulations. For each value of 𝐿 𝑛 the sheath is negative below the curves and become positive above. The same electron temperature is used in the nearanode region and in the Bohm velocity 𝑐 𝑆 on the cathode side.

𝐵 = (

T e 𝐿 𝑛 𝑚𝜈 𝑒𝑛 𝑒𝛼𝑐 𝑆 )

1/2

The magnetic field at the transition point is plotted in Figure 3 as a function of ion mass for 𝑇 𝑒 =5 eV and 𝜈 𝑒𝑛 = 5 × 10 6 s -1 , for 𝐿 𝑛 =1 cm and 𝐿 𝑛 =2cm ( is taken equal to ½).

In the case of a Hall thruster, the 𝛼 coefficient (the ratio of the plasma density at the sonic point to the maximum plasma density) cannot be easily estimated. The sonic point is located at the beginning of the acceleration region and the 𝛼 coefficient depends on the characteristics of a given thruster (magnetic field, anomalous transport, mass flow rate, wall losses, etc…). Note also that, since there is a large decrease of the electron temperature in the ionization region, the electron temperature that should be used in the Bohm velocity 𝑐 𝑆 for a real device is actually larger than the electron temperature 𝑇 𝑒 in the near-anode region.

We conclude that although the expression above is based on crude approximations, it provides a useful scaling of the magnetic field at the transition point, as a function of collision frequency, gradient length, and ion mass.

Finally, in the case of a non-uniform magnetic field in the near-anode region, the expression above for the magnetic field at the transition point can be modified to:

√ 𝐵 2 ̅̅̅̅ = ( 𝑇 𝑒 𝐿 𝑛 𝑚𝜈 𝑒𝑛 𝑒𝛼𝑐 𝑆 ) 1/2
where 𝐵 2 ̅̅̅̅ is the spatial average of 𝐵 2 in the near-anode region, i.e. over the length 𝐿 𝑛 .

B. 1D PIC-MCC simulation of the transition

A 1D PIC-MCC simulation has been used to check the scaling of the analytical relation derived above. The simulation is based on the assumptions defined at the beginning of this section. A voltage 𝑈 𝐷 = 200 𝑉 is applied between cathode and anode. Electrons are supposed to undergo only elastic collisions with a constant momentum transfer frequency 𝜈 𝑒𝑛 = 5 × 10 6 s -1 (corresponding to pressure on the order of 1 mtorr, typical of the near-anode region of Hall thrusters). Ions are collisionless. A constant magnetic field 𝐵 is applied perpendicular to the discharge axis. Electron-ion pairs are injected between x 1 and x 2 according to a cosine distribution, as in Ref. 43 :

𝑆(𝑥) = 𝑆 0 𝑐𝑜𝑠 (𝜋 𝑥-𝑥 𝑀 𝑥 2 -𝑥 1 ) for 𝑥 1 ≤ 𝑥 ≤ 𝑥 2 with 𝑥 𝑀 = 𝑥 2 +𝑥 1 2 𝑆(𝑥) = 0 for 𝑥 < 𝑥 1 or 𝑥 > 𝑥 2
In the simulations presented here the discharge gap d, and 𝑥 1 and 𝑥 2 are taken as: 𝑑 = 1.5 cm, 𝑥 1 = 0.7𝑑 = 1.05 cm, 𝑥 2 = 0.9𝑑 = 1.35 cm Electrons and ions are injected in the ionization region with Maxwellian velocity distribution functions with temperature 𝑇 𝑒 = 5 eV for electrons and 𝑇 𝑖 = 0.1 eV for ions. To maintain a constant electron temperature in the ionization region, the velocity of electrons in this region is changed randomly according to a Maxwellian distribution at a constant frequency 𝜈 𝑀𝑎𝑥 = 5 × 10 7 𝑠 -1 . These "fictitious" collisions affect electron transport in the magnetized plasma of the ionization region but we are mainly interested in electron and ion transport in the region between the ionization region and the anode and the detailed physics of the ionization region is out of the scope of this paper. The applied voltage is set to 𝑈 𝐷 = 200 𝑉 in all the simulations below. The value of this voltage has little effects on the results because the electron temperature is fixed independently (in contrast with experiments where the electron temperature is controlled by the applied voltage). The simulations were performed until steady state was reached, for increasing values of the magnetic field until a transition from negative to positive anode sheath was observed. To check the analytical relation derived above, the simulations were performed for different gas masses (keeping the same electron collision frequency) corresponding to helium, argon, and xenon (see the symbols in Figure 3). 𝑆 0 is chosen so that the maximum plasma density at steady state allows to run the simulation with 250 grid points (i.e. the electron Debye length is larger than d/250=60 m. 𝑆 0 defines an injection current density:

𝐽 𝑠 = 𝑒(|Γ 𝑒𝑎 | + |Γ 𝑒𝑐 |) = 𝑒(|Γ 𝑖𝑎 | + |Γ 𝑖𝑐 |) = 𝑒 ∫ 𝑆(𝑥)𝑑𝑥 𝑑 0 = 𝑒𝑆 0 2 𝜋 (𝑥 2 -𝑥 1 )
where Γ is a charged particle flux. The indexes e and i refer to electrons and ions, and the indexes a and c refer to anode and cathode. Since |Γ 𝑖𝑐 | ≫ |Γ 𝑖𝑎 | and Γ 𝑖𝑐 = 𝛼𝑛 0 𝑐 𝑆 we see that 𝑛 0 is proportional to 𝑆 0 /𝑐 𝑆 . Therefore, in order to keep similar values of the plasma density and number of grid points for the simulations of gases with different mass 𝑀, the value of 𝑆 0 was varied roughly as 𝑀 -1/2 , i.e. the discharge current density in the simulations was smaller for larger ion mass. The values of 𝑆 0 in the simulations were such that the total injection current density 𝐽 𝑠 was equal to 30 A/m 2 for helium, 10 A/m 2 for argon, and 3 A/m 2 for xenon. The number of particles per cell in these 1D simulations was on the order of 10 4 . The symbols in Figure 3 show the value of the magnetic field corresponding to the transition from negative to positive anode sheath as observed in the PIC-MCC simulations for the mass of helium, argon and xenon (everything else except 𝑆 0 being kept constant: electron temperature, discharge voltage and electron collision frequency). These points are relatively close to the curve obtained with the analytical expression of the transition magnetic field for 𝐿 𝑛 = 1 cm .

)

n e B=18 G (a) B=21 G (b) B=21.5 G  (V)
x (cm) B=18 G The transition from negative to positive sheath predicted by the PIC-MCC simulations is illustrated in Figure 4 where the axial variations of the electron and ion density profiles (Figure 4a) and of the electric potential in the near-anode region (Figure 4b) are displayed. The potential profiles indicate that the transition takes place for B slightly below 21 G in this example (helium). Note that the plasma density is larger for larger magnetic fields, as expected. For magnetic fields above the transition, the relation between potential and ion density becomes Boltzmann (ions are collisionless and are repelled by the anode sheath), and since the plasma density decay is almost linear, the potential varies logarithmically with position (see the effect of the ion temperature on the plasma potential in Figure 4b).

When the magnetic field is increased above the transition point, the 1D PIC simulations predict a significant increase of the anode sheath length with a large potential increase in the sheath. Electrons are extracted as in a Child Langmuir magnetized sheath. We will see in the sections below that for magnetic fields around or above the transition point the solution is actually not stable so that the 1D model is no longer meaningful. Instabilities and rotating structures are present around and above the transition magnetic field. The onset of instabilities and the formation of coherent structures in the direction perpendicular to the discharge axis and magnetic field are studied on the basis of dispersion relations in section IV and are analyzed with 2D PIC-MCC simulations in section V.

IV. DISPERSION RELATION IN THE NEAR-ANODE REGION

A very complete overview of the theory of instabilities in partially magnetized plasmas of E B discharges that includes the effects of density gradients, magnetic field gradients, electron gyroviscosity and collisions has been recently given by A. Smolyakov and his colleagues 9,10,20,[START_REF] Romadanov | [END_REF]48 .

We consider here a purely azimuthal instability (i.e. in the direction perpendicular to the cathode-anode discharge and to the magnetic field) and neglect the component parallel to the magnetic field and the axial component of the wave vector. The dispersion relation in that case, neglecting charge separation effects and without magnetic field gradient writes (formula (31) of Ref. 9 ):

𝜔 * + 𝑘 2 𝜌 𝑒 2 (𝜔 -𝜔 0 + 𝑖𝜈 𝑒𝑛 ) 𝜔 -𝜔 0 + 𝑘 2 𝜌 𝑒 2 (𝜔 -𝜔 0 + 𝑖𝜈 𝑒𝑛 ) = 𝑘 2 𝑐 𝑠 2 𝜔 2
where 𝑘 is the wave number in the azimuthal direction, 𝜌 𝑒 = (𝑒T e 𝑚 ⁄ ) 1/2 𝜔 𝑐𝑒 ⁄ is the electron Larmor radius, 𝜔 * = 𝑘𝑣 * , with 𝑣 * = T e (BL n ) ⁄ the diamagnetic drift velocity, 𝜔 0 = 𝑘𝑣 0 , with 𝑣 0 = 𝐸 0 𝐵 ⁄ the drift velocity, 𝑐 𝑠 = (𝑒T e 𝑀 ⁄ ) 1/2 is the ion sound speed and 𝜈 𝑒𝑛 is the electron-neutral collision frequency.

This dispersion relation can be considered as a high frequency or small wavelength extension of the Simon-Hoh instability 9 . The 𝑘 2 𝜌 𝑒 2 terms in this general dispersion relation comes from the inclusion of the electron inertia and gyroviscosity in the linearized equations 9 . In the limit of long wavelengths,i.e. when 𝑘𝜌 𝑒 ≪ 1, the dispersion relation takes the form of the classical, collisionless Simon-Hoh dispersion relation. In the limit of short wavelengths, 𝑘𝜌 𝑒 ≫ 1, the dispersion relation results in the ion sound mode, 𝜔 2 = 𝑘 2 𝑐 𝑠 2 that can be destabilized by E B drift or diamagnetic drift, and collisions. This equation can be written as follows, in normalized units:

Ω 3 + Ω 2 ( 1 Κ V * -Κ𝑉 0 + 𝑖Ν 𝑒 ) -Ω(1 + 𝐾 2 ) + Κ𝑉 0 (1 + 𝐾 2 )
-𝑖Ν 𝑒 𝐾 2 = 0

where Ω = The electric field and density gradient are in the same direction. The electron-neutral collision frequency is set to 𝜈 𝑒𝑛 = 5 × 10 6 s -1 (full lines) or to 10 6 𝑠 -1 (dotted lines, (a) only). The dashed line in (b) corresponds to the ion acoustic instability dispersion relation 𝜔 = 𝑘𝑐 𝑠 (i.e. 𝜔 𝜔 𝐿𝐻 ⁄ = 𝑘𝜌 𝑒 ). In normalized units, the three cases correspond to 𝑉 0 = 0, 6.5, and 26, with 𝑉 * =65. The dispersion relation Ω(𝐾) above has three independent normalized parameters, V * , 𝑉 0 and Ν 𝑒 respectively related to diamagnetic and E B drift normalized to the ion sound speed, and to the electron collision frequency normalized to the lower hybrid angular frequency.

E 0 =2 V/cm E 0 =0.5 V/cm  LH E 0 =0 (a) E 0 =2 V/cm E 0 =0.5 V/cm  LH k e E 0 =0 (b)
L n =5 cm L n =2 cm L n =5 cm  LH L n = (a) L n =2 cm L n =  LH k e (b)
The real and imaginary parts of the of the solutions of the dispersion relation give, respectively, the angular frequency 𝜔 and the growth rate 𝛾 of the instability. The plasma is not stable to perturbations for positive growth rates, i.e. 𝛾 > 0. Figure 5 shows three examples of dispersion curves Ω(𝐾). corresponding to xenon, with a magnetic field B=40 Gauss, 𝑇 𝑒 = 5 eV, L n = 1 cm ( T e /L n =5 V/cm), for three values of 𝐸 0 : 0, 0.5, and 2 V/cm, and with 𝜈 𝑒𝑛 = 5 × 10 6 s -1 . The growth rate of the instability exhibits some well-defined maxima that shift to lower wavenumbers, i.e. larger wavelengths as the plasma electric field 𝐸 0 increases. For large values of 𝑘𝜌 𝑒 the instability tends to the ion sound mode (dashed line of Figure 5b). For 𝐸 0 = 0 the instability corresponds to the destabilization of the ion-sound mode by density gradient and collisions. The growth rate of this mode increases with the collision frequency. For larger values of the plasma electric field, the growth rate is less sensitive to the collision frequency as can be seen in Figure 5a (compare the full and dotted lines for 𝐸 0 = 0.5 V/cm and 𝐸 0 = 2 V/cm), and the maximum growth rate decreases with increasing collision frequency. These modes are typical of the Simon-Hoh instability which occurs in plasmas where density gradient and electric field are in the same direction. Figure 6 shows dispersion curves for different values of the density gradient and for fixed values of 𝐸 0 , B, 𝑇 e , and 𝜈 𝑒 , respectively 2 V/cm, G, 5 eV, and 5× 10 6 s -1 , in xenon. We see that the growth rate is maximum for smaller wavenumbers, i.e. larger wavelengths, for decreasing plasma density gradients (increasing 𝐿 𝑛 ).

A situation symmetric to the zero field case, 𝐸 0 = 0, T e 𝐿 𝑛 ⁄ ≠ 0, and 𝜈 𝑒𝑛 ≠ 0 described above is the zerodensity gradient case (𝐿 𝑛 ∞ ), where 𝐸 0 ≠ 0, T e 𝐿 𝑛 ⁄ = 0 and 𝜈 𝑒𝑛 ≠ 0. This mode corresponds to a lower hybrid mode destabilized by E B drift and collisions (termed as resistive mode in Ref. 59 ). We see in Figure 6 that the mode frequency in the zero-density gradient case is close to the lower hybrid frequency at the maximum growth rate. In the limit 𝜔 < 𝜔 0 and 𝜈 𝑒𝑛 < 𝜔 0 the growth rate increases with the collision frequency as 9,59 : 7 presents a useful overview of the characteristics of the instability as a function of density gradient and plasma electric field. The maximum growth rate, the angular frequency, and the wavelength 𝜆 = 2𝜋/𝑘 of the instability at the maximum growth rate are plotted in normalized units in Figure 7a, b, andc ⁄ for a particular case (xenon, 𝐵 = 40 G and 𝑇 𝑒 = 5 eV, giving speed 𝑐 𝑠 = 1.9 × 10 3 m/s and 𝜔 𝐿𝐻 =1.43 × 10 6 rad/s) is also provided for clarity. We see that the growth rate is on the order or larger than the lower hybrid frequency in most of the considered domain of variations of 𝑉 0 and V * , i.e. several 10 6 s -1 for xenon at 40 G, and increases with increasing density gradient.

ω 𝜔 𝐿𝐻 = 1 + 𝑖 𝜈 𝑒𝑛 2𝜔 0 = 1 + 𝑖 1 2𝑘𝜌 𝑒 𝜈 𝑒𝑛 𝜔 𝐿𝐻 1 𝑉 0 Figure
The wavelength at maximum growth rate becomes very small around 𝑉 0 ⋍ 1 , i.e. when the 𝐸 0 𝐵 ⁄ drift velocity is small and on the order of the ion acoustic velocity (white region of the contour plot in Figure 7c) . A sharp increase of the angular frequency can be seen in this region (Figure 7b), corresponding to the ion sound mode where the angular frequency increases with increasing wave number (decreasing wavelength). In this region the dispersion relation should actually take into account charge separation (i.e. a 𝑘 2 𝜆 𝐷 2 term should be added in the left-hand side of the dispersion relation above, see, e.g. eq. ( 19) of Ref. 10 ) to properly describe the ion acoustic wave, and the wavelength of the instability at the maximum growth rate should depend on the Debye length for lower values of the plasma density.

The wavelength of the instability at maximum growth rate increases for increasing 𝑉 0 (i.e. plasma electric field) and decreasing 𝑉 * (density gradient). For large enough 𝑉 0 the wavelength at maximum growth rate increases sharply with decreasing 𝑉 * when the angular frequency becomes on the order of the lower hybrid frequency. The instability growth rate in this region also decreases with decreasing 𝑉 * and is below the lower hybrid frequency.

Finally, Figure 8 shows the velocity of the azimuthal wave at the maximum growth rate as a function of 𝑉 0 and 𝑉 * in the conditions of Figure 7. ⁄ in the same conditions as Figure 7.

The phase velocity of the wave is close to the ion sound velocity 𝑐 𝑠 in a large part of the domain, is smaller than 𝑐 𝑠 below 𝑉 0 = 1, and increases above 𝑐 𝑠 in the resistive, lower hybrid mode.

We summarize below the main conclusions of this section, which are valid for partially magnetized plasmas typical of the near-anode region of Hall thrusters or magnetrons under the assumptions of constant electron temperature, constant magnetic field, no ionization, and collisionless ions: -The development of instabilities is controlled by three main parameters related to the plasma density gradient, plasma electric field, and collision frequency, namely V * = -The Simon-Hoh instability plays a dominant role in a large part of the (V * , V 0 ) domain. -The maximum growth rate of the instability is on the order of the lower hybrid frequency and slightly increases with V * from 𝜔 𝐿𝐻 at V * ≈ 10 to 4𝜔 𝐿𝐻 at V * ≈ 60. The angular frequency at maximum growth rate varies from 𝜔 𝐿𝐻 at V * ≈ 10 to several 𝜔 𝐿𝐻 when V 0 decreases, with a sharp increase in the ion sound wave region around V 0 = 1. -The wavelength at maximum growth rate increases from values on the order of the electron Larmor radius 𝜌 𝑒 for V 0 around 1 to several 𝜌 𝑒 when V * decreases and/or V 0 increases. For low values of the density gradient, i.e. low V * , the wavelength increases sharply -The domains of destabilization by density gradient and collisions or by electric field and collisions are limited to small regions of the (V * , V 0 ) domain corresponding respectively to low values of V 0 (V 0 < 1) or low values of V * (V * ≾ 0.3 V 0 ). -The phase velocity of the wave at the maximum growth rate is close to or slightly above the ion acoustic velocity in most of the (V * , V 0 ) domain. In the results presented in this paper, the dimensions are taken as: 𝑑 = 1.5 cm, 𝑥 1 = 0.7𝑑 = 1.05 cm, 𝑥 2 = 0.9𝑑 = 1.35 cm. The results presented here have been obtained for a periodic length 𝑤 = 6 cm. The number of grid points is 128 × 512 and most of the simulations are performed with about 100 particles per cell.

V. 2D PIC MCC SIMULATION OF THE NEAR-ANODE REGION

In order to study the development of instabilities in the E direction (azimuthal direction of E plasma devices), 2D PIC-MCC simulations were performed in a 2D domain perpendicular to the constant magnetic field and defined by the axial and azimuthal directions, as in Ref. 43 . As in this reference and in the 1D model of section III.B of the present paper, ionization is not treated self-consistently and its profile (cosine in the axial direction, as in III.B, and uniform in the azimuthal direction) is imposed.

The domain is Cartesian and is periodic in the azimuthal direction, as shown in Figure 9. As in section III.B, only elastic electron-neutral collisions are taken into account, with a constant frequency of 5 × 10 6 s -1 i.e. roughly corresponding to a pressure around 1 mtorr. The cathodeanode distance, d, is 1.5 cm and the position of the ionization region is as in section III.B. The distance between the end of the ionization region at x 1 , and the anode (i.e. the length of the near-anode region) is about 1 cm. The width of the domain in the (periodic) azimuthal direction is 6 cm in the simulations below. The magnitude of the ionization source term is fixed as in the 1D case by imposing an injection current density 𝐽 𝑆 = 𝑒 ∫ 𝑆(𝑥)𝑑𝑥 𝑑 0

. Electrons and ions are injected in the ionization region with Maxwellian velocity distribution functions at temperatures 𝑇 e = 5 eV and 𝑇 i = 0.1 eV respectively. The electron temperature is maintained at 𝑇 e in the ionization region by changing the particle velocities according to a Maxwellian distribution at a constant rate of 5 × 10 7 s -1 . In the results presented here the ion masses are those of helium or xenon (convergence is faster for lighter ions since the growth rate of the instabilities scales as 𝜔 𝐿𝐻 ).

A. Formation of the instability

Figure 10 shows the axial-azimuthal distribution of the ion density at different times during the growth of the instability in helium for a magnetic field B=30 G. The initial condition is the 1D steady state solution (which corresponds to a positive anode sheath in helium at 30 G, see section 3). The 1D solution is clearly not stable and instabilities in the azimuthal direction quickly develop in the form of plasma filaments in a few hundreds of ns and propagate in the negative y (P B) direction. The wavelength of the instability, 𝜆, is initially small, approximately 2 mm at t=600 ns, i.e. on the order of the electron gyroradius at 30 G. The wavelength then increases in time (see Figure 10 at t=2 s) while the filaments tend to merge and form larger structures. At t=6 s macroscopic structures of size in the cm range can be clearly seen but small size fluctuations in the mm range are still present.

The azimuthal electric field, shown in Figure 11, exhibits fluctuations with amplitudes on the order of 5 × 10 3 V/m and the time evolution of the structures observed on the azimuthal field is consistent with that of the ion density. Note the increase of the axial component of the wavevector (neglected in the dispersion relations of section IV) in the non-linear stage of the instability.

Figure 12 shows the growth as a function of time, of the ion density relative azimuthal fluctuations at a given axial position 𝑥 0 = 0.5𝑑 , and for different values of the magnetic field, defined by the ratio:

𝛿𝑛 𝑟𝑚𝑠 𝑛 0 (𝑥 0 , t) = [∫ 𝑛 𝑖 2 (𝑥 0 , 𝑦, 𝑡) 𝑑𝑦 𝑤 ⁄ 𝑤 0 -(∫ 𝑛 𝑖 (𝑥 0 , 𝑦, 𝑡) 𝑑𝑦 𝑤 ⁄ 𝑤 0 ) 2 ]
1/2 ∫ 𝑛 𝑖 (𝑥 0 , 𝑦, 𝑡) 𝑑𝑦 𝑤 ⁄ 𝑤 0 where 𝑤 is the length of the domain in the azimuthal direction and 𝑛 0 = ∫ 𝑛 𝑖 (𝑥 0 , 𝑦, 𝑡) 𝑑𝑦 𝑤 ⁄ 𝑤 0

is the ion density at the axial position 𝑥 0 and time t, averaged over the azimuthal direction.

The results of Figure 12 show that the characteristic time of the instability growth in the linear stage deduced from the simulations is in the 100 ns range. This is consistent with a growth rate on the order of the lower hybrid angular frequency 𝜔 𝐿𝐻 (for B=30 G in helium, 𝜔 𝐿𝐻 ≈ 6 × 10 6 𝑠 -1 ), as predicted by the theory (see Figure 7a). The dashed lines in Figure 12b correspond to a growth rate equal to 𝜔 𝐿𝐻 . To compare more precisely the theoretical growth rate with the predictions of the simulations, we can estimate V * and V 0 in the simulations. For the initial conditions of the simulations the plasma density gradient length is on the order of 1 cm so V * = T e (𝐿 𝑛 𝐵𝑐 𝑠 ) ⁄ ≈ 16 . The parameter V 0 = 𝐸 0 (𝐵𝑐 𝑠 ) ⁄ is more difficult to estimate in these conditions but can be considered smaller than 3 (the potential increase around 𝑥 0 at the onset of the instability growth is less than 1 V over 1 cm). For these values of V 0 and V * , the maximum growth rate γ predicted by the theory (Figure 7a) is between 𝜔 𝐿𝐻 and 1.5𝜔 𝐿𝐻 , which is consistent with the simulations of Figure 12. The growth rate increases with increasing magnetic field as expected from the 𝜔 𝐿𝐻 scaling. For magnetic fields below 20 G, the sheath becomes negative in helium in our conditions (see section III) and, as can be seen in Figure 12 the plasma density fluctuations become weaker and they are no longer significant below B=15G. The theoretical wavelength at the maximum growth rate displayed in Figure 7c, is on the order of one to two electron Larmor radii for the above values of V 0 and V * , which is also consistent with the results of the PIC simulations shown in Figure 10 and Figure 11 during the linear growth. Finally, the expected wave phase velocity in this range of values of V 0 and V * , should be very close to the ion sound velocity, which is also in agreement with the PIC simulations.

B. Saturation -Fine structures and macroscopic structures

The instability saturates after a few s in the example above (in helium) and, after the linear growth, the small filamentary structures tend to merge to form larger structures. The evolution in time of the instability after the linear phase can be understood as follows. In the filamentary structure of Figure 10 at t=1 s, the plasma density gradient tends to be smaller (increase of 𝐿 𝑛 ) along the plasma filaments and larger between filaments, which means that V * decreases in a filament. According to Figure 7 this corresponds to a shift of the maximum growth rate toward values where the wavelength at maximum growth rate is larger. This can explain the formation of larger structures in the azimuthal direction as the plasma expands in the near-anode region, and the evolution from microinstabilities to rotating spokes. The typical shape of the azimuthally non-uniform plasma structure during the non-linear phase is represented in Figure 13a. The azimuthal size of the larger scale non-uniformities is on the order of 3 cm in these conditions, i.e. half of the azimuthal period of the simulation domain. The finite azimuthal period of the simulation domain clearly defines an upper limit for the wavelength of these structures. It is interesting to see on Figure 13a that small structures, of size on the order of the electron gyroradius coexist with more macroscopic structures. It seems that the small structures grow at the edges of the larger scale structures, i.e. where the density gradients are larger. Figure 13b displays the axial profiles of the ion density at two different azimuthal positions of the contour plot of Figure 13a. We see that the axial plasma density gradient varies along the azimuthal direction and that the finer structures tend to grow in the region of larger density gradients.

Although the regime of Figure 13 is quite non-linear, we note that the formation of fine structures at large density gradients, i.e. for smaller values of 𝐿 𝑛 and larger values of 𝑉 * , is consistent with the fact that the wavelengths of the instabilities at maximum growth rate tends to decrease with increasing V * .

C. Propagation velocity of the instability

The phase velocity of the instability in the azimuthal direction, estimated from the time evolution of the ion density structures, is slightly below the ion acoustic velocity in helium for 𝑇 e = 5 eV, i.e. 1.1 × 10 4 m/s. We did not notice a significant increase of this velocity with increasing magnetic field between 20 and 40 G (such behavior was observed by Powis et al. 39 with PIC simulation in a cylindrical domain but for larger values of the magnetic field).

A simple way to measure the propagation velocity of the wave is to look at the current collected on a small portion of the anode. Figure 14 shows the total electron current per unit length at the anode and the current collected on a 1.2 length segment of the anode. The total electron current per unit length reaches a constant value on a short time scale, on the order of 1-2 s. This current is consistent with the imposed injection current per unit length 𝐽 𝑆,𝑤 = 𝑒 ∫ 𝑑𝑦 𝑤 0 ∫ 𝑆(𝑥)𝑑𝑥 𝑑 0 = 1.3 A m ⁄ (the ion current to the anode is small so the electron current to the anode is close to the injection current density, i.e. 𝐽 𝑒𝑎,𝑤 ≈ 𝐽 𝑆,𝑤 ). The electron current per unit length collected by the anode segment shown in Figure 14 presents strong and regular oscillations at a frequency around 230 kHz after time t=10 s. Since the azimuthal dimension of the macroscopic structure or rotating spoke is on the order of 3 cm (Figure 13), its velocity is about 7 × 10 3 m/s i.e. slightly lower than the ion sound speed in helium for an electron temperature of 5 eV.

D. Results in xenon

Very similar results have been obtained in xenon. As expected, the time scale of the evolution of the instability is significantly increased with respect to helium because the growth rate scales as the lower hybrid frequency 𝜔 𝐿𝐻 , i.e. as the inverse square root of the ion mass for a given magnetic field. Because of the differences in the ion sound velocity, the parameters V * and V 0 are larger by a factor of about 6 in xenon compared with their values in helium for the same values of 𝐵, T e 𝐿 𝑛 ⁄ , and 𝐸 0 . The time evolution of the rms fluctuations of the ion density in the mid-gap for different values of the magnetic field intensity are shown in Figure 15. The growth rate roughly scales with the lower hybrid angular frequency for magnetic fields larger than 40 G. For 𝐵 = 40 G , the growth rate is lower than 𝜔 𝐿𝐻 , and for 𝐵 = 30 G the instability is not significant.

The fluctuations of the ion density in the non-linear phase can be seen in Figure 16 for different values of the magnetic field intensity. Some fluctuations are present for 𝐵 = 40 G although the value of the magnetic field is very close to the magnetic field corresponding to the transition from negative to positive anode sheath in the 1D model (see section III). The amplitude of the fluctuations is however relatively small and the instability remains in a short wavelength, filamentary structure even for large simulation times. When the magnetic field is increased the filaments tend to merge and form larger structures. The number of these larger structures decreases as the magnetic field increases, and only one macroscopic structure is present for B=100 G (note however that the structures become less coherent as the magnetic field is increased).

The time variations of the total electron current on the anode and on a segment of the anode are plotted in Figure 17 for the 𝐵 = 80 G case. The fluctuations of the current on the anode segment correspond to the motion of the plasma structures in the azimuthal direction. The apparent increase of the time interval between current peaks is a signature of the formation of larger structures as time evolves. From the last two peaks, at 80 and 110 s, and since the size of the plasma structures around these times is about half the width of the simulation domain as can be seen in Figure 16, one can infer that the velocity of the instability in the azimuthal direction is in the 1000 m/s range. The ion acoustic velocity of xenon ions for 𝑇 𝑒 = 5 eV is 1.9 × 10 3 m/s.

VI. DISCUSSION

A. Assumptions and relevance to the near-anode region of magnetized plasma devices

The model presented in this paper is based on two important simplifying assumptions: 1) negligible ionization in the near-anode region, and 2), constant transverse magnetic field.

The approximation of negligible ionization (and constant electron temperature) was very useful for studying the onset azimuthal instabilities and for comparing with the predictions of the dispersion relation. With respect to applications and to the particular case of the near-anode region of a Hall thruster we believe that this approximation is not good in a global regime of oscillations of a Hall thruster, i.e. under conditions where "breathing oscillations" are important. In such regime it is expected that the location of the ionization region is modulated in time and a non-negligible part of the overall ionization may take place close to the anode. Under conditions of small or non-existent breathing oscillations, ionization takes place closer to the exhaust plane and ionization in the near-anode region becomes less important. The exhaustive measurements and simulation results presented for example by the JPL group on a 6 kW thruster show that the electron temperature in the first half length of the channel from the anode is less than 5 eV 60 in this thruster. The influence of ionization is this rotating regime remains to be clarified.

The assumption of constant magnetic field allows simple comparisons with the predictions of dispersion relations, but is clearly not realistic in the case of a Hall thruster where the radial magnetic field decreases from values on the order of 150-200 G in the channel exhaust plane to values of 10 G or less at the anode. The results presented here give a useful estimate of the average magnetic field (more precisely root mean square magnetic field) in the near-anode region below which electron transport can be considered as classical and above which azimuthal instabilities leading to anomalous transport are likely to develop. The effect of the magnetic field gradient present in the near-anode region of Hall thrusters will be explored in a future work and the results from PIC simulations will be compared with the predictions from the theory 9,10,20,48 .

B. Anomalous transport induced by the instability

The anomalous transport induced by the instability in the near-anode region can be quantified by calculating an effective electron collision frequency 𝜈 𝑒𝑓𝑓 , cross-field electron mobility, 𝜇 𝑒𝑓𝑓 or effective Hall parameter h 𝑒𝑓𝑓 . This is done by averaging the drift diffusion expression of the electron flux in time and space. The time integration is made in the non-linear steady state regime over a long enough time  ( is on the order or larger than 200 s in the results below). The integration in space is made over the near-anode region and in the azimuthal direction. Since the electron current density to the anode, 𝐽 𝑆 , is fixed and is directly related to the ionization source term 𝑆, we can write:

𝐽 𝑆 𝑒 = ∫ 𝑆(𝑥)𝑑𝑥 𝑑 0 =𝜇 𝑒𝑓𝑓 1 𝜏 1 𝑤 ∫ 𝑑𝑥𝑑𝑦 ∫(𝑛∇ϕ -T e ∇𝑛) 𝑑𝑡
This gives:

ℎ 𝑒𝑓𝑓 = 1 𝜇 𝑒𝑓𝑓 𝐵 = 𝜔 𝑐𝑒 𝜈 𝑒𝑓𝑓 ≈ 𝑒 𝐽 𝑆 1 𝐵 ( T e 𝐿 𝑛 𝑛 0 + 𝑛𝐸)
where 𝑛 0 is the time averaged plasma density at the end of the ionization region (𝑥 = 𝑥 1 , see Figure 2), and 𝑛𝐸 is an average in time and space over the near-anode region (between 𝑥 = 0 and 𝑥 = 𝑥 1 and 𝑦 = 0 and 𝑦 = 𝑤).

The effective Hall parameter ℎ 𝑒𝑓𝑓 can also be obtained from the ratio of the time and space average azimuthal and axial electron current in the near anode region. The effective Hall parameters calculated with these two methods are in good agreement, with differences smaller than 20%.

Figure 18 shows the variations of the effective Hall parameter and electron mobility deduced from the simulations, as a function of magnetic field, for helium and xenon. In the considered range of magnetic fields above the negativepositive anode sheath transition, the effective Hall parameter tends to decrease with increasing magnetic field, while the effective mobility is practically constant.

The effective Hall parameter in helium is below the Bohm value (𝜇 𝐵𝑜ℎ𝑚 = 1 16𝐵 ⁄ , i.e. ℎ 𝐵𝑜ℎ𝑚 = 16), is on the order of 15 at 20 G and decreases to 5 at 80 G. In xenon, ℎ 𝑒𝑓𝑓 is larger, on the order of 40 at 40 G, and decreases to 14 at 100 G. Note that these values of the effective Hall parameters are significantly smaller than the values (more than 100) obtained with PIC simulations in the region downstream of the maximum magnetic field, i.e. in the acceleration region of a Hall thruster, where the Electron Cyclotron Drift Instability evolving toward a modified ion acoustic instability was found to be dominant 43 . This means that anomalous transport can be large in the near-anode region, and can be responsible for rotating structures if the local magnetic field is large enough.

C. Short wavelength and long wavelength modes

The short wavelength mode is present in the simulation during the linear growth of the instability (first few s in Figure 10 and Figure 11 in helium, for B=30 G), and can be present even in a steady state regime for low values of the magnetic field (Figure 16 in xenon for B=40 G). B=30 G is above the transition field defined in section III in helium and B=40 G is around the transition field in xenon.

For magnetic fields close to the transition field, i.e. around B=40 G in xenon, the simulations shows that the time averaged electric field in the near-anode region is close to zero and the time averaged plasma potential is above the anode potential, as can be seen in Figure 19. The short wavelength instability that is seen in the B=40 G case of Figure 16 is therefore not a Simon-Hoh instability but corresponds to an ion acoustic wave destabilized by density gradient and collisions (i.e. the 𝑉 0 = 0 part of Figure 7). The wavelength in the B=40 G case of Figure 16 is on the order of 2 mm, close to the electron Larmor radius 𝜌 𝑒 = (𝑒𝑇 𝑒 𝑚 ⁄ ) 1/2 𝜔 𝑐𝑒 ≈ 1.33 mm. ⁄ The dispersion relation for these conditions corresponds to the 𝐸 0 = 0 case of Figure 5, which predicts a maximum growth rate for 𝑘𝜌 𝑒 ≈ 6.2, i.e. for a wavelength 𝜆 = 2𝜋 𝑘 ⁄ ≈ 1.33 mm. It is possible to show that when charge separation effects are taken into account in the dispersion relation (Debye length term), the theoretical wavelength at maximum growth rate is larger in our simulation case with relatively low plasma density, and is closer to the 2 mm wavelength observed in the B=40 G case of Figure 16.

Simulations with different ion masses show that the acoustic wave observed around the transition point is very coherent and exhibit ion trapping effects. Figure 20 shows a (y-v y ) typical phase space plot of the ions at a given time of the simulation in the steady state regime for helium at 20 G (shorter simulation domain in the azimuthal direction, w=1.5 cm), and for xenon, 40 G, w= 6 cm (Figure 16). For the helium, B=20 G case, the electron Larmor radius is 𝜌 𝑒 ≈ 2.66 mm and we see in Figure 20 that the wavelength is close the electron Larmor radius.

For a larger magnetic field, B=100 G, in xenon, the time averaged plasma potential and plasma density are quite different from the B=40 G case (Figure 19). 16). The length of the simulation domain in the azimuthal direction in the helium case was shorter in this example (w= 1.5 cm).

The plasma potential at 100 G slightly increases in the direction of the anode in part of the anode region. The average electric field is therefore in the same direction as the plasma density gradient and the instability observed the B=100 G case of Figure 16 is a Simon-Hoh instability. This is also true for the 60 and 80 G cases of Figure 16 in xenon, and for the 30 G case of Figure 10 in helium.

VII. SUMMARY AND CONCLUSION

The onset of instabilities in a low temperature partially magnetized plasma under conditions typical of the nearanode region of Hall thrusters or magnetron discharges has been studied with a simplified model. The model assumes a prescribed ionization source on the cathode side of a magnetized discharge with a given applied voltage, and describes electron and ion transport from the ionization region to the anode through the "near-anode region". In this region, the magnetic field is relatively low, but the Hall parameter due to electron-neutral collisions is still large, on the order of 100 or larger. Only elastic collisions are considered for electrons, with a constant collision frequency, while ions are supposed to be collisionless. In this first approach, the magnetic field is supposed to be constant and ionization in the near-anode region is neglected.

The transition from a negative to a positive anode sheath is observed with a 1D PIC-MCC simulation of the near-anode region, when the magnetic field is increased. A simple expression has been derived for the magnetic field at the transition between negative and positive anode sheath.

The magnetic field at the transition point can be expressed as a function of electron collision frequency, length of the near-anode region, electron temperature and ion mass and this expression is in reasonable agreement with predictions from the 1D PIC-MCC simulation.

A 2D PIC-MCC simulation has been used to study the development of instabilities in the axial-azimuthal plane perpendicular to the magnetic field. The results show that azimuthal instabilities develop for magnetic fields above the transition point, i.e. for values of the magnetic field where the 1D model predicts the formation of a positive anode sheath. Using 1D axial solutions as initial conditions of the 2D model (azimuthally uniform conditions), short wavelength (close to the electron gyroradius) instabilities start to form during the linear phase and rotate in the azimuthal direction. If the magnetic field is close the transition point, the wavelength of the instability remains short, on the order of the electron gyroradius, over times much longer than all the relevant characteristic times (e.g. electron and ion transit times in the near-anode region). The plasma appears as composed of filaments extending from the ionization region to the anode. The plasma density gradient tends to decrease along a filament. When the magnetic field is increased, some of the filaments merge, forming larger structures rotating in the azimuthal direction. The short wavelength structures can coexist with the larger structures and seem to grow from the sharp end (i.e. large density gradient regions) of the larger structures. For large enough magnetic field a single rotating structure is observed along the length of the simulation domain in the azimuthal direction

The development of the instability is consistent with fluid dispersion relations that take into account the effect of E B electron drift, diamagnetic drift, electron inertia and gyroviscosity terms, and collisions 9 . For values of the magnetic field close to the transition field, the short wavelength regime present at steady state corresponds to an ion sound mode destabilized by density gradient and collisions. For larger magnetic fields the instability evolves toward a longer wavelength mode typical of a classical collisionless Simon-Hoh regime.

The conditions of the simulations presented here (fixed electron temperature, no ionization) are oversimplified and the results cannot be directly compared to specific experiments. However, recent measurements of electron fluctuations via an adapted Thomson scattering diagnostic have shown the coexistence of millimeter scale fluctuations and rotating spokes in the plasma of a planar magnetron 28 . The authors attributed the small-scale fluctuations to the Electron Cyclotron Drift Instability (ECDI). Since the measurements seem to have been performed in the low electric field of the bulk plasma, on the anode side of the maximum plasma density (to be confirmed), the results presented here suggest that the observed fluctuations could be due to an ion acoustic wave destabilized by density gradient and collisions.

The anomalous transport induced by the instability and deduced from the simulations in the near-anode region is large even at moderate magnetic fields and the effective Hall parameter is significantly smaller than the value predicted by PIC-MCC simulations in the acceleration region of Hall thrusters and attributed to an Electron Cyclotron Drift Instability. For example, for a near-anode magnetic field of 50 G, the effective Hall parameter deduced from the simulations presented here is on the order of 40 in xenon and 10 in helium, corresponding to effective collision frequencies of 2 × 10 7 s -1 and 8 × 10 7 s -1 respectively. The Hall parameter decreases with magnetic field in the considered range, leading to a quasi-constant effective electron mobility (the effective electron collision frequency increases as 𝐵 2 ). Although the assumptions of constant magnetic field and of no-ionization in the near-anode region are restrictive in terms of applications to real situations, this work give some useful estimates of the conditions for instabilities in the near-anode region and some information on the nature and properties of these instabilities. It shows that instabilities can be present without ionization and with azimuthally uniform neutral atom densities. In the relatively low magnetic field conditions considered here, the frequency and growth rate of the instabilities were on the order of the lower hybrid frequency and the phase velocity was close to the ion sound speed.

Figure 1 :

 1 Figure 1: Typical axial profiles of the applied magnetic field, ionization rate, and electric field (a) in a Hall thruster, and (b), in a planar magnetron. The grey region characterized by a low electric field, negligible ionization and moderate magnetic field is called the "near-anode region" in this paper. The anode sheaths are negative in these two examples since the sheath electric field repels electrons and attracts ions.

Figure 1

 1 Figure 1 also displays the ionization region and the region of large electric field, downstream of the ionization region.In a Hall thruster (Figure1a) this large field region, called the acceleration region, is located in the quasineutral plasma. The ion beam is extracted on the right of the acceleration region, and an emissive cathode allows neutralization of the ion beam as well as providing a source to sustain the discharge. In a magnetron discharge (Figure1b) the large field region is the ion sheath next to the cathode. In both discharges, electrons accelerated in the large field region release their energy and ionize the gas in the ionization region. The thickness of the ionization region is relatively small because of the large magnetic field and is on the order of 1 cm in the conditions of Hall thrusters[START_REF] Goebel | Fundamentals of Electric Propulsion: Ion and Hall Thrusters[END_REF][START_REF] Boeuf | [END_REF] .In "modern" Hall thrusters the radial magnetic field is maximum around the channel exhaust plane (Figure1a) and the near-anode region is relatively short (less than 2 cm). Proper operation of a Hall thruster is achieved when most of the ionization takes place immediately upstream of the acceleration region and we will assume here, as said above, that ionization is negligible in the near-anode region. This supposes a significant decrease of the electron temperature

Figure 2 :

 2 Figure 2 : (a) Profiles of the magnetic field and ionization rate used in the model of the negative to positive sheath transition; ionization is neglected in the near-anode region between the anode at 𝑥 = 0 and the end of the ionization region at 𝑥 = 𝑥 1 ; (b) typical expected profiles of the electric field (case of a negative sheath) and electron and ion densities. The dimensions in the simulations are: 𝑑 = 1.5 cm, 𝑥 1 = 0.7𝑑 = 1.05 cm, 𝑥 2 = 0.9𝑑 = 1.35 cm

Figure

  Figure3: Magnetic field at the negative-positive sheath transition point, as a function of ion mass for an electron temperature 𝑇 𝑒 =5eV, an electron-neutral collision frequency of 5 × 10 6 s -1 , and for two values of the density gradient length, with 𝛼 = 0.5. The colored full lines are from the analytical expression (see text) and the circle symbols are obtained from the PIC-MCC simulations. For each value of 𝐿 𝑛 the sheath is negative below the curves and become positive above. The same electron temperature is used in the nearanode region and in the Bohm velocity 𝑐 𝑆 on the cathode side.

Figure 4 :

 4 Figure 4 : 1D PIC-MCC simulations in helium showing the transition from negative to positive anode sheath: (a) axial profiles of the electron and ion densities for two values of the magnetic field: 18 G (negative anode sheath) and 21 G (positive anode sheath); (b) profiles of the potential in the near-anode region for three values for the magnetic field. Conditions: applied voltage 𝑈 𝐷 = 200 𝑉, gap 𝑑 = 1.5 cm, given (cosine) profile of the ionization source term 𝑆(𝑥) between 𝑥 1 = 1.05 cm and 𝑥 2 = 1.45 cm, with 𝐽 𝑠 = ∫ 𝑒𝑆(𝑥)𝑑𝑥 𝑑 0 = 30 A m 2⁄ , electron and ion temperatures 𝑇 𝑒 = 5 eV, 𝑇 𝑖 = 0.1 eV (full lines) or 𝑇 𝑖 = 0.5 eV (dashed lines), constant electron-neutral momentum exchange collision frequency 𝜈 𝑒𝑛 = 5 × 10 6 𝑠 -1 (no inelastic collisions), collisionless ions.

Figure 5 :

 5 Figure 5 : (a) Growth rate and (b) angular frequency, normalized to the lower hybrid angular frequency as a function of azimuthal wavenumber times electron Larmor radius for three values of the axial electric field E0 in the plasma for xenon, B=40 G, (i.e. 𝜔 𝐿𝐻 =1.43 × 10 6 rad/s), 𝑇 𝑒 = 5 eV , Ln= 1 cm ( T e /Ln =5 V/cm).The electric field and density gradient are in the same direction. The electron-neutral collision frequency is set to 𝜈 𝑒𝑛 = 5 × 10 6 s -1 (full lines) or to 10 6 𝑠 -1 (dotted lines, (a) only). The dashed line in (b) corresponds to the ion acoustic instability dispersion relation 𝜔 = 𝑘𝑐 𝑠 (i.e. 𝜔 𝜔 𝐿𝐻 ⁄ = 𝑘𝜌 𝑒 ). In normalized units, the three cases correspond to 𝑉 0 = 0, 6.5, and 26, with 𝑉 * =65.

Figure 6 :

 6 Figure 6 : : (a) Growth rate and (b) angular frequency, normalized to the lower hybrid angular frequency as a function of azimuthal wavenumber times electron Larmor radius for three values of the density gradient length Ln , for xenon, B=40 G, (i.e. 𝜔 𝐿𝐻 =1.43 × 10 6 rad/s), 𝑇 𝑒 = 5 eV. The plasma electric field E0 is set to 2V/cm. The electric field and density gradient are in the same direction. The electron-neutral collision frequency is set to 𝜈 𝑒𝑛 = 5 × 10 6 s -1 The dashed line in (b) corresponds to the ion acoustic instability dispersion relation 𝜔 = 𝑘𝑐 𝑠 (i.e. 𝜔 𝜔 𝐿𝐻 ⁄ = 𝑘𝜌 𝑒 ).In normalized units, the three cases correspond to 𝑉 * = 0, 13, and 32, with 𝑉 0 =13.

Figure 7 : 3 . 5 .

 735 Figure 7 : Contours of constant, normalized (a) maximum growth rate, (b) angular frequency at the maximum growth rate, and (c) wavelength at the maximum growth rate, as a function of 𝑉 0 = 𝐸 0 (𝐵𝑐 𝑠 ) ⁄ and 𝑉 * = T e (𝐿 𝑛 𝐵𝑐 𝑠 ) ⁄ for 𝜈 𝑒 𝜔 𝐿𝐻 ⁄ = 3.5. The axis on the top and on the right of the plots give the corresponding values of 𝐸 0 and 𝑇 e 𝐿 𝑛 ⁄ for the particular case of xenon with 𝐵 = 40 G, 𝑇 𝑒 = 5 eV, and 𝜈 𝑒𝑛 =5× 10 6 s -1 . The value of 𝜔 𝐿𝐻 in that case is 1.43 × 10 6 rad/s. The growth rate and angular frequency are normalized to the lower hybrid angular frequency and the wavelength of the instability is normalized to the electron Larmor radius.

  respectively, in the form of contour plots as a function of 𝑉 0 = 𝐸 0 𝐵𝑐 𝑠 and V * = T e 𝐿 𝑛 𝐵𝑐 𝑠 . A scale showing the corresponding values of 𝐸 0 and T e 𝐿 𝑛

Figure 8 :

 8 Figure 8 : Phase velocity 𝑉 w of the azimuthal wave normalized to the ion sound velocity 𝑐 𝑠 , as a function of 𝑉 0 = 𝐸 0 (𝐵𝑐 𝑠 ) ⁄ and 𝑉 * = T e (𝐿 𝑛 𝐵𝑐 𝑠 ) ⁄ in the same conditions as Figure 7.

Figure 9 :

 9 Figure 9 : axial-azimuthal domain of the 2D PIC-MCC simulations. In the results presented in this paper, the dimensions are taken as: 𝑑 = 1.5 cm, 𝑥 1 = 0.7𝑑 = 1.05 cm, 𝑥 2 = 0.9𝑑 = 1.35 cm.The results presented here have been obtained for a periodic length 𝑤 = 6 cm. The number of grid points is 128 × 512 and most of the simulations are performed with about 100 particles per cell.

Figure 10 :

 10 Figure 10: Axial-azimuthal (x,y) distributions of the ion density at four different times; helium mass, 𝐵 = 30 G, 𝜈 𝑒𝑛 = 5 × 10 6 . Scale: 1.8, 1.3, 1.4, 1.5x10 16 m -3 for t=0.6, 1, 2, and 3 s respectively. The current per unit length 𝐽 𝑆,𝑤 (see text) is equal to 1.3 A/m. The wave propagates in the negative y direction.

Figure 11 :

 11 Figure 11 : Axial-azimuthal (x,y) distributions of the azimuthal component of the electric field in the conditions of Figure 10. The units, in kV/m, are [-6.3,+4.7], [-6.7,+5.2], [-5.9,+5.4], [-5.2,+4.3] for 0.6, 1, 2, and 3 s, respectively.

Figure 12 :

 12 Figure 12 : Time variations of the relative azimuthal rms fluctuations of the ion density at the axial position 𝑥 0 = 0.5𝑑 on two different time scales for different values of the magnetic field in helium (constant electron-neutral momentum exchange frequency of 5 × 10 6 s -1 in the conditions of Figure 10 and Figure 11. The dashed lines in (b) correspond to a growth rate equal to 𝜔 𝐿𝐻 for each value of the magnetic field.

Figure 13 :Figure 14 :

 1314 Figure 13 : Axial-azimuthal distribution of the ion density (top) and profiles of the ion density along the axial position at two azimuthal positions, y=0.38 cm and y=2.13 cm at time t=35 s, same conditions as Figure 10, helium mass, B=30 G, electron momentum exchange frequency 5 × 10 6 s -1 . MIN MAX

Figure 15 :

 15 Figure 15 : Time variations of the relative azimuthal rms fluctuations of the ion density at the axial position 𝑥 0 = 0.5𝑑 on two different time scales for different values of the magnetic field for xenon (constant electron-neutral momentum exchange frequency of 5 × 10 6 s -1 ) in in the same geometry as above. The dashed lines in (b) are shown for comparison and correspond to a growth rate equal to 𝜔 𝐿𝐻 for each value of the magnetic field.

Figure 16 :

 16 Figure 16 : Axial-azimuthal (x,y) distributions of the ion density for four values of the magnetic field, 40, 60, 80, and 100 G; xenon mass, 𝜈 𝑒𝑛 = 5 × 10 6 s -1 . Scale: 1.8, 1.3, 1.4, 1.5 × 10 16 m -3 . The plots correspond to simulation times larger than 100 s for all magnetic fields. The imposed current per unit length 𝐽 𝑆,𝑤 is equal to 0.16 A/m.

Figure 17 :

 17 Figure 17 : Total electron current per unit length on the 6 cm anode and current per unit length on a 1.2 cm segment of the anode, showing the motion of the instability in the azimuthal direction. Same conditions as Figure 16, xenon mass, B=80 G, electron momentum exchange frequency 5 × 10 6 s -1 .

Figure 18 :

 18 Figure18: Effective Hall parameter (symbols) and electron mobility (dashed line) in the near-anode region as a function of magnetic field, for xenon (red) and helium (blue).

Figure 19 :

 19 Figure 19 : Axial profiles of the time averaged (a), ion densities, and (b), plasma potential for B=40 G and B=100 G in xenon in the conditions of Figure 16 .

Figure 20 :

 20 Figure 20 : Azimuthal, (y,vy) , phase plots of the ions in helium for B=20 G (top) and xenon for B=40 G (bottom, conditions of Figure16). The length of the simulation domain in the azimuthal direction in the helium case was shorter in this example (w= 1.5 cm).

16 Future work will include the effects of magnetic field gradients and ionization.


Acikgoz Degenerate

Poisson Distribution

Ugur Duran
email: ugur.duran@iste.edu.tr

Mehmet Acikgoz
email: acikgoz@gantep.edu.tr

Degenerate Poisson Distribution

Keywords: 2010 Mathematics Subject Classi…cation. Primary 60E05, Secondary 11B73, 33B10 Poisson distribution, raw moments, Bell polynomials, degenerate exponential function, unsigned Stirling number of the …rst kind

come

Introduction

Probability distributions possess many importance and have been commonly used in statistics and probability theory. They can be given by a table or an equation which links each outcome of a statistical experiment with its probability of occurrence, see [START_REF] Evans | Statistical distributions, 2nd Edition[END_REF][START_REF] Frank | Handbook of the Poisson Distribution[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF] for more information about the topic on probability distribution. One of the most useful and applicable distribution is the Poisson probability distribution (cf. [START_REF] Evans | Statistical distributions, 2nd Edition[END_REF][START_REF] Frank | Handbook of the Poisson Distribution[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF]). The Poisson distribution is a discrete probability distribution that expresses the probability of a given number of events occurring in a …xed interval of time or space if these events occur with a known constant rate and independently of the time since the last event. The Poisson distribution can be also utilized for the number of events in other speci…ed intervals such as volume, area or distance, cf. [START_REF] Evans | Statistical distributions, 2nd Edition[END_REF][START_REF] Frank | Handbook of the Poisson Distribution[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF] and see also the references cited therein.

Several generalizations of the Poisson distribution such as the Polya process, truncation of the domain, stuttering Poisson distribution, linear combinations of Poisson variables, particle counting, collective risk theory, mixed Poisson distributions, Erlang process and Poisson mixing distributions have been deeply discussed in the book [START_REF] Frank | Handbook of the Poisson Distribution[END_REF] in 1967. In recent years, some scientists have been considered several new extensions of the Poisson distribution (cf. [START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF]). For instance, Hermann [START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF] considered a generalization of the Poisson distribution based on the generalized Mittag-Le-er function and calculated the corresponding raw moments algebraically in terms of the familiar Bell polynomials. Mahmoudi [START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF] introduced the exponentiated Weibull-Poisson (EWP) distribution which includes four parameters as increasing, decreasing, bathtubshaped and unimodal failure rate and satis…ed diverse properties for this distribution covering moments, quantiles, its probability density function and its reliability and failure rate functions. Porwal [START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF] de…ned Mittag-Le-er type Poisson distribution and then provided several properties of this distribution. In this paper, we consider a novel extension for the Poisson distribution which includes the degenerate exponential function and we call it the degenerate Poisson distribution. Then, we investigate multifarious properties of the degenerate Poisson distribution covering the …rst and the second raw moments. We also attained the skewness and the kurtosis for this distribution. Furthermore, we derive its moment generating function by which we consider the degenerate Bell polynomials and present a connection for the mentioned polynomials associated with the unsigned Stirling numbers of the …rst kind.

Here we provide some basic information taken from the references [START_REF] Evans | Statistical distributions, 2nd Edition[END_REF][START_REF] Frank | Handbook of the Poisson Distribution[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF].

De…nition 1. The n-th moment of a discrete probability distribution about X = 0 is de…ned by n = E (X n ) and the moment generating function (m.g.f.) of a random variable X is denoted by M X (t) and de…ned by

M X (t) = E e tX :
The probability density function of Poisson distribution with parameter is de…ned as follows

p x ( ) = P (X = x) = 1 N x
x! for x = 0; 1; 2; : : : ;

(1.1)
where N is the normalization constant and is determined by the requirement of normalizability of the distribution, which coincides with the zeroth raw moment 0 :

0 = 1 X x=0 p x ( ) = 1 X x=0 1 N x x! = 1 N 1 X x=0 x x! = 1 N e = 1 (1.2)
which implies N = e . The …rst and the second raw moments for the Poisson distribution are given below:

0 = E (X) = and 1 = E X 2 = + 2 ; (1.3)
therefore, we obtain the mean E (X) and the variance V (X) of the Poisson distribution (Mean) E (X) = and (Variance) V (X) = :

(1.4)

The moment generating function of the Poisson distribution is as follows:

E e xt = 1 X x=0 e xt p x ( ) = e 1 X x=0 (e t ) x x! = e (e t 1) = 1 X n=0 B n ( ) t n n! ; (1.5)
which is the generating function of the classical Bell polynomials B n ( ) (cf . [START_REF] Bell | Exponential polynomials[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Kim | Some identities of Bell polynomials[END_REF][START_REF] Kim | On degenerate Bell numbers and polynomials[END_REF][START_REF] Kim | On partially degenerate Bell numbers and polynomials[END_REF]). So, for n 0, we write the higher raw moments by means of the Bell polynomials as follows:

n = E (X n ) = B n ( ) .
(1.6)

A more detailed information about the topics mentioned above can be found in [START_REF] Evans | Statistical distributions, 2nd Edition[END_REF][START_REF] Frank | Handbook of the Poisson Distribution[END_REF][START_REF] Hermann | Generalization of the fractional poisson distribution[END_REF][START_REF] Mahmoudi | Exponentiated Weibull-Poisson distribution: Model, properties and applications[END_REF][START_REF] Porwal | On Mittag-Le-er type Poisson distribution[END_REF] and each of the related references cited therein.

Degenerate Poisson Distribution

In this section, we consider an extension for the Poisson distribution via the degenerate exponential function e x (t) and the -falling factorial (x) n; . We call it the degenerate Poisson distribution. We then perform to obtain the …rst few raw moments and the moment generating function for the aforesaid distribution. By means of this moment generating function, we de…ne the degenerate Bell polynomials and …nd the properties for these polynomials.

For being a real number, the -falling factorial (x) n; is de…ned by (see [START_REF] Dattoli | Finite sums and generalized forms of Bernoulli polynomials[END_REF][START_REF] Duran | On Gould-Hopper-based fully degenerate poly-Bernoulli polynomials with a q-parameter[END_REF])

(x) n; = x(x )(x 2 ) (x (n 1
) ); n = 1; 2; : : :

1; n = 0: (2.1)
The degenerate exponential function e x (t) for a real number is given by (cf. [START_REF] Duran | On Gould-Hopper-based fully degenerate poly-Bernoulli polynomials with a q-parameter[END_REF][START_REF] Kim | On degenerate Bell numbers and polynomials[END_REF][START_REF] Kim | On partially degenerate Bell numbers and polynomials[END_REF]) It is readily seen that lim !0 e x (t) = e xt (cf. [START_REF] Duran | On Gould-Hopper-based fully degenerate poly-Bernoulli polynomials with a q-parameter[END_REF][START_REF] Kim | On degenerate Bell numbers and polynomials[END_REF][START_REF] Kim | On partially degenerate Bell numbers and polynomials[END_REF]). From (2.1) and (2.2), the degenerate exponential function possesses the following series expansion

e x (t) = (1 + t)
e x (t) = 1 X n=0 (x) n; t n n! . (2.3)
We now introduce a generalization of the Poisson distribution via the -falling factorial (x) n;r as follows:

p x ( ; ) = P (X = x) = 1 N ( ) x;
x! for x = 0; 1; 2; : : : ;

(2.4) the normalization constant N follows from the requirement of normalizability of the distribution:

0 = 1 X x=0 p x ( ; ) = 1 N 1 X x=0 ( ) x; x! = 1 N e (1) = 1; (2.5)
and so N = e (1) = (1 + ) . Thus, we state our main de…nition as follows.

De…nition 2. The degenerate Poisson distribution with two parameters and is de…ned by the following formula p x ( ; ) = P (X = x) = e (1) ( ) x;

x! for x = 0; 1; 2; : : : :

We now investigate some properties of the degenerate Poisson distribution. We …rstly give the mean and variance of the mentioned distribution as follows.

Theorem 1. The mean E (X) and the variance V (X) of the degenerate Poisson distribution are given by

E (X) = 1 + and V (X) = (1 + ) 2 .
(2.7)

Proof. By De…nition 2, we acquire

1 = E (X) = 1 X x=0 xp x ( ; ) = 1 X x=0 xe (1) ( ) x;
x! = e (1)

1 X x=1 ( ) x;
(x 1)! = e (1)

1 X x=0 ( ) x+1;
x! = e (1) which is the asserted mean in (2.7) and

1 X x=0 ( )( 2 ) ( (x 1) )( x ) x! = e (1)
2 = E X 2 = 1 X x=0 x 2 p x ( ; ) = e (1)
1 X x=0 x 2 ( ) x; x! = e (1) 1 X x=1 (x 1 + 1) ( ) x;
(x 1)! = e (1)

1 X x=2 ( ) x; (x 2)! + e (1) 1 X x=1 ( ) x;
(x 1)! = e (1)

1 X x=0 ( )( 2 ) ( (x 1) )( x ) x! + e (1) 1 X x=0 ( ) x;
x! = e (1) ( )

1 X x=0 ( 2 ) x;
x! + e (1)

1 X x=0 ( ) x;
x! = e (1) ( )e 2 (1) + e (1) e

= ( )

(1 + ) 2 + 1 + = 2 + (1 + ) 2 :
Therefore, we get the variance

V (X) = E X 2 (E (X)) 2 = 2 + (1 + ) 2 2 (1 + ) 2 = (1 + ) 2 ;
which is the desired variance in (2.7).

Remark 1. As ! 0, the mean and variance in (2.7) reduce to the mean and variance for familiar Poisson distribution in (1.3).

The third raw moment (skewness) and fourth raw moment (kurtosis) are given by the following theorem.

Theorem 2. The skewness and the kurtosis for the degenerate Poisson distribution are given below

3 = E X 3 = ( ) ( 2 )
(1 + )

3 + 3 ( ) (1 + ) 2 + 1 + (2.8) and 4 = E X 4 = ( ) ( 2 )( 3 )
(1 + )

4 + 6 ( ) ( 2 )
(1 + )

3 + 7 ( ) (1 + ) 2 + 1 + . (2.9)
Proof. Utilizing the similar methods used in the proof of the Theorem 1, the proof of Theorem 2 is based on the following expansions x 3 = x (x 1) (x 2) + 3x (x 1) + x and x 4 = x (x 1) (x 2) (x 3) + 6x (x 1) (x 2) + 7x (x 1) + x:

So, we omit the proof. The moment generating function E e tX of the degenerate Poisson distribution is obtained by the following procedure:

Preprints

E e tX = e (1) 1 X x=0 e tx ( ) x; (x 1)! = 1 X n=0 e (1) 1 X x=0 x n ( ) x; (x 1)! ! t n n! . (2.10)
In usual, the moment generating function of the familiar Poisson distribution coincides the generating function of the usual Bell polynomials, see the formula (1.5). By this motivation, we de…ne the degenerate Bell polynomials, di¤erent from Kim's degenerate Bell polynomials in [START_REF] Kim | On degenerate Bell numbers and polynomials[END_REF][START_REF] Kim | On partially degenerate Bell numbers and polynomials[END_REF], as follows:

1 X n=0 B n; ( ) t n n! = e (1) 1 X x=0 e tx ( ) x; x! . (2.11)
In view of (2.10) and (2.11), for n 0, the degenerate Bell polynomials are given by the following series:

B n; ( ) = e (1) 1 X x=0 x n ( ) x; x! . (2.12)
Remark 2. As tends to zero, the formula given in (2.12) reduce to the well known formula for usual Bell polynomials (cf. [START_REF] Kim | Some identities of Bell polynomials[END_REF][START_REF] Kim | On degenerate Bell numbers and polynomials[END_REF][START_REF] Kim | On partially degenerate Bell numbers and polynomials[END_REF]):

B n ( ) = e 1 X x=0 x n x x! .
With the help of Theorems 1 and 2, the …rst few degenerate Bell polynomials are

B 0; ( ) = 1 B 1; ( ) = 1+ B 2; ( ) = ( ) (1+ ) 2 + 1+ B 3; ( ) = ( )( 2 ) (1+ ) 3 + 3 ( ) (1+ ) 2 + 1+ B 4; ( ) = ( )( 2 )( 3 ) (1+ ) 4 + 6 ( )( 2 ) (1+ ) 3 + 7 ( ) (1+ ) 2 + 1+ B 5; ( ) = ( ) ( 2 )( 3 )( 4 )
(1+ ) 5 + 10 ( ) ( 2 )( 3 )
(1+ ) 4 + 25 ( )( 2 ) (1+ ) 3 + 15 ( ) (1+ ) 2 + 1+ B 6; ( ) = ( ) ( 2 )( 3 )( 4 )( 5 )
(1+ ) 6 + 15 ( ) ( 2 )( 3 )( 4 )
(1+ ) 5 +65 ( ) ( 2 )( 3 )
(1+ ) 4 + 90 ( )( 2 ) (1+ ) 3 + 35 ( ) (1+ ) 2 + 1+ :
The di¤erence operator is de…ned by (see [START_REF] Dattoli | Finite sums and generalized forms of Bernoulli polynomials[END_REF])

f (x) = f (x + ) f (x) ; 6 = 0: (2.13)
Note that as goes to zero, the di¤erence operator reduce to the usual derivative operator:

lim !0 f (x) = lim !0 f (x + ) f (x) = f 0 (x) = df (x) dx :
The following di¤erence rules hold true (cf. [START_REF] Duran | On Gould-Hopper-based fully degenerate poly-Bernoulli polynomials with a q-parameter[END_REF]): x! = ( )( 2 ) ( (x 1) ) 1:2:::

k (x) n; = n! (n k)! (x) n k; ; 0 k n (2.
(x 1) x = x ( ) ( 2 )
( (x 1) ) 1:2:::

(x 1) = x ( ) 2 x 1 = x x 1 Y k=1 k . (2.19)
The unsigned Stirling numbers of the …rst kind is de…ned by the rising factorial given below (cf. [START_REF] Araci | On weighted q-Daehee polynomials with their applications[END_REF])

x n = x (x + 1) (x + n 1) = n X k=0 S 1 (n; k) x k ; (2.20)
which satis…es the following recurrence relation

S 1 (n + 1; k) = nS 1 (n; k) + S 1 (n; k 1) (2.21)
and the following relations for n being a natural number. We now give the following theorem.

Theorem 4. For n 1, we have

n Y k=1 k = 1 n! n X k=0 S 1 (n + 1; k + 1) k ( ) n k : (2.23)
Proof. We prove this theorem by induction: For n = 1, we get

( ) = 1 1! 1 X k=0 S 1 (2; k + 1) k ( ) 1 k = S 1 (2; 1) ( ) + S 1 (2; 2) ;
which is correct because of the values S 1 (2; 1) = 1 = S 1 (2; 2). We assume that the formula in To prove for n = m + 1, we compute that

m+1 Y k=1 k = m + 1 m Y k=1 k = (m + 1) m + 1 1 m! m X k=0 S 1 (m + 1; k + 1) k ( ) m k = 1 (m + 1)! m X k=0 S 1 (m + 1; k + 1) k+1 ( ) m k + 1 m! m X k=0 S 1 (m + 1; k + 1) k ( ) m+1 k = S 1 (m + 1; m + 1) m+1 (m + 1)! + 1 m! S 1 (m + 1; 1) ( ) m+1 + 1 (m + 1)! m 1 X k=0 S 1 (m + 1; k + 1) k+1 ( ) m k + 1 m! m X k=1 S 1 (m + 1; k + 1) k ( ) m+1 k = m+1 (m + 1)! + ( ) m+1 + 1 (m + 1)! m X k=1 S 1 (m + 1; k) k ( ) m+1 k + 1 m! m X k=1 S 1 (m + 1; k + 1) k ( ) m+1 k = m+1 (m + 1)! + ( ) m+1 + 1 (m + 1)! m X k=1 [S 1 (m + 1; k) + (m + 1) S 1 (m + 1; k + 1)] k ( ) m+1 k = m+1 (m + 1)! + ( ) m+1 + 1 (m + 1)! m X k=1 S 1 (m + 2; k + 1) k ( ) m+1 k = 1 (m + 1)! m+1 X k=1 S 1 (m + 2; k + 1) k ( ) m+1 k ;
where we use the formulas (2.21) and (2.22). Therefore, the proof is completed by induction.

By means of Theorem 4 and formula (2.19), we give

( ) x; = x 1 X k=0 S 1 (x; k + 1) k+1 ( ) x 1 k : (2.24)
Using (2.24), we rewrite (2.6) and (2.12) with explicit representations as p x ( ; ) = e (1)

x!

x 1 X k=0 S 1 (x; k + 1) k+1 ( )

x 1 k = 1 X m=0 x 1 X k=0 ( ) m;
m!x! S 1 (x; k + 1) k+1 ( ) x n x!

x 1 X k=0 S 1 (x; k + 1) k+1 ( )

x 1 k = 1 X m=0 1 X x=0 x 1 X k=0 ( ) m;
m!x! x n S 1 (x; k + 1) k+1 ( ) x 1 k :

Conclusion

In the present paper, the degenerate Poisson distribution covering the degenerate exponential function and the -falling factorial have been de…ned and investigated. Various properties of the degenerate Poisson distribution such as the …rst to the fourth raw moments and di¤erence rule have been provided with their proofs. The corresponding moment generating function has been investigated, by which the degenerate Bell polynomials have been considered and their properties have been developed. The results obtained in this paper reduce to the known results for the Poisson distribution and Bell polynomials as approaches zero.

x and e 1

 1 (t) = e (t) . (2.2) Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 20 August 2019 doi:10.20944/preprints201908.0214.v1

  (www.preprints.org) | NOT PEER-REVIEWED | Posted: 20 August 2019 doi:10.20944/preprints201908.0214.v1

S 1 (

 1 n; n) = 1, S 1 (n + 1; 1) = n! and S 1 (n; n 1) = n (

S 1

 1 (2.23) holds for n = m, namely m (m + 1; k + 1) k ( ) m k : Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 20 August 2019 doi:10.20944/preprints201908.0214.v1

preprints.org) | NOT PEER-REVIEWED | Posted: 20 August 2019 doi:10.20944/preprints201908.0214.v1

  In order to get more properties about the degenerate Poisson distribution, we consider that

	( ) x;

14) and e x (t) = te x (t) : (2.15) From (2.14) and (2.15), we obtain e (1) = e (1) and ( ) x; = x ( ) x 1; . (2.16)

We give a product rule of the di¤erence operator in a special case as follows.

Preprints (www.

Con ‡ict of Interest: The authors declare that they have no con ‡ict of interest.

Lemma 1. For m; x being two natural numbers, we then have ( ) m; ( ) x; = m ( ) m 1; ( ) x; + x ( + ) m; ( ) x 1; .

(2.17)

Proof. In terms of (2.13) and (2.16), we attain

which is the claimed result (2.17).

We now improve the properties of the degenerate Poisson distribution. We give the following theorem. Proof. By means of the De…nition 2, the de…nition of the di¤erence operator (2.13) and Lemma 1, we investigate

(m 1)! + ( ) x 1;

(x 2)!