Electrospray propulsion: Modeling of the beams of droplets and ions of highly conducting propellants

Electrospray thrusters are electrostatic accelerators of charged droplets and ions generated by electrosprays.¹ The low thrust associated with a single emitter (of the order of 1 $μ$N) combined with its high stability enables precision spacecraft positioning applications. This technology has been demonstrated by the Disturbance Reduction System-Space Technology 7 mission (DRS-ST7), a precursor of the Laser Interferometer Space Antenna mission (LISA).^2,3 Furthermore, higher thrust can be achieved by using arrays of electrosprays,^4–7 enabling efficient primary propulsion and attitude control at the low power levels available to SmallSats. Overspraying of the extractor and accelerator electrodes is a major lifetime-limiting mechanism in electrospray propulsion.⁸ The progressive accumulation of fluid on these electrodes leads to the formation of electrically conducting films and the catastrophic shorting of power supplies; it may also lead to backspraying between the grids and the emitter electrode and erratic and progressive worsening of the performance of the thruster.⁹ Having a detailed first-principles understanding of the physics of the beam and the capability to model, it is important to properly design the electrodes and to predict potential impingement of the beam during the long periods of operation, often several years, typical of electric thrusters.¹⁰ 

Electrosprays operating in the cone-jet mode¹¹ atomize the liquid propellant into submicrometric droplets¹² with narrow distributions of diameters and charge-to-mass ratios. This is the emission mode commonly used in electrospray propulsion. In an experimental setting, the liquid is typically fed through a capillary tube at a higher potential than a nearby plate, and the electric field shapes the liquid into a conical meniscus as it flows out of the tip. A thin jet forms at the vertex of the so-called Taylor cone and, through a breakup caused by capillary instability, generates charged droplets. Molecular ions are also emitted from the droplets and possibly from the surface of the jet as well.¹³ The electric field accelerating the droplets and ions is induced by both the electric potentials set at the electrodes and the charges carried by the particles in the beam. The contribution of this space charge is key to the radial expansion of the beam.¹⁴ Particle-on-particle Coulombic repulsion dominates the initial region of the beam following the emission point.¹⁵ This initial region is small compared to other characteristic lengths of the electrodes such as the gaps between the emitter, extractor, and accelerator electrodes and the diameters of the emitter and the orifices of the extractor and accelerator.

The sprays of the charged droplets produced by electrosprays have been modeled mostly under atmospheric conditions.¹⁶ It was first modeled by Gañán-Calvo et al.¹⁷ using a Lagrangian approach, but the set of ordinary differential equations could not be solved for a large enough number of droplets to represent the space charge effect. By the end of the 2000s, simplified Lagrangian models and more computational power enabled the simulation of electrosprays emitting thousands of droplets,^18,19 and the models were used to predict deposition patterns.²⁰ The computational time scales with the square of the number of droplets due to their mutual interaction. Therefore, simulating electrosprays of submicrometric droplets is challenging due to the very large number of droplets present in typical regions of interest. This problem can be alleviated by using different integration time steps throughout the computational domain.²¹ An Eulerian model²² can provide insights and realistic results at a fraction of the cost of a Lagrangian simulation but cannot fully capture the coupling between moving droplets. Gamero-Castaño¹⁵ has developed an Eulerian model for electrospray beams under vacuum conditions.

The model presented in this paper combines two different approaches: a Lagrangian model for the inner region near the jet breakup where the trajectories of all individual charged droplets are integrated simultaneously and where Coulomb repulsion between droplets is fully captured while integrating the equations of motion and an Eulerian model for the outer region where the structure of the beam is obtained by computing the envelopes of beamlets in which the population of droplets is subdivided. The two regions are coupled by the initial conditions of the envelopes, which are calculated with the solution of the inner problem and by the influence of the space charge of each region on the electric field. This strategy greatly reduces the computational time while retaining the important particle-on-particle Coulomb interaction where is significant. Finally, the numerical solution is validated with experimental results obtained with the ionic liquid 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide, EMI-Im,²³ which is the propellant used by the electrospray thruster recently demonstrated by the DRS-ST7 mission and whose beams have been characterized in detail.^24–28 EMI-Im produces a particular type of beam consisting of a mixture of charged droplets and molecular ions without significant spatial separation. There exist other prototypical cases such as ion beams devoid of charged droplets,²⁹ beams of main and satellite droplets that are spatially separated,¹⁵ electrosprays at atmospheric pressure,³⁰ etc., that can be studied with the present model by using the appropriate initial conditions and distributions for the constituent particles.

The model considers a steady and axisymmetric electrospray. Figure 1(a) shows a sketch of the computational domain in cylindrical coordinates. The emitter electrode has a diameter of 0.36 mm and a conical tip with a 45$°$ half-angle. The extractor electrode is located 0.73 mm from the vertex of the emitter, has a thickness of 0.45 mm, and has an orifice with a diameter of 1.78 mm. This geometry matches an electrospray source previously used to characterize beams of the ionic liquid EMI-Im,²⁵ and these experiments will help validating the model. Beams of EMI-Im have also been measured by Thuppul et al.²⁷ The red grid is used to discretize the space charge and evaluate the electric field at the nodes and only needs to extend over the region occupied by the beam. The breakup of the jet produces charged droplets of varying diameters $D$ and charge-to-mass ratios $ξ$⁠, with a joint density distribution $f(D,ξ)$⁠. The vertex of the grid is placed at the nominal point where the jet breaks into droplets. The position $xi$ of the $i$th droplet can be computed by integrating its equation of motion,

where the electric field $E$ is induced by the distribution of charges on the surfaces of the electrodes, on the cone-jet and throughout the beam. Although the breakup of the jet is time-dependent, we consider constant nominal values for its position $zJ$ and for the velocity of fluid $vJ$⁠. Moreover, for the highly conducting liquids of interest to electrospray propulsion, the jet has a substantial amount of charge on its surface, making it susceptible to lateral oscillations.²⁶ Its effect on the initial conditions of the droplets is implemented by assigning values that randomly deviate from the nominal breakup position and jet velocity,

where $αi$ and $βi$ are small offsets and $iz$ is the unit vector along the axial direction. The droplets are inserted sequentially in the computational domain at times $to,i$ fixed by the conservation of mass

where $Di$ and $DJ$ stand for the diameters of the $i$th droplet and the jet at the breakup, respectively. Finally, the solution of the Poisson equation for the electric potential $ϕ$

yields the electric field, $E=−∇ϕ$⁠. $q$ is the charge of a droplet, $δ$ is the Dirac-delta function, and $ε0$ is the vacuum permittivity. As shown in Fig. 1(a), to solve the Poisson equation, we use as boundary conditions the potential of the emitter electrode, zero normal electric field at the radial-azimuthal plane upstream of the emitter, and zero potential at the extractor and all other boundaries.

Although (1)–(5) provide a full description of the beam of droplets, the integration of these equations is impractical due to the large number of droplets, the many time steps needed in most situations of practical interest, and the need to solve the Poisson equation at every time step. Instead, we use a model that separates the computational domain in two regions: an inner, smaller zone surrounding the emission point, characterized by a high droplet density, and where (1)–(5) are solved for a large sequence of droplets (i.e., in this inner region we retain individual particle-on-particle Coulomb interactions) and an outer, larger region where the lower droplet density and the reduced importance of Coulomb scattering are compatible with a continuous model for the droplet velocity and charge density. The two solutions are matched at the common boundary, where the inner solution provides initial conditions for the outer problem. Figure 1(b) shows the extent of the inner region in a typical simulation. Its characteristic length is 10 $μ$m, about two orders of magnitude smaller than the distance between the emitter and the extractor. This characteristic length can be evaluated with the solution by finding the distance from the emission point where the radial component of the electric field induced by droplet-on-droplet Coulomb repulsion is a small fraction of the total value.

To solve the continuum, outer problem, we separate the droplets into groups with constant charge-to-mass ratio $ξk$ and compute the volumetric charge density $ρk$ and velocity $Vk$ fields for each group. These fields fulfill conservation of charge and momentum

where the Poisson equation for the electric potential

yields the electric field. The space charge term includes contributions from all droplet groups in the outer region, $∑ρk$⁠, as well as the space charge in the inner region, $ρin$⁠. We use the method of characteristics to integrate (7) for a given electric field, where the trajectory of any droplet with charge to mass ratio $ξk$ is a characteristic curve of this partial differential equation. Once a set of characteristics is computed, the axial $Vz$ and radial $Vr$ components of the velocity field are obtained by integrating the ordinary differential equations

along the characteristics. The characteristics themselves are determined by integrating (1), starting at the boundary between the inner and outer regions. The required initial velocities are obtained by averaging the inner solution, which also yields the charge flux along the boundary. For example, the initial velocity $Vo,k$ of a characteristic with a charge-to-mass ratio $ξk$ at a particular position of the boundary is computed as the charge-weighed average of the velocities of all $N$ individual droplets with that charge-to-mass ratio exiting within a small neighborhood of that position,

With this information, we compute the two characteristics enveloping each droplet group, as well as additional characteristics in between. Once the velocity field $Vk$ is determined along the characteristics, Eq. (6) readily yields the charge density field.

The sequence of droplets needed for the inner solution must be representative of the variation of sizes and charge-to-mass ratios produced by the jet breakup. Although the resulting distributions have been measured for micrometer-sized and slightly smaller droplets,¹⁵ only partial experimental information is available for highly conducting propellants. In particular, the charge-to-mass ratios can be measured with the time-of-flight technique, but the diameters of such small droplets cannot be fully characterized. In the absence of detailed information, we assume the distribution

where $g(ξ)$ is obtained from time-of-flight experiments and the mean of the normal distribution depends on the charge-to-mass ratio, $Dc(ξ)$⁠. In the calculations, we use a typical standard deviation $σ$ of 10% of $Dc$⁠.³¹ 

When characterizing a beam, both the time-of-flight $τtof$ and retarding potential $ϕRP$ of a droplet must be measured to obtain its charge-to-mass ratio,

Here, $Ltof$ is the length of the time-of-flight path of the detector. $τtof$ and $ϕRP$ can be measured simultaneously with analyzers operating in tandem, but such configuration only samples a small region of the beam.²⁶ Instead, we use a time-of-flight measurement of the full beam, together with the average value of the retarding potential, to estimate $g(ξ)$⁠. As an example, Fig. 2(a) shows the time-of-flight curve of an EMI-Im beam with a current of 305 nA.²⁵ The droplet and ion populations are easily separated by their very different velocities and times-of-flight, the ions representing 15.4% of the total beam current. The average retarding potential of the droplet population is 1471 V, and the time-of-flight path is 0.147 m. When using these value in (13), the time-of-flight curve for the droplets is directly converted into the charge-to-mass ratio cumulative distribution shown in Fig. 2(b) and into the density distribution $g(ξ)$ upon differentiation. For the purpose of numerical calculations, the distribution is discretized into a large number of bins with the centers of the bins shown as red dots on the distribution curves.

We use linear instability analysis of capillary breakup to estimate $Dc(ξ)$⁠.²⁶ The small electrical relaxation time of EMI-Im compared to characteristic breakup times makes the process quasi-electrostatic. Furthermore, these jets have very high values of the Ohnesorge and Taylor numbers (viscous diffusion is very efficient, and the electrostatic pressure is always large, in some cases exceeding the capillary pressure). Under such conditions, the diameter $D∗$ of the most likely droplet produced by the breakup (i.e., the diameter with maximum growth rate derived from the instability analysis) is approximately given by Ref. 26 

where the Taylor number is

$IB$⁠, $ς$⁠, and $γ$ stand for the beam current, surface charge density of the jet, and the surface tension of the liquid, respectively. The charge to mass ratio of a droplet scales with the square of its diameter in an equipotential breakup. In addition, since conservation of charge and mass requires the average charge-to-mass ratio of the most likely droplet to be similar to that of the jet, we can approximate the sought mean diameter by

Once the joint $f(D,ξ)$ distribution is determined, the sequence of droplets is constructed in two steps. First, we generate a large pool of droplets (typically, 100 000) by sampling $f(D,ξ)$⁠. Second, subsets of droplets (typically, 15) are randomly extracted from this pool, subject to the condition that their total charge is within a small fraction of that in a section of the jet with the same volume. Subsets that fulfill this condition are added to the sequence and removed from the pool. This ensures that the sequence of droplets fulfills conservation of charge at frequencies comparable to that of droplet emission. Since the net charge in the jet is distributed on its surface, and the electric current transported by the jet is the advection of this surface charge, the net charge contained in a section of the jet of length $L$ is equal to $IBL/vJ$⁠.³² 

To estimate the initial position and velocity of the droplets we use a technique developed by Gamero-Castaño for determining the nominal velocity and potential of the jet at the breakup.^15,33 The technique is based on the natural dispersion of the droplets’ charge-to-mass ratios and assumes that the variations of the potential and the velocity of the droplets at emission are much smaller than the voltage drop along the cone-jet and the velocity gained by the propellant along the jet. Under these conditions, all droplets are emitted at approximately the same potential $ϕJ$ and velocity $vJ$⁠, which are obtained from a regression of the droplets’ retarding potentials and charge-to-mass ratios. Reference 26 tabulates $vJ$ and the voltage drop along the cone-jet, $ϕE−ϕJ$⁠, for cone-jets of EMI-Im operated at 21 $°$C. $ϕE$ and $ϕJ$ are the potential of the emitter and the potential of the jet at the breakup, respectively. The experimental values range between 547 and 431 m/s and 196 and 455 V for beam currents between 230 and 450 nA. The nominal position of the breakup $zJ$ is estimated with $ϕJ$ and a numerical solution of the electric potential along the axis. The values for a cone-jet with a current of 300 nA are $vJ=501$ m/s, $ϕE−ϕJ=274$ V, and $zJ=10.2μ$m. For reference, the diameters of the jet and the most likely droplet are $DJ=27.7$ nm and $D∗=45.0$ nm. All numerical solutions presented in this article are done for this 300 nA EMI-Im beam.

Finally, to incorporate the jet’s lateral oscillations on the initial conditions of the droplets, we randomly set their initial positions within a circle of radius $ro$ centered on the axis (⁠$ro$ should have a value of the order of the diameter of the jet) and use a random oscillation with maximum angle $φo$ to distribute $vJ$ among the three components of the initial velocity

where $Ai$ and $Bi$ are numbers picked randomly from a continuous uniform distribution between 0 and 1 and $θi$ is an azimuthal angle picked randomly from a continuous uniform distribution between 0 and $2π$⁠. $αi$ and $βi$ were introduced in Eqs. (2) and (3).

We use the continuum approach (6)–(8) to model the expansion of the ion beam. Retarding potential measurements show that the ions present in electrosprays of EMI-Im are emitted from droplets in flight within and immediately downstream of the breakup region, and possibly from the surface of the jet.²⁶ This bounds the initial positions of the trajectories. We divide the ion current into groups, each one simulating emission from a stationary droplet. These sources of ions are placed at regular intervals inside the beam of droplets and, based on retarding potential measurements, within a few hundred volts from the breakup. We define $m×n$ characteristics for each group, all starting at the position of the droplet source and with initial velocities simulating spherical and homogeneous emission

The modulus of the initial velocity is derived from the self-potential of the droplet, $ve=2ξionϕself$⁠, with $ϕself=q/(2πε0D)$⁠. Since the breakup is nearly equipotential, the emission velocity is to first approximation independent of the diameter of the droplet. For an electrospray of EMI-Im with a current of 300 nA, our estimates of the droplet self-potential and the ion emission velocity are $ϕself=19$ V and $ve=2688$ m/s. Ion emission from the jet can be simulated as an additional group with characteristics starting at the axis immediately upstream of the breakup point, and with a radial velocity derived from the self-potential of the jet.

We use the boundary element method with linear elements to solve the Poisson equation and evaluate the electric field. Brebbia and Dominguez³⁴ provide an excellent description of the technique, and Bakr³⁵ focuses on its application to axisymmetric geometries such as the one considered here. For a given charge density field, the potential and the electric field at any point inside the simulation domain are computed with boundary integrals involving the potential and its gradient normal to the boundary, and the position of the point. The evaluation of the integral requires the discretization of the boundary into $NΓ$ nodes: of the $2×NΓ$ values of the potential and the normal gradient needed to compute the integral, $NΓ$ values are specified as boundary conditions, while the remaining $NΓ$ values must be computed by considering each node as the load point in the boundary integral. This results in a system of $NΓ$ linear algebraic equations for the unknown potentials and gradients. Using $u$ to denote the vector of unknowns (a subset of the potentials and gradients at the nodes), $w$ for the vector of specified boundary conditions, and $q$ for a vector with the net charge inside the cells of the grid (see Fig. 1), the system of algebraic equations and its solution can be written as

After separating $u$ and $w$ into vectors $ϕΓ$ and $eΓ$ with the potentials and gradients at the boundary nodes, the potential field at any point $x$ is computed explicitly as

with similar equations for the components of the electric field. Note that the potential is two dimensional (axisymmetric). The elements of the matrices $A$⁠, $B$⁠, and $M$ and the vectors $g$⁠, $h$⁠, and $m$ are geometric factors that only depend on the discretization of the boundaries, the mesh used to discretize the space charge and, in the case of $g$⁠, $h$⁠, and $m$⁠, on the evaluation point. Therefore, these matrices and vectors only need to be computed once in an initial step of the calculations.

Equations (19) and (20) assume a continuous space charge field which, combined with the axisymmetric geometry, result in a smooth and axisymmetric electric field. This is appropriate for the outer region of the beam model. However, droplet-on-droplet Coulomb interactions are retained in the inner region, resulting in a three-dimensional electric field. To account for this, we define a spherical neighborhood around each droplet with a cut-off radius $rcf$ and compute the contribution to the field by neighboring droplets as a direct sum of pair potentials

When computing the Coulomb field induced by neighboring droplets we subtract their contribution to the axisymmetric field, in order to avoid double-counting. In summary, we use the smooth electric field (20) when integrating the characteristics of the momentum equation in the outer region and the three-dimensional field (21) when integrating the trajectories of individual droplets in the inner region.

The flow chart in Fig. 3 shows the main steps in the algorithm used for solving the model. In an initial block, we define basic inputs such as the beam current, emitter potential, temperature, domain geometry, and boundary conditions; compute the dimensionless numbers present in the equations and the sequence of droplets to fly; define the grids for discretizing the space charge and evaluating the axisymmetric electric field; and compute the matrices and vectors of constant coefficients needed for evaluating the electric field. This is followed by a main loop where the inner solution is computed and processed to determine the initial conditions for the characteristics of the outer region; and an inner loop is executed for computing the outer solution. In both the main and inner loops, the space charge is updated with the latest available solution, and the electric field is reevaluated with this information. The solution converges rapidly, and the main and inner loops can be stopped after a few iterations.

The simultaneous integration of the trajectories of droplets in the inner region is done with the Störmer–Verlet method. We employ a time step $Δt$ based on the emission period of the critical droplet, $τ∗=2D∗3/(3DJ2vJ)$⁠, typically $Δt=0.01τ∗$⁠, which is sufficient to produce a maximum error in the exit angle of the droplets’ trajectories smaller than 1%. The trajectories are computed in three-dimensional cartesian coordinates using the electric field (21). The characteristics of (7) for both the droplet groups in the outer region and the ion groups are integrated in cylindrical coordinates using a fourth order Runge–Kutta scheme with a variable time step.

Figure 4(a) shows a snapshot of droplets in the inner region, depicted as circles with sizes proportional to their diameters. For simplicity, we only plot the radial and axial components of the positions in cylindrical coordinates. Once the droplets cross the boundary between the inner and outer regions, they are removed from the integration, and their exit positions and velocities are averaged to obtain the initial conditions for the outer problem. The average exit velocities of six representative droplet groups are illustrated with vectors spanning the full exit range of each group (11 vectors for each group). A maximum jet oscillation angle $φmax=14.75°$ and a emission circle with $ro/DJ=5.7$ are used to set the initial conditions (Sec. III C will show that these choices yield the best comparison with experiments). The trajectories of most droplets rapidly separate from the axis within a couple of micrometers (otherwise, all droplets would be below the 14.75$°$ line) due to Coulomb particle–pair repulsion with nearby droplets. There is also a segregation of droplets by sizes, with smaller droplets preferentially moving farther away from the axis. Figure 4(b) shows a quantitative analysis of the average velocities and the distribution of droplets at the exit boundary of the inner region, using as abscissa the polar angle along the boundary. The accumulated current data on top shows the fraction of the current that exits the inner region below a given polar angle, for each droplet group. For comparison, the polar angles at which the velocity vectors in Fig. 4(a) are calculated coincide with the abscissas of the dots in the accumulated current curves of Fig. 4(b). Although the segregation of droplets by charge-to-mass ratios is evident (the higher the charge-to-mass ratio the larger the exit polar angle at a current fraction of 50%), the droplet groups overlap significantly. In fact, a small fraction of the droplets with the highest charge-to-mass ratio also appear near the axis of the beam (see distribution for $ξ=8.06$⁠). The bottom of Fig. 4(b) shows the angle formed by the components of the average exit velocities. This angle depends on the exit position and is independent of the charge-to-mass ratio, matching the polar angle up to approximately $φ=28°$ and falling below for higher values.

The electric field determines the trajectories of the droplets. Figure 5 analyzes its radial component in the inner region, which is key for the expansion of the beam. The axial component is mostly determined by the electrodes and is approximately constant within this small region. Figure 5(a) shows the fraction of the radial component due only to direct contributions from neighboring droplets, while Fig. 5(b) shows the fraction of the radial component induced by the space charge everywhere. The direct contribution from neighboring droplets is important only within a small zone near the emission point and decays rapidly away from it, while the contribution of the space charge remains significant further downstream. Thus, although repulsion between charged droplets is key for the expansion of the beam, direct droplet-on-droplet interactions only need to be retained within a few micrometers from the emission point. Elsewhere, the continuum treatment of the space charge is an excellent approximation. Figure 5(c) shows the fraction of the radial component induced by the space charge in the inner region. The comparison with Fig. 5(b) indicates that the coupling between the inner and outer regions is important, and the space charge in the outer region needs to be accounted for when computing the electric field in the inner region. Finally, Fig. 5(d) shows the value of the radial electric field. In the absence of space charge, the field would decrease monotonically to zero toward the axis, but when here is a beam of charged particles the space charge induces a radial electric field that peaks near the axis and is most intense near the emission region.

Figure 6(a) shows three characteristic curves for each droplet group featured in Fig. 4. The characteristics represent the lower envelope of each beamlet, the top envelope, and the envelope enclosing 50% of the current of the group. The characteristics do not intersect because the electric field in the outer region, mostly the result of the potential difference and the geometry of the electrodes, does not have a strong focusing effect. Thus, the angular separation of droplets according to charge-to-mass ratios observed in the inner region extends to the outer region. This could be modified by implementing a different electrode geometry. Figure 6(b) shows the local angle of the characteristics: all characteristics are pushed toward the axis by the electric field, but this effect is more intense in characteristics starting at higher polar angles [see Fig. 4(b)]. Note also that this effect is stronger at increasing charge-to-mass ratio, e.g., compare the characteristics for $ξ/ξJ=2.04$ and 0.6 starting with an angle near 18$°$ because the retarding potential of a droplet decreases at increasing charge-to-mass ratio, $ϕRP≅vJ2/2ξ+ϕJ$⁠, effectively increasing the bending efficiency of a given electric field.

Figure 7 shows the characteristics of the ion beam. We consider emission from droplets in flight, as well as from the surface of the jet upstream of the breakup. We also plot for reference the top, the bottom, and the 50% current envelopes of the droplet group with highest charge-to-mass ratio, $ξ/ξJ=19.3$⁠. Two results are most significant: first, if ions were emitted from the jet just upstream of the breakup, they would be narrowly distributed and concentrated in the outer region of the beam of droplets, enveloping it. Thus, if a significant fraction of the ion current were emitted from the jet, the beam current profiles would display a strong and narrow local maximum at the largest polar angles of its range. This is not observed in experiments with EMI-Im, which instead show that ions appear throughout the beam, down to its axis.^24,25,27 On the basis of this comparison, we conclude that there is not significant emission of ions from the surface of the jet in EMI-Im electrosprays. Second, the characteristics of ions emitted from droplets in flight overlap with the droplet beam, in agreement with experiments. This may seem to contradict the segregation of droplets by charge-to-mass ratios, which if extrapolated to the ions would make them appear fully separated from the droplets, enveloping them. However, the emission of ions and the subsequent interaction with charged particles is markedly different: droplets are emitted with a large kinetic energy per unit charge, mostly directed along the axis, and near the emission region repel away from other droplets with similar velocities. In these scattering events, the droplets with lower inertia and higher specific charge undergo the higher deflections, giving rise to angular segregation. On the other hand ions are evaporated from droplets in flight in every direction (including toward the axis and backwards toward the emitter), with similar initial speeds. Furthermore, Coulomb collisions between ions and droplets are not only unlikely (the droplet number density decreases rapidly downstream of the breakup) but also weak. In addition to the initial emission speed (associated with the droplet’s self-potential of few tens of volts), the ion gains significant kinetic energy during its flight in the potential field, before potentially flying near a second droplet. Thus, if the impact parameter is smaller than the radius of the droplet, the ion simply penetrates the droplet disappearing from the beam; otherwise, its trajectory only suffers a slight deflection. In summary, an ion only undergoes a significant “scattering event,” i.e., the repulsion by the droplet from which it is emitted, which sends it into the beam with a significant velocity and in any possible direction. From this point on the interaction with other charged particles in the beam is well approximated by the interaction with a continuous space charge field.

In addition to the velocity of the fluid and the position of the jet breakup obtained from experiments, small radial displacements and radial velocity components must be specified to fully define the initial conditions of droplets. In the absence of a numerical solution of the breakup dynamics, we rely on (17) which randomly sets these values as functions of two parameters, $ro$ and $φo$⁠. We next compare experimental measurements²⁵ with numerical solutions for different values of $ro$ and $φo$ to determine the effect of these two parameters on the expansion of the beam and to find optimum values. Figure 8 plots the root mean square of the difference between the experimental and computed accumulated current profiles

for values of $ro/DJ$ and $φo$ between 1 and 6.57 and 2$°$ and 15.25$°$⁠. The accumulated current $Iacc(φ)$ is the current transported by the beam between the axis and the polar angle $φ$⁠. The difference is maximum for the smallest values of $ro$ and $φo$ and minimum for $ro/DJ=5.7$ and $φo=14.7°$⁠. It is worth noting that the optimum values are not unphysical: $ro/DJ=5.7$ means that droplets are emitted from within a circle with a radius that is 5.2 times the radius of the critical droplet, while $φo=14.7°$ indicates that the section of the jet that is oscillating is approximately 0.6 $μ$m long, a fraction of the estimated total length of 10.2 $μ$m.

Figure 9 compares the experimental current profile with model results for several values of $ro$ and $φo$⁠. The profiles are shown in both accumulated current form, Figs. 9(a) and 9(b), and in current density $J$ form, Figs. 9(c) and 9(d),

We plot $r2J$ for easier comparison at arbitrary distance from the emitter. The effect of the size of the emission area on the spreading of the beam is stronger at intermediate and large angles, and it does not significantly change the profile near the axis [compare curves with $φo=5°$ and emission radii of 3 and 1.5 in Figs. 9(a) and 9(c)]. The oscillation angle has the opposite effect, changes on $φo$ only affect significantly the profiles at small and intermediate angles [compare curves with $ro/DJ=3$ and oscillation angles of 10$°$ and 5$°$ in Figs. 9(a) and 9(c)]. These dependencies of the beam structure are related to the angular segregation of droplets by size: the smaller the emission area the higher the initial Coulomb repulsion between droplets, which has a larger effect on small droplets (they have less inertia and their higher charge-to-mass ratios make them more sensitive to the electric field) and preferentially pushes them outward. On the other hand, the larger droplets are less affected by Coulomb repulsion, and need an initial radial velocity component to move away from the axis. Figures 9(b) and 9(d) show how the optimal values $ro/DJ=5.7$ and $φo=14.7°$ produce a solution that matches well the measured profile. Figures 9(b) and 9(d) also show the profiles of the droplet and ion beams for this optimal case, as well as the profile for values of the emission area and oscillation angle far from the optimum (⁠$ro/DJ=1$ and $φo=2°$⁠). In the simulations, all ions are emitted from droplets in flight, resulting in a beam with an angular range similar to that of the droplet beam.

The algorithm for solving the model relies on integrating trajectories with an approximate electric field, using these trajectories to update the space charge and the electric field, recomputing the trajectories, and iterating until convergence. Because the space charge in both the outer and inner regions affect the solution in each one, the algorithm uses a nested loop for computing the outer region, within the body of a main loop that computes the overall solution. The convergence of the algorithm is illustrated in Fig. 10 as a function of the global iteration number $n=(n1−1)×3+n2$⁠, where $n1$ is the iteration step in the main loop, $n1=1,2,3$⁠, and $n2$ is the step in the inner loop, $n2=1,2,3$⁠. Figure 10(a) plots two measures of the error of the solution, the root mean square of the accumulated current profiles with respect to the last iteration

and the root-mean-square of the exit angles of all droplet characteristics with respect to the last iteration

The error in the first iteration (⁠$n1=1$⁠, $n2=1$⁠) is 6.0% in the accumulated current profile and 2.1$°$ in the exit angle. The error rapidly reduces to 0.26% and 0.27$°$ in the second iteration (⁠$n1=1$⁠, $n2=2$⁠), i.e., the first time that the outer solution is used to update the electric field in the outer region without recomputing the inner region. An additional iteration of the outer region without recalculating the inner region (third iteration, $n1=1$⁠, $n2=3$⁠) does not significantly improve the solution. The solution improves again when the outer solution is used to recalculate the inner region (fourth iteration, $n1=2$⁠, $n2=1$⁠). However, at this point the error is already small, surely smaller than differences in the solution due to uncertainties in the initial conditions. Figure 10(b) illustrates this by plotting the accumulated current profiles obtained in the first, second, and ninth iterations. The solutions for the second and ninth iterations practically coincide when plotted in this scale, the difference being much smaller than between the solution and the experimental profile (see Fig. 9). The algorithm converges fast, executing two iterations of the inner loop, followed by one final iteration to recompute the inner and outer regions, is sufficient to obtain an accurate result.

We present a model for electrospray beams that solves this problem as a continuum except near the emission region where, due to the high number density of droplets, retaining Coulomb pair-interactions between droplets is essential and a Lagrangian approach must be employed. At the boundary between the Lagrangian and Eulerian regions, the former provides the initial conditions needed for solving the latter. Furthermore, because the space charge in both regions determines the electric field in either one, the algorithm for solving the problem uses an iterative scheme in which the different regions are solved sequentially, to provide estimates of the electric field of increasing accuracy. The computed current density profiles match well the experimental results. The model explains the paradoxical coincidence between the angular ranges of droplets and ions observed in experiments with the ionic liquid EMI-Im. The model shows that this is a consequence of the emission of ions from droplets in flight. If ions were emitted from the surface of the jet they would envelop the droplet beam, in contradiction with experiments.

The model requires the joint distribution of charge-to-mass ratio and diameters of the droplets. It also requires the initial velocity and position of the droplets upon emission from the jet. Time-of-flight and retarding potential diagnostics provide this key information. In addition, the droplets are naturally emitted with small radial offsets in their positions and velocities. We do not have the ability to accurately predict these offsets and instead model them with a random distribution based on two parameters, the maximum radius of the emission area and the maximum oscillation angle of the jet. The comparison between the experimental and calculated current density profiles guides the selection of these two parameters, which are found to have physical values.

Four important problems benefit from the beam model and provide opportunities for extending this research: (a) the design of the electrodes of electrospray thrusters for minimizing the interception of propellant by the extractor and accelerator. Over time, this creates conducting films that short the high voltage sources powering the electrodes, and the catastrophic failure of the electrospray thruster. Note that the presence of an accelerator or additional electrodes downstream of the extractor has a negligible effect on the inner region, but contributes to the shape of the beam in the outer region; (b) investigating the interaction of the plumes of electrospray thrusters with the low density plasma that surrounds a spacecraft, and with its surfaces; (c) the model can be extended to include viscous drag on the droplets, in order to study the deposition and coating of electrosprayed particles under atmospheric conditions and reduced pressures; (d) in addition to beams similar to those of EMI-Im, the model can be used to study other emission modes such as the pure ion emission regime, mixed ion-droplet regimes with ions emitted from both droplets in flight and the base of the cone, electrosprays of droplets from liquids with lower electrical conductivities, etc. Each emission mode has characteristic initial conditions and distributions for the charged particles produced by the electrospray (these are also inputs to the model), which in turn depend on the physical properties and the flow rate of the liquid being electrosprayed.

This work was funded by the Jet Propulsion Laboratory, Award No. 1620723 and the Air Force Office of Scientific Research, Award No. FA9550-21-1-0200. We thank the monitors of the program, Dr. John Ziemer and Dr. Colleen Marrese-Reading for their support.

The authors have no conflicts to disclose.

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