The Inﬂuence of the Magnetic Field Line Curvature on Wall Erosion near the Hall Thruster Exit Plane

: One of the main factors that limit the lifetime of the Hall effect Thrusters (HETs) is the erosion of the acceleration channel caused by the ﬂux of energetic ions. The magnetic ﬁeld that is curved and convex towards the anode has been widely used in HETs because of its role in reducing the divergence angle of the ion beam and the channel wall erosion. However, the mechanism of the inﬂuence of the magnetic ﬁeld line curvature on the wall erosion is not clear. Therefore, in this paper, a 2D3V numerical model based on the immersed-ﬁnite-element and particle-in-cell (IFE-PIC) method is established to simulate the radial-azimuthal plane near the exit of the Hall thruster. The effect of the tilt angle of the magnetic ﬁeld line on the wall sputtering erosion rate is analyzed. The results show that compared to the case with the electric ﬁeld E perpendicular to the magnetic ﬁeld B , the energy of the ions hitting the channel wall is smaller and the wall erosion is weaker when the magnetic ﬁeld lines are convex to the anode. As the tilt angle of the magnetic ﬁeld lines increases from 0 ◦ to 60 ◦ , the erosion rate is reduced by 90%. Conversely, when the magnetic ﬁeld lines are convex to the exit plane of the channel, the wall erosion is much more serious compared to the case with the orthogonal electric ﬁeld E and the magnetic ﬁeld B . As the tilt angle of the magnetic ﬁeld line changes from 0 ◦ to 60 ◦ , the erosion rate is enhanced by 171%. The results in this paper are instructive for the design and optimization of the magnetic ﬁeld of the HETs.


Introduction
Hall-effect thrusters (HETs) have become the most widely used type of thruster device due to their combined advantages of high thrust density, high specific impulse, system simplicity, and high reliability [1,2].Typical HETs consist of an annular ceramic channel.The cross-section and working schematic of the HETs are shown in Figure 1.The anode located upstream in the channel works as the neutral propellant distributor and the external cathode placed in the near-field plume region works as an electron emitter.A large potential difference was applied between the anode and the hollow cathode, creating an electric field along the axial direction that can accelerate the ions produced by ionization and generate thrust.A magnetic field dominated by the radial component along the channel can trap electrons and increase their residence time within the discharge to promote the efficient ionization of neutral propellant gas [3].With the complexity of spacecraft missions escalating, the service life of the thruster devices becomes more demanding.Therefore, as one of the most important factors limiting the life of HETs, the sputtering wall erosion caused by the flux of energetic ions can lead to the mass loss of the discharge channel and removes HETs from the working order.Much experimental and theoretical research has been carried out to investigate the mechanism of sputtering wall erosion [4][5][6][7].
and theoretical research has been carried out to investigate the mechanism of spu wall erosion [4][5][6][7].The magnetic field topology has been found to be a central design issue for because the main working principle of HETs is to use the orthogonal electric and ma fields to sustain the ionization process and accelerate ions.Currently, the commonl magnetic field configuration is the type convex to the anode proposed by Morozo [8] for its proven ability to reduce the divergence angle of the ion beam and impro discharge performance.The characteristics of the configuration are that the main c nent is along the radial direction and the magnetic field lines near the wall are dens that in the center of the channel, forming the so-called "plasma lens" [9,10].A num studies have been carried out to look for insight into this lens-type magnetic field uration.Hofer et al. [11,12] investigated the effect of the magnetic field on the plu vergence angle and plasma oscillations by controlling the magnetic field configu with additional trim coils.It was found that the curvature of the magnetic field li axial gradient of the main radial magnetic field had a great influence on the perfor of the thruster.Fruchtman et al. [13] theoretically analyzed the relationship betwee netic focusing and the plume divergence angle and confirmed the focusing effect curved magnetic field line on the plasma.Garrigues [14] and Liu [15] also confirm effect of the magnetic field topology on the process of sputtering erosion by study subject of ion divergence and wall erosion conditions under the influence of di magnetic field configurations.Liu et al. [16] studied the ion bombardment and w sion in a Hall thruster with magnetic lens design using the PIC method and foun the wall erosion rate as well as the length of the erosion zone decrease monotonical the increase in the curvature of the lens-type magnetic field.
These previous results suggest that the curved magnetic field configuration is h in converging the ion beam and reducing the sputtering erosion caused by energet The magnetic field topology has been found to be a central design issue for HETs because the main working principle of HETs is to use the orthogonal electric and magnetic fields to sustain the ionization process and accelerate ions.Currently, the commonly used magnetic field configuration is the type convex to the anode proposed by Morozov et al. [8] for its proven ability to reduce the divergence angle of the ion beam and improve the discharge performance.The characteristics of the configuration are that the main component is along the radial direction and the magnetic field lines near the wall are denser than that in the center of the channel, forming the so-called "plasma lens" [9,10].A number of studies have been carried out to look for insight into this lens-type magnetic field configuration.Hofer et al. [11,12] investigated the effect of the magnetic field on the plume divergence angle and plasma oscillations by controlling the magnetic field configuration with additional trim coils.It was found that the curvature of the magnetic field line and axial gradient of the main radial magnetic field had a great influence on the performance of the thruster.Fruchtman et al. [13] theoretically analyzed the relationship between magnetic focusing and the plume divergence angle and confirmed the focusing effect of the curved magnetic field line on the plasma.Garrigues [14] and Liu [15] also confirmed the effect of the magnetic field topology on the process of sputtering erosion by studying the subject of ion divergence and wall erosion conditions under the influence of different magnetic field configurations.Liu et al. [16] studied the ion bombardment and wall erosion in a Hall thruster with magnetic lens design using the PIC method and found that the wall erosion rate as well as the length of the erosion zone decrease monotonically with the increase in the curvature of the lens-type magnetic field.
These previous results suggest that the curved magnetic field configuration is helpful in converging the ion beam and reducing the sputtering erosion caused by energetic ions.However, most of the previous studies were conducted using several fixed magnetic field configurations, which failed to reach a universal conclusion on the influence of magnetic field line tilt angle on sputtering erosion.Moreover, the effects of the magnetic field line tilt angle on the sheath formation, the particle motion near the wall, and the wall sputtering process have not been well analyzed, and the influencing mechanism is also unclear.Therefore, in this paper, a 2D3V IFE-PIC numerical model was established to simulate the near-exit radial-azimuthal (r-θ) plane of HETs.By changing the direction and magnitude of the magnetic field axial component B z , different tilt angles of the magnetic field line are realized.The wall erosion rate and relevant parameters are analyzed under different tilt angles.The study presented in this paper has important implications for the design and optimization of the magnetic field configuration.
The organization of the paper is as follows.In Section 2, the relevant model selection, magnetic field setup, and wall erosion calculation method are presented.In Section 3, the ion beam properties and wall erosion under different magnetic field line tilt angles are compared and the associated physical mechanisms are analyzed.The paper is then summarized in Section 4.

PIC Simulation Model
In this paper, a 2D3V IFE-PIC numerical model is established to simulate the radialazimuthal plane near the exit of the Hall thruster.The model is based on the Cartesian coordinate system, as shown in Figure 2. The simulated domain consists of the dielectric wall and a portion of the discharge channel.Its upper part is set to be quasi-neutral, corresponding to the quasi-neutral mainstream area of the discharge channel, and the lower part is considered BN dielectric wall.For clarity of presentation, the upper surface of the wall is taken as the starting of the radial direction of the simulated region, which is denoted as y = 0.The size of the simulation where λ De denotes the Debye length.The model ignores the curvature of the wall in the azimuth and has a finite length in both Ox and Oz directions.A magnetic field B y along the radial direction Oy and a constant electric field E z along the axial direction Oz are applied in the numerical model.Thus, the electrons will drift in the azimuthal direction of Ox (i.e., E x × B y ).tilt angle on the sheath formation, the particle motion near the wall, and the wall sputtering process have not been well analyzed, and the influencing mechanism is also unclear.Therefore, in this paper, a 2D3V IFE-PIC numerical model was established to simulate the near-exit radial-azimuthal (r-θ) plane of HETs.By changing the direction and magnitude of the magnetic field axial component Bz, different tilt angles of the magnetic field line are realized.The wall erosion rate and relevant parameters are analyzed under different tilt angles.The study presented in this paper has important implications for the design and optimization of the magnetic field configuration.The organization of the paper is as follows.In Section 2, the relevant model selection, magnetic field setup, and wall erosion calculation method are presented.In Section 3, the ion beam properties and wall erosion under different magnetic field line tilt angles are compared and the associated physical mechanisms are analyzed.The paper is then summarized in Section 4.

PIC Simulation Model
In this paper, a 2D3V IFE-PIC numerical model is established to simulate the radialazimuthal plane near the exit of the Hall thruster.The model is based on the Cartesian coordinate system, as shown in Figure 2. The simulated domain consists of the dielectric wall and a portion of the discharge channel.Its upper part is set to be quasi-neutral, corresponding to the quasi-neutral mainstream area of the discharge channel, and the lower part is considered BN dielectric wall.For clarity of presentation, the upper surface of the wall is taken as the starting of the radial direction of the simulated region, which is denoted as y = 0.The size of the simulation domain is , where λDe denotes the Debye length.The model ignores the curvature of the wall in the azimuth and has a finite length in both Ox and Oz directions.A magnetic field By along the radial direction Oy and a constant electric field Ez along the axial direction Oz are applied in the numerical model.Thus, the electrons will drift in the azimuthal direction of Ox (i.e., Ex × By).In the numerical calculation, only singly charged electrons and ions are simulated.The secondary electron emission (SEE) and the collisions between particles are not considered.The propellant type used here is argon.A quasi-neutral region of width 5λDe is set near the incident boundary [17].At each time step, a fixed number of ions (determined by the ion density n0 of the quasi-neutral region) enter the simulation region from the incident boundary, and then the numbers of electrons and ions in the quasi-neutral region are counted.If the number of electrons is less than that of the ions, then the corresponding electrons are replenished to make them consistent.If the electron number is greater than In the numerical calculation, only singly charged electrons and ions are simulated.The secondary electron emission (SEE) and the collisions between particles are not considered.The propellant type used here is argon.A quasi-neutral region of width 5λ De is set near the incident boundary [17].At each time step, a fixed number of ions (determined by the ion density n 0 of the quasi-neutral region) enter the simulation region from the incident boundary, and then the numbers of electrons and ions in the quasi-neutral region are counted.If the number of electrons is less than that of the ions, then the corresponding electrons are replenished to make them consistent.If the electron number is greater than that of the ions, no replenishment is required.In addition, the velocity of the incident electrons obeys a Maxwellian distribution, while the ions are considered to be cold, and the incident velocity is chosen equal to the Bohm velocity.
The electrons and ions are treated as PIC particles in the simulation area.The position and velocity of the particles are tracked in all three directions, while the electric field is calculated only in the (Ox-Oy) plane by solving Poisson's equation, which can be expressed as where ε is the vacuum permittivity, Φ is the electric potential, e is the elementary charge, and n e and n i are the density of electrons and ions.The charge densities n i and n e are calculated by weighting all particles to the 2D radial-azimuthal (Ox-Oy) simulation area.
The movements of particles are determined from the integration of the equation of motion, which can be described as [18]: where m and v are the mass and velocity of the particles, B is the magnetic field, and E is the electric field determined by E = −∇Φ .Periodic boundary conditions are applied for Poisson's equation and for the particles at the azimuthal boundaries along the Ox direction.At the positions of x = 0 and x = L x , the electric potential satisfies.If the particles exceed the boundaries in this direction, they are allowed to re-inject into the simulated region from the opposite side.The potential at the incident boundary (i.e., y = L y ) is supposed to be ϕ = 0.The electrons and ions that run out of this boundary are removed from the simulation domain.At the position of y = 0, the dielectric wall is adopted.The particles that hit the wall are removed but their charge will be accumulated on the wall surface [19], and the wall potential can be calculated using the IFE method based on the charge density on the wall.Poisson's equation is not solved in the Oz direction.When particles run out of the boundary along the Oz direction, they will be re-injected from the opposite side while keeping the position and velocity in the (Ox-Oy) plane unchanged.However, to prevent particles from being infinitely accelerated along the Oz direction, the velocities of the re-injected electrons need to be redistributed with the initial temperature T e , and the ion velocities need to be reset as 0 because they are considered to be cold.
In order to maintain stability and reduce numerical heating, the PIC simulations require a time step that satisfies [20]: where ω pe = (e 2 n e )/(ε 0 m) is the electron plasma frequency.At the same time, the grid spacing also needs to be limited as follows [20]: where λ De = (ε 0 T e )/(|e|n e ) is the Debye length with T e the electron temperature.In order to reduce numerical heating and statistical noise, PIC simulations usually require approximately 100 particles per cell.The parameters involved in this paper are given in Table 1.The reliability and accuracy of the model used in this paper have been verified in Ref. [21].

Model Settings
The commonly used magnetic field configuration of Hall thrusters is the type convex to the anode.For example, a typical magnetic field configuration applied to spt100 is shown in Figure 3a, from which it can be seen that the magnetic field lines are curved in the channel and almost perpendicular to the wall near the exit.In order to simplify the calculation, it is assumed that the magnetic field line is axially symmetric along the center line of the channel, as was performed in Ref. [16].The angle α between the lines and the thruster's radial direction can be introduced to describe the characteristics of the magnetic field configuration near the channel exit, which can be simply expressed by:

Model Settings
The commonly used magnetic field configuration of Hall thrusters is the type convex to the anode.For example, a typical magnetic field configuration applied to spt100 is shown in Figure 3a, from which it can be seen that the magnetic field lines are curved in the channel and almost perpendicular to the wall near the exit.In order to simplify the calculation, it is assumed that the magnetic field line is axially symmetric along the center line of the channel, as was performed in Ref. [16].The angle α between the lines and the thruster's radial direction can be introduced to describe the characteristics of the magnetic field configuration near the channel exit, which can be simply expressed by: arctan( ) The α is considered to be zero when the magnetic field line is perpendicular to the wall (i.e., E⊥B), α < 0 when the magnetic field line is convex to the anode, and α > 0 when the magnetic field lines convex to the exit plane of the channel.The schematic of the simplified magnetic field configuration and the angle of α is presented in Figure 3b.In order to ensure that the ECDI caused by the radial magnetic field and axial electric field (Ez × By) does not affect the simulation results as the magnetic field changes, the axial electric field Ez and the radial magnetic field By remain unchanged in all cases.Thus, the magnetic field B is only determined by the magnitude and direction of the axial magnetic field component Bz, as shown as follows: The α is considered to be zero when the magnetic field line is perpendicular to the wall (i.e., E⊥B), α < 0 when the magnetic field line is convex to the anode, and α > 0 when the magnetic field lines convex to the exit plane of the channel.The schematic of the simplified magnetic field configuration and the angle of α is presented in Figure 3b.
In order to ensure that the ECDI caused by the radial magnetic field and axial electric field (E z × B y ) does not affect the simulation results as the magnetic field changes, the axial electric field E z and the radial magnetic field B y remain unchanged in all cases.Thus, the magnetic field B is only determined by the magnitude and direction of the axial magnetic field component B z , as shown as follows: By changing the direction and value of the axial magnetic field B z , the value of α is set as ±15 • , ±30 • , ±45 • , and ±60 • , respectively.Once the value of α is determined, the magnetic field keeps unchanged during the numerical simulation.Because the parameter variations on the outer wall are the same as those on the inner wall, only the results on the inner wall are presented in this paper [16].

Sputtering Model
The erosion of the channel wall primarily comes from the impingement of the energetic ions.The erosion rate q of the channel wall is related to the sputtering yield of the ceramic wall as well as the velocity and angle distribution of the incident particles, and it can be described as [22]: where J i is the ion incident flux density, N is the atomic density of the sputtered wall material, θ i is the angle between the ion incident direction and the surface normal direction, and Y(E ni , θ i ) is the sputtering yield.For the erosion of the channel wall, the sputtering yield represents the number of eroded particles per incident ion at energy E ni and at normal incidence θ i , and the influence of E ni and θ i is independent [23].Therefore, the expression of Y(E ni , θ i ) can be written as [24]: Here, S(E ni ) is the ion energy sputtering coefficient, and it is related to the sputtering threshold E th of the wall material.In this paper, the value of E th is set as 30 eV and the universal semi-empirical formula for S(E ni ) based on the Sigmand theory is adopted [24,25].Y (θ i ) is the ion incident angle-dependent sputtering yield.It is based on a Yamamura empirical formula, which is summarized by experimental results on a large number of materials sputtering [26].
In the simulation, the change of wall morphology caused by the sputtering erosion is not considered.All the results presented in this paper are obtained after convergence.The time at which the simulation reaches the steady state is denoted as t = 0, and all statistics mentioned below start from this moment.

Wall Erosion When the Value of α Is 0
When the electric field E is perpendicular to the magnetic field B, the values of the tilt angle α and the axial magnetic field component B z are 0. In this orthogonal-field configuration, the electrons form a stable E × B drift in the azimuthal direction.The drift velocity can be expressed as: [27]: In HETs, these drift velocities of electrons are as high as 10 6 m/s [20].However, the ions in the channel can be considered unmagnetized due to their large mass relative to the electrons, and their drift velocity along the azimuth direction is almost zero.The large velocity difference between the electrons and the ions produces a strong instability that displays ion acoustic-like characteristics along the azimuthal direction, which is also known as electron cyclotron drift instability (ECDI).The ECDI is a short wavelength, highfrequency kinetic instability that manifests as large potential and density oscillations [27].
As shown in Figure 4, obvious periodic modulations of electric potential and electron density driven by ECDI are observed within the quasi-neutral region of the near-sheath layer.The wavelength and the frequency of the instability are 1.4 mm and 57 MHz, respectively.Moreover, the sheath with an electric potential of −37 V is formed on the wall in Figure 4 because the electrons move faster than the ions.Moreover, the sheath with an electric potential of −37 V is formed on the wall in Figure 4 because the electrons move faster than the ions.Due to the potential gradient between the azimuthal periodic modulations, a local electric field formed in the Ox direction.As shown in Figure 5a, the value of the electric field along the azimuth direction alternates between positive and negative.Therefore, this local electric field is defined as  x E , where "~" indicates that the electric field is oscillating.Figure 5b shows the distributions of the electric field along the azimuthal direction at the distances of 20λDe and 40λDe from the wall.It can be seen that there is a phase difference between the two electric fields due to the inclination angle of the wave structure of ECDI along the traveling direction.It should be noted that the average value of   For HETs, the presence of the near-wall sheath increases the velocity of ions to the wall and prevents the approach of electrons, which induces the radial divergence of the ion beam and intensifies the wall erosion.By collecting the sputtering parameters of energetic ions that hit the wall at different moments, the erosion rate can be obtained.Figure 6 shows the sputtering erosion rate along the azimuthal direction in different time periods with the orthogonal electric fields E and magnetic fields B. Here, the erosion rate was Due to the potential gradient between the azimuthal periodic modulations, a local electric field formed in the Ox direction.As shown in Figure 5a, the value of the electric field along the azimuth direction alternates between positive and negative.Therefore, this local electric field is defined as E x , where "~" indicates that the electric field is oscillating.Figure 5b shows the distributions of the electric field along the azimuthal direction at the distances of 20λ De and 40λ De from the wall.It can be seen that there is a phase difference between the two electric fields due to the inclination angle of the wave structure of ECDI along the traveling direction.It should be noted that the average value of E x is not zero, and the average value of the electric field on the two lines is very close and both are approximately 1.2 × 10 3 V/m.Moreover, the sheath with an electric potential of −37 V is formed on the wall in Figure 4 because the electrons move faster than the ions.Due to the potential gradient between the azimuthal periodic modulations, a local electric field formed in the Ox direction.As shown in Figure 5a, the value of the electric field along the azimuth direction alternates between positive and negative.Therefore, this local electric field is defined as  x E , where "~" indicates that the electric field is oscillating.Figure 5b shows the distributions of the electric field along the azimuthal direction at the distances of 20λDe and 40λDe from the wall.It can be seen that there is a phase difference between the two electric fields due to the inclination angle of the wave structure of ECDI along the traveling direction.It should be noted that the average value of   For HETs, the presence of the near-wall sheath increases the velocity of ions to the wall and prevents the approach of electrons, which induces the radial divergence of the ion beam and intensifies the wall erosion.By collecting the sputtering parameters of energetic ions that hit the wall at different moments, the erosion rate can be obtained.Figure 6 shows the sputtering erosion rate along the azimuthal direction in different time periods with the orthogonal electric fields E and magnetic fields B. Here, the erosion rate was For HETs, the presence of the near-wall sheath increases the velocity of ions to the wall and prevents the approach of electrons, which induces the radial divergence of the ion beam and intensifies the wall erosion.By collecting the sputtering parameters of energetic ions that hit the wall at different moments, the erosion rate can be obtained.Figure 6 shows the sputtering erosion rate along the azimuthal direction in different time periods with the orthogonal electric fields E and magnetic fields B. Here, the erosion rate was averaged in ∆T = 3.2 µs for accuracy.It can be found from Figure 6 that the sputtering erosion rate of the wall was basically stable after the simulation had converged, and the average value of the sputtering erosion rate is approximately 2.9 × 10 −10 m/s.
Sci. 2023, 13, x FOR PEER REVIEW averaged in ΔT = 3.2 μs for accuracy.It can be found from Figu erosion rate of the wall was basically stable after the simulation average value of the sputtering erosion rate is approximately 2.9 0.0 x (mm) Erosion rate (m/s)

Wall Erosion When the Value of α Is Positive
When the value of α is positive, the magnetic field line is con will lead to a change in the electric potential distribution compa electric field E perpendicular to the magnetic field B. Since By and unchanged in all cases, the periodic instability structure caused by region can always be observed and only the wavelength in the sponding to different α varies, as shown in Figure 7.When α = 15° and exerts little impact on the ECDI.Thus, the wavelength in th same as the situation where the electric field E is perpendicular However, with the increase in α, the wavelength of ECDI decrease seen by comparing Figures 4 and 7, the potential on the wall with that of the case with E⊥B, and the wall potential decreases as the

Wall Erosion When the Value of α Is Positive
When the value of α is positive, the magnetic field line is convex to the anode, which will lead to a change in the electric potential distribution compared to the case with the electric field E perpendicular to the magnetic field B. Since B y and E z are preset and remain unchanged in all cases, the periodic instability structure caused by B y × E z in the simulated region can always be observed and only the wavelength in the simulation area corresponding to different α varies, as shown in Figure 7.When α = 15 • , the value of B z is small and exerts little impact on the ECDI.Thus, the wavelength in the simulated area is the same as the situation where the electric field E is perpendicular to the magnetic field B. However, with the increase in α, the wavelength of ECDI decreases.In addition, as can be seen by comparing Figures 4 and 7, the potential on the wall with positive α is lower than that of the case with E⊥B, and the wall potential decreases as the value of α increases.The behavior of electrons could be the explanation for the above phenomenon.The periodic distribution of the potential driven by ECDI leads to an electric field   E B will be generated.The new orthogonal field induces a new drifting motion perpendicular to the wall surface.If the value of α is positive, the new drifting motion of the electrons is towards the wall.This movement increases the velocity of the electrons tending to the wall.Therefore, the electron flux that can reach the wall becomes larger, and the potential drop in the sheath increases.As the value of the positive α increases, the radial drift velocity increases.More electrons can reach the wall and make the wall potential decreases.Therefore, it can be seen from the electron density distribution shown in Figure 8 that as the value of positive α increases, the wall potential decreases and the low-density region of the electrons near the wall expands.The behavior of electrons could be the explanation for the above phenomenon.The periodic distribution of the potential driven by ECDI leads to an electric field E x along the azimuthal direction, as shown in Figure 5.The average value of E x is non-zero and positive.When the magnetic field component B z is applied in the axial direction, a new orthogonal field E x × B z will be generated.The new orthogonal field induces a new drifting motion perpendicular to the wall surface.If the value of α is positive, the new drifting motion of the electrons is towards the wall.This movement increases the velocity of the electrons tending to the wall.Therefore, the electron flux that can reach the wall becomes larger, and the potential drop in the sheath increases.As the value of the positive α increases, the radial drift velocity increases.More electrons can reach the wall and make the wall potential decreases.Therefore, it can be seen from the electron density distribution shown in Figure 8 that as the value of positive α increases, the wall potential decreases and the low-density region of the electrons near the wall expands.Although the axial component of the curved magnetic field does not directly affect the ion distribution, the ion behavior will be affected by the electric potential dominated by the electron motion.When the electric field E is perpendicular to the magnetic field B, most of the ions are captured by the "ion trapping" phenomenon induced by ECDI and move primarily in the azimuthal direction within the channel.When the magnetic field line is convex to the exit plane of HETs, the potential of the wall decreases.The potential gradient along the radial direction allows a greater acceleration of the ions moving toward the wall.Figure 9 shows the ion velocity distribution function for the radial component with different values of the positive α.In order to eliminate the interference of the quasineutral region on the particle velocity, all ions within 10λDe from the upper boundary are neglected in the statistics.It can be seen that as the value of the positive α increases, the average velocity of ions toward the wall increases, and the distribution of the ion velocity becomes more uniform.Although the axial component of the curved magnetic field does not directly affect the ion distribution, the ion behavior will be affected by the electric potential dominated by the electron motion.When the electric field E is perpendicular to the magnetic field B, most of the ions are captured by the "ion trapping" phenomenon induced by ECDI and move primarily in the azimuthal direction within the channel.When the magnetic field line is convex to the exit plane of HETs, the potential of the wall decreases.The potential gradient along the radial direction allows a greater acceleration of the ions moving toward the wall.Figure 9 shows the ion velocity distribution function for the radial component with different values of the positive α.In order to eliminate the interference of the quasi-neutral region on the particle velocity, all ions within 10λ De from the upper boundary are neglected in the statistics.It can be seen that as the value of the positive α increases, the average velocity of ions toward the wall increases, and the distribution of the ion velocity becomes more uniform.Figure 10 shows the average sputtering erosion rate of the wall corresponding to different positive α, which is obtained by averaging the erosion rate along the Ox direction after simulation convergence.It can be seen that compared with the case in which the electric field E is perpendicular to the magnetic field B, the applied axial magnetic field component Bz results in a larger wall erosion rate.The increasing potential drop along the radial direction promotes the divergence of the ion beam in the discharge channel and increases the energy of the ions hitting the wall.Thus, the sputtering erosion of the ions on the wall is aggravated.Moreover, with the increasing value of the positive α, the radial velocity of the ions increases, which makes the wall erosion more serious.Compared with the case with the electric field E perpendicular to the magnetic field B, when the value of the angle α is 15°, 30°, 45°, and 60°, the wall erosion rate is enhanced by 4.9%, 31%, 70%, and 171%, respectively.

α (°)
Erosion rate (m/s) Figure 10 shows the average sputtering erosion rate of the wall corresponding to different positive α, which is obtained by averaging the erosion rate along the Ox direction after simulation convergence.It can be seen that compared with the case in which the electric field E is perpendicular to the magnetic field B, the applied axial magnetic field component B z results in a larger wall erosion rate.The increasing potential drop along the radial direction promotes the divergence of the ion beam in the discharge channel and increases the energy of the ions hitting the wall.Thus, the sputtering erosion of the ions on the wall is aggravated.Moreover, with the increasing value of the positive α, the radial velocity of the ions increases, which makes the wall erosion more serious.Compared with the case with the electric field E perpendicular to the magnetic field B, when the value of the angle α is 15 • , 30 • , 45 • , and 60 • , the wall erosion rate is enhanced by 4.9%, 31%, 70%, and 171%, respectively.Figure 10 shows the average sputtering erosion rate of the wall corresponding to different positive α, which is obtained by averaging the erosion rate along the Ox direction after simulation convergence.It can be seen that compared with the case in which the electric field E is perpendicular to the magnetic field B, the applied axial magnetic field component Bz results in a larger wall erosion rate.The increasing potential drop along the radial direction promotes the divergence of the ion beam in the discharge channel and increases the energy of the ions hitting the wall.Thus, the sputtering erosion of the ions on the wall is aggravated.Moreover, with the increasing value of the positive α, the radial velocity of the ions increases, which makes the wall erosion more serious.Compared with the case with the electric field E perpendicular to the magnetic field B, when the value of the angle α is 15°, 30°, 45°, and 60°, the wall erosion rate is enhanced by 4.9%, 31%, 70%, and 171%, respectively.

Wall Erosion under Negative α
When the value of α is negative, the magnetic field line is convex to the channel exit.Similar to what is presented in Section 3.2, the drift formed by the electromagnetic field E x × B z causes the electron to have a radial velocity.The difference is that the radial drift velocity here is away from the wall.Under the influence of this radial drift, some electrons with small incident velocity may directly leave the simulation area from the incident boundary, resulting in the reduction of the electron density within the simulation region.Another section of electrons that used to be able to enter the simulation area is forced to slow down and stay near the entrance, which leads to an accumulation of electrons near the entrance.The difference in the behavior of electrons causes the potential distribution under negative α to show a significant difference from the condition of α > 0. As shown in Figure 11, the wavelength of ECDI becomes shorter due to the reduction of electron density in the simulation area.In addition, the potentials of the dielectric wall are generally higher than that in Figure 7 and even local positive potentials appear in the simulation area.Moreover, as the value of α decreases, both the wall potential and the local positive potential generated by the ion accumulation in the simulation area increase.

Wall Erosion under Negative α
When the value of α is negative, the magnetic field line is convex to the channel exit.Similar to what is presented in Section 3.2, the drift formed by the electromagnetic field   E B causes the electron to have a radial velocity.The difference is that the radial drift velocity here is away from the wall.Under the influence of this radial drift, some electrons with small incident velocity may directly leave the simulation area from the incident boundary, resulting in the reduction of the electron density within the simulation region.Another section of electrons that used to be able to enter the simulation area is forced to slow down and stay near the entrance, which leads to an accumulation of electrons near the entrance.The difference in the behavior of electrons causes the potential distribution under negative α to show a significant difference from the condition of α > 0. As shown in Figure 11, the wavelength of ECDI becomes shorter due to the reduction of electron density in the simulation area.In addition, the potentials of the dielectric wall are generally higher than that in Figure 7 and even local positive potentials appear in the simulation area.Moreover, as the value of α decreases, both the wall potential and the local positive potential generated by the ion accumulation in the simulation area increase.These phenomena can be explained by the more detailed behavior of electrons.As the direction of the radial drift velocity is away from the wall, the movement of the These phenomena can be explained by the more detailed behavior of electrons.As the direction of the radial drift velocity is away from the wall, the movement of the electrons toward the wall slows down or is even reverse accelerated.Therefore, whether electrons will accumulate near the entrance or escape directly from the incident boundary depends on the magnitude of the radial drift velocity.Figure 12 shows the distribution of electron density corresponding to different negative α.It can be seen that there is a large density gradient in the radial direction.When the value of α is −15 • and −30 • , the magnitude of the drift velocity away from the wall is small.It makes it difficult for the electrons to enter the simulation region, so the density accumulation can be observed near the entrance in Figure 12a,b.However, when the value of α is −45 • or −60 • , the drift velocity becomes larger, making it impossible for electrons to remain near the entrance.Most of the electrons escape from the simulation region.As a result, the average electron density near the upper boundary decreases, as shown in Figure 12c,d.In addition, as the value of the negative α decreases, the velocity of the electrons away from the wall increases, resulting in a decreased number of electrons that can reach the wall.Therefore, as shown in Figure 11, the potential of the dielectric wall increases with the decreasing value of the negative α.When the value of α is 60 • , there is even a positive wall potential.However, the ions are able to enter the simulation area for the given incident conditions because they are not bound by the magnetic field, thus resulting in a local positive potential in the simulation area.Meanwhile, as the electron density decreases due to the change in α, a higher value of the positive potential is presented due to the growth of the difference in n e and n i .
Appl.Sci.2023, 13, x FOR PEER REVIEW 13 of 18 electrons toward the wall slows down or is even reverse accelerated.Therefore, whether electrons will accumulate near the entrance or escape directly from the incident boundary depends on the magnitude of the radial drift velocity.Figure 12 shows the distribution of electron density corresponding to different negative α.It can be seen that there is a large density gradient in the radial direction.When the value of α is −15° and −30°, the magnitude of the drift velocity away from the wall is small.It makes it difficult for the electrons to enter the simulation region, so the density accumulation can be observed near the entrance in Figure 12a,b.However, when the value of α is −45° or −60°, the drift velocity becomes larger, making it impossible for electrons to remain near the entrance.Most of the electrons escape from the simulation region.As a result, the average electron density near the upper boundary decreases, as shown in Figure 12c,d.In addition, as the value of the negative α decreases, the velocity of the electrons away from the wall increases, resulting in a decreased number of electrons that can reach the wall.Therefore, as shown in Figure 11, the potential of the dielectric wall increases with the decreasing value of the negative α.When the value of α is 60°, there is even a positive wall potential.However, the ions are able to enter the simulation area for the given incident conditions because they are not bound by the magnetic field, thus resulting in a local positive potential in the simulation area.Meanwhile, as the electron density decreases due to the change in α, a higher value of the positive potential is presented due to the growth of the difference in ne and ni.The electron motion and wall potential distribution have important effects on the behavior of the ions when the value of α is negative.The decrease in the potential drop between the dielectric wall and the quasi-neutral region makes it difficult for ions to reach The electron motion and wall potential distribution have important effects on the behavior of the ions when the value of α is negative.The decrease in the potential drop between the dielectric wall and the quasi-neutral region makes it difficult for ions to reach the wall.More electrons are trapped in the channel, creating local positive potentials, as shown in Figure 11.The resulting potential gradient along the radial direction further prevents the ions from approaching the walls, even causing them to show acceleration in the direction away from the wall.In addition, the ion velocity induced by the "ion trapping" phenomenon of ECDI is also in the direction away from the wall.Therefore, different from the case with positive α shown in Figure 9, some of the ion velocities when α is negative are along the direction far away from the wall, as shown in Figure 13.Moreover, as the value of negative α decreases, the number of the ions whose velocity is away from the wall increases, and the maximum value of the velocity in this direction rises.Thus, the overall distribution of IVDF becomes flat.
the wall.More electrons are trapped in the channel, creating local positive pote shown in Figure 11.The resulting potential gradient along the radial directio prevents the ions from approaching the walls, even causing them to show accele the direction away from the wall.In addition, the ion velocity induced by the " ping" phenomenon of ECDI is also in the direction away from the wall.Therefo ent from the case with positive α shown in Figure 9, some of the ion velocities w negative are along the direction far away from the wall, as shown in Figure 13.M as the value of negative α decreases, the number of the ions whose velocity is aw the wall increases, and the maximum value of the velocity in this direction rises.overall distribution of IVDF becomes flat.The difference in IVDF for radial components shown in Figures 9 and 13 different degrees of sputtering erosion on the dielectric wall.When the value of ative, the ions entering the simulation area are slowed down or even reversely acc resulting in a reduction in both the number and energy of ions that can hit the erosion rate values in Figure 14 are generally smaller than that in the cases of α > 0. Meanwhile, as the value of the negative α decreases, the average velocity o away from the wall increases, and the erosion rate of the dielectric wall induce ion sputtering decreases.Therefore, compared with the cases of E⊥B and α > magnetic field lines convex to the anode direction enables better focusing of the and reduction of the erosion rate.Compared with the case with the electric fie pendicular to the magnetic field B, when the value of the angle α is −15°, −30°, −60°, the erosion rate decreases by 5.5%, 59%, 83%, and 90%, respectively.The difference in IVDF for radial components shown in Figures 9 and 13 results in different degrees of sputtering erosion on the dielectric wall.When the value of α is negative, the ions entering the simulation area are slowed down or even reversely accelerated, resulting in a reduction in both the number and energy of ions that can hit the wall.The erosion rate values in Figure 14 are generally smaller than that in the cases of α = 0 and α > 0. Meanwhile, as the value of the negative α decreases, the average velocity of the ions away from the wall increases, and the erosion rate of the dielectric wall induced by the ion sputtering decreases.Therefore, compared with the cases of E⊥B and α > 0, having magnetic field lines convex to the anode direction enables better focusing of the ion beam and reduction of the erosion rate.Compared with the case with the electric field E perpendicular to the magnetic field B, when the value of the angle α is −15 • , −30 • , −45 • , and −60 • , the erosion rate decreases by 5.5%, 59%, 83%, and 90%, respectively.According to the above conclusions, several relevant phenomena can be explained in the Hall thrusters.Firstly, the magnetic field configuration that is convex toward the anode is always chosen for application in Hall thrusters because this configuration corresponds to the case with negative α in this paper.According to the numerical results, this type of magnetic field configuration benefits the convergence of ion beams and can reduce the sputtering erosion rate on the wall.Secondly, the erosion near the exit of the thruster is more serious than that inside the channel because the value of the negative α near the exit is smaller than that near the anode.According to the result shown in Figure 14, the larger value of α corresponds to a larger erosion rate when α is negative.Therefore, in the channel with a smaller value of the negative α, the effect of the magnetic field on shielding ion sputtering is more effective, and the wall erosion is less severe.

Conclusions
In this paper, a 2D3V IFE-PIC numerical model was established to simulate a portion of the radial-azimuthal plane near the exit of the HETs.The effect of the magnetic field line curvature on the sheath formation and dielectric wall erosion was analyzed under different tilt angles of the magnetic field.
When the electric field E is perpendicular to the magnetic field B, the erosion rate caused by ion sputtering on the dielectric wall remains constant after simulation convergence with a value of 2.9 × 10 −10 m/s.If the magnetic field line is convex to the anode, the number of electrons that can hit the wall increases, which causes the wall potential to decrease.The increased potential gradient between the wall and the quasi-neutral region accelerates the ions towards the wall, leading to a more intense degree of wall erosion compared to the case where the electric field E is perpendicular to the magnetic field B. Furthermore, as α increases, the erosion rate becomes larger, and the value of the erosion rate increased by 4.9%, 31%, 70%, and 171% as α changed from 0° to 15°, 30°, 45°, and 60°, respectively.If the magnetic field line is convex to the exit plane of HETs, the number of electrons that can reach the wall decreases, resulting in an increase in the wall potential.The decreased potential drop along the radial direction reduces the energy of the ions hitting the wall.Therefore, the wall sputtering erosion rate is smaller than that with orthogonal electromagnetic fields.As α varies from 0° to −15°, −30°, −45°, and −60°, the value of the erosion rate decreases by 5.5%, 59%, 83%, and 90%, respectively.According to the above conclusions, several relevant phenomena can be explained in the Hall thrusters.Firstly, the magnetic field configuration that is convex toward the anode is always chosen for application in Hall thrusters because this configuration corresponds to the case with negative α in this paper.According to the numerical results, this type of magnetic field configuration benefits the convergence of ion beams and can reduce the sputtering erosion rate on the wall.Secondly, the erosion near the exit of the thruster is more serious than that inside the channel because the value of the negative α near the exit is smaller than that near the anode.According to the result shown in Figure 14, the larger value of α corresponds to a larger erosion rate when α is negative.Therefore, in the channel with a smaller value of the negative α, the effect of the magnetic field on shielding ion sputtering is more effective, and the wall erosion is less severe.

Conclusions
In this paper, a 2D3V IFE-PIC numerical model was established to simulate a portion of the radial-azimuthal plane near the exit of the HETs.The effect of the magnetic field line curvature on the sheath formation and dielectric wall erosion was analyzed under different tilt angles of the magnetic field.
When the electric field E is perpendicular to the magnetic field B, the erosion rate caused by ion sputtering on the dielectric wall remains constant after simulation convergence with a value of 2.9 × 10 −10 m/s.If the magnetic field line is convex to the anode, the number of electrons that can hit the wall increases, which causes the wall potential to decrease.The increased potential gradient between the wall and the quasi-neutral region accelerates the ions towards the wall, leading to a more intense degree of wall erosion compared to the case where the electric field E is perpendicular to the magnetic field B. Furthermore, as α increases, the erosion rate becomes larger, and the value of the erosion rate increased by 4.9%, 31%, 70%, and 171% as α changed from 0 • to 15 • , 30 • , 45 • , and 60 • , respectively.If the magnetic field line is convex to the exit plane of HETs, the number of electrons that can reach the wall decreases, resulting in an increase in the wall potential.The decreased potential drop along the radial direction reduces the energy of the ions hitting the wall.Therefore, the wall sputtering erosion rate is smaller than that with orthogonal electromagnetic fields.As α varies from 0 • to −15

Figure 1 .
Figure 1.Geometry structure of the PIC model in radial-azimuthal plane.

Figure 1 .
Figure 1.Geometry structure of the PIC model in radial-azimuthal plane.

Figure 2 .
Figure 2. Geometry structure of the PIC model in azimuthal-radial plane.

Figure 2 .
Figure 2. Geometry structure of the PIC model in azimuthal-radial plane.

Figure 3 .
Figure 3. Schematic diagram of (a) typical magnetic field topology of SPT-100; (b) the simplified magnetic field configuration and tilt angle of magnetic field lines.

Figure 3 .
Figure 3. Schematic diagram of (a) typical magnetic field topology of SPT-100; (b) the simplified magnetic field configuration and tilt angle of magnetic field lines.

Figure 4 .
Figure 4.The 2D spatial distribution of (a) electric potential and (b) electron density with E⊥B.

xE
is not zero, and the average value of the electric field on the two lines is very close and both are approximately 1.2 × 10 3 V/m.

Figure 5 .
Figure 5. Distribution of electric field in (a) two-dimensional and (b) one-dimensional space along azimuth at y = 0.84 mm and y = 1.68 mm with E⊥B.

Figure 4 .
Figure 4.The 2D spatial distribution of (a) electric potential and (b) electron density with E⊥B.

Figure 4 .
Figure 4.The 2D spatial distribution of (a) electric potential and (b) electron density with E⊥B.

xE
is not zero, and the average value of the electric field on the two lines is very close and both are approximately 1.2 × 10 3 V/m.

Figure 5 .
Figure 5. Distribution of electric field in (a) two-dimensional and (b) one-dimensional space along azimuth at y = 0.84 mm and y = 1.68 mm with E⊥B.

Figure 5 .
Figure 5. Distribution of electric field in (a) two-dimensional and (b) one-dimensional space along azimuth at y = 0.84 mm and y = 1.68 mm with E⊥B.

Figure 6 .
Figure 6.The erosion rate of the dielectric wall along the azimuthal direc ods with E⊥B.

Figure 6 .
Figure 6.The erosion rate of the dielectric wall along the azimuthal direction in different time periods with E⊥B.

Figure 7 .
Figure 7.The 2D spatial distributions of electric potential with different values of the positive α.

E
azimuthal direction, as shown in Figure 5.The average value of x is non-zero and positive.When the magnetic field component Bz is applied in the axial direction, a new orthogonal field x z

Figure 7 .
Figure 7.The 2D spatial distributions of electric potential with different values of the positive α.

Figure 8 .
Figure 8.The 2D spatial distribution of electron density with different values of the positive α.

Figure 8 .
Figure 8.The 2D spatial distribution of electron density with different values of the positive α.

Figure 9 .
Figure 9. Ion velocity distribution function for radial component with different values of the positive α.

Figure 9 .
Figure 9. Ion velocity distribution function for radial component with different values of the positive α.

Figure 9 .
Figure 9. Ion velocity distribution function for radial component with different values of the positive α.

Figure 10 .
Figure 10.The average erosion rate of the dielectric wall with different values of the positive α.

18 Figure 10 .
Figure 10.The average erosion rate of the dielectric wall with different values of the positive α.
Figure 11.The 2D distributions of electric potential with different values of the negative α.

Figure 11 .
Figure 11.The 2D distributions of electric potential with different values of the negative α.

Figure 12 .
Figure 12.The 2D distributions of electron density with different values of the negative α.

Figure 12 .
Figure 12.The 2D distributions of electron density with different values of the negative α.

Figure 13 .
Figure 13.Ion velocity distribution function for radial component with different values o ative α.

Figure 13 .
Figure 13.Ion velocity distribution function for radial component with different values of the negative α.

Figure 14 .
Figure 14.The averaged erosion rate of the dielectric wall corresponding to different values of the negative α.

Figure 14 .
Figure 14.The averaged erosion rate of the dielectric wall corresponding to different values of the negative α.

Table 1 .
Numerical parameters used in the 2D3V PIC simulation.