A multi-scale approach to electrospray ion source modeling has been developed. The evolution of a single-emitter electrospray plume in a pure ionic regime is simulated with a combination of electrohydrodynamic fluids and n-body particle modeling. Simulations are performed for the ionic liquid, EMI-BF4, firing in a positive pure-ion mode. The metastable nature of ion clusters is captured using an ion fragmentation model informed by molecular dynamics simulations and experimental data. Results are generated for three operating points (120, 324, and 440 nA) and are used to predict performance relevant properties, such as the divergence angle and the extractor surface impingement rate. Comparisons to experimental data recorded at similar operating points are provided.

## NOMENCLATURE

- $a$
tip linear eccentricity

- $a$
particle acceleration vector

- $cp$
liquid heat capacity

- $d$
distance between the tip and the electrode

- $El$
background Laplace electric field

- $Ep$
particle-induced Coulombic (Poisson) electric field

- $e$
elementary charge

- $h$
Planck constant

- $I$
total current emitted

- $j$
current density vector

- $kT$
thermal conductivity

- $mi$
ion mass

- $n$
normal vector to the interface

- $p$
pressure

- $Q$
flow rate

- $qi$
ion charge

- $Rc$
tip curvature radius

- $r$
radial cylindrical coordinate

- $r$
particle position vector

- $rc$
fluid channel radius

- $rh$
extractor aperture radius

- $rsim$
radius of the simulation domain

- $rtip$
tip base radius

- $T$
liquid temperature

- $Te$
electric stress tensor

- $Tf$
fluid stress tensor

- $t$
tangential vector to the interface

- $v$
particle/fluid velocity vector

- $Z$
hydraulic impedance coefficient

- $z$
vertical cylindrical coordinate

- $zel$
extractor thickness

- $zsim$
vertical distance between the tip base and the extractor

- $zt(r)$
tip profile as a function of $r$

- $zup$
vertical extension of the domain above the extractor

### Greek

- $\Delta G$
ion free energy of solvation

- $\eta 0$
non-dimensional tip asymptote parameter

- $\Gamma E$
extractor electrode

- $\Gamma Es$
extractor aperture

- $\Gamma EXT$
domain external boundary

- $\Gamma M$
meniscus interfacial domain

- $\Gamma T$
tip profile boundary

- $\gamma $
surface tension coefficient

- $\kappa (T)$
temperature-dependent electrical conductivity

- $\Omega l$
simulation bulk liquid domain

- $\Omega v$
simulation vacuum domain

- $\mu $
fluid viscosity

- $\varphi $
electric potential

- $\psi $
azimuthal cylindrical coordinate

- $\rho sc$
volumetric charge density in the liquid

- $\sigma $
interfacial charge density

- $\tau $
species mean lifetime

## I. INTRODUCTION

The use of field emission techniques for ion beam formation from a liquid reservoir has wide applications, including lithography,^{1} mass spectrometry,^{2} and high efficiency space propulsion.^{3} For space propulsion, a high velocity ion plume transfers momentum to the vehicle with minimal propellant usage. Electrospray architectures, such as porous needle arrays^{4} and capillary-fed^{5} ionic liquid, and liquid metal^{6} ion sources have all been implemented in flight-tested propulsion systems. In these systems, anywhere from tens to thousands of individual ion emitters are clustered together to generate the momentum required to maneuver the vehicle.

In each case, the system performance depends on the evolution of the individual ion beams from the point of formation at the meniscus-plume (i.e., liquid-vacuum) boundary. The processes of interest for performance span spatial scales ranging from nanometers to centimeters and temporal scales ranging from picoseconds to hours. The prediction of device-level performance and time-integrated phenomena, such as surface interactions, requires starting with a clear picture of the evolution of an individual ion plume from the source.

The evolution of ion plumes in field emission emitters is determined by the structure of the electrostatic fields created on the specific geometry of the emitter. In its most common configuration, a field emission emitter forms a liquid meniscus on top of the needle. When exposed to an external electric field, the meniscus interface adopts a conical shape, generally in equilibrium or in balance between electric, hydrodynamic, and surface tension stresses. Near the apex of the meniscus, the electric fields are strong enough to trigger the field emission process. The position of the apex and initial velocities of the ions are given by the hydrodynamic and stress conditions of the meniscus. The interface of the meniscus departs from the classical Taylor cone structure, generally due to the existence of a finite pressure drop sustaining the flow of ions and the fact that the meniscus is held in a confined tip geometry, where the Taylor solution does not apply. Such modified equilibrium shapes (at least those conserving their axial symmetry) exist under a limited range of conditions, namely, a narrow range of external fields, a minimum hydraulic impedance, and particularly small meniscus sizes on the order of 3–5 $\mu $m.^{7}

This micrometer-scale order of magnitude lies close to the diffraction limit for traditional optical microscopy, which has precluded direct observation of these ionic liquid menisci. Attempts for a direct observation of these menisci have been performed with electron microscopy with limited success since electrons were shown to produce abnormal solidification of the ionic liquid during the observation process.^{8} Consequently, any experimental characterization of these menisci has been relegated to indirect measurements of the total current emitted or plume characterization. The integration of an electrohydrodynamic source model with a discrete plume model presents new opportunities for validation of the meniscus solution. Furthermore, it offers promise to bridge the gap between the macroscale design parameter space (emitter geometry, liquid properties) and the plume evolution via realistic initial conditions in ways that cannot be captured through molecular dynamics or simulations alone,^{9} limited to domain size, or constant axial flux,^{10} which are usually constrained to a limited set of operating conditions.

Electrospray ion sources have been well-characterized experimentally in both single emitter^{11–13} and array configurations.^{4,14} Yet, there is still a gap in fundamental understanding of how the micro-scale processes occurring near the emission region and the electric field structure affect these observations. In particular, information about the dynamics of fragmentation is of particular interest to performance.^{15} As in plasma devices, it is necessary to understand fundamental processes in order to design for maximum lifetime and efficiency.

The characteristics of electrospray plumes (currents per emitter below 1 mA and electrode gaps below 1 mm)^{11,16} make it possible to kinetically track individual particles from the point of emission from the liquid meniscus to beyond the extractor electrode. Recent work has produced equilibrium shapes of the pure ion^{17} and the cone-jet^{18} emission modes. Davis *et al*.^{19} have shown that an n-body approach to tracking electrostatic forces in the plume is feasible for individual emitters. For example, for ionic liquid ion sources emitting currents about $I\u2248300$ nA, a reference for the density of particles $n$ in an acceleration region of $L\u2248100$ $\mu $m length, beam cross section at 50 $\mu $m from the extractor $S\u2248300$ $\mu m2$ (considering a beam of $15\xb0$), potential acceleration at 50 $\mu $m $\varphi \u2248900$ V, and an average mass to charge ratio of singly charged particles of $em=106$ C/kg yields about $n\u2248IeS2em\varphi \u22481017particles/m3$. For a plume simulation volume extending three times the distance between the tip and the extractor $V\u22483LS$, this yields about $N=105$ particles in the simulation volume, which lies on the limit of affordable computational times per time step in direct $n$-body approaches in CPUs.^{20} While the ions can be treated as discrete particles, it has long been known that molecular effects, namely, the fragmentation of metastable ion clusters, have important implications on plume evolution.^{11} Relevant work by Miller and Lozano^{13} and Prince *et al*.^{21} have used experimental and numerical methods to capture the fragmentation rates of solvated ion clusters for various ionic liquids.

The goal of this work is to develop a multiscale model of an individual electrospray ion plume capturing evolution from the sub-micrometer scale emission region to the steady-state region hundreds of micrometers downstream of the extractor. Unique contributions of this approach include (1) a steady-state fluid model for initial conditions, (2) an electric field self-consistent with emission conditions, and (3) a fragmentation model consistent with the accelerating field. This model of the individual emitter plume will serve as a starting point for expanding to array-level investigations to capture processes, such as beamlet mixing and bi-polar ion neutralization.

## II. MODEL DESCRIPTION

### A. Approach

The simulation framework developed through this work spans approximately three orders of magnitude in spatial evolution and contains both continuum fluid and discrete particle components. The emission properties are informed by a continuum electrohydrodynamic (EHD) meniscus model. This model resolves the 2D axially symmetric equilibrium shapes of a meniscus in steady pure-ion evaporation, as a function of the geometry of the emitter and a set of physical parameters of the ionic liquid (density, dielectric constant, surface tension coefficient, conductivity). The model is based on the work by Coffman *et al*.^{22} and Gallud and Lozano^{7} and adapted to handle curved electrode geometries. This adaptation is only geometrical and maintains the physics described in the cited models. The geometrical adaptations extend the emitter electrode from a planar sheet to a needle structure and include an aperture in the extractor electrode to permit the ion flow, as usually encountered in experimental setups. The boundary conditions of the updated geometry are discussed in Secs. II B and II C. The EHD model provides the total current and current density of ions emitted along the meniscus interface in equilibrium, and the structure of the background electric field accelerates the free ions. Details of the EHD model are provided in Sec. II C. Once emitted from the liquid surface, the ions are tracked as discrete particles using an n-body approach to capture inter-particle electrostatic interactions. Details of the n-body model are provided in Secs. III A and III C. This simulation focuses on the pure-ion emission mode of an electrospray source and contains only the three dominant ion types: monomers, dimers, and trimers. The metastable nature of dimer and trimer ion-neutral clusters is captured by a rate model informed by both molecular dynamics simulations and published experimental measurements. Details of the fragmentation model are provided in Sec. III D.

### B. Domain description

The domain considered in this study is shown in Fig. 1. It consists of an axisymmetric hyperboloidal tip ($\Gamma T$) of base radius $rtip=75\mu m$ attached to a flat plate of radius $rsim=300\mu m$ and at a distance $zsim=250\mu m$ of a flat extractor electrode $\Gamma E$ of thickness $zel=30\mu m$. The distance, $d$, between the emission site and the extractor is $100\mu m$ and is treated as vacuum ($\Omega v$). The extractor electrode has an aperture hole of radius $rh=150\mu m$ to permit the flow of particles across it. The simulation domain extends $zup=570\mu m$ beyond the end of the extractor electrode.

The hyperboloidal electrode tip profile [$r,z(r)$] is given by the hyperboloid equation

where $a$ and $\eta 0$ are geometrical parameters. The parameter $a$ is the linear eccentricity of the hyperboloid,

The geometry is specified according to the distance $d$ from the electrode tip to the center of the plane where the extractor plate starts ($\Gamma E$) and the curvature radius $Rc$ at the apex of the tip $\Gamma T$. The parameter $\eta 0$ is a non-dimensional constant governing the asymptote of the tip hyperboloid,

The values of $Rc$ and $d$ used in this simulation are 11 and $100\mu m$, respectively.

Figure 2 shows the vicinity of the tip apex (the asterisk region in Fig. 1). The tip is truncated to accommodate a fluid channel ($\Omega l$) of radius $rc=5\mu m$, where the meniscus experiencing pure-ion evaporation anchors. The meniscus interface is identified as $\Gamma M$. The inlet of the channel $\Gamma I$ is situated at a distance equal to the channel radius $rc$ from the rim of the fluid channel. Similar to Gallud and Lozano^{7} and Coffman *et al*.,^{22} the channel is considered to be very long compared to the meniscus size, and only a small portion of length equal to $rc$ is treated computationally. The axes are non-dimensionalized by the channel length $rc$.

### C. Emission conditions: The EHD model

The EHD model finds the equilibrium shapes of an ionic liquid meniscus experiencing steady pure-ion evaporation under the electrode geometry shown in Fig. 1. Solutions are obtained iteratively by sequentially solving the electrohydrodynamic equilibrium inside and outside the meniscus until residuals are reduced to a negligible value.^{22} In particular, the Poisson equation is solved in $\Omega l$ and the Laplace equation in $\Omega v$ to obtain the Maxwell electric stress distribution along the meniscus interface $\Gamma M$ in the vacuum side ($\tau ev$) and in the liquid side ($\tau el$),

where $\rho sc$ is the space charge in the liquid. The breakup of quasi-neutrality in the ionic liquid is due to the presence of conductivity gradients originated by the heating of the meniscus.^{7} The ion space charge is neglected in the vacuum for ionic liquids ($\rho sc\u22480$ on $\Omega v$) since zeroth order of magnitude analysis shows that it does not modify the electric field that the meniscus observes to a sufficient extent.^{22} In the bulk and meniscus interface, convective charge transport can be neglected as well, and an Ohmic model for the conductivity is considered: $j=\kappa (T)E$. A potential drop $V$ is applied between the tip $\Gamma T$ and the extractor $\Gamma E$,

The normal component of the electric field vanishes on the external boundaries of the simulated domain,

The charge conservation equation is solved in $\Omega l$ to account for the current density distribution $j\u22c5n$ perpendicular to the meniscus interface $\Gamma M$,

where the Ohmic model for the current density is used and $\kappa (T)$ is the conductivity of the ionic liquid, which depends on temperature. The model considers that the evaporated current density follows a phenomenological activated process law,^{23}

where $\sigma $ is the surface charge density, $T$ is the temperature of interface, $kB$ is the Boltzmann constant, $h$ is the Planck constant, $\Delta G$ is the energy of solvation of the ions, $E\u2217\u22484\pi \epsilon 0q3\Delta G2$ is the critical field of emission, $q$ is the charge, and $\epsilon 0$ is the electric permittivity of vacuum.

For the ionic liquids considered in this work, the critical field of emission is about $109$ $Vm$.^{24} This yields an approximate value of $\Delta G\u22481.3$ eV used in the simulation for singly charged ions.

The incompressible Navier–Stokes equations are solved in the liquid domain ($\Omega l$) to obtain the distribution of viscous stresses along $\Gamma M$,

where $\rho $ is the density of the fluid, $v$ is the velocity of the fluid, $Tf=\u2212pI+\mu (\u2207v+\u2207vT)$ is the fluid stress tensor, $p$ is the pressure, and $\mu $ is the viscosity of the fluid. The solver assumes that the flow is fully developed and paraboloidal at the beginning of the computational domain ($\Gamma I$), with a constant pressure $p$ equal to the drop originating from the friction of the upstream flow with the walls. This pressure drop is given by the pressure at the propellant reservoir $pr$ (zero in vacuum conditions) minus the product of a hydraulic impedance coefficient $Z$ with the fluid volumetric flow rate $Q$,

As in most electrokinetic flows, collisions between charge carriers originate oscillations that are dissipated as thermal conduction in the bulk (Joule heating) and dominate the increase of the temperature in the ionic liquid. Parameters, such as viscosity or conductivity, depend strongly on the temperature of the liquid. These changes affect the structure of the equilibrium shape in the meniscus^{7} and, consequently, the initial conditions of the ions. These effects are captured in the energy conservation equation, which provides temperature distributions in $\Omega l$,

where $cp$ is the heat capacity, $\kappa T$ is the thermal conductivity, and $\Phi $ is the viscous dissipation power per unit volume. Equation (13) is subject to an insulating condition along the interface and a fixed temperature of the electrode wall,

The meniscus shape $\Gamma M$ is perturbed along the iterative process until the equilibrium condition is met. The equilibrium condition is such that

where $Tev,Tel$ are the electric stress tensors in the vacuum and liquid, respectively, $Tf$ is the fluid stress tensor, and $\gamma $ is the surface tension coefficient. The superscript $T$ indicates a transpose vector. Initial particle positions, velocities, and probability of emission are given by the equilibrium solutions or when the difference between the electrical and viscous stresses balances the surface tension stress in the normal direction [Eq. (15)] and vanishes in the tangential direction [Eq. (16)] along all the points in the interface $\Gamma M$.

The iterative process starts with a guess of the initial equilibrium shape that is used at a first stage to solve the Poisson equation in the liquid (4), the Laplace equation in the vacuum (5), charge conservation (8), and kinetic law for ion evaporation (9) subject to the boundary conditions in Eqs. (6) and (7). This first stage yields the potential $\varphi $ surface charge distribution in the interface $\sigma $ and in the liquid $\rho sc$. The solution of these variables shows the distribution of electrical stresses along the meniscus interface in the vacuum $Tev$ and liquid $Tel$ together with a distribution of current density along the interface [Eq. (9) can be obtained]. The current emitted can also be obtained from the electric problem by integrating the solution of the current density along the meniscus interface as

The solution of the Navier–Stokes equations is solved after the electric problem [Eqs. (10) and (11)] subject to the outlet condition in Eq. (12) together with mixed Dirichlet and Neumann boundary conditions at the interface. The Dirichlet boundary condition fixes the normal velocity distribution along the interface, which should be consistent with the mass conservation of the current evaporated,

The Neumann boundary condition provides the value of the tangential fluid stress through Eq. (16). The third stage of the solver uses the energy equation (13) to compute the temperature in the meniscus interface, subject to the boundary conditions in Eq. (14) and using the velocity $v$ and current density $j$ postprocessed from the potential solution of the electric problem. After the different physics are resolved (electric, fluid, and thermal subproblems), the meniscus will not generally be in equilibrium. Equation (15) will be fulfilled up to a certain residual value $R$ that has to be very small in the equilibrium condition,

In the last stage of the solver, the residual distribution itself is used to update the equilibrium interface guess for the next iteration until the residual is small enough to consider that an equilibrium shape is found. Details of this surface update mechanism are out of the scope of this paper and can be found in the full description of the model.^{7}

The operational conditions used in the simulations are

Constant . | Value . |
---|---|

1337 V (120 nA) | |

Voltage (V) | 1823 V (324 nA) |

2127 V (440 nA) | |

Hydraulic impedance ($Z$) | $3.66\xd71019Pas/m3$ |

Electrode temperature ($Tw$) | 296 (K) |

Average-charge-to-mass ratio ($q/m$) | $5.9\xd7105$ C/kg |

Reservoir pressure ($pr$) | 0 Pa |

Constant . | Value . |
---|---|

1337 V (120 nA) | |

Voltage (V) | 1823 V (324 nA) |

2127 V (440 nA) | |

Hydraulic impedance ($Z$) | $3.66\xd71019Pas/m3$ |

Electrode temperature ($Tw$) | 296 (K) |

Average-charge-to-mass ratio ($q/m$) | $5.9\xd7105$ C/kg |

Reservoir pressure ($pr$) | 0 Pa |

Once the equilibrium condition is met, the equilibrium fluid velocity fields $v$ and positions $z(r)$ along the interface are used as initial conditions for ions in the n-body trajectory propagator.

Figure 2 shows the equilibrium shapes obtained for the different currents tested in this paper. It can be noticed that higher equilibrium currents correspond to meniscus shapes with a higher curvature. This result is due to the fact that the magnitude of the hydraulic impedance pressure drop is proportional to the current emitted, and it balances the electric stress to a greater extent.

For the simulation results presented in this paper, the equilibrium shapes were computed using the physical parameters similar to EMI-BF4,^{22} and a value of the hydraulic impedance coefficient $Z=4\xd71019$ $Pam3$, which is reference for electrospray tips firing in the pure-ion regime.^{12} The specific values used are as follows:

Constant . | Value . |
---|---|

Electrical conductivity ($\kappa $) | 1.36 S/m |

Viscosity ($\mu $) | 0.038 Pa s |

Surface tension ($\gamma $) | 0.0452 N m |

Density ($\rho $) | 1.24 $g/cm3$ |

Specific heat capacity ($cp$) | 1555.5 J/(kg K) |

Thermal conductivity ($\kappa T$) | 0.195 W/(m K) |

Constant . | Value . |
---|---|

Electrical conductivity ($\kappa $) | 1.36 S/m |

Viscosity ($\mu $) | 0.038 Pa s |

Surface tension ($\gamma $) | 0.0452 N m |

Density ($\rho $) | 1.24 $g/cm3$ |

Specific heat capacity ($cp$) | 1555.5 J/(kg K) |

Thermal conductivity ($\kappa T$) | 0.195 W/(m K) |

## III. SIMULATION MECHANICS

### A. Computational approach

The n-body simulation approach allows for the explicit time integration of individual ion trajectories. The initial conditions for field-evaporated ions are informed by the static solution to the EHD model given in Sec. II C. The process for assigning initial position ($r0$) and velocity ($v0$) values for injected particles is described in Sec. III B. A different EHD solution is used for each operating point (i.e., current or voltage specified). The EHD result provides both the initial conditions for the particles and the background electric field through which they are accelerated.

The trajectory for each particle is integrated in 3D using the following algorithm:

The n-body force calculation step is parallelized using CPU multithreading for computational efficiency. The system specifications of the cluster used and the computing time can be found here. (The simulation was run on a server with 64 Intel Xeon CPU cores of 2.27 GHz, with computations distributed over 32 cores and across 640 total threads. Under these conditions, the simulation time for the 440 nA simulation was 15 h, 2 min of computation excluding the initialization time of the particles.)

### B. Injection

At each time step, a particle is injected at a randomized initial position $(r0,z0,\psi 0)$. The coordinate $r0$ is determined with a differential probability density $dP\u221dj\u22c5ndAe$, where $j\u22c5n$ is the emitted current density at each point along the meniscus surface obtained from Eq. (9) and $dAe$ is the differential area of emission. If the meniscus equilibrium interface is parameterized by $\Gamma M(r,z(r))$, then the differential area of emission yields

where $d\u2113$ is the axisymmetric meniscus interface arc-length differential. Figure 3 shows both the non-dimensional emitted current density $j\u22c5n$ and the normalized probability density for the three currents considered in this paper.

It should be noted that the expected emission probability exactly at the apex is zero due to the occurrence of a singularity since the emission area collapses $(dAe=0)$.

The coordinate $z0$ is chosen in accordance to the meniscus parameterization [$z0=z(r0)$]. The azimuthal coordinate $\psi 0$ is chosen randomly, with uniform probability between $\psi 0\u2208[0,2\pi ]$.

The initial velocity of each ion is informed by the EHD solution. The magnitude is equal to the fluid velocity just inside the meniscus boundary and the vector is aligned with the surface normal,

Thus, ions are emitted from the surface with relatively low initial velocities ($v0<1$ m/s) but are immediately accelerated by the high electric fields immediately outside the meniscus.

### C. Particle trajectory integration

The particle trajectory integrator dynamics are defined by Newton’s second law for a charged particle of mass, $mi$, and charge, $qi$, in the presence of other particles with charges, $qj$,

where the Laplace field $EL(ri)$ is defined as the background Laplace field originated by the potential difference between the electrodes $\Gamma T$ and $\Gamma E$, and the space charge contribution to the electric field, $EC$, is summed over all particles in the domain.

The trajectories of the particles are integrated using a leapfrog scheme, which is second order accurate, symplectic and can be written into a time-centered formulation. The integration is performed in a three step process, and the code employs the “kick-drift-kick” scheme,

where $vin$ and $ain$ are the velocity and acceleration of the particle $i$ at the step $n$ and $\Delta t$ is the time step for the integration. Even for rather large time steps and a long simulation time, the leapfrog scheme maintains qualitatively correct behavior without long-term secular trends seen in other integration methods, such as Runge–Kutta.

The electric fields accelerating the particles in Eq. (22) are computed separately via interpolation of the EHD solution. The EHD model provides the values of the Laplace electric field in an unstructured Delaunay mesh. The Laplace field is interpolated using a Delaunay point search with the barycentric coordinates of the mesh elements.^{25} The space charge field is computed by summing over the interactions between the particles as given by the second term in Eq. (22). The simulation domain is divided into three subdomains, which differ in their electric field characteristics and time step used to iterate the particle motion, as shown in Fig. 4.

Region I is a hemisphere of $5\mu $m starting at the meniscus tip. In this region, electric fields are very strong, especially near the meniscus tip, where the ions are extracted due to field evaporation. The large gradients in field require a sufficiently small time step to ensure the accuracy of the trajectory computation.

The upper bound on the required time step is found heuristically as shown in Fig. 5. The plume is evolved to a steady state producing a Poisson potential distribution. A test particle is situated at a particular position on the meniscus interface, and its trajectory is integrated with multiple time-step options. The meniscus size is very small compared to the scale of the trajectory plotted in Fig. 5; therefore, it is not shown. The relationship between time step and accuracy is shown in the subaxis plot in Fig. 5, where the positions at $z$ are compared with respect to a very small time step of $\Delta t=10\u221216s$. The time steps used in this paper below $\Delta t=1.4x10\u221212s$ are found to resolve the trajectory to within 0.1% accuracy in the region near the meniscus. The convergence plot shows near second-order convergence, as predicted by the leapfrog scheme.

The number of particles that are injected at each time step, $dN$, can, in general, be specified by taking into account the current from the EHD model $I$ and the selected time step,

In all cases analyzed in this work, the time step is set such that exactly one particle is injected per iteration. For the currents analyzed (120–440 nA), the resulting time step is below the $10\u221211$ s cutoff.

In this region, electrostatic forces based on Coulomb’s law are calculated every time step. This region also contains the highest rates of solvated ion cluster fragmentation due to the high background fields. The probability of fragmentation for each solvated particle is updated as a function of the local electric field at each time step. The process for calculating fragmentation probabilities is described in Sec. III D.

In Region II, the same routine for electric field interpolation is used, but the time step for inter-particle force calculations is relaxed. In this region, the space charge contribution to the electric field is updated at a rate that is $1/10$th of the Region I rate. In the leapfrog formulation [Eq. (24)], the acceleration due to the Laplace field is updated at each time step, whereas the acceleration due to Coulombic forces is updated every $10$th iteration. This approach is justified as inter-particle forces have a negligible contribution to the overall electrostatic force ($EC<0.01EL$) beyond a radial distance of $5\mu $m, as will be shown in the results.

In Region III, the particles have experienced $>95%$ of the total potential drop established by the accelerating electrodes. This region marks the simulation boundary for this study and the edge of the region of interest for single-emitter analysis. Beyond this boundary, particles are linearly propagated at constant velocity, fragmentation follows a constant rate law model, and inter-particle Coulombic forces are no longer captured.

### D. Fragmentation

Ionic liquid cluster fragmentation is the process by which the molecules in a cluster dissociate, typically producing a single ion and a neutral cluster. These neutral clusters can collect on the extractor grid or interact with other portions of the ion beam. Fragmentation in the electric field region also affects the energy distribution of the beam, as different parts of the fragmented clusters reach different final kinetic energies.^{13} This can have a detrimental impact on electrospray propulsive efficiency.^{13,26–28} Additionally, these energetic neutrals could generate secondary electrons at detectors, affecting experimental measurements if not properly suppressed.^{29,30} For these reasons, it is critical to include fragmentation in electrospray plume modeling.

Fragmentation occurs both due to the high internal energy of ionic liquid clusters and the force of the electric field pulling the clusters apart.^{10,13} Previous work by Miller and Lozano indicates that the fragmentation process for a population of clusters is an exponential decay process with an electric field-dependent constant rate that occurs on the picosecond to microsecond timescale in a typical electrospray thruster.^{13,31} In the n-body model, ionic liquid cluster fragmentation is modeled as discrete events occurring each time step in the simulation. The probability of a given ion cluster fragmenting during a time step of duration $\delta t$ is given by the integration of the constant rate over the length of the time step,

where $\tau $ is the mean lifetime of the cluster. At the beginning of each time step, the probability of fragmentation is calculated from the mean lifetime of the cluster and the electric field at the location of the cluster. Thus, mean lifetimes for populations of clusters under the same conditions can be used to govern fragmentation behavior of individual clusters via this fragmentation probability. Provided we know the mean lifetime of clusters as a function of electric field, we can perform this process at each time step to model fragmentation with varying timescales over the timescale of the beam simulation. The rest of this section details how the mean lifetimes were determined for each cluster species as a function of the applied electric field. Some error in the fragmentation model is introduced by calculating the probability of fragmentation using the electric field experienced by the cluster at the beginning of each time step. This error results in changes of less than 1% in the final beam energy distribution when the time step is smaller than 2 ps.

Previous work indicates that the mean lifetime of an ion cluster is dependent on the cluster size, the total internal energy of the cluster, and the strength of the applied electric field.^{13,21,32,33} Thus, as the clusters move through the electric field between the emission site and the extractor, the rate at which they fragment, and thus the probability that they fragment, changes. The internal energy of the clusters is assumed to remain constant after emission. Molecular dynamics (MD) results from the previous work provide the dependence of cluster mean lifetime on total internal energy and electric field strength.^{34} Details of the simulation procedures and results are provided in the Appendix. In short, the simulations followed the procedures developed by Coles and Lozano using the Canongia Lopes and Padua force field in LAMMPS.^{27,35–37} This force field is based on quantum mechanical modeling and has been used previously to characterize fragmentation in ionic liquids.^{27,37} While there are limitations to any force parameterization for simulation, a detailed comparison of the available force fields is beyond the scope of this work. $EMI\u2212BF4$ dimers were initialized with random atom positions and velocities, and a Nosé–Hoover thermostat was used to bring the dimer to the desired temperature.^{38} The electric field was applied, and the dimer was simulated until fragmentation occurred. Mean lifetimes were determined by averaging the time it took 5000 independently simulated dimers to fragment for varying initial internal energies and electric field strength conditions. The results of this previous work indicate that the molecular dynamics simulation temperature parameter can be used interchangeably with the total internal energy of the ion cluster.^{34} Fragmentation rates are referenced here using the temperature of the MD simulations for ease of comparison with the previous work.

The n-body simulation accounts for fragmentation in dimers, trimers, and dimers formed from the fragmentation of trimers, also called trimer products. All clusters of one type are assumed to have the same temperature. Further work is required to incorporate nonuniform temperature distributions for each cluster type. Trimers and trimer products are assumed to have the same temperature and electric field dependence as given by the MD results for dimers. Further work is required to determine fragmentation rates for other cluster types using molecular dynamics.

The exact percentage of each cluster species in the beam and the internal energy of each species for a given emitted current level and emitter geometry are not known.^{13} These initial conditions in the beam determine the fragmentation rates of the clusters and the total amount of fragmentation occurring in the emitter. One way to determine these initial conditions is to compare simulation outputs with experimental energy distribution data.^{10,33} The energy distribution of the beam is analyzed using retarding potential analysis (RPA) curves.^{13} The shape of an RPA curve is determined almost entirely by the percentage of each species in the beam and the fragmentation rates of each species as a function of electric field strength.^{13} The effect of spherical spreading on the experimental RPA curves is not accounted for in this work and is the subject of ongoing research.

Iteration was used to determine the percentage of each cluster species in the beam and the internal energy of those species that best matched available experimental data. n-body simulations were repeated using different temperatures for each cluster type and different beam species compositions. For each iteration, the dependence of fragmentation rates on electric field strength was given by the molecular dynamics results for the internal energy specified for that iteration. The fragmentation rates used can be found in the Appendix. Simulated RPA curves were compared to experimental data with the same emitted current to determine appropriate temperatures for each cluster type and the portion of the beam composed of each cluster type.

For each set of possible initial conditions, a single RPA curve generated from 20 ns of simulation was used for comparison to experimental data. This simulation time allowed for trimers, the largest ion species, to reach the end of the electric field region. Ions in the simulation that did not reach the end of the electric field region were not included in the RPA curve. Experimental data are typically collected by averaging beam signals over several minutes of firing. This results from the need to sweep a retarding voltage across a range of 1–2 kV and operate within the response time of the electronics. However, simulated fragmentation processes will come to a steady state far earlier than RPA detection electronics can capture. The fragmentation rates are dominated by the Laplacian electric field contribution. Thus, after all ion species reach the end of the electric field region, the fragmentation processes will have reached their steady-state behavior. Assuming a constant emitter and liquid meniscus geometry, the beam energy distribution will reach a steady state at this time. Moreover, a large number of ions of each species ensure that even low probability fragmentation events are represented in the RPA curve. Thus, it is appropriate to compare the 20 ns simulated RPA curve with experimental data taken over a much longer timescale as has been done in recent work.^{10}

RPA curves are generated by calculating the kinetic energy of each ion cluster when they reach a plane that is designated as the detector. The distance to the detector plane is 8.85 cm, which is the same as the distance to the detector for the experimental curves used for comparison. At the end of the simulation, all ion clusters that traveled through the extractor were propagated to the detector plane assuming negligible external forces. The probability of fragmentation during the time for the clusters to travel to the detector plane was calculated and used to determine how many of the clusters would fragment in this region. Although the electric field extends beyond the extractor, the mean lifetime past the extractor was taken to be the mean lifetime of EMI-BF4 dimers with no electric field or 1.49 $\mu s$.^{13}

Figure 6 shows the mean lifetimes of the different cluster types originally generated by MD simulations that result in the most accurate representation of the experimental RPA curve when used in the n-body simulation. The temperatures of the ion species were 1000, 1350, and 810 K for dimers, trimers, and trimer products, respectively. This simulation was composed of 40% monomers, 40% dimers, and 20% trimers. Figure 7 shows the similarity between the simulated RPA curve and the experimental RPA curve at 323 and 327 nA of emitted current, respectively.

There are two notable differences between the simulated and experimental RPA curves that can be explained by the fragmentation and RPA curve implementations. First, there is a discrepancy at the lowest voltages. This region of the RPA curve corresponds to very low energy particles that likely result from multiple fragmentations of clusters larger than trimers, which are not included in the n-body simulation. Second, the experimental curve vertical steps are rounded more than the numerical curve. This is likely a result of the energy spreading that results from the source not being located exactly at the center of curvature of the spherical RPA detector.^{13} Overall, the fragmentation model matches well with experimental data and reproduces a simulated energy spectrum close to the empirical results.

## IV. SIMULATION CASES

Multiple case studies were performed using the n-body code for a range of stable menisci configurations from the EHD model.

All of the simulations utilize the ionic liquid 1-ethyl-3-methyl-imidazolium-tetrafluoroborate (EMI-BF4). The physical coefficient values used in this simulation correspond to the parameters of EMI-BF4.^{22} Molecules emitted are either monomers, dimers, or trimers. Higher order particles and droplets are excluded due to their very high probability of breaking up directly after emission and their negligible contribution to measured mass spectra when operating in the pure-ion emission mode.^{11} All neutrals included in the simulation are a result of the fragmentation of solvated ion clusters.

The three cases analyzed are a low current stable meniscus emitting 120 nA, a medium current stable meniscus emitting 324 nA, and a high current meniscus emitting 440 nA, which operates near the static stability boundary. The EHD results are generated for a hydraulic impedance, which is representative of electrospray emitters firing in the pure-ion regime. Simulations are run until the ion plume reaches a steady state in the simulation domain according to the parameters presented below.

The parameters for the simulations are summarized in Table I. The number of steps is equal to the number of injected particles since the code is set up to inject one particle each time step ($dN=1$). The number of particles in the table is all of the particles injected into the domain over the full simulation time.

Current (nA) . | Voltage (V) . | Time (ns) . | Time step (ps) . | Number of particles . |
---|---|---|---|---|

20 776 | ||||

120 | 1337 | 20 | 1.34 | Charged: 14 949 |

Neutrals: 5827 | ||||

56 261 | ||||

324 | 1823 | 20 | 0.50 | Charged: 40 392 |

Neutrals: 15 869 | ||||

59 690 | ||||

440 | 2127 | 16^{a} | 0.37 | Charged: 43 857 |

Neutrals: 15 833 |

Current (nA) . | Voltage (V) . | Time (ns) . | Time step (ps) . | Number of particles . |
---|---|---|---|---|

20 776 | ||||

120 | 1337 | 20 | 1.34 | Charged: 14 949 |

Neutrals: 5827 | ||||

56 261 | ||||

324 | 1823 | 20 | 0.50 | Charged: 40 392 |

Neutrals: 15 869 | ||||

59 690 | ||||

440 | 2127 | 16^{a} | 0.37 | Charged: 43 857 |

Neutrals: 15 833 |

^{a}

Shorter simulation time than the other two cases due to computational burden at a smaller time step.

Figures 8 displays the plume parameters at the final iteration of the simulation for the 440 nA condition. Monomers are indicated by blue markers, dimers by red, trimers by green, and neutrals by yellow. The coordinate $z=0$ corresponds to the tip electrode apex. The tip profile is not shown in this figure and would elongate through $z<0$. The extractor electrode profile ($\Gamma E$ in Fig. 1) is shown in gray. The extractor has a circular aperture of $150\mu $m radius and a thickness of $30\mu $m.

The final state for each simulated operating point reveals a relatively well-defined ion plume surrounded by a more widely diverging neutral population yet to fully evolve on the ion timescale. While not captured in the simulation time, it is observed that the lowest energy neutrals on the periphery of the plume have trajectories that will lead to electrode impact at longer timescales. The displayed results are analyzed quantitatively in Sec. V to evaluate the divergence angle, the fragmentation fraction, and the neutral flux.

The fragmentation fraction is displayed in each figure next to the total number of neutrals at the end of simulation. It is used as a metric to evaluate the validity of the simulations with respect to observed experimental conditions. The fragmentation fraction is defined as

using the number of dimers, trimers, and neutrals at the last time step.

The fragmentation fraction is reported for each simulation and amounts to between 60% and 65% for all three cases, which is in close agreement with experimental measurements for similar operating points.^{13}

### A. Energy conservation

As the n-body approach only considers energy-conserving Coulomb collisions, comparing the change in potential and kinetic energy can be used as a tool to assess the accuracy of this the model. In this case, the change in kinetic energy of the particles beyond their injection energy should be approximately equal to the electrostatic potential difference traversed. To assess conservation of energy, the change in kinetic energy and the change in electric potential are calculated for each unfragmented ion that traverses the $250\mu $m simulation domain in the simulation time. The leapfrog integration scheme has second-order accuracy; thus, the increase in the kinetic energy of the particles overall should be within 1% of the change in potential energy. Coulomb collisions are energy conserving and thus would result in an energy spread but not a net increase or decrease in kinetic energy when summed over the population of particles.

Figure 9 shows a histogram of the ratio of the change in kinetic energy to the change in electrostatic potential between the injection values and the values at the final time step. Results are shown for each simulation case. The results indicate that energy conservation is satisfied to within 0.1–0.5% for these simulations. These results also confirm that space charge effects, particularly the Boersch effect, which would be captured in the energy spread, are minimal.

## V. PLUME ANALYSIS AND EXPERIMENTAL COMPARISON

### A. Divergence angle

The divergence angle is calculated by integrating over the solid angle ($0.5\xb0$ steps) and solving for the angle at which 80% and 99% of all particles are contained. This calculation is performed at both the extractor location and the simulation boundary of $250\mu $m. The calculated angle is equivalent at both locations; thus, only the $250\mu $m result is reported. The current density is determined from the number of particles, all singly charged, crossing the boundary per unit area and time. To ensure a statistically relevant number of particles pass through the hemisphere, the divergence angle is calculated by averaging over the last 2–4 ns when the plume has reached a steady state. Results are shown for all currents in Fig. 10. The divergence angle shown represents the half-angle as azimuthal symmetry is assumed for the plume. Two discretizations for the angles are used in the calculations: $3\xb0$ for the plot and $0.5\xb0$ for the values reported in Table II. This achieves a higher resolution for the numerical value and less noise in the plot.

Current (nA) . | 80% (250 μm)
. | 99% (250 μm)
. |
---|---|---|

120 | $7.0\xb0$ | $12.0\xb0$ |

324 | $8.0\xb0$ | $12.0\xb0$ |

440 | $8.0\xb0$ | $12.0\xb0$ |

Current (nA) . | 80% (250 μm)
. | 99% (250 μm)
. |
---|---|---|

120 | $7.0\xb0$ | $12.0\xb0$ |

324 | $8.0\xb0$ | $12.0\xb0$ |

440 | $8.0\xb0$ | $12.0\xb0$ |

The angles containing 80% and 99% of the charged species are presented in Table II. The resulting divergence half-angles in the simulated plumes exhibit only little dependence on the operating current and are within 7–12$\xb0$ for all cases considered.

These simulations also predict the evolution of a band of neutral molecules produced via fragmentation around the ion plume. The divergence angle of the neutral population is shown in Table III. For the neutral particles, an increase in current leads to an increase in the divergence. Unlike the ions, there is a significant difference in the angle, which contains 80% and 99% of the neutral particles. This indicates that they follow a different angular distribution profile and their spreading is amplified at higher currents.

Current (nA) . | 80% (250 μm)
. | 99% (250 μm)
. |
---|---|---|

120 | $9.5\xb0$ | $22.5\xb0$ |

324 | $10.0\xb0$ | $27.5\xb0$ |

440 | $10.0\xb0$ | $29.0\xb0$ |

Current (nA) . | 80% (250 μm)
. | 99% (250 μm)
. |
---|---|---|

120 | $9.5\xb0$ | $22.5\xb0$ |

324 | $10.0\xb0$ | $27.5\xb0$ |

440 | $10.0\xb0$ | $29.0\xb0$ |

The profile of an electrospray beam has been experimentally observed to be something between parabolic and super-Gaussian. The super-Gaussian shape is given by

where $\theta $ is the dependent variable, $A=1$ is the magnitude due to normalization, $\theta t=0$ the tilt in the angle as we assume perfect on axis firing, $\theta 0=8.812\xb0$$[8.742,8.881]$, and $n=2.342$$[2.215,2.469]$ the sharpness. The model was proposed by Thuppul *et al.*,^{39} and the values for $n$ and $\theta 0$ are determined by using a least square method curve fit, where the values in the square parentheses are the 95% confidence interval values of the curve fit. The same procedure is used to fit the experimental data with the same model yielding an estimate of $\theta 0=13.18\xb0$$[12.46,13.89]$ and $n=2.3$$[1.517,3.083]$.

Overall, the simulation results for ion beam divergence are in close agreement with the experimental results in shape as can be seen by the curve fit parameters. The difference in the angle may be due to multiple factors, which are discussed below.

The divergence angle is tightly related to the radius of curvature $rt$ of the ion trajectory, which can be obtained with an expression for the centripetal force $miac=qE\u22c5n=miv2rt$,

where $vi0$ and $\varphi 0$ are the initial velocity and potential energy of the ion at the beginning of flight on the meniscus interface.

It is hard to predict general trends on the divergence angle of the trajectories in electrospraying since the elements of Eq. (30), namely, $vi0,\varphi 0,\varphi (r),E(r)$, are itself a result of the equilibrium solution of the emitting meniscus and the electric field structure affected by the geometry of the electrodes. In the assumption that most of the particles arrive at a straight trajectory ($rt\u2192\u221e$) very close to the meniscus apex, some qualitative observations can be made. For example, it can be deduced that in the absence of fragmentation, the trajectories of the ions will mostly depend on the geometrical structure of the field regardless of their charge-to-mass ratio if neglecting the initial velocity term in Eq. (30). However, when dimers and trimer species fragment under acceleration by the (mostly axial) electric field, the radius of curvature of their trajectory will decrease from their unfragmented counterparts due to the increase in the charge-to-mass ratio after fragmentation. Thus, the relationship between fragmentation and the electric field structure is critical to accurately predicting beam divergence.

The structure of the electric field depends on aspects, such as the geometrical configuration of the electrodes and the meniscus interface or the presence of the space charge. In this regard, a tip located closer to the extracting electrode could enhance the intensity of the radial electric field compared to the axial field and intensify the spreading effect on the divergence angles. This could explain the differences of the simulated ion divergence angle with the experimental data presented in Fig. 10, where the emitter is placed flush with the extractor grid and has an aperture of $1400\mu $m instead of $300\mu $m.

The shape of the meniscus near the emission region could also have an important impact on the divergence angle, especially when menisci exhibit cusped equilibrium shapes that strengthen the radial component of the electric field near the emission region. In the literature of liquid metal ion sources, it has been observed that these cusped equilibrium shapes are favored by shorter emission region length scales $r\u2217$ and higher currents. The emission region length-scale $r\u2217$ is itself a by-product of solving the electrohydrodynamic equations presented in Sec. II C but can be approximated to $r\u2217\u2248\gamma E\u22172$ by using a rough order of magnitude approximation of Eq. (15) if considering an emission region similar to a spherical cap. Consequently, higher critical fields and smaller surface tension coefficients will enhance the cusped shape of the meniscus near the emission region and possibly incur in higher ion divergence angles. It is worth mentioning that no precise data of $E\u2217$ exist for ionic liquids. Simulation tools, such as the one presented in this paper, could enable physics-based correlations of $E\u2217$ to metrics, such as the divergence angle, which is an experimentally observable quantity.

In addition to increasing the radial component of the electric field, space charge effects are shown to elongate these cusps.^{40} Interestingly, space charge effects are strongly related to the lower charge-to-mass ratio of the ions $qi/mi$. While no direct effects on the ion dynamics are expected from the variability of $qi/mi$ except those derived from fragmentation, lower $qi/mi$ ions have longer residence time near the meniscus interface and, therefore, more screening effects of the field that the meniscus sees. The effect of this additional screening is an elongated cusped interface, which, as mentioned before, favors the strength of the radial field compared to the axial field. The dependence of the divergence angle with $qi/mi$ through an enhanced meniscus cusp length is a known experimental observation for liquid metal ion sources.^{41}

While it is beyond the scope of this study to provide a definitive explanation for the discrepancy, the integrated fluid-discrete particle simulation framework will enable future explorations of this parameter space. The ability to resolve the physics on length and temporal scales spanning emission dynamics to the far field structure of the acceleration region is required. A contribution of this preliminary work is showing the relative unimportance of space charge in beam spreading for currents in the range of 100–500 nA. This result leaves the electric field structure, the meniscus geometry (including its coupled effect with the plumes), and molecular interactions near the emission site (see Sec. V D) as factors to be explored in future work.

### B. Laplace vs Poisson field

The relative importance of the Poisson space charge effects vs the Laplace extracting field depends on the concentration of ions at a particular point in the domain. These effects tend to be more significant near the vicinity of the emission region, where the plume has not expanded yet and the ion density is maximum. For reference, particle densities can be seen in Fig. 18 as a function of the axial coordinates for the maximum current simulated. Still, the characteristic low currents of pure-ion evaporation in ionic liquids originate space charge concentrations of such a limited strength that the Laplace field contributes $85%$ or more of the total field near the emission region. This can be seen on the right side of Figs. 11–13. They show $|EC||EC|+|EL|$, the axisymmetric average of the magnitude of the Poisson electric field over the sum of both Poisson and Laplace field magnitudes near the emission region. This average fraction expands with the current density, as expected by the increase of ion concentration.

The diagrams shown capture the last simulation step to ensure the steady state of the plume. On the left side, the effect of this balance between $EL$ and $EC$ on the ion trajectories is plotted. The solid lines show the trajectories of a monomer that experiences the averaged axisymmetric Laplace field only, and the dashed lines include the effect of space charge. Figures 11–13 show the emission region of the meniscus (about 250 nm width), with its apex drawn in bold black. It can be observed that the space charge effects on the trajectories are most significant at the highest currents and for particles starting at regions slightly displaced from the apex of the meniscus tip. This result may be an effect of the increased particle concentration in these regions. At such a close distance to the emission region, the particle distribution still resembles the probability density function shown in Fig. 3.

### C. Fragmentation and neutral fluxes

The incorporation of fragmentation modeling reveals the importance of the neutral particle population. As the fragmentation of solvated clusters occurs across the domain, neutral pairs are created with a wide range of energies and trajectories.

Figure 14 shows the distribution of fragmentation rates across the simulation domain, which represents the convolution of the probability of fragmentation and residence time of ions in regions of varying field strength. The largest electric fields (beyond $108$ V/m) exist within a few micrometers of the emission region. The fragmentation events influencing the energy distribution of ions occur in the acceleration region, defined as Regions I and II in the domain. While fragmentation continues to occur in the field free region (up to 50% of the events, according to experimental findings^{13}), these events are unimportant to the evolution of the beam in the present domain volume.

The fragmentation rate plot reveals populations of neutrals with vastly different characteristics. The neutrals created near the emission site have low energies and wide divergence angles. This population of neutrals is expected to impinge on the bottom of the extractor grid within the microsecond timescale. While beyond the steady-state timescale of the ions in the domain, the neutral flux to the extractor electrode can be estimated from the 16–20 ns results assuming ballistic trajectories.

The neutral flux is estimated by linearly propagating the neutrals, which originated from fragmentation processes during a 16–20 ns simulation until they reach the extractor surface and averaging over the simulation time. To predict flux gradients, the surface is divided into $3\mu $m square bins. The number of neutrals that fall within each bin is summed and then spline interpolated to produce the impingement map shown in Fig. 15. Azimuthally averaged flux and energy predictions are shown in Figs. 16 and 17, respectively. Figure 17 shows both the median and max energy of impinging neutrals.

The result of impingement depends on the energies of the neutrals, as they may deposit on the surface intact, reflect or break apart upon impact, or remove material from the surface. Neutrals originating from fragmentation processes close to the emitter tip have much lower kinetic energies than those which fragmented further into the acceleration region. Figure 17 shows that peak neutral energies are expected to be near the extraction potential at the middle of the aperture and quickly decay to 20 eV or less within the first $50\mu $m radial band for all currents studied. The average energy of impinging neutrals is below 10 eV for radial distances of $150\mu $m or greater where the extractor electrode surface is defined. The noise in predicted values for large radii is likely due to the relatively small number of particles generated in the 16–20 ns simulation time. Further investigations using molecular dynamics and interfacial physics are required to predict the resulting surface evolution.

### D. Collisions

The only particle interaction inherently captured in this model is Coulombic repulsion between charged species. In this case, the neutrals formed from fragmentation events do not interact with other particles in the plume due to their net zero charge. A detailed model of ion-neutral interactions and the timescales required to accurately capture them are beyond the scope of this study. Furthermore, the simulation times are not long enough for slow-moving neutral fragments to fill the simulation domain. The longest simulation time in this work was 20 ns, which sets a lower bound for all detectable collisions at 50 MHz.

Despite this limitation, the initial plume evolution can be used to estimate collision frequencies between different particles and inform longer time-scale investigations. Analysis was performed in post-processing to assess the rates at which various particle pairs might come within close enough contact to interact at the molecular scale. First, the densities and velocities of each species are averaged across the cross section of the beam as a function of axial distance from the emission site. This information is then used to make estimates of axially varying collision rates, which are beyond the limit of detection of the simulation.

The collision rates are estimated by

where $n1$ and $n2$ are the densities of the particles (Fig. 18), $\sigma $ is the collisional cross section, $v1$ and $v2$ are the average particle velocities (Fig. 19), and $vol$ is the volume of the cross-sectional slice of the conically bounded plume at a given axial location. A collision cross section of $3\xd710\u221218m2$ is used for all calculations as an estimate for the boundary at which solvated clusters definitively attach or detach based on previous molecular dynamics studies.^{27}

Estimated collision rates for monomers–trimers, monomers–dimers, and monomers–neutrals were calculated from the highest current simulation as shown in Fig. 20. The monomer–neutral collisions are only evaluated up to $10\mu $m from emission as the slow-moving neutral population has only fully propagated this far. Results are shown for other particle pairs up to an axial distance of $25\mu $m to capture the region of highest collisionality. A comparison of these projected interaction rates (>MHz) with the required operating time for the ion source (seconds to hours) indicates that these types of interactions will be important for long-term plume evolution.

## VI. CONCLUSIONS AND FUTURE WORK

This work presents the framework and results for multiscale, multiphysics modeling of ion emission from an electrically stressed liquid meniscus. The steady-state emission conditions, including the meniscus topology, are informed by an electrohydrodynamic (EHD) model of pure-ion emission. The EHD solution determines the emitted current, the ion initial conditions, and the Laplacian electric field solution for the domain. The evolution of the ion plume is captured using an n-body model for interparticle electrostatic forces. Fragmentation effects are taken into account based on an experimentally and numerically informed model with dependence on the surrounding electric field.

While the simulation times are limited to 20 ns or less in this work, results can be used to predict longer timescale phenomena as the ion species reach a steady-state condition at this time. The steady-state energy spectra predicted by this model accurately reproduce experimentally measured results. Steady-state density and velocity profiles for various species are used to predict collision rates, indicating that simulations at the microsecond scale or longer will be needed to capture collisional effects. The products of the predicted ion-neutral and ion-ion collisions should be further investigated using molecular dynamics and fed back into the simulation to predict effects on plume structure and composition. Similarly, the flux of low energy neutrals to the periphery of the extractor electrode is expected to reach a steady-state $\u22481\mu $s after the onset of emission. For the geometry modeled in this work, neutral fragments impinging on the grid have energies below 10 eV or approximately 0.5% of the applied potential.

While the simulation times are limited to 20 ns or less in this work, results can be used to predict longer timescale phenomena as the ion species reach a steady-state condition in this time. The steady-state energy spectra predicted by this model accurately reproduces experimentally measured results. Steady-state density and velocity profiles for various species are used to predict collision rates, indicating that simulations at the microsecond scale or longer will be needed to capture collisional effects. Similarly, the flux of low energy neutrals to the periphery of the extractor electrode is expected to reach steady-state approximately $1\mu $s after the onset of emission. For the geometry modeled in this work, neutral fragments impinging on the grid have energies below 10 eV or approximately 0.5% of the applied potential.

The multiscale model developed here captures the properties and evolution of an ionic liquid ion beam operating in a pure-ion emission mode. The resulting plume properties can be used as the input for simulations extended in both spatial and temporal dimensions. Examples of these applications include the interactions of beams in multiplexed emitter arrays and molecular phenomena, such as surface interactions and inter-particle collisions.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge support from the NASA Early Stage Innovations (Grant No. 80NSSC19K0211) and the NY Space Grant Fellowship. E. M. Petro would like to thank Summer Hoss for contributing to the energy and collisionality analysis. M. Schroeder would like to thank Kyle Sonandres and David Hernandez for their contribution to the fragmentation simulations. X. Gallud would like to thank Steven Liu and Aditya Mehrotra for their contributions improving the computational efficiency of the code. S. K. Hampl would like to thank the German Academic Scholarship Foundation for supporting his research stay at MIT. In addition, P. C. Lozano would like to thank the M. Alemán-Velasco Foundation for its support.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

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

### APPENDIX: MOLECULAR DYNAMICS DETAILS

Molecular dynamics simulations were performed to determine the fragmentation rates of ionic liquid clusters in an electric field. 5000 dimers were simulated individually at different internal energies and electric field conditions, and the mean lifetime was calculated for each of these conditions from the results. This appendix details the simulation process and the key results used to determine the fragmentation rates that best matched experimental data in the n-body model.

The simulation process is based on the previous work by Coles *et al.* and Prince *et al.*^{21,26,27,33} The simulations used LAMMPS, an open-source MD software that has been used previously for ionic liquid simulations.^{10,26,27,42,43} The simulations use the Lopes and Padua force field.^{35} The time step for the simulations was 0.5 fs. The procedure for simulating a single cluster was as follows:

Determine coordinates of atoms in the cluster from stable emitted clusters. These data were taken from emission simulations performed by Coles

*et al*.^{26,27}Randomize the velocities of each atom of the cluster at 298 K.

Apply a Nosé–Hoover thermostat to equilibrate the temperature of the cluster to the desired temperature.

^{38}The damping parameter of the thermostat was 1000 fs.Apply a constant electric field strength until the separation between the molecules of the cluster reaches the fragmentation threshold $\delta f(sim)$. More about $\delta f(sim)$ is described below.

The goal of the thermostat in the second step is to set the internal energy of the cluster. Unfortunately, current MD methods do not provide a simple way to set internal energy; therefore, the temperature of the cluster was used instead. However, results show that the standard deviation in initial internal energy for each group of samples was much lower than the difference in energy between the simulated temperatures. Fragmentation rates are referenced here using the temperature of the MD simulations.

Ionic liquid clusters simulated with electric fields often develop large separations between the molecules in the cluster without fragmenting.^{26} In order to determine when fragmentation occurs, $\delta f(sim)$ must be larger than the largest maximum separation that could occur without fragmentation. A $\delta f(sim)$ value of 40 Å was selected in agreement with the previous work.^{10,27} Each cluster was simulated until this value of $\delta f(sim)$ was reached. While a separation of 40 Å is appropriate for determining whether a cluster has fragmented, the time at which that separation is reached does not necessarily represent the time at which the cluster began to fragment. Particularly for the high electric field and high temperature cases, the time it takes for the maximum separation to reach 40 Åis large compared to the time it takes for the cluster to begin to fragment. To account for this, the time at fragmentation was determined during post-processing using the value $\delta f(post)$. The final chosen $\delta f(post)$ value for dimers was 20 Å, which is in agreement with the recent work.^{10} The mean lifetime, $\tau $, was determined for a set of samples simulated at the same conditions by averaging the lifetimes of all of the samples. Table IV shows the mean lifetime results for positive EMI-$BF4$ dimers.

E (V/A) | 0.8 | 0.3 | 0.15 | 0.1 | 0.05 | 0.02 | 0.01 | 0.005 | 0.0005 | 0.00005 |

T (K) | ||||||||||

300 | 0.582 | 4.739 | … | … | … | … | … | … | … | … |

600 | 0.571 | 2.644 | 301.638 | … | … | … | … | … | … | … |

1000 | 0.555 | 1.665 | 20.867 | 87.001 | … | … | … | … | … | … |

1500 | 0.538 | 1.289 | 5.950 | 15.521 | 54.426 | 162.737 | 244.584 | 333.153 | 433.944 | … |

2000 | 0.518 | 1.103 | 3.304 | 6.414 | 16.050 | 33.355 | 41.604 | 48.082 | 56.912 | 56.983 |

2500 | 0.498 | 0.986 | 2.234 | 3.202 | 7.969 | 13.663 | 16.540 | 18.992 | 20.923 | 21.236 |

3000 | 0.477 | 0.890 | 1.776 | 2.598 | 4.881 | 7.719 | 9.084 | 10.208 | 11.105 | 11.120 |

E (V/A) | 0.8 | 0.3 | 0.15 | 0.1 | 0.05 | 0.02 | 0.01 | 0.005 | 0.0005 | 0.00005 |

T (K) | ||||||||||

300 | 0.582 | 4.739 | … | … | … | … | … | … | … | … |

600 | 0.571 | 2.644 | 301.638 | … | … | … | … | … | … | … |

1000 | 0.555 | 1.665 | 20.867 | 87.001 | … | … | … | … | … | … |

1500 | 0.538 | 1.289 | 5.950 | 15.521 | 54.426 | 162.737 | 244.584 | 333.153 | 433.944 | … |

2000 | 0.518 | 1.103 | 3.304 | 6.414 | 16.050 | 33.355 | 41.604 | 48.082 | 56.912 | 56.983 |

2500 | 0.498 | 0.986 | 2.234 | 3.202 | 7.969 | 13.663 | 16.540 | 18.992 | 20.923 | 21.236 |

3000 | 0.477 | 0.890 | 1.776 | 2.598 | 4.881 | 7.719 | 9.084 | 10.208 | 11.105 | 11.120 |

## REFERENCES

*Proceedings of the 2001 Particle Accelerator Conference*(IEEE, 2001), Vol. 1, pp. 76–80.

*AIAA Propulsion and Energy 2020 Forum*(AIAA, 2020), p. 3612.

*AIAA Propulsion and Energy 2020 Forum*(AIAA, 2020), p. 3609.

*2019 International Electric Propulsion Conference*(IEPC, 2019).

*Proceedings of the 36th International Electric Propulsion Conference, IEPC-2019, Vienna, Austria*(IEPC, 2019), Vol. 590.

*GPU Computing Gems Emerald Edition*(Morgan Kaufmann, 2011), pp. 113–132.

*AIAA Propulsion and Energy Forum and Exposition 2019*(AIAA, 2019).

*48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit 2012*(AIAA, 2012).

*49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference*(AIAA, 2013), Vol. 1, Part F.

*The 36th International Electric Propulsion Conference*(IEPC, 2019), pp. 1–10.

*53rd AIAA/SAE/ASEE Joint Propulsion Conference 2017*(AIAA, 2017), pp. 1–17.

*44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit*(AIAA, 2008).

*45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit*(AIAA, 2009).