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.

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

ΔG

ion free energy of solvation

η0

non-dimensional tip asymptote parameter

ΓE

extractor electrode

ΓEs

extractor aperture

ΓEXT

domain external boundary

ΓM

meniscus interfacial domain

ΓT

tip profile boundary

γ

surface tension coefficient

κ(T)

temperature-dependent electrical conductivity

Ωl

simulation bulk liquid domain

Ωv

simulation vacuum domain

μ

fluid viscosity

ϕ

electric potential

ψ

azimuthal cylindrical coordinate

ρsc

volumetric charge density in the liquid

σ

interfacial charge density

τ

species mean lifetime

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 arrays4 and capillary-fed5 ionic liquid, and liquid metal6 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 μ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 emitter11–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 ion17 and the cone-jet18 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 I300 nA, a reference for the density of particles n in an acceleration region of L100μm length, beam cross section at 50 μm from the extractor S300μm2 (considering a beam of 15°), potential acceleration at 50 μm ϕ900 V, and an average mass to charge ratio of singly charged particles of em=106 C/kg yields about nIeS2emϕ1017particles/m3. For a plume simulation volume extending three times the distance between the tip and the extractor V3LS, 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 Lozano13 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.

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 Lozano7 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.

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

FIG. 1.

Computational domain.

FIG. 1.

Computational domain.

Close modal

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

z(r)=η0a24+r21η02,
(1)

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

a=2d1+Rcd.
(2)

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

η0=(1+Rcd)12.
(3)

The values of Rc and d used in this simulation are 11 and 100μ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 (Ωl) of radius rc=5μm, where the meniscus experiencing pure-ion evaporation anchors. The meniscus interface is identified as ΓM. The inlet of the channel ΓI is situated at a distance equal to the channel radius rc from the rim of the fluid channel. Similar to Gallud and Lozano7 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.

FIG. 2.

Equilibrium shapes of the ionic liquid meniscus for the simulated currents. The liquid domain is referenced as Ωl, the meniscus interface as ΓM, the inlet boundary as ΓI, and the wetted part of the tip boundary as ΓT. Length scales are non-dimensionalized by the channel radius rc. The top left subaxis contains a detailed structure of the meniscus interface near the emission region, where a diagram of the current density profile distribution normal to the interface is shown. The current density data correspond to the I = 120 nA case and is not drawn up to scale. Numerical data of the current densities for the three currents simulated can be found in Fig. 3.

FIG. 2.

Equilibrium shapes of the ionic liquid meniscus for the simulated currents. The liquid domain is referenced as Ωl, the meniscus interface as ΓM, the inlet boundary as ΓI, and the wetted part of the tip boundary as ΓT. Length scales are non-dimensionalized by the channel radius rc. The top left subaxis contains a detailed structure of the meniscus interface near the emission region, where a diagram of the current density profile distribution normal to the interface is shown. The current density data correspond to the I = 120 nA case and is not drawn up to scale. Numerical data of the current densities for the three currents simulated can be found in Fig. 3.

Close modal

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 Ωl and the Laplace equation in Ωv to obtain the Maxwell electric stress distribution along the meniscus interface ΓM in the vacuum side (τev) and in the liquid side (τel),

2ϕ=ρscε0inΩl,
(4)
2ϕ=0inΩv,
(5)

where ρ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 (ρsc0 on Ω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=κ(T)E. A potential drop V is applied between the tip ΓT and the extractor ΓE,

ϕ=VonΓTΓI,ϕ=0onΓE.
(6)

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

ϕn=En=0onΓEXT.
(7)

The charge conservation equation is solved in Ωl to account for the current density distribution jn perpendicular to the meniscus interface ΓM,

j=(κE)=0inΩl,
(8)

where the Ohmic model for the current density is used and κ(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 

jn=σkBThexp(ΔGEnEkBT),
(9)

where σ is the surface charge density, T is the temperature of interface, kB is the Boltzmann constant, h is the Planck constant, ΔG is the energy of solvation of the ions, E4πε0q3ΔG2 is the critical field of emission, q is the charge, and ε0 is the electric permittivity of vacuum.

For the ionic liquids considered in this work, the critical field of emission is about 109Vm.24 This yields an approximate value of ΔG1.3 eV used in the simulation for singly charged ions.

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

ρ(v)v=T+ρscEonΩl,
(10)
v=onΩl,
(11)

where ρ is the density of the fluid, v is the velocity of the fluid, Tf=pI+μ(v+vT) is the fluid stress tensor, p is the pressure, and μ is the viscosity of the fluid. The solver assumes that the flow is fully developed and paraboloidal at the beginning of the computational domain (Γ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,

p=prQZ=prIρqmZ.
(12)

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 meniscus7 and, consequently, the initial conditions of the ions. These effects are captured in the energy conservation equation, which provides temperature distributions in Ωl,

ρcpvT=κT2T+jjκ(T)+ΦinΩl,
(13)

where cp is the heat capacity, κT is the thermal conductivity, and Φ 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,

Tn=0onΓM,T=TwonΓTΓI.
(14)

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

nT(TevTelTf)n=γnonΓM,
(15)
tT(TevTelTf)n=0onΓM,
(16)

where Tev,Tel are the electric stress tensors in the vacuum and liquid, respectively, Tf is the fluid stress tensor, and γ 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 Γ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 ϕ surface charge distribution in the interface σ and in the liquid ρ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

I=ΓMjndΓM.
(17)

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,

vn=qmjnonΓM.
(18)

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,

R=nT(TevTelTf)nγnonΓM.
(19)

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

ConstantValue
 1337 V (120 nA) 
Voltage (V) 1823 V (324 nA) 
 2127 V (440 nA) 
Hydraulic impedance (Z3.66×1019Pas/m3 
Electrode temperature (Tw296 (K) 
Average-charge-to-mass ratio (q/m5.9×105 C/kg 
Reservoir pressure (pr0 Pa 
ConstantValue
 1337 V (120 nA) 
Voltage (V) 1823 V (324 nA) 
 2127 V (440 nA) 
Hydraulic impedance (Z3.66×1019Pas/m3 
Electrode temperature (Tw296 (K) 
Average-charge-to-mass ratio (q/m5.9×105 C/kg 
Reservoir pressure (pr0 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×1019Pam3, which is reference for electrospray tips firing in the pure-ion regime.12 The specific values used are as follows:

ConstantValue
Electrical conductivity (κ1.36 S/m 
Viscosity (μ0.038 Pa s 
Surface tension (γ0.0452 N m 
Density (ρ1.24 g/cm3 
Specific heat capacity (cp1555.5 J/(kg K) 
Thermal conductivity (κT0.195 W/(m K) 
ConstantValue
Electrical conductivity (κ1.36 S/m 
Viscosity (μ0.038 Pa s 
Surface tension (γ0.0452 N m 
Density (ρ1.24 g/cm3 
Specific heat capacity (cp1555.5 J/(kg K) 
Thermal conductivity (κT0.195 W/(m K) 

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:

  1. Interpolate to find a local electric field.

  2. Perform n-body Coulombic force summation.

  3. Perform kick–drift–kick leapfrog integration [Eqs. (23)–(25)].

  4. Account for probabilistic fragmentation (Sec. III D).

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.)

At each time step, a particle is injected at a randomized initial position (r0,z0,ψ0). The coordinate r0 is determined with a differential probability density dPjndAe, where jn 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 ΓM(r,z(r)), then the differential area of emission yields

dAe=2πrd=2πr1+(dzdr)2dr,
(20)

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

FIG. 3.

Radial probability of injection is shown in black. Non-dimensional current density emitted along the meniscus interface is shown in blue (see the subaxis in Fig. 2). The reference current density for non-dimensionalized is j=6.78Am2 The x axis represents the distance of a meniscus interface point to the axis of symmetry in non-dimensional coordinates, where rrc=1 represents the point where the meniscus attaches to the tip electrode. Results shown for the three currents in the legend.

FIG. 3.

Radial probability of injection is shown in black. Non-dimensional current density emitted along the meniscus interface is shown in blue (see the subaxis in Fig. 2). The reference current density for non-dimensionalized is j=6.78Am2 The x axis represents the distance of a meniscus interface point to the axis of symmetry in non-dimensional coordinates, where rrc=1 represents the point where the meniscus attaches to the tip electrode. Results shown for the three currents in the legend.

Close modal

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 ψ0 is chosen randomly, with uniform probability between ψ0[0,2π].

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,

v0=(vinn)n^.
(21)

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.

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,

miai=qi(EL+EC)=qiEl(ri)+jqiqj(rirj)4πε0|rirj|3,
(22)

where the Laplace field EL(ri) is defined as the background Laplace field originated by the potential difference between the electrodes ΓT and Γ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,

vin+12=vin+ai(rin)Δt2,
(23)
rin+12=rin+vin+12Δt2,
(24)
vin+1=vin+12+ai(rin+12)Δt2,
(25)

where vin and ain are the velocity and acceleration of the particle i at the step n and Δ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.

FIG. 4.

Schematic of the different computational regions used in the simulation (not to scale).

FIG. 4.

Schematic of the different computational regions used in the simulation (not to scale).

Close modal

Region I is a hemisphere of 5μ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 Δt=1016s. The time steps used in this paper below Δt=1.4x1012s 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.

FIG. 5.

Particle trajectories at different orders of magnitude time steps. Trajectory convergence analysis is shown in the subaxis plots with respect to the solution at time step Δt=1016. The slope of the error convergence plot is about 1.72 in logarithmic scale dlogϵdlogΔt.

FIG. 5.

Particle trajectories at different orders of magnitude time steps. Trajectory convergence analysis is shown in the subaxis plots with respect to the solution at time step Δt=1016. The slope of the error convergence plot is about 1.72 in logarithmic scale dlogϵdlogΔt.

Close modal

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,

dN=IqΔt.
(26)

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 1011 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/10th 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 10th 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μ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.

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 δt is given by the integration of the constant rate over the length of the time step,

p=1eδtτ,
(27)

where τ 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. EMIBF4 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 μ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.

FIG. 6.

Optimal mean lifetimes as a function of electric field for different cluster species for 327 nA. Cluster temperatures are 1000, 1350, and 810 K for dimers, trimers, and trimer products, respectively.

FIG. 6.

Optimal mean lifetimes as a function of electric field for different cluster species for 327 nA. Cluster temperatures are 1000, 1350, and 810 K for dimers, trimers, and trimer products, respectively.

Close modal
FIG. 7.

Comparison of numerical and experimental RPA curves for 324 and 327 nA, respectively. This simulation was composed of 40% monomers, 40% dimers, and 20% trimers. Experimental data collected by Dr. Catherine Miller.13 

FIG. 7.

Comparison of numerical and experimental RPA curves for 324 and 327 nA, respectively. This simulation was composed of 40% monomers, 40% dimers, and 20% trimers. Experimental data collected by Dr. Catherine Miller.13 

Close modal

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.

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.

TABLE I.

Parameters for the simulation cases.

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 16a 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 16a 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 (ΓE in Fig. 1) is shown in gray. The extractor has a circular aperture of 150μm radius and a thickness of 30μm.

FIG. 8.

Plot of the plume after 16 ns simulation time at 440 nA.

FIG. 8.

Plot of the plume after 16 ns simulation time at 440 nA.

Close modal

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

frac=nneutralsndimers+ntrimers100
(28)

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 

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μ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.

FIG. 9.

Change in kinetic vs electrostatic energy.

FIG. 9.

Change in kinetic vs electrostatic energy.

Close modal

The divergence angle is calculated by integrating over the solid angle (0.5° 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μm. The calculated angle is equivalent at both locations; thus, only the 250μ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° for the plot and 0.5° for the values reported in Table II. This achieves a higher resolution for the numerical value and less noise in the plot.

FIG. 10.

Divergence angles vs beam intensity (normalized current density) for all three simulations at 250μm from the emitter tip. Rotational symmetry is assumed. A 3° discretization is used for the graph. The simulation data are compared to experimental data from Lozano,11 and both respective data sets are curve fit using a super-Gaussian model.

FIG. 10.

Divergence angles vs beam intensity (normalized current density) for all three simulations at 250μm from the emitter tip. Rotational symmetry is assumed. A 3° discretization is used for the graph. The simulation data are compared to experimental data from Lozano,11 and both respective data sets are curve fit using a super-Gaussian model.

Close modal
TABLE II.

Divergence half-angles for the two confidence intervals and at the simulation boundary. The angles are in 0.5 steps due to the discretization.

Current (nA)80% (250 μm)99% (250 μm)
120 7.0° 12.0° 
324 8.0° 12.0° 
440 8.0° 12.0° 
Current (nA)80% (250 μm)99% (250 μm)
120 7.0° 12.0° 
324 8.0° 12.0° 
440 8.0° 12.0° 

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° 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.

TABLE III.

Divergence half-angles for the neutral particles and the two confidence intervals at the simulation boundary. Same discretization applies as in Table II.

Current (nA)80% (250 μm)99% (250 μm)
120 9.5° 22.5° 
324 10.0° 27.5° 
440 10.0° 29.0° 
Current (nA)80% (250 μm)99% (250 μm)
120 9.5° 22.5° 
324 10.0° 27.5° 
440 10.0° 29.0° 

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

f(θ)=Aexp[((θθt)2θ02)n],
(29)

where θ is the dependent variable, A=1 is the magnitude due to normalization, θt=0 the tilt in the angle as we assume perfect on axis firing, θ0=8.812°[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 θ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 θ0=13.18°[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=qEn=miv2rt,

rt=mivi2qiE(r)n=vi02+2(qi/mi)(ϕ0ϕ(r))(qi/mi)E(r)n,
(30)

where vi0 and ϕ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,ϕ0,ϕ(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) 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μm instead of 300μ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 and higher currents. The emission region length-scale r is itself a by-product of solving the electrohydrodynamic equations presented in Sec. II C but can be approximated to rγE2 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 exist for ionic liquids. Simulation tools, such as the one presented in this paper, could enable physics-based correlations of E 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.

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.

FIG. 11.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for low current simulation: 120 nA.

FIG. 11.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for low current simulation: 120 nA.

Close modal
FIG. 12.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for medium current simulation: 324 nA.

FIG. 12.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for medium current simulation: 324 nA.

Close modal
FIG. 13.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for high current simulation: 440 nA. The higher error when excluding space charge effects at higher currents is visible.

FIG. 13.

LHS: particle trajectories of a monomer, including space charge effects (dashed line), and only considering the Laplace field; RHS: Ratio of Coulomb to Poisson electric field for high current simulation: 440 nA. The higher error when excluding space charge effects at higher currents is visible.

Close modal

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.

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 findings13), these events are unimportant to the evolution of the beam in the present domain volume.

FIG. 14.

(a) Distribution of fragmentation events per nanosecond as a function of the electric field. (b) Approximate axial location for each simulation over the magnitude of the electric field.

FIG. 14.

(a) Distribution of fragmentation events per nanosecond as a function of the electric field. (b) Approximate axial location for each simulation over the magnitude of the electric field.

Close modal

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μ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.

FIG. 15.

Interpolated number flux of neutrals at 100μm (the location of the extractor electrode). Data are shown for the high current (440 nA) simulation. The black dots represent the (x,y) projection of the neutral location at z=100μm.

FIG. 15.

Interpolated number flux of neutrals at 100μm (the location of the extractor electrode). Data are shown for the high current (440 nA) simulation. The black dots represent the (x,y) projection of the neutral location at z=100μm.

Close modal
FIG. 16.

Neutral flux at different radii.

FIG. 16.

Neutral flux at different radii.

Close modal
FIG. 17.

Impact energies of neutrals at different radii.

FIG. 17.

Impact energies of neutrals at different radii.

Close modal
FIG. 18.

Densities of all particle species vs axial location for the 440 nA simulation.

FIG. 18.

Densities of all particle species vs axial location for the 440 nA simulation.

Close modal

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μm radial band for all currents studied. The average energy of impinging neutrals is below 10 eV for radial distances of 150μ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.

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

R=n1n2σ|v1v2|vol,
(31)

where n1 and n2 are the densities of the particles (Fig. 18), σ 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×1018m2 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 

FIG. 19.

Average velocities of all particle species vs axial location for the 440 nA simulation.

FIG. 19.

Average velocities of all particle species vs axial location for the 440 nA simulation.

Close modal

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μ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μ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.

FIG. 20.

Estimated collision rates for all species with monomers in the region close to the tip for the 440 nA simulation.

FIG. 20.

Estimated collision rates for all species with monomers in the region close to the tip for the 440 nA simulation.

Close modal

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 1μ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μ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.

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.

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.

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:

  1. 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

  2. Randomize the velocities of each atom of the cluster at 298 K.

  3. 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.

  4. Apply a constant electric field strength until the separation between the molecules of the cluster reaches the fragmentation threshold δf(sim). More about δ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, δf(sim) must be larger than the largest maximum separation that could occur without fragmentation. A δf(sim) value of 40 Å was selected in agreement with the previous work.10,27 Each cluster was simulated until this value of δ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 δf(post). The final chosen δf(post) value for dimers was 20 Å, which is in agreement with the recent work.10 The mean lifetime, τ, 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.

TABLE IV.

EMI-BF4 positive dimer mean lifetimes in units of picoseconds. Ellipses indicate that simulations were not run at these conditions due to the limitations of the computational resources used.

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 
1.
J.
Melngailis
, “Ion sources for nanofabrication and high resolution lithography,” in Proceedings of the 2001 Particle Accelerator Conference (IEEE, 2001), Vol. 1, pp. 76–80.
2.
J. B.
Fenn
,
M.
Mann
,
C. K.
Meng
, and
S. F.
Wong
, “
Electrospray ionization—Principles and practice
,”
Mass Spectrom. Rev.
9
,
37
70
(
1990
).
3.
P. C.
Lozano
,
B. L.
Wardle
, and
P.
Moloney
, “
Nanoengineered thrusters for the next giant leap in space exploration
,”
MRS Bull.
40
,
842
849
(
2015
).
4.
E.
Petro
,
A.
Bruno
,
P.
Lozano
,
L. E.
Perna
, and
D.
Freeman
, “Characterization of the TILE electrospray emitters,” in AIAA Propulsion and Energy 2020 Forum (AIAA, 2020), p. 3612.
5.
J. K.
Ziemer
,
C.
Marrese-Reading
,
S.
Arestie
,
N. R.
Demmons
,
R. E.
Wirz
,
A.
Collins
, and
M.
Gamero
, “Progress on developing LISA microthruster technology,” in AIAA Propulsion and Energy 2020 Forum (AIAA, 2020), p. 3609.
6.
D.
Krejci
,
A.
Reissner
,
T.
Schönherr
,
B.
Seifert
,
Z.
Saleem
, and
R.
Alejos
, “Recent flight data from IFM nano thrusters in a low earth orbit,” in 2019 International Electric Propulsion Conference (IEPC, 2019).
7.
X.
Gallud
and
P. C.
Lozano
, “
The emission properties, structure and stability of ionic liquid menisci undergoing electrically assisted ion evaporation
,”
J. Fluid Mech.
933
,
A43
(
2022
).
8.
K. J.
Terhune
,
L. B.
King
,
K.
He
, and
J.
Cumings
, “
Radiation-induced solidification of ionic liquid under extreme electric field
,”
Nanotechnology
27
(
37
),
375701
(
2016
).
9.
J.
Asher
,
Z.
Huang
,
C.
Cui
, and
J.
Wang
, “
Multi-scale modeling of ionic electrospray emission
,”
J. Appl. Phys.
131
(
1
),
014902
(
2022
).
10.
N.
Nuwal
,
V. A.
Azevedo
,
M. R.
Klosterman
,
S.
Budaraju
,
D. A.
Levin
, and
J. L.
Rovey
, “
Multiscale modeling of fragmentation in an electrospray plume
,”
J. Appl. Phys.
130
(
18
),
184903
(
2021
).
11.
P. C.
Lozano
, “
Energy properties of an EMI-Im ionic liquid ion source
,”
J. Phys. D: Appl. Phys.
39
(
1
),
126
(
2006
).
12.
C.
Pérez-Martínez
and
P. C.
Lozano
, “
Ion field-evaporation from ionic liquids infusing carbon xerogel microtips
,”
Appl. Phys. Lett.
107
(
4
),
043501
(
2015
).
13.
C. E.
Miller
and
P. C.
Lozano
, “
Measurement of the dissociation rates of ion clusters in ionic liquid ion sources
,”
Appl. Phys. Lett.
116
(
25
),
254101
(
2020
).
14.
D.
Krejci
,
F.
Mier-Hicks
,
R.
Thomas
,
T.
Haag
, and
P.
Lozano
, “
Emission characteristics of passively fed electrospray microthrusters with propellant reservoirs
,”
J. Spacecr. Rockets
54
(
2
),
447
458
(
2017
).
15.
N. M.
Uchizono
,
A. L.
Collins
,
A.
Thuppul
,
P. L.
Wright
,
D. Q.
Eckhardt
,
J.
Ziemer
, and
R. E.
Wirz
, “
Emission modes in electrospray thrusters operating with high conductivity ionic liquids
,”
Aerospace
7
(
10
),
141
(
2020
).
16.
M.
Tajmar
,
A.
Genovese
, and
W.
Steiger
, “
Indium field emission electric propulsion microthruster experimental characterization
,”
J. Propul. Power
20
(
2
),
211
218
(
2004
).
17.
C.
Coffman
,
M.
Martínez-Sánchez
,
F. J.
Higuera
, and
P. C.
Lozano
, “
Structure of the menisci of leaky dielectric liquids during electrically-assisted evaporation of ions
,”
Appl. Phys. Lett.
109
(
23
),
231602
(
2016
).
18.
M.
Gamero-Castaño
and
M.
Magnani
, “
Numerical simulation of electrospraying in the cone-jet mode
,”
J. Fluid Mech.
859
,
247
267
(
2019
).
19.
M. J.
Davis
,
A. L.
Collins
, and
R. E.
Wirz
, “Electrospray plume evolution via discrete simulations,” in Proceedings of the 36th International Electric Propulsion Conference, IEPC-2019, Vienna, Austria (IEPC, 2019), Vol. 590.
20.
R.
Yokota
and
L. A.
Barba
, “Treecode and fast multipole method for N-body simulation with CUDA,” in GPU Computing Gems Emerald Edition (Morgan Kaufmann, 2011), pp. 113–132.
21.
B. D.
Prince
,
C. J.
Annesley
,
R. J.
Bemish
, and
S.
Hunt
, “Solvated ion cluster dissociation rates for ionic liquid electrospray propellants,” in AIAA Propulsion and Energy Forum and Exposition 2019 (AIAA, 2019).
22.
C. S.
Coffman
,
M.
Martínez-Sánchez
, and
P. C.
Lozano
, “
Electrohydrodynamics of an ionic liquid meniscus during evaporation of ions in a regime of high electric field
,”
Phys. Rev. E
99
(
6
),
063108
(
2019
).
23.
J. V.
Iribarne
, “
On the evaporation of small ions from charged droplets
,”
J. Chem. Phys.
64
(
6
),
2287
(
1976
).
24.
C. S.
Coffman
, “Electrically-assisted evaporation of charged fluids: Fundamental modeling and studies on ionic liquids,” Ph.D. thesis (Massachusetts Institute of Technology, 2016).
25.
B.
Joe
, “
GEOMPACK—A software package for the generation of meshes using geometric algorithms
,”
Adv. Eng. Softw. Workst.
13
(
5–6
),
325
331
(
1991
).
26.
T. M.
Coles
,
T. P.
Fedkiw
, and
P. C.
Lozano
, “Investigating ion fragmentation in electrospray thruster beams,” in 48th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit 2012 (AIAA, 2012).
27.
T. M.
Coles
and
P. C.
Lozano
, “Investigating efficiency losses from solvated ion fragmentation in electrospray thruster beams,” in 49th AIAA/ASME/SAE/ASEE Joint Propulsion Conference (AIAA, 2013), Vol. 1, Part F.
28.
P. C.
Lozano
, “
Energy properties of an EMI-Im ionic liquid ion source
,”
J. Phys. D: Appl. Phys.
39
(
1
),
126
134
(
2006
).
29.
N.
Uchizono
,
A.
Collins
,
C.
Marrese-Reading
,
S.
Arestie
,
J.
Ziemer
, and
R.
Wirz
, “
The role of secondary species emission in vacuum facility effects for electrospray thrusters
,”
J. Appl. Phys.
130
,
3301
(
2021
).
30.
J.
Magnusson
,
A.
Collins
, and
R.
Wirz
, “
Polyatomic ion-induced electron emission (IIEE) in electrospray thrusters
,”
Aerosp. Sci. J.
7
,
153
(
2020
).
31.
E. M.
Petro
,
C. E.
Miller
,
J.
Schmidt
, and
P.
Lozano
, “Development of an electrospray fragmentation model for kinetic plume modeling,” in The 36th International Electric Propulsion Conference (IEPC, 2019), pp. 1–10.
32.
B. D.
Prince
,
P.
Tiruppathi
,
R. J.
Bemish
,
Y. H.
Chiu
, and
E. J.
Maginn
, “
Molecular dynamics simulations of 1-ethyl-3-methylimidazolium bis[(trifluoromethyl)sulfonyl]imide clusters and nanodrops
,”
J. Phys. Chem. A
119
(
2
),
352
368
(
2015
).
33.
B. D.
Prince
,
A. L.
Patrick
,
C. J.
Annesley
,
R. J.
Bemish
,
S. W.
Miller
,
M. L.
Hause
, and
K. M.
Vogelhuber
, “A combined experimental and theoretical treatment of ionic liquid thermal dissociation,” in 53rd AIAA/SAE/ASEE Joint Propulsion Conference 2017 (AIAA, 2017), pp. 1–17.
34.
M.
Schroeder
, “Numerical characterization of fragmentation in ionic liquid clusters,” M.S. thesis (Massachusetts Institute of Technology, 2021), arXiv:2111.08448.
35.
J. N.
Canongia Lopes
,
J.
Deschamps
, and
A. A. H.
Pádua
, “
Modeling ionic liquids using a systematic all-atom force field
,”
J. Phys. Chem. B
108
(
6
),
2038
2047
(
2004
).
36.
A. P.
Thompson
,
H.
Metin Aktulga
,
R.
Berger
,
D. S.
Bolintineanu
,
W.
Michael Brown
,
P. S.
Crozier
,
P. J.
in ’t Veld
,
A.
Kohlmeyer
,
S. G.
Moore
,
T. D.
Nguyen
,
R.
Shan
,
M. J.
Stevens
,
J.
Tranchida
,
C.
Trott
, and
S. J.
Plimpton
, “
LAMMPS—A flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales
,”
Comput. Phys. Commun.
271
,
108171
(
2022
).
37.
J. N.
Canongia Lopes
and
A.
Pádua
, “
CL&P: A generic and systematic force field for ionic liquids modeling
,”
Theor. Chem. Acc.
131
,
3330
(
2012
).
38.
D. J.
Evans
and
B. L.
Holian
, “
The Nose–Hoover thermostat
,”
J. Chem. Phys.
83
,
4069
(
1985
).
39.
A.
Thuppul
,
A. L.
Collins
,
P. L.
Wright
,
N. M.
Uchizono
, and
R. E.
Wirz
, “
Mass flux and current density distributions of electrospray plumes
,”
J. Appl. Phys.
130
(
10
),
103301
(
2021
).
40.
R. G.
Forbes
and
N. N.
Ljepojevic
, “
Liquid-metal ion source theory: Electrohydrodynamics and emitter shape
,”
Surf. Sci.
266
(
1–3
),
170
175
(
1992
).
41.
L. W.
Swanson
,
G. A.
Schwind
,
A. E.
Bell
, and
J. E.
Brady
, “
Emission characteristics of gallium and bismuth liquid metal field ion sources
,”
J. Vac. Sci. Technol.
16
,
1864
(
1979
).
42.
N.
Takahashi
and
P. C.
Lozano
, “Computational investigation of molecular ion evaporation in electrospray thrusters,” in 44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit (AIAA, 2008).
43.
N.
Takahashi
and
P. C.
Lozano
, “Atomistic numerical approach to ion evaporation from a tungsten surface for electrospray thrusters,” in 45th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit (AIAA, 2009).