Molecular Dynamics Electrospray Simulations of Coarse-Grained Ethylammonium Nitrate ( EAN ) and 1-Ethyl-3-Methylimidazolium Tetrafluoroborate ( EMIM-BF 4 )

In this work, the use of molecular dynamics as a predictive tool for modeling the atomistic behavior of electrospray propulsion is discussed. 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF4) and ethylammonium nitrate (EAN) were considered as two limits of ionic liquid (IL) propellants that tend to operate in an ion versus a droplet mode. The emission modes were found to depend on the electro-chemical properties of the IL propellant. The aprotic EMIM-BF4-based electrospray emitted primarily monomers and trimers as the dominant species and only small quantities of droplets. In contrast, trimers were the dominant emitted species in the protic EAN emissions with a significantly large contribution from droplets towards the total emission current, suggesting that EMIM-BF4-based colloid thrusters operate in ion mode and EAN-based devices operate in the droplet mode. Furthermore, the formation of the Taylor cone was found to depend on the mass flow rate and the external electric field strength. This paper provides a framework that can be extended for use to simulate any other ILs or their combinations.


Introduction
Colloid thrusters are a subset of electrostatic propulsion devices that use relatively benign non-metallic ionic liquids (ILs), which are defined as salts having a melting point less than room temperature (373 K).ILs are used as propellants in colloid thrusters with the advantage of producing higher thrust-to-power ratios compared to the traditional liquid cesium-based field emission electric propulsion (FEEP) [1] thrusters.When a capillary fed by an ionic liquid is placed in the presence of an external electric field, a fine jet of ions is emitted from a Taylor [2] cone formed at the head of the capillary, and the emitted ions are further accelerated by the electric field.The emission from the Taylor cone is also known as an electrospray and its applications are not just limited to colloid thrusters.Electrosprays have been used extensively in the field of micro-fabrication because they allow for controlled deposition [3] for applications involving microfilm deposition [4], microcircuit manufacturing [5] and ion beam lithography [6].The performance and operability of colloid thrusters have already been tested in space on the NASA ST-7 LISA Pathfinder mission [7,8].
The operational mode of the colloid thruster is dependent on whether the ILs breakdown into either individual ion species or sub-micron-sized droplets (colloids) that are accelerated by an electric field, generating high exhaust velocities.Ionic liquids are classified into either aprotic or protic types.The imidazolium-based ILs, such as 1-ethyl-3-methylimidazolium tetrafluoroborate (EMIM-BF 4 ), fall under the aprotic class of ILs, and as shown by Borner et al. [9], EMIM-BF 4 operates in the ion mode [10].On the other hand, ethylammonium nitrate (EAN) is classified as a protic IL, which are known to be environmentally more benign than other organic solvents [11] and easier to synthesize compared to aprotic ILs [12].EAN has been reported to emit in the droplet mode, producing a jet of nanodroplets [13].The electro-chemical properties of the ILs used as propellants and their mass flow rates determine the mode of operation of the colloid thruster.Extensive research has been performed to determine the emission modes of the colloid thrusters.Since future propellants may use mixtures of protic and aprotic types of ILs, molecular dynamics (MD) electrospray simulations can be used to analyze the differences in the emission behavior specific to these two types of ILs.
A large number of ILs can be used as propellants for electrosprays [14,15], and the availability of a large number of inter-atomic potentials also makes it possible to model virtually any of the ILs.However, MD simulations scale on the order of O(N 2 ), where, N is the number of atoms, making MD very expensive to model large simulations (>10,000 atoms) typically used to study electrosprays.In this work, effective field coarse graining method [16,17] was used to derive coarse-grained potential of EAN, and its efficacy for the electrospray processes was verified by successfully performing MD simulations of Taylor cone formation and subsequent emission of ion-species.Previous works [16][17][18] have derived coarse-grained (CG) potentials to verify bulk IL physical properties or have used approximations to modify the CG potential to perform electrospray simulations [19,20].
Prior to designing colloid thrusters, it is necessary to have a simulation-based ability to predict thrust, required voltages and the mode of operation.MD is a particle-based method that allows one to understand the nature of electrosprays at an atomistic level [21].In this work, a bottom-up approach was used to analyze the correlation between the electro-chemical properties, such as the Coulomb interaction of the ILs, and the emission of ion species from the capillary, the required extrusion voltage and the formation time of the Taylor cone.This work provides key insight into the differences in the spatial characteristics of monomer and droplet emission, which has previously only been postulated [22].
Although MD is a mature computational technique, there are several challenges that need to be addressed when investigating the properties of ionic liquid electrosprays.The stability of the Taylor cone also depends on the externally-applied electric field.Colloid thruster experiments provide important information regarding the species-specific emission currents [23] and the thrust generated for the applied potential across the device, but they lack the specificities of the shape and strength of the electric field generated around the mouth of the capillary.It should also be noted that while the experiments are conducted over length and time scales of micro-meters and micro-seconds, the current computational capability only allows one to simulate electrosprays at length and timescales of nano-meters and nano-seconds.Thus, it becomes incumbent for the MD simulations to use external boundary conditions that recreate the Taylor cone formation and emission characteristics similar to the experiments without the knowledge of the spatial and temporal evolution of the electric fields generated during the experiments.Therefore, a discussion of the appropriate electric field boundary conditions is also included.
The remainder of the paper is as follows.In Section 2, the methodology used to model the coarse-grained potential of EAN is presented.The simulation domain and the boundary conditions required to generate electrospray emission using the MD simulations are discussed in Section 3. The differences between EAN and EMIM-BF 4 in terms of their Coulomb interactions are analyzed and presented in Section 4. The formation of the Taylor cone and the emission currents obtained from the EAN and EMIM-BF 4 electrospray simulations are discussed in Section 5. Conclusions and remaining challenges are outlined in Section 6.

Coarse-Graining of EAN Using EFCG
Since a large number of molecules is required to simulate non-equilibrium processes such as electrosprays in MD, these simulations are computationally very expensive.Researchers such as Borner et al. [24] have overcome this problem by using coarse-grained potentials to model electrospray simulations of EMIM-BF 4 .The coarse-grained potential of EMIM-BF 4 was obtained by modifying the CG potential of EMIM-NO 3 derived by Wang et al. [16][17][18] from the all-atom potential using the effective field coarse-graining (EFCG) method.The EFCG [16,17] method works on the principle of simplifying non-essential group of atoms and their respective degrees-of-freedom to a single coarse-grained site [18].For EAN electrospray simulations, it was postulated that modeling the relative dynamics of hydrogen or oxygen atoms may not be vital.Therefore, for the CG model of EAN, it was proposed that the atoms forming CH 3 , CH 2 , NH 3 and NO 3 should be simplified to CG sites E2, E1, N1 and NN, respectively, as shown in Figure 1a,b, respectively.This allows one to simulate EAN using only four CG sites as opposed to 15 atoms required for the all-atom simulations.The EFCG method lays down the rules to obtain a CG potential from the existing all-atom potentials.The location of each CG site is taken as the center-of-mass of its constituent atoms, thus satisfying the linear relationship criteria between the location of the CG sites with its constituent atoms.Constituent atoms forming one CG site cannot belong to any other CG site, and each CG site should have at least one constituent atom.Furthermore, the force acting on one CG site should be equal to the sum of forces acting on constituent atoms forming that CG site.The mass and partial charge of the CG site is calculated as the sum of masses and partial charges of the constituent atoms forming that site, as discussed in other previous work [25].This conserves kinetic energy and electric charges when transitioning from the all-atom to the CG model.

Determination of CG Bonded Interaction Parameters
The coarse-graining procedure requires separate treatments to derive bonded and non-bonded interactions from all-atom potentials.Bonded interactions account for the energy of covalent bonds, angles and dihedral terms.The non-bonded interactions have contributions from the van der Waals and Coulomb terms.To obtain the bonded interaction parameters, an all-atom simulation of 64 EAN ion-pairs was energy minimized and equilibrated using a grand-canonical ensemble (NPT) followed by a canonical ensemble (NVT) with a 0.5-fs time step, for a simulation time of 2.5 ns, each at a temperature of 295 K.After performing equilibration, 4000 samples of atomic positions were collected by running the system with a micro-canonical ensemble (NVE) over three million time steps.Each sample represents a snapshot of the position of all atoms of 64 EAN ion-pairs every 750 time steps (0.375 ps).For each sample, 64 molecules per sample, using a post-processing code, bond-lengths between groups forming methylene(E2)-methyl(E1) and methyl(E1)-ammonium(N1) were obtained along with the values of the angle formed by the methylene-methyl-ammonium group.Therefore, for each of the bonds and covalent angle, 256,000 values were obtained from which probability distribution functions (PDFs) were generated for the two bonds (between E1-E2 and N1-E1) and one angle (E1-E2-N1).
The generated PDFs were then approximated using a curve fit of the form, where µ is the median of the sampled data and σ is the its standard deviation, both of which were obtained from the all-atom MD simulations.This form of the normal distribution is comparable to a Boltzmann distribution, where the bonded interaction is given by V b (r), k B is the Boltzmann constant and T = 295 K is the canonical temperature.The potential energy for the bonded interactions due to bond stretching is calculated using a harmonic function of the form, and for the valence angle potential interactions, where k r and k θ are the harmonic spring constants of the CG covalent bonds and angles, respectively.Substituting the harmonic functions from Equations ( 3) and (4) into Equation ( 2) and comparing the resulting expression with the normal distribution from Equation (1), it is observed that r and θ are the bond distance and angle formed by the three center of masses, respectively.The r 0 and θ 0 represent the equilibrium values for the bond and angles and are equal to the median µ of the bonds and angles from the sampled data.

Determination of CG Non-Bonded Interaction Parameters
The non-bonded potential has contributions from the van der Waals dispersion and Coulomb electrostatic interactions, The van der Waals contribution is simulated using the Lennard-Jones potential, where n and m are the total number of atoms on the two molecules under consideration, σ is the distance at which the energy is zero between the atom pair and is the lowest energy between the atom pair.The values of σ and are calculated from ab initio calculations [26], and r ij is the distance between the atom pair.
In MD, the electrostatic energy is calculated using Coulomb's relationship, where r ij is the distance between atoms i and j, C is Coulomb's constant, q i , q j are the user-defined partial charges on atoms i and j, which are obtained from ab initio calculations [26], respectively.When the CG assumptions are made, the information regarding the van der Waals and the intra-molecular Coulomb interactions is lost, but the ability to model the inter-molecular Coulomb interactions is retained since the partial charges of each CG site is approximated as the sum of partial charges of the constituent atoms forming that CG site.In the proposed EAN CG model, there are four CG sites resulting in 10 distinct pairs (e.g., E1-E1, E1-E2, E1-N1, N1-NN, etc).
Combining a group of constituent atoms that occupy a finite volume in the simulation domain to a single point leads to losses in the spatial variations that exist in the all-atom MD domain.This loss of spatial variation is important because it influences the intra-molecular component of non-bonded interaction forces.To account for different spatial configurations observed in the all-atom EAN system, a Monte Carlo approach using 8000 different molecular orientations was implemented by placing the constituent atoms forming one of the CG pairs at a fixed point and placing the constituent atoms of the second CG in a random orientation, radially 2.0 Å away from the center of mass of the first group of atoms.Keeping the first group of atoms fixed, the second group is moved radially 0.02 Å outwards along the line joining the center of masses of these two molecules for all 8000 orientations.At every increment, the non-bonded van der Waals and Coulomb interactions are recorded as a function of radial distance for each CG pair combination.This process is repeated up to a radial distance of 20 Å, and the values obtained from the 8000 samples are averaged.From the MD simulations, the van der Waals and Coulomb force exerted between the constituent atoms forming the two CG sites are obtained.However, this Coulomb force is the sum of inter-, as well intra-molecular electrostatic interactions and is referred to as the total Coulomb force, F total .
The inter-molecular Coulomb force is then subtracted from the sum of van der Waals and F total to obtain the effective force, F e f f ective , where the inter-molecular Coulomb force is, such that Q 1 and Q 2 are the sum of the partial charges on the constituent atoms forming the first and second CG site and R 12 is the radial distance between the center of masses of the two CG sites.Consider two CG sites, E1 from the cation and NN, which forms the anion.The Coulomb and van der Waals forces averaged over 8000 configurations are shown in Figure 2a,b.The total Coulomb force, F total , obtained from Equation (8), and the inter-molecular Coulomb force, F inter-molecular , obtained using Equation (10), are shown in Figure 2a.The F inter-molecular is then subtracted from the sum of van der Waals and F total , provided as Equation ( 9), which is shown in Figure 2b.This effective force and its corresponding potential is inputted to the MD simulations as a tabulated potential to model the CG non-bonded interaction potential.Using the bonded and non-bonded CG potentials, 4000 samples of atom positions were collected, similar to that from the all-atom simulations.The probability distribution function (PDFs) of bonds connecting sites E1-N1, E1-E2 and angle N1-E1-E2 were generated and compared with those from the all-atom simulations.It was found that while the distribution of angles agreed well, the median value for both CG bonds was 0.0075 Å less than the specified equilibrium bond distances.Thus, a correction was applied to the bonded CG model, such that the equilibrium bond distances for both CG bonds were increased by 0.0075 Å, and the final values are shown in Table 1.After implementing the second correction, equilibration and NVE simulations were repeated.Using the simulations performed with the corrected bonded CG potential, the comparison between the CG bond and angle distributions are shown in Figure 3a-c.It can be observed that the distribution of the CG bonds between sites E1-E2, E1-N1 and angle N1-E1-E2 are in excellent agreement with the distributions obtained from the all-atom simulations.To validate the non-bonded interactions of the derived CG potential, 64 CG EAN ion-pairs were placed in a periodic cell of size 25 × 25 × 25 Å, and the system was energy minimized.This was followed by 2.5 ns of NPT and NVT ensemble equilibration simulations, both with a time step of 0.5 fs, similar to the all-atom MD simulations.From these equilibrium simulations, the mass density of the CG system was found to be 1080 kg/m 3 , which was lower than the desired mass density of 1216 kg/m 3 because the non-bonded interactions were too repulsive at larger than expected inter-CG site distances.To correct this, the non-bonded interaction force values for all atom pairs were shifted by 0.25 Å.The system was re-equilibrated after this correction, and the desired mass density of 1219 kg/m 3 was obtained.To verify the accuracy of the CG potential, RDFs from the CG simulations were compared with those obtained from the MD all-atom simulations, as shown in Figure 4. To obtain the all-atom site RDFs, the center of mass of the constituent atoms forming that group were considered.The RDFs are typically used to ensure the accuracy of long-range structural order of any liquid system MD simulations, and therefore, the focus is on the comparison of the first and second co-ordination shell radii and not their magnitudes.For all the CG pairs, the first and second co-ordination shells were generated with reasonable accuracy in comparison with the all-atom simulation, which validates the non-bonded interaction parameters for EAN CG model.The accuracy of the CG EAN potential was also verified by reproducing the electrical conductivity of 2.39 S/m, which is in agreement with the experimental value of 2.47 [13].The tabulated potential files contain the details of the non-bonded interactions for the 10 possible CG combinations and are provided as Supplementary Materials.The prime objective of coarse graining was to achieve a lower computational cost of the electrospray MD simulations.For a 64-molecule periodic cell, the number of simulated particles was reduced from 960-256.Using 8 Xeon E5-2680 @ 2.7-GHz processors on the STAMPEDE system, a speed up of approximately five-times was obtained.Since the electrospray simulations typically have more than 19,000 molecules, the potential speed up and the savings in the computational cost are beneficial.
Similar to the coarse-graining performed for EAN, Wang et al. [17] have reduced the non-essential degrees-of-freedom in the EMIM-NO 3 ion-pairs.Using CG, the number of required particles to simulate one molecule was reduced from 27 atoms to just five CG sites.The same CG potential was modified by replacing the anion from NO 3 to BF 4 to obtain the CG potential for EMIM-BF 4 .The two CH 3 groups were simplified to two single sites 'M1' and 'M2', respectively, as shown in Figure 1c,d, respectively.The largest computational savings is observed by reducing the large imidazolium group to a single CG site.The CH 2 group in the cation and the BF 4 anion were simplified to sites 'MR' and 'BF4', respectively.Equilibration and electrospray simulations of CG EMIM-BF 4 were performed for this work using the potential derived by Wang et al. [16].Borner et al. [24], also modeled the electrical conductivity of CG EMIM-BF 4 in MD simulations and found good agreement with the experimental values.

Numerical Parameters and Boundary Conditions Used for the MD Electrospray Simulations
The specification of the boundary conditions necessary to generate appropriate applied external electric fields and the MD simulation setup are presented next.

Boundary Conditions Used to Obtain External Electric Field
The external electric fields play an important role in determining the emitted ion species and the general shape of the structure of Taylor cone at the mouth of the capillary.The external electric field implemented on the electrospray simulations was calculated using the potential, φ l , derived by solving Laplace's equation, The boundary conditions used to solve Equation ( 11) are shown in Figure 5a.A three-dimensional schematic of the simulation domain is shown in Figure 5b.The entire outer surface of the capillary was grounded using a Dirichlet boundary condition of φ = 0 V.The simulation domain was non-periodic in all three dimension, and a Neumann boundary condition was applied on all three boundaries, as shown in Figure 5a.Therefore, the total electric potential felt by the i-th CG site would be, such that Q i and Q j are the partial charges on the CG sites, n is the total number of CG sites and R c = 20 Å is the cut-off radius used for the simulations.Equation (11) was solved using the generalized minimal residual (GMRES) method [27] with a grid-size (h x , h y , and h z ) equal to the Coulomb cut-off radii, R c = 20 Å.Note that the contribution from the second term of Equation ( 12) to the total electric potential on a CG site only needs to be computed once, at the zeroth time step, whereas the first term of Equation ( 12) was computed every time step in the MD simulation.Grid convergence for Equation (11) was tested by reducing the grid size from 20-5.0 Å.It was observed that the solution to Equation ( 11), using the boundary setup and conditions shown in Figure 5 and an extraction voltage of −17 V, generated an external electric field with normal, as well as radial components, as shown in Figure 6a,b, respectively.Both, the normal and radial electric field components are strongest at the meniscus or mouth of the capillary and become gradually weaker away from the mouth of the capillary.It was found that the presence of the radial electric field near the capillary meniscus is necessary to generate the Taylor cone shape [28].The unique shape of the electric field also effects the stability of the Taylor cone and the type of ions species that are emitted from it.

Details of the MD Simulation Domain
The accuracy of the long-range Coulomb interaction models is vital to the prediction of the Taylor cone formation and emission characteristics [29].Various long-range Coulomb models, such as the direct Coulomb (DC) [28], shifted force Coulomb sum (SFCS) [30] and the particle-particle particle-mesh (PPPM) [31] methods, can be used to simulate ionic liquid systems.However, as shown in other recent work [28], the highest accuracy is obtained only when the DC method with a relatively large cut-off radius of 20 Å is used.Therefore, all simulations were performed in LAMMPS [32] using the DC method to calculate all electrostatic non-bonded interactions with a cut-off radius of R c = 20 Å.
The simulations results discussed in this paper were performed using the boundary conditions and domain dimensions shown in Figure 5.The extraction voltages used for the EAN and EMIM-BF 4 simulations were φ = −60 and −17 V, respectively.The extraction voltage of −17 V was used for EMIM-BF 4 as it generates potential gradients comparable to those that occur over the length scale of the experiments.It was observed that an extraction voltage of −17 V did not cause emission from EAN.Instead, by trial and error, it was found that a starting extraction voltage of −60 V generated emission from EAN.The capillary of radius 56 Å and height 275 Å is first filled with 9455 ion-pairs for EMIM-BF 4 simulations and 19,810 ion-pairs for the EAN simulations.The ion-pairs in the capillary were first energy minimized using the system potential energy as the criteria to remove any non-physical CG site or bond overlaps that may exist.The ion-pairs were then equilibrated to a temperature of 297 K using the canonical (NVT) ensemble.A repulsive wall was used to generate a mass flow rate of 1.22 × 10 −12 and 1.14 × 10 −11 kg/s for the EMIM-BF 4 and EAN electrospray simulations, respectively.All extrusion simulations were performed using a time step of 0.5 fs.The emitted ion-pairs were then sampled at the extractor ring, and the emission current was calculated.The emitted positive ion species were classified based on the number of cations present in the aggregate.Using the definition, (EMIM − BF 4 ) n × EMIM + , if n was zero, the aggregate was termed a monomer, and if n = 1 or 2, the aggregate was defined as a dimer or trimer, respectively.The emission currents are shown as cumulative moving averages, where each data set represents the data obtained over 20 ps.

Comparison of the Cation-Anion Separation Energy of EAN vs. EMIM-BF 4
As stated previously, EMIM-BF 4 belongs to the aprotic class of ionic liquids, whereas EAN belongs to the protic class of ionic liquids.These two ILs show significant differences in their molecular structures and electro-chemical behaviors.Ionic liquids belonging to the protic ILs (PIL) subset are formed by proton transfer due to stoichiometric combination of a Brønsted acid and base [11,33].EAN is reported to "self-assemble into a sponge-like nanostructure in bulk" [34] and form separate polar and apolar domains.The polar and apolar domains consist of charged ammonium and nitrate groups and long chains of alkyl groups in the EAN bulk, respectively.EMIM-BF 4 in comparison is a larger molecule than EAN and does not show the tendency to form separate polar and apolar domains in the bulk.The average number of anions surrounding a cation, also known as the coordination number of EAN, was observed to be 1.9, and that of EMIM-BF 4 was 1.5, suggesting that the bulk structure of EAN is more closely packed than of EMIM-BF 4 .In other previous work [35], it was observed that EAN contains an extensive network of hydrogen bonds, which is not present in EMIM-BF 4 , possibly making the emission of ion-pairs from the EAN bulk more difficult to achieve than EMIM-BF 4 for the same extraction voltage.
As mentioned in the Introduction, EMIM-BF 4 operates in the ion emission mode during the electrospray process.In contrast, Alonso-Matilla et al. [13] has shown that ethylammonium nitrate (EAN) produces stable emission of fine droplets as small as 1 nm in diameter, which is preferable in certain electrospray operations because it results in a better charge-to-mass ratio, leading to more efficient colloid thrusters.The physical and chemical properties of EMIM-BF 4 are different from those of EAN.With respect to mass density, surface tension and conductivity, the values for EMIM-BF 4 compared to EAN are 1240 vs. 1216 kg/m 3 , 0.054 vs. 0.047 N/m, and 1.31 vs. 2.47 S/m, respectively.Borrajo-Pelaez and Gamero-Castano [36,37] have observed that the charge-to-mass ratio of EAN droplets was approximately two times higher than that for EMIM-BF 4 .When the same acceleration potential (applied voltage difference between the electrodes) was applied, the average droplet velocity of EAN was observed to be higher than that of EMIM-BF 4 (10 vs. 7.6 km/s respectively).However, it has also been observed that EAN requires roughly three-times higher extraction potential to generate emission compared to EMIM-BF 4 .
The external electric field generates the force required to separate ions from the IL bulk by acting on the partial charges of the individual CG sites.The Coulomb attraction between the cation and anion determines the energy required to separate and move the ion species in the presence of their respective electric field.The Coulomb attraction of the system was analyzed by simplifying the extrusion simulation as a two-step process.First, the applied electric field has to overcome the energy of the ionic bond and separate the cation and anion and, second, move the cation away from the bulk towards the cathode.The bulk simulations, schematics of which are shown in Figure 7, were used to understand the change in Coulomb energy during both steps of the extrusion process.The height, h, represents the distance of the meniscus from the platinum base on which 1000 IL ion-pairs are resting.The value of h was 53 and 79 Å, for EAN and EMIM-BF 4 , respectively.The Coulomb interaction measurements were found to be sensitive to inter-particle distances as small as 0.01 Å.
To simulate single ion-pair/ion emission simulations, a cation-anion ion-pair at the center of an IL bulk of 1000 ion-pair was selected and artificially moved through the IL bulk towards the cathode in the presence of an applied electric field.The ion-pair was moved 0.5 Å every time step, normal to the meniscus, in the direction of the applied electric field, and after every movement at each time step, the system was energy minimized to avoid atom and bond overlaps due to the artificial displacement of the ion-pair.The potential, total (sum of potential and kinetic energy) and Coulomb energy of the entire system were monitored to find the changes in the interaction energy due to the motion of the ion-pair first, through the bulk and then through the externally-applied electric field.This simulation was repeated by artificially moving only the cation of the ion-pair towards the meniscus and then towards the cathode using the same procedure.These simulations were performed for EAN and EMIM-BF 4 using different constant electric field strengths of 8.0 and 3.0 V/nm, respectively.The amount of energy required to separate the cation-anion pair can be estimated from the difference in the Coulomb energy to move just the cation versus that necessary to move the ion-pair, measured at the meniscus, as shown in Figure 8 (it was observed that the contributions from the covalent bond and angles and van der Waals energy were not significant, and the potential energy of the system was primarily only due to Coulomb interactions).The solid and dashed lines represent the energy values for movement of the entire ion-pair and of only the cation, respectively.Large spatial fluctuations are observed in both EAN and EMIM-BF 4 .These fluctuations occur when the artificially-displaced test ion-pair disrupts the long-range order of the liquid.However, it should be noted that the smaller oscillations between 30 and 50 Å in Figure 8a are not noise, but occur when the artificially-displaced test ion-pair encounters an alkyl aggregate within the bulk in EAN.These smaller fluctuations were observed in the region closer to the liquid-vacuum interface rather than deeper inside the EAN bulk and were not witnessed in EMIM-BF 4 .Comparing the differences in the energy due to the motion of the ion-pair versus the motion of only the cation, there is a difference of 4.0 and 1.5 eV for EAN and EMIM-BF 4 , respectively.This explains the fundamental differences that exist in the emission characteristics of EAN and EMIM-BF 4 ILs.The extensive network of hydrogen bonds existing in EAN provides resistance to the separation of cations and anions [35].However, it should be noted that only the normal electric field was used in these bulk simulations to separate the cation-anion pair, and the presence of both the normal and radial fields will make the emission of ions easier.

Evolution of Taylor Cone Structure of EAN vs. EMIM-BF 4
The evolution of Taylor cone formation generated by EMIM-BF 4 ion-pair emissions is shown in Figure 9, using snapshots obtained from the MD simulations.Up to 400 ps, only a few loose ions near the meniscus are emitted, as shown in Figure 9a, because the combination of low mass flow rate and low extrusion voltage create a low emission rate out of the capillary.Even at 400 ps, the apex shape formation is evident, but is not developed.At 500 ps, the initial stages of the formation of a Taylor cone structure can be seen at the mouth of the capillary.The combination of the strongest electric fields, both radial and normal, near the mouth of the capillary allows the Taylor cone to become fully formed, but the radial electric field also leads to emission of loose ions from the sides of the cone structure, as shown in Figure 9b.The cone structure is fully developed after 600 ps, as shown in Figure 9c.Once the Taylor cone structure is fully formed, the emission of ion-pairs occurs primarily from the apex, which is in a region of lower electric field strength.This allows the emission of larger ion-species to occur beyond a 600-ps simulation time.
Similarly, the evolution of the Taylor cone structure from the EAN simulation is shown in Figure 10.Due to a much higher flow rate compared to the EMIM-BF 4 simulations, a large amount of IL bulk is seen emitted outside the capillary even at 200 ps.However, as shown in Figure 10a, due to the stronger Coulomb potential in the EAN bulk, the emitted IL does not yet have a cone-like structure.
At the same time, loose ion-pairs are emitted from the top of the IL bulk, which is present outside the capillary.However, at 300 ps, the radial electric fields at the mouth of capillary start generating emissions in the lateral direction, gradually forming the Taylor cone shape, as shown in Figure 10b.Again, the formation of large droplets can be observed near the apex of this cone-like structure due to the presence of weaker electric fields there.The Taylor cone becomes fully formed at 400 ps, as shown in Figure 10c, with a steady emission of large droplets observed from its apex.

Comparison of Species-Specific Electrospray Emission Currents
The emission currents were sampled at the extractor for each emitted ion-species from the EMIM-BF 4 and EAN electrospray simulations.For EMIM-BF 4 , the smaller ion species, such as the monomers and dimers, emitted during the initial stages of the Taylor cone formation are sampled after 500 ps, as shown in Figure 11a.However, once the Taylor cone is fully formed after 600 ps, larger species such as the trimers are emitted followed by much larger droplets.It should also be noted that some of the larger droplets undergo further breakdown during their travel from the apex of the Taylor cone towards the extractor plane due to the weaker, but nonetheless effective electric field.
The larger droplets mostly break apart into monomers and trimers.Thus, the droplets have the lowest contribution towards the total emission currents.The monomers and trimers form the most dominant emitted species, and the droplets are the least dominant emitted species in EMIM-BF 4 , thus confirming that EMIM-BF 4 operates in the ion instead of in the droplet regime.The emission of monomers from EAN occurs much earlier, at 200 ps, at the extractor plane.In contrast to the EMIM-BF 4 simulation, the dimers and trimers are observed at approximately the same time, after 350 ps, rather than at later time steps.The distinct cone formation, shown in Figure 10c, allows a continuous emission of trimers, which are continually sampled at the extractor plane.It was hypothesized that dimers are the species in lowest concentration since EAN ion-pairs are either emitted as monomers from the sides of the Taylor cone or as larger ion-species from the tip of the Taylor cone.This was verified by generating velocity distributions of the monomers and droplets.A sampling plane located at 100 Å above the Taylor cone was selected because at that location, the emitted droplets and monomers do not undergo further acceleration and maintain their emitted velocities.This location also ensures that sufficient droplet statistics are obtained, since the droplets break down into smaller emission species as they travel towards the extractor plane.The tangential and axial velocity distributions of the droplets are shown in Figure 12a,b.Tangential velocities are generated because of the interaction of the radial electric field with the emitted ion species.The shape of the tangential velocity distribution of droplets has the form of a normal distribution with a mean value close to zero, suggesting that the majority of droplets emitted from the capillary have negligible tangential velocities.Since the radial electric fields are weakest near the tip of the capillary, Figure 12a can be interpreted to show that the droplets lack tangential velocities because they are emitted from the tip of the capillary.The distribution of droplet axial velocities, shown in Figure 12b, suggests that the majority of the droplets are emitted with an axial velocity between 40 and 80 km/s.In contrast, a similar analysis of monomer velocities does not give this result.The monomer tangential velocity distribution, shown in Figure 12c, does not have a maximum probability close to zero, but instead exhibits a noticeable number of samples at non-zero tangential velocities.This result suggests that emitted monomers interact with the radial electric field, and since the field is strongest near the sides of the Taylor cone, the majority are emitted from the sides of the Taylor cone, as was concluded by Lozano [22] in his experimental studies on the EMIM-BF 4 electrosprays.Finally, because the mass of monomers is comparatively much lower than that of droplets, they are strongly accelerated in the axial direction with no clear maxima, as shown in Figure 12d.
This type of information is vital to understanding the modes of operation of a colloid thruster for a particular ionic liquid.Furthermore, comparing the fraction of ion-to-droplet emission allows one to predict the efficiency of such a device.Unlike EMIM-BF 4 , the larger EAN ion-species, such as trimers and droplets, do not undergo further breakdown while traveling towards the extractor plane.The Coulomb interactions in the larger emitted species would require a much stronger electric field to break down the droplets into monomers and dimers.Therefore, for EAN, this results in trimers being the dominant emitted species, followed by monomers and droplets.In fact, the droplet emission current contributes approximately a quarter of the total emitted current.

Figure 2 .
Figure 2. Non-bonded interaction forces between CG sites E1-NN as a function of radial distance.(a) The total Coulomb force Equation (8) and inter-molecular Coulomb force Equation (10); (b) The van der Waals and the effective Equation (9) non-bonded force.

Figure 3 .
Figure 3. Probability distribution comparison of the sampled bond distance and angle data from the all-atom (AA) and coarse-grained (CG) MD simulations.(a) Bond distance E1-E2 from the coarse-grained model; (b) Bond distance E2-N1 from the coarse-grained model; (c) Angle N1-E1-E2 from the coarse-grained model.

Figure 5 .
Figure 5. Boundary conditions and isometric projection of the simulation domain used for the MD electrospray simulations.(a) Boundary conditions; (b) Isometric projection.

Figure 6 .
Figure 6.The electric fields generated by solving Equation (11) with the boundary conditions shown in Figure 5 and extraction voltage of −17 V. (a) Normal electric field; (b) Radial electric field.

Figure 7 .
Figure 7. Schematic of the simulation domain used for the bulk simulations to investigate cation-anion separation energy.

Figure 8 .
Figure 8.The Coulomb energy of 1000 ion pairs is shown, such that the solid and dashed lines represent the case when the entire cation-anion and when only the cation were moved towards the cathode, respectively.An extrusion electric field of 8 and 3 V/nm was assumed for EAN and EMIM-BF 4 , respectively.(a) EAN; (b) EMIM-BF 4 .

Figure 12 .
Figure 12.The tangential and axial velocity distributions of the monomers and droplets emitted from an EAN electrospray.The total number of droplets and monomers sampled is 390 and 321, respectively.(a) Tangential (V x and V y ) velocities of droplets; (b) Axial (V z , towards the extractor plane) velocities of droplets; (c) Tangential (V x and V y ) velocities of monomers; (d) Axial (V z , towards the extractor plane) velocities of monomers.

Table 1 .
Parameters for the CG EAN, shown in Figure1b, bonded interactions.