Ion emission from nanodroplets undergoing Coulomb explosions: a continuum numerical study

Abstract

A droplet charged above the Rayleigh limit is unstable. In the resulting dynamical process, referred to as a Coulomb explosion, smaller droplets with higher charge-to-mass ratios are ejected, reducing the charge of the parent droplet below the Rayleigh limit. Furthermore, if the droplet is sufficiently small, the electric field on its surface can promote ion field emission. Ion emission can lower the charge of a spherical droplet below its Rayleigh limit, keeping it stable, or reduce the charge of a deforming droplet, changing its dynamics and potentially preventing the Coulomb explosion. This article develops a continuum phase field electrohydrodynamic model to study the interplay between Coulomb explosions and ion emission, using charged nanodroplets of the ionic liquid 1-ethyl-3-methylimidazolium bis(trifluoromethylsulfonyl)imide (EMI-Im) as a case study. In small droplets (diameter $D \lesssim 20$ nm for EMI-Im), the electric field is strong enough to emit ions in the early phase of the droplet's evolution, suppressing the Coulomb explosion. For $20 \lesssim D \lesssim 45$ nm, the electric field on the EMI-Im droplet may not promote significant ion emission; however, as the unstable droplet becomes ellipsoidal, ions are emitted from its vertices, ultimately suppressing the Coulomb explosion while shedding $20\unicode{x2013}40\,\%$ of the initial charge. For larger EMI-Im droplets, $45 \lesssim D \lesssim 100$ nm, the evolution typical of a Coulomb explosion is observed, accompanied by ion emission which is however insufficient to prevent the Coulomb explosion. Ion emission and the smaller progeny droplets account for $24\,\%$ and $16\,\%$ of the initial charge, respectively.

1. Introduction

Electrosprays operating in the cone-jet mode (Zeleny 1917; Taylor 1964; Cloupeau & Prunet-Foch 1989; Fernández De La Mora & Loscertales 1994; Collins et al. 2008) have the unique ability for producing sprays of charged droplets with narrow size distributions, down to radii of a few nanometres (De Juan & Fernández De La Mora 1997; Rosell-Llompart, Grifoll & Loscertales 2018). The fate of these droplets is governed by the amount of charge they carry relative to their size which, when sufficiently high, can lead to the disintegration of the droplet and/or the shedding of charge by ion field emission. Rayleigh (1882) predicted that a spherical droplet charged beyond a critical value, referred to as the Rayleigh limit, is unstable and will disintegrate. This process is called the Rayleigh fission or Coulomb explosion of the droplet. The pressures inside and outside of the droplet are equal at the Rayleigh limit, i.e. the mechanical and electrical contributions to the surface tension cancel out one another (Rayleigh 1882; Taylor 1964; Basaran & Scriven 1989). This condition can be written in terms of the charge $Q_{Ray}$ carried by the droplet or the strength of the electric field $E_{Ray}$ on its surface:

(1.1a,b)\begin{equation} Q_{Ray}=8{\rm \pi}\sqrt{\varepsilon_{o} \gamma R^{3}}, \quad E_{Ray}=\sqrt{\frac{4\gamma}{\varepsilon_{o}R}}, \end{equation}

where $\gamma$ and $R$ are the surface tension and the radius of the droplet, respectively, and $\varepsilon _{o}$ is the permittivity of vacuum. A droplet charged at or above the Rayleigh limit will undergo a Coulomb explosion, shedding a fraction of its charge and mass in the form of smaller progeny droplets. Interested readers may refer to the photographs of Coulomb explosions taken by Gomez & Tang (1994), Duft et al. (2003) and Giglio et al. (2008).

The prediction of Rayleigh motivated several experimental studies. Doyle, Moffett & Vonnegut (1964) used Millikan's oil drop technique to study the Coulomb explosion of evaporating aniline and water droplets. They estimated a charge loss of roughly $30\,\%$ in the form of smaller progeny droplets, with very low mass loss. Schweizer & Hanson (1971) reported a charge and mass loss of $25\,\%$ and $5\,\%$, respectively. Duft et al. (2003) used an electrohydrodynamic levitation technique coupled with high-resolution time-elapse photography to isolate a charged droplet and capture the transient behaviour of the Coulomb explosion; they reported a charge loss of $33\,\%$ and a mass loss smaller than $1\,\%$. Giglio et al. (2008) conducted similar experiments and also performed numerical simulations. They also reported an experimental charge loss of $33\,\%$. Other studies report charge losses between $20\,\%\unicode{x2013}40\,\%$ (Roulleau & Desbois 1972; Taflin, Ward & Davis 1989), with the exception of sulphuric acid droplets which lose $50\,\%$ of their original charge (Richardson, Pigg & Hightower 1989). The images captured by Duft et al. (2003) and Giglio et al. (2008) have encouraged analytical and numerical computations of the transient behaviour of exploding droplets under different Reynolds numbers (Betelú et al. 2006; Burton & Taborek 2011; Collins et al. 2013; Radcliffe 2013; Garzon, Gray & Sethian 2014; Gawande, Mayya & Thaokar 2017, 2020). The experimental and numerical evidence show that a Coulomb explosion can be described as a three-step process: first, the slightly perturbed spherical droplet elongates into an ellipsoid of increasing aspect ratio, which eventually develops conical tips at the vertices; second, a fine jet is issued from each cusp and accelerated by the axial electric field; third, the jets, which are naturally unstable, break into progeny droplets which may undergo further Coulomb explosions depending on their electrification level.

Ion field emission is a second mechanism by which a droplet can shed charge. It is a kinetic process in which ions evaporating from the surface of a liquid must overcome an energy barrier (the ion solvation energy), which is lowered by the electric field. Early experiments demonstrated that ion field emission occurs in charged nanodroplets (Iribarne & Thomson 1976; Thomson & Iribarne 1979; Katta, Rockwood & Vestal 1991). Loscertales & Fernández De La Mora (1995) analysed the solid residues left behind by charged nanodroplets after complete evaporation of the liquid phase to quantify the electric field required for ion field emission. Labowsky, Fenn & Fernández De La Mora (2000) developed a continuum ion evaporation model that was in good agreement with the experimental findings of Gamero-Castaño & Fernández De La Mora (2000b) and Hogan & Fernández De La Mora (2009). These studies predict the electric field $E^{*}$ necessary for ion emission to be in the range $0.8\unicode{x2013}2\,{\rm V}\,{\rm nm}^{-1}$, which is only possible in highly charged droplets with diameters of tens of nanometres or smaller. These droplets can originate from much larger droplets in aerosols, where large residence times combined with solvent evaporation lead to a chain of Coulomb explosions producing increasingly smaller droplets. Alternatively, highly charged nanodroplets can be directly produced by electrospraying liquids with high electrical conductivities (Gamero-Castaño & Cisquella-Serra 2020; Miller et al. 2021; Perez-Lorenzo & Fernández De La Mora 2022).

The nanodroplets produced by electrospraying highly conducting liquids are often charged above the Rayleigh limit (Gamero-Castaño & Cisquella-Serra 2020; Miller et al. 2021; Perez-Lorenzo & Fernández De La Mora 2022). Moreover, Gamero-Castaño & Cisquella-Serra (2020) and Perez-Lorenzo & Fernández De La Mora (2022) observe that a significant fraction of the total current in these electrosprays is carried by ions, e.g. 20–25 % of the total current in electrosprays of EMI-Im and 1-ethyl-3-methylimidazolium tris(perfluoroethyl)trifluorophosphate (EMI-FAP). Both studies find that the ions must be evaporating from droplets in flight and propose that, in some cases, the emission occurs from droplets undergoing Coulomb explosions. The general consensus is that ions will evaporate from the surface of the droplet if the field emission limit is reached before the Rayleigh limit, i.e. if $E_{Ray} > E^{*}$ (Iribarne & Thomson 1976; Labowsky 1998; Labowsky et al. 2000). However, a droplet charged above the Rayleigh limit and having an electric field smaller than $E^{*}$ may develop areas with larger electric fields during the Coulomb explosion and emit ions. This shedding of charge may reduce substantially the net charge of the parent droplet, modify the dynamics of the Coulomb explosion and even prevent it.

The goal of this article is to study the interaction between a Coulomb explosion and ion field emission in highly charged nanodroplets. We develop an electrohydrodynamic, phase field model to compute the deformation of a droplet charged above the Rayleigh limit, while including ion field emission. We use this model to study EMI-Im nanodroplets in the diameter range 10–100 nm, because experiments show the simultaneous presence of Coulomb explosions and ion field emission under these conditions (Gamero-Castaño & Cisquella-Serra 2020). We address two specific questions: the extent to which ion emission suppresses Coulomb explosions and the magnitude of the charge emitted from typical nanodroplets. The remainder of this article is organized as follows: § 2 describes the numerical model; § 3.1 analyses the Coulomb explosion in the absence of ion emission, with the goals of setting a baseline and validating the numerical model; § 3.2 analyses the effects of ion emission in droplets of varying size charged above the Rayleigh limit; and § 4 summarizes the findings and recommends future work.

2. Problem formulation and numerical set-up

2.1. Problem formulation

Figure 1 shows the computational domain in cylindrical coordinates $\{r,z\}$. There are two subdomains, $\varOmega _{l}$ and $\varOmega _{g}$, occupied by an ionic liquid and ambient gas, respectively. The ionic liquid forms a spherical droplet of radius $R$ centred in the origin of coordinates, with a net charge exceeding the Rayleigh limit, $Q_o = \eta Q_{Ray}$ with $\eta > 1$. The droplet is slightly deformed to make it spheroidal while keeping its volume, such that the major axis (along $z$) is $1.04$ times the minor axis (Burton & Taborek 2011; Giglio et al. 2020). The droplet is then allowed to evolve from this starting configuration. The problem exhibits axial and planar ($z=0$) symmetries, and therefore the governing equations are solved in the quadrant bounded by segments $\varGamma _{b}$, $\varGamma _{r}$, $\varGamma _{u}$ and $\varGamma _{l}$. The viscosity, density, relative permittivity and electrical conductivity of the media are $\mu _i$, $\rho _i$, $\varepsilon _i$ and $K_i$, respectively, where the subscript $i$ denotes the ambient gas ($g$) or the ionic liquid ($l$).

Figure 1. Schematic of the problem and computational domain.

Electrohydrodynamic problems involving a free surface are usually solved using the leaky dielectric model (Melcher & Taylor 1969; Saville 1997), which treats the interface between the liquid and the surrounding gas as a surface. All excess charge is assumed to reside on the surface (surface charge). Conservation of charge is fulfilled by imposing a conservation equation on the surface, together with an Ohms law and the assumption of negligible volumetric charge in the bulk of the liquid. Drawbacks include the difficult calculation of the position of the free surface and the failure of model assumptions in ultra-fast liquid disintegration (Ganán-Calvo et al. 2016; Pillai et al. 2016). To avoid these problems, we use the phase field method (Anderson, Mcfadden & Wheeler 1998; Jacqmin 1999; Yue et al. 2004). This method regards the interface as a thin region of finite thickness, which can be easily tracked, and a continuous distribution of all independent variables throughout the media. It has the additional advantage of eliminating the assumption of negligible volumetric charge in the bulk of the liquid. The phase field method has been proven to be accurate in a variety of electrohydrodynamic problems (Tomar et al. 2007; López-Herrera, Popinet & Herrada 2011; López-Herrera et al. 2015; Mandal et al. 2015; Ganán-Calvo et al. 2016; Pillai et al. 2016).

We have previously developed an electrohydrodynamic phase field model for electrified jets and validated it with existing experimental and numerical results (Misra & Gamero-Castaño 2022). The ionic liquid and the surrounding gas are modelled as a continuum by defining a scalar phase variable $\phi$ that varies uniformly across the two media. Here, $\phi$ is tracked using the Cahn–Hilliard equations (Anderson et al. 1998), which is assigned values of $-1$ and $1$ in the gas and the liquid far from the interface, respectively. The surface separating the gas and liquid phases is defined as the loci where $\phi =0$. The physical properties are defined as continuous functions of $\phi$ throughout the media. Specifically, $\rho$, $\mu$ and $\varepsilon$ are the weighted arithmetic means of the phase variable:

(2.1a–c)$$\begin{gather} \rho=\rho_{g}\left(\frac{1-\phi}{2}\right)+\rho_{l}\left(\frac{1+\phi}{2}\right), \quad \mu=\mu_{g}\left(\frac{1-\phi}{2}\right)+\mu_{l}\left(\frac{1+\phi}{2}\right), \nonumber\\ \varepsilon=\varepsilon_{g}\left(\frac{1-\phi}{2}\right)+\varepsilon_{l}\left(\frac{1+\phi}{2}\right), \end{gather}$$

whereas, to mitigate unphysical charge leakage, the electrical conductivity is defined as

(2.2)\begin{equation} \frac{1}{K}=\frac{1}{K_{g}}\left(\frac{1-\phi}{2}\right)+ \frac{1}{K_{l}}\left(\frac{1+\phi}{2}\right) \end{equation}

(Tomar et al. 2007; López-Herrera et al. 2011; Roghair et al. 2015; Huh & Wirz 2022). Presently, the model of Misra & Gamero-Castaño (2022) is supplemented with the standard transport equation for ion field emission (Iribarne & Thomson 1976; Loscertales & Fernández De La Mora 1995):

(2.3)\begin{equation} J_{i}=\frac{k_{B}T}{h}\sigma \exp\left(-\frac{\Delta G - G(E_{n}^{g})}{k_{B}T}\right), \end{equation}

where $J_{i}$ is the ion current density evaporated from the surface, $k_B$ is the Boltzmann constant, $T$ is the temperature of the liquid, $\sigma$ is the surface charge density and $h$ is Planck's constant. Here, $\Delta G -G(E_{n}^{g})$ is the energy barrier that the evaporating ion must overcome. Additionally, $\Delta G$ is the ion solvation energy, which is not accurately known for ionic liquids although it is estimated to be in the range of 1.4–2 eV (Iribarne & Thomson 1976; Loscertales & Fernández De La Mora 1995; Gamero-Castaño & Fernández De La Mora 2000a,c). Also, $G(E_{n}^{g})$ is the reduction of the energy barrier due to the normal component of the electric field on the gas side of the surface $E_{n}^{g}$. We use the model by Iribarne & Thomson (1976) to evaluate the reduction of the energy barrier:

(2.4)\begin{equation} G(E_{n}^{g})=\left(\frac{e^{3}E_{n}^{g}}{4{\rm \pi}\varepsilon_{o}}\right)^{1/2}, \end{equation}

where $e$ stands for the elementary unit charge. Equation (2.3) suggests that a meaningful ion current can only occur when $\Delta G - G(E_{n}^{g})= {O}(k_{B}T)$. Since $k_{B}T \ll \Delta G$ at room conditions, the characteristic electric field for ion field emission is given by (Coffman et al. 2016; Coffman, Martínez-Sánchez & Lozano 2019; Gallud & Lozano 2022):

(2.5)\begin{equation} E^{*}=\frac{4{\rm \pi}\varepsilon_{o}{\Delta G}^{2}}{e^3}. \end{equation}

The dependent variables velocity $\boldsymbol {u}$, pressure $p$, volumetric charge density $\rho _e$, electric potential $V$ and the phase variable $\phi$ fulfil the equations of conservation of mass, momentum and charge, a modified Poisson's equation and the Cahn–Hillard equation throughout the computational domain:

(2.6)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}= 0, \end{gather}

(2.7)\begin{gather} \frac{\partial(\rho \boldsymbol{u})}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho \boldsymbol{u}\boldsymbol{u})={-}\boldsymbol{\nabla}p +\boldsymbol{\nabla} \boldsymbol{\cdot} (\mu (\boldsymbol{\nabla u}+\boldsymbol{\nabla u}^{\rm T}))+ \boldsymbol{F}_{es} + \boldsymbol{F}_{st}, \end{gather}

(2.8)\begin{gather} \frac{\partial\rho_e}{\partial t} +\boldsymbol{\nabla} \boldsymbol{\cdot} (\rho_e \boldsymbol{u})={-}\boldsymbol{\nabla} \boldsymbol{\cdot} (K\boldsymbol{E})-\frac{k_{B}T}{h}\rho_{e}\exp\left(-\frac{\Delta G - G(E_{n}^{g})}{k_{B}T}\right), \end{gather}

(2.9)\begin{gather} \varepsilon \nabla^{2}V + \boldsymbol{\nabla}V \boldsymbol{\cdot} \boldsymbol{\nabla}\varepsilon={-}\frac{\rho_e}{\varepsilon_{o}}, \end{gather}

(2.10a,b)\begin{gather} \frac{\partial \phi}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \phi=\zeta \gamma \nabla^{2} \psi , \quad \psi=\frac{1}{\xi}(\phi^{2}-1)\phi-\xi\nabla^{2}\phi . \end{gather}

Here, $\boldsymbol {F}_{es}$ and $\boldsymbol {F}_{st}$ are the electrostatic and surface tension volumetric forces:

(2.11)\begin{gather} \boldsymbol{F}_{es}= \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{T}_{e} = \boldsymbol{\nabla} \boldsymbol{\cdot} \varepsilon \varepsilon_o\left(\boldsymbol{EE}-\frac{1}{2}\boldsymbol{I}|E|^{2}\right) = \rho_e\boldsymbol{E}-\frac{1}{2}\varepsilon_o \boldsymbol{E} \boldsymbol{\cdot} \boldsymbol{E} \boldsymbol{\nabla}\varepsilon, \end{gather}

(2.12)\begin{gather} \boldsymbol{F}_{st}=\gamma\psi\boldsymbol{\nabla}\phi. \end{gather}

The volumetric charge conservation equation (2.8) can be derived from the general Poisson–Nernst–Planck (PNP) equation for charged species (Saville 1997; Herrada et al. 2012; Gañán-Calvo et al. 2018), see Appendix B for further details. The equation includes charge convection and conduction, and a term accounting for the loss of charge due to ion field emission. Equation (2.10a,b) is the Cahn–Hillard equation for the phase variable $\phi$ (Anderson et al. 1998; Jacqmin 1999, 2000; Yue et al. 2004); $\zeta$ is the phase field mobility parameter (a constant) and $\xi$ is the diffuse interface thickness which is indicative of the sharpness of the artificial interface separating the two phases. Brackbill, Kothe & Zemach (1992), Jacqmin (2000), Anderson et al. (1998), Jacqmin (1999) and Yue et al. (2004) provide derivations of the surface tension force, (2.12). It is worth noting that in the sharp interface limit, $\xi \rightarrow 0$, the phase field formulation transforms into the conventional treatment based on the use of a surface (Brackbill et al. 1992; Yue et al. 2004). We employ the following boundary conditions:

(2.13)\begin{gather} \boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{E} = 0, \quad \boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{u} = 0, \quad \boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{\nabla}\psi=\boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi = 0, \quad \text{on}\ \varGamma_{b}, \end{gather}

(2.14)\begin{gather} V=0, \quad \boldsymbol{u} = 0, \quad \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\nabla}\psi= \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi = 0, \quad \rho_{e}=0, \quad p=0, \quad \text{on}\ \varGamma_{r}, \end{gather}

(2.15)\begin{gather} V=0, \quad \boldsymbol{u} = 0, \quad \boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{\nabla}\psi= \boldsymbol{e}_{z} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi = 0, \quad \rho_{e}=0, \quad p=0, \quad \text{on} \varGamma_{u}, \end{gather}

(2.16)\begin{gather} \boldsymbol{e}_{r}\boldsymbol{\cdot} \boldsymbol{E}=0, \quad \frac{\partial (\boldsymbol{e}_{z}\boldsymbol{\cdot}\boldsymbol{u})}{\partial r}=\boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{u}= 0, \quad \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\nabla}\psi= \boldsymbol{e}_{r} \boldsymbol{\cdot} \boldsymbol{\nabla}\phi = 0 ,\quad \text{on}\ \varGamma_{l}. \end{gather}

We write the system of equations in dimensionless form using $l_c = R$, $t_c= \mu _{l}l_c/\gamma$, $u_c = l_c / t_c$, $p_c = \gamma / l_{c}$, $E_c = Q_{o}/(4 {\rm \pi}\varepsilon _{o} l_{c}^{2})$ and $\rho _{e,c}=\varepsilon _{o} E_c / l_{c}$ as the characteristic scales for length, time, velocity, pressure, electric field and volumetric charge, respectively. In particular, the dimensionless forms of the governing equations (2.6)–(2.10a,b) read

(2.17)\begin{gather} \boldsymbol{\widetilde{\nabla}} \boldsymbol{\cdot} \boldsymbol{\tilde{u}}= 0, \end{gather}

(2.18)\begin{gather} \frac{1}{{{Oh}^2}}\frac{\partial(\frac{\rho}{\rho_l} \boldsymbol{\tilde{u}})}{\partial \tilde{t}}+ \frac{1}{{{Oh}^2}}\boldsymbol{\widetilde{\nabla}} \boldsymbol{\cdot} \left(\frac{\rho}{\rho_l} \boldsymbol{\tilde{u}}\boldsymbol{\tilde{u}}\right)={-}\boldsymbol{\widetilde{\nabla}}\tilde{p} +\boldsymbol{\widetilde{\nabla}} \boldsymbol{\cdot} \left(\frac{\mu}{\mu_l} \left(\boldsymbol{\widetilde{\nabla} \tilde{u}}+\boldsymbol{\widetilde{\nabla} \tilde{u}}^{\rm T}\right)\right)+ \varGamma \boldsymbol{\tilde{F}}_{es} + \boldsymbol{\tilde{F}}_{st}, \end{gather}

(2.19)\begin{gather} {\frac{\varPi_{1}}{\varepsilon_l}} \frac{\partial\tilde{\rho}_e}{\partial \tilde{t}} + {\frac{\varPi_{1}}{\varepsilon_l}} \boldsymbol{\widetilde{\nabla}} \boldsymbol{\cdot} \left(\tilde{\rho}_e \boldsymbol{\tilde{u}}\right)={-}\boldsymbol{\widetilde{\nabla}} \boldsymbol{\cdot} \left(\frac{K}{K_l}\boldsymbol{\tilde{E}}\right)-{\frac{\varPi_{2}}{\varepsilon_l}}\tilde{\rho}_{e} \exp\left(-\varPi_{3}\left(1- \sqrt{\tilde{E}_{n}^{g}/ \varPi_{4}} \right)\right), \end{gather}

(2.20)\begin{gather} \varepsilon \widetilde{\nabla}^{2}\tilde{V} + \boldsymbol{\widetilde{\nabla}}\tilde{V} \boldsymbol{\cdot} \boldsymbol{\widetilde{\nabla}}\varepsilon={-}\tilde{\rho}_e, \end{gather}

(2.21a,b)\begin{gather} \frac{\partial \phi}{\partial \tilde{t}} + \boldsymbol{\tilde{u}} \boldsymbol{\cdot} \boldsymbol{\widetilde{\nabla}} \phi= \frac{1}{{Pe}}\widetilde{\nabla}^{2} \tilde{\psi} , \quad \tilde{\psi}=\frac{1}{\tilde{\xi}}(\phi^{2}-1)\phi-\tilde{\xi}\widetilde{\nabla}^{2}\phi, \end{gather}

where dimensionless variables are designated with an overtilde. Equations (2.17)–(2.21a,b) include eight dimensionless numbers: ${Oh}$, $\varGamma$, $\varPi _1$, $\varPi _2$, $\varPi _3$, $\varPi _4$, $Pe$ and $\varepsilon _l$. In addition, the ratios $\rho _g / \rho _l$, $\mu _g/\mu _l$ and $K_g/K_l$ are nearly zero and $\varepsilon _g \cong 1$ for all liquid/gas combinations. The Ohnesorge number is the ratio between the viscous time scale $t_c$ and the inertial time scale $\sqrt {\rho _{l} R^{3}/\gamma }$:

(2.22)\begin{equation} {Oh}=\frac{\mu_l}{\sqrt{\gamma \rho_{l}R}}, \end{equation}

and measures the relative importance of viscous and inertial forces. The Taylor number:

(2.23)\begin{equation} \varGamma= \frac{\varepsilon_{o}E_{c}^{2}R}{\gamma} \end{equation}

measures the relative importance between the electrostatic and capillary stresses. In particular, $\varGamma = 4$ indicates that the droplet is charged at the Rayleigh limit. The $\varPi _1$ is the ratio between the electrical relaxation time $t_e = \varepsilon _l \varepsilon _o/K_l$ and the characteristic time scale:

(2.24)\begin{equation} \varPi_1 = \frac{t_e}{t_c}= \frac{ \varepsilon_l \varepsilon_o \gamma }{K_l\mu_{l} R}. \end{equation}

The $\varPi _1$ is indicative of the speed at which the charge migrates from the bulk to the surface in an attempt to make the liquid phase equipotential. The dimensionless groups $\varPi _2$, $\varPi _3$ and $\varPi _4$ determine the importance of ion field emission:

(2.25a–c)\begin{equation} \varPi_2 = \frac{\varepsilon_l \varepsilon_o k_{B}T }{ K_l h}, \quad \varPi_3 = \frac{\Delta G }{k_{B}T}, \quad \varPi_4 = \frac{4{\rm \pi}\varepsilon_{o}{\Delta G}^{2}}{e^3 E_{c}}. \end{equation}

Here, $\varPi _2$ is the ratio between the electrical relaxation time and the molecular evaporation time $h/k_B T$, $\varPi _3$ is the ratio between the ion solvation energy and the molecular thermal energy and $\varPi _4$ is the dimensionless characteristic electric field for ion emission, $\tilde {E}^{*}$. Finally, the Péclet number measures the relative importance of convection and diffusion in the Cahn–Hilliard equation:

(2.26)\begin{equation} {Pe}=\frac{R^3}{\zeta\gamma t_c}. \end{equation}

For all the numerical cases considered in this study, we make $\zeta =(R\xi )^{2}/p_{c}t_c$ (Ding, Gilani & Spelt 2010; Mandal et al. 2015), i.e. ${Pe}=1/\tilde {\xi }^2$.

2.2. Numerical implementation

Our goal is to study the role of ion field evaporation in the formation of electrosprayed nanodroplets of highly conducting liquids, in particular ionic liquids. To the best of our knowledge, this is the only practical situation in which droplets directly form both above the Rayleigh limit and can field emit ions. The central problem is understanding how ion field emission modifies the sprays of a given liquid at varying flow rate, or equivalently as a function of the droplet radius. We focus on the electrosprays of the ionic liquid EMI-Im because it has been thoroughly studied, e.g. it is known that ion emitted from droplets in flight is a key element of the physics, and that some droplets undergo Coulomb explosions while others emit ions and do not fragment (Gamero-Castaño & Cisquella-Serra 2020). Thus, in all calculations, we use the physical properties of EMI-Im, $\varepsilon _{l}=12.2$, $\rho _{l}=1520\,{\rm kg}\,{\rm m}^{-3}$, $\mu _{l}=0.032\,{\rm Pa}\,{\rm s}$, $K_{l}=0.74\,{\rm S}\,{\rm m}^{-1}$ and $\gamma =0.0349\,{\rm N}\,{\rm m}^{-1}$. For the gas phase, we use $\varepsilon _{g}=1$, $\rho _{g}=1\,{\rm kg}\,{\rm m}^{-3}$, $\mu _{g}= 10^{-4}\,{\rm Pa}\,{\rm s}$ and $K_{g}=10^{-9}\,{\rm S}\,{\rm m}^{-1}$. A precise value of the ion solvation energy is not available and we use the estimate $\Delta G = 1.62$ eV (Iribarne & Thomson 1976; Loscertales & Fernández De La Mora 1995; Gamero-Castaño & Fernández De La Mora 2000a). In all simulations, we consider an initial charge $4\,\%$ above the Rayleigh limit, $Q_o=1.04Q_{Ray}$.

Table 1 shows typical flow rates, beam currents and droplet radii of EMI-Im electrosprays (Gamero-Castaño & Cisquella-Serra 2020), together with the values of size-dependent dimensionless numbers (we assume a charge $4\,\%$ above the Rayleigh limit). The values of the dimensionless numbers that do not depend on the radius of the droplet are $\varGamma = 4.016$, $\varPi _2 = 906$, $\varPi _3 = 63.1$, $\varepsilon _{l}=12.2$ and ${Pe} = 10^4$ (we use $\xi = R/100$ in all calculations). Note that ion field emission from droplets of highly conducting liquids is characterized by high Ohnesorge numbers (all ionic liquids have high viscosities, comparable or larger to that of EMI-Im, and the droplet radii must be at most a few tens of nanometres to sustain the electric fields needed for significant ion emission) and small $\varPi _1$ (the charge is, for all purposes, relaxed on the surface). Therefore, varying the droplet radius for a given ionic liquid at constant excess charge over the Rayleigh limit is equivalent to varying $\varPi _4$ while keeping $\varGamma$, $\varPi _2$, $\varPi _3$, ${Pe}$ and $\varepsilon _{l}$ constant, in the limit ${Oh}^2 \gg 1$ and $\varPi _1 \ll 1$.

Table 1. Flow rates, beam currents and droplet radii of EMI-Im electrosprays (Gamero-Castaño & Cisquella-Serra 2020), together with size-dependent dimensionless numbers for droplets charged $4\,\%$ above the Rayleigh limit.

We use COMSOL Multiphysics Software to solve the system of equations (Comsol, Inc. 2019). We employ for most equations COMSOL's built-in laminar flow, electrostatics and phase field interface, which uses a finite element solver. The volumetric charge conservation equation (2.8) cannot be incorporated with built-in interfaces so we use instead the partial differential equation interface. Additionally, we include the electrostatic (2.11) and surface tension (2.12) volumetric forces in (2.7) as forcing terms. The equations are solved in COMSOL's weak formulation framework. The phase variable $\phi$ is discretized using a cubic-order Lagrange element; $\boldsymbol {u}$, $V$ and $\rho _{e}$ are discretized using quadratic-order Lagrange elements; and $p$ is discretized using a linear Lagrange element. We use the parallel sparse solver MUMPS for marching the solution in time. MUMPS uses a second-order backward differential formulation scheme with variable time step, computed using the Courant–Friedrichs–Lewy condition (Courant, Friedrichs & Lewy 1967). The maximum time step is set at $\Delta t_{max}/t_{c}=0.01$, while typical time steps at the start of the simulation are of the order of $10^{-5}$. When needed (e.g. to analyse the results), the position of the surface is computed as the loci where $\phi =0$ by interpolation. For spatial discretization, we use a triangular grid with grid size $h$. The computational domain is divided into two regions, as depicted in figure 2. Region A contains the droplet and its grid size is $h=R/64$, while region B has a grid size $h=R/30$. The variation of the phase variable $\phi$ is depicted in figure 2 along with details of the triangular grid. The diffuse interface thickness is fixed at $\xi =R/100$ (Misra & Gamero-Castaño 2022). The post-processing of the results is performed with the COMSOL–MATLAB interface.

Figure 2. A section of the computational domain along with details of the grid and the variation of the phase variable $\phi$.

Initially, we impose a homogeneous charge density $\rho _{e,o}=3 \eta Q_{Ray}/(4 {\rm \pi}R^{3})$ in $\varOmega _l$, with $\eta =1.04$. The charge is allowed to relax at zero fluid velocity, migrating to the interface to make the droplet equipotential (Misra & Gamero-Castaño 2022). Once the charge relaxes, we compute the deformation of the droplet over time. In most calculations, we use $H=L=10R$ and have verified that these values do not affect the solution. When calculating the elongation of the jet, we use a longer domain with $H=14R$. Appendix B provides details about the charge relaxation step and Appendix C presents a grid independence study.

3. Results and discussions

3.1. The Coulomb explosion of a droplet in the absence of ion field emission

This problem has been thoroughly investigated (Burton & Taborek 2011; Collins et al. 2013; Garzon et al. 2014; Gawande et al. 2017) and the comparison with our solution makes it possible to validate the numerical model. In the remainder of the paper, this case is referred to as ‘pure Rayleigh fission’, PRF. Physically, this scenario corresponds to large droplets with electric fields substantially lower than $E^{*}$. For the calculations, we use a droplet with $D=50$ nm having the physical properties of EMI-Im, with the exception of using a lower electrical conductivity, $K=0.095\,{\rm S}\,{\rm m}^{-1}$, to simplify capturing the evolution. This leads to $Oh\sim 28$ and $\varPi _{1}\sim 0.05$.

Figure 3(a) shows the evolution of a droplet with $Q_{o}=1.04Q_{Ray}$. Initially nearly spherical, the droplet becomes an ellipsoid of increasing aspect ratio. The vertices eventually transition into cusps, which elongate to form jets. We quantify the evolution of the droplet with the aspect ratio, ${\it AR}$, defined as the quotient between the spans of the droplet along the axial and radial directions. Figure 3(b) shows the aspect ratio and the electric field on the vertices $\tilde {E}_{tip}$, as a function of time. The shape of the droplet just prior to the formation of jets is referred to as the critical shape. Its aspect ratio, ${\it AR} \cong 3.82$, is close to the experimental value of $3.7\unicode{x2013}3.85$ for droplets in the Stokes limit (Duft et al. 2003; Achtzehn et al. 2005; Giglio et al. 2008); it also compares well with the values of previous numerical calculations, $3.85\unicode{x2013}3.87$ (Gawande et al. 2017, 2020; Giglio et al. 2020). Here, $\tilde {E}_{tip}$ reaches its maximum value for an aspect ratio slightly larger than that of the critical shape. The aspect ratio grows rapidly once the jet is formed. Appendix A provides additional comparison with experimental data.

Figure 3. Coulomb explosion of a droplet without ion emission (PRF), with initial charge $1.04Q_{Ray}$, $Oh=28$ and $\varPi _{1}=0.05$: (a) evolution of the droplet's shape with time. The aspect ratios (${\it AR}$) are $1.36$, $1.89$, $2.43$, $3.20$, $3.79$, $4.79$ and $8.41$ at increasing time; (b) evolution of the aspect ratio and the electric field on vertex $\tilde {E}_{tip}$; (c) evolution of the droplet starting right before the critical shape is formed (${\it AR}=3.68$), when the total charge is reduced to (c i) $Q_{o}=0.64Q_{Ray}$ and (c ii) $Q_{o}=0.65Q_{Ray}$.

Prior experiments and calculations have shown that the jets are unstable and generate small droplets. The main elongated drop left after the detachment of the progenies relaxes back to a spherical shape, with a charge that is substantially smaller than the original and most of the original mass (Duft et al. 2003; Giglio et al. 2008). An accurate calculation of the dynamics of the jet breakup and the formation of droplet progenies requires a higher grid resolution than employed in this article (the capture of radial features is limited by our choice of interface thickness, $\xi =R/100$), but doing so increases the computational cost dramatically. Although we do not capture the entire fission process, we can estimate the amount of charge that a parent droplet must shed to become stable. Starting with the shape at ${\it AR}=3.68$ (just prior to the formation of the critical shape), we progressively reduce the charge and observe the evolution to discern whether the droplet develops the jets or goes back to the spherical shape. Figure 3(c) illustrates this approach: starting with a shape ${\it AR}=3.68$, if the charge is made to be $Q_{o}=0.65Q_{Ray}$, the droplet still forms a jet and is unstable. Conversely, if the charge is slightly reduced to $Q_{o}=0.64Q_{Ray}$, the aspect ratio of the droplet decreases and it becomes spheroidal. This suggests that a droplet charged at the Rayleigh limit will shed $36\,\%$ of its charge in a Coulomb explosion. Duft et al. (2003) and Giglio et al. (2008) report an experimental charge loss of $33\,\%$, which is in good agreement with our prediction. A recent equipotential model predicts a charge loss of roughly $39\,\%$ (Gawande et al. 2017). There is a variation of the fraction of lost charge reported in the literature, which is likely due to viscosity and conductivity effects. For example, for low viscous de-ionized water droplets ($Oh=0.023$), the critical shape is achieved at ${\it AR}=2.7\unicode{x2013}2.9$ and a charge loss of roughly $20\,\%$ has been measured (Giglio et al. 2020), whereas experiments in the Stokes limit ($Oh\gg 1$) report a charge loss of $33\,\%$ (Duft et al. 2003; Achtzehn et al. 2005; Giglio et al. 2008).

3.2. Droplet charged above the Rayleigh limit with ion field emission

The simulation of a droplet charged above the Rayleigh limit including ion field emission (labelled as WFE, ‘with field emission’) is motivated by the experimental studies of Gamero-Castaño & Cisquella-Serra (2020), Miller et al. (2021) and Perez-Lorenzo & Fernández De La Mora (2022), who report substantial ion currents in electrosprays of several ionic liquids. Retarding potential measurements indicate that these ions are emitted from droplets likely to be charged above the Rayleigh limit. Our simulations reveal three distinct outcomes depending on the initial size of the droplet.

3.2.1. Small radius regime, region $I$

Sufficiently small droplets charged above the Rayleigh limit emit ions from the spherical and spheroidal shapes. Figure 4 shows the evolution of droplets with diameters of 10 and 20 nm charged $4\,\%$ above the Rayleigh limit. The 10 nm droplet undergoes marginal deformation while emitting ions, before becoming a stable sphere. The $20$ nm droplet follows a similar path, but undergoes a larger deformation. In either case, the droplet does not develop the cusps preceding the formation of the jets and droplet progenies, i.e. the shedding of charge by ion field emission prevents the Coulomb explosion.

Figure 4. Evolution of EMI-Im droplets emitting ions, with an initial charge $Q_{o}= 1.04Q_{Ray}$: (a) $D = 10\,{\rm nm}$; (b) $D = 20\,{\rm nm}$. The times $t/t_{c}$ for each snapshot are shown above the drops.

Figure 5(a) shows the charge of three ion-emitting droplets with diameters of $10$, $15$ and $20$ nm as a function of time, while figure 5(b) shows their aspect ratios together with that of a $50$ nm droplet in the PRF regime (without ion emission). The initial charge of all droplets is $4\,\%$ above the Rayleigh limit. The characteristic time for the shedding of charge scales with $t_c$, while the fraction of the initial charge that is evaporated increases with decreasing droplet diameter. This is to be expected from the dependency of the electric field on the diameter of a droplet at the Rayleigh limit, (1.1a,b). Note that as the droplets emit ions, their aspect ratios increase. The droplets remain quasi-equipotential during the deformation, and the electric field on the surface scales with the inverse of the square of the local radius of curvature. Thus, as the aspect ratio of the droplet increases, ion emission preferentially takes place from the vertices of the ellipsoid. The aspect ratios of the droplets plateau once ion emission becomes negligible, a larger droplet exhibits a larger deformation. From this state of maximum deformation, the droplet falls back to the spherical shape, now stable with a constant charge that can be significantly smaller than the Rayleigh limit. The 10, 15 and 20 nm droplets lose $38\,\%$, $28\,\%$ and $20\,\%$ of their initial charge, respectively. Given the substantial charge loss, the electric stress at the vertices of the ellipsoid is insufficient to form the cusps and jets characteristic of a Coulomb explosion.

Figure 5. EMI-Im droplets with different diameters emitting ions (initial charge $Q_{o}= 1.04Q_{Ray}$): (a) evolution of the charge held by the droplet; (b) evolution of the aspect ratio and comparison with that of a droplet in the PRF regime (without ion emission).

Figure 6 shows the evolution of the normal component of the electric field at the vertex of the droplet. In all cases, the electric field levels off near $E_{n}^{g}/E^{*} \approx 0.76$ when ion emission is significant, i.e. $0.76\times E^{*}$ is the effective value of the electric field during ion emission: when the electric field is slightly above $0.76\times E^{*}$, the high intensity of the ion flux (see exponential dependency of the emitted current on the electric field) rapidly reduces the charge of the droplet, bringing down the electric field to this effective value for ion emission; and ion emission is negligible if the electric field is slightly below $0.76\times E^{*}$. The initial electric field for the two smaller droplets is above the effective value and therefore both of them are characterized by an initial brief period of intense ion emission. Conversely, the electric field for the larger droplet, $D = 20$ nm, is initially below the effective value and, although the droplet emits ions from its initial state, it needs to deform to reach the effective electric field and increase ion emission. The total charge evaporated during the time in which the electric field is approximately constant is significant for all three droplets (see figure 5a). The freezing of the electric field during ion emission was previously observed in the experiments by Loscertales & Fernández De La Mora (1995). Finally, note that the evolution of the electric field in figure 6 may exhibit, at times, an unphysical wiggle. This artefact is due to how we evaluate the electric field for the purpose of making the plot: $E_{n}^{g}$ is defined at the surface, which in general does not coincide with the position of the nodes in the computational grid. Thus, we compute $E_{n}^{g}$ by interpolating the values of the electric field between two nodes, in a region where the electric field and the volumetric charge exhibit large gradients. This causes the wiggle. A more accurate way of computing $E_{n}^{g}$ would consist of using more nodes in the grid to fit the values of the electric field to appropriate analytical functions and determining the value at the surface from the fitting. However, note that this is only important for plotting purposes, because $E_{n}^{g}$ is not a variable used in the calculations (the surface and therefore the electric field at the surface are not part of the calculations in the phase field method).

Figure 6. Evolution of the outer electric field on the vertices of EMI-Im droplets emitting ions, with an initial charge $Q_{o} = 1.04Q_{Ray}$.

3.2.2. Medium radius regime, region $II$

In this case, the initial electric field on the surface of the spherical droplet is insufficient to emit ions. Similarly to the PRF case, the unstable droplet deforms into an ellipsoid of increasing aspect ratio at constant charge. As the electric field near the vertices of the ellipsoid increases, it eventually becomes high enough to promote significant ion emission, discharging the droplet and preventing the Coulomb explosion.

Figure 7 displays the evolution of a $D = 30$ nm droplet initially charged $4\,\%$ above the Rayleigh limit. Figure 7(a) shows the droplet deforming into an ellipsoid of increasing aspect ratio, reaching a maximum deformation, and going back to a sphere without developing the cusps and jets typical of a Coulomb explosion. Figure 7(b) shows the aspect ratio, the charge and the electric field at the vertices of the droplet as functions of time. Initially, the trend for the aspect ratio coincides with that for the PRF case because the electric field near the vertices is insufficient to support ion emission. However, as the electric field at the vertices reach a critical value at $t/t_{c}\sim 10$, the droplet begins to shed charge and the aspect ratios of the PRF and WFE droplets start to separate. The aspect ratio of the WFE droplet, rather than accelerating, reaches a maximum at $t/t_{c}\sim 38$. By this time, the droplet has lost $17.3\,\%$ of its charge by ion field emission, which takes place at nearly constant electric field, $E_{n}^{g}/E^{*} \sim 0.76\unicode{x2013}0.77$. Most of the ion emission occurs up to the maximum of the aspect ratio. Minor ion emission continues to take place as the droplet relaxes to the stable spherical configuration, being negligible for $t/t_{c}\gtrsim 54$. Because of the reduced electrification, the cusps and jets never develop, i.e. the Coulomb explosion is suppressed.

Figure 7. Results for a $D = 30$ nm droplet with initial charge $Q_{o}=1.04Q_{Ray}$, including ion emission: (a) evolution of the droplet's shape ($t/t_{c}$ is shown above the drops); (b) evolution of the aspect ratio (for comparison, we also include the aspect ratio of a PRF droplet), charge and outer electric field at the vertices of the droplet.

Figure 8 shows the evolution for a larger droplet, $D=40$ nm, also charged $4\,\%$ above the Rayleigh limit. Although, in this case, ion emission is triggered at a later time due to the lower initial electric field, it also suppresses the Coulomb explosion. The droplet follows the PRF trend for a longer time due to the later onset of ion emission. At $t/t_{c}\sim 14$, the electric field at the vertices is sufficiently high to trigger ion emission ($E_{n}^{g}/E^{*}\sim 0.76\unicode{x2013}0.77$). The droplet starts to emit charge and the aspect ratio separates from the PRF droplet. Note that the vertices of the ellipsoid are almost transitioning into cusps when the aspect ratio plateaus at $t/t_{c}\sim 38.2$. Interestingly, the aspect ratio exceeds $3.8$, which is the critical aspect ratio for the PRF case. The larger 40 nm droplet sheds a higher fraction of its initial charge than the $30$ nm droplet, $\sim 24\,\%$ and $\sim 19\,\%$, respectively.

Figure 8. Same as figure 7 for a larger droplet, $D = 40$ nm.

These findings (regions I and II) provide an alternative scenario for ion emission to that proposed by previous authors who assume that it takes place from a spherical droplet only when its initial electric field $E_o$ can promote significant ion emission while being lower than that associated with the Rayleigh limit, $E^{*}\leq E_o< E_{Ray}$ (Iribarne & Thomson 1976; Labowsky 1998; Labowsky et al. 2000). Furthermore, they confirm the qualitative arguments put forward by Gamero-Castaño & Fernández de la Mora (2000) and Gamero-Castaño & Cisquella-Serra (2020), namely that droplets undergoing Coulomb explosions may emit ions from areas of high electric fields, suppressing the fission process.

3.2.3. Large radius regime, region $III$

Lastly, we discuss relatively large droplets charged above the Rayleigh limit, specifically in the diameter range $45\unicode{x2013}100$ nm. The unstable droplet evolves into an ellipsoid, developing cusps and jets. Ion emission may take place from the vertices and the jets, but cannot prevent the Coulomb explosion. Both ion emission and droplet progenies reduce the charge of the droplet.

Figure 9 shows a simulation for a $D = 45$ nm droplet initially charged $4\,\%$ above the Rayleigh limit. The evolution of the shape is qualitatively similar to the PRF case: the spherical droplet becomes an ellipsoid of increasing aspect ratio; it develops cusps at the vertices; and the cusps move away from the centre forming jet-like extensions. Ion emission is negligible until $t/t_{c}=15.4$. Beyond this time, the electric field at the vertices is high enough to emit ions, lowering the charge of the droplet and slowing the speed at which the aspect ratio increases. However, the shedding of charge does not impede the formation of cusps and jets, depicted in the droplet profiles at $t/t_{c}= 33.9, 37.8$ and 40.7. It is interesting that electric fields high enough to evaporate ions only develop at the vertices of the spheroid and at the cusps. At $t/t_{c}\gtrsim 38$, the electric field at the cusps is not sufficient to continue the emission of ions and the net charge of the droplet remains constant. The jets continue to elongate, and we assume that natural instability will lead to their breakup and the formation of droplet progenies. Overall, the droplet sheds roughly $28\unicode{x2013}29\,\%$ of its charge by ion emission. Profiles of the electric field along the surface at different times are shown in figure 9(c), $t/t_{c}=35.3$, and figure 9(d), $t/t_{c}=41.2$. At the earlier time, the electric field increases towards the cusps, where it reaches a value consistent with significant ion evaporation, $E_{n}^{g}/E^{*}=0.76\unicode{x2013}0.77$. The profile for $t/t_{c}=41.2$, i.e. for a situation with developed jets, shows two distinct maxima. The lower maximum occurs near the transition region of the jet $z/R\sim 2.6\unicode{x2013}2.7$, whereas the larger electric field occurs at the tip of the jet. The value of the latter is $E_{n}^{g}/E^{*} \sim 0.67$, i.e. insufficient to trigger significant ion emission, as reflected by the constancy of the charge carried by the droplet at this time.

Figure 9. Simulation of a $D = 45$ nm droplet with initial charge $Q_{o}=1.04Q_{Ray}$, undergoing Coulomb explosion and ion field emission (WFE): (a) droplet shape at different times; (b) evolution of the charge and aspect ratio of the droplet; (c) profile of the droplet's surface and outer electric field at $t/t_{c}=35.3$; (d) profile of the droplet's surface and outer electric field at $t/t_{c}=41.2$.

Figure 10 shows that larger droplets (diameters of 50, 80 and 100 nm) behave similarly to the droplet with a diameter of 45 nm: ion emission is absent from the initial spherical configuration; the deformation follows the PRF path until the electric field at the vertices is sufficient to evaporate ions; and ion emission does not prevent the formation of jets and the Coulomb explosion. The aspect ratio curves match the PRF trend for most of the evolution; however, the three curves depart from the PRF trend before the critical shape forms. This is due to the onset of ion emission from the vertices before the formation of cusps. Larger droplets follow the PRF path more closely and we expect that ion emission never occurs for sufficiently large droplets. The three droplets loose a similar fraction of their charge by ion emission, approximately $24\,\%$, but the loss is more gradual and starts at an earlier time for the smaller droplet. As in the case of the $D = 45$ nm droplet, ion emission eventually ceases at a time when the jet is still elongating, suggesting that the droplets will undergo Coulomb explosions and produce droplet progenies.

Figure 10. Results for droplets with diameters of 50, 80 and 100 nm: (a) evolution of the aspect ratio and comparison with that of a PRF droplet; (b) evolution of the droplet's charge.

In our calculations, the large droplets that emit ions during a Coulomb explosion lose approximately $24\,\%$ of their charge by ion emission. Since the total charge loss of PRF droplets is approximately $36\,\%$, we estimate that the droplets in region III lose approximately $16\,\%$ of the charge in the form of droplet progenies (considering that $Q_{o}=1.04Q_{Ray}$). The rationale behind this estimate is that the shape of the ion-emitting droplet at the time the jets form is similar to the critical shape of the PRF droplet and therefore they both must lose a similar amount of charge to fall back to the stable spherical configuration (this is the argument made by figure 3c). The comparison of the droplet shapes shown in figure 11 supports this estimate: the shape of the PRF droplet right before the critical point (${\it AR} = 3.68$, see figure 3a) matches the profile of the $D = 100$ nm droplet well, with slight differences near the vertices. Our estimate that large droplets lose approximately $24\,\%$ of their charge by ion emission agrees with the findings of Perez-Lorenzo & Fernández De La Mora (2022), who hypothesize that the majority of ions observed in the electrosprays of EMI-FAP are emitted from nanodroplets undergoing Coulomb explosions. The experiments of Perez-Lorenzo & Fernández De La Mora (2022) show that approximately $20\unicode{x2013}25\,\%$ of the total current is emitted in the form of ions, a value not far from our prediction; furthermore, the average diameter of the primary droplets of EMI-FAP in these experiments is estimated to be near 90 nm, i.e. within regime III.

Figure 11. Profiles of a $D = 100$ nm droplet undergoing Coulomb explosion with ion emission (${\it AR}=3.80$ and ${\it AR}=7.10$) and comparison with a PRF droplet near the critical shape (${\it AR}=3.68$).

Figure 12(a) shows the total charge emitted from EMI-Im droplets in regions I and II by ion emission as a function of the diameter of the droplet. Figure 12(b) plots the electric field on the surface of a spherical droplet charged 4 % above the Rayleigh limit, as a function of the diameter of the droplet (the electric field is normalized with the characteristic value for ion emission $E^{*}$); figure 12(b) also shows the domains of regions I, II and III. Here, $E_{c}/E^{*} \gtrsim 0.77$ in region I and, therefore, ion emission is significant from the initial spherical configuration. Although the droplet deforms into an ellipsoid of increasing aspect ratio, ion emission prevents the formation of cusps and the Coulomb explosion and, once enough charge is shed, the droplet goes back to a spherical shape with a charge below the Rayleigh limit. Region II is characterized by droplets with initial electric fields in the range $0.47 \lesssim E_{c}/E^{*} \lesssim 0.77$. The initial electric field is not sufficient to trigger ion emission; however, as the droplet deforms and the curvature at the vertices becomes large enough, ions are emitted from this small region, discharging the droplet and preventing the formation of cusps and the Coulomb explosion. The droplet subsequently retracts to the stable spherical configuration below the Rayleigh limit. Region III is characterized by initial electric fields $E_{c}/E^{*} \lesssim 0.47$. These droplets develop cusps and jets, and hence undergo Coulomb explosions yielding droplet progenies. These droplets lose a significant fraction of their charge by ion emission from their vertices and jets; however, the reduction of charge does not prevent the Coulomb explosion and the formation of droplet progenies.

Figure 12. (a) Charge lost by ion emission ($Q_{i}$) for droplets in regions I and II; (b) electric field on the surface of a spherical droplet charged 4 % above the Rayleigh limit, as a function of the diameter of the droplet and ranges of regions I, II and III.

We can estimate the maximum diameter of the droplet, $D_{max}$, for which ion emission is significant, i.e. the diameter separating region III from the PRF regime. Figure 3(b) shows that, for EMI-Im droplets in the PRF regime charged near the Rayleigh limit, the maximum value of the electric field, $E_{tip}/E_{Ray} \cong 2.83$, occurs at the cusps when the aspect ratio is slightly larger than that of the critical shape. Furthermore, ion emission is significant only when $E_{max}/E^{*} \geq 0.76$. These two relations can be combined to obtain $D_{max}$ (considering droplets $4\,\%$ above the Rayleigh limit):

(3.1)\begin{equation} D_{max}=\frac{\gamma}{2 \varepsilon_{o}^3} \left( \frac{e^3}{0.26 {\rm \pi}\Delta G^2} \right) ^2. \end{equation}

Here, $D_{max} = 140$ nm for the conditions used in our numerical calculations. We expect droplets with larger diameters to be in the PRF regime.

We have used an initial charge $4\,\%$ over the Rayleigh limit or equivalently, $\varGamma = 4.016$, in all simulations. A higher initial charge excess should have a minor effect on the size windows of regions I, II and III because although ions would be emitted earlier during the deformation of the droplet, the exponential nature of the ion emission law fixes the maximum electric field at a precise fraction of $E^*$. The amount of charge emitted from a droplet of a given radius with initial charge $\eta Q_{Ray}$ will increase with the initial charge and should approximately be that of the droplet charged $4\,\%$ over the Rayleigh limit plus $(\eta -1.04)Q_{Ray}$. Furthermore, using a small charge excess over the Rayleigh limit is representative of practical conditions: for example, charged droplets in aerosols reach the Rayleigh limit by evaporation of the liquid phase at constant charge, i.e. $Q_o/Q_{Ray}$ increases towards 1; and the charged nanodroplets produced during the jet breakup of an electrospray are, at most, slightly charged above the Rayleigh limit, because the breakup process distributes the charge and mass of the jet among neighbouring droplets to minimize the excess charge (Misra & Gamero-Castaño 2022).

4. Conclusions

We have developed an electrohydrodynamic phase field model to study the coupling of ion emission and Coulomb explosion in charged liquid droplets. The model accounts for the finite electric conductivity of the liquid. We draw the following major conclusions from the current study.

1. (i) Highly viscous droplets charged slightly above the Rayleigh limit emit roughly $36\,\%$ of their charge in the form of smaller droplet progenies. This analysis does not take into account ion field emission and, therefore, is valid for relatively large droplets where the electric field on the surface is not sufficient to support ion emission. For the EMI-Im droplets and conditions considered in the current article, we expect droplets with diameter $D > 140$ nm to follow this trend.

2. (ii) Ion field emission suppresses Coulomb explosions in droplets that are charged above the Rayleigh limit and are sufficiently small ($D < 45$ nm in our simulations with EMI-Im droplets). Between $20\text { and }40\,\%$ of the initial charge of the droplet is shed by ion emission (we assume an initial charge $4\,\%$ above the Rayleigh limit).

3. (iii) There is an intermediate size range ($45 \lesssim D \lesssim 100$ nm in our simulations with EMI-Im droplets) in which droplets charged above the Rayleigh limit undergo Coulomb explosions, and also emit ions from the vertices and jets of the deforming droplet. We find that approximately $24\,\%$ of the initial charge of the droplet is shed by ion emission and estimate that roughly $15\unicode{x2013}16\,\%$ is lost in the form of droplet progenies. This prediction matches well the experiments and analysis of Perez-Lorenzo & Fernández De La Mora (2022).

4. (iv) Ion emission freezes the value of the electric field in areas where emission is significant. We quantify this value at $E_{n}^{g}/E^{*}\sim 0.76\unicode{x2013}0.77$.

The continuum approach and the leaky dielectric model have been successful in reproducing Coulomb explosions in the past (Giglio et al. 2008; Gawande et al. 2017, 2020; Giglio et al. 2020). However, these prior studies compare numerical results with experimental data for droplets 10–50 $\mathrm {\mu }{\rm m}$ in diameter (Duft et al. 2003; Achtzehn et al. 2005; Giglio et al. 2008). Ion emission is relevant only in nanometric droplets and, in these systems, the continuum approximation as well as key assumptions of the leaky dielectric model such as constant electrical conductivity may be less accurate. However, we note that molecular dynamic simulations of ion emission and Coulomb explosions in droplets with diameter less than $10$ nm (Luedtke et al. 2008) show good qualitative agreement with our calculations. Moreover, recent continuum studies of ion emission from Taylor cones provide confidence in our treatment of this problem (Higuera 2008; Coffman et al. 2016, 2019; Gallud & Lozano 2022). Regarding the validity of the leaky dielectric model, it has been shown that this formulation accurately reproduces the similar problem of cone-jets (Gamero-Castaño & Magnani 2019). However, we plan to improve the present model by using an electrokinetic formulation, thus relaxing the assumption of a constant electrical conductivity.