The influence of operating parameters on the dynamic characteristics of minimized electron cyclotron resonance ion thrusters for space gravitational wave detection

Abstract

Electron cyclotron resonance ion thrusters (ECRITs) have the potential to be used for space gravitational wave detection due to their wide thrust range. However, an unclear understanding of dynamic processes of ECRITs with strongly coupled multi-operating parameters limits further improvements on thrust noise and response velocity by feedback control systems. An integrative mathematical model considering the non-Maxwell electron energy distribution function for ECRITs is validated by experiments and used to study the influence of operating parameters on the dynamic processes of thrusters, which provides a new simplified grid model. Simulation results show the response processes with microwave (MW) power can be divided into two stages. The characteristic times of the first and second stages are respectively several microseconds and 10 ms, which are respectively dominated by plasma motion and the volume effect. The overshoot of screen grid (SG) current decreases, while its response time remains unchanged when the response time of MW power is prolonged. The response time of SG current with a step increase of flow rate is approximately 10 ms, consistent with the volume effect. The SG current decreases with rise of flow rate for high flow rate operations due to the small increment of ion density limited by low electron temperature, the decrease of ion Bohm velocity and reduction of sheath extraction area. The influence of grid voltage on the dynamic process of the SG current depends on variation ranges of extraction capabilities. When variations of sheath extraction area are limited, the response time is 5 μs, consistent with plasma response time. It is prolonged to 0.5 ms if sheath extraction area variations are large because they cause obvious variations of plasma parameters in the discharge chamber. These dynamic results can not only facilitate designing feedback controllers of micro-propulsion systems for high-precision space missions, but also provide guidance for ion sources to generate highly stable or rapid-response ion beam.

1. Introduction

Rapid advances of astronautic technology give high-precision scientific experiments the potential to be carried out in space due to the advantages of low environment noise and vastness. These high-accuracy tasks put forward strict demands on the stability of satellites, which can be realized by micro-propulsion systems to offset the non-conservative forces on satellites (Maghami et al. 2017). Space gravitational wave (GW) detection is a representative, which aims to detect GWs in the low frequency band by laser interference measuring the optical length variations between two proof masses (Huang et al. 2017). There are two main orbit plans for space GW detection. The first one is the geocentric orbit, which is put forward by the TianQin plan. The detection of the TianQin plan relies on three identical spacecrafts, which are placed on geocentric orbits with a semi-major axis of 10⁵ km and form a nearly equilateral triangle (Luo et al. 2016). The second one is the heliocentric orbit, which is put forward by the laser interferometer space antenna (LISA). It comprises a cluster of three spacecrafts, located at the corners of an equilateral triangle of 5 million kms. The cluster is in a heliocentric orbit trailing the earth by 20° (Danzmann 2000). Although the orbit plans for the TianQin plan and LISA are different, the requirements of thrusters for these two plans are similar because the main non-conservative forces on satellites are analogous. The propulsion system developed for space GW detection needs to satisfy wide thrust range (5∼100 μN), low thrust noise (0.1 μN Hz^−0.5), short response time (50 ms), high thrust resolution (0.1 μN) and long lifetime (10 000 h) (Liu et al. 2022). The thrust response time refers to the time when the thrust is stabilized after its command is adjusted.

Some micro-thrusters have been tested in orbit to validate their feasibility for space GW detection. The colloid micro-Newton thrusters (CMNTs) and cold gas thrusters (CGTs) are launched on LISA Pathfinder. Test results of CMNTs show the thrust range only covers 5∼30 μN and the reliability is seriously weakened by the complex structures. A CMNT in orbit fails due to obstruction of the emitter (Anderson et al. 2018). Compared with CMNTs, the reliabilities of CGTs are highly improved. However, their response processes are too slow and their low specific impulses result in propellant mass burden (Armano et al. 2019). The CGTs in Tianqin-1 show similar problems (Luo et al. 2020). The Hall micro-thrusters (HMTs) and radio-frequency ion thrusters (RITs) are tested on Taiji-1. Test results show the thrust noise of the HMT does not satisfy requirements and is hard to suppress because of its low thrust estimation accuracy (Xu et al. 2021). The thrust noise of the RIT is 0.2 μN Hz^−0.5. The thrust minimum of the RIT is close to the lower bound of requirements, which easily causes flameout when the thrust is adjusted near its lower boundary (He, Duan & Qi 2021). The electron cyclotron resonance ion thrusters (ECRITs) are effective electric thrusters. They utilize microwaves to heat electrons based on electron cyclotron resonance (ECR) and enhance ionization intensity. They accelerate ions by the strong electric field in the grid system and produce thrust. Compared with other thrusters, many advantages of ECRITs give them potential for the space GW detection task. Firstly, the thrust estimation accuracy is high, which means the thrust noise can be effectively suppressed by feedback control. Besides, the ECRITs can work in very low pressure due to the resonance heating mechanism, which means the thrust minimum can be further reduced (Lieberman & Lichtenberg 2005). Furthermore, their specific impulses are high, which can effectively reduce propellant mass. The reliabilities of ECRITs have been validated by many space missions, like the μ10 launched in the Falcon and Falcon II asteroid detectors (Yamamoto et al. 2006; Nishiyama et al. 2020) and the μ1 in HODOYOSHI-4 and PROCYON (Funase, Koizumi & Nakasuka 2014; Hiroki 2015). Faced with space GW detection, a minimized ECRIT is designed, whose thrust range covers 1∼100 μN (Meng, Zhu & Yu 2023). Its thrust noise, resolution and response time are closely related to the design of the feedback control system.

Test results show that thrust noise is higher in the low frequency band. Besides, the thrust response velocity is hard to improve when adjusting the flow rate under an open-loop condition. On one hand, the feedback control is advantageous to reducing thrust noise in the low frequency band. On the other hand, feedback control can improve the thrust response velocity by automatically generating closed-loop microwave power commands while adjusting the flow rate. Therefore, a high-precision multi-parameter coordination feedback control system is crucial for space GW detection, whose design is based on the clear understanding of the dynamic characteristics of ECRITs with operating parameters. This is mainly because these dynamic characteristics, especially the characteristic time and amplitude of the overshoot, can directly affect the loop gain design, stability and performances of a feedback controller. However, the strong coupling of multi-operating parameters for ECRITs makes their dynamic characteristics complex. The flow rate and microwave (MW) power influence not only plasma parameters upstream of the screen grid (SG) and the ion extraction properties, but also ion acceleration characteristics via the SG sheath potential. The grid voltages affect not only the acceleration voltage, but also ion extraction properties via the sheath shape and the extraction area. Experiments are hard to solve this problem. On one hand, the dynamic characteristics of ECRITs are hard to be totally reflected by experiments due to the slow response of the flow rate controller. On the other hand, the mechanism of the dynamic process is closely related to the SG sheath and hard to analyse, limited by diagnostic methods. Compared with experiments, simulation is a cost-effective way to study the physical process of ECRITs. The particle in cell (PIC) model and global model are two common simulation tools. The PIC model has been widely used to study physical processes like the sheath structure (Liu, Niu & Yu 2019) and the interaction between plasma and walls (Roy & Pandey 2002), which is significant for optimization of the thruster (Ding et al. 2017; Fede et al. 2021). Fu et al. have studied plasma response processes with the adjustment of grid voltages by PIC simulation and concluded the characteristic time of ECRITs is 1.5 μs (Fu et al. 2022). The plasma response time refers to the characteristic time of ions moving in the discharge chamber. It is of the order of several micro-seconds according to the size of the discharge chamber and ion velocity. Besides, for the plasma motion, the adjustment of operating parameters may result in the diffusion motion of neutral gas, which is called the volume effect. Its characteristic time is usually much longer than the plasma response time. Therefore, on the basis of Fu's work, the influence of the operating parameters on the dynamic process of ECRITs considering the volume effect is crucial for clear understanding of thruster dynamic characteristics, design of a feedback controller and improvement of thruster performances. Although many simplified ways (Szabo 2001) and parallel strategies have been used (Duras et al. 2017), computational cost for the PIC model is still high, especially considering the response time of the flow rate. Compared with the PIC model, the computation cost of the global model is low. However, it can only simulate average plasma parameters in the discharge chamber and cannot reflect sheath properties and ion beamlet divergence.

In this paper, an integrative model is built based on the 1 cm ECRIT designed for space GW detection. The discharge chamber is simulated by a global model. The sheath and ion beamlet are respectively described by an elliptic paraboloid and rotating hyperboloid. The integrative model can reflect the sheath shape, divergence of ion beamlet and intersection rate variations with operating parameters. Besides, it agrees well with experimental results. Although the model is built based on micro-Newton ECRITs, the modelling methods and conclusions are still applicable to ECRITs with other thrust levels. Only some key parameters like the size of the discharge chamber and the effective loss areas of electrons and ions should be revised. Besides, the simplification method also provides guidance for other ion sources which adopt ion optics to extract and control ion beam. The structure of this paper is shown as follows. In § 2, the research objective, experimental set-up and simulation model are introduced; in § 3, the dynamic results of state parameters and thruster performances with the adjustment of MW powers, flow rates and grid voltages are analysed. The state parameters include electron density, electron temperature, atom density and atom temperature; the performances mainly include SG current, acceleration grid (AG) current and thrust. In § 4, the conclusion is summarized.

2. Simulation model and experimental set-up

2.1. Research objective

In this paper, the research objective is the ECRIT developed by the Harbin Institute of Technology faced with the space GW detection task. The structure and a photo of the thruster are illustrated as figure 1. The propellant Xe is fed into the discharge chamber through a flow tube and a gas distributor. The MW is injected into the thruster through an antenna. The opposite direct current (DC) voltages are applied to the SG and AG, which can form a strong electric field between the grids and accelerate ions to produce thrust. Besides, it adopts two-stage permanent magnets to form the cusped magnetic field and confine the electrons. On one hand, the magnetic configuration can enhance ionization intensity and improve specific impulses. On the other hand, the magnetic configuration can reduce the interaction between plasma and walls, which can increase thruster lifetime. The diameter of the discharge chamber is 1 cm. There are 85 holes of grids. The diameters of SG holes and AG holes are respectively 0.72 mm and 0.4 mm. The thicknesses of SG and AG are respectively 0.2 mm and 0.3 mm. The distance between two grids is 0.25 mm. Experimental results show its thrust range is 1∼100 μN by adjusting operating parameters. The ERCIT should be ignited before the detection task starts. Under the detection mode, the drag-free controller generates thrust commands according to the non-conservative forces on satellites. The operating parameters of ECRITs are automatically and continuously adjusted by the feedback controller.

Figure 1. The illustration of ECRIT. (a) Structure and (b) photo.

2.2. Experimental set-up

Experiments are conducted to obtain necessary data for the simulation model and its validation. The experimental platform is illustrated as figure 2. The vacuum chamber is cuboid, whose size is 0.8 × 0.8 × 0.7 m. Two molecular pumps and one mechanical pump work to maintain the pressure below 10⁻⁴ Pa. The mass rate of propellants is controlled by an Alicat flow meter. The DC sources are developed by HSPY, whose ranges are 0∼1500 V. The MW power sources can realize a 0∼20 W range. These devices satisfy the experimental requirements.

Figure 2. Experimental platform.

There are two experiments which are conducted in the paper. Firstly, the electron energy distribution needs to be tested by a Langmuir probe because it deviates from a Maxwell distribution as it is affected by the complex cusped magnetic field and many inelastic collision processes. The real electron energy distribution is critical for computing the rate coefficients of the ionization and excitation reactions. The Langmuir probe consists of a cylindrical tantalum and a ceramic tube, whose diameter is 0.7 mm. Considering the strongly magnetized plasma in the discharge chamber, the planar probe is utilized to avoid some problems of cylindrical probes like radical sheath expansion and electron non-saturation. The flat end of the Langmuir probe is located at the ECR layer, which is 4.5 mm from the SG axially and 2 mm from the central axis radically. The scanning range of probe voltage is −35–150 V. The electron energy distribution is obtained by differentiating the transition region of the I–V characteristic curve, which is shown in § 2.3. Secondly, we measure SG current and AG current under different operations to validate the model accuracy. These electric parameters are recorded by a YOKOGAWA wave recorder. These experimental results are shown at the beginning of §§ 3.1–3.3.

2.3. Simulation model

The integrative model is based on some assumptions listed as follows. The plasma parameters in the discharge chamber are uniform; the sheath profile upstream of the SG is approximately elliptic paraboloid (Yan et al. 2020); the ion beamlet shape is approximately a rotating hyperboloid considering the Coulomb repulsion effect of ions in the grid system (Shang et al. 2019; Li et al. 2022). The main structure of the model is shown in figure 3. The plasma parameters in the discharge chamber are simulated by a global model. The discharge chamber model is connected with the sheath model by plasma parameters upstream of the SG and grid voltage. The ion beamlet model is connected with the sheath model by sheath profile parameters.

Figure 3. The structure of the model.

It can be seen from figure 3 that the influences of the flow rate and MW power on the ion acceleration characteristics are considered via the sheath potential in this model. Besides, the influence of grid voltage on the ion extraction properties is considered via variations of the sheath shape in the model. Compared with traditional global model, the integrative model built in this paper can not only analyse plasma parameters in the discharge chamber, but also analyse variations of the sheath shape with operating parameters. Besides, this model considers the divergence of the ion beamlet and the intersection effect of the AG, which are helpful for analysing thruster dynamic processes with adjustment of operating parameters. The main equations and parameters can be referenced in Niu, Liu & Yu (2024), which are introduced briefly as follows. The discharge chamber model includes the ion number conservation equation (Lafleur et al. 2022), the atom number conservation equation (Magaldi et al. 2022), the electron energy conservation equation (Lucken et al. 2019; Nakagawa et al. 2020) and the atom energy conservation equation (Chabert et al. 2012; Grodein et al. 2016). These equations are solved to update the ion density, electron temperature, atom density and atom temperature.

The rate coefficients of the ionization reaction and excitation reaction are two key parameters to calculate the plasma parameters in the discharge chamber, which are determined by the cross-sections of collisions and the electron energy distribution function. The cross-sections of the ionization and excitation reaction are as in Niu et al. (2024). However, the electron energy distribution function deviates from a Maxwell distribution affected by the complex cusped magnetic field and many inelastic collision processes (especially the ionization and excitation process). Therefore, the ECRIT is diagnosed by a Langmuir probe in the ECR region to obtain its electron energy distribution function.

The diagnostic result of the electron energy distribution is obtained by differentiating the transition region of I–V characteristic curves. Several typical electron energy distribution functions with parameters, shown as (2.1), are combined to approach the measurement results. Here, ${T_e}$ is equivalent electron temperature. The computations of ${c_{1i}}$ and ${c_{2i}}$ are shown as (2.2). Parameter i is the characteristic parameter of the distribution function. We can obtain a Maxwell distribution when i equals to 1 and a Druyvesteyn distribution when i equals to 2. Here, $\varGamma$ is gamma function. The parameters are fitted by the least square method and then revised to satisfy the definition, shown as (2.3). The diagnostic and fitted results are shown in figure 4(a). The rate coefficients of ionization and excitation reactions are calculated by (2.4) and shown in figure 4(b). Here, Q and $\sigma (E)$ are respectively the rate coefficient and cross-section of corresponding collision process and m is the electron mass

(2.1)\begin{gather}f(E) = \sum\limits_{i = 1}^n {{c_{1i}}T_e^{ - 1.5}\sqrt E \exp ( - {c_{2i}}{{(E/{T_e})}^i})} ,\end{gather}

(2.2)\begin{gather}\left. {\begin{array}{c@{}} {{c_{1i}} = i{{\left( {\dfrac{2}{3}} \right)}^{3/2}}\dfrac{{{{[\varGamma (5/2i)]}^{3/2}}}}{{{{[\varGamma (3/2i)]}^{5/2}}}}}\\ {{c_{2i}} = {{\left( {\dfrac{2}{3}} \right)}^{3/2}}{{\left[ {\dfrac{{\varGamma (5/2i)}}{{\varGamma (3/2i)}}} \right]}^i}} \end{array}} \right\},\end{gather}

(2.3)\begin{gather}\int_0^\infty {f(E)\,\textrm{d}E} = 1,\end{gather}

(2.4)\begin{gather}Q = \int_0^\infty {\sigma (E)} \sqrt {\frac{{2E}}{m}} f(E)\,\textrm{d}E.\end{gather}

Figure 4. Electron energy distribution function and the calculated rate coefficient. (a) Electron energy distribution function and (b) reaction rate coefficient.

The absorbed MW power is shown as (2.5) and (2.6) (Budden 1966). To reflect the saturation effect of MW power absorption, the model sets the upper limit of absorbed MW power when the plasma density reaches its critical value, shown as (2.7). The relationship between the SG sheath potential and electron temperature upstream of the SG is shown as (2.8) (Goebel & Katz 2008)

(2.5)\begin{gather}{P_{\textrm{abs}}} = {P_{\textrm{input}}}{\eta _p},\end{gather}

(2.6)\begin{gather}{\eta _p} = (1 - {\textrm{e}^{ - \mathrm{\pi }\eta }}),\end{gather}

(2.7)\begin{gather}\eta = \frac{{\omega _{\textrm{pe}}^2}}{{\omega c\alpha }};\quad {\omega _{\textrm{pe}}} = \sqrt {\frac{{{n_{\textrm{ecr}}}{e^2}}}{{\varepsilon m}}} ;\quad \alpha = \frac{1}{B}\frac{{\textrm{d}B}}{{\textrm{d}z}},\end{gather}

(2.8)\begin{gather}{n_c} = {({f_p}/9)^2},\end{gather}

(2.9)\begin{gather}{\phi _e} = {T_e}\left\{ {\ln \left( {\sqrt {\frac{M}{{2\mathrm{\pi }m}}} } \right) - \ln \left( {\frac{{{A_i}}}{{{A_e}}}} \right)} \right\}.\end{gather}

Here,${P_{\textrm{input}}}$ is the input MW power; ${\eta _p}$ is the absorbed efficiency of plasma for MW power; ${f_p}$ is the MW frequency, 4.2 GHz; c is the light speed; $\alpha$ is the normalized magnetic induction intensity gradient in the ECR region; ${\omega _{\textrm{pe}}}$ is the plasma oscillation angular frequency; ${n_{\textrm{ecr}}}$ is the average ion density in the ECR region, satisfying ${n_{\textrm{ecr}}} = {f_{\textrm{ecr}}}n$, n is the average electron density in the discharge chamber; $\varepsilon$ is the permittivity of free space; ${n_c}$ is the critical plasma density considering the saturation effect of MW power absorption by plasma. Also, M and m are respectively the Xe atom and electron masses and ${A_i}$ and ${A_e}$ are respectively the ion and electron loss areas on the SG.

We adopt the elliptic paraboloid to describe the sheath profile upstream of the SG considering the radial and angular symmetry of ECRITs, shown as (2.9). The origin of the coordinate system is located at the intersection point of the sheath and the central axis of the thruster. The z-axis is the direction of the central axis of the thruster. The x-axis and y-axis are two orthogonal axes perpendicular to the z-axis, satisfying the right-hand spiral rule. The ion beamlet shape is described by a rotating hyperboloid, shown as (2.10). Symbol ‘b ₁’ is the minimum cross-sectional area of the ion beamlet and ‘a’ is its z-direction location; symbol ‘d ₀’ describes the bending degree of the ion beamlet. The relationships between these parameters with sheath profiles are shown in Niu et al. (2024)

(2.10)\begin{gather}{x^2} + {y^2} = 2pz,\end{gather}

(2.11)\begin{gather}{x^2} + {y^2} = {d_0}{(z - a)^2} + {b_1}.\end{gather}

The reference value of the sheath thickness is the distance from the front side of the SG to the centre of the AG, shown as (2.11). Symbols ${l_s}$ and ${l_a}$ are respectively the thicknesses of the SG and the AG. Symbol ${l_{sa}}$ is the distance between the SG and the AG. The sheath thickness L_s is the distance from the origin to the centre of the AG, shown as (2.12) and related to the grid voltages, sheath potential and plasma density upstream of the SG. Symbols ${f_{\textrm{ls}}}$ and $\zeta$ are parameters of the sheath thickness. Symbol ${\lambda _d}$ is the Debye length of plasma upstream of the SG. The sheath profile is determined by the relationship between L_s and L_{s 0}, which is parallel to the SG when L_s equals to L_{s 0}

(2.12)\begin{gather}{L_{s0}} = {l_s} + {l_{sa}} + 0.5{l_a},\end{gather}

(2.13)\begin{gather}{L_s} = {f_{\textrm{ls}}}{\lambda _d}{\left( {\frac{{e({\phi_{\textrm{es}}} + {U_s} - {U_a})}}{{k{T_{\textrm{es}}}}}} \right)^\zeta }.\end{gather}

We choose a single-aperture model to study the sheath and ion beamlet characteristics with the assumption that the plasma density is uniform upstream of the SG. The simulation object is a cylinder, whose radius and height are respectively the radius of the SG aperture and the distance from the front sidewall of the SG to the back sidewall of the AG. The front sidewall refers to the surface close to the discharge chamber. The ion number conservation equation in the grids is shown as (2.13). The term on the left side of (2.13) is the temporal gradient of the average ion density. The terms on the right side of (2.13) are respectively the ion density change rates caused by ions entering the SG, extracted from the AG and impinging on the AG (sidewall and inner wall of the AG aperture). The calculation of parameters in (2.13) are shown in Niu et al. (2024)

(2.14)\begin{equation}\frac{{\textrm{d}{n_v}}}{{\textrm{d}t}} = \frac{{{n_s}{A_s}{T_{s2}}\sqrt {\frac{{e{T_{\textrm{es}}}}}{M}} }}{{N{V_{\textrm{grid}}}}} - \frac{{{f_{v1}}{n_v}{\eta _{f1}}{A_{\textrm{ah1}}}{v_{i1}}}}{{{V_{\textrm{grid}}}}} - \frac{{({f_{v2}}{\eta _{f2}}{A_{\textrm{ah2}}} + {f_{v3}}{\eta _{f3}}{A_{\textrm{ah3}}}){n_v}{v_{i2}}}}{{{V_{\textrm{grid}}}}}.\end{equation}

Here, ${n_v}$ is the average ion density for the single-aperture body; N is the SG aperture number; ${V_{\textrm{grid}}}$ is the volume of single-aperture body; ${f_{v1}}$, ${f_{v2}}$ and ${f_{v3}}$ are respectively the ratio of ion density on the AG exit, AG sidewall, AG aperture inner wall and ${n_v}$; ${\eta _{f1}}$, ${\eta _{f2}}$ and ${\eta _{f3}}$ are the area integrals of the cosine of the angle respectively between the extracted ion acceleration direction, the ion (impinging on the AG sidewall) acceleration direction, the ion (impinging on the AG aperture inner wall) acceleration direction and the z direction; ${A_{\textrm{ah1}}}$ is the extraction area of the ion beamlet; ${A_{\textrm{ah2}}}$ is the impingement area on the AG sidewall; ${A_{\textrm{ah3}}}$ is the impingement area on the AG aperture inner wall; ${v_{i1}}$ and ${v_{i2}}$ are the speeds of ions respectively extracted from the AG and impinging on the AG.

The SG current, AG current, AG current from charge exchange (CEX) ions and thrust are illustrated as (2.14)∼(2.17). The AG current includes three parts, namely the primary ions impinging on the AG sidewall and AG aperture inner wall from the discharge chamber and the CEX ions impinging on the AG sidewall from downstream of the AG (Lingwei et al. 2010). The CEX part is shown as (2.16), whose parameters can be referenced in Niu et al. (2024). The thrust includes two parts, which are respectively produced by ions and neutral gas. The atom speed is shown as (2.18). Symbols ${v_a}$ and ${T_g}$ are connected by the ideal gas state equation. The unit of ${T_g}$ is eV, which is updated by the atom energy conservation equation

(2.15)\begin{gather}{I_b} = {n_s}e{A_s}{T_{s2}}\sqrt {\frac{{k{T_{\textrm{es}}}}}{M}} ,\end{gather}

(2.16)\begin{gather}{I_a} = {f_{v2}}{n_v}e{\eta _{f2}}{N_s}{A_{\textrm{ah2}}}{v_{i2}} + {f_{v3}}{n_v}e{\eta _{f3}}{N_s}{A_{\textrm{ah3}}}{v_{i2}} + {I_{\textrm{cex}}},\end{gather}

(2.17)\begin{gather}{I_{\textrm{cex}}} = \int_0^{{z_0}} {{f_I}e{A_a}\frac{{\textrm{d}{n_{\textrm{cex}}}}}{{\textrm{d}t}}\,\textrm{d}z} ,\end{gather}

(2.18)\begin{gather}F = {\eta _{f1}}\frac{{({I_b} - {I_a})}}{e}\sqrt {\frac{{2M({\phi _{\textrm{es}}} + {U_s})}}{e}} + \frac{1}{4}{n_g}\sqrt {\frac{{8e{T_g}}}{{\mathrm{\pi }M}}} {A_a}{T_a}{\eta _c}M{v_a},\end{gather}

(2.19)\begin{gather}{v_a} = \sqrt {3e{T_g}/M} .\end{gather}

Here, ${I_b}$ is the SG current; ${n_s}$ is the plasma density upstream of the SG; ${A_s}$ is the SG area; ${T_{s2}}$ is the effective transparency of the SG; ${T_{\textrm{es}}}$ is the electron temperature upstream of the SG; ${I_a}$ is the AG current; ${\phi _{\textrm{es}}}$ is the sheath potential upstream of the SG; ${U_s}$ is the SG voltage; ${T_a}$ is the optical transparency of the AG; ${\eta _c}$ is the Clausing factor of the grid system; ${I_{\textrm{cex}}}$ is the AG current due to CEX reaction downstream of the AG; ${n_{\textrm{cex}}}$ is the density of CEX ions; ${f_I}$ is the factor describing the proportion of CEX ions hitting the AG; ${v_a}$ is the atom speed; ${T_g}$ is the neutral gas temperature.

In this paper, the integrative model is built to study the influences of operating parameters on the dynamic process of ECRITs. Therefore, the model accuracy needs to be validated. Firstly, the stable accuracies of the model with variations of MW power, flow rate and grid voltages are all validated in § 3 by experiments. Besides, the dynamic accuracy of the model is validated by adjusting MW power. This is mainly because the response time of the MW source is approximately 200 μs, much shorter than that of flow rate controller and DC sources used in the experiments. Therefore, the dynamic characteristics of ECRITs can be reflected from the grid current by adjusting MW power. In § 3, there is an analysis about the response time of the plasma parameters which refers to their stable time after operations are adjusted.

3. Thruster dynamic characteristics

3.1. Influence of MW power on thruster dynamic characteristics

When the flow rate and SG/AG voltages are set as 0.03 sccm and 600 V/−150 V, the SG current under different MW powers is shown in figure 5(a). Sccm is the unit of flow rate representing standard millilitres per minute. It can be seen that the simulation results agree well with experiments, with an error less than 5 %. Then the MW power is adjusted from 0.15 to 0.6 W by the MW source control system while keeping the flow rate and grid voltages invariant. The MW power is measured by the wave detector and its response time is approximately 200 μs. We choose this process as the model input and obtain the corresponding simulation results, which agree well with experiments, shown in figure 5(b). It shows that the response time is approximately 10 ms, which is much longer than the response process of the MW power and can reflect the dynamic characteristics of ECRITs with the adjustment of MW power.

Figure 5. Comparison of experiments and model results with the variation of MW power. (a) Stable characteristics and (b) dynamic characteristics.

Experimental results show that the characteristic time of the ECRIT dynamic process is of the order of 10 ms if the MW power can be adjusted much more rapidly. Actually, the characteristic time of the MW transfer velocity and electron acceleration by cyclotron resonance determines that the energy can be absorbed by electrons instantaneously after the adjustment of the MW power. However, rapid increase of the electron temperature breaks the ionization balance, which results in the motion of neutral gas and redistribution of atom density. Considering the volume of the discharge chamber and atom motion velocity, the recovery time of ionization balance is approximately 10 ms. More detailed analysis is listed as follows.

The flow rate and SG/AG voltages are respectively set as 0.03 sccm and 600 V/−150 V. The state parameters in the discharge chamber are shown as figure 6 when the MW power steps from 0.15 to 0.45 W with 0.15 W intervals. In the figures showing the dynamic processes of thruster parameters with the adjustment of operating parameters, we use the blue dashed line to represent the adjustment of operating parameters (MW power, flow rate and DC voltage) and use the red solid line to represent the dynamic processes of thruster parameters. It can be seen from figure 6(a) that the electron temperature firstly increases with enhancement of the MW electric field after the MW power rises. The increase of high-energy electrons can enhance the ionization process and improve electron density. However, the ionization process is a strong energy dissipation process. Therefore, the promotion of ionization intensity costs more electron energy and the electron temperature falls with the increase of electron density. The response time of this stage is tens of microseconds, which is consistent with the plasma characteristic time.

Figure 6. Simulation results of state parameters with step response of MW. (a) Electron temperature, (b) electron density, (c) atom density, (d) atom temperature, (e) source term of collisions of ions and atoms and (f) source term of collisions of ions and wall.

The response time of the next stage is approximately 10 ms due to the volume effect. The enhancement of the ionization process causes a drop of the atom density locally, especially in the ECR region, and further results in the diffusion motion of neutral gas. The motion of neutral gas causes the decline of average atom density of the discharge chamber, including the main ionization region (the cusp region). The decrease of the average atom density further results in the fall of ionization intensity, which can be reflected from the electron density at this stage. Therefore, the average electron density firstly increases and then decreases, which forms an overshoot with the step increase of MW power. The actual volume of the discharge chamber determines the time scale of neutral gas filling and it is of the order of tens of milliseconds. Therefore, the volume effect can be reflected from the atom density in figure 6(c).

It can be seen from figure 6(d) that the response process of atom temperature can also be divided into two stages. The response time of the first stage is approximately 0.2 ms, longer than the plasma response time with the increase of MW power. The main processes influencing atom temperature include the collision process between ions and neutral gas and the collision process between neutral gas and the discharge chamber wall. The source terms of these two processes are shown as figures 6(e) and 6(f). It can be seen that the characteristic time of the atom temperature variation due to collisions between ions and neutral gas is tens of microseconds, consistent with the plasma response. The characteristic time of the atom temperature variation due to collisions between neutral gas and the wall is 0.2 ms, which is the main factor determining the response time of atom temperature in the first stage.

The SG current and thrust with a step increase of MW power are shown as figure 7. Under the adjustment process, the thruster is optimal perveance and the AG current is very low. Therefore, the thrust is mainly determined by the SG current with grid voltages unchanged. The SG current firstly increases after MW power increases, which characteristic time is tens of microseconds caused by the enhancement of the ionization process. The characteristic time of the next stage is approximately 10 ms, which is the characteristic time of neutral gas diffusion motion resulting from the broken ionization balance. The overshoots of the SG current and thrust occur because of the decrease of average atom density and ionization intensity due to the neutral gas diffusion motion after the ionization balance is broken.

Figure 7. Simulation results of thruster performance with step response of MW. (a) Screen grid current and (b) thrust.

Preliminary experiments have validated that thrust noise can be obviously reduced by feedback control when the MW power and ion beam current are respectively chosen as the controlled variable and feedback parameter (Zhou et al. 2022). Therefore, four first-order inertia response processes of MW power, shown as (3.1)∼(3.4), are used to further study their influences on the SG current. Their response times are respectively 1 ms, 5 ms, 10 ms and 15 ms, which are chosen considering that the characteristic time of the SG current with the step variation of MW power is approximately 10 ms. To further study the influence of MW power overshoot characteristics on the SG current, the second-order under-damped processes with the same response time are used as the model input, shown as (3.5)∼(3.8)

(3.1)\begin{gather}{G_p} = \frac{{5000}}{{s + 5000}},\end{gather}

(3.2)\begin{gather}{G_p} = \frac{{1000}}{{s + 1000}},\end{gather}

(3.3)\begin{gather}{G_p} = \frac{{500}}{{s + 500}},\end{gather}

(3.4)\begin{gather}{G_p} = \frac{{333.33}}{{s + 333.33}},\end{gather}

(3.5)\begin{gather}{G_p} = \frac{{3.6 \times {{10}^7}}}{{{s^2} + 6000s + 3.6 \times {{10}^7}}},\end{gather}

(3.6)\begin{gather}{G_p} = \frac{{4 \times {{10}^6}}}{{{s^2} + 2000s + 4 \times {{10}^6}}},\end{gather}

(3.7)\begin{gather}{G_p} = \frac{{{{10}^6}}}{{{s^2} + 1000s + {{10}^6}}},\end{gather}

(3.8)\begin{gather}{G_p} = \frac{{4.444 \times {{10}^5}}}{{{s^2} + 666.7s + 4.444 \times {{10}^5}}}.\end{gather}

The SG current under various MW power response processes is shown in figure 8. It can be seen from figure 8(a) that the SG current shows overshoot characteristics although there is no overshoot for MW power under the first-order inertia process. The overshoot of the SG current declines as the MW power response process is prolonged. It can be concluded from the above analysis that the average atom density decreases due to the neutral gas diffusion after the ionization balance is broken by the rise of MW power. The characteristic time is approximately 10 ms. When the response time of MW power is less than 10 ms, the MW power changes little during the time interval from the response time of the MW to that of the thruster (10 ms). However, the neutral gas density and ionization intensity still fall in this time interval, which results in the decline of the SG current. It is the main reason why there is an overshoot for the SG current and its response time is almost unchanged under different MW response processes when the response time of MW power is less than 10 ms. Besides, the overshoot of the SG current declines as the MW power response process is prolonged because the continuous increase of MW power can weaken the decrease degree of the SG current caused by the fall of the average atom density. When the response time of MW power exceeds the characteristic time of the volume effect (10 ms), the increase of MW power can make the SG current increase continuously and therefore the overshoot phenomenon of the SG current gradually disappears.

Figure 8. Simulation results of screen current under different responses of MW. (a) First-order inertia process and (b) second-order under-damped process.

It can be seen from figure 8(b) that there are always overshoots for SG current under the second-order under-damped MW power response processes because the output MW power itself has overshoots. Under this condition, the response process of SG current shows similar characteristics to that under the first-order inertia process. When the response time of MW power is shorter than the characteristic time of the volume effect, there is a stage where the SG current still falls after the MW power becomes stable. This phenomenon gradually disappears when the response process of MW power is prolonged. Besides, it can be seen from figure 8(b) that the overshoots of the SG current decrease under slower MW power response processes, although the overshoots of MW powers are the same under these second-order under-damped response processes. This phenomenon reflects that the overshoots of the SG current depend on not only the overshoots of MW power, but also the response times of MW power, which are closely related to the characteristics of ECRITs.

3.2. Influence of flow rate on thruster dynamic characteristics

When the MW power and SG/AG voltage are set as 0.3 W and 500 V/−100 V, the SG currents under different flow rates is shown as figure 9, which shows that the simulation error is less than 5 %.

Figure 9. Comparison of experiments and model results with the variation of flow rate.

The state parameters in the discharge chamber are shown in figure 10 when the flow rate increases from 0.03 to 0.1 sccm in two steps. It can be seen from figure 10 that the average atom density increases with the step rise of the flow rate, which can prompt the ionization intensity. Therefore, electron density rises with the enhancement of ionization intensity, which consumes more electron energy and causes the drop of electron temperature. Besides, the electron density increment caused by the increase of flow rate from 0.03 to 0.06 sccm is more obvious than that from 0.06 to 0.1 sccm. Correspondingly, the atom density increment caused by the increase of flow rate from 0.03 to 0.06 sccm is less obvious than that from 0.06 to 0.1 sccm. It is mainly because the average electron temperature is higher under the first adjustment process, which causes that the enhancement of ionization intensity under the first adjustment process is more obvious than that under the second adjustment process. The response time of state parameters is 10 ms, consistent with the volume effect.

Figure 10. Simulation results of state parameters with step response of flow rate. (a) Atom density, (b) electron density, (c) electron temperature and (d) atom temperature.

It can be seen from figure 10(d) that the atom temperature firstly decreases and then rises after the step increase of the flow rate. The response time of the decrease process is 0.1 ms. It is mainly because the temperature of neutral gas flowing into the discharge chamber is below that previously in the discharge chamber. Accompanied by this process, the ion density increases with the enhancement of the ionization process. The collisions between ions and neutral gas heat the neutral gas and result in the increase of atom temperature. The response time of this process is approximately 10 ms. When the flow rate increases from 0.03 to 0.06 sccm, the ion density increment caused by the enhancement of ionization intensity is large, which leads to a strong heat effect of ions on atoms. Therefore, the atom temperature stays at a higher level after the flow rate increases. However, when the flow rate increases from 0.06 to 0.1 sccm, the heat effect of ions on atoms is weak, which causes the atom temperature to decrease due to inflow of low-temperature neutral gas. The response processes of the SG current and thrust with the step increase of flow rate are shown in figure 11.

Figure 11. Simulation results of thruster performance with step response of flow rate. (a) Screen grid current, (b) thrust and (c) sheath thickness.

It can be seen that the response characteristics of the thrust are similar to those of the SG current. It is mainly because the perveance condition of the thruster is under optimal perveance and the AG current is negligible compared with the SG current. Besides, the grid voltage is invariant during the adjustment process. Therefore, the thrust is mainly determined by the SG current. It can be seen that the response time of the SG current is approximately 10 ms, which is consistent with the characteristic time of the volume effect. Besides, it can be seen that the SG current decreases when the flow rate increases from 0.06 to 0.1 sccm. The SG current depends on the ion density upstream of the SG, ion Bohm velocity and ion extraction area. Although the increase of flow rate directly results in the rise of ion density upstream of the SG, the increment is small because the ionization intensity is limited by the low electron temperature. Meanwhile, the increase of flow rate results in the drop of electron temperature, which reduces the ion Bohm velocity upstream of the SG. Besides, it can be seen from figure 11(c) that the sheath thickness decreases with the increase of ion density upstream of the SG, which further results in the decline of the ion extraction area. Therefore, the SG current decreases when the flow rate rises from 0.06 to 0.1 sccm. When the flow rate increases from 0.03 to 0.06 sccm, the ion density increases at large scale. Therefore, the SG current still increases although the sheath thickness and Bohm velocity decrease under this adjustment process.

3.3. Influence of grid voltage on thruster dynamic characteristics

When the MW power and flow rate are set as 0.6 W and 0.06 sccm, the SG current under different SG voltages is shown as figure 12. The corresponding AG voltages are respectively −50 V, −50 V, −100 V, −150 V, −150 V and −200 V with the increase of SG voltages. It can be seen that the simulation results agree well with experiments, with an error of less than 5 %.

Figure 12. Comparison of experiments and model results with the variation of grid voltage.

The flow rate and MW power are respectively set as 0.06 sccm and 0.6 W. The state parameters in the discharge chamber are shown in figure 13 when the SG voltage steps from 200 to 500 V with 150 V intervals. The AG voltage remains unchanged during the first adjustment process and steps from −50 to −100 V during the second adjustment process. It can be seen that the electron density firstly decreases with the step increase of SG voltage because of the stronger extraction effect on the ions upstream of the SG by the grid system. The drop of electron density reduces the ionization intensity and less electron energy is consumed by the ionization process, which effects lead to the rise of the electron temperature. Ions in the discharge chamber move with the Bohm velocity, which is estimated to be 3 × 10³ m s⁻¹ if we assume the average electron temperature is 8 eV. The channel length of the thruster is 15 mm, which means that the characteristic time of this process is approximately 5 μs. The main drop of electron density occurs on a time scale that is of the order of this time interval. However, the decline of ionization intensity results in the increase of atom density. This response process is determined by the diffusion motion of neutral gas. Simulation results show that this process lasts approximately 0.5 ms. It can be seen that the response process of atom density due to step increase of the SG voltage is much faster than that due to the step increase of the MW power (figure 6c). It is mainly because the influence of grid voltages on state parameters in the discharge chamber is much weaker than that of MW power, which can be validated by the comparison between figures 6 and 13. Increase of atom density enhances the ionization intensity and results in the overshoot of electron density and electron temperature. Meanwhile, the response time of electron density and electron temperature is also prolonged to 0.5 ms because of the broken balance of ionization process and diffusion motion of neutral gas.

Figure 13. Simulation results of state parameters with step response of grid voltage. (a) Electron density, (b) electron temperature, (c) atom density, (d) atom temperature, (e) sheath thickness and (f) extraction area.

It can be seen from figure 13(d) that the variation characteristic of atom temperature is similar to that of ion density. It is mainly because the heat effect of ions on atoms by collisions dominates when the flow rate is stable. Besides, it can be seen the variation range of state parameters during the second adjustment process is larger than that during the first adjustment process, although the variations of SG voltages are the same. Analysis shows that it is mainly caused by the variation range of the extraction area. It can be seen from figure 13(e) that variations of sheath thickness are similar during the two adjustments. However, the sheath is thicker during the second adjustment process. Therefore, a similar variation of sheath thickness can result in a larger variation of the extraction area during the second adjustment, which can be validated from figure 13(f). Besides, it can be seen from figure 13(e) that the sheath thickness firstly increases and then falls, which forms an overshoot. The sheath thickness firstly increases mainly because of the stronger extraction capability and decline of ion density upstream of the SG with the step increase of grid voltages. Then the sheath thickness decreases mainly because the ion density upstream of the SG increases while the grid voltage remains invariant in the latter part of the response process.

The response processes of thruster performances are shown as figure 14. It can be seen that the response time of the SG current for the first and second adjustment processes are respectively 5 μs and 0.5 ms. The response time of the SG current during the second adjustment agrees well with that of state parameters, which is mainly influenced by the diffusion motion of neutral gas. However, the response time of the SG current during the first adjustment is consistent with the time interval when ions move through the discharge chamber with the Bohm velocity. The motion of neutral gas seems to have negligible influence on the response process of the SG current, which is caused by two reasons. On the one hand, the increments of state parameters are smaller under the first adjustment process compared with those under the second adjustment process. On the other hand, it can be seen from figure 13 that the variation direction of the ion density upstream of the SG is opposite to that of the ion Bohm velocity and ion extraction area. The SG current is determined by the product of these three parts. Therefore, variations of state parameters have negligible influence on the response of SG current when their ranges are inconspicuous.

Figure 14. Simulation results of thruster performance with step response of grid voltage. (a) Screen grid current, (b) AG current, (c) ion beam current and (d) thrust.

It can be seen from figure 14(b) that the AG current decreases to a very low level when the SG voltage increases from 200 to 500 V. When the SG and AG voltages are respectively 200 V and −50 V, the SG sheath is very thin and under the over-perveance condition. Under this condition, the ions ejected from the SG sheath tend to move into the sidewall of the AG and cause the AG current to be high. With an increase of grid voltages, the thruster gradually changes from over-perveance to optimal perveance and most ions ejected from the SG sheath can escape from the AG holes, which directly leads to the decrease of the AG current. The response time of the AG current for the first adjustment is approximately 5 μs, which is consistent with that of the SG current. There is an overshoot for the AG current, which is much higher than that of the SG current. The ion density in the grid system firstly decreases obviously with the increase of grid voltages. Then ions move from the discharge chamber to the grids, which contributes to the rise of ion density in the grid system and results in an overshoot of the AG current. There are no overshoots of the AG current for the second adjustment. This is mainly because most ions can escape from the AG holes when the SG sheath is under the optimal-perveance condition. The AG current mainly comes from the slow ions produced by CEX collisions downstream of the AG. Therefore, the ion density variation has negligible influence on the AG current for the second adjustment. It can be seen that the ion beam current and thrust have similar response processes with the step increase of grid voltages. The overshoot is obvious and its response time is approximately 5 μs for the first response process. The overshoot is inconspicuous and its response time is approximately 0.5 ms for the second response process. It is mainly determined by the response processes of the SG current and AG current.

4. Conclusion

In this paper, an integrative model of a minimized ECRIT is built to study the influence of operating parameters on the dynamic process of a thruster. The stable and dynamic model accuracy is validated by experiments. The main conclusions are listed as follows. These conclusions can not only facilitate designing the feedback controller of micro-propulsion systems applied in high-precision space missions, but also provide guidance for ion sources to generate a highly stable or rapid-response ion beam.

1. (i) The electron temperature increases sharply with the step rise of MW power and then decreases with the increase of electron density because more electron energy is consumed with the enhancement of ionization intensity. The response time of this process is tens of microseconds, consistent with the plasma response time. Then the broken ionization balance results in diffusion motion of neutral gas, which further causes a decrease of electron density. Therefore, there are overshoots of the electron density and SG current. The response time of this process is 10 ms, consistent with the volume effect. Further simulation results show the overshoot of the SG current decreases under a slower MW response process, although the response time of the SG current remains 10 ms when the response time of MW power is shorter than 10 ms.

2. (ii) The atom density increases gradually with the step increase of flow rate, which strengthens the ionization intensity and leads to an increase of electron density and drop of electron temperature. During this process, the atom temperature firstly decreases due to inflow of low-temperature neutral gas and then increases because of the heat effect of ion–atom collisions. The response time of state parameters is the same as that of the SG current, namely approximately 10 ms, which is consistent with the volume effect. Besides, simulation results show that the SG current decreases with the increase of flow rate for high flow rate operations. This is mainly because the increment of ion density upstream of the SG by the increase of flow rate is limited by the low electron temperature when the flow rate is high. Under this condition, the decreases of ion Bohm velocity and sheath extraction area cause the decline of the SG current.

3. (iii) The step increase of SG voltage directly results in the decrease of electron density and ionization intensity in the discharge chamber due to the stronger extraction capability of grids. This process breaks the previous ionization balance and causes the increase of atom density and an overshoot of electron density due to the diffusion motion of neutral gas. The response time of state parameters is approximately 0.5 ms, which is slower than the response time caused by a step rise of MW power. This is mainly because the influence of grid voltages on plasma parameters of the discharge chamber is very limited. Besides, simulation results show that the response time of the SG current is approximately 5 μs with the step increase of grid voltage when the variations of extraction capability are limited, which is consistent with plasma response time. With the increase of extraction capability variation range, the response time of the SG current is prolonged to 0.5 ms because larger variations of grid voltages have stronger influences on the plasma parameters of the discharge chamber.