Simple model of multi-scale and multi-site emissions for porous ionic liquid electrospray thrusters

Electrospray thrusters, which extract ions directly from a liquid propellant without generating plasma, have attracted considerable attention for small spacecraft systems because of their high propulsion efficiency, even in low-power operations (e.g., a few watts). Electrospray thrusters were first developed as colloid thrusters,¹ which achieve large thrust but low specific impulses. Recently, ionic liquids (ILs), which are molten salt at room temperature, have been developed.² If thrusters emit ions from the IL of the propellant at higher speeds, they achieve higher specific impulses with smaller thrusts.³ This operation mode, which is called “purely ionic regime,” is useful for long-term operation.^4–6 These IL electrospray thruster technologies have introduced the possibility of large increments of spacecraft velocity for maneuver even for small and low-power spacecraft.⁷ 

Three basic types of IL electrospray thrusters with different emitter structures have been proposed: capillary,^8–10 externally wetted,^3,11,12 and porous.^13,14 The porous type can be easily manufactured via mechanical machining^15–18 and allows the propellant feeding system to be simple and small. It can also easily achieve a purely ionic regime.^19,20 In addition, it can utilize multi-site emission, which enhances the thrust. The emission region is large on the needle tip, with many emission sites. The multiple menisci may be generated by periodic internal pressure owing to various pores or instability of ILs on the porous surface.^21–23 Thus, the relatively large tip curvature of the emitter needle tip potentially induces multi-site emission and increases the thrust, rather than a completely sharp emitter geometry.^18,21

There are large differences in the physics scale in porous IL emission. Figure 1 shows the multi-scale of porous IL electrospray thrusters: (a) an entire thruster head (⁠1 cm) has many porous emitters, (b) one porous needle typically has a size less than 1 mm, and (c) IL menisci with sizes of – appear on the porous needle tip, leading to the multi-site emissions. There are 10³- to 10⁶-order size differences between the emitter needle and menisci.

An empirical model,²⁴ a statistical prediction model,²⁵ and a combination of numerical simulations^26,27 have been used to describe the emissions. However, the physical modeling of multi-site emission is unclear. Accurately predicting the current emitted from porous emitters and reproducing the current vs voltage (⁠–⁠) characteristics have been challenging owing to the complex geometry of IL menisci on the porous emitter needle. Optical observation of the menisci is close to or below the diffraction limit, and electron microscopy also resulted in limited success because high-energy electrons caused the solidification of the IL.²⁸ 

In this study, the proposed model predicted the current emitted from the porous emitter using the surface electric field on the needle, which can be calculated as the solution of Laplace’s equation. The current–voltage (⁠–⁠) characteristics are reproduced using the model, avoiding multi-scale complexity. For the modeling, we focus on the assumption that ionic conductivity determines ion emission rather than ion evaporation. References 29 and 30 have reported that conductivity significantly affects ion emissions.

The proposed model for describing emitted current density incorporates ionic conductivity suppressed by porous media (model parameter) with the surface electric field on the porous needle and the threshold surface electric field for emission. The total emitted current can be calculated via the area integration of the emitted current density ⁠. The conductivity was evaluated using the calculated surface electric field and the measured total emitted current ⁠. The threshold electric field was also obtained via experiments. The detailed model equations are described in Sec. II.

The model successfully simulated the – characteristics using a constant ionic conductivity for 7 m class porous emitters, where the expansion of the emission area resulted in the curve shape corresponding to the multi-site emission. The proposed physical model also explains the role of porous media in emission. The ion–pore wall friction stress restricts the ionic conductivity in the porous media, and the conductivity determines the emitted current.

Figure 3 shows the typical – plots produced by the aforementioned model [Eqs. (16) and (17)] when varying the parameter / in Eq. (16). The plots employ an idealized geometric configuration, listed in Table I, with the tip of the emitter needle located in the extractor grid plane (emitter–extractor gap = 0). Here, is a constant (⁠ = 8.91 for EMI-DCA), whereas is increased from 10–100 S/m in increments of 10 S/m. Additionally, these plots assume a certain fixed threshold V/m because emission starts when ⁠. This ion-emission threshold is discussed in Sec. IV. As the surface electric field increases with an increasing voltage, this model indicates that the expansion of the emission area (⁠⁠) influences the shape of the – curve. Moreover, a decrease in ⁠, which represents the ion-transport suppression, leads to a reduction in the emitted current, as shown in Fig. 3, where was determined using the experimental data presented in Sec. IV.

The thruster head consists of a porous emitter and an extractor grid, where porous materials of two different internal scales are used. The first is a 50 nm class porous material manufactured via a standard acid-leaching process,⁵³ and the second is a 7 m class porous material manufactured via a sintering process (by Nippon Electric Glass Co., Ltd.). Figure 4(a) shows the surface geometry of the emitter array measured with a laser scanning microscope (Keyence, VK-X3000) and design parameters, where represents the height of the emitter needles, represents the radius of the needle tips, represents the half-apex angle of the emitter needle, and represents the radius of the extractor holes. Figure 4(b) shows the photograph taken downstream from one emitter needle and an aligned extractor grid. The alignment error between the needle tip and the center of the extractor grid holes is within 0.03 mm.

The values of these design parameters are presented in Table I. We fabricated multiple needle shapes with a height of 1.0 mm and spacing of 1.0 mm via computer numerical control machining. One emitter array has 76 needles in 1.0  1.0 cm porous media. The grid gap is adjusted using shim rings between the thruster housing and the extractor grid.

The – characteristics were measured to determine the total current from the porous emitter. Figure 5 shows the experimental setup used to apply voltages and measure currents. The porous emitter was connected to a source meter (Keithley, 2657A) through a 1 MΩ resistor via the IL, and the extractor grid was connected to another source meter (Keithley, 2612B) through a 100 kΩ resistor. The extractor voltage was fixed at 0 V, and a 1.0 Hz bipolar pulse voltage was applied to the emitter to measure the – characteristics. The emitter current was obtained by averaging 60 points with a measurement interval of 10 ms at each applied voltage.

An IL of 1-ethyl-3-methylimidazolium dicyanamide (EMI-DCA) with a molecular weight of 177.23 g/mol was used as the propellant. All measurements were conducted in a vacuum chamber at 5.0 × 10^-3 Pa or less. Notably, we only focused on the emitter current to validate the proposed model, but it was also confirmed that the maximum extractor current was less than 10% of the emitter current.

Figures 6 and 7 show the – characteristics for grid gap with 7 μm class and 50 nm class emitters, respectively. The scattered plots are experimental, and the error bars represent the standard deviation of repeated consecutive data. The emitter current for the 50 nm class and 7 m class emitters has an error of approximately 10%–20% in consecutive operation.

The line plots in Figs. 6 and 7 are calculated using the proposed model [Eqs. (16) and (17)]. The parameter / in Eq. (16) was adjusted to fit the calculated – characteristics to the experimental data, thereby allowing the determination of through the least squares method. Note that (EMI-DCA) is constant, and the parameter estimation is discussed in Sec. IV C. The value of the least squares method was evaluated, indicating the validity of the fitting. All the values were larger than 0.95.

In the proposed model, the emitted current density was proportional to according to Eq. (16). Thus, if the emission area is constant, the emitted current is a linear function of emitter voltage ⁠, and the – characteristics should be linear, where is proportional to ⁠. However, the emission area expands with an increase in ⁠, and the emitted current is determined via the area integration of the expanding emission area according to Eq. (17). This emission-area expansion increases the increment rate of the emitted current with an increase in ⁠. Consequently, the expanding emission area in the increase in the emitter voltage results in the curve shape of the – characteristics.

In previous research, two types of – characteristics have been reported. In Ref. 18, which considered a geometry similar to that of the emitters examined in this study, the – characteristics exhibited a gradual curve shape. The emission area expands with an increase in the emitter voltage. This thruster has a relatively large emitted current owing to multi-site emissions. In Ref. 54, the curve shape of the – characteristics changes to a straight line shape in the middle of the range of the increasing emitter voltage. This transition of the – characteristics shape indicates that the emission area expands up to the middle of the range of the increasing emitter voltage; however, emission-area expansion is restricted, and the emission area becomes constant as the emitter voltage increases further.

For the 7 m class emitter, the shape of the – characteristics is successfully reproduced by the proposed model, and the shape is curved with an increase in ⁠, as shown in Fig. 6. This shape of the – characteristics of the 7 m class emitter corresponds to the case of Ref. 18.

For the 50 nm class emitter, a curve shape is observed up to the middle of the range of the increasing and a slightly straighter shape is observed thereafter, as shown in Fig. 7. The shape of the – characteristics of the 50 nm class emitter corresponds to the case of Ref. 54. The emission-area expansion for the 50 nm class emitter may be restricted in the middle of the range of the increasing ⁠.

Figure 8 shows the experimental results of the current vs electric field (⁠–⁠) characteristics with 7 μm and 50 nm class emitters for the grid gap. The error magnitude of the repeated experiment is approximately equivalent to the size of the marker [Fig. 8(a) 1 A and Fig. 8(b) 0.1 A]. The horizontal axis indicates the maximum value of the surface electric field at the needle tip in the axial -direction, and the vertical axis indicates the emitter current per needle ⁠. The emission starts when a certain electric field is applied, and the emitted current increases with ⁠. The starting point of the curve is evaluated as ⁠. Figure 8 indicates that is determined to be 20 V/μm regardless of the polarity and porosity. This value was used for the data for the proposed model [Eq. (17)].

Figure 9 shows the two-dimensional distributions of the magnitude of the electric field near the needle tip and the emission area for different grid gaps and emitter voltages. The emission area is defined as the location where the local exceeds in the proposed model [Eq. (16)], where the direction of is perpendicular to the surface (see Fig. 11 in the  Appendix). The emission area increases with an increase in the emitter voltage and a reduction in the grid gap, as shown in Fig. 9. In Fig. 8(a), the emitted current increases with a reduction in the grid gap. Thus, as the gap becomes smaller, the emission area increases, and a large current is emitted from the large emission area (i.e., multi-site emission). If currents are emitted from only one point at the tip (i.e., single-site emission), the – characteristics should be identical among all the grid gaps.

In Fig. 8(b), the – characteristics are almost invariable among all the grid gaps. When a 50 nm class emitter is used, the expansion of the emission area may be restricted in the middle of the range of the increasing ⁠, as mentioned in Sec. IV A. It is considered that the ion transport is insufficient upstream from the porous needle surface in the case of the 50 nm class emitter.

Table II presents the estimated parameter ⁠. This conductivity was obtained via the fitting process described in Sec. IV A. Thus, the line plots shown in Figs. 6 and 7 represent the values obtained using the numerical model [Eqs. (16) and (17)], employing the parameters listed in Table II.

In addition, the ion–pore wall friction coefficient can be evaluated because and are described with as using Eqs. (6) and (8). Table II presents the estimated (as ⁠).

For the 7 m class emitter, the of anions and cations can be considered constant within 14% and 6% from their mean values in ⁠. The maximum error of consecutive experiments is approximately 20% of the mean values. Thus, is considered constant and independent of the grid gap ⁠. This result of the constant implies that there is sufficient ion transport to the porous surface as global flow in the porous emitter and that the emission area expands with an increase in ⁠, as discussed in Secs. IV A and IV B.

Table II indicates that the ion–pore wall friction is two or three orders of magnitude larger than the ion–ion friction ⁠. Thus, the ion–pore wall friction determines the conductivity, and the conductivity determines the emitted current. Here, in particular, ions are current carriers in ILs, and conductivity is based on the mobility of ions as Eq. (8). Moreover, the suppression of conductivity by ion–pore wall friction may play an essential role in ion evaporation (purely ionic regime) because this model indicates that conductivity is related to the flow rate as ⁠, where is a cross section of the flow in a porous emitter, and is defined as ⁠.^39,42 The current emitted as the droplet is extinguished and evaporated determines the emitted current when the flow rate is sufficiently small.^42,55 Thus, a trade-off may exist: a larger porous scale increases the emitted current owing to the large but may induce a droplet regime. This trade-off implies the existence of an optimal porosity for the porous emitter for the large emitted current and the purely ionic regime.

In the case of the 50 nm class emitter, depends on experimental conditions. This non-constant implies that the emission-area expansion may be restricted in the middle of the range of the increasing because of insufficient ion transport to a porous surface, as discussed in Secs. IV A and IV B. The insufficient ion transport as global flow in porous media arises from the porous structure.⁵⁶ The non-constant for the various values is explained as follows. The emitted current is calculated via area integration on using Eq. (17). However, the actual emission area of the 50 nm class emitter is restricted and smaller than from the middle of the range of the increasing ⁠. Thus, when is small and is calculated to be larger than the actual emission area, the model parameter is underestimated as a small value using the measured current ⁠. In addition, is overestimated for small values.

The emitted current of the 7 m class emitter exceeds that of the 50 nm class emitter in Figs. 6 and 7, which is explained by the smaller in the former case. Additionally, the emitted current of the anions is smaller than that of the cations in Figs. 6 and 7. The difficulty in emitting anions has also been reported.¹⁸ This smaller emission current of anions compared with that of cations is also explained by the larger of anions compared with cations. However, the anion emission may be affected by the field emission effect. When a strong negative electric field is applied (⁠⁠), the electron cloud is formed around menisci.⁵⁷ The electron cloud shields the electric field and the surface charge and weakens the electric field in the IL, reducing the emitted current density. Thus, is evaluated as a smaller value and is evaluated as a larger value in anion emission compared with the cation emission.

The proposed model was also applied to experimental data from previous studies. Figure 10 shows its application to the experimental data in Ref. 22, including the plot of a Bayesian estimation model from Ref. 25 for comparison. By analyzing the emitter geometry employed in Ref. 22, we found that the tip curvature is not axially symmetric and varies in the azimuthal directions, that of tip 1 ranging from 17 to 53 m and tip 2 ranging from 15 to 17 m. However, to calculate the electric fields for the analysis, we assumed the following tip curvatures:  m for tip 1 and  m for tip 2. Moreover, the error bars obtained via repeated experiments and the manufacturing accuracy of the emitter geometry were uncertain in Ref. 22. Thus, the evaluation values obtained from Fig. 10 should be considered as approximations (tip 1: and 216 S/m for the anion and cation , respectively; tips 2: and 1192 S/m for the anion and cation, respectively). Thus, the differences in the tip curvatures of tips 1 and 2 resulted in significant differences in the values. As shown in Fig. 10, the proposed model, despite its simplicity, reproduced the – shapes as effectively as the Bayesian model.²⁵ 

A simple model for multi-site emission [Eq. (16)] was proposed based on momentum conservation of transported ions through porous media, with the assumption that conductivity controls ion emission.²⁹ Using the proposed model, the – characteristics of experiments were reproduced successfully with a 7 m class emitter. It is difficult to apply the model to the 50 nm class emitter because ion transport may be insufficient, and emission-area expansion is restricted with the increasing emitter voltage.

The proposed model indicates that the expansion of the emission area (i.e., multi-site emission) causes the unique curve shape of – characteristics.

The ion–ion friction and ion–pore wall friction determine ⁠. Moreover, the ion–pore friction is the primary factor determining because is two orders of magnitude larger than ⁠. The larger emitted current for the 7 m class emitter compared with that for the 50 nm class emitter is explained by the smaller ion–pore wall friction.

The proposed model assumes that the surface charge converges on the order of magnitude of owing to the quasi-neutral region and the non-neutral sheathe (double layer) when a strong electric field is applied to the IL (⁠⁠). The boundary conditions considering the non-neutral region should be investigated in future research.

This work was partially supported by JSPS KAKENHI under Grant No. JP21H01530, the Canon Foundation, the Advisory Committee for Space Engineering of Japan Aerospace Exploration Agency, Grant-in-Aid for Exploratory Research in the area of human space technology and space exploration, and JST FOREST Program under Grant No. JPMJFR2129. The authors acknowledge Nippon Electric Glass Co., Ltd. for offering porous materials, thank Mr. Takase for fabricating thruster prototypes at the Advanced Machining Technology Group of Japan Aerospace Exploration Agency, and thank Dr. Kazuma Emoto for his insightful reviews and valuable suggestions regarding this paper.

The authors have no conflicts to disclose.

Koki Takagi: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Project administration (equal); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Yusuke Yamashita: Conceptualization (equal); Formal analysis (supporting); Methodology (supporting); Software (supporting); Validation (supporting); Writing – review & editing (equal). Ryudo Tsukizaki: Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (equal); Writing – review & editing (equal). Kazutaka Nishiyama: Project administration (equal); Resources (equal); Supervision (equal). Yoshinori Takao: Funding acquisition (equal); Project administration (equal); Resources (equal); Software (supporting); Supervision (equal); Validation (supporting); Writing – original draft (supporting); Writing – review & editing (equal).

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

Equation (19) is discretized using second-order central difference methods applicable to the non-surface region, as depicted in Fig. 11(a). A cut-cell method⁵⁸ was applied to calculate the surface electric field on the exact boundary in the surface region depicted in Fig. 11(a). Figure 11(b) shows the control volume for the cut-cell method. ⁠, ⁠, ⁠, and represent the electric fields defined on the left, top, bottom, and right interfaces of the control volume, respectively; represents the distance between cell centers; and represents the distance between the exact emitter surface and the cell centers in the surface region.