Importance of stepwise ionization from the metastable state in electron cyclotron resonance ion thrusters

In electron cyclotron resonance (ECR) thrusters, the plasma mode transition is a critical phenomenon because it determines the maximum thrust performance. In ECR ion thrusters, ionization generally occurs in the magnetic confinement region, where electrons are continuously heated by ECR and confined by magnetic mirrors. However, as the flow rate increases, ionization is also observed outside the magnetic confinement region, and this induces the plasma mode transition. In our previous work, two-photon absorption laser-induced fluorescence (TALIF) analysis revealed that the stepwise ionization from the metastable state plays an important role in the ionization process. However, the distribution of the stepwise ionization has not yet been revealed because of the long lifetime of the metastable state. In this study, this distribution was investigated using one experimental and two numerical approaches. First, TALIF was applied to two types of gas injection with clear differences in thrust performance and ground-state neutral density distribution. In the first simulation, the metastable state particle simulation was used to estimate the excitation rate distribution. In the second study, simulations of the electric field of microwaves were used to estimate the contribution of the stepwise ionization to the plasma density. The experimental and numerical results revealed that the stepwise ionization spreads outside the magnetic confinement region because of the diffusion of metastable particles, and this spread induces the plasma mode transition, explaining the difference between the two types of gas injection.

However, measurement of the local excitation and ionization rates is still challenging. Generally, the ionization and excitation rates can be written as k i, j = σ i, j v e n e n j , where k i, j is i-th production rate from the lower state j (e.g., the ground state), σ i, j v e is the rate coefficient from state i to state j, n e is electron density, and n j is the number density at the lower state j. To evaluate σ i, j v e , the cross-section σ i, j and electron velocity distribution function are necessary. Of course, this function can be theoretically evaluated using the Langmuir probe and LTS; however, it is very difficult to measure non-Maxwellian electrons because of the low signal ratio. Additionally, in low-temperature plasmas, not only direct ionization from the ground state, but also ionization from the lower excited state, i.e., stepwise ionization, can be important because relatively low-energy electrons can ionize the excited neutral particles. For example, plasma simulation models including these excited states have been proposed by Boeuf [21], Hagelaar [22], and Hara [23]. Chaplin suggests that the stepwise ionization from the metastable state is very important in Hall effect thrusters by collisional radiative model (CRM) [24].
Hence, an experimental evaluation of the ionization process that includes direct and stepwise ionization is strongly desirable. In our group, understanding the ionization process of the ECR ion thrusters used in the Hayabusa [25] and Hayabusa2 [26] missions is one of the most important issues. In the thrusters, the plasma mode transition occurs as the input flow rate increases, which restricts the maximum thrust performance. Generally, in ECR ion thrusters, ionization occurs in the magnetic confinement region where electrons are continuously heated by the magnetic mirror confinement. However, we previously reported that ionization is also observed outside the magnetic confinement region as the flow rate increases, which induces plasma mode transition [27]. Therefore, to investigate the ionization process, we proposed a method for investigating excitation process by the ground-state neutral density and spontaneous emission intensity at that point [28]. Then, an analysis of the excitation process was used to estimate the ionization process. The measurement revealed that stepwise ionization from metastable neutral particles 1s 5 (Paschen notation) plays an important role in the ionization process [28].
However, the distribution of the stepwise ionization is not yet clear because the metastable state particles have a long natural lifetime [29]. In this study, to investigate this distribution, one experimental and two numerical approaches were used. In the experiment, the ground-state neutral density and spontaneous emission intensity were measured using two-photon absorption laser-induced fluorescence (TALIF) for two types of gas injections that are substantially different not only in thrust performance [30], but also in metastable neutral density [31], ground-state neutral density [32], and spontaneous emission intensity [33]. These differences can be useful for revealing the distribution of stepwise ionization. In one simulation, a metastable state particle simulation considering de-excitation due to collision and wall diffusion was performed to specify the excitation rate distribution in the metastable state. Then, the simulation results were validated using the density distribution obtained by laser absorption spectroscopy (LAS) [30,31]. Second, to estimate the plasma density distribution, the electric field of the microwaves was simulated because it strongly depends on the plasma density (e.g., cut-off phenomenon). Here, the plasma density profile was calculated using the plasma diffusion equation [34], thus revealing the contribution of the stepwise ionization. The simulated electric field was also validated using experimental data measured by an electro-optic (EO) probe [35]. These numerical approaches are useful for estimating unobserved parameters (e.g., the excitation rate distribution and stepwise ionization distribution) and the correlation between each pair of parameters.
The rest of this paper is organized as follows. The configuration and an explanation of the plasma mode transition are presented in "ECR ion thruster" section. The TALIF measurements of the ground-state neutral density and spontaneous emission are given in "Ground-state neutral density and spontaneous emission intensity measurements" section. The numerical models are described in "Numerical simulation" section. The experimental results, simulation results, and discussion are presented in "Results" section. Finally, the conclusions are summarized in "Conclusion" section.

ECR ion thruster
Configuration A schematic of the ECR ion thruster is shown in Fig. 1. Microwaves are transmitted into the discharge chamber from the antenna at the bottom of the waveguide. In the discharge chamber, two magnet rings generate a mirror magnetic field. Larmor radius of electrons in the mirror magnetic field is the order of 0.001-0.01 cm. Hence, electrons are strongly magnetized and confined by the mirror magnetic field. At that time, Fig. 1 Schematic of a gridded microwave discharge ion thruster. WG and DC injections were evaluated. The ECR region is where the electron cyclotron frequency and microwave frequency match electrons are continuously heated by ECR in the mirror magnetic field, which is called the magnetic confinement region. Then, the electrons heated by ECR ionize neutral particles, and then plasma is generated. Larmor radius of ions is the order of 0.1-1.0 cm, which indicates that ions are weakly magnetized. As a results, ions can be transported to the grid and ions are accelerated by a high voltage between the screen grid and accelerator grid.
In this study, two types of gas injection were used: waveguide (WG) injection and discharge chamber (DC) injection. In WG injection, the propellant is injected from the bottom of the waveguide, whereas in DC injection, the propellant is injected between the two magnet rings. The microwave power was fixed at 34 W, and the flow rate was varied from 1 to 4 standard cubic centimeters per minute (sccm). The screen and accelerator voltages were fixed at 1500 and − 350 V, respectively. The background pressure was 2.0 × 10 − 3 Pa and contained 3.0 sccm of xenon. The vacuum facilities are described in detail elsewhere [27].

Plasma mode-transition
There are two main types of magnetic field in ECR plasma sources: one mirror magnetic field (our study [28] and magnetic nozzle [36]), and multipole type in azimuthal direction [37,38]. Plasma mode transition has been observed in both types of magnetic field [36][37][38]. The mode transition in an ECR ion thruster is shown in Fig. 2 [27]. Here, when the flow rate is optimal, the thrust (ion beam current) is maximized. At a low flow rate (below the optimal flow rate), the plasma is concentrated in the magnetic confinement region. In contrast, as the flow rate increases, plasma is also generated outside the magnetic confinement region. At the optimal flow rate, red luminescence is generated inside the boundary of the waveguide. At rates higher than the optimal flow rate, the red luminescence further increases. This transition indicates that the extraction rate distribution changes drastically as the flow rate increases. Additionally, because the ionization process strongly correlates with the excitation process, the ionization distribution should also change.

Ground-state neutral density and spontaneous emission intensity measurements
Principle Using the TALIF experimental setup, the ground-state neutral density and spontaneous emission intensity were measured at the same point, as shown in Fig. 3. Our idea is based on the fact that the wavelengths of fluorescence and spontaneous emission by electron-impact excitation are the same. For instance, in normal laser-induced fluorescence on metastable ions, the wavelength of the fluorescence is different from that of spontaneous emission [13,14]. Thus, in TALIF measurement, the main background noise is the spontaneous emission at the same wavelength, and hence both parameters can be measured using the same setup. Generally, the spontaneous emission intensity I j can be expressed as follows [39].
Here, I 0 is a constant and n e is the electron density. The subscripts i, j, and g are the upper, lower, and ground states, respectively. This equation states that the spontaneous emission is the sum of the direct excitation from ground state g and the stepwise excitation from lower excited state i. Because the spontaneous emission intensity is proportional to n g , the stepwise excitation can be evaluated by the correlation between n g and spontaneous emission I j . Noted that Eq. (1) assumes optically thin, that is, the radiation trapping is neglected [24]. Additionally, because the ionization process is strongly related to the excitation process, the ionization process is estimated using an analysis of the excitation process.

Determination of excitation process
Generally, there are many selections of two-excitation processes as summarized in [40]. We propose that the excitation process should be selected following three perspectives.
(1). Trade-off between fluorescence intensity and spontaneous emission intensity (2). Lower state is resonant state (3). Difference between energy of excitation and singly-charged ionization Because the spontaneous emission at the same wavelength of fluorescence is the inevitable noise, (I) is required for measuring the ground-state neutral density with high signal to noise (S/N) ratio as pointed in C. Eichhorn [41]. (II) is required for minimizing the uncertainly of optically depth. Generally, if the optical depth is not sufficiently thin, the spontaneous emission intensity decreases because the excited neutral density at lower state reabsorbs in the spontaneous emission [24]. Hence, the lower state should be resonant state because the density is generally lower than that of metastable state. Last, (III) is required for minimizing the difference of energy between the excitation and singly-charged ionization. In this thruster, the singly charged ions account for more than 90% of all ions [9], and hence we focus on the ionization process of singly-charged ions. In this experiment, the ionization process is investigated by using the analysis of the excitation process; thus, the discrepancy of electron population that can react excitation and ionization should be small.
Based on the above criteria, we recommend the excitation process of 1s 0 to 2p 3 as shown in Fig. 3. In this excitation process, the fluorescence at 834.68 nm emits from the transition 2p 3 to 1s 2 . From the previous result of spectra [32], the wavelength of 834.68 nm is the best choice in view of point (I). Then, because lower excited state 1s 2 is resonant state, (II) is satisfied. Actually, we confirmed that the assumption of optically thin establishes within 1% of the spontaneous emission intensity under the typical electron density of 10 16 − 10 17 m −3 [23,24,42]. In point (III), the excitation energy is 11.1 eV and the difference in energy between excitation and singly-charged ionization is only 1 eV [43], which indicates that if the stepwise excitation is non-negligible, the stepwise ionization of singly-charged ions is also non-negligible [28].

Experimental setup
Details of the experimental setup can be found in [28,32]. Hence, only a short explanation is provided here. A Nd:YAG laser (NL310, EKSPLA) at a wavelength of 532 nm that emitted 300 mJ was utilized as the pump laser. Then, the 532 nm laser was first converted to 672.87 nm by a dye (LDS-698) and its third harmonic (672.87/3 = 224.29 nm) was extracted using a dye laser (LiopStar-E-N-D2400, LIOP-TEC). Two quartz windows with diameters of 4 mm were opened for laser injection and detection while the thruster performance was maintained. The laser was focused on the center of the waveguide. To detect the fluorescence and spontaneous emission at 834.68 nm, a band-pass filter and a photomultiplier tube were used. The full width at half maximum of the band-pass filter (Sigma Koki Co., Ltd.) is less than 2 nm to eliminate the emissions at 828.0 nm and 840.9 nm [28]. To reduce the noise due to fluctuations in spontaneous emission, 3000 fluorescence signals were averaged. When measuring the spontaneous emission intensity, the laser was stopped and the signal was recorded. The calibration to obtain the ground-state neutral density is given in [28,44]. The ground-state neutral density and spontaneous emission were measured at the three distances from the screen grid: − 7, − 11.5, and − 16 cm.

Metastable state particle simulation
Metastable neutral particles 1s 5 have a long natural lifetime (40 s) [29], but their lifetime can be reduced by collision and wall diffusion significantly. The de-excitation frequency ν de is the sum of the collision with other particles ν coll and wall diffusion ν wall , that is, ν de = ν coll + ν wall . Considering thruster length L ≈ 10 cm and thermal velocity of metastable neutral particles v th ≈ 10 2 m/s, ν wall can be estimated as ν wall ≈ v th /L ≈ 10 3 Hz. Two types of collisions are possible: collisions with ground-state neutral particles ν g and collisions with electrons ν e . However, ν g can be negligible compared to ν e [45] thus, ν coll ≈ ν e . ν e can be expressed as ν e = σ total v e n e , where σ total is the total cross-section of de-excitation, v e is the electron particle velocity, and n e is the electron density. The total cross-section can be expressed as follows [39,46].
The sum over i occurs for the transition from the upper states (2p i , 3d i , 2s i , 3p i ) to 1s 4 , which is optically coupled to the ground state, and intra-transition from 1s 5 and 1s i (i = 1, 2, 3, 4). The sum over k occurs in stepwise ionization including singly charged and multiply charged ionization. From the probe measurement and full kinetic simulation [42], the electron temperature in the discharge chamber and waveguide was varied from 1 to 20 eV. From the simulated cross-section using the relativistic distorted wave approach [20,47], we have σ 1, k ≈ 10 −20 − 10 −19 m 2 in the range of T e = 1-20 eV. Then, from experimental measurements, the stepwise ionizations are σ 2, k ≈ 10 −20 − 10 −19 m 2 in the range of T e = 1-20 eV [48]. Additional factors must also be considered. First, there is the possibility that the de-excited particles will reproduce, but this reproduction is not dominant [46]. Other reactions such as Penning ionization and three-body collision are negligible because of the low pressure [21]. As a result, σ total ≈ 10 −19 − 10 −18 m 2 is valid in the range of T e = 1-20 eV, which is consistent with the results of Karabadzhak [39].
From the full-kinetic simulation [42], the electron density is n e ≈ 10 16 − 10 17 m −3 inside the discharge chamber. Thus, considering the electron velocity v e ≈ 10 6 m/s, ν e ≈ 10 3 − 10 5 Hz was used. However, inside the waveguide, the electron density can be less than that in the discharge chamber, especially near the bottom of the waveguide. If electron density n e ≤ 10 16 m −3 , ν e can be smaller than ν wall , which indicates that wall diffusion is dominant. Some studies have reported that the wall diffusion model is ν wall = (2.405 2 D x /r 2 ) −1 + (ν 0 /r) −1 [19,34]. However, in the thruster, it is difficult to model the diffusion between the waveguide and discharge chamber. Thus, in this study, wall diffusion was directly evaluated by particle simulation. Specifically, the deexcitation due to wall diffusion was simulated by deleting the particles when they reached the wall. The Monte Carlo approach was used to simulate the collision [49]. Simulated particles were deleted using the following probability of collision:P c = 1 − exp(−ν coll Δt), where Δt is the simulation time step. To keep the calculations within a numerical error of 1%, ν coll Δt = 1/50 was selected [49].
Finally, the excitation rate distribution at the metastable state K ex was modeled. According to the plasma mode transition, the excitation rate distribution in the metastable state can be expressed using driving parameters η w and L w . First, η w is defined as a linear interpolation of the excitation production rate at waveguide K ex, wave and magnetic confinement region K ex, mag as K ex, total = η w K ex, wave + (1 − η w )K ex, mag . Here, K ex, total is the total excitation production rate. Second, to express the axial distribution, the generation length L w inside the waveguide was introduced, defined as shown in Fig. 4. For simplicity, the metastable neutral particles were generated randomly in each magnetic confinement region and the waveguide. To specify the magnetic field line, the vector potential was calculated using finite element method magnetics [50].

Effect of stepwise ionization on microwave propagation
To estimate the contribution of the stepwise ionization to plasma density, the electric field of the microwaves was simulated. For simplicity, a quasi-1D simulation model was used. By specifying the electromagnetic field and current density as E y ðz; tÞ ≈ e E y ðzÞe iω m t and J y ðz; tÞ ≈ e J y ðzÞe iω m t , the following equation can be obtained [51].
Here, the electric field is assumed to be only an electromagnetic field, i.e., the electrostatic field is neglected: ∇ ! Á E ! ¼ 0. Moreover, ω m is the frequency of the microwaves and ∇ 2 t is the transverse part of the Laplacian. To simulate the electric field in the central axis, we employed quasi-1D approximation [52]. The detailed derivation and the justification were written in Appendix. Using this approximation, the transverse part can be approximated as ∇ 2 t e E y ðzÞ ≈ −ðρ 1;1 =RÞ 2 e E y ðzÞ, where ρ 1, 1 ≈ 1.84 is the first posi- Fig. 4 Model of metastable generation using driving parameters η w and L w . Parameter η w is the ratio of the ionization rates in the waveguide and the discharge chamber, and L w is the distance of the ionization area from the screen grid tive root of ∂J B1 /∂r = 0 (J B1 is first Bessel function) and R is the radius and ρ 1, 1 corresponds to the eigenvalue of transverse electric (TE) 11 mode [35]. Additionally, the current density was estimated by cold-electron approximation, i.e., e J y ðzÞ ≈ σ p ðzÞ e E y ðzÞ, where σ p (z) is the complex plasma conductivity, expressed as follows where, ν m is the momentum frequency and ω ce is electron cyclotron frequency. Noted that the conductivity of ions is neglected because ions plasma frequency are much larger than microwave one. Additionally, because the electron cyclotron frequency is one order smaller than the microwave frequency and the momentum collision frequency is three order small than the microwave frequency in the central axis [28], the plasma conductivity is approximated in the non-magnetized case as follows.
Using these approximations, one Helmholtz equation is obtained.
where, ω pe is electron plasma frequency. Finally, to determine the plasma density, the 1D plasma diffusion equation is solved as follows [28].
where ν iz, g and ν iz, m are the collision frequencies of the ground-state and metastable neutral particles, respectively. In this study, to investigate the contribution of stepwise ionization to the plasma density, the ratio of ground-state neutral density to metastable neutral density n m /n g was varied.
To solve Eqs. (6) and (7) numerically, a tridiagonal matrix algorithm was utilized. To solve Eq. (6), the perfect electrical conductivity condition e E y ¼ 0 was employed as the boundary condition at the grid and the bottom of the waveguide. To solve Eq. (7), the Bohm condition −D∂n p /∂x = Γ b was employed at the grid. Here, Γ b is the Bohm flux, where Γ b = n p v b exp(−1/2)η sc , v b is the Bohm velocity, and η sc ≈ 0.8 is the effective transparency of the screen grid. At the bottom of the waveguide, an artificial boundary condition of n p = 0 from a certain length z < L p inside the waveguide was employed. The boundary condition was derived from the radial transport of electrons, i.e., the loss of the side wall at the waveguide. The magnetic field is not parallel to any axial direction except for the central axis (Fig. 1). Because the plasma is strongly restricted along the magnetic field, the plasma is lost at the wall. A length L p = 9 cm was empirically determined by agreement in the experimental results of the microwaves.

Ground-state neutral density versus spontaneous emission intensity
The ion beam currents with respect to the propellant flow rate are shown in Fig. 5. The red points indicate the results at the optimal flow rate. In the WG injection, the maximum ion beam current is 138-139 mA at the optimal flow rate of 1.8 sccm. In the DC injection, the maximum ion beam current is 168-175 mA at the optimal flow rates of 2.9 sccm or 3.0 sccm. The spontaneous emission intensity with respect to the groundstate neutral density is shown in Fig. 6. In all cases, the spontaneous emission intensity rapidly increases with respect to the ground-state neutral density around the optimal flow rates (red points). From Eq. (1), the sharp increase can be caused by three components: stepwise excitation, an increase in electron density, and an increase in electron temperature. In a previous study, by eliminating the ECR layer inside the waveguide, the magnitude of the sharp increase was decreased by approximately 50%, but it was still present, and the plasma mode transition was not eliminated [28]. Therefore, the sharp increase cannot be explained only by the increase in electron temperature. Considering the rate coefficient and the density ratio of metastable neutrals in the range of T e = 3 to 5 eV measured by Langmuir probe [53], the spontaneous emission from metastable neutral particles 1s 5 accounts for 20%-250% of the spontaneous emission from the ground state [28]. In addition, the rate coefficient suggests that stepwise ionization occurs in the same order as that of stepwise excitation, which increases electron density [28]. Comparing WG injection with DC injection, the ground-state neutral density of WG injection is higher than that of DC injection, which is consistent with the direct simulation Monte Carlo (DSMC) results [54]. In WG injection, the neutral density increased towards the bottom of the waveguide. By contrast, in DC injection, the density at 3.5-4.0 sccm decreases towards the bottom of the waveguide. Thus, the absolute value of the density distribution of the ground state differs in WG and DC injection, which leads to a difference in the optimal flow rate. Focusing on the axial dependence of the spontaneous emission intensity, the maximum value of spontaneous emission decreases by only approximately 48% from z = − 7 cm to z = − 16 cm in WG injection, whereas it decreases by approximately 25% in DC injection. In "Spatial structure of plasma modetransition" section, we discuss how the difference in ground-state neutral density and spontaneous emission generates the difference in thrust performance between WG and DC injection based on the simulation results.

Estimation of excitation rate distribution at metastable state
As explained in "Metastable state particle simulation" section, the de-excitation frequency ν de can be varied in the range ν wall ≈ 10 3 Hz ≤ ν de ≤ ν wall + ν coll ≈ 10 5 Hz. Thus, to evaluate the sensitivity of ν de with respect to the metastable neutral density distribution, four cases, ν coll = 0, 10 3 , 10 4 , 10 5 Hz, were simulated. As an example, Fig. 7 shows the normalized 2D metastable neutral density distribution for two cases: (a) η w = 0% and (b) η w = 100%, where L w = 10 cm. The simulation results indicate that in both cases, ν coll is very sensitive to the density distribution when ν coll ≥ 10 4 Hz because ν wall is on the order of 10 3 Hz.
The simulated metastable neutral density was compared with the experimental one in the central axis obtained by LAS [31,33]. In numerical simulations, it is difficult to accurately match the absolute value of excitation rate with the experiment. Thus, in this study, each absolute value of metastable neutral density is matched with the experimental result. Five cases for L w , that is, L w = 6, 7, 8, 9, 10 cm, and four cases for η w , that is, η w = 0, 0.25, 0.5, 0.75, 1.0, were simulated. Finally, to obtain the variation of de- excitation frequency due to collision, ν coll = 0, 10 3 , 10 4 , 10 5 Hz were simulated for each case. In total, 5 × 5 × 4 = 100 cases were simulated, and the simulation results that minimized the discrepancies in the experiment al results were then selected. Figure 8 shows the axial distribution of metastable state 1s 5 density. Each simulation result includes the minimum and maximum values due to the variation in ν coll = 0 -10 5 Hz.
The most notable result is that for at Fig. 8(a), where WG 1 sccm, and Fig. 8(d), where DC 2 sccm, the results for η w = 0% are in the best agreement with the experimental results. This indicates that the excitation rate distribution in the metastable state is concentrated in the magnetic confinement region. In Fig. 8(b) and (c), the results for η w = 75% (L w = 10 cm) are in the best agreement with the experiment. This indicates that 75% of the metastable neutrals were generated in the waveguide. At z=− 10 to − 15 cm, the density gradient is shown, which indicates that the particles diffuse from the production area at z > − 10 cm. Because the maximum value of the simulation  [31,33], red and blue lines: simulation results. Simulation results include the minimum and maximum values due to variations in ν coll = 0~10 5 Hz agrees with that of the experiment, it is suspected that the electron density is very low in this region. In Fig. 8(e), η w = 25 % (L w = 7 cm) cm is in the best agreement, and in Fig. 8(f), η w = 25 % (L w = 10 cm) cm is in the best agreement. Compared to WG injection, η w inside the waveguide was reduced. Comparing Fig. 8(e) with Fig. 8(f), L w is larger, which suggests that the production area is greater at the bottom of the waveguide.
Using the comparison between simulation and experiment, the excitation rate distribution in the metastable state can be estimated. At a low flow rate, the excitation is concentrated in the magnetic confinement region, and then the excitation transitions to the waveguide as the flow rate increases. Additionally, the particle simulation firstly revealed that the transition of the excitation distribution corresponds to the diffusion of metastable neutral particles, which suggests that the stepwise ionization can broad by the diffusion of metastable state neutral particles. The ratio of stepwise ionization will be investigated by numerical simulation of electromagnetic field of microwaves.

Contribution of the stepwise ionization to plasma density
Here, an electron temperature of 3 eV was assumed, and hence ν iz, m /ν iz, g = 20 was used [23]. The ground-state neutral density was estimated as ν iz, g ≈ 10 3 Hz from the order estimation of n e ≈ 10 17 m −3 , v e ≈ 10 6 m/s, and σ iz ≈ 10 −20 m 2 . Finally, the groundstate neutral density was fixed at n g ≈ 10 19 m −3 from the experimental results near the optimal flow rate (Fig. 6). In these situations, to investigate the contribution of stepwise ionization, n m /n g were varied from n m /n g = 1%-10%. Fig. 9 shows the plasma density and electric field of microwaves in four cases: n m /n g = 0 % , 1 % , 5 % , 10%. The experimental data is the electric field of the microwave near the optimal flow rate (2 sccm) in WG injection [35]. First, we note that the simulated electric field can reproduce the experiment. Thus, because the simulation case n m /n g = 5% is in the best agreement with the experiment, this suggests that the experimental plasma density can be estimated to be on the order of n e ≈ 10 17 m −3 . The microwaves cannot be transported to the discharge chamber when n m /n g > 5%. Of course, even if stepwise ionization is not considered, the electric field can be reproduced when there is a large ground-state neutral density. In this case, if n Ã g ¼ n g þ ν iz;m =ν iz;g Â n m ¼ 2 Â10 19 m −3 is used, the source term (right-hand side of Eq. (7)) has the same value; Fig. 9 Results for (a) plasma density and (b) electric field of microwaves for various ratios n m /n g . The experimental electric field was measured by EO probe near the optimal flow rate (2 sccm) in WG injection [55] thus, the same electric field can be obtained. However, we emphasize that near the optimal flow rate, such a case did not occur in the experiment. As shown in Figs. 6 and 8, the metastable neutral density increases with respect to the flow rate. By contrast, in WG injection, the ground-state neutral density decreases near the optimal flow rate (WG: 1.5-1.8 sccm), and in the DC injection, the density is approximately constant around n g ≈ 1 × 10 19 m −3 (DC: 2.7-3.1 sccm). Thus, the change in the metastable and ground state neutral density indicates that n m /n g increases near the optimal flow rate. Under such conditions, the increase in plasma density near the optimal flow rate cannot be explained without the stepwise ionization.

Spatial structure of plasma mode-transition
Based on the simulation and experimental results, we propose the spatial structure for the plasma mode transition shown in Fig. 10. Specifically, in this spatial structure, plasma mode transition consists of the following three steps: (i). At low flow rates (WG: 1sccm, DC: 2 sccm), direct ionization and excitation occur in the magnetic confinement region. Then, the generated metastable neutral particles are transported to the inside of the waveguide and its exit as shown in Fig. 7. (ii). Stepwise ionization and excitation occur in and at the exit of the waveguide by the transportation of the metastable neutral particles to the waveguide, and more particles are transported into the waveguide. (iii). When the plasma density inside the waveguide becomes high enough to prevent the propagation of microwaves, plasma mode transition occurs.
Step (i) is justified by the good agreement in the LAS measurement results at low flow rates ( Fig. 8 (a) and 9(d)).
Step (ii) is also consistent with the comparison between simulation and experiment ( Fig. 8 (b) and 9(e)), and can be explained by the sharp increase in the spontaneous emission intensity around the optimal flow rate (Fig. 6).
Step (iii) is consistent with the comparison between simulation and experiment ( Fig. 8 (c) and 9(f)), and is consistent with the electric field of the microwaves measured by the EO probe [35].
Additionally, using steps (i)-(iii), we discuss the difference between WG and DC injection. WG injection enhances steps (ii) and (iii) because of the larger ground-state neutral density (Fig. 6), which causes plasma mode transition at lower flow rates. Then, focusing on the axial dependence of the spontaneous emission intensity, we find that the ground-state neutral density distribution in WG injection induces a large ratio in excitation rate η w as shown in Fig. 8. Because of the higher metastable neutral density, the maximum value of spontaneous emission decreases by only approximately 48% from z = − 7 cm to − 16 cm in the WG injection, whereas it decreases by approximately 25% in the DC injection. Therefore, steps (i)-(iii) can explain the difference between WG and DC injection in terms of ground state neutral density, spontaneous emission intensity, and metastable neutral density.
Finally, we discuss whether the proposed mechanism of plasma mode transition can be applied to different types of magnetic field or other applications. Because the diffusion of metastable state neutral particles from magnetic confinement region occurs in any magnetic field, it is considered that proposed mechanism of plasma modetransition can be applied to other different types of magnetic field, such as multipole ECR plasma sources [37,38]. Another thing that some ECR plasma sources are used to generate multiply charged ions [55,56]. Then, some paper reported that stepwise (step by step) ionization from lower charged ions is important to generate multiply charged ions [55,56]. In this situation, it is considered that the diffusion of lower charged ions, i.e., ion leakage, is more important than that of metastable state neutral particles. In summary, our proposed mechanism should be effective under any magnetic field in ECR plasma sources where single-charged ions are dominant.

Conclusion
In this study, to investigate the stepwise ionization distribution from the metastable state that can induce plasma mode transition, one experimental and two numerical Fig. 10 Spatial structure of plasma mode transition. (i) At below optimal flow rates, direct ionization and excitation occurs in the magnetic confinement region. (ii) The generated metastable neutral particles are transported to the inside and exit of the waveguide. From the transported metastable neutral particles, stepwise excitation and ionization occurs in and at the exit of the waveguide, and more particles are transported inside the waveguide. (iii) When the plasma density inside the waveguide is high enough to prevent the propagation of microwaves, the plasma mode transition occurs approaches were performed. First, to investigate the ionization process in two types of injection, the ground-state neutral density and spontaneous emission intensity were measured. Then, a metastable neutral particle simulation including de-excitation was used to estimate the excitation rate distribution in the metastable state. From this estimation, it was found that at a low flow rate, the excitation is concentrated in the magnetic confinement region, and the excitation then transits to the waveguide as the flow rate increases. Second, the electric field of the microwaves was simulated to estimate the contribution of stepwise ionization to plasma density. It was revealed that the simulation results can reproduce the measured data. The experimental results for n m /n g demonstrate that the increase in plasma density near the optimal flow rate cannot be explained without stepwise ionization.
These experimental and simulation results indicate that the location of the stepwise ionization transits to outside of the magnetic confinement region because of the diffusion of metastable particles. Additionally, this transition enhances the plasma density inside the waveguide, resulting in plasma mode transition. This finding can explain the difference between the two types of gas injection from the ground-state neutral density distribution and the axial dependence of spontaneous emission intensity. We believe that this finding will lead to the development of future performance improvements, and our proposed mechanism should be effective under any magnetic field in ECR plasma sources where single-charged ions are dominant.
Here, k 2 c ¼ k 2 z −k 2 m , where k 2 m ¼ μ 0 ε 0 ω 2 m and then boundary condition; ∂ g H z0 ðzÞ=∂r ¼ 0 must be satisfied at r = R. The analytical solution of Eq. (11) is Bessel function and then it is satisfied with following condition.
where ρ m, n is the n-th positive root of ∂B s, m (ρ n )/∂r = 0, where B s, m is m-th the Bessel function. By using this approximation, the transverse part of Laplacian can be expressed as follows.
Thus, proposed approximation means that the transverse part of Laplacian is estimated under approximation of Eq. (9) and then e E y ðzÞ is calculated by the Eq. (14).
Thus, if the change of amplitude on the electric field is large, the approximation generates discrepancy. In the absence of plasma ( e J y ¼ 0Þ, the electromagnetic field of microwaves in central axis is shown in Fig. 11. The electric field of the microwaves in the Fig. 11 Parallel component E y of the electric field of the microwaves in the central axis in the absence of plasma. The quasi-1D approximation are compared with the results of the 3D FDTD simulation and EO probe from [35]