# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

Eric Giglio, J. Rangama, S. Guillous, T. Le Cornu

► To cite this version:

Eric Giglio, J. Rangama, S. Guillous, T. Le Cornu. Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets. Physical Review E, 2020, 101 (1), 10.1103/Phys-RevE.101.013105. [https://doi.org/10.1103/Phys-RevE.101.013105](https://doi.org/10.1103/Phys-RevE.101.013105)

HAL Id: hal-03046756

https://hal.science/hal-03046756

Submitted on 21 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

Eric Giglio, J. Rangama, S. Guillous, T. Le Cornu

► To cite this version:

Eric Giglio, J. Rangama, S. Guillous, T. Le Cornu. Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets. Physical Review E, 2020, 101 (1), 10.1103/Phys-RevE.101.013105. [https://doi.org/10.1103/Phys-RevE.101.013105](https://doi.org/10.1103/Phys-RevE.101.013105)

HAL Id: hal-03046756

https://hal.science/hal-03046756v1

Submitted on 21 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

Article in PHYSICAL REVIEW E · January 2020

DOI: 10.1103/PhysRevE.101.013105

CITATIONS

10

READS

215

4 authors, including:

<figure><img src="image_2.png" /><figcaption>Eric Giglio</figcaption></figure>

Centre de Recherche sur les Ions, les Matériaux et la Photonique (CIMAP)

63 PUBLICATIONS

648 CITATIONS

SEE PROFILE

<figure><img src="image_3.png" /><figcaption>Jimmy Rangama</figcaption></figure>

Centre de Recherche sur les Ions, les Matériaux et la Photonique (CIMAP)

129 PUBLICATIONS

1,734 CITATIONS

SEE PROFILE

<figure><img src="image_4.png" /><figcaption>Stéphane Guillous</figcaption></figure>

Atomic Energy and Alternative Energies Commission

45 PUBLICATIONS

483 CITATIONS

SEE PROFILE
# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

Article in PHYSICAL REVIEW E · January 2020

DOI: 10.1103/PhysRevE.101.013105

CITATIONS

10

4 authors, including:

<figure><img src="image_2.png" /><figcaption>Eric Giglio</figcaption></figure>

Centre de Recherche sur les Ions, les Matériaux et la Photonique (CIMAP)

63 PUBLICATIONS

648 CITATIONS

SEE PROFILE

<figure><img src="image_3.png" /><figcaption>Stéphane Guillous</figcaption></figure>

Atomic Energy and Alternative Energies Commission

45 PUBLICATIONS

483 CITATIONS

SEE PROFILE

READS

215

<figure><img src="image_4.png" /><figcaption>Jimmy Rangama</figcaption></figure>

Centre de Recherche sur les Ions, les Matériaux et la Photonique (CIMAP)

129 PUBLICATIONS

1,734 CITATIONS

SEE PROFILE
# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

E. Giglio, J. Rangama, S. Guillous, and T. Le Cornu

Centre de Recherche sur les Ions, les Matériaux et la Photonique,
Université de Caen Normandie, ENSICAEN, CEA, CNRS, 14000 Caen, France

In this work, we model and simulate the shape evolution of critically charged droplets, from the initial spherical shape to the charge emission and back to the spherical shape. The shape deformation is described using the Viscous Correction Viscous Potential Flow (VCVPF) model, which is a potential flow approximation of the Navier-Stokes equation for incompressible Newtonian fluids. The simulated shapes are compared to snapshots of experimentally observed drop deformations. We highlight the influence of the dimensionless viscosity and charge carrier mobility of the liquid on the shape evolution of droplets and discuss the observed trends. We give an explanation as to why the observed deformation pathways of positively and negatively charged pure water droplets differ and give a hint as to why negatively charged water droplets emit more charge during charge break-up than positively charged ones.

PACS numbers: 47.20.Dr, 47.20.Ma

## I. INTRODUCTION

In 1882, Lord Rayleigh predicted that charged droplets become unstable against an infinitesimal spheroidal deformation if the Coulomb pressure at the surface balances the capillary pressure [1]. In a pioneering work, Duft et al. [2, 3] confirmed the Rayleigh prediction of the Coulomb instability of charged droplet by observing the spontaneous shape evolution of critically charged ethylene glycol droplets. By trapping charged micro-droplets in a Paul trap, they succeeded in taking snapshots of the ultra-fast droplet deformation, from the initial spherical shape to the moment jets are emitted from the tips of the elongated drop and back to the spherical shape. Later, the same group also observed the shape evolution of supercooled (-5° C) pure water droplets [4] and heated glycerin droplets [5]. As was already pointed out in [4] and in the PhD work of T. Müller [6], liquid droplets having different viscosities and/or electrical conductivities show a noticeable different shape evolution. The latter is supported by our recent observations of the shape evolution of critically charged de-ionized (pure) water droplets at room temperature. Surprisingly, we found that positively and negatively charged pure water droplets show slightly different shape deformations. That they behave differently is also corroborated by the data found in [6], where negatively charged droplets at room temperature were found to emit almost twice as much charge as positively charged ones.

In this work, we want to understand in more detail how the viscosity and electrical conductivity influence the shape deformation of critically charged droplets. Some simulations of critically charged droplet deformations have already be performed by several authors. We may cite the work of Belelú et al. [7] and Gawande et al. [8] who presented a study of the shape deformation of a charged viscous perfectly conducting drop suspended in air. Using the boundary-element method to solve the

Stokes flow equation (creeping flow approximation), they found that the viscous drop attains an aspect ratio of 3.86 at charge breakup, close to the observed value of 3.85 [2, 3]. However the Stokes flow approximation may break down for the ultra-fast low-viscosity droplet deformation, for which the advective inertial forces may dominate the viscous forces. The case of an inviscid perfectly conducting liquid droplets was discussed in [4] where the simulated pathway was compared to glycol and supercooled water droplets. The authors concluded that viscous forces or finite conductivity are necessary in order to simulate accurately the different experimentally observed droplet deformations. Burton and Taborek [12] extended the study of inviscid liquid drops by considering a finite Ohmic bulk and surface conductivity and showed the influence of both conductivities on the elongation of the drop and on the tip formation. Radcliffe [10, 11] used finite element methods (FEM) to solve the Navier-Stokes (NS) equation and to simulate the deformation of perfectly conducting viscous droplets. His preliminary results showed a strong influence of the viscosity on the elongation of the droplet at the moment pointed tips are formed. Also related to our study is the work of Nganguia et al. [13], who studied the inertia effects on the electrode formation of a viscous drop under a DC electric field as well the nice work of Collins et al. [14] who simulated the jet formation of electrified liquids and discussed the finesse of the jet as a function of the liquid properties. In a later work, their simulations unequivocally show that electrospray daughter droplets are coulombically stable at the instant they are created and that there exists a universal scaling law for the daughter droplets charge [15]. Alternatively, Gañán-Calvo et al. used the promising Volume of Fluid (VoF) method for tracking and locating the free surface (fluid-gas interface) [16]. They succeeded in simulating with high precision the first drop emitted from the jet of an electrified parent droplet and give universal scaling laws for the diameter and electric charge of the first issued daughter droplet,
# Influence of the viscosity and charge mobility on the shape deformation of critically charged droplets

E. Giglio, J. Rangama, S. Guillous, and T. Le Cornu

Centre de Recherche sur les Ions, les Matériaux et la Photonique,
Université de Caen Normandie, ENSICAEN, CEA, CNRS, 14000 Caen, France

In this work, we model and simulate the shape evolution of critically charged droplets, from the initial spherical shape to the charge emission and back to the spherical shape. The shape deformation is described using the Viscous Correction Viscous Potential Flow (VCVPF) model, which is a potential flow approximation of the Navier-Stokes equation for incompressible Newtonian fluids. The simulated shapes are compared to snapshots of experimentally observed drop deformations. We highlight the influence of the dimensionless viscosity and charge carrier mobility of the liquid on the shape evolution of droplets and discuss the observed trends. We give an explanation as to why the observed deformation pathways of positively and negatively charged pure water droplets differ and give a hint as to why negatively charged water droplets emit more charge during charge break-up than positively charged ones.

PACS numbers: 47.20.Dr, 47.20.Ma

## I. INTRODUCTION

In 1882, Lord Rayleigh predicted that charged droplets become unstable against an infinitesimal spheroidal deformation if the Coulomb pressure at the surface balances the capillary pressure [1]. In a pioneering work, Duft et al. [2, 3] confirmed the Rayleigh prediction of the Coulomb instability of charged droplet by observing the spontaneous shape evolution of critically charged ethylene glycol droplets. By trapping charged micro-droplets in a Paul trap, they succeeded in taking snapshots of the ultra-fast droplet deformation, from the initial spherical shape to the moment jets are emitted from the tips of the elongated drop and back to the spherical shape. Later, the same group also observed the shape evolution of supercooled (-5° C) pure water droplets [4] and heated glycerin droplets [5]. As was already pointed out in [4] and in the PhD work of T. Müller [6], liquid droplets having different viscosities and/or electrical conductivities show a noticeable different shape evolution. The latter is supported by our recent observations of the shape evolution of critically charged de-ionized (pure) water droplets at room temperature. Surprisingly, we found that positively and negatively charged pure water droplets show slightly different shape deformations. That they behave differently is also corroborated by the data found in [6], where negatively charged droplets at room temperature were found to emit almost twice as much charge as positively charged ones.

In this work, we want to understand in more detail how the viscosity and electrical conductivity influence the shape deformation of critically charged droplets. Some simulations of critically charged droplet deformations have already be performed by several authors. We may cite the work of Belelú et al. [7] and Gawande et al. [8] who presented a study of the shape deformation of a charged viscous perfectly conducting drop suspended in air. Using the boundary-element method to solve the

Stokes flow equation (creeping flow approximation), they found that the viscous drop attains an aspect ratio of 3.86 at charge breakup, close to the observed value of 3.85 [2, 3]. However the Stokes flow approximation may break down for the ultra-fast low-viscosity droplet deformation, for which the advective inertial forces may dominate the viscous forces. The case of an inviscid perfectly conducting liquid droplets was discussed in [4] where the simulated pathway was compared to glycol and supercooled water droplets. The authors concluded that viscous forces or finite conductivity are necessary in order to simulate accurately the different experimentally observed droplet deformations. Burton and Taborek [12] extended the study of inviscid liquid drops by considering a finite Ohmic bulk and surface conductivity and showed the influence of both conductivities on the elongation of the drop and on the tip formation. Radcliffe [10, 11] used finite element methods (FEM) to solve the Navier-Stokes (NS) equation and to simulate the deformation of perfectly conducting viscous droplets. His preliminary results showed a strong influence of the viscosity on the elongation of the droplet at the moment pointed tips are formed. Also related to our study is the work of Nganguia et al. [13], who studied the inertia effects on the electrode formation of a viscous drop under a DC electric field as well the nice work of Collins et al. [14] who simulated the jet formation of electrified liquids and discussed the finesse of the jet as a function of the liquid properties. In a later work, their simulations unequivocally show that electrospray daughter droplets are coulombically stable at the instant they are created and that there exists a universal scaling law for the daughter droplets charge [15]. Alternatively, Gañán-Calvo et al. used the promising Volume of Fluid (VoF) method for tracking and locating the free surface (fluid-gas interface) [16]. They succeeded in simulating with high precision the first drop emitted from the jet of an electrified parent droplet and give universal scaling laws for the diameter and electric charge of the first issued daughter droplet,
which were validated both numerically and experimentally.

We propose here a numerical model based on the potential flow approximation of the NS equation to simulate the shape evolution of a critically charged droplet. The model accounts for the viscosity and electrical properties of the liquid. In particular, it accounts for the charge mobility of the excess charge carriers at the interface and shows how the liquid properties influence the shape deformation of a critically charged droplet. In section II we describe and discuss our model and highlight the dimensionless quantities that control the shape deformation. In section III we show our experimentally observed snapshots of positively charged pure water droplets and compare the snapshots to simulated ones. We also recall the superellipsoid fit, already introduced in [4], that is used to define the deformation pathway of critically charged droplets. We discuss the influence of the viscosity and charge carrier mobility on the deformation pathway and explain the observed trends. Finally we explain why positively and negatively charged droplets deform differently. It should be emphasized that the interest in charge breakup of droplets is not purely academic as the production of micro-sized or nano-sized jets is of present considerable importance in the charge separation in electrified clouds and in modern technologies such as electro-spray ionization, ink jet printers and fuel injection that require a fine control of the size of the jets and their decay droplets [17].

## II. DYNAMICS OF THE DROPLET DEFORMATION

We want to simulate the nonlinear time-evolution of the droplet shape that develops after the onset of the Coulomb instability of an initially spherical charged drop of radius $R$ in a gas $g$ atmosphere. A liquid of dynamic viscosity $\mu$ and density $\rho$ is supposed electrically conductive. As a result, the charges accumulate at the droplet surface and the electric field is divergence-free $\vec{\nabla} \cdot \vec{E} = 0$ inside and outside the droplet volume. This corresponds to the Taylor-Melcher leaky dielectric model [18, 20]. The droplet deformation is assumed to be axisymmetric and the liquid of the droplet incompressible. The latter condition implies that the velocity vector field of the fluid is divergence-free (solenoidal),

$$ \vec{\nabla} \cdot \vec{u} = 0 \qquad (1) $$

ensuring that the time evolution of the droplet occurs at constant volume. The Navier-Stokes equation for an incompressible (Newtonian) fluid in the droplet volume $\Omega$ is expressed by the equation

$$ \frac{d\vec{u}}{dt} = \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \vec{\nabla} \vec{u} = \frac{1}{\rho} \vec{\nabla} \cdot T \qquad (2) $$

where $T$ is the stress tensor, defined by

$$ T = -pI + 2\mu D + M \qquad (3) $$

with $I$ being the identity matrix, $p$ the pressure in the liquid, $D$ the strain-rate tensor,

$$ D = \frac{1}{2} (\vec{\nabla} \vec{u} + (\vec{\nabla} \vec{u})^t) \qquad (4) $$

and $M$ the Maxwell stress tensor

$$ M_{i,j} = \varepsilon_0 \varepsilon_r E_i E_j - \frac{\varepsilon_0 \varepsilon_r}{2} \vec{E}^2 \delta_{i,j} , \qquad (5) $$

with $\varepsilon_r$ being the dielectric constant of the fluid. Note that $T$ is symmetric, $T = T^t$. As there are no volume charges, the electric force per unit volume is zero, $\vec{\nabla} \cdot M = 0$. Thus, in the leaky dielectric model, the Maxwell tensor is divergence-free.

### Boundary conditions for viscous charged droplets

In the following sections, the equations are defined at the droplet interface, if not specified otherwise. The pressure $p$ in the liquid is obtained via the boundary condition; the jump ($||$) of the normal stress across the liquid-gas interface of the droplet is balanced by the capillary (curvature) force per unit area $P_\gamma$,

$$ T_{nn} = \vec{n} \cdot (T^g - T) \cdot \vec{n} = T_{nn}^g - T_{nn} = P_\gamma \qquad (6) $$

The latter is due to the phenomenon of surface tension $\gamma$, and is given by the Young-Laplace equation

$$ P_{\gamma} = \gamma(\vec{\nabla} \cdot \vec{n}) = -2\gamma H , \qquad (7) $$

where $H$ is the mean curvature at a surface point and $\gamma$ is considered constant. Having defined the normal vector $\vec{n}$ at the drop interface pointing outward, the mean curvature of a sphere of radius $R$ is negative, $H = -1/R$, so that $P_\gamma > 0$. This definition ensures that the pressure in the liquid is higher than the pressure of the atmosphere. We suppose that the viscosity of the outer gas atmosphere is negligible $\mu^g \ll \mu$ and that the pressure of the outer atmosphere $p^g$ is constant (and set to zero), so at the drop interface $\partial\Omega$ one has

$$ \begin{aligned} P_{\gamma} &= T_{nn}^g - T_{nn} \\ &= -(p^g - p) + 2(\mu^g - \mu)D_{nn} + M_{nn}^g - M_{nn} \\ &\approx p - 2\mu D_{nn} + M_{nn}^g - M_{nn} \end{aligned} \qquad (8) $$

We introduce the Coulomb pressure $P_c$, defined as the jump of the normal component of the Maxwell tensor across the interface

$$ \begin{aligned} P_c &= M_{nn}^g - M_{nn} \\ &= \frac{\varepsilon_0}{2} ((E_n^g)^2 - \varepsilon_r E_n^2 + (\varepsilon_r - 1)E_\tau^2) \end{aligned} \qquad (9) $$

where $\vec{\tau}$ stands for the tangential vector at the surface along the polar angle. Injecting the surface charge density $\sigma = \varepsilon_0(E_n^g - \varepsilon_r E_n)$ into (9) allows expressing the Coulomb pressure as a function of $\sigma$ and the electric field
which were validated both numerically and experimentally.

We propose here a numerical model based on the potential flow approximation of the NS equation to simulate the shape evolution of a critically charged droplet. The model accounts for the viscosity and electrical properties of the liquid. In particular, it accounts for the charge mobility of the excess charge carriers at the interface and shows how the liquid properties influence the shape deformation of a critically charged droplet. In section II we describe and discuss our model and highlight the dimensionless quantities that control the shape deformation. In section III we show our experimentally observed snapshots of positively charged pure water droplets and compare the snapshots to simulated ones. We also recall the superellipsoid fit, already introduced in [4], that is used to define the deformation pathway of critically charged droplets. We discuss the influence of the viscosity and charge carrier mobility on the deformation pathway and explain the observed trends. Finally we explain why positively and negatively charged droplets deform differently. It should be emphasized that the interest in charge breakup of droplets is not purely academic as the production of micro-sized or nano-sized jets is of present considerable importance in the charge separation in electrified clouds and in modern technologies such as electro-spray ionization, ink jet printers and fuel injection that require a fine control of the size of the jets and their decay droplets [17].

## II. DYNAMICS OF THE DROPLET DEFORMATION

We want to simulate the nonlinear time-evolution of the droplet shape that develops after the onset of the Coulomb instability of an initially spherical charged drop of radius $R$ in a gas $g$ atmosphere. A liquid of dynamic viscosity $\mu$ and density $\rho$ is supposed electrically conductive. As a result, the charges accumulate at the droplet surface and the electric field is divergence-free $\vec{\nabla} \cdot \vec{E} = 0$ inside and outside the droplet volume. This corresponds to the Taylor-Melcher leaky dielectric model [18, 20]. The droplet deformation is assumed to be axisymmetric and the liquid of the droplet incompressible. The latter condition implies that the velocity vector field of the fluid is divergence-free (solenoidal),

$$ \vec{\nabla} \cdot \vec{u} = 0 \qquad (1) $$

ensuring that the time evolution of the droplet occurs at constant volume. The Navier-Stokes equation for an incompressible (Newtonian) fluid in the droplet volume $\Omega$ is expressed by the equation

$$ \frac{d\vec{u}}{dt} = \frac{\partial \vec{u}}{\partial t} + \vec{u} \cdot \vec{\nabla} \vec{u} = \frac{1}{\rho} \vec{\nabla} \cdot T \qquad (2) $$

where $T$ is the stress tensor, defined by

$$ T = -pI + 2\mu D + M \qquad (3) $$

with $I$ being the identity matrix, $p$ the pressure in the liquid, $D$ the strain-rate tensor,

$$ D = \frac{1}{2} (\vec{\nabla} \vec{u} + (\vec{\nabla} \vec{u})^t) \qquad (4) $$

and $M$ the Maxwell stress tensor

$$ M_{i,j} = \varepsilon_0 \varepsilon_r E_i E_j - \frac{\varepsilon_0 \varepsilon_r}{2} \vec{E}^2 \delta_{i,j} , \qquad (5) $$

with $\varepsilon_r$ being the dielectric constant of the fluid. Note that $T$ is symmetric, $T = T^t$. As there are no volume charges, the electric force per unit volume is zero, $\vec{\nabla} \cdot M = 0$. Thus, in the leaky dielectric model, the Maxwell tensor is divergence-free.

### Boundary conditions for viscous charged droplets

In the following sections, the equations are defined at the droplet interface, if not specified otherwise. The pressure $p$ in the liquid is obtained via the boundary condition; the jump ($||$) of the normal stress across the liquid-gas interface of the droplet is balanced by the capillary (curvature) force per unit area $P_\gamma$,

$$ T_{nn} = \vec{n} \cdot (T^g - T) \cdot \vec{n} = T_{nn}^g - T_{nn} = P_\gamma \qquad (6) $$

The latter is due to the phenomenon of surface tension $\gamma$, and is given by the Young-Laplace equation

$$ P_\gamma = \gamma(\vec{\nabla} \cdot \vec{n}) = -2\gamma H , \qquad (7) $$

where $H$ is the mean curvature at a surface point and $\gamma$ is considered constant. Having defined the normal vector $\vec{n}$ at the drop interface pointing outward, the mean curvature of a sphere of radius $R$ is negative, $H = -1/R$, so that $P_\gamma > 0$. This definition ensures that the pressure in the liquid is higher than the pressure of the atmosphere. We suppose that the viscosity of the outer gas atmosphere is negligible $\mu^g \ll \mu$ and that the pressure of the outer atmosphere $p^g$ is constant (and set to zero), so at the drop interface $\partial\Omega$ one has

$$ \begin{aligned} P_\gamma &= T_{nn}^g - T_{nn} \\ &= -(p^g - p) + 2(\mu^g - \mu)D_{nn} + M_{nn}^g - M_{nn} \\ &\approx p - 2\mu D_{nn} + M_{nn}^g - M_{nn} \end{aligned} \qquad (8) $$

We introduce the Coulomb pressure $P_c$, defined as the jump of the normal component of the Maxwell tensor across the interface

$$ \begin{aligned} P_c &= M_{nn}^g - M_{nn} \\ &= \frac{\varepsilon_0}{2} ((E_n^g)^2 - \varepsilon_r E_n^2 + (\varepsilon_r - 1)E_\tau^2) \end{aligned} \qquad (9) $$

where $\vec{\tau}$ stands for the tangential vector at the surface along the polar angle. Injecting the surface charge density $\sigma = \varepsilon_0(E_n^g - \varepsilon_r E_n)$ into (9) allows expressing the Coulomb pressure as a function of $\sigma$ and the electric field
$(E_n, E_τ)$, evaluated at the liquid side of the liquid-gas interface of the droplet.

$$P_c = \frac{\epsilon_0}{2} \left( \left( \frac{\sigma}{\epsilon_0} + \epsilon_r E_n \right)^2 - \epsilon_r E_n^2 + (\epsilon_r - 1) E_τ^2 \right) \quad (10)$$

Finally rearranging (8) we obtain an equation for the liquid pressure at the interface,

$$p = P_γ - P_c + 2μD_{nn} \quad . \quad (11)$$

The shear stress $T_{τn}$, in contrast, is continuous across the free interface,

$$T_{τn} \parallel \vec{\tau} \cdot (T^g - T) \cdot \vec{n} = 0 \quad (12)$$

$$= 2μ^g D_{τn}^g + M_{τn}^g - (2μD_{τn} + M_{τn}) = 0 \quad (13)$$

We deduce from (13), by considering the viscous shear stress at the gas side of the interface to be negligible, that the viscous shear stress in the liquid must be equal to the Maxwell shear stress across the interface,

$$μ2D_{τn} = M_τ^g - M_{τn} = σE_τ \quad , \quad (14)$$

where the latter equality is deduced from (5). Note that for a perfectly conducting liquid, $E_τ = 0$, yields the more usual condition that the viscous shear stress is zero at the free liquid-gas interface. Finally the prohibition of mass transfer across the interface is guaranteed by the kinematic boundary condition,

$$\left( \frac{d\vec{s}}{dt} - \vec{u}(\vec{s}) \right) \cdot \vec{n} = 0 \quad . \quad (15)$$

where $\vec{u}(\vec{s})$ is the fluid velocity at the boundary represented by $\vec{s}$ and $d\vec{s}/dt$ the velocity of the boundary. Equation (15) describes thus the shape deformation of the droplet.

### Charge dynamics at the interface

We suppose that all excess charges are located within a thin layer (skin) at the interface so that the surface charge approximation holds. The charge dynamics at the interface is thus governed by the surface charge transport equation on a closed deforming surface. The latter is moving with the fluid velocity evaluated at the surface (15). The surface charges are dragged by the liquid flow and pushed by the electric field along the surface. The surface current

$$\vec{j}_s = σ(u_τ + λE_τ)\vec{\tau} \quad (16)$$

has thus two contributions. The first is the velocity field of the fluid tangent to the surface $u_τ$ and the second is the tangent component of the electric field $E_τ$ times the charge mobility $λ$. We may safely neglect the diffusion term which is usually 2 orders of magnitude smaller than

the two former terms. The surface charge evolution eventually reads [19–21],

$$\frac{Dσ}{Dt} - 2Hσu_n + \vec{\nabla}_s \cdot (\vec{\tau}σ(u_τ + λE_τ)) = κE_n \quad (17)$$

where the convective time derivative means the derivative when we follow the surface along a direction normal to itself, that is $\frac{Dσ}{Dt} = \frac{\partialσ}{\partial t} + u_n \frac{\partialσ}{\partial n}$. The second left hand term (LHT), $2Hσu_n$, accounts for the deforming surface. The third LHT involves the surface divergence operator $\vec{\nabla}_s$, which is defined in the supplementary material (36) and accounts for the charge transport along the surface. The last term, proportional to the bulk conductivity $κ$, plays the role of a source term, which moves charges from the bulk to the surface or moves charges from the surface into the bulk. Eventually, the divergence-free electric field inside, $\vec{E}$, and outside, $\vec{E}^{(g)}$, the droplet are deduced by requiring that $\vec{\tau} \cdot (\vec{E}^{(g)} - \vec{E}) = 0$ and $\vec{n} \cdot (\vec{E}^{(g)} - \epsilon_r \vec{E}) = σ/ε_0$ at the interface. Note that (17) differs from the one used in [12], where the drag of charges by the fluid velocity $u_τ$ was omitted and where a constant surface conductivity was used instead of $λσ$.

## A. Viscous Corrected Viscous Potential Flow approximation (VCVPF)

### The VPF approach

We assume in the following that the droplet deformation is well described using the potential flow approximation. The latter implies that the velocity field $\vec{u}$ of the fluid can be written as the gradient of the scalar velocity potential $Ψ$,

$$\vec{u} = \vec{\nabla}Ψ \quad \text{in } Ω \quad . \quad (18)$$

Because of (1), $Ψ$ satisfies the Laplace equation in the droplet volume $Ω$, $\vec{\nabla}^2Ψ = 0$. Within these assumptions, the axisymmetric NS equation simplifies to the 1D Euler equation, giving the time-evolution of the velocity potential $Ψ$ on the interface as a function of the liquid pressure $p$.

$$\frac{\partialΨ}{\partial t} + \frac{\vec{u}^2}{2} + \frac{p}{ρ} = 0 \quad , \quad (19)$$

where the liquid pressure at the interface reads (11),

$$p = P_γ - P_c + 2μD_{nn} \quad , \quad (20)$$

with the difference that here the strain rate tensor $D$ is computed using only the irrotational velocity field $\vec{\nabla}Ψ$. Equations (19) and (20) yield the viscous potential flow model (VPF). Compared to the inviscid liquid model used in [4], the pressure has now an additional term which is proportional to the viscosity $μ$. Note that, while for
$(E_n, E_τ)$, evaluated at the liquid side of the liquid-gas interface of the droplet.

$$P_c = \frac{\varepsilon_0}{2} \left( \left( \frac{\sigma}{\varepsilon_0} + \varepsilon_r E_n \right)^2 - \varepsilon_r E_n^2 + (\varepsilon_r - 1) E_\tau^2 \right) \quad (10)$$

Finally rearranging (8) we obtain an equation for the liquid pressure at the interface,

$$p = P_γ - P_c + 2μD_{nn} \quad . \quad (11)$$

The shear stress $T_{τn}$, in contrast, is continuous across the free interface,

$$T_{τn} = \vec{\tau} \cdot (T^g - T) \cdot \vec{n} = 0 \quad (12)$$

$$= 2μ^g D_{τn}^g + M_{τn}^g - (2μD_{τn} + M_{τn}) = 0 \quad (13)$$

We deduce from (13), by considering the viscous shear stress at the gas side of the interface to be negligible, that the viscous shear stress in the liquid must be equal to the Maxwell shear stress across the interface,

$$μ2D_{τn} = M_τ^g - M_{τn} = σE_τ \quad , \quad (14)$$

where the latter equality is deduced from (5). Note that for a perfectly conducting liquid, $E_τ = 0$, yields the more usual condition that the viscous shear stress is zero at the free liquid-gas interface. Finally the prohibition of mass transfer across the interface is guaranteed by the kinematic boundary condition,

$$\left( \frac{d\vec{s}}{dt} - \vec{u}(\vec{s}) \right) \cdot \vec{n} = 0 \quad . \quad (15)$$

where $\vec{u}(\vec{s})$ is the fluid velocity at the boundary represented by $\vec{s}$ and $d\vec{s}/dt$ the velocity of the boundary. Equation (15) describes thus the shape deformation of the droplet.

### Charge dynamics at the interface

We suppose that all excess charges are located within a thin layer (skin) at the interface so that the surface charge approximation holds. The charge dynamics at the interface is thus governed by the surface charge transport equation on a closed deforming surface. The latter is moving with the fluid velocity evaluated at the surface (15). The surface charges are dragged by the liquid flow and pushed by the electric field along the surface. The surface current

$$\vec{j}_s = σ(u_τ + λE_τ)\vec{\tau} \quad (16)$$

has thus two contributions. The first is the velocity field of the fluid tangent to the surface $u_τ$ and the second is the tangent component of the electric field $E_τ$ times the charge mobility $λ$. We may safely neglect the diffusion term which is usually 2 orders of magnitude smaller than

the two former terms. The surface charge evolution eventually reads [19–21],

$$\frac{D\sigma}{Dt} - 2H\sigma u_n + \vec{\nabla}_s \cdot (\vec{\tau}\sigma(u_τ + \lambda E_τ)) = \kappa E_n \quad (17)$$

where the convective time derivative means the derivative when we follow the surface along a direction normal to itself, that is $\frac{D\sigma}{Dt} = \frac{\partial\sigma}{\partial t} + u_n \frac{\partial\sigma}{\partial n}$. The second left hand term (LHT), $2H\sigma u_n$, accounts for the deforming surface. The third LHT involves the surface divergence operator $\vec{\nabla}_s$, which is defined in the supplementary material (36) and accounts for the charge transport along the surface. The last term, proportional to the bulk conductivity $κ$, plays the role of a source term, which moves charges from the bulk to the surface or moves charges from the surface into the bulk. Eventually, the divergence-free electric field inside, $\vec{E}$, and outside, $\vec{E}^{(g)}$, the droplet are deduced by requiring that $\vec{\tau} \cdot (\vec{E}^{(g)} - \vec{E}) = 0$ and $\vec{n} \cdot (\vec{E}^{(g)} - \varepsilon_r \vec{E}) = \sigma/\varepsilon_0$ at the interface. Note that (17) differs from the one used in [12], where the drag of charges by the fluid velocity $u_τ$ was omitted and where a constant surface conductivity was used instead of $λσ$.

## A. Viscous Corrected Viscous Potential Flow approximation (VCVPF)

### The VPF approach

We assume in the following that the droplet deformation is well described using the potential flow approximation. The latter implies that the velocity field $\vec{u}$ of the fluid can be written as the gradient of the scalar velocity potential $Ψ$,

$$\vec{u} = \vec{\nabla} \Psi \quad \text{in } \Omega \quad . \quad (18)$$

Because of (1), $Ψ$ satisfies the Laplace equation in the droplet volume $Ω$, $\vec{\nabla}^2Ψ = 0$. Within these assumptions, the axisymmetric NS equation simplifies to the 1D Euler equation, giving the time-evolution of the velocity potential $Ψ$ on the interface as a function of the liquid pressure $p$.

$$\frac{\partial Ψ}{\partial t} + \frac{\vec{u}^2}{2} + \frac{p}{\rho} = 0 \quad , \quad (19)$$

where the liquid pressure at the interface reads (11),

$$p = P_γ - P_c + 2μD_{nn} \quad , \quad (20)$$

with the difference that here the strain rate tensor $D$ is computed using only the irrotational velocity field $\vec{\nabla}Ψ$. Equations (19) and (20) yield the viscous potential flow model (VPF). Compared to the inviscid liquid model used in [4], the pressure has now an additional term which is proportional to the viscosity $μ$. Note that, while for
an irrotational flow the viscous force per unit volume is zero inside the drop, $\vec{\nabla} \cdot D = 0$, its contribution to the surface pressure is generally non-zero

$$D_{nn} \neq 0 \text{ and } D_{\tau n} \neq 0 . \qquad (21)$$

As a result, the jump of the shear stress $T_{\tau n}$ across the interface violates the zero shear stress boundary condition at a free surface,

$$T_{\tau n} - M_{\tau n} = \mu 2 D_{\tau n} \neq \sigma E_{\tau} . \qquad (22)$$

Indeed, VPF does not have enough degrees of freedom to enforce that the shear strain rate $D_{\tau n}$ equals the shear stress of the Maxwell tensor. The rotational part of the velocity field is missing here to fulfill the condition (22). In the particular case of a perfectly conducting liquid, where the shear stress of the Maxwell tensor is zero, the shear strain rate $D_{\tau n}$ should vanish at the interface, which is in contradiction with (21). Compared to a rotational flow where the condition (22) can be enforced, VPF results in an incorrect dispersion relationship for small surface oscillations, i.e. the damping of surface modes is underestimated for low viscosities [22].

## The VCVPF model

In order to compensate for the non-zero shear stress boundary condition in VPF, Wang and co-workers [23–25] derived a viscous correction formulation for the irrotational pressure. They assumed that a viscous correction of the pressure at the gas-liquid interface is required to compensate for the non-vanishing irrotational shear stress. The idea is that the pressure correction is a real viscous pressure which varies from the pressure in the irrotational flow outside a narrow vorticity layer near the gas-liquid surface, to the required value at the interface. They called this solution the viscous correction of the VPF model. They proposed thus to add a viscous pressure $p_{\mu}$ to the liquid pressure $p$ of the VPF model,

$$p = P_{\gamma} - P_{c} + 2\mu D_{nn} - p_{\mu} . \qquad (23)$$

The viscous pressure $p_{\mu}$ satisfies the Laplace equation in the drop volume $\Omega$ and is defined such as the power of the viscous pressure balances the power of the non-zero shear stress across the boundary for each surface element $ds^2$,

$$u_n(-p_{\mu})ds^2 + u_{\tau}(-2\mu D_{\tau n} + \sigma E_{\tau})ds^2 = 0 . \qquad (24)$$

The viscous pressure being assumed harmonic, Eq. (24) defines uniquely the viscous potential $p_{\mu}$ in the droplet volume. Combining (19), (23) and (24) yields the time evolution of the velocity potential in the viscous correction for VPF (VCVPF) model. The viscous pressure $p_{\mu}$ adds dissipation at the interface, that originates from the shear stress boundary condition. However, energy dissipation within the bulk is still absent. In the case of

finite conductivity, the tangential component of the electric field may be non-zero, and $p_{\mu}$ accounts also for the traction originating from the non-vanishing shear stress of the Maxwell tensor $\sigma E_{\tau}$ at the interface. The technicalities for the evaluation of the viscous pressure $p_{\mu}$ are detailed in the Appendix Sec. 2. We want to emphasize that the VCVPF model can also be derived from the dissipation model, where no ad hoc viscous pressure is introduced, but where the zero shear stress boundary condition is implicitly satisfied [22].

## Numerical approach

As $\Psi$ satisfies the Laplace equation in the droplet volume, $\vec{\nabla}^2 \Psi = 0$, $\Psi$ can be written as a linear combination of harmonic basis functions, having the symmetry of the problem and the correct behavior close to the origin. We choose to solve the Laplace equation in the prolate spheroidal coordinate system as it is well adapted to our boundary problem, while the harmonic basis functions are still separable. Knowing $\Psi$ on the boundary, the expansion coefficients are readily obtained [4]. This is the same technique we used for evaluating the viscous pressure (see the Appendix Sec. 2). We expand $\Psi$ typically on the first $\leq 100$ harmonic basis functions $f_m$ and the boundary $\vec{s}_i$ is discretized using $i = 300$ grid points distributed so as to cluster close to the tips. The jet formation and emission are not simulated here, as the underlying numerical spectral method, based on harmonic function expansion in prolate coordinate system, is not well suited to solve the Laplace equation for shapes with almost singular tips. Noting that we have also $\vec{\nabla}^2 V = 0$, the harmonic basis expansion technique is also used to evaluate the electric field inside and outside the droplet.

## Validity of VCVPF

We note that VCVPF only approaches the power of dissipation of the NS equation [23], mainly because the rotational part of the velocity field is missing, which is also necessary to satisfy explicitly the shear stress condition (13) at the free surface. Padrino et al. [22] studied the dispersion relationships, giving the damping of small amplitude surface oscillations as a function of the surface mode, for both models, VPF and VCVPF and compared them to the exact NS solution. They conclude that, while VPF underestimates the damping of the surface oscillations, VCVPF tends to overestimate the damping with respect to the exact NS solution. However, for low viscous liquids, $\mu < 10^{-2}\sqrt{\rho\gamma R}$ (cf. 27), VCVPF does recover the exact dispersion relationship, while VPF underestimated the damping by at least a factor 2 compared to VCVPF. We prefer thus VCVPF over VPF for the study of shape deformation of low viscous liquid drops. We also want to highlight that unlike the Stokes equation, which was used in [7, 8] to study the shape deformation
an irrotational flow the viscous force per unit volume is zero inside the drop, $\vec{\nabla} \cdot D = 0$, its contribution to the surface pressure is generally non-zero

$$D_{nn} \neq 0 \quad \text{and} \quad D_{\tau n} \neq 0 . \qquad (21)$$

As a result, the jump of the shear stress $T_{\tau n}$ across the interface violates the zero shear stress boundary condition at a free surface,

$$T_{\tau n} - M_{\tau n} = \mu 2 D_{\tau n} \neq \sigma E_{\tau} . \qquad (22)$$

Indeed, VPF does not have enough degrees of freedom to enforce that the shear strain rate $D_{\tau n}$ equals the shear stress of the Maxwell tensor. The rotational part of the velocity field is missing here to fulfill the condition (22). In the particular case of a perfectly conducting liquid, where the shear stress of the Maxwell tensor is zero, the shear strain rate $D_{\tau n}$ should vanish at the interface, which is in contradiction with (21). Compared to a rotational flow where the condition (22) can be enforced, VPF results in an incorrect dispersion relationship for small surface oscillations, i.e. the damping of surface modes is underestimated for low viscosities [22].

### The VCVPF model

In order to compensate for the non-zero shear stress boundary condition in VPF, Wang and co-workers [23–25] derived a viscous correction formulation for the irrotational pressure. They assumed that a viscous correction of the pressure at the gas-liquid interface is required to compensate for the non-vanishing irrotational shear stress. The idea is that the pressure correction is a real viscous pressure which varies from the pressure in the irrotational flow outside a narrow vorticity layer near the gas-liquid surface, to the required value at the interface. They called this solution the viscous correction of the VPF model. They proposed thus to add a viscous pressure $p_{\mu}$ to the liquid pressure $p$ of the VPF model,

$$p = P_{\gamma} - P_{c} + 2\mu D_{nn} - p_{\mu} . \qquad (23)$$

The viscous pressure $p_{\mu}$ satisfies the Laplace equation in the drop volume $\Omega$ and is defined such as the power of the viscous pressure balances the power of the non-zero shear stress across the boundary for each surface element $ds^2$,

$$u_n(-p_{\mu})ds^2 + u_{\tau}(-2\mu D_{\tau n} + \sigma E_{\tau})ds^2 = 0 . \qquad (24)$$

The viscous pressure being assumed harmonic, Eq. (24) defines uniquely the viscous potential $p_{\mu}$ in the droplet volume. Combining (19), (23) and (24) yields the time evolution of the velocity potential in the viscous correction for VPF (VCVPF) model. The viscous pressure $p_{\mu}$ adds dissipation at the interface, that originates from the shear stress boundary condition. However, energy dissipation within the bulk is still absent. In the case of

finite conductivity, the tangential component of the electric field may be non-zero, and $p_{\mu}$ accounts also for the traction originating from the non-vanishing shear stress of the Maxwell tensor $\sigma E_{\tau}$ at the interface. The technicalities for the evaluation of the viscous pressure $p_{\mu}$ are detailed in the Appendix Sec. 2. We want to emphasize that the VCVPF model can also be derived from the dissipation model, where no ad hoc viscous pressure is introduced, but where the zero shear stress boundary condition is implicitly satisfied [22].

### Numerical approach

As $\Psi$ satisfies the Laplace equation in the droplet volume, $\vec{\nabla}^2\Psi = 0$, $\Psi$ can be written as a linear combination of harmonic basis functions, having the symmetry of the problem and the correct behavior close to the origin. We choose to solve the Laplace equation in the prolate spheroidal coordinate system as it is well adapted to our boundary problem, while the harmonic basis functions are still separable. Knowing $\Psi$ on the boundary, the expansion coefficients are readily obtained [4]. This is the same technique we used for evaluating the viscous pressure (see the Appendix Sec. 2). We expand $\Psi$ typically on the first $\leq 100$ harmonic basis functions $f_m$ and the boundary $\vec{s}_i$ is discretized using $i = 300$ grid points distributed so as to cluster close to the tips. The jet formation and emission are not simulated here, as the underlying numerical spectral method, based on harmonic function expansion in prolate coordinate system, is not well suited to solve the Laplace equation for shapes with almost singular tips. Noting that we have also $\vec{\nabla}^2V = 0$, the harmonic basis expansion technique is also used to evaluate the electric field inside and outside the droplet.

### Validity of VCVPF

We note that VCVPF only approaches the power of dissipation of the NS equation [23], mainly because the rotational part of the velocity field is missing, which is also necessary to satisfy explicitly the shear stress condition (13) at the free surface. Padrino et al. [22] studied the dispersion relationships, giving the damping of small amplitude surface oscillations as a function of the surface mode, for both models, VPF and VCVPF and compared them to the exact NS solution. They conclude that, while VPF underestimates the damping of the surface oscillations, VCVPF tends to overestimate the damping with respect to the exact NS solution. However, for low viscous liquids, $\mu < 10^{-2}\sqrt{\rho\gamma R}$ (cf. 27), VCVPF does recover the exact dispersion relationship, while VPF underestimated the damping by at least a factor 2 compared to VCVPF. We prefer thus VCVPF over VPF for the study of shape deformation of low viscous liquid drops. We also want to highlight that unlike the Stokes equation, which was used in [7, 8] to study the shape deformation
of critically charged droplets, the VPF and VCVPF models include both the advective term $\rho \frac{\ddot{u}^2}{2}$, accounting thus also for inertial effects.

## B. Dimensionless quantities

We introduce the characteristic time $\tau_0$, giving the timescale of the deformation, and the characteristic capillary pressure $P_\gamma^\circ$ of the spherical drop,

$$ \tau_0 = \sqrt{\frac{\rho R^3}{\gamma}} \quad , \quad P_\gamma^\circ = 2\gamma/R \quad . \quad (25) $$

Defining the dimensionless time $\tilde{t} = t/\tau_0$, velocity potential $\tilde{\Psi} = \Psi \frac{\tau_0}{R^2}$, capillary pressure $-\tilde{H} = P_\gamma/P_\gamma^\circ$ and viscous pressure $\tilde{p}_\mu = p_\mu/P_\gamma^\circ$, allows putting the time evolution of the velocity potential $\Psi$ in VCVPF into a dimensionless from

$$ \frac{\partial \tilde{\Psi}}{\partial \tilde{t}} + \frac{\tilde{u}^2}{2} + 2(-\tilde{H} - X\tilde{P}_c + \text{Oh}\tilde{D}_{nn} - \tilde{p}_\mu) = 0, \quad (26) $$

where the quantities with a tilde are dimensionless. We introduced in (26) the dimensionless viscosity $\text{Oh}$ (Ohne-sorge number),

$$ \text{Oh} = \frac{\mu}{\sqrt{\rho\gamma R}} \quad (27) $$

which relates the viscous forces to inertial and surface tension forces. Larger $\text{Oh}$ indicate a greater influence of the viscosity. Note that $\tilde{p}_\mu$ depends on $\text{Oh}$ via (24). We further introduced the fissility parameter

$$ X = \frac{P_c^\circ}{P_\gamma^\circ} = \frac{Q^2}{64\pi^2\varepsilon_0\gamma R^3} \quad , \quad (28) $$

defined as the ratio of the Coulomb pressure over the surface pressure for the spherical droplet shape. If $X < 1$ the spherical shape of droplet is stable. If $X = 1$, the spherical droplet is unstable and deforms spontaneously. The dimensionless Coulomb pressure $\tilde{P}_c$ found in (26) is defined by $P_c = \tilde{P}_c P_c^\circ = X \tilde{P}_c P_\gamma^\circ$, yielding

$$ \tilde{P}_c = (\tilde{\sigma} + \varepsilon_r \tilde{E}_n)^2 - \varepsilon_r \tilde{E}_n^2 + (\varepsilon_r - 1) \tilde{E}_\tau \quad , \quad (29) $$

where $\tilde{\sigma}$ and $\tilde{E}$ are given by the relations $\sigma = \tilde{\sigma}\sqrt{4X\varepsilon_0\gamma/R}$ and $E = \tilde{E}\sqrt{4X\gamma/(\varepsilon_0 R)}$.

The surface charge dynamics at the interface is governed by equation (17). Introducing the surface and bulk relaxation rates

$$ \tau_s^{-1} = \lambda \sqrt{\frac{X\gamma}{R^3\varepsilon_0}} \quad , \quad \tau_b^{-1} = \frac{\kappa}{\varepsilon_0\varepsilon_r} \quad (30) $$

allows to define the dimensionless charge carrier mobility $\tilde{\lambda}$ and bulk conductivity $\tilde{\kappa}$

$$ \tilde{\lambda} = \frac{\tau_0}{\tau_s} = \lambda \sqrt{\frac{4X\rho}{\varepsilon_0}} \quad (31) $$

$$ \tilde{\kappa} = \frac{\tau_0}{\tau_b} = \kappa \sqrt{\frac{\rho R^3}{\gamma \varepsilon_0^2 \varepsilon_r^2}} \quad (32) $$

Note that $\tilde{\lambda}$ is proportional to $\sqrt{X}$ and thus to the charge $|Q|$ carried by the droplet. Equation (17) can now be put into the dimensionless form

$$ \frac{D\tilde{\sigma}}{D\tilde{t}} - 2\tilde{H}\tilde{\sigma}\tilde{u}_n + \vec{\nabla}_{\tilde{s}} \cdot (\vec{\tau}\tilde{\sigma}(\tilde{u}_\tau + \tilde{\lambda}\tilde{E}_\tau)) = \tilde{\kappa}\varepsilon_r\tilde{E}_n \quad . \quad (33) $$

Note that the factor $\varepsilon_r$, found in the last term, simply originates from the definition of the bulk relaxation time $\tau_b$ (30). The charge dynamics at the interface is controlled by two parameters namely $\tilde{\lambda}$ and $\tilde{\kappa}$. If $\tilde{\lambda} \gg 1$ or $\tilde{\kappa} \gg 1$, then the charge relaxation dominates the rate of motion of the fluid and the surface charge distribution is close to equilibrium, meaning that the electric potential at the surface is an equipotential. We recover thus the case of a perfectly conducting liquid, where the surface charge distribution is always at equilibrium during the droplet deformation, $\tilde{\sigma} = \tilde{\sigma}_{eq}$, so that $\tilde{P}_c = \tilde{\sigma}_{eq}^2$ depends merely on the droplet shape $\tilde{s}$. From equation (26) we deduce that, in the particular case of a perfectly conducting liquid, the shape deformation depends only on the dimensionless viscosity $\text{Oh}$ and on the fissility $X$. If the droplet is critically charged, $X = 1$, the shape deformation depends only on the dimensionless viscosity $\text{Oh}$.

However, for liquid drops with sufficiently low charge relaxation rates, so that $\tilde{\lambda} \lesssim 1$ and bulk conductivity $\tilde{\kappa} \lesssim 1$, the surface charge distribution may not be at equilibrium during the ultra-fast droplet deformation (in the sense that the potential at the surface is no longer an equipotential). This is because the surface charges can now be dragged by the liquid flow out of equilibrium. As a result, the inner electric field is non-zero at the interface and the dimensionless Coulomb pressure $\tilde{P}_c$ deviates from the Coulomb pressure of a perfect conductor. Hence, in the case of low conducting liquids, the shape deformation of the droplet depends not only on $\text{Oh}$ and $X$ but also indirectly on the parameters $\tilde{\lambda}$ and $\tilde{\kappa}$.

## III. SHAPE DEFORMATION OF CRITICALLY CHARGED DROPLETS

We took snapshots of the ultrafast shape deformation of pure water droplets at the Rayleigh limit. The experimental setup used for trapping and taking snapshots of the droplet is similar to the one presented in [3–5]. A charged droplet is injected into a Paul trap, where it evaporates until it reaches the critical fissility $X = 1$. The onset of the Coulomb instability is monitored by a photo-diode which triggers the flash lamp and the CCD
of critically charged droplets, the VPF and VCVPF models include both the advective term $\rho \frac{\ddot{u}^2}{2}$, accounting thus also for inertial effects.

## B. Dimensionless quantities

We introduce the characteristic time $\tau_0$, giving the timescale of the deformation, and the characteristic capillary pressure $P_\gamma^\circ$ of the spherical drop,

$$ \tau_0 = \sqrt{\frac{\rho R^3}{\gamma}} \quad , \quad P_\gamma^\circ = 2\gamma/R \quad . \quad (25) $$

Defining the dimensionless time $\tilde{t} = t/\tau_0$, velocity potential $\tilde{\Psi} = \Psi \frac{\tau_0}{R^2}$, capillary pressure $-\tilde{H} = P_\gamma/P_\gamma^\circ$ and viscous pressure $\tilde{p}_\mu = p_\mu/P_\gamma^\circ$, allows putting the time evolution of the velocity potential $\Psi$ in VCVPF into a dimensionless from

$$ \frac{\partial \tilde{\Psi}}{\partial \tilde{t}} + \frac{\tilde{u}^2}{2} + 2(-\tilde{H} - X \tilde{P}_c + \text{Oh} \tilde{D}_{nn} - \tilde{p}_\mu) = 0, \quad (26) $$

where the quantities with a tilde are dimensionless. We introduced in (26) the dimensionless viscosity $\text{Oh}$ (Ohne-sorge number),

$$ \text{Oh} = \frac{\mu}{\sqrt{\rho\gamma R}} \quad (27) $$

which relates the viscous forces to inertial and surface tension forces. Larger $\text{Oh}$ indicate a greater influence of the viscosity. Note that $\tilde{p}_\mu$ depends on $\text{Oh}$ via (24). We further introduced the fissility parameter

$$ X = \frac{P_c^\circ}{P_\gamma^\circ} = \frac{Q^2}{64\pi^2\varepsilon_0\gamma R^3} \quad , \quad (28) $$

defined as the ratio of the Coulomb pressure over the surface pressure for the spherical droplet shape. If $X < 1$ the spherical shape of droplet is stable. If $X = 1$, the spherical droplet is unstable and deforms spontaneously. The dimensionless Coulomb pressure $\tilde{P}_c$ found in (26) is defined by $P_c = \tilde{P}_c P_c^\circ = X \tilde{P}_c P_\gamma^\circ$, yielding

$$ \tilde{P}_c = (\tilde{\sigma} + \varepsilon_r \tilde{E}_n)^2 - \varepsilon_r \tilde{E}_n^2 + (\varepsilon_r - 1) \tilde{E}_\tau \quad , \quad (29) $$

where $\tilde{\sigma}$ and $\tilde{E}$ are given by the relations $\sigma = \tilde{\sigma}\sqrt{4X\varepsilon_0\gamma/R}$ and $E = \tilde{E}\sqrt{4X\gamma/(\varepsilon_0 R)}$.

The surface charge dynamics at the interface is governed by equation (17). Introducing the surface and bulk relaxation rates

$$ \tau_s^{-1} = \lambda \sqrt{\frac{X\gamma}{R^3\varepsilon_0}} \quad , \quad \tau_b^{-1} = \frac{\kappa}{\varepsilon_0\varepsilon_r} \quad (30) $$

allows to define the dimensionless charge carrier mobility $\tilde{\lambda}$ and bulk conductivity $\tilde{\kappa}$

$$ \tilde{\lambda} = \frac{\tau_0}{\tau_s} = \lambda \sqrt{\frac{4X\rho}{\varepsilon_0}} \quad (31) $$

$$ \tilde{\kappa} = \frac{\tau_0}{\tau_b} = \kappa \sqrt{\frac{\rho R^3}{\gamma \varepsilon_0^2 \varepsilon_r^2}} \quad (32) $$

Note that $\tilde{\lambda}$ is proportional to $\sqrt{X}$ and thus to the charge $|Q|$ carried by the droplet. Equation (17) can now be put into the dimensionless form

$$ \frac{D\tilde{\sigma}}{D\tilde{t}} - 2\tilde{H}\tilde{\sigma}\tilde{u}_n + \vec{\nabla}_{\tilde{s}} \cdot (\vec{\tau}\tilde{\sigma}(\tilde{u}_\tau + \tilde{\lambda}\tilde{E}_\tau)) = \tilde{\kappa}\varepsilon_r\tilde{E}_n \quad . \quad (33) $$

Note that the factor $\varepsilon_r$, found in the last term, simply originates from the definition of the bulk relaxation time $\tau_b$ (30). The charge dynamics at the interface is controlled by two parameters namely $\tilde{\lambda}$ and $\tilde{\kappa}$. If $\tilde{\lambda} \gg 1$ or $\tilde{\kappa} \gg 1$, then the charge relaxation dominates the rate of motion of the fluid and the surface charge distribution is close to equilibrium, meaning that the electric potential at the surface is an equipotential. We recover thus the case of a perfectly conducting liquid, where the surface charge distribution is always at equilibrium during the droplet deformation, $\tilde{\sigma} = \tilde{\sigma}_{eq}$, so that $\tilde{P}_c = \tilde{\sigma}_{eq}^2$ depends merely on the droplet shape $\tilde{s}$. From equation (26) we deduce that, in the particular case of a perfectly conducting liquid, the shape deformation depends only on the dimensionless viscosity $\text{Oh}$ and on the fissility $X$. If the droplet is critically charged, $X = 1$, the shape deformation depends only on the dimensionless viscosity $\text{Oh}$.

However, for liquid drops with sufficiently low charge relaxation rates, so that $\tilde{\lambda} \lesssim 1$ and bulk conductivity $\tilde{\kappa} \lesssim 1$, the surface charge distribution may not be at equilibrium during the ultra-fast droplet deformation (in the sense that the potential at the surface is no longer an equipotential). This is because the surface charges can now be dragged by the liquid flow out of equilibrium. As a result, the inner electric field is non-zero at the interface and the dimensionless Coulomb pressure $\tilde{P}_c$ deviates from the Coulomb pressure of a perfect conductor. Hence, in the case of low conducting liquids, the shape deformation of the droplet depends not only on $\text{Oh}$ and $X$ but also indirectly on the parameters $\tilde{\lambda}$ and $\tilde{\kappa}$.

## III. SHAPE DEFORMATION OF CRITICALLY CHARGED DROPLETS

We took snapshots of the ultrafast shape deformation of pure water droplets at the Rayleigh limit. The experimental setup used for trapping and taking snapshots of the droplet is similar to the one presented in [3–5]. A charged droplet is injected into a Paul trap, where it evaporates until it reaches the critical fissility $X = 1$. The onset of the Coulomb instability is monitored by a photo-diode which triggers the flash lamp and the CCD
<figure><img src="image_5.png" /><figcaption>FIG. 1: Comparison between simulated (top) and observed (bottom) snapshots of $\emptyset = 52\mu m$ deionized water droplets at room temperature. Initially, the red droplets are charged at the Rayleigh limit $X = 1$. At charge breakup, 20 % of the charge is removed and the fissility for the green droplet is $X = 0.64$. Shape labels: (a) lemon, (b) peanut, (c) diamond.</figcaption></figure>

camera. Delaying the trigger allows taking snapshots of the time evolution of the shape deformation. At the onset of the instability, the trapped water droplets have a $\emptyset = 52\mu m$. Assuming that water at room temperature has a viscosity of 1 cP, a density of $1000 \text{ kg/m}^3$ and a surface tension of $0.072 \text{ N/m}$, the droplets are characterized by a low dimensionless viscosity of $\text{Oh}=0.023$ (see table I), and are thus particularly good candidates for the simulation with our VCVPF code. Surprisingly no jets were detected on the snapshots, neither for positively nor for negatively charged water droplets at room temperature. We took more than 30 snapshots of the moment where the elongated droplet shapes have pointed ends, but found no trace of small daughter droplets or jets on the snapshots, unlike for glycol and supercooled water [3, 4]. We concluded that the jet was too thin to be observed and estimated that the jet radius, if any, must be smaller than the resolution of the snapshots of $0.36 \mu m$. These observations are corroborated by the experimental results presented in the Ph.D work of R. Müller [6]. The lack of observed jets made us assume that the kinetic energy loss during charge emission should be negligible. We also assumed that the charge is emitted in an amount of time short enough to consider the droplet deformation frozen during charge emission. In our numerical code, the charge emission is simulated thus by simply removing instantaneously a fraction of the charge as soon as pointed tips are formed. The kinetic energy accumulated during the first part of the deformation is conserved during the charge emission. In the example shown in Fig. 1 the deionized water droplets were positively charged. The measurements in [6] indicate that positively charged pure water droplets emit at room temperature about 20

<table><caption>TABLE I: Parameters of deionized (pure) water droplets at room temperature.</caption><thead><tr><th>R (μm)</th><th>μ (cP)</th><th>ρ (kg/m³)</th><th>γ (mN/m)</th><th>Oh</th><th>τ₀ (μs)</th></tr></thead><tbody><tr><td>26</td><td>1</td><td>1000</td><td>72</td><td>0.023</td><td>16</td></tr></tbody></table>

% of the initial charge at charge breakup.

We simulated with our numerical code, based on the VCVPF model, the complete shape deformation of an initially critically charged droplet, including the charge breakup and the "return" to the stable sub-critically charged spherical drop. We used in our simulations a bulk conductivity of $\kappa_b = 5.5 \times 10^{-6} \text{ S/m}$ and a charge mobility of $\lambda = 37 \times 10^{-8} \text{ m}^2/\text{s/V}$ (see table II). We compare in Fig. 1 the simulated shapes to the observed ones. The red simulated shapes in Fig. 1 are critically charged droplets ($X = 1$) and show the ultra-fast shape deformation after the onset of the instability. Once the pointed tips are formed, we removed, as suggested by [6], 20% of the charge and the fissility reduces instantaneously to $X = 0.64$. At that moment, the shape fits remarkably well the experimentally found lemon-like shape. But even more remarkably is the agreement between simulated and observed shapes during the sub-critical shape deformation (green shapes). Especially, the peanut and diamond shape have been found by the simulation. The general agreement validates our approach for describing low-viscous liquid droplets.

**Superellipsoid:** In order to compare easily the deformation pathways of different droplets, the shapes are fitted with the mathematical expression of a superellipsoid of the form [30]

$$ \left| \frac{z}{a} \right|^n + \left| \frac{\sqrt{x^2 + y^2}}{b} \right|^n = 1 \quad . \qquad (34) $$

The shapes are fitted with the constraint of constant volume so that the fit had only two independent parameters, namely $a/b$ and $n$. Expression (34) can fit either capsules, spheroids or pointed shapes, depending on $n$. The aspect ratio $a/b$ gives the elongation of the droplet. Adjusting $(a/b, n)$ for the droplet shape for each snapshot yields a number of points that can be plotted on a graph with $a/b$ as the abscissa and $n$ as the ordinate, so that the complete droplet deformation path can be represented reasonably well on a 2D graph. Fitting the shapes by a superellipsoid allows us to compare the deformation pathway of droplets having different liquid properties. Such a graphical representation of the deformation pathway has already been used in [4] to highlight the differences in the shape deformation of glycol and supercooled water droplets or in [10, 11] to show the influence of the viscosity on the deformation pathway. Below, we give examples of the fit parameters $(a/b, n)$ found for various observed shapes. We took $\sim 270$ snapshots of the fast deforming water droplet and fitted the shapes on the snapshots by the superellipsoid function. The resulting set of
<figure><img src="image_5.png" /><figcaption>FIG. 1: Comparison between simulated (top) and observed (bottom) snapshots of $\emptyset = 52\mu m$ deionized water droplets at room temperature. Initially, the red droplets are charged at the Rayleigh limit $X = 1$. At charge breakup, 20 % of the charge is removed and the fissility for the green droplet is $X = 0.64$. Shape labels: (a) lemon, (b) peanut, (c) diamond.</figcaption></figure>

camera. Delaying the trigger allows taking snapshots of the time evolution of the shape deformation. At the onset of the instability, the trapped water droplets have a $\emptyset = 52\mu m$. Assuming that water at room temperature has a viscosity of 1 cP, a density of $1000 \text{ kg/m}^3$ and a surface tension of $0.072 \text{ N/m}$, the droplets are characterized by a low dimensionless viscosity of Oh=0.023 (see table I), and are thus particularly good candidates for the simulation with our VCVPF code. Surprisingly no jets were detected on the snapshots, neither for positively nor for negatively charged water droplets at room temperature. We took more than 30 snapshots of the moment where the elongated droplet shapes have pointed ends, but found no trace of small daughter droplets or jets on the snapshots, unlike for glycol and supercooled water [3, 4]. We concluded that the jet was too thin to be observed and estimated that the jet radius, if any, must be smaller than the resolution of the snapshots of 0.36 $\mu m$. These observations are corroborated by the experimental results presented in the Ph.D work of R. Müller [6]. The lack of observed jets made us assume that the kinetic energy loss during charge emission should be negligible. We also assumed that the charge is emitted in an amount of time short enough to consider the droplet deformation frozen during charge emission. In our numerical code, the charge emission is simulated thus by simply removing instantaneously a fraction of the charge as soon as pointed tips are formed. The kinetic energy accumulated during the first part of the deformation is conserved during the charge emission. In the example shown in Fig. 1 the deionized water droplets were positively charged. The measurements in [6] indicate that positively charged pure water droplets emit at room temperature about 20

<table><caption>TABLE I: Parameters of deionized (pure) water droplets at room temperature.</caption><thead><tr><th>R (μm)</th><th>μ (cP)</th><th>ρ (kg/m³)</th><th>γ (mN/m)</th><th>Oh</th><th>τ₀ (μs)</th></tr></thead><tbody><tr><td>26</td><td>1</td><td>1000</td><td>72</td><td>0.023</td><td>16</td></tr></tbody></table>

% of the initial charge at charge breakup.

We simulated with our numerical code, based on the VCVPF model, the complete shape deformation of an initially critically charged droplet, including the charge breakup and the "return" to the stable sub-critically charged spherical drop. We used in our simulations a bulk conductivity of $\kappa_b = 5.5 \times 10^{-6} \text{ S/m}$ and a charge mobility of $\lambda = 37 \times 10^{-8} \text{ m}^2/\text{s/V}$ (see table II). We compare in Fig. 1 the simulated shapes to the observed ones. The red simulated shapes in Fig. 1 are critically charged droplets ($X = 1$) and show the ultra-fast shape deformation after the onset of the instability. Once the pointed tips are formed, we removed, as suggested by [6], 20% of the charge and the fissility reduces instantaneously to $X = 0.64$. At that moment, the shape fits remarkably well the experimentally found lemon-like shape. But even more remarkably is the agreement between simulated and observed shapes during the sub-critical shape deformation (green shapes). Especially, the peanut and diamond shape have been found by the simulation. The general agreement validates our approach for describing low-viscous liquid droplets.

**Superellipsoid:** In order to compare easily the deformation pathways of different droplets, the shapes are fitted with the mathematical expression of a superellipsoid of the form [30]

$$ \left| \frac{z}{a} \right|^n + \left| \frac{\sqrt{x^2 + y^2}}{b} \right|^n = 1 \quad . \qquad (34) $$

The shapes are fitted with the constraint of constant volume so that the fit had only two independent parameters, namely $a/b$ and $n$. Expression (34) can fit either capsules, spheroids or pointed shapes, depending on $n$. The aspect ratio $a/b$ gives the elongation of the droplet. Adjusting $(a/b, n)$ for the droplet shape for each snapshot yields a number of points that can be plotted on a graph with $a/b$ as the abscissa and $n$ as the ordinate, so that the complete droplet deformation path can be represented reasonably well on a 2D graph. Fitting the shapes by a superellipsoid allows us to compare the deformation pathway of droplets having different liquid properties. Such a graphical representation of the deformation pathway has already been used in [4] to highlight the differences in the shape deformation of glycol and supercooled water droplets or in [10, 11] to show the influence of the viscosity on the deformation pathway. Below, we give examples of the fit parameters $(a/b, n)$ found for various observed shapes. We took $\sim 270$ snapshots of the fast deforming water droplet and fitted the shapes on the snapshots by the superellipsoid function. The resulting set of
<figure><img src="image_6.png" /><figcaption>FIG. 2: Fit parameters (a/b, n) found for various observed shapes using the superellipsoid equation (34). Red contours show the corresponding superellipsoid shape.</figcaption></figure>

(a/b, n) was reported on a graph, yielding the pathway of the droplet deformation, see Fig. 3. The arrows give the direction of the time evolution of the path. It should be noted that the scattering in the experimental data shown in the (a/b,n) graphs results mainly from the qualitative fit of the snapshots by the superellipsoid equation. The fit was done by software and then the values (a/b,n) were verified for consistency by the authors. Similar shapes can give some dispersion in the fitted parameters especially for n > 2.5 and n < 1.7. In the following, we will use the deformation pathway as an analysis tool to study the influence of the liquid properties on the droplet shape deformation.

<figure><img src="image_7.png" /><figcaption>FIG. 3: Graph (a/b, n) showing the deformation pathway of a critically charged deionized water droplet. Starting with the spherical shape (1, 2), arrows give the time evolution of the path.</figcaption></figure>

### A. Influence of the viscosity

In a first step, we focus on the influence of the dimensionless viscosity Oh on the deformation pathway. We assume here that the fluid is perfectly conducting, so that Oh is the only free parameter controlling the deformation

pathway. We simulated for different Ohnesorge numbers Oh, ranging from 0.001 to 0.8, the shape deformation of critically charged droplets ($X = 1$). The pathways for several values of Oh are shown in Fig. 4. In all simulations, the initial shape was a spheroid $n = 2$ with a ratio $a/b = 1.03$ (almost spherical). First, we note that the general deformation pathway is similar for all cases: the instability starts with $a/b = 1.03$ and $n = 2$. Then $a/b$ increases, while $n$ remains nearly constant at $\simeq 2$, meaning that the initial deformation remains mainly prolate spheroidal. With ongoing elongation of the drop, $n$ starts to decrease and the shape becomes lemon-like with pointed tips.

<figure><img src="image_8.png" /><figcaption>FIG. 4: Deformation pathway of a critically charged droplet as a function of the Ohnesorge number Oh. At the pointed tip appearance ($n = 1.5$), the fissility is reduced from $X = 1$ to $X = 0.64$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5: Elongation $a/b$ of the droplet at charge breakup, (taken at $n = 1.5$) as a function of the dimensionless viscosity Oh. Purple squares stand for the VCVPF model and green empty circles for the VPF model. The horizontal dotted black line is $a/b = 3.86$ and gives the elongation of the droplet in the case of the Stokes flow approximation [7, 8].</figcaption></figure>

When looking in more detail, we find that with in-
<figure><img src="image_6.png" /><figcaption>FIG. 2: Fit parameters (a/b, n) found for various observed shapes using the superellipsoid equation (34). Red contours show the corresponding superellipsoid shape.</figcaption></figure>

(a/b, n) was reported on a graph, yielding the pathway of the droplet deformation, see Fig. 3. The arrows give the direction of the time evolution of the path. It should be noted that the scattering in the experimental data shown in the (a/b,n) graphs results mainly from the qualitative fit of the snapshots by the superellipsoid equation. The fit was done by software and then the values (a/b,n) were verified for consistency by the authors. Similar shapes can give some dispersion in the fitted parameters especially for n > 2.5 and n < 1.7. In the following, we will use the deformation pathway as an analysis tool to study the influence of the liquid properties on the droplet shape deformation.

<figure><img src="image_7.png" /><figcaption>FIG. 3: Graph (a/b, n) showing the deformation pathway of a critically charged deionized water droplet. Starting with the spherical shape (1, 2), arrows give the time evolution of the path.</figcaption></figure>

### A. Influence of the viscosity

In a first step, we focus on the influence of the dimensionless viscosity Oh on the deformation pathway. We assume here that the fluid is perfectly conducting, so that Oh is the only free parameter controlling the deformation

pathway. We simulated for different Ohnesorge numbers Oh, ranging from 0.001 to 0.8, the shape deformation of critically charged droplets ($X = 1$). The pathways for several values of Oh are shown in Fig. 4. In all simulations, the initial shape was a spheroid $n = 2$ with a ratio $a/b = 1.03$ (almost spherical). First, we note that the general deformation pathway is similar for all cases: the instability starts with $a/b = 1.03$ and $n = 2$. Then $a/b$ increases, while $n$ remains nearly constant at $\simeq 2$, meaning that the initial deformation remains mainly prolate spheroidal. With ongoing elongation of the drop, $n$ starts to decrease and the shape becomes lemon-like with pointed tips.

<figure><img src="image_8.png" /><figcaption>FIG. 4: Deformation pathway of a critically charged droplet as a function of the Ohnesorge number Oh. At the pointed tip appearance ($n = 1.5$), the fissility is reduced from $X = 1$ to $X = 0.64$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. 5: Elongation $a/b$ of the droplet at charge breakup, (taken at $n = 1.5$) as a function of the dimensionless viscosity Oh. Purple squares stand for the VCVPF model and green empty circles for the VPF model. The horizontal dotted black line is $a/b = 3.86$ and gives the elongation of the droplet in the case of the Stokes flow approximation [7, 8].</figcaption></figure>

When looking in more detail, we find that with in-
creasing Oh and thus with increasing viscosity of the liquid, the droplet adopts a more and more elongated shape before pointed tips appear at the ends, typically when $n \sim 1.5$. Indeed, the viscous forces at the interface are stronger for fast moving/deforming boundaries, or alternatively, higher surface modes are more damped than lower modes. As a result, for larger Oh, the tip formation is retarded during the droplet elongation and consequently, tips appear for shapes with higher aspect ratios. For $Oh=10^{-3}$, the maximal elongation is characterized by the ratio $a/b = 2.66$, close to the value of 2.65 found for inviscid fluids [4]. At the other extremity, for $Oh=0.4$, the aspect ratio $a/b$ is 3.72, close to the factor 3.86 found by [7, 8] in the case of Stokes flow (creeping flow with vanishing inertia forces) corresponding to large $Oh \gg 1$. The trend is shown in Fig. 5, where we give the aspect ratio $a/b$ of the droplet at the moment $n = 1.5$ as a function of $Oh \in [0, 1.4]$. For comparison, we also add the results obtained for the VPF model. In the low viscosity domain $Oh< 0.01$, when the dimensionless viscous pressures $\tilde{D}_{nn}$ and $\tilde{p}_{\mu}$ are negligible compared to the difference between the capillary and the Coulomb pressure, VPF and VCPPF give similar pathways and final aspect ratios. However because VCVPF recovers in that domain the correct dispersion relationship, as discussed in a previous section, VCVPF is preferable over VPF also for low viscous liquids. For larger dimensionless viscosities $Oh> 0.8$, the aspect ratio in VCVPF tends to the one obtained in the Stokes flow approximation, while VPF clearly overestimates the elongation of high viscous drops confirming that VCVPF may be better suited for a larger range of viscosities. For water droplets, characterized by $Oh=0.023$, the differences in the pathways between VPF and VCVPF are indeed small but nevertheless notable. Typically, a VCVPF calculation with $Oh=0.023$ is similar to a VPF calculation with $Oh=0.046$.

For the simulated pathways in Fig. 4, once the pointed tips have formed ($n = 1.5$), 20% of the initial charge is suddenly removed and the fissility parameter drops instantaneously to $X = 0.64$. At the same time, we reduced the kinetic energy by 80%, by simply multiplying $\psi$ by the factor $\sqrt{0.2}$. The latter is an arbitrary choice and motivated here so as to allow to simulate also very low viscous drops $Oh = 10^{-3}$, for which the low energy dissipation could otherwise prevent the numerically stable computation of the return path of such low-viscous drops. After charge emission, the pressure due to the surface tension dominates the Coulomb pressure at the surface and the pointed tips retract so as to flatten the ends. This is monitored by the increase of the parameter $n$ from 1.5 to up to 3.0 to fit an elongated capsule. The maximal value of $n$ during the return path depends on the Oh number. Because of viscous forces at the interface, the accumulated kinetic energy is dissipated and the drops return finally to the spherical shape, unlike the simulations done with the inviscid liquid [4, 12].

The shape deformation of viscous drops shows a similar behaviour than the damped forced oscillator. For

$Oh< 0.2$, the dynamics is under-damped and the return to the sphere goes through oblate shapes as indicated by the ratio $a/b < 1$. In the $(a/b, n)$ graph, the droplet oscillates with decreasing amplitude, between the prolate and oblate shapes and between pointed (diamonds) and flattened (capsule) tip shapes. We note that $Oh \simeq 0.2$ corresponds to the value where the droplet is critically damped, meaning that it no longer overshoots when returning to the asymptotic spherical shape. For $Oh> 0.2$, the motion is overdamped and the droplet returns slowly to the stable spherical shape. In the particular case of $Oh \ge 0.4$, the return path is even characterized by $n \le 2$, meaning that the tips stay pointed for a while even after charge emission. This is a surprising result and the prediction may not be meaningful as VCVPF calculations with $Oh \ge 0.4$ are expected to deviate significantly from full NS calculations, as was shown in [22] in their study of small drop surface oscillations.

## B. Influence of the charge carrier mobility

<figure><img src="image_10.png" /><figcaption>FIG. 6: Deformation pathway of a droplet with $Oh=0.023$ (typically a $\emptyset = 50\mu m$ water droplet at room temperature) for various dimensionless charge carrier mobility $\tilde{\lambda}$, ranging from 0.1 to 10 and for the perfect conductor (p.c.).</figcaption></figure>

In the second step, we focus on the influence of the charge mobility of excess charge carrier on the deformation pathway. We consider the particular case where the dimensionless bulk conductivity $\tilde{\kappa} \ll 1$ is sufficiently small so that the charge dynamics are dominated by the surface current $j_s$ and depends mainly on the dimensionless charge mobility $\tilde{\lambda}$. The simulations were performed with $\tilde{\kappa} = 0.1$ and with a dimensionless viscosity of $Oh=0.023$ (corresponding to a $\emptyset = 52\mu m$ droplet of deionized water at room temperature, see table I and II). In figure 6, we show that the first part of the pathway (corresponding to $X=1$) depends notably on the charge mobility $\tilde{\lambda}$. Considering that the charge mobility affects mainly the "fine-structure" of the Coulomb pressure, its influence on the pathway is rather remarkable. We see that with decreasing mobility, the aspect ratio $a/b$ of the droplet at tip formation ($n \sim 1.5$) decreases from 2.85
creasing Oh and thus with increasing viscosity of the liquid, the droplet adopts a more and more elongated shape before pointed tips appear at the ends, typically when $n \sim 1.5$. Indeed, the viscous forces at the interface are stronger for fast moving/deforming boundaries, or alternatively, higher surface modes are more damped than lower modes. As a result, for larger Oh, the tip formation is retarded during the droplet elongation and consequently, tips appear for shapes with higher aspect ratios. For Oh=10⁻³, the maximal elongation is characterized by the ratio $a/b = 2.66$, close to the value of 2.65 found for inviscid fluids [4]. At the other extremity, for Oh=0.4, the aspect ratio $a/b$ is 3.72, close to the factor 3.86 found by [7, 8] in the case of Stokes flow (creeping flow with vanishing inertia forces) corresponding to large Oh $\gg 1$. The trend is shown in Fig. 5, where we give the aspect ratio $a/b$ of the droplet at the moment $n = 1.5$ as a function of Oh $\in [0, 1.4]$. For comparison, we also add the results obtained for the VPF model. In the low viscosity domain Oh< 0.01, when the dimensionless viscous pressures $\tilde{D}_{nn}$ and $\tilde{p}_\mu$ are negligible compared to the difference between the capillary and the Coulomb pressure, VPF and VCPPF give similar pathways and final aspect ratios. However because VCVPF recovers in that domain the correct dispersion relationship, as discussed in a previous section, VCVPF is preferable over VPF also for low viscous liquids. For larger dimensionless viscosities Oh> 0.8, the aspect ratio in VCVPF tends to the one obtained in the Stokes flow approximation, while VPF clearly overestimates the elongation of high viscous drops confirming that VCVPF may be better suited for a larger range of viscosities. For water droplets, characterized by Oh=0.023, the differences in the pathways between VPF and VCVPF are indeed small but nevertheless notable. Typically, a VCVPF calculation with Oh=0.023 is similar to a VPF calculation with Oh=0.046.

For the simulated pathways in Fig. 4, once the pointed tips have formed ($n = 1.5$), 20% of the initial charge is suddenly removed and the fissility parameter drops instantaneously to $X = 0.64$. At the same time, we reduced the kinetic energy by 80%, by simply multiplying $\psi$ by the factor $\sqrt{0.2}$. The latter is an arbitrary choice and motivated here so as to allow to simulate also very low viscous drops Oh = 10⁻³, for which the low energy dissipation could otherwise prevent the numerically stable computation of the return path of such low-viscous drops. After charge emission, the pressure due to the surface tension dominates the Coulomb pressure at the surface and the pointed tips retract so as to flatten the ends. This is monitored by the increase of the parameter $n$ from 1.5 to up to 3.0 to fit an elongated capsule. The maximal value of $n$ during the return path depends on the Oh number. Because of viscous forces at the interface, the accumulated kinetic energy is dissipated and the drops return finally to the spherical shape, unlike the simulations done with the inviscid liquid [4, 12].

The shape deformation of viscous drops shows a similar behaviour than the damped forced oscillator. For

Oh< 0.2, the dynamics is under-damped and the return to the sphere goes through oblate shapes as indicated by the ratio $a/b < 1$. In the $(a/b, n)$ graph, the droplet oscillates with decreasing amplitude, between the prolate and oblate shapes and between pointed (diamonds) and flattened (capsule) tip shapes. We note that Oh $\simeq 0.2$ corresponds to the value where the droplet is critically damped, meaning that it no longer overshoots when returning to the asymptotic spherical shape. For Oh> 0.2, the motion is overdamped and the droplet returns slowly to the stable spherical shape. In the particular case of Oh≥ 0.4, the return path is even characterized by $n \le 2$, meaning that the tips stay pointed for a while even after charge emission. This is a surprising result and the prediction may not be meaningful as VCVPF calculations with Oh≥ 0.4 are expected to deviate significantly from full NS calculations, as was shown in [22] in their study of small drop surface oscillations.

## B. Influence of the charge carrier mobility

<figure><img src="image_10.png" /><figcaption>FIG. 6: Deformation pathway of a droplet with Oh=0.023 (typically a $\emptyset = 50\mu m$ water droplet at room temperature) for various dimensionless charge carrier mobility $\tilde{\lambda}$, ranging from 0.1 to 10 and for the perfect conductor (p.c.).</figcaption></figure>

In the second step, we focus on the influence of the charge mobility of excess charge carrier on the deformation pathway. We consider the particular case where the dimensionless bulk conductivity $\tilde{\kappa} \ll 1$ is sufficiently small so that the charge dynamics are dominated by the surface current $j_s$ and depends mainly on the dimensionless charge mobility $\tilde{\lambda}$. The simulations were performed with $\tilde{\kappa} = 0.1$ and with a dimensionless viscosity of Oh=0.023 (corresponding to a $\emptyset = 52\mu m$ droplet of deionized water at room temperature, see table I and II). In figure 6, we show that the first part of the pathway (corresponding to $X=1$) depends notably on the charge mobility $\tilde{\lambda}$. Considering that the charge mobility affects mainly the "fine-structure" of the Coulomb pressure, its influence on the pathway is rather remarkable. We see that with decreasing mobility, the aspect ratio $a/b$ of the droplet at tip formation ($n \sim 1.5$) decreases from 2.85
to 2.6. Taking $\tilde{\lambda} \gg 10$ gives a pathway which is indistinguishable from the perfect conductor (p.c) case. At $n = 1.5$ the charge is reduced by 20%. The pathway is now dominated by the capillary pressure because the Coulomb pressure has been reduced significantly due to the emitted charge ($X = 0.64$). The effect of the charge mobility on the sub-critical part of the pathways was indeed found negligible and the return pathways are all similar. They differ only because they start with different initial conditions (different $a/b$ ratios).

The reason why smaller charge mobilities $\tilde{\lambda}$ tend to accelerate the tip formation is explained by the competition between the charges being dragged by the fluid and the charges being pushed by the tangent component of the electric field at the surface. For low-viscous drops, $0 < Oh < 0.03$, the velocity at the surface may be so high that the flow drags a large amount of charges to the tips. The large concentration of charges at the tips tends to accelerate the formation of pointed ends. If the accumulated charge density at the tips exceeds the equilibrium charge density $\sigma_{eq}$ (the equilibrium charge corresponds to the charge distribution that results in an equipotential at the surface) it generates an electric field directed from the tip to the equator of the droplet, acting against the drag of the flow. The larger the charge mobility, the larger the surface charge relaxation rate $\tau_s^{-1}$ and the closer the charge distribution stays to the equilibrium charge distribution $\sigma_{eq}$. That is why the deformation pathway tends to the p.c. case for $\tilde{\lambda} \gg 1$. On contrary, for $\tilde{\lambda} \ll 1$ the surface charges move mainly with the fluid to the ends, where they accumulate. The finite charge mobility has thus a non-negligible effect on the shape deformation and should be taken into account when looking at the pathway in detail. We will see below that the charge mobility allows to discriminate between positively and negatively charged pure water droplets.

### C. Influence of the charge polarity: Case of positively and negatively water droplets

By changing the sign of the polarization ring of the injector, negatively or positively charged droplets can be injected and trapped. From the recorded snapshots, we reconstructed in Fig. 7 the deformation pathway of positively (green squares) and negatively (black circles) charged drops. Surprisingly, the pathways differ and the positively charged droplets show at charge breakup lemon like shape with an aspect ratio of $a/b = 2.9$ larger than the negatively charged ones ($a/b = 2.7$). As the radius was 26 µm for both cases and as the viscosity and density are not expected to depend on the sign of the charge, we supposed that the different behavior comes from the different mobilities of the charge carriers in water. Indeed, both protons $H^+$ and hydroxides $OH^-$ in aqueous solution have anomalously large diffusion coefficients, but with the mobility of $H^+$ almost twice as large as that of $OH^-$ anions [29, 31].

<figure><img src="image_11.png" /><figcaption>FIG. 7: Deformation pathway for positively and negatively charged water droplets having a radius of 26 µm. Solid and dashed lines are simulations, black closed circles and green open squares are experimental results of negatively and positively charged deionized water droplets, respectively. In the legend, $\lambda$ stands for the dimensionless charge mobility (tilde has been omitted) and $\Delta Q/Q$ is the relative amount of emitted charge. The red dashed line is added to show the effect on the return path for $\tilde{\lambda} = 2$, in the case where 20% of the charge was removed instead of 40%.</figcaption></figure>

<table><caption>TABLE II: Electrical properties of pure water depending on the polarity of the carried charge. Proton and hydroxide mobilities at room tempertaure are found in [29, 31].</caption><thead><tr><th>Polarity</th><th>κ (μS/m)</th><th>λ (m²/s/V)</th><th>ε<sub>r</sub></th><th>κ̃</th><th>λ̃</th></tr></thead><tbody><tr><td>pos</td><td>5.5</td><td>37 × 10<sup>-8</sup></td><td>80</td><td>0.1</td><td>4.0</td></tr><tr><td>neg</td><td>5.5</td><td>20 × 10<sup>-8</sup></td><td>80</td><td>0.1</td><td>2.1</td></tr></tbody></table>

We simulated the deformation pathway using the VCVPF approach, where the charge dynamics at the surface are explicitly taken into account by (17). The latter depends on the mobility of the charge carriers via the dimensionless mobility $\tilde{\lambda}$ and on the bulk conductivity via the dimensionless factor $\tilde{\kappa}$. In table I and II, we give the relevant parameters of a deionized water droplet used in the simulations. Positively and negatively charged droplets distinguish only by the charge carrier mobility and, in both cases, the dimensionless charge mobility dominates the dimensionless bulk conductivity $\tilde{\lambda} \gg \tilde{\kappa}$. We are thus in the same configuration presented in section III B and the deformation pathway is expected to depend on $\tilde{\lambda}$, with the bulk conductivity having merely a negligible effect on the deformation dynamics. The simulated path is superimposed in Fig. 7 to the experimental data. We found a remarkable agreement between the simulated and observed pathway. We focus our analysis on the first part ($X = 1$) of the pathway for which the charge mobility has a non-negligible influence. We clearly see that positively (green open squares) and neg-
to 2.6. Taking $\tilde{\lambda} \gg 10$ gives a pathway which is indistinguishable from the perfect conductor (p.c) case. At $n = 1.5$ the charge is reduced by 20%. The pathway is now dominated by the capillary pressure because the Coulomb pressure has been reduced significantly due to the emitted charge ($X = 0.64$). The effect of the charge mobility on the sub-critical part of the pathways was indeed found negligible and the return pathways are all similar. They differ only because they start with different initial conditions (different $a/b$ ratios).

The reason why smaller charge mobilities $\tilde{\lambda}$ tend to accelerate the tip formation is explained by the competition between the charges being dragged by the fluid and the charges being pushed by the tangent component of the electric field at the surface. For low-viscous drops, $0 < Oh < 0.03$, the velocity at the surface may be so high that the flow drags a large amount of charges to the tips. The large concentration of charges at the tips tends to accelerate the formation of pointed ends. If the accumulated charge density at the tips exceeds the equilibrium charge density $\sigma_{eq}$ (the equilibrium charge corresponds to the charge distribution that results in an equipotential at the surface) it generates an electric field directed from the tip to the equator of the droplet, acting against the drag of the flow. The larger the charge mobility, the larger the surface charge relaxation rate $\tau_s^{-1}$ and the closer the charge distribution stays to the equilibrium charge distribution $\sigma_{eq}$. That is why the deformation pathway tends to the p.c. case for $\tilde{\lambda} \gg 1$. On contrary, for $\tilde{\lambda} \ll 1$ the surface charges move mainly with the fluid to the ends, where they accumulate. The finite charge mobility has thus a non-negligible effect on the shape deformation and should be taken into account when looking at the pathway in detail. We will see below that the charge mobility allows to discriminate between positively and negatively charged pure water droplets.

### C. Influence of the charge polarity: Case of positively and negatively water droplets

By changing the sign of the polarization ring of the injector, negatively or positively charged droplets can be injected and trapped. From the recorded snapshots, we reconstructed in Fig. 7 the deformation pathway of positively (green squares) and negatively (black circles) charged drops. Surprisingly, the pathways differ and the positively charged droplets show at charge breakup lemon like shape with an aspect ratio of $a/b = 2.9$ larger than the negatively charged ones ($a/b = 2.7$). As the radius was 26 µm for both cases and as the viscosity and density are not expected to depend on the sign of the charge, we supposed that the different behavior comes from the different mobilities of the charge carriers in water. Indeed, both protons $H^+$ and hydroxides $OH^-$ in aqueous solution have anomalously large diffusion coefficients, but with the mobility of $H^+$ almost twice as large as that of $OH^-$ anions [29, 31].

<figure><img src="image_11.png" /><figcaption>FIG. 7: Deformation pathway for positively and negatively charged water droplets having a radius of 26 µm. Solid and dashed lines are simulations, black closed circles and green open squares are experimental results of negatively and positively charged deionized water droplets, respectively. In the legend, $\lambda$ stands for the dimensionless charge mobility (tilde has been omitted) and $\Delta Q/Q$ is the relative amount of emitted charge. The red dashed line is added to show the effect on the return path for $\tilde{\lambda} = 2$, in the case where 20% of the charge was removed instead of 40%.</figcaption></figure>

<table><caption>TABLE II: Electrical properties of pure water depending on the polarity of the carried charge. Proton and hydroxide mobilities at room tempertaure are found in [29, 31].</caption><thead><tr><th>Polarity</th><th>κ (µS/m)</th><th>λ (m²/s/V)</th><th>ε<sub>r</sub></th><th>κ̃</th><th>λ̃</th></tr></thead><tbody><tr><td>pos</td><td>5.5</td><td>37 × 10<sup>-8</sup></td><td>80</td><td>0.1</td><td>4.0</td></tr><tr><td>neg</td><td>5.5</td><td>20 × 10<sup>-8</sup></td><td>80</td><td>0.1</td><td>2.1</td></tr></tbody></table>

We simulated the deformation pathway using the VCVPF approach, where the charge dynamics at the surface are explicitly taken into account by (17). The latter depends on the mobility of the charge carriers via the dimensionless mobility $\tilde{\lambda}$ and on the bulk conductivity via the dimensionless factor $\tilde{\kappa}$. In table I and II, we give the relevant parameters of a deionized water droplet used in the simulations. Positively and negatively charged droplets distinguish only by the charge carrier mobility and, in both cases, the dimensionless charge mobility dominates the dimensionless bulk conductivity $\tilde{\lambda} \gg \tilde{\kappa}$. We are thus in the same configuration presented in section III B and the deformation pathway is expected to depend on $\tilde{\lambda}$, with the bulk conductivity having merely a negligible effect on the deformation dynamics. The simulated path is superimposed in Fig. 7 to the experimental data. We found a remarkable agreement between the simulated and observed pathway. We focus our analysis on the first part ($X = 1$) of the pathway for which the charge mobility has a non-negligible influence. We clearly see that positively (green open squares) and neg-
atively (black closed circles) charged water droplets can already be distinguished by the first part of their respective pathway. The VCVPF calculation with $\tilde{\lambda} = 4$ follows rather the path of the positively charged droplets (green squares) with a final $a/b$ ratio of about 2.8 for $n = 1.5$. On the other hand, the simulation with $\tilde{\lambda} = 2$ agrees better with the measurements of negatively (black points) charged droplets. The comparison corroborates that the observed difference in the ($X = 1$) part of the pathway can be explained by the difference in the charge mobility of $H^+$ and $OH^-$ in water.

After charge breakup, the subcritical pathway is dominated by the amount of emitted charge. In particular, the slope of the path after charge emission in the graph ($a/b, n$) increases with increasing emitted charge, as the larger difference between the capillary pressure and Coulomb pressure accelerates the flattening of the tips. In the case of negatively charged droplets, we found that removing 40% of the initial charge gives the best "visual fit" for the return pathway. Removing only 20% of the charge (see red dashed line in Fig. 7), shifts the return path by $\simeq 0.2$ to the right on the $a/b$ axis, resulting in a less satisfying visual fit with the black closed circles. In the case of positively charged droplets the best visual fit with the green squares was obtained for a relative emitted charge $\Delta Q/Q$ of 20%. These findings are corroborated by the observations reported by R. Müller [6], who found that negatively charged pure water droplets at room temperature having $\emptyset = 62$ µm emit about 40% of the initial charge at charge breakup while positively emit merely 20%. We note however that the large scattering of the green squares and black dots in Fig. 7 does not allow us to clearly discriminate between 20% and 40% of the emitted charge.

It seems thus that the reason to why positively charged pure water droplets have a aspect ratio $a/b$ larger than negatively charged ones is because $H^+$ has a higher charge mobility than $OH^-$. For such low-viscous droplets, Oh= 0.023, the drag of the charges by the fluid velocity is in competition with the push by the electric field. It happens that for $OH^-$ the latter effect is less efficient so that more charges accumulate at the tips, accelerating the tip formation. The charge velocity tangent to the surface, given as the sum of the drag velocity $\tilde{u}_\tau$ and "electric field" velocity $\tilde{\lambda}\tilde{E}_\tau$ is illustrated in the upper panel of Fig.8 for three different values of $\tilde{\lambda}$. The quantities in Fig.8 are given as a function of the dimensionless distance $r/R$, corresponding to the distance from the symmetry axis to a surface point, such that $r \to 0$ corresponds to the tips (see inset). We note that close to the tips ($r/R \simeq 0.06$), the velocity of the charges is reduced by a factor 2 when increasing $\tilde{\lambda} = 1$ to $\tilde{\lambda} = 4$. The accumulated charge density at the tips (lower panel) for $\tilde{\lambda} = 1$ was also found to be larger about by a factor 2 than for $\tilde{\lambda} = 4$, accelerating thus the tip formation. The higher charge density may also explain why negatively charged pure water droplets emit more charge than positively charged pure water droplets, as was observed in [6].

<figure><img src="image_12.png" /><figcaption>FIG. 8: (a) Dimensionless velocity of the charges tangent to the surface, $\tilde{u}_\tau + \tilde{\lambda}\tilde{E}_\tau$, as a function of the dimensionless distance $r/R$ defined in the inset. Positive values indicate velocities directed towards the tips. (b) Dimensionless surface charge density $\tilde{\sigma}$ as a function of $r/R$. The quantities are calculated for three different charge mobilities $\tilde{\lambda}$. The blue line with pluses gives the equilibrium surface charge density $\tilde{\sigma}_{eq}$ in the case of a perfectly conducting liquid (p.c.).</figcaption></figure>

If instead of pure water we had used sea water, having a bulk conductivity of 5 S/m so that $\tilde{\kappa} \simeq 10^4 \gg \tilde{\lambda}_\alpha$ for all excess charge carrier $\alpha$, the charge dynamics would be dominated by the bulk current $\tilde{\kappa}\tilde{E}_n$ and the liquid well approached by a perfect conductor. We expect thus positively and negatively charge droplets of sea water to yield the same deformation path and emit the same amount of charge at charge break-up. The latter prediction may be verified in future experiments.

## IV. CONCLUSION

We conclude this work by highlighting the remarkable results. On the experimental front, we trapped charged pure water droplets in a Paul trap and took snapshots of the ultra-fast shape deformation after the onset of the Coulomb instability. This was done for positively and negatively charged droplets at room temperature. We noted that, just before charge emission, negatively charged droplets are more elongated than positively ones. On the numerical front, we developed a numerical code, based on the VCVPF approach, which includes viscous effects at the interface. The fluid dynamics is coupled to the surface charge dynamics, so that the Coulomb pressure depends also on the electrical properties of the liquid. For the particular case of pure water droplets,
atively (black closed circles) charged water droplets can already be distinguished by the first part of their respective pathway. The VCVPF calculation with $\tilde{\lambda} = 4$ follows rather the path of the positively charged droplets (green squares) with a final $a/b$ ratio of about 2.8 for $n = 1.5$. On the other hand, the simulation with $\tilde{\lambda} = 2$ agrees better with the measurements of negatively (black points) charged droplets. The comparison corroborates that the observed difference in the ($X = 1$) part of the pathway can be explained by the difference in the charge mobility of $H^+$ and $OH^-$ in water.

After charge breakup, the subcritical pathway is dominated by the amount of emitted charge. In particular, the slope of the path after charge emission in the graph ($a/b, n$) increases with increasing emitted charge, as the larger difference between the capillary pressure and Coulomb pressure accelerates the flattening of the tips. In the case of negatively charged droplets, we found that removing 40% of the initial charge gives the best "visual fit" for the return pathway. Removing only 20% of the charge (see red dashed line in Fig. 7), shifts the return path by $\simeq 0.2$ to the right on the $a/b$ axis, resulting in a less satisfying visual fit with the black closed circles. In the case of positively charged droplets the best visual fit with the green squares was obtained for a relative emitted charge $\Delta Q/Q$ of 20%. These findings are corroborated by the observations reported by R. Müller [6], who found that negatively charged pure water droplets at room temperature having $\emptyset = 62$ µm emit about 40% of the initial charge at charge breakup while positively emit merely 20%. We note however that the large scattering of the green squares and black dots in Fig. 7 does not allow us to clearly discriminate between 20% and 40% of the emitted charge.

It seems thus that the reason to why positively charged pure water droplets have a aspect ratio $a/b$ larger than negatively charged ones is because $H^+$ has a higher charge mobility than $OH^-$. For such low-viscous droplets, $\text{Oh}=0.023$, the drag of the charges by the fluid velocity is in competition with the push by the electric field. It happens that for $OH^-$ the latter effect is less efficient so that more charges accumulate at the tips, accelerating the tip formation. The charge velocity tangent to the surface, given as the sum of the drag velocity $\tilde{u}_\tau$ and "electric field" velocity $\tilde{\lambda}\tilde{E}_\tau$ is illustrated in the upper panel of Fig.8 for three different values of $\tilde{\lambda}$. The quantities in Fig.8 are given as a function of the dimensionless distance $r/R$, corresponding to the distance from the symmetry axis to a surface point, such that $r \to 0$ corresponds to the tips (see inset). We note that close to the tips ($r/R \simeq 0.06$), the velocity of the charges is reduced by a factor 2 when increasing $\tilde{\lambda} = 1$ to $\tilde{\lambda} = 4$. The accumulated charge density at the tips (lower panel) for $\tilde{\lambda} = 1$ was also found to be larger about by a factor 2 than for $\tilde{\lambda} = 4$, accelerating thus the tip formation. The higher charge density may also explain why negatively charged pure water droplets emit more charge than positively charged pure water droplets, as was observed in [6].

<figure><img src="image_12.png" /><figcaption>FIG. 8: (a) Dimensionless velocity of the charges tangent to the surface, $\tilde{u}_\tau + \tilde{\lambda}\tilde{E}_\tau$, as a function of the dimensionless distance $r/R$ defined in the inset. Positive values indicate velocities directed towards the tips. (b) Dimensionless surface charge density $\tilde{\sigma}$ as a function of $r/R$. The quantities are calculated for three different charge mobilities $\tilde{\lambda}$. The blue line with pluses gives the equilibrium surface charge density $\tilde{\sigma}_{eq}$ in the case of a perfectly conducting liquid (p.c.).</figcaption></figure>

If instead of pure water we had used sea water, having a bulk conductivity of 5 S/m so that $\tilde{\kappa} \simeq 10^4 \gg \tilde{\lambda}_\alpha$ for all excess charge carrier $\alpha$, the charge dynamics would be dominated by the bulk current $\tilde{\kappa}\tilde{E}_n$ and the liquid well approached by a perfect conductor. We expect thus positively and negatively charge droplets of sea water to yield the same deformation path and emit the same amount of charge at charge break-up. The latter prediction may be verified in future experiments.

## IV. CONCLUSION

We conclude this work by highlighting the remarkable results. On the experimental front, we trapped charged pure water droplets in a Paul trap and took snapshots of the ultra-fast shape deformation after the onset of the Coulomb instability. This was done for positively and negatively charged droplets at room temperature. We noted that, just before charge emission, negatively charged droplets are more elongated than positively ones. On the numerical front, we developed a numerical code, based on the VCVPF approach, which includes viscous effects at the interface. The fluid dynamics is coupled to the surface charge dynamics, so that the Coulomb pressure depends also on the electrical properties of the liquid. For the particular case of pure water droplets,
we showed that the VCVPF model describes accurately all steps of the shape deformation. The pathway from the initial spherical shape to the lemon shape, the instantaneous charge emission with negligible kinetic energy loss and the pathway back to the sphere by going through exactly the same exotic shapes are well reproduced. The excellent agreement gives much confidence in the underlying model and numerical algorithms of the code. We found that the VCVPF model predicts that with increasing viscosity, or more precisely, with increasing Oh, the final shape before charge break-up is more and more elongated. We found that, at charge emission, the ratio a/b lies between 2.66 for inviscid liquids and 3.86 in the limiting case of a Stokes flows. We found that, if the dimensionless charge mobility $\tilde{\lambda}$ dominates the dimensionless bulk conductivity $\tilde{\kappa}$ of the liquid while being of the order of 1, the droplet deformation is sensitive to the mobility of the charge carriers. We showed that for $\tilde{\lambda} = 1$ the charge dynamics is dominated by the fluid velocity at the interface which drags the charge carriers to the tips, accelerating the formation of pointed tips. On the other hand, for $\tilde{\lambda} = 4$ the electric field is efficient enough to avoid a rapid accumulation of charge at the tips, delaying the formation of pointed tips. Finally, our simulations explained that, even though the charge mobility of protons is higher only by a factor 2 than the one of OH$^-$, positively and negatively charged droplets show a notable different aspect ratio at charge breakup. The lower mobility of OH$^-$ may also explain why more charge (by factor 2) is released for negatively charged droplets at charge breakup.

## Acknowledgments

This work was funded by the French National Research Agency (ANR), program ANR-07-JCJC-0062-01 (DynaMIC).

# SUPPLEMENTARY MATERIAL

## 1. Surface divergence

Let $\vec{\tau}$ be a vector tangent to the surface defined by the vector $\vec{s}$. The surface divergence operator

$$ \vec{\nabla}_s \cdot (\Psi(\vec{s}) \vec{\tau}(\vec{s})) \quad (35) $$

for arbitrary function $\Psi(\vec{s})$ defined on the surface, may not be well known to the reader, so we give a computationally useful expression of the surface divergence found in (17). Let ($\xi > 1, -1 \le \eta \le 1, 0 \le \varphi \le 2\pi$) be curvilinear orthogonal coordinates and $\vec{r}(\xi, \eta, \varphi)$ an arbitrary vector. Because the droplet deformation and charge dynamics are axisymmetric, the function $\Psi(\xi, \eta) = \sigma(u_\tau + \lambda E_\tau)$ does not depend on the azimuthal angle $\varphi$. The surface vector $\vec{s}(\eta)$ depends thus only on $\eta$. The general expression simplifies in that case to

$$ \begin{aligned} \vec{\nabla}_s \cdot (\vec{\tau}(\vec{s}) \Psi(\vec{s})) &= \vec{n} \cdot (\vec{\nabla} \wedge (\Psi(\vec{r})\vec{e}_\varphi)) \Big|_{\vec{r}=\vec{s}} \\ &= \vec{n} \cdot \vec{e}_\xi \frac{1}{h_\eta h_\varphi} \frac{d}{d\eta} (h_\varphi(\vec{s})\Psi(\vec{s})) \quad (36) \end{aligned} $$

where $\vec{\tau} = \vec{e}_\varphi \wedge \vec{n}$ is the tangent vector and $h_\eta$ and $h_\varphi$ the scale factors associated to the curvilinear coordinate system and evaluated at the surface. The term $d(\cdot)/d\eta$ in (36) can be accurately computed using for example a Gauss-Lobatto interpolation scheme, which is appropriate if the endpoints are known (here $h_\varphi(\vec{s})\Psi(\vec{s}) = 0$ for $\eta = \pm 1$).

## 2. Numerical method for computing the viscous potential

The VCVPF model is only useful if one can compute the viscous pressure $p_\mu$, which enters the normal stress term. We propose here a numerical method to compute the viscous pressure in VCVPF. Making use of the fact that the Laplace operator is separable in prolate spheroid coordinate system ($\xi \ge 1, -1 \le \eta \le 1, 0 \le \varphi \le 2\pi$), we get a simple expression for the harmonic functions $f_m$ that satisfy the Laplace equation $\vec{\nabla}^2 f_m = 0$. As we want $f_m(\xi, \eta) = f_m(\xi, -\eta)$ to be axisymmetric and even and defined at the origin $f(1, 0) < \infty$, we have

$$ f_m(\xi, \eta) = P_{2m}(\xi)P_{2m}(\eta) \quad (37) $$

where $P_m(\cdot)$ are Legendre polynomials of order $m$. Let $\{f_m\}$ be a basis of harmonic function. We may now expand the viscous pressure on the harmonic basis $f_m$,

$$ p_\mu(\xi, \eta) = \sum_{m=1}^{M} b_m f_m(\xi, \eta) \quad (38) $$

The expansion coefficients $\{b_m\}$ are deduced using the boundary condition (24),

$$ \sum_{m=1}^{M} b_m \underbrace{f_m(\vec{s}_i) u_n(\vec{s}_i) ds_i}_{F_{m,i}} = \underbrace{u_\tau(\vec{s}_i)(2\mu D_{\tau n}(\vec{s}_i) ds_i)}_{d_i} - \underbrace{u_\tau(\vec{s}_i) \sigma_i(\vec{s}_i) E_\tau(\vec{s}_i)) ds_i}_{e_i} \quad (39) $$

where $\vec{s}_i \le N$ is an ensemble of $N$ points at the surface. Equation (39) is a system of $N$ equations and $M$ unknowns of the form

$$ \sum_{m=1}^{M} b_m F_{m,i} = d_i - e_i, \quad 1 \le i \le N \quad (40) $$

If $N \ge M$ the system is overdetermined and can be solved using a singular value decomposition (SVD), which finds the coefficient $\{b_m\}$ as least squares solutions.
we showed that the VCVPF model describes accurately all steps of the shape deformation. The pathway from the initial spherical shape to the lemon shape, the instantaneous charge emission with negligible kinetic energy loss and the pathway back to the sphere by going through exactly the same exotic shapes are well reproduced. The excellent agreement gives much confidence in the underlying model and numerical algorithms of the code. We found that the VCVPF model predicts that with increasing viscosity, or more precisely, with increasing Oh, the final shape before charge break-up is more and more elongated. We found that, at charge emission, the ratio a/b lies between 2.66 for inviscid liquids and 3.86 in the limiting case of a Stokes flows. We found that, if the dimensionless charge mobility $\tilde{\lambda}$ dominates the dimensionless bulk conductivity $\tilde{\kappa}$ of the liquid while being of the order of 1, the droplet deformation is sensitive to the mobility of the charge carriers. We showed that for $\tilde{\lambda} = 1$ the charge dynamics is dominated by the fluid velocity at the interface which drags the charge carriers to the tips, accelerating the formation of pointed tips. On the other hand, for $\tilde{\lambda} = 4$ the electric field is efficient enough to avoid a rapid accumulation of charge at the tips, delaying the formation of pointed tips. Finally, our simulations explained that, even though the charge mobility of protons is higher only by a factor 2 than the one of OH$^-$, positively and negatively charged droplets show a notable different aspect ratio at charge breakup. The lower mobility of OH$^-$ may also explain why more charge (by factor 2) is released for negatively charged droplets at charge breakup.

## Acknowledgments

This work was funded by the French National Research Agency (ANR), program ANR-07-JCJC-0062-01 (DynaMIC).

# SUPPLEMENTARY MATERIAL

## 1. Surface divergence

Let $\vec{\tau}$ be a vector tangent to the surface defined by the vector $\vec{s}$. The surface divergence operator

$$ \vec{\nabla}_s \cdot (\Psi(\vec{s}) \vec{\tau}(\vec{s})) \quad (35) $$

for arbitrary function $\Psi(\vec{s})$ defined on the surface, may not be well known to the reader, so we give a computationally useful expression of the surface divergence found in (17). Let ($\xi > 1, -1 \le \eta \le 1, 0 \le \varphi \le 2\pi$) be curvilinear orthogonal coordinates and $\vec{r}(\xi, \eta, \varphi)$ an arbitrary vector. Because the droplet deformation and charge dynamics are axisymmetric, the function $\Psi(\xi, \eta) = \sigma(u_\tau + \lambda E_\tau)$ does not depend on the azimuthal angle $\varphi$. The surface vector $\vec{s}(\eta)$ depends thus only on $\eta$. The general expression simplifies in that case to

$$ \begin{aligned} \vec{\nabla}_s \cdot (\vec{\tau}(\vec{s}) \Psi(\vec{s})) &= \vec{n} \cdot (\vec{\nabla} \wedge (\Psi(\vec{r}) \vec{e}_\varphi)) \Big|_{\vec{r}=\vec{s}} \\ &= \vec{n} \cdot \vec{e}_\xi \frac{1}{h_\eta h_\varphi} \frac{d}{d\eta} (h_\varphi(\vec{s}) \Psi(\vec{s})) \quad (36) \end{aligned} $$

where $\vec{\tau} = \vec{e}_\varphi \wedge \vec{n}$ is the tangent vector and $h_\eta$ and $h_\varphi$ the scale factors associated to the curvilinear coordinate system and evaluated at the surface. The term $d(\cdot)/d\eta$ in (36) can be accurately computed using for example a Gauss-Lobatto interpolation scheme, which is appropriate if the endpoints are known (here $h_\varphi(\vec{s}) \Psi(\vec{s}) = 0$ for $\eta = \pm 1$).

## 2. Numerical method for computing the viscous potential

The VCVPF model is only useful if one can compute the viscous pressure $p_\mu$, which enters the normal stress term. We propose here a numerical method to compute the viscous pressure in VCVPF. Making use of the fact that the Laplace operator is separable in prolate spheroid coordinate system ($\xi \ge 1, -1 \le \eta \le 1, 0 \le \varphi \le 2\pi$), we get a simple expression for the harmonic functions $f_m$ that satisfy the Laplace equation $\vec{\nabla}^2 f_m = 0$. As we want $f_m(\xi, \eta) = f_m(\xi, -\eta)$ to be axisymmetric and even and defined at the origin $f(1, 0) < \infty$, we have

$$ f_m(\xi, \eta) = P_{2m}(\xi) P_{2m}(\eta) \quad (37) $$

where $P_m(\cdot)$ are Legendre polynomials of order $m$. Let $\{f_m\}$ be a basis of harmonic function. We may now expand the viscous pressure on the harmonic basis $f_m$,

$$ p_\mu(\xi, \eta) = \sum_{m=1}^{M} b_m f_m(\xi, \eta) \quad (38) $$

The expansion coefficients $\{b_m\}$ are deduced using the boundary condition (24),

$$ \sum_{m=1}^{M} b_m \underbrace{f_m(\vec{s}_i) u_n(\vec{s}_i) ds_i}_{F_{m,i}} = \underbrace{u_\tau(\vec{s}_i)(2\mu D_{\tau n}(\vec{s}_i) ds_i)}_{d_i} - \underbrace{u_\tau(\vec{s}_i) \sigma_i(\vec{s}_i) E_\tau(\vec{s}_i)) ds_i}_{e_i} \quad (39) $$

where $\vec{s}_i \le N$ is an ensemble of $N$ points at the surface. Equation (39) is a system of $N$ equations and $M$ unknowns of the form

$$ \sum_{m=1}^{M} b_m F_{m,i} = d_i - e_i, \quad 1 \le i \le N \quad (40) $$

If $N \ge M$ the system is overdetermined and can be solved using a singular value decomposition (SVD), which finds the coefficient $\{b_m\}$ as least squares solutions.
[1] "On the Equilibrium of Liquid Conducting Masses Charged with Electricity", Lord Rayleigh, Philosophical Magazine, Vol. 14, No. 87, pp. 184-186 (1882)

[2] D.Duft, H. Lebius, B.A. Huber, C. Guet, T. Leisner, Phys. Rev. Lett. 89 0845031 (2002)

[3] D.Duft, T. Achtzehn, R. Müller, B.A. Huber, and T. Leisner, Nature (London) 421, 128 (2003)

[4] E. Giglio, B. Gervais, J. Rangama, B. Manil, B.A. Huber, D. Duft, R. Müller, T. Leisner, C. Guet Phys Rev E 77 p 036319 (2008)

[5] T. Achtzehn, R. Müller, D. Duft and T. Leisner, Eur. Phys. J. D 34 p 311 (2005)

[6] Dissertation by René. Müller (2010), Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen, Fakultät für Mathematik und Naturwissenschaften der Technischen Universität Ilmenau. https://www.db-thueringen.de/receive/dbt_mods_00015683

[7] S.I. Betelu, M.A. Fontelos, U. Kindelan and O. Vantzos, Phys. Fluids 18 051706 (2006)

[8] Neha Gawande, Y. S. Mayya, and Rochish Thaokar, Phys. Rev. Fluids 2, 113603 (2017)

[9] J. F. de la Mora, Annu. Rev. Fluid Mech. 39, 217 (2007).

[10] A. Radcliffe, International Journal for Numerical Methods in Fluids 71, Issue 2, p 249-268 (2013), https://doi.org/10.1002/ftd.3667

[11] A. Radcliffe, Int J Model Simul., Vol. 7 No. 3 (2016) 1650016

[12] J. C. Burton and P. Taborek, Phys. Rev. Lett. 106 p 144501 (2011)

[13] H. Nganguya, Y.-N. Young, A. T. Layton, M.-C. Lai, W.-F. Hu "Electrohydrodynamics of a viscous drop with inertia" PHYSICAL REVIEW E 93, 053114 (2016)

[14] R.T. Collins, J.J. Jones, M.T. Harris and O.A. Basaran, Nature Physics Vol 4, Issue 2, 149 (2008)

[15] R.T. Collins, K. Sambath, M.T. Harris and O.A. Basaran, Proc. Natl. Acad. Sci. March 26, (2013) 110 (13) 4905-4910; https://doi.org/10.1073/pnas.1213708110
Proc Natl Acad Sci U S A. 2013 Mar 26;110(13):4905-10. doi: 10.1073/pnas.1213708110

[16] A. M. Gañán-Calvo, J. M. López-Herrera, N. Rebollo-Muñoz and J. M. Montanero, Scientific Reports, 6:32357 (2016), https://doi.org/10.1038/srep32357

[17] A.G. Bailey, Atomisation and Spray Technology (ISSN 0266-3481), vol. 2, no. 2, (1986), p. 95-134.

[18] J. R. Melcher and G. I. Taylor, Annu. Rev. Fluid Mech. 1, p 111. (1969)

[19] I.W. McAllister, J. Phys. D 23 359 (1990)

[20] D.A. Saville "The Taylor-Melcher Leaky Dielectric Model", Annual Review of Fluid Mechanics, Vol. 29, p 27-64 (1997) https://doi.org/10.1146/annurev.fluid.29.1.27

[21] A. Castellanos "Electrohydrodynamics", section 1.4, p 13, Springer

[22] J.C.Padrino T.Funada D.D.Joseph, International Journal of Multiphase Flow Volume 34, Issue 1, p 61-75 (2008)

[23] D.D. Joseph and J. Wang, J. Fluid. Mech. vol 505, p 365 (2004)

[24] J. WANG, D.D. JOSEPH and T. FUNADA, Journal of Fluid Mechanics 552 p 383 (2005)

[25] J. Wang and D.D. Joseph, J. Fluid. Mech. vol 557, p 145 (2006)

[26] C. Pantano, A. M. Gañán-Calvo, and A. Barrero, J. Aerosol Sci. 25, 1065 (1994).

[27] A. M. Gañán-Calvo, Phys. Rev. Lett. 79 p 217 (1997)

[28] A. N. Zharov, A. I. Grigor'ev and S. O. Shiryaeva, Tech. Phys. (2005) Vol 50, p 1124, https://doi.org/10.1134/1.2051449

[29] P. W. Atkins and Julio De Paula, Publisher: Oxford ; New York : Oxford University Press, (2006) ISBN 0198700725 (table 21.5 page 763)

[30] E. Weisstein, "Superellipsoid." From MathWorld - A Wolfram Web Resource. http://mathworld.wolfram.com/Superellipsoid.html

[31] S. Hi Lee and J. C. Rasaiah, The Journal of Chemical Physics, 135, 124505 (2011)
[1] "On the Equilibrium of Liquid Conducting Masses Charged with Electricity", Lord Rayleigh, Philosophical Magazine, Vol. 14, No. 87, pp. 184-186 (1882)

[2] D.Duft, H. Lebius, B.A. Huber, C. Guet, T. Leisner, Phys. Rev. Lett. 89 0845031 (2002)

[3] D.Duft, T. Achtzehn, R. Müller, B.A. Huber, and T. Leisner, Nature (London) 421, 128 (2003)

[4] E. Giglio, B. Gervais, J. Rangama, B. Manil, B.A. Huber, D. Duft, R. Müller, T. Leisner, C. Guet Phys Rev E 77 p 036319 (2008)

[5] T. Achtzehn, R. Müller, D. Duft and T. Leisner, Eur. Phys. J. D 34 p 311 (2005)

[6] Dissertation by René. Müller (2010), Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen, Fakultät für Mathematik und Naturwissenschaften der Technischen Universität Ilmenau. https://www.db-thueringen.de/receive/dbt_mods_00015683

[7] S.I. Betelu, M.A. Fontelos, U. Kindelan and O. Vantzos, Phys. Fluids 18 051706 (2006)

[8] Neha Gawande, Y. S. Mayya, and Rochish Thaokar, Phys. Rev. Fluids 2, 113603 (2017)

[9] J. F. de la Mora, Annu. Rev. Fluid Mech. 39, 217 (2007).

[10] A. Radcliffe, International Journal for Numerical Methods in Fluids 71, Issue 2, p 249-268 (2013), https://doi.org/10.1002/ftd.3667

[11] A. Radcliffe, Int J Model Simul., Vol. 7 No. 3 (2016) 1650016

[12] J. C. Burton and P. Taborek, Phys. Rev. Lett. 106 p 144501 (2011)

[13] H. Nganguya, Y.-N. Young, A. T. Layton, M.-C. Lai, W.-F. Hu "Electrohydrodynamics of a viscous drop with inertia" PHYSICAL REVIEW E 93, 053114 (2016)

[14] R.T. Collins, J.J. Jones, M.T. Harris and O.A. Basaran, Nature Physics Vol 4, Issue 2, 149 (2008)

[15] R.T. Collins, K. Sambath, M.T. Harris and O.A. Basaran, Proc. Natl. Acad. Sci. March 26, (2013) 110 (13) 4905-4910; https://doi.org/10.1073/pnas.1213708110
Proc Natl Acad Sci U S A. 2013 Mar 26;110(13):4905-10. doi: 10.1073/pnas.1213708110

[16] A. M. Gañán-Calvo, J. M. López-Herrera, N. Rebollo-Muñoz and J. M. Montanero, Scientific Reports, 6:32357 (2016), https://doi.org/10.1038/srep32357

[17] A.G. Bailey, Atomisation and Spray Technology (ISSN 0266-3481), vol. 2, no. 2, (1986), p. 95-134.

[18] J. R. Melcher and G. I. Taylor, Annu. Rev. Fluid Mech. 1, p 111. (1969)

[19] I.W. McAllister, J. Phys. D 23 359 (1990)

[20] D.A. Saville "The Taylor-Melcher Leaky Dielectric Model", Annual Review of Fluid Mechanics, Vol. 29, p 27-64 (1997) https://doi.org/10.1146/annurev.fluid.29.1.27

[21] A. Castellanos "Electrohydrodynamics", section 1.4, p 13, Springer

[22] J.C.Padrino T.Funada D.D.Joseph, International Journal of Multiphase Flow Volume 34, Issue 1, p 61-75 (2008)

[23] D.D. Joseph and J. Wang, J. Fluid. Mech. vol 505, p 365 (2004)

[24] J. WANG, D.D. JOSEPH and T. FUNADA, Journal of Fluid Mechanics 552 p 383 (2005)

[25] J. Wang and D.D. Joseph, J. Fluid. Mech. vol 557, p 145 (2006)

[26] C. Pantano, A. M. Gañán-Calvo, and A. Barrero, J. Aerosol Sci. 25, 1065 (1994).

[27] A. M. Gañán-Calvo, Phys. Rev. Lett. 79 p 217 (1997)

[28] A. N. Zharov, A. I. Grigor'ev and S. O. Shiryaeva, Tech. Phys. (2005) Vol 50, p 1124, https://doi.org/10.1134/1.2051449

[29] P. W. Atkins and Julio De Paula, Publisher: Oxford ; New York : Oxford University Press, (2006) ISBN 0198700725 (table 21.5 page 763)

[30] E. Weisstein, "Superellipsoid." From MathWorld - A Wolfram Web Resource.
http://mathworld.wolfram.com/Superellipsoid.html

[31] S. Hi Lee and J. C. Rasaiah, The Journal of Chemical Physics, 135, 124505 (2011)