2.1. Computational domain

Figure 1 shows the computational domain. We use cylindrical coordinates $\{x,r\}$. The emitter is modelled as an equipotential cylinder with surfaces $\varSigma _{1e}$ and $\varSigma _{2e}$, supporting a Taylor cone with a free surface $\varSigma _{12}\cup \varSigma _{34}$. Here $R(x)$ is the radius of the cross-section of the cone. The domain is enclosed by an outer boundary formed by a plane $\varSigma _{1c}$ perpendicular to the emitter, and the surface $\varSigma _{1o}$ placed far enough from the meniscus so that the particular choice of position does not significantly affect the solution. The bulk of the Taylor cone is divided into regions 2 and 3: the ion emission regime is characterized by very small flow rates and the fluid motion and dissipation effects are expected to be significant only near the tip of the meniscus, region 3. Similarly, the space surrounding the Taylor cone is divided into regions 1 and 4 to take advantage of the sharp decrease of the space charge of the beam away from the emission area (the space charge density is only retained in region 4). The use of these regions reduces the computational time because regions 1 and 2 can be resolved with simplified equations and computational grids with lower resolution than in regions 3 and 4. The model equations are highly nonlinear and must be solved iteratively. The extent of regions 3 and 4 is determined in each iteration by monitoring the ion emission area and the volumetric charge density in the beam.

Figure 1. Sketch of the model domain.

2.2. Governing equations

The model is based on the leaky-dielectric formulation (Melcher & Taylor Reference Melcher and Taylor1969; Saville Reference Saville1997), extended to include ion-field emission, dissipation and the ion beam. The velocity $\boldsymbol {V}$ and pressure $P$ fields in the liquid fulfil the continuity and conservation of momentum equations,

(2.1)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{V}=0, \end{gather}

(2.2)\begin{gather} \rho\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{V}={-}\boldsymbol{\nabla} P+\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}_\mu, \end{gather}

where $\rho$ and $\mu$ are the density and viscosity of the liquid and $\boldsymbol {\tau }_\mu$ is the viscous stress tensor,

(2.3)\begin{equation} \boldsymbol{\tau}_\mu=\mu\left(\boldsymbol{\nabla}\boldsymbol{V}+{\boldsymbol{\nabla}\boldsymbol{V}}^{\rm T}\right). \end{equation}

Ohmic and viscous dissipation are expected to be important, and the equation of conservation of energy is used to compute the associated variation in temperature $T$,

(2.4)\begin{equation} \rho c\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla} T=k\nabla^2T+\boldsymbol{\tau}_\mu\boldsymbol{:}\boldsymbol{\nabla}\boldsymbol{V}+\boldsymbol{E}^i\boldsymbol{\cdot}\boldsymbol{J}^i, \end{equation}

where $c$ and $k$ are the heat capacity and thermal conductivity of the liquid, while $\boldsymbol {E}$ and $\boldsymbol {J}$ stand for the electric field and the current density. When needed, superscripts $i$ and $o$ denote inside and outside the Taylor cone, respectively. The electrical conductivity $K$ and the viscosity are strong functions of temperature, and these dependencies are modelled with exponential relations valid for ionic liquids (Okoturo & VanderNoot Reference Okoturo and VanderNoot2004; Leys et al. Reference Leys, Wübbenhorst, Preethy Menon, Rajesh, Thoen, Glorieux, Nockemann, Thijs, Binnemans and Longuemart2008),

(2.5)\begin{gather} K(T)=Y_K{\rm e}^{{B_K}/({T-T_K})}, \end{gather}

(2.6)\begin{gather} \mu(T)=Y_\mu {\rm e}^{{B_\mu}/({T-T_\mu})}, \end{gather}

where the constants $Y_K$, $B_K$, $T_K$, $Y_\mu$, $B_\mu$ and $T_\mu$ are liquid specific. All other liquid properties are regarded constant, independent of temperature. Following the leaky-dielectric formulation, the volumetric charge density inside the liquid is assumed to be negligible and Ohm's law is used to express the current density, $\boldsymbol {J}^i = K\boldsymbol {E}^i$. Conservation of charge and the irrotational nature of the electric field then provide an equation for the electric potential $\varPhi ^i$ inside the liquid

(2.7)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} (K\boldsymbol{E}^i) = 0\quad {\rm and} \quad \boldsymbol{E} ={-}\boldsymbol{\nabla} \varPhi \quad \rightarrow \quad K\nabla^2\varPhi^i+\boldsymbol{\nabla}\varPhi^i\boldsymbol{\nabla} K=0. \end{equation}

Due to the assumption of zero volumetric charge, the model cannot require the fulfilment of both conservation of charge and Gauss’ law for the electric field. To fulfil both principles, the model would need to consider a non-zero volumetric charge $\rho _e$ and a current density $\boldsymbol {J}^i = K\boldsymbol {E}^i + \rho _e \boldsymbol {V}$. In this case, the Navier–Stokes equation (2.2) would also include the volumetric force $\rho _e \boldsymbol {E}^i$. We would like to point out that, similarly to Gallud & Lozano (Reference Gallud and Lozano2022), we had also solved the problem including this volumetric force by using a fictitious volumetric charge density $\rho _e = ({\varepsilon \varepsilon _0 \boldsymbol {\nabla }\varPhi ^i\boldsymbol {\nabla } K})/{K}$. The resulting solution is nearly identical to that of the present model. However, we have decided not to use this fictitious volumetric charge because there is no physical or mathematical justification for it.

In the vacuum surrounding the Taylor cone the electric potential $\varPhi ^o$ fulfils Poisson equation,

(2.8)\begin{equation} \nabla^2\varPhi^o={-}\frac{\chi\omega}{\varepsilon_0}. \end{equation}

Here $\chi$ is the charge-to-mass ratio of the ions, $\omega$ is the ion mass density and $\varepsilon _0$ is the permittivity of vacuum. Following the leaky-dielectric formulation, net charge is only present on the surface of the Taylor cone, and it is modelled as a surface charge density $\sigma$ fulfilling the conservation equation

(2.9)\begin{equation} \frac{{\rm d}}{{\rm d} x}\left(r\sigma\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{t}\right)= R (1+{R^{'}}^2)^{1/2} (KE_n^i-J_\omega), \end{equation}

where $\boldsymbol {t}$ is the unit vector tangential to the surface and $J_\omega$ is the field-emitted ion current density. Here $E_n$ is the component of the electric field normal to the surface. Equation (2.9) balances the charge convected on the surface, the charge injected from the bulk, and the field-emitted charge. Additional equations that must be fulfilled at the surface $\varSigma _{12}\cup \varSigma _{34}$ include the jump of the normal components of the electric field, the balance of normal and tangential stresses, and the kinematic velocity constraint (the velocity component normal to the surface is coupled to the emitted current density),

(2.10)\begin{gather} E_n^o-\varepsilon E_n^i=\frac{\sigma}{\varepsilon_0}, \end{gather}

(2.11)\begin{gather} \boldsymbol{t}\boldsymbol{\cdot}\boldsymbol{\tau}_\mu\boldsymbol{n}=\boldsymbol{t}\boldsymbol{\cdot}\left(\boldsymbol{\tau}_M^o-\boldsymbol{\tau}_M^i\right)\boldsymbol{n}, \end{gather}

(2.12)\begin{gather} \gamma\boldsymbol{n}\boldsymbol{\cdot}\left(\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n}\boldsymbol{\mathsf{I}}\right)\boldsymbol{n}-P+\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\tau}_\mu\boldsymbol{n}=\boldsymbol{n}\boldsymbol{\cdot}\left(\boldsymbol{\tau}_M^o-\boldsymbol{\tau}_M^i\right)\boldsymbol{n}, \end{gather}

(2.13)\begin{gather} \boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{n}=\frac{J_\omega}{\chi\rho}. \end{gather}

Here $\varepsilon$ and $\boldsymbol {n}$ are the dielectric constant of the liquid and the unit vector normal to the surface, respectively. The Maxwell stress tensor $\boldsymbol {\tau }_M$ is given by

(2.14)\begin{equation} \boldsymbol{\tau}_M=\varepsilon\varepsilon_0\boldsymbol{E}\boldsymbol{E}-\varepsilon_0\frac{\boldsymbol{E}\boldsymbol{\cdot}\boldsymbol{E}}{2}\left[\varepsilon-\rho\left(\frac{\partial\varepsilon}{\partial\rho}\right)_T\right], \end{equation}

where the term $\rho (\partial \varepsilon /\partial \rho )_T$ is zero due to the incompressibility assumption. The emitted ions are accelerated away from the surface by the electric field. They are modelled as a continuum where the only significant force is that induced by the average electric field, i.e. we neglect ion–ion collisions. The continuity and conservation of momentum equations for the ion density and ion velocity $\boldsymbol {v}$ fields are

(2.15)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\omega\boldsymbol{v}\right)=0, \end{gather}

(2.16)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\omega\boldsymbol{v}\boldsymbol{v}\right)=\chi\omega\boldsymbol{E}^o. \end{gather}

Ion evaporation is modelled with the kinetic law proposed by Iribarne & Thomson (Reference Iribarne and Thomson1976). This equation relates the ion current density with the energy barrier that emitted ions must overcome

(2.17)\begin{equation} J_\omega=\frac{k_BT}{h}\sigma\exp\left(-\frac{{\Delta G}_S-G_E}{k_BT}\right), \end{equation}

where $k_B$ and $h$ are the Boltzmann and the Planck constants, ${\Delta G}_S$ is the ion solvation energy and $G_E$ is the electrostatic contribution to the energy barrier. A common approach for determining $G_E$ consists of computing the electric field acting on the emitted ion as the superposition of the fields induced by its image charge and the surface charge (Müller Reference Müller1941, Reference Müller1956; Iribarne & Thomson Reference Iribarne and Thomson1976). In this article we use an improved formulation that takes into account the dielectric properties and the local curvature of the emitting surface (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2022),

(2.18)\begin{equation} G_E=\left(1-\frac{R_c}{r^{*}}\right)q R_c E_n^o+\frac{q^2(\varepsilon-1)}{8{\rm \pi}\varepsilon_0 R_c}\sum_{k=0}^\infty{\frac{1}{\varepsilon+1+\dfrac{1}{k}}\left(\frac{R_c}{r^{*}}\right)^{2k+2}}, \end{equation}

where $R_c$ is the local average curvature of the surface, $r^*$ is the position of the energy barrier respect to the surface and $q$ is the charge of the ion, equal to the fundamental charge $e$ for singly charged ions. Here $r^*$ is obtained numerically (Magnani & Gamero-Castaño Reference Magnani and Gamero-Castaño2022).

2.3. Characteristic scales and dimensionless equations

We estimate the characteristic length $L_c$ of the ion-emitting region with a simplified model for ion emission from an ideal Taylor cone. The electric field on the surface of the ideal Taylor cone (Taylor Reference Taylor1964) is given by

(2.19)\begin{equation} E_T(r)=\sqrt{\frac{2\gamma}{\varepsilon_0}\frac{\cos\theta_T}{r}}, \end{equation}

where $r$ is the radial cylindrical coordinate and $\theta _T \cong 49.29^\circ$ is the cone semiangle (Taylor Reference Taylor1964). The electric field is singular at the vertex and will induce ion emission. Our simplified model uses (2.19) to approximate $E_n^o$, and combines (2.9), (2.10) and (2.17) to eliminate $E_n^i$ and $\sigma$ and obtain an explicit expression for the current density of emitted ions (we assume negligible surface velocity),

(2.20)\begin{equation} J_\omega(r)=\dfrac{E_T(r)}{\dfrac{\varepsilon}{K}+\dfrac{h}{\varepsilon_0 k_B T}\exp\left(\dfrac{{\Delta G}_S-G_E}{k_B T}\right)}. \end{equation}

This expression resembles Ohms law with two resistances in series and driven by the electric field: the resistance $\varepsilon /K$ is associated with the injection of charge on the surface from the bulk, and $({h}/{\varepsilon _0 k_B T})\exp (({{\Delta G}_S-G_E})/{k_B T})$ is associated with the evaporation of ions from the surface. Thus, the current density exhibits two limiting behaviours depending on the dominance of each mechanism: far from the vertex the electric field is small, the energy barrier is large and ion evaporation restricts the transport

(2.21)\begin{equation} J_\omega(r\rightarrow\infty) \cong \frac{\varepsilon_0 k_B T}{h}\exp\left(-\frac{{\Delta G}_S-G_E}{k_B T}\right) E_T, \end{equation}

while near the vertex the electric field is large, the energy barrier is negligible, and the emission is restricted by the availability of charge injected from the bulk

(2.22)\begin{equation} J_\omega(r\rightarrow 0) \cong \frac{K}{\varepsilon}E_T. \end{equation}

In the first limit the surface charge is in near equipotential balance with the outer electric field, $\sigma \cong \varepsilon _0E_T$, while in the second limit ion emission depletes the surface charge exponentially, $\sigma \cong [({k_B\varepsilon _0}/{hK})\varepsilon T\exp (-(({{\Delta G}_S-G_E})/{k_B T}))]^{-1} \varepsilon _0 E_T$. Note also that $r J_\omega (r)$ tends to zero towards the vertex, and therefore the total emitted current $({2{\rm \pi} }/{\sin \theta _T})\int _{0}^{\infty } r J_\omega (r) \,{\rm d}r$ remains finite. This result is non-trivial given the exponential dependence of the ion emission law on the electric field (2.17), and the singular behaviour of the field at the vertex. However, this analysis comes with the caveat that our approximation of the electric field, $E_T$, may be inaccurate near the vertex due to the depletion of the surface charge, and will need to be confirmed with the numerical solution of the first-principles model. Figure 2 plots (2.20) and its two limiting behaviours. The current density is largest near the axis and falls sharply soon after, clearly defining a region from which most of the emission takes place. We define the characteristic length of the ion-emitting region as the value of the radial coordinate where both limiting behaviours are equal, which can be written explicitly when using the energy barrier for a planar conductor, $G_E = \sqrt {q^3E_n^o / (4{\rm \pi} \varepsilon _0)}$ (Iribarne & Thomson Reference Iribarne and Thomson1976),

(2.23)\begin{equation} L_c=\frac{\cos\theta_T q^6 \gamma}{8{\rm \pi}^2{\varepsilon_0}^3{{\Delta G}_S}^4}\left[1-\frac{k_B T}{{\Delta G}_S} \ln\left(\frac{\varepsilon\varepsilon_0 k_B T}{h K}\right)\right]^{{-}4} = \alpha \frac{\cos\theta_T q^6 \gamma }{8{\rm \pi}^2{\varepsilon_0}^3 {{\Delta G}_S}^4}, \end{equation}

where $({k_B T}/{{\Delta G}_S})\ln ({\varepsilon \varepsilon _0 k_B T}/{h K})\ll 1$ for typical ion emission conditions, e.g. it has a value of 0.101 for ${\Delta G}_S=1.6$ eV, $T=25\,^\circ$C, $\varepsilon =10$ and $K=1\ {\rm S}\ {\rm m}^{-1}$. The size of the ion-emitting region is a strong function of the ion solvation energy, proportional to the surface tension, and a weak function of the dielectric constant, the electrical conductivity and temperature.

Figure 2. Estimation of the ion current density emitted from a Taylor cone, (2.20), computed with $T= 25\,^\circ {\rm C}$, $\varepsilon = 10$, ${\Delta G}_S = 1.6$ eV, $\gamma = 0.05\ {\rm N}\ {\rm m}^{-1}$ and $K= 1\ {\rm S}\ {\rm m}^{-1}$.

Table 1 shows the characteristic scales employed in the non-dimensionalization of the equations. The characteristic viscosity and conductivity are the values at the upstream temperature $T_0$. The characteristic electric field is equal to the electric field on the surface of an ideal Taylor cone at $r=L_c$. The characteristic fluid velocity is obtained using the relation between the total emitted current and the flow rate $Q$, which in the ion emission regime are proportional,

(2.24)\begin{equation} Q=\frac{2{\rm \pi}}{\chi\rho}\int_S{J_\omega}\,{\rm d}S. \end{equation}

From this equation the characteristic flow rate is defined as

(2.25)\begin{equation} Q_c=\frac{2{\rm \pi}\varPhi_c K_c L_c}{\chi\rho}, \end{equation}

yielding the characteristic fluid velocity $Q_c/{\rm \pi} {L_c}^2$. We list below the governing equations in dimensionless form, with dimensionless variables written without additional markings,

(2.26)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{V}=0, \end{gather}

(2.27)\begin{gather} \boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{V}={-}\frac{\varPi_2}{2\cos\theta_T}\boldsymbol{\nabla} P+\frac{1}{{{Re}}}\left[\mu\nabla^2\boldsymbol{V}+\left(\boldsymbol{\nabla}\boldsymbol{V}+{\boldsymbol{\nabla}\boldsymbol{V}}^{\rm T}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}\mu\right], \end{gather}

(2.28)\begin{gather} Pe\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{\nabla} T-\nabla^2T=\frac{2\varPi_3}{{{Re}}\varPi_2}\mu\left(\boldsymbol{\nabla}\boldsymbol{V}+{\boldsymbol{\nabla}\boldsymbol{V}}^{\rm T}\right)\boldsymbol{:}\boldsymbol{\nabla}\boldsymbol{V}+K\boldsymbol{\nabla}\varPhi^i\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi^i, \end{gather}

(2.29)\begin{gather} K\nabla^2\varPhi^i={-}\boldsymbol{\nabla}\varPhi^i\boldsymbol{\cdot}\boldsymbol{\nabla} K, \end{gather}

(2.30)\begin{gather} \nabla^2\varPhi^o={-}\omega, \end{gather}

(2.31)\begin{gather} E_n^o-\varepsilon E_n^i=\frac{\sigma}{2\varPi_3}, \end{gather}

(2.32)\begin{gather} \frac{{\rm d}}{{\rm d} x}\left(r\sigma\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{t}\right)=R (1+{R^{'}}^2)^{1/2} (KE_n^i-J_\omega), \end{gather}

(2.33)\begin{gather} J_\omega=\varPi_1 T\sigma\exp\left(-\frac{\varPi_G-G_E}{T}\right), \end{gather}

(2.34)\begin{gather} G_E=\varPi_4\left[\left(1-\frac{R}{r^{*}}\right)R E_n^o+\varPi_5\frac{\varepsilon-1}{R}\sum_{k=0}^\infty{\frac{1}{\varepsilon+1+\dfrac{1}{k}}\left(\frac{R}{r^{*}}\right)^{2k+2}}\right], \end{gather}

(2.35)\begin{gather} \boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{n}=\frac{J_\omega}{2}, \end{gather}

(2.36)\begin{gather} \omega\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v}+\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\omega=0, \end{gather}

(2.37)\begin{gather} \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{v}={-}\varPi_2\varPi_3\boldsymbol{\nabla}\varPhi^o, \end{gather}

(2.38)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{n}-P+\frac{4\cos\theta_T}{{{Re}}\varPi_2}\mu\boldsymbol{n}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{V}\boldsymbol{n}=\cos\theta_T\left[{E_n^o}^2-\varepsilon{E_n^i}^2+\left(\varepsilon-1\right){E_t}^2\right], \end{gather}

(2.39)\begin{gather} E_t\sigma=\frac{2\varPi_3}{{{Re}}\varPi_2}\mu\boldsymbol{t}\boldsymbol{\cdot}\left(\boldsymbol{\nabla}\boldsymbol{V}+{\boldsymbol{\nabla}\boldsymbol{V}}^{\rm T}\right)\boldsymbol{n}. \end{gather}

The equations include the 10 dimensionless numbers shown in table 2.

Table 1. Characteristic scales used to non-dimensionalize the equations.

Table 2. Dimensionless numbers parametrizing the solution.

3.1. Electrostatic solution

The electric field, the surface charge and the flux of field-evaporated ions are computed by solving (2.29)–(2.33) which, after algebraic manipulations and discretization, are written as a single linear system of algebraic equations. First, we note that the right-hand sides of (2.29) and (2.30) are significant only near the tip of the meniscus, i.e. in regions 3 and 4, respectively. In these regions we solve the equations using finite differences in orthogonal grids (see Appendix A) using a second-order central difference for the space derivatives. In regions 1 and 2 these equations are well approximated by the Laplace equation, and solved using the boundary element method (BEM) (Brebbia & Dominguez Reference Brebbia and Dominguez1994; Brebbia, Telles & Wrobel Reference Brebbia, Telles and Wrobel2012; Bakr Reference Bakr2013). This method discretizes and solves the Laplace equation only on the boundary of the domain by converting it into the linear system $\boldsymbol{\mathsf{H}}\boldsymbol {\varPhi }=\boldsymbol{\mathsf{G}}\boldsymbol {\varPhi }'$, where the potential and its normal derivative at the nodes of the discretized boundary are related by the $\boldsymbol{\mathsf{H}}$ and $\boldsymbol{\mathsf{G}}$ matrices computed with standard BEM techniques. Furthermore, on the surface of the meniscus we combine (2.31)–(2.33) to eliminate $J_\omega$ and $\sigma$ and obtain a single equation relating the normal components of the electric field on both sides of the surface, $E_n^i$ and $E_n^o$, as follows:

(3.1)\begin{align} E_n^i&=\dfrac{E_n^o}{\varepsilon+\dfrac{K}{2\varPi_3\left[\varPi_1 T\exp\left(\dfrac{G_E-\varPi_G}{T}\right)+\dfrac{t_r}{r}V_t+\dfrac{{\rm d}V_t}{{\rm d}t}\right]}}\nonumber\\ &\quad +\dfrac{\dfrac{{\rm d}\sigma}{{\rm d}t}V_t}{K+2\varepsilon\varPi_3\left[\varPi_1T\exp\left(\dfrac{G_E-\varPi_G}{T}\right)+\dfrac{t_r}{r}V_t+\dfrac{{\rm d}V_t}{{\rm d}t}\right]}. \end{align}

Here $V_t$ is the surface tangential velocity and $t_r$ is the radial component of the tangential vector. This equation is written for simplicity as

(3.2)\begin{equation} E_n^i=B_E E_n^o+B_\sigma, \end{equation}

where the coefficients $B_E$ and $B_{\sigma }$ are evaluated with the solution from the previous iteration.

The BEM applied to region 1 is written in matrix form as

(3.3)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{H}}_{12} & \boldsymbol{\mathsf{H}}_{1o} & \boldsymbol{\mathsf{H}}_{14} & \boldsymbol{\mathsf{G}}_{12} & -\boldsymbol{\mathsf{G}}_{1c} & \boldsymbol{\mathsf{G}}_{14}\end{bmatrix}\begin{bmatrix} \boldsymbol{\varPhi}_{12}\\ \boldsymbol{\varPhi}_{1o}\\ \boldsymbol{\varPhi}_{14}\\ \boldsymbol{E}_n^o\\ \boldsymbol{E}_{1e}\\ \boldsymbol{E}_{1c}\\ \boldsymbol{E}_{14}\end{bmatrix}={-}\begin{bmatrix} \boldsymbol{\mathsf{H}}_{1e} & \boldsymbol{\mathsf{H}}_{1c} & \boldsymbol{\mathsf{G}}_{1o}\end{bmatrix}\begin{bmatrix} \boldsymbol{\varPhi}_{1e}\\ \boldsymbol{\varPhi}_{1c}\\ \boldsymbol{E}_{1o}\end{bmatrix}, \end{equation}

where the subindexes in all elements refer to the labelling of the boundaries in figure 1. The $\boldsymbol {\varPhi }_{1e}$, $\boldsymbol {\varPhi }_{1c}$ and $\boldsymbol {E}_{1o}$ vectors are boundary conditions,

(3.4a–c)\begin{equation} \boldsymbol{\varPhi}_{1e}=\varPi_\varPhi, \quad \boldsymbol{\varPhi}_{1c}=0, \quad \boldsymbol{E}_{1o} = 0. \end{equation}

The BEM applied to region 2 is written in matrix form as

(3.5)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{H}}_{21} & \boldsymbol{\mathsf{H}}_{23} & \boldsymbol{\mathsf{G}}_{21} & -\boldsymbol{\mathsf{G}}_{2e} & -\boldsymbol{\mathsf{G}}_{23}\end{bmatrix}\begin{bmatrix} \boldsymbol{\varPhi}_{21}\\ \boldsymbol{\varPhi}_{23}\\ \boldsymbol{E}_n^i\\ \boldsymbol{E}_{2e}\\ \boldsymbol{E}_{23}\end{bmatrix}={-}\left[\boldsymbol{\mathsf{H}}_{2e}\right]\left[\boldsymbol{\varPhi}_{2e}\right], \end{equation}

with boundary condition

(3.6)\begin{equation} \boldsymbol{\varPhi}_{2e}=\varPi_\varPhi. \end{equation}

In region 3 equation (2.29) is written in the $\{\xi _3,\eta _3\}$ orthogonal coordinate system using the relations defined in Appendix A:

(3.7)\begin{align} &\left(\dfrac{\partial\xi_3}{\partial x}^2+\dfrac{\partial\xi_3}{\partial r}^2\right)\dfrac{\partial^2\varPhi^i}{\partial{\xi_3}^2}+ \left(\dfrac{\partial\eta_3}{\partial x}^2+\dfrac{\partial\eta_3}{\partial r}^2\right)\dfrac{\partial^2\varPhi^i}{\partial{\eta_3}^2}+ \left(\dfrac{\partial^2\xi_3}{\partial x^2}+\dfrac{\partial^2\xi_3}{\partial r^2}+\dfrac{1}{r}\dfrac{\partial\xi_3}{\partial r}\right)\dfrac{\partial\varPhi^i}{\partial\xi_3}\nonumber\\ &\qquad +\left(\dfrac{\partial^2\eta_3}{\partial x^2}+\dfrac{\partial^2\eta_3}{\partial r^2}+\dfrac{1}{r}\dfrac{\partial\eta_3}{\partial r}\right)\dfrac{\partial\varPhi^i}{\partial\eta_3}={-}\dfrac{1}{K}\left[ \dfrac{\partial K}{\partial\xi_3}\left(\dfrac{\partial\xi_3}{\partial x}^2+\dfrac{\partial\xi_3}{\partial r}^2\right)\dfrac{\partial\varPhi^i}{\partial\xi_3}\right.\nonumber\\ &\qquad + \left. \dfrac{\partial K}{\partial\eta_3}\left(\dfrac{\partial\eta_3}{\partial x}^2+\dfrac{\partial\eta_3}{\partial r}^2\right)\dfrac{\partial\varPhi^i}{\partial\eta_3} \right]. \end{align}

This equation is discretized using finite differences and written in matrix form

(3.8)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{M}}_3 & \boldsymbol{\mathsf{B}}_{34} & \boldsymbol{\mathsf{B}}_{23}\end{bmatrix}\begin{bmatrix} \boldsymbol{\varPhi}_{3}\\ \boldsymbol{\varPhi}_{34}\\ \boldsymbol{\varPhi}_{23}\end{bmatrix}=\left[\boldsymbol{F}_3\right], \end{equation}

where $\boldsymbol{\mathsf{M}}_3$ and $\boldsymbol {F}_3$ are the matrix and the forcing term resulting from the discretization, while $\boldsymbol{\mathsf{B}}_{23}$ and $\boldsymbol{\mathsf{B}}_{34}$ are matrices associated with the matching of boundary conditions in the interfaces with regions 2 and 4, i.e. the electric potential and the normal component of the electric field are continuous across these interfaces. Matrix $\boldsymbol{\mathsf{M}}_3$ also contains the symmetry condition at the axis,

(3.9)\begin{equation} \frac{\partial\varPhi^3}{\partial r}(\xi_3,0)=0. \end{equation}

Similarly, (2.30) is solved in region 4 in an orthogonal grid and, after discretization, written as

(3.10)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{M}}_4 & \boldsymbol{\mathsf{B}}_{43} & \boldsymbol{\mathsf{B}}_{41}\end{bmatrix}\begin{bmatrix} \boldsymbol{\varPhi}_{4}\\ \boldsymbol{\varPhi}_{34}\\ \boldsymbol{\varPhi}_{14}\end{bmatrix}={-}\left[\boldsymbol{\omega}\right], \end{equation}

where the matrix $\boldsymbol{\mathsf{M}}_4$ contains the discretization of the equation and the boundary condition at the axis

(3.11)\begin{equation} \frac{\partial\varPhi^4}{\partial r}(\xi_4,0)=0 \end{equation}

while matrices $\boldsymbol{\mathsf{B}}_{43}$ and $\boldsymbol{\mathsf{B}}_{41}$ contain the matching of boundary conditions between regions.

Equations (3.2), (3.3), (3.5), (3.8) and (3.10) are merged into a single linear system of algebraic equations to solve for the electric potential in regions 3 and 4, $\boldsymbol {\varPhi }^3$ and $\boldsymbol {\varPhi }^4$, the electric potential on the surface of the meniscus, $\boldsymbol {\varPhi }_{12}$ and $\boldsymbol {\varPhi }_{34}$, and the normal components of the electric field on either side of the surface, $E_n^i$ and $E_n^o$. With this solution the tangential component of the electric field is computed by differentiating the electric potential along the surface, while the surface charge density is obtained from (2.31).

3.2. Fluid dynamic solution

The continuity and momentum equations, (2.26) and (2.27), are solved in the stream function-vorticity formulation, which in axisymmetric problems reduces the system of three coupled equations for the two velocity components and the pressure to a system of two coupled equations for the stream function $\varPsi$ and the vorticity $\varOmega$,

(3.12)\begin{equation} \frac{\partial^2\varPsi}{\partial r^2}+ \frac{\partial^2\varPsi}{\partial x^2}-\frac{1}{r} \frac{\partial\varPsi}{\partial r}={-}r\varOmega, \end{equation}

(3.13)\begin{align} & \frac{\partial^2\varOmega}{\partial r^2}+\frac{\partial^2\varOmega}{\partial x^2}+\frac{1}{r}\frac{\partial\varOmega}{\partial r}-\frac{\varOmega}{r^2}=\frac{{{Re}}}{\mu r}\left(\frac{\partial\varPsi}{\partial r}\frac{\partial\varOmega}{\partial x}-\frac{\partial\varPsi}{\partial x}\frac{\partial\varOmega}{\partial r}+\frac{1}{r}\frac{\partial\varPsi}{\partial x}\varOmega\right)\nonumber\\ &\qquad -\frac{1}{\mu}\left[\frac{1}{r}\left(\frac{\partial^2\varPsi}{\partial x^2}-\frac{\partial^2\varPsi}{\partial r^2}+\frac{1}{r}\frac{\partial\varPsi}{\partial r}\right)\left(\frac{\partial^2\mu}{\partial r^2}-\frac{\partial^2\mu}{\partial x^2}\right)\right.\nonumber\\ &\qquad + \left.\frac{2}{r}\left(\frac{1}{r}\frac{\partial\varPsi}{\partial x}-2\frac{\partial^2\varPsi}{\partial x\partial r}\right)\frac{\partial^2\mu}{\partial x\partial r}+\frac{\partial\varOmega}{\partial x}\frac{\partial\mu}{\partial x}+\frac{\partial\varOmega}{\partial r}\frac{\partial\mu}{\partial r}\right], \end{align}

with

(3.14a,b)\begin{gather} U=\frac{1}{r}\frac{\partial\varPsi}{\partial r}, \quad V={-}\frac{1}{r}\frac{\partial\varPsi}{\partial x}, \end{gather}

(3.15)\begin{gather} \boldsymbol{\varOmega}=\boldsymbol{\nabla}\times\boldsymbol{V}, \end{gather}

where $U$ and $V$ are the axial and radial components of the velocity. Equations (3.12) and (3.15) are solved only in region 3, where the liquid velocity is expected to be significant, using the $\{\xi _3,\eta _3\}$ orthogonal coordinates defined in Appendix A,

(3.16)$$\begin{gather} \left(\frac{\partial\xi_3}{\partial x}^2+\frac{\partial\xi_3}{\partial r}^2\right)\frac{\partial^2\varPsi}{\partial{\xi_3}^2}+ \left(\frac{\partial\eta_3}{\partial x}^2+\frac{\partial\eta_3}{\partial r}^2\right)\frac{\partial^2\varPsi}{\partial{\eta_3}^2}+ \left(\frac{\partial^2\xi_3}{\partial x^2}+\frac{\partial^2\xi_3}{\partial r^2}-\frac{1}{r}\frac{\partial\xi_3}{\partial r}\right)\frac{\partial\varPsi}{\partial\xi_3}\nonumber\\ +\left(\frac{\partial^2\eta_3}{\partial x^2}+\frac{\partial^2\eta_3}{\partial r^2}-\frac{1}{r}\frac{\partial\eta_3}{\partial r}\right)\frac{\partial\varPsi}{\partial\eta_3}={-}r\varOmega, \end{gather}$$

(3.17) \begin{align} &\left(\frac{\partial\xi_3}{\partial x}^2+\frac{\partial\xi_3}{\partial r}^2\right)\frac{\partial^2\varOmega}{\partial{\xi_3}^2}+ \left(\frac{\partial\eta_3}{\partial x}^2+\frac{\partial\eta_3}{\partial r}^2\right)\frac{\partial^2\varOmega}{\partial{\eta_3}^2}+ \left(\frac{\partial^2\xi_3}{\partial x^2}+\frac{\partial^2\xi_3}{\partial r^2}+\frac{1}{r}\frac{\partial\xi_3}{\partial r}\right)\frac{\partial\varOmega}{\partial\xi_3}\nonumber\\ &\quad +\left(\frac{\partial^2\eta_3}{\partial x^2}+\frac{\partial^2\eta_3}{\partial r^2}+\frac{1}{r}\frac{\partial\eta_3}{\partial r}\right)\frac{\partial\varOmega}{\partial\eta_3}=\frac{{{Re}}}{r\mu}\left[\left(\frac{\partial\eta_3}{\partial x}\frac{\partial\xi_3}{\partial r}-\frac{\partial\xi_3}{\partial x}\frac{\partial\eta_3}{\partial r}\right)\left(\frac{\partial\varPsi}{\partial\xi_3}\frac{\partial\varOmega}{\partial\eta_3}-\frac{\partial\varPsi}{\partial\eta_3}\frac{\partial\varOmega}{\partial\xi_3}\right)\right.\nonumber\\ &\quad +\left.\left(\frac{\partial\xi_3}{\partial x}\frac{\partial\varPsi}{\partial\xi_3}-\frac{\partial\eta_3}{\partial x}\frac{\partial\varPsi}{\partial\eta}\right)\frac{\varOmega}{r}\right] +\frac{1}{\mu}\left[\frac{1}{r}\left(\frac{\partial^2\varPsi}{\partial x^2}-\frac{\partial^2\varPsi}{\partial r^2}+\frac{1}{r}\frac{\partial\varPsi}{\partial r}\right)\left(\frac{\partial^2\mu}{\partial x^2}-\frac{\partial^2\mu}{\partial r^2}\right)\right.\nonumber\\ &\quad + \left.\frac{2}{r}\left(\frac{2\partial^2\varPsi}{\partial x\partial r}-\frac{1}{r}\frac{\partial\varPsi}{\partial x}\right)\frac{\partial^2\mu}{\partial x\partial r}- \left(\frac{\partial\xi_3}{\partial x}^2+\frac{\partial\xi_3}{\partial r}^2\right)\frac{\partial\varOmega}{\partial\xi_3}\frac{\partial\mu}{\partial\xi_3}+ \left(\frac{\partial\eta_3}{\partial x}^2+\frac{\partial\eta_3}{\partial r}^2\right)\frac{\partial\varOmega}{\partial\eta_3}\frac{\partial\mu}{\partial\eta_3}\right], \end{align}

where the last terms are given in $\{x,r\}$ coordinates for brevity. The equations, discretized using finite differences, fulfil the boundary conditions

(3.18)\begin{gather} \left. \begin{gathered} \varPsi(\xi_3,0)=0 \quad \frac{\partial\varPsi}{\partial\xi_3}(\xi_3,1)=\dfrac{-rJ_\omega}{2\sqrt{\dfrac{\partial\xi_3}{\partial x}^2+\dfrac{\partial\xi_3}{\partial r}^2}}\\ \dfrac{\partial\varPsi}{\partial\xi_3}(0,\eta_3)=0 \quad \varPsi(1,\eta_3)=0 \end{gathered} \right\}, \end{gather}

(3.19)\begin{gather} \left. \begin{gathered} \varOmega(\xi_3,0)=0 \quad \varOmega(\xi_3,1)=\dfrac{{{Re}}\varPi_2}{2\varPi_3}\dfrac{E_t\sigma}{\mu}-\dfrac{2}{r}\left[\dfrac{1}{t_x}\dfrac{\partial t_r}{\partial t}\dfrac{\partial\varPsi}{\partial n}-\dfrac{t_r}{r}\dfrac{\partial\varPsi}{\partial t}+\dfrac{\partial^2\varPsi}{\partial t^2}\right]\\ \frac{\partial\varOmega}{\partial\xi_3}(0,\eta_3)=0 \quad \varOmega(1,\eta_3)=0 \end{gathered} \right\}. \end{gather}

The linear systems of algebraic equations for the stream function and the vorticity resulting from the discretization and application of boundary conditions are solved separately to improve numerical stability. The boundary condition for the vorticity on the surface of the meniscus is derived from (2.39), and links the fluid dynamic and the electrostatic problems. The balance of tangential stresses on the surface is the main driver of the liquid flow, while the outward velocity due to ion emission has a smaller contribution (the flow rates of ion-emitting Taylor cones are low).

3.3. Fluid temperature and the variation of physical properties

The temperature can be written as $T(x, r)=T_0 + T_1(x, r)$, where $T_0$ is the upstream temperature. We solve (2.28) for $T_1(x,r)$ in region 3, using the $\{\xi _3,\eta _3\}$ orthogonal coordinates

(3.20)\begin{align} &\left(\frac{\partial\xi_3}{\partial x}^2+\frac{\partial\xi_3}{\partial r}^2\right)\frac{\partial^2T_1}{\partial{\xi_3}^2}+ \left(\frac{\partial\eta_3}{\partial x}^2+\frac{\partial\eta_3}{\partial r}^2\right)\frac{\partial^2T_1}{\partial{\eta_3}^2}+ \left(\frac{\partial^2\xi_3}{\partial x^2}+\frac{\partial^2\xi_3}{\partial r^2}+\frac{1}{r}\frac{\partial\xi_3}{\partial r}\right)\frac{\partial T_1}{\partial\xi_3}\nonumber\\ &\quad +\left(\frac{\partial^2\eta_3}{\partial x^2}+\frac{\partial^2\eta_3}{\partial r^2}+\frac{1}{r}\frac{\partial\eta_3}{\partial r}\right)\frac{\partial T_1}{\partial\eta_3}={-}\frac{Pe}{r}\left(\frac{\partial\xi_3}{\partial x}\frac{\partial\eta_3}{\partial r}+\frac{\partial\eta_3}{\partial x}\frac{\partial\xi_3}{\partial r}\right)\left(\frac{\partial\varPsi}{\partial\xi}\frac{\partial T_1}{\partial\eta}-\frac{\partial\varPsi}{\partial\eta}\frac{\partial T_1}{\partial\xi}\right)\nonumber\\ &\quad -\frac{2\varPi_3}{{{Re}}\varPi_2}\frac{\mu}{r^2}\left[2\frac{\partial^2\varPsi}{\partial x\partial r}^2+\left(\frac{\partial^2\varPsi}{\partial r^2}-\frac{\partial^2\varPsi}{\partial x^2}-\frac{1}{r}\frac{\partial\varPsi}{\partial r}\right)^2+2\left(\frac{\partial^2\varPsi}{\partial x\partial r}-\frac{1}{r}\frac{\partial\varPsi}{\partial x}\right)^2\right]\nonumber\\ &\quad -K\left[\left(\frac{\partial\xi_3}{\partial x}^2+\frac{\partial\xi_3}{\partial r}^2\right)\frac{\partial\varPhi^i}{\partial\xi_3}^2+\left(\frac{\partial\eta_3}{\partial x}^2+\frac{\partial\eta_3}{\partial r}^2\right)\frac{\partial\varPhi^i}{\partial\eta_3}^2\right], \end{align}

where several stream function terms in the right-hand side are written in $\{x,r\}$ coordinates for brevity. The boundary conditions for the temperature field are

(3.21)\begin{equation} \left. \begin{gathered} \frac{\partial T_1}{\partial\eta_3}(\xi_3,0)=0 \quad \frac{\partial T_1}{\partial\eta_3}(\xi_3,1)=0\\ T_1(0,\eta_3)=0 \quad \frac{\partial T_1}{\partial\xi_3}(1,\eta_3)=0 \end{gathered} \right\}. \end{equation}

The temperature equation (3.20) is discretized using finite differences, resulting in a linear system of algebraic equations that incorporates the boundary conditions. The solution for the temperature field is then used to compute the viscosity and conductivity near the tip of the meniscus, using (2.5) and (2.6).

3.4. Ion beam

We compute the density and velocity fields of the ion beam by integrating the equations of conservation of mass (2.36) and momentum (2.37),

(3.22)\begin{gather} u\frac{\partial\omega}{\partial x}+v\frac{\partial\omega}{\partial r}+\omega\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial r}+\frac{v}{r}\right)=0, \end{gather}

(3.23)\begin{gather} \left. \begin{aligned} u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial r} & ={-}\varPi_2\varPi_3\frac{\partial\varPhi^o}{\partial x},\\ u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial r} & ={-}\varPi_2\varPi_3\frac{\partial\varPhi^o}{\partial r}, \end{aligned} \right\} \end{gather}

using the method of characteristics. The distance along the trajectory of an ion is defined as the new coordinate of integration,

(3.24)\begin{equation} \frac{{\rm d}}{{\rm d}s}=\frac{\partial x}{\partial s}\frac{\partial}{\partial x}+\frac{\partial r}{\partial s}\frac{\partial}{\partial r}, \end{equation}

and, using this new coordinate, the equations are rewritten as

(3.25)\begin{equation} \left. \begin{gathered} \frac{{\rm d}u}{{\rm d}s}={-}\varPi_2\varPi_3\frac{\partial\varPhi^o}{\partial x} \quad \frac{{\rm d}v}{{\rm d}s}={-}\varPi_2\varPi_3\frac{\partial\varPhi^o}{\partial r},\\ \frac{{\rm d}\omega}{{\rm d}s}=\omega\left[\frac{\varPi_2\varPi_3}{u^2+v^2}\left(u\frac{\partial\varPhi^o}{\partial x}+v\frac{\partial\varPhi^o}{\partial r}\right)-\frac{1}{u^2+v^2}\left(u\frac{\partial v}{\partial p}-v\frac{\partial u}{\partial p}\right)-\frac{v}{r}\right], \end{gathered} \right\} \end{equation}

where $p$ is the coordinate normal to $s$. Characteristic lines are obtained by integrating

(3.26a,b)\begin{equation} \frac{\partial x}{\partial s}=u \quad {\rm and} \quad \frac{\partial r}{\partial s}=v. \end{equation}

The initial positions of the characteristic lines coincide with the surface of the meniscus, while the initial value of the ion density is obtained from the emitted current density and the initial ion velocity,

(3.27)\begin{equation} \omega_0=\frac{J_\omega}{\chi\left\|\boldsymbol{v}_0\right\|}. \end{equation}

We estimate the initial ion velocity $\boldsymbol {v}_0$ as the average velocity of the ions that can escape the energy barrier,

(3.28)\begin{equation} v_{avg}=\frac{\displaystyle\int_{\Delta G}^\infty{v f_M(v)}\,{\rm d}v}{\displaystyle\int_{\Delta G}^\infty{f_M(v)}\,{\rm d}v}, \end{equation}

where $\Delta G={\Delta G}_S-G_E$ is the local value of the energy barrier. The ion population in the meniscus is in thermal equilibrium and therefore characterized by a Maxwellian distribution $f_M(v)$. A good estimate for the initial velocity of the ions is that associated with (3.28), reduced by the climb of the energy barrier. In dimensionless quantities, the resulting initial ion density and velocity are

(3.29)\begin{equation} \left. \begin{gathered} \omega_0=\dfrac{J_\omega}{2\chi\left\|\boldsymbol{v}_0\right\|},\\ \left\|\boldsymbol{v}_0\right\|=\sqrt{\dfrac{T}{m}\left[\dfrac{2\left(\dfrac{\Delta G}{T}+1\right)}{2\sqrt{\dfrac{\Delta G}{T}}+\sqrt{\rm \pi}\exp\left(\dfrac{\Delta G}{T}\right)\operatorname{erfc}\sqrt{\dfrac{\Delta G}{T}}}\right]^2-\dfrac{2\Delta G}{m}}. \end{gathered} \right\} \end{equation}

3.5. Pressure in the meniscus

The pressure field in the meniscus is needed to enforce the equilibrium of normal stresses, (2.38). The pressure along the surface is computed by integrating the momentum equation using the latest solution for the velocity field. Equation (2.27) is first multiplied by the vector $\boldsymbol {t}$ tangential to the surface

(3.30)\begin{align} \frac{\partial}{\partial t}\left(\frac{\cos\theta_T}{\varPi_2}\boldsymbol{V}\boldsymbol{\cdot}\boldsymbol{V}+P\right)=\frac{2\cos\theta_T}{\varPi_2}\left[\boldsymbol{V}\times\boldsymbol{\varOmega}+\frac{\mu\nabla^2\boldsymbol{V}+\left(\boldsymbol{\nabla}\boldsymbol{V}+{\boldsymbol{\nabla}\boldsymbol{V}}^{\rm T}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}\mu}{{{Re}}}\right]\boldsymbol{\cdot}\boldsymbol{t}. \end{align}

This equation is then converted into the surface reference frame $\{t,n\}$ defined in Appendix A, and integrated. All terms containing ${\partial \mu }/{\partial n}$ are zero due to the boundary condition for the temperature field on the surface. The final expression for the pressure along the surface is

(3.31)\begin{align} P&=P_{tip}+\frac{2\cos\theta_T}{\varPi_2}\left[\int_0^t{\frac{\partial\varPsi}{\partial t}\frac{\varOmega}{r}}\,{\rm d}t -\frac{1}{2r^2}\left(\frac{\partial\varPsi}{\partial t}^2+\frac{\partial\varPsi}{\partial n}^2\right)\right]\nonumber\\ &\quad +\frac{2\cos\theta_T}{{{Re}}\varPi_2}\int_0^t{\mu\left(\frac{\partial\varOmega}{\partial n}-\frac{\varOmega t_x}{r}\right)+\frac{2}{r}\frac{\partial\mu}{\partial t}\left(\frac{1}{t_x}\frac{\partial t_r}{\partial t}\frac{\partial\varPsi}{\partial t}+\frac{t_r}{r}\frac{\partial\varPsi}{\partial n}-\frac{\partial^2\varPsi}{\partial t\partial n}\right)}\,{\rm d}t, \end{align}

where $P_{tip}$ is the pressure at the tip of the meniscus. The surface coordinate $t$ starts at the tip of the meniscus and progresses along the surface up to the emitter. Once the integration is performed, $P_{tip}$ is adjusted so that the pressure at the emitter matches the upstream condition $P=P_0$ (we set $P_0 = 0$).

3.6. Surface optimization

The position of the surface of the meniscus is computed in each iteration by minimizing the error in the balance of normal stresses. To do so (2.38) is written as

(3.32)$$\begin{gather} \frac{1}{\sqrt{\varrho^2+{\varrho'}^2}}\left(\frac{\varrho\varrho''-{\varrho'}^2}{\varrho^2 +{\varrho'}^2}+\frac{\varrho'}{\varrho\tan\theta}-2\right)={-}P-\cos\theta_T\left[{E_n^o}^2-\varepsilon{E_n^i}^2+\left(\varepsilon-1\right){E_t}^2\right]\nonumber\\ +\, \frac{4\cos\theta_T}{{{Re}}\varPi_2}\frac{\mu}{r}\left(\frac{1}{n_r}\frac{\partial n_r}{\partial n}\frac{\partial\varPsi}{\partial n}+\frac{t_x}{r}\frac{\partial\varPsi}{\partial t}+\frac{\partial^2\varPsi}{\partial t\partial n}\right), \end{gather}$$

where the left-hand side is the surface tension stress written in terms of the position of the surface in polar coordinates $\varrho =\varrho (\theta )$. The right-hand side of (3.32) is evaluated with the latest solution. This formulation in polar coordinates (originating at the tip of the emitter, see Appendix B) avoids numerical instabilities both at the vertex of the meniscus ($\theta =0$) and at the anchoring point, $\theta ={\rm \pi} /2$. The use of Cartesian coordinates would result in very large or infinite derivatives at these points. The apparent singularity at $\theta =0$ in the $\varrho '/\varrho \tan \theta$ term is resolved with a Taylor expansion of $\varrho '$ in the vicinity of the vertex, which shows the term to be finite for a tip with a finite radius of curvature,

(3.33)\begin{equation} \lim_{\theta\rightarrow 0}{\frac{\varrho'}{\varrho\tan\theta}}=\lim_{\theta\rightarrow 0}{\frac{\varrho'(0)+\varrho''(0)\theta}{\varrho\theta}}=\frac{\varrho''(0)}{\varrho(0)}. \end{equation}

Equation (3.32) is used to optimize the surface by applying the integral least-square technique. The residual along the surface is defined as

(3.34)\begin{gather} \mathcal{E}=\frac{1}{2}\int_0^{{\rm \pi}/2}f^2(\varrho(\theta),y(\theta),\theta)+\left(\varrho'-y\right)^2 \,{\rm d}\theta, \end{gather}

(3.35)$$\begin{gather} f(\varrho,y,\theta)=\frac{\varrho y'-y^2}{\varrho^2+y^2}+\frac{y}{\varrho\tan\theta}-2+\sqrt{\varrho^2+y^2}\left\{\cos\theta_T\left[{E_n^o}^2-\varepsilon{E_n^i}^2+\left(\varepsilon-1\right){E_t}^2\right]\right.\nonumber\\ +\left. P-\frac{4\cos\theta_T}{{{Re}}\varPi_2}\frac{\mu}{r}\left(\frac{1}{n_r}\frac{\partial n_r}{\partial n}\frac{\partial\varPsi}{\partial n}+\frac{t_x}{r}\frac{\partial\varPsi}{\partial t}+\frac{\partial^2\varPsi}{\partial t\partial n}\right)\right\}, \end{gather}$$

where the first derivative $\varrho '$ has been converted into a second optimization variable $y$ to improve the stability of the algorithm. To minimize the residual $\mathcal {E}$, the surface is discretized with the same set of nodes used for the electrostatic problem, and the resulting integration provides two equations at each point $j$ for the surface position $\varrho _j$ and its first derivative $y_j$,

(3.36)\begin{gather} \int_0^{{\rm \pi}/2}f\frac{\partial f}{\partial\varrho_j}+\left(\varrho'-y\right)\frac{\partial\varrho'}{\partial\varrho_j}\,{\rm d}\theta=0, \end{gather}

(3.37)\begin{gather} \int_0^{{\rm \pi}/2}f\frac{\partial f}{\partial y_i}+\left(\varrho'-y\right)\frac{\partial\varrho'}{\partial y_j}\,{\rm d}\theta=0. \end{gather}

This system of algebraic equations is solved using the Newton's method. Regions 3 and 4 are redefined in each iteration with the updated surface. The transition between regions 2 and 3 is set at the point in the surface where the ion current density drops to $10^{-12}$ of its maximum which, due to the exponential nature of the ion emission law, is always near the vertex of the meniscus; and the transition between regions 1 and 4 is set at the axial point where the ion density drops to $10^{-8}$ of its maximum.

3.7. Breakdown of the continuum hypothesis near the vertex of the meniscus

Ion emission is localized in a small region surrounding the tip of the meniscus. The continuum hypothesis breaks down within a distance from the vertex of a few molecular radii (for the ionic liquids of interest in this study, the molecular radius if approximately one half of a nanometre), and we do not expect nor require the model equations to be accurate in this molecular-size region. Instead, we define a Knudsen number in terms of the distance $\delta$ from the vertex of the meniscus

(3.38)\begin{equation} K_n=\frac{\lambda_m}{2\delta}=\frac{1}{\delta}\left(\frac{3 m}{4{\rm \pi}\rho}\right)^{1/3}, \end{equation}

where the mean free path $\lambda _m$ has been expressed in terms of the molecular mass $m$ and the liquid density, use a critical value of the Knudsen number to define the extent of the region where the continuum hypothesis breaks down, and extrapolate the continuum solution into this region by using a patching function that preserves axial symmetry

(3.39)\begin{equation} \mathcal{F}(\delta<\delta_K)=\mathcal{F}(\delta_K)+\frac{\mathcal{F}'(\delta_K)}{2}\left(\frac{\delta^2}{\delta_K}-\delta_K\right). \end{equation}

Here $\delta _K$ is the position on the surface where the Knudsen number reaches the critical value. We typically adopt a critical value of 0.1 for the Knudsen number, or equivalently $\delta _K=3.93$ nm, and use (3.39) to evaluate the stream function, the vorticity, the pressure and the viscous stress near the vertex.

3.8. Solving scheme

The input parameters in the model are the radius of the emitter, the distance and the potential difference between the emitter and the extractor electrode, the physical properties of the liquid and the upstream temperature. The pressure jump $P_0$ across the surface is set to zero at the anchor point on the emitter. The solution of the model provides the shape of the meniscus; the velocity, pressure and temperature fields in the liquid phase; the electric field and surface charge; and the density and velocity fields of the ion beam. Figure 3 shows a flow chart of the algorithm used to solve the system of equations. First, we make a guess for the shape of the meniscus, assume a zero velocity field and a constant temperature throughout the fluid. With these values we start the iterative procedure: the sizes of regions 3 and 4 are determined (in the first iteration we use $10\times L_c$ as initial guess); the boundaries of the four regions are discretized and the orthogonal grids of regions 3 and 4 computed. The grid discretizations of regions 3 and 4 have $100\times 50$ and $50\times 100$ points, respectively, with higher point density near the tip of the meniscus and close to the surface of the liquid (see Appendix A). The boundary of region 1 is discretized using 1400 points, while region 2 uses 1100 points. One thousand of these points are shared between the two regions, and are used to discretize the portion of the surface not included in the orthogonal grids of regions 3 and 4. The boundary elements of this portion of the surface have variable size for a smooth transition between the spacing of the orthogonal grids and the rest of the domain discretization. The electric problem, the fluid dynamic problem and the temperature equation are solved in this order; the viscosity and electrical conductivity are updated with the computed temperature field; and the ion beam is computed. These steps are iterated until convergence of the total ion current. At this point the residual of the surface normal stresses is evaluated and if it is below a specified threshold, typically $10^{-5}$, the solution is accepted as final, otherwise the surface is modified using (3.36) and (3.37) and, with the updated shape of the meniscus, a new iteration is started.

Figure 3. Flow chart of the algorithm for solving the system of equations.