# Charge dynamics at the interface separating two dielectric media.

Eric Giglio

► To cite this version:

Eric Giglio. Charge dynamics at the interface separating two dielectric media.: Two cases: (1) Shape deformation of droplets charged at the Rayleigh limit. (2) Ion transport through insulating capillaries. Fluid mechanics [physics.class-ph]. Normandie Université, 2019. tel-04520514

HAL Id: tel-04520514

https://hal.science/tel-04520514

Submitted on 27 Mar 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>
Pour obtenir le diplôme d'Habilitation à Diriger des Recherches

Spécialité PHYSIQUE

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

Boulevard Henri Becquerel, BP5133, 14070 Caen cedex 5

Charge dynamics at the interface of two dielectric media.

Two cases :

(i) shape deformation of droplets charged at the Rayleigh limit.

(ii) ion transport through insulating capillaries

Présentée et soutenue par

Eric GIGLIO

Thèse soutenue publiquement le 3 décembre 2019
devant le jury composé de

<table><tr><td>Mme. Francesca GULMINELLI</td><td>Professeur de l'Université CAEN-NORMANDIE<br>Lab. de Physique Corpusculaire (LPC)</td><td>Rapporteur</td></tr><tr><td>M. Christoph LEMELL</td><td>Professeur de l'Univ. of Technology of Vienna<br>Institut for Theoretical Physics</td><td>Rapporteur</td></tr><tr><td>M. Bruno MANIL</td><td>Professeur de l'Université Paris 13<br>Laboratoire de Physique des Lasers (LPL)</td><td>Rapporteur</td></tr><tr><td>M. Jose Maria MONTANERO</td><td>Professeur de l'Université d'Extremadura (Espagne)<br>Department of Mechanic Engineering, Energy and Materials,</td><td>Examinateur</td></tr><tr><td>M. Philippe RONCIN</td><td>Directeur de Recherche au CNRS<br>Institut des Sciences Moléculaires d'Orsay, Paris-11, (ISMO)</td><td>Examinateur</td></tr><tr><td>M. Emmanuel BALANZAT</td><td>Directeur de Recherche au CNRS<br>Centre de Rech. sur les Ions, les Matériaux et la Photonique (CIMAP)</td><td>Examinateur, Garant</td></tr></table>
<figure><img src="image_2.png" /></figure>
Charge dynamics at the interface separating two
dielectric media.

(1) Shape deformation of droplets charged at the
Rayleigh limit

(2) Ion transport through insulating capillaries

Giglio Eric

giglio@ganil.fr

https://cv.archives-ouvertes.fr/eric-giglio

https://www.researchgate.net/profile/Eric_Giglio/research

Centre de Recherche sur les Ions, les Matériaux et la Photonique
(CIMAP), Normandie Univ, ENSICAEN, UNICAEN, CEA, CNRS

Boulevard Henri Becquerel BP5133 14070 Caen cedex 5, France
<figure><img src="image_3.png" /></figure>
Short note about of the choice of two topics presented in
the manuscript

In this short note, the author wants two highlight the historical link between the seemingly independent topics presented in chapter one and chapter two. While both topics study different objects and physical phenomena, they have in common the charge dynamics of free charges at the interface separating two dielectric media. From 2007 to 2013, the author studied the spontaneous droplet deformation of critically charged droplets. During that time, the author developed, among other, a model for the surface charge dynamics at the droplet surface, i.e., at the interface separating a charged conducting liquid from a gas interface with negligible viscosity. In 2013, the author was asked for theoretical support in the transport of an ion beam through insulating capillaries. The knowledge gathered from the previous work allowed the author to imagine a (CPU efficient) approach for the charge dynamics in the capillary wall and evaluation of the electric field. The numerical model developed for the computation of the beam trajectories through the capillary in chapter two in thus based on the charge dynamics model derived from the droplet topic presented in chapter one. From 2013 up to now, the numerical model is used successfully to describe and predict experimental results in the field of beam transport through insulating capillaries.
<figure><img src="image_4.png" /></figure>
# Contents

<table><tr><td>1</td><td>Shape deformation of critically charged microdroplets</td><td>1</td></tr><tr><td>1.1</td><td>Introduction</td><td>1</td></tr><tr><td>1.1.1</td><td>The Rayleigh limit of charged droplets</td><td>2</td></tr><tr><td>1.1.2</td><td>How does the deformation develops?</td><td>10</td></tr><tr><td>1.1.3</td><td>Motivation</td><td>14</td></tr><tr><td>1.2</td><td>Experimental setup</td><td>16</td></tr><tr><td></td><td>Injector</td><td>18</td></tr><tr><td></td><td>Taking Snapshots</td><td>19</td></tr><tr><td>1.2.1</td><td>Deformation pathway</td><td>19</td></tr><tr><td>1.3</td><td>Experimental Results</td><td>21</td></tr><tr><td></td><td>Mercury droplets</td><td>21</td></tr><tr><td></td><td>Glycol-water droplets</td><td>23</td></tr><tr><td></td><td>Pure water droplets</td><td>24</td></tr><tr><td>1.4</td><td>Theoretical models and simulations</td><td>27</td></tr><tr><td>1.4.1</td><td>Navier-Stokes equation for Newtonian incompressible fluids</td><td>28</td></tr><tr><td></td><td>Boundary conditions for viscous charged drops</td><td>30</td></tr><tr><td>1.4.2</td><td>Charge dynamics at the interface</td><td>32</td></tr><tr><td>1.4.3</td><td>Kinetic energy dissipation</td><td>33</td></tr><tr><td>1.4.4</td><td>Potential Flow approximation: Case of a perfectly conducting inviscid liquid drop</td><td>34</td></tr><tr><td></td><td>Surface charge distribution for a perfectly conducting liquid</td><td>35</td></tr><tr><td></td><td>Dimensionless equations:</td><td>36</td></tr><tr><td></td><td>Simulated results using the PF approach</td><td>37</td></tr><tr><td>1.4.5</td><td>Viscous Potential Flow (VPF) approximation</td><td>42</td></tr><tr><td>1.4.6</td><td>Viscous Corrected Viscous Potential Flow (VCVPF) approach</td><td>43</td></tr><tr><td></td><td>Dimensionless time evolution of the velocity potential Ψ and surface charge σ.</td><td>44</td></tr><tr><td></td><td>Decay rate: VPF vs VCVPF</td><td>46</td></tr><tr><td></td><td>Charge break up</td><td>47</td></tr><tr><td>1.5</td><td>Simulated Results with VCVPF</td><td>49</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>1.5.1</td>
      <td>Influence of the dimensionless viscosity Oh on the deformation pathway</td>
      <td>49</td>
    </tr>
    <tr>
      <td>1.5.2</td>
      <td>Comparison of VCVPF calculations with experimental results</td>
      <td>52</td>
    </tr>
    <tr>
      <td></td>
      <td>Glycol and supercooled water droplets</td>
      <td>52</td>
    </tr>
    <tr>
      <td></td>
      <td>Glycol-Water mixture</td>
      <td>54</td>
    </tr>
    <tr>
      <td></td>
      <td>Pure water droplets</td>
      <td>55</td>
    </tr>
    <tr>
      <td>1.5.3</td>
      <td>Influence of the charge carrier mobility</td>
      <td>56</td>
    </tr>
    <tr>
      <td></td>
      <td>Positive vs negative drops</td>
      <td>58</td>
    </tr>
    <tr>
      <td>1.5.4</td>
      <td>Boundary element method</td>
      <td>62</td>
    </tr>
    <tr>
      <td></td>
      <td>Simulated Results</td>
      <td>63</td>
    </tr>
    <tr>
      <td>1.5.5</td>
      <td>Finite element method</td>
      <td>64</td>
    </tr>
    <tr>
      <td></td>
      <td>Simulated Results</td>
      <td>66</td>
    </tr>
    <tr>
      <td>1.6</td>
      <td>Summary and Conclusion</td>
      <td>69</td>
    </tr>
    <tr>
      <td>1.7</td>
      <td>Rotationally Corrected Viscous Potential Flow (RCVPF)</td>
      <td>72</td>
    </tr>
    <tr>
      <td></td>
      <td>Choice of the rotational basis set</td>
      <td>75</td>
    </tr>
    <tr>
      <td></td>
      <td>First simulations with RCVPF: Damped surface oscillations</td>
      <td>76</td>
    </tr>
    <tr>
      <td colspan="2"><strong>Appendices</strong></td>
      <td><strong>79</strong></td>
    </tr>
    <tr>
      <td colspan="2"><strong>A</strong></td>
      <td><strong>80</strong></td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Numerical method of the Potential Flow (PF) model</td>
      <td>80</td>
    </tr>
    <tr>
      <td>A.2</td>
      <td>Time propagator (valid also for VCVPF and RCVPF)</td>
      <td>84</td>
    </tr>
    <tr>
      <td>A.3</td>
      <td>Electric field in the case of a non-equilibrium charge distribution</td>
      <td>85</td>
    </tr>
    <tr>
      <td>A.4</td>
      <td>Numerical method for the evaluation of p<sub>μ</sub> in VCVPF</td>
      <td>86</td>
    </tr>
    <tr>
      <td>A.5</td>
      <td>Derivation of the RCVPF model</td>
      <td>87</td>
    </tr>
    <tr>
      <td>A.6</td>
      <td>Numerical method of RCVPF</td>
      <td>88</td>
    </tr>
    <tr>
      <td>2</td>
      <td><strong>Guiding and focusing power of insulating capillaries</strong></td>
      <td><strong>95</strong></td>
    </tr>
    <tr>
      <td></td>
      <td>Acknowledgments</td>
      <td>95</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Electric conductivity in glasses</td>
      <td>100</td>
    </tr>
    <tr>
      <td>2.1.1</td>
      <td>Bulk conductivity</td>
      <td>100</td>
    </tr>
    <tr>
      <td></td>
      <td>Enhanced conductivity due to excess charge carriers:</td>
      <td>102</td>
    </tr>
    <tr>
      <td>2.1.2</td>
      <td>Surface conductivity</td>
      <td>103</td>
    </tr>
    <tr>
      <td>2.1.3</td>
      <td>Non-linear conductivity</td>
      <td>103</td>
    </tr>
    <tr>
      <td>2.1.4</td>
      <td>Interface charge carrier injection</td>
      <td>104</td>
    </tr>
    <tr>
      <td></td>
      <td>Inner surface</td>
      <td>104</td>
    </tr>
    <tr>
      <td></td>
      <td>Outer surface</td>
      <td>105</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Charge dynamics in the capillary walls</td>
      <td>107</td>
    </tr>
    <tr>
      <td>2.2.1</td>
      <td>Excess charge carrier density</td>
      <td>108</td>
    </tr>
    <tr>
      <td>2.2.2</td>
      <td>Surface charge approximation</td>
      <td>109</td>
    </tr>
    <tr>
      <td>2.2.3</td>
      <td>Surface charge dynamics at the interfaces</td>
      <td>110</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>2.2.4</td>
      <td>Evaluation of the electric field at the interfaces</td>
      <td>111</td>
    </tr>
    <tr>
      <td></td>
      <td>Polarization charges in dielectrics</td>
      <td>111</td>
    </tr>
    <tr>
      <td></td>
      <td>Matrix representation of the surface integrals</td>
      <td>113</td>
    </tr>
    <tr>
      <td>2.2.5</td>
      <td>Image charge potential of projectile</td>
      <td>114</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Equation of motion of the projectile</td>
      <td>115</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Charge relaxation rates</td>
      <td>115</td>
    </tr>
    <tr>
      <td></td>
      <td>Examples</td>
      <td>120</td>
    </tr>
    <tr>
      <td></td>
      <td>Dimensionless analysis</td>
      <td>121</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Simulations</td>
      <td>123</td>
    </tr>
    <tr>
      <td>2.5.1</td>
      <td>Numerical model: InCa4D</td>
      <td>123</td>
    </tr>
    <tr>
      <td>2.5.2</td>
      <td>Features of the Numerical model InCa4D</td>
      <td>124</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Most significant simulated results</td>
      <td>130</td>
    </tr>
    <tr>
      <td>2.6.1</td>
      <td>Transverse self-organized focusing</td>
      <td>130</td>
    </tr>
    <tr>
      <td></td>
      <td>The mechanism</td>
      <td>130</td>
    </tr>
    <tr>
      <td></td>
      <td>Explaining experimental results</td>
      <td>130</td>
    </tr>
    <tr>
      <td>2.6.2</td>
      <td>Radial self organized Focusing</td>
      <td>132</td>
    </tr>
    <tr>
      <td>2.6.3</td>
      <td>Influence of Secondary Electrons on high intensity beam transmission</td>
      <td>136</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Experimental Setup</td>
      <td>140</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Remarkable experimental results</td>
      <td>145</td>
    </tr>
    <tr>
      <td>2.8.1</td>
      <td>Depletion (discharge) rate</td>
      <td>145</td>
    </tr>
    <tr>
      <td>2.8.2</td>
      <td>Experimental evidence of self-organized focusing</td>
      <td>146</td>
    </tr>
    <tr>
      <td>2.9</td>
      <td>Conclusion</td>
      <td>150</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Perspectives</td>
      <td>153</td>
    </tr>
    <tr>
      <td>2.10.1</td>
      <td>Stabilization of the focused ion beam</td>
      <td>153</td>
    </tr>
    <tr>
      <td>2.10.2</td>
      <td>Transmission through curved insulating capillaries</td>
      <td>156</td>
    </tr>
    <tr>
      <td>2.10.3</td>
      <td>Transmission of injected electron beam</td>
      <td>158</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td>Appendix : Numerical methods</td>
      <td>159</td>
    </tr>
    <tr>
      <td>2.11.1</td>
      <td>Injected charge at the inner surface</td>
      <td>159</td>
    </tr>
    <tr>
      <td>2.11.2</td>
      <td>Accelerating the charging up</td>
      <td>160</td>
    </tr>
    <tr>
      <td>2.11.3</td>
      <td>Evaluation of the potential and electric field inside the capillary</td>
      <td>161</td>
    </tr>
    <tr>
      <td></td>
      <td>The S matrix</td>
      <td>163</td>
    </tr>
    <tr>
      <td>2.11.4</td>
      <td>Induced surface charges</td>
      <td>165</td>
    </tr>
  </tbody>
</table>
List of Figures

<table><tr><td>1.1 Soap film</td><td>4</td></tr><tr><td>1.2 Ring electrode Paul trap</td><td>5</td></tr><tr><td>1.3 Scheme of experimental setup used by Taflin</td><td>6</td></tr><tr><td>1.4 Size and charge of an evaporating glycol droplet</td><td>7</td></tr><tr><td>1.5 Forced oscillation of an evaporating droplet</td><td>8</td></tr><tr><td>1.6 Eigenfrequencies f<sub>n</sub> of the quadrupole mode</td><td>9</td></tr><tr><td>1.7 Deformation energy of a critically charged droplet</td><td>10</td></tr><tr><td>1.8 super-critically charged droplet deformations</td><td>11</td></tr><tr><td>1.9 Snapshots of the shape evolution of a glycol droplet</td><td>13</td></tr><tr><td>1.10 Snapshots of a supercooled water</td><td>14</td></tr><tr><td>1.11 Photo of the setup in CIMAP</td><td>16</td></tr><tr><td>1.12 Schematic experimental setup</td><td>17</td></tr><tr><td>1.13 Time resolved light scattering intensity</td><td>18</td></tr><tr><td>1.14 Scheme of a piezo-electric injector</td><td>19</td></tr><tr><td>1.15 Superellipsoids</td><td>20</td></tr><tr><td>1.16 Deformation pathway for glycol and supercooled water droplets</td><td>20</td></tr><tr><td>1.17 Snapshots of 60% water and 40% glycol droplets</td><td>22</td></tr><tr><td>1.18 Deformation pathway for mixed liquid droplet</td><td>24</td></tr><tr><td>1.19 Pathway for positively and negatively charged water droplets</td><td>25</td></tr><tr><td>1.20 Snapshots of the instability of positively charged water droplets</td><td>26</td></tr><tr><td>1.21 Calculated shape sequence of inviscid liquid</td><td>37</td></tr><tr><td>1.22 Local fissility and and tip formation</td><td>40</td></tr><tr><td>1.23 Pathway of an inviscid perfectly conducting drop</td><td>41</td></tr><tr><td>1.24 Decay rate VPF vs VCVPF</td><td>46</td></tr><tr><td>1.25 Deformation pathways with different Oh using VCVPF</td><td>50</td></tr><tr><td>1.26 Aspect ratio a/b vs Oh</td><td>51</td></tr><tr><td>1.27 Simulated pathway in VCVPF for Glycol and supercooled water</td><td>52</td></tr><tr><td>1.28 Snapshots with the VCVPF approach</td><td>53</td></tr><tr><td>1.29 Simulated deformation pathway of mixed liquid drop</td><td>54</td></tr><tr><td>1.30 Comparison between experimental and observed water snapshots</td><td>55</td></tr><tr><td>1.31 Influence of charge mobility</td><td>57</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>1.32</td>
      <td>Positively vs negatively charged drops</td>
      <td>59</td>
    </tr>
    <tr>
      <td>1.33</td>
      <td>surface charges as a function of mobility</td>
      <td>60</td>
    </tr>
    <tr>
      <td>1.34</td>
      <td>Representation of the sphere by truncated conical panels.</td>
      <td>63</td>
    </tr>
    <tr>
      <td>1.35</td>
      <td>Simulated time-evolution of the tip shape. x and z are dimension-less, normalized to the initial droplet radius.</td>
      <td>64</td>
    </tr>
    <tr>
      <td>1.36</td>
      <td>FEM discretization</td>
      <td>65</td>
    </tr>
    <tr>
      <td>1.37</td>
      <td>Oscillatory relaxation with FEM.</td>
      <td>67</td>
    </tr>
    <tr>
      <td>1.38</td>
      <td>Pathway using NS and FEM</td>
      <td>68</td>
    </tr>
    <tr>
      <td>1.39</td>
      <td>N vs Oh</td>
      <td>76</td>
    </tr>
    <tr>
      <td>1.40</td>
      <td>Quadupole oscillations in RCVPF</td>
      <td>77</td>
    </tr>
    <tr>
      <td>A.1</td>
      <td>Prolate coordinate system</td>
      <td>81</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>nano- vs macro-capillaries</td>
      <td>96</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Scheme of nano-capillaries</td>
      <td>97</td>
    </tr>
    <tr>
      <td>2.3</td>
      <td>Scheme of tapered-capillaries.</td>
      <td>98</td>
    </tr>
    <tr>
      <td>2.4</td>
      <td>Borosilicate</td>
      <td>100</td>
    </tr>
    <tr>
      <td>2.5</td>
      <td>Excess charge distribution in the capillary wall</td>
      <td>110</td>
    </tr>
    <tr>
      <td>2.6</td>
      <td>Scheme of charges at 3 interfaces</td>
      <td>112</td>
    </tr>
    <tr>
      <td>2.7</td>
      <td>Scheme of the surface current</td>
      <td>116</td>
    </tr>
    <tr>
      <td>2.8</td>
      <td>Weight factors</td>
      <td>119</td>
    </tr>
    <tr>
      <td>2.10</td>
      <td>Secondary electrons</td>
      <td>127</td>
    </tr>
    <tr>
      <td>2.11</td>
      <td>Focusing by external electrode</td>
      <td>129</td>
    </tr>
    <tr>
      <td>2.12</td>
      <td>Transverse field compression</td>
      <td>131</td>
    </tr>
    <tr>
      <td>2.13</td>
      <td>Transmission through tapered capillary as a function of tilt angle</td>
      <td>132</td>
    </tr>
    <tr>
      <td>2.14</td>
      <td>Snapshots of focused trajectories</td>
      <td>134</td>
    </tr>
    <tr>
      <td>2.15</td>
      <td>Radial focusing through short capillary</td>
      <td>136</td>
    </tr>
    <tr>
      <td>2.16</td>
      <td>Ion trajectories through tilted straight capillary</td>
      <td>137</td>
    </tr>
    <tr>
      <td>2.17</td>
      <td>Effect of SE on transmission</td>
      <td>138</td>
    </tr>
    <tr>
      <td>2.18</td>
      <td>Transmission with SE.</td>
      <td>139</td>
    </tr>
    <tr>
      <td>2.19</td>
      <td>Photos of capillary holder</td>
      <td>141</td>
    </tr>
    <tr>
      <td>2.20</td>
      <td>Scheme of capillary holder</td>
      <td>142</td>
    </tr>
    <tr>
      <td>2.21</td>
      <td>Emittance meter</td>
      <td>143</td>
    </tr>
    <tr>
      <td>2.22</td>
      <td>Capillary puller</td>
      <td>144</td>
    </tr>
    <tr>
      <td>2.23</td>
      <td>Simulated discharge curve of total charge in a capillary</td>
      <td>144</td>
    </tr>
    <tr>
      <td>2.24</td>
      <td>Fitted (κ<sub>b</sub>, κ<sub>s</sub>) for some measured discharge rate</td>
      <td>146</td>
    </tr>
    <tr>
      <td>2.25</td>
      <td>Beam spot on the PSD showing the focusing on the injected beam</td>
      <td>147</td>
    </tr>
    <tr>
      <td>2.26</td>
      <td>Simulated and experimental time evolution of the transmitted fraction</td>
      <td>148</td>
    </tr>
    <tr>
      <td>2.27</td>
      <td>Focusing electric field inside the capillary</td>
      <td>150</td>
    </tr>
    <tr>
      <td>2.28</td>
      <td>Influence of the emittance on the transmission</td>
      <td>153</td>
    </tr>
    <tr>
      <td>2.29</td>
      <td>Influence of the grounded tip on the electric potential</td>
      <td>154</td>
    </tr>
  </tbody>
</table>
<table><tr><td>2.30 curved capillary</td><td>156</td></tr><tr><td>2.31 Toroidal coordinate system</td><td>157</td></tr><tr><td>2.32 2D-grid for charge deposition</td><td>159</td></tr><tr><td>2.33 2D-grid for mapping the interior of the capillary</td><td>162</td></tr><tr><td>2.34 Screening of surface charges</td><td>166</td></tr></table>
# Chapter 1

# Shape deformation of critically charged microdroplets

**Collaborations (CIMAP):** J. Rangama, A. Radcliffe*, S. Guillous, B. Manil, B.A. Huber

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

* Post-doctoral position (12 months)

**Collaborations (Ilmenau):** D. Duft, R. Müller, T. Leisner, M. Brickman, C. Wender, T. Achtzehn

*Universität Ilmenau D; Contrat : PROCOP (00353XC, 05927XC)*

**Founding:** ANR "Dynamic" (ANR-07-JCJC-0062-01) 2008-2011

## 1.1 Introduction

If electric charges are present in liquids that are pulverized, the formed droplets may be charged. Such pulverization processes are for example found in oceans or waterfalls and are encountered in the atmosphere. Electrified droplets are also involved in metallic powder production, ink-jet printing, microspray lubrication, emulsification, spray painting and occur in a variety of industrial processes [19]. The droplets range typically from the micrometer to the millimetre. It was found advantageous to create charge droplets that can be guided by electric fields, allowing to optimize the covering of the surface of a given object. If the droplets are pulverized in the presence of an electric field, the amount of charge (and sign) carried by the droplet can be controlled by the electric field. Electroscope ionization, a technique that is revolutionizing the mass-spectrometric analysis of large biomolecules is another application of electrified droplets. Now, if a droplet is
charged to a critical value so that the repulsive electrostatic force overcomes the attractive surface tension force, the droplet breaks up, sending out more charged daughter droplets and leaving behind a stable residual drop. The charge is expected to be emitted from fine jets, triggered by the spontaneous shape deformation of the critically charged drop.

In this manuscript the author will describe the spontaneous shape deformation of critically charged droplets after the onset of the charge instability. The work will briefly introduce the notion of the Rayleigh limit of charged droplets, discuss the experimental setup and observed data sets. In the main part of this first chapter, a detailed discussion of the modelling and simulation of the droplet deformation, from the initially spherical shape to the intermediate shape just before the excess charge is ejected and back to spherical sub-critically charged drop, is proposed. We are in particular interested in the influence of the viscosity and conductivity of liquid on the shape deformation and charge break-up. The project was supported by the ANR "DYNAMIC" 2007, with the author being the project leader.

### 1.1.1 The Rayleigh limit of charged droplets

At equilibrium, the total charge carried by the spherical droplet is uniformly distributed at the surface (similar to a plain spherical conductor), as the ions are free to move in the liquid. This generates a force at the surface that acts outwards, counteracting the inwards force due to the surface tension. The stability of charged droplets has been studied theoretically as early as 1882 by Lord Rayleigh [2]. Considering small deformation amplitudes, he expanded the axial surface deformation $R(\theta, t)$ on Legendre polynomials $P_n$, where $a_m$ is the amplitude of the surfaces mode $m \ge 2$,

$$R(\theta, t) = R_0 + a_2(t)P_2(\cos\theta) + a_4(t)P_4(\cos\theta) + \dots, \quad a_m \ll R_0 \quad (1.1)$$

and calculated the eigenfrequency $f_m$ of each mode $m \ge 2$. He concluded that, in an incompressible charged liquid droplet of radius $R$, the quadupole oscillation becomes unstable against small perturbations as soon as the disruptive Coulomb pressure $P_c = \sigma^2/2\epsilon_0$ due to the surface charges $\sigma$ equals the cohesive capillary pressure $P_s = 2\gamma/R$ due to the surface tension $\gamma$ of the liquid. This result is expressed by the equation of natural frequency $f_m$ of the $m^{th}$ surface mode of a spherical droplet,

$$f_m = \frac{1}{2\pi} \left( \frac{\gamma m(m-1)}{\rho R^3} [(m+2) - 4X] \right)^{1/2} \quad (1.2)$$
where $\rho$ is the density of the liquid. The fissility parameter $X$ stands for the ratio between the two pressures and is defined as <fnref n="1" />,

$$X = \frac{P_c}{P_s} = \frac{\sigma^2/2\varepsilon_0}{2\gamma/R} = \frac{Q^2}{64\pi^2\gamma\varepsilon_0 R^3} \quad (1.3)$$

where we introduced the total charge $Q = 4\pi\sigma R^2$ of the droplet. The last term of (1.3) gives the explicit dependence of $X$ on the droplet radius and shows that the fissility increases with decreasing radius $R$. The frequency equation (1.2) indicates that for a sub-critical fissility, $X < 1$, the modes ($m > 2$) have all a positive frequency and the charged droplet is stable against small perturbations. At the Rayleigh limit, when the fissility becomes critical, $X \ge 1$, the frequency of the quadrupole moment ($m = 2$) becomes imaginary and its mode is unstable against small perturbations, resulting into a droplet that deforms spontaneously. While Lord Rayleigh predicted in his paper [2] the formation of fine jets, he gave no further details about the shape deformation after the onset of the "Rayleigh instability" of a charged droplet.

## Taylor Cone

By charging to several thousand Volt a drop of liquid standing at the lower end of a glass tube, Zeleny observed in 1917 for the first time the disintegration of charged drops [3]. He measured the potential at which they disintegrated owing to the formation of a pointed end from which issued a narrow jet. Working with isolated charged bubbles, Sir G. Taylor (1964) was the first to point out that the drop elongates until it quickly develops an apparently conical end and that a narrow jet appears at the vertex [8]. Moreover Taylor realized that a conical point could exist in equilibrium and gave the condition for its steady state. Namely, the exact balance of capillary pressure and electrostatic pressure on a perfectly conducting conical fluid occurs only for an opening half-angle of 49.3°. Indeed, formation of conical liquid meniscus and associated jet emission are complex fluid dynamical phenomena and received much attention in recent years. Cloupeau et al. [24] (1989) determined the domains of operation in this cone-jet mode, after achieving good reproducibility of phenomena. How does a cone-jet form and how does the thin jet that emanates from the cone's tip break up into small drop was theoretically (numerically) investigated among other by Pantano et al. [30] (1994), by Gañán [33] (1997), by Jens Egger [34] (1997), by de la Mora et al. [48] (2007)

<footnote n="1">¹The fissility, expressed as the ratio Coulomb energy $E_c$ over twice the surface energy $E_s$ of the droplet, $X = E_c/2E_s$, is often used in the literature to express the Rayleigh limit of the droplet. While it is a compact notation, it is not really helpful to understand the origin of the instability and should be depreciated in favour of (1.3).</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1.1: Soap film, microsecond exposures of successive stages (a) of jet formation; (b), (c), (d) subsequent collapse (Taylor 1964).</figcaption></figure>

and by Collins *et al.* [50] (2008). The latter report simulations and experiments that enable a comprehensive picture of the mechanisms of cone formation, jet emission and break-up that occurred during electrodynamic tip streaming from a liquid film of finite conductivity. They propose a scaling law to predict the size of the daughter droplets produced from jet-break-up. More recently Collins *et al.* [57] (2013) gave universal scaling laws for the disintegration of electrified drops, that is, neutral drops in a uniform strong electric field. They used theory and simulation, to show that conductivity can be tuned to yield three scaling regimes for radius and charge of the daughter droplets. Similarly, Gañán-Calvo *et al.* (2016) proposed universal scaling laws of the first ejected daughter droplet from an electrified drop [61]. It ought to be emphasized that the interest is not purely academic as the production of microsized or nanosized jets is of present considerable importance in modern technologies that require a fine control of the size of jets and their decay droplets, including electrospray ionization mass spectrometry (ESI-MS) of large biomolecules [23] as well as printing and coating processes [47]. Electrified droplets were also found to influence the raindrops formation in thunderclouds [13].

Observation of charge emission near the Rayleigh limit

Several experimental studies have attempted to verify Lord Rayleigh's prediction. Hendricks in 1962 produced charged oil drops and accelerated them through a large potential to measure size and charge [6]. In nearly all cases, the charge on the droplet was near or below the Rayleigh limit. Observations of exploding single droplets were made by Doyle *et al.* [7] (1964) using a Millikan-type apparatus.
They examined the stability of evaporating charged Aniline and water droplets. During the evaporation, only neutral molecules are emitted<fnref n="2" />. As a result, the initial charge is conserved until the droplet reach a critical radius, for which the droplet becomes unstable. They reported that at breakup, the parents drop ejected a cloud of very small droplets, which carried away about 30% of its charge but relatively little mass. Schweizer and Hanson [10] (1971) experimented with a charged droplet

<figure><img src="image_6.png" /><figcaption>Figure 1.2: Ring electrode Paul trap</figcaption></figure>

in a quadrupole trap (type Paul trap). This quadupole trap used superposed ac and dc fields to maintain the droplet at the center of the ac ring electrode. The dc field generates a vertical force on the droplet that balances its weight, allowing to extract the mass of the droplet. They reported confirmation of the Rayleigh limit (with a spread of ±4%) for droplets with a radius between 7.5 and 20 µm. The droplets lost 23 ± 5% of their initial charge but only 5 ± 5% of their mass in the explosion. Several other experimental studies reported a range of charge and mass losses. For example, using levitation techniques, Roulleau and Desbois [11] estimated a charge loss in the range of 16% - 40%, but could not detect the mass loss. Rhim and co-workers [20] (1987) at the Jet Propulsion Laboratory (JPL) used a hybrid electrostatic-acoustic technique to levitate charged droplets with diameter of the order of 1 mm. They found that water droplets burst prior to reaching the Rayleigh limit. They were not able to quantify their results, but they speculated that the drop became resonant with the quadrupole field before reaching the Rayleigh limit. Their video camera system was not sensitive enough to detect mass loss during the explosion. In the work of Taflin et al. [22] (1989), micro-droplets were trapped in Paul traps and the size monitored precisely with a laser and a photomultiplier tube mounted a right angle. The "optical resonance spectrum" was used to obtain the size of the droplet to within ±0.01%. They too found that droplets burst prior to reaching the Rayleigh limit with a critical charge only 0.8 times the Rayleigh charge. They investigated if the premature explosion was due

<footnote n="2">²Due to the strong image-charge interaction of charged particles with a dielectric surface, ions are tighter bound to the droplet surface than neutral molecules.</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 1.3: Scheme of experimental setup used by Taflin (1989) and Duft (2003). The photomultiplier at right angle with the laser is not shown.</figcaption></figure>

to thermal fluctuations on the surface, but concluded otherwise. They found it
also not likely that a resonance associated with the ac frequency and the droplet
oscillation could lead to premature fissioning of the droplet as the eigen-frequency
of the droplet with only 80% of the Rayleigh charge would be much higher than the
60 Hz of the ac trapping field. Their preferred explanation is that contamination
of the surface could lead to regions of low surface tension, which would reach the
Rayleigh limit before a pure component droplet would do so. They also noted
that with each explosion, the mass loss was within 1% – 2.3%, while the charge
loss ranged between 10% and 18%. They found however no correlation between
mass and charge loss in the reported data. Using electrodynamic balance for drop
levitation, Richardson *et al.* [1] observed a 15% charge loss for dioctyl-phthalate
droplets and around 50% for sulphuric acid drops. The mass loss observed in the
case of both liquids was less than 2.5%. Gomez and Tang [29] measured the charge
and size of heptane droplets generated by electrostatic sprays and observed rupture
of droplets, when they had charging levels between 70% and 80% of the Rayleigh
limit. On the whole, experimental results estimate the charge loss in the range
of 15% - 50% with less than 5% mass loss. The discrepancy observed in different
experiments may be attributed to the errors encountered in the measurements due
to the very fast dynamics of the process as well as the influences of other parameters
such as applied electric field, acceleration due to gravity, and aerodynamic drag.
There is however an agreement that the droplets seem to burst prior to reaching
the Rayleigh limit with a critical charge only 0.8 times the Rayleigh charge.

More recently, Duft *et al.* [38] (2002) report direct observation of the breakup of
<figure><img src="image_8.png" /><figcaption>Figure 1.4: Size and charge of an evaporating glycol droplet. The dotted region around the first Coulomb instability is enlarged in the inset (Duft 2002).</figcaption></figure>

charged micro-droplets where the amplitude and phase of the forced quadrupole shape oscillation of the droplets are taken as probes for their fissility. Indeed, at resonance (X=1), the phase between the droplet oscillation and the external quadupole AC field is shifted by $\phi = \pi/2$. They found that for the glycol droplet investigated, the breakup takes place for $\phi \approx 1.0$ rad which corresponds to $X \approx 0.97$ [14], independent of the size of the droplet at the moment of instability. This was the first clear indication that for evaporating highly charged microdroplets of glycol, the instability occurs indeed close to the Rayleigh limit, confirming Rayleigh prediction. This is in contrast to the previous investigations [6, 7, 10, 20, 22, 29] which report a critical fissility between 0.5 and 0.8, when using the tabulated surface tension of the droplet liquid. However, using the tabulated surface tension of glycol, Duft et al. found an apparent critical fissility of $X \approx 0.7$ as well. They noticed that the apparent fissility decreases with subsequent instabilities of the same droplet and therefore tentatively attribute the reduced surface tension to surface active contaminants, which become more and more concentrated as the droplet shrinks<fnref n="3" />.

This prediction is also corroborated by our simulations of a low viscous charged

<footnote n="3">³At this point, one may also ask in how far the ions that accumulate at the surface influence the surface tension of the liquid droplet.</footnote>
<figure><img src="image_9.png" /><figcaption>Figure 1.5: Simulated forced oscillation of an evaporating droplet (given by the aspect ratio, major axis over minor axis, (a/b) as a function of time in a quadrupole AC field. The time-dependence of the fissility $X(t)$ is due to the decreasing droplet radius in time. The amplitude and mean value of the driving AC field are scaled and shifted for better readability of the phase. The spontaneous spheroidal deformation of the droplet takes place for $X \approx 0.997$ with a phase shift of about $\pi/2$, as highlighted by the vertical black line. The simulation used a quadrupole AC field of 200 Hz having a potential amplitude at the droplet surface of 40 mV.</figcaption></figure>

droplets in an external oscillating quadrupole field of the Paul trap. Simulations show that the forced oscillation in a Paul trap will not trigger the spontaneous deformation until the fissility of the droplet gets close to the Rayleigh limit. In figure 1.5, we simulated the droplet oscillation in an AC quadrupole field. The frequency of the AC field and its potential at the droplet's surface were chosen to reproduce the experimental conditions of the setup used by Duft (2003). The electric potential of the external AC field drives the mobile charges at the droplet's surface, inducing a quadruple Coulomb pressure at the surface, making the quadrupole mode of the surface oscillate at the frequency of the AC field. Details about the simulations will be given in section 1.4. Due to evaporation, the radius $R(t)$ of the droplet decreases with time. As a result the fissility $X(t)$ of the droplet increases (1.3), while the eigenfrequency of the quadrupole moment $f_2(t)$ decreases. Using
<figure><img src="image_10.png" /><figcaption>Figure 1.6: Eigenfrequencies $f_n(X)$ of the quadrupole mode (violet) and octupole mode (green) as a function of the fissility $X$, normalized with respect to the quadrupole frequency $f_2(0) = 8.5$ kHz of a neutral ($X = 0$) water droplet of $R_0 = 100$ µm. The frequency of the quadrupole field of the trap is typically about 200 Hz, thus $f_{AC}/f_2(0) = 0.0176$. In this case, the quadruple eigenfrequency $f_2(X)$ of the given water droplet crosses the AC frequency at $X = 0.997$.</figcaption></figure>

equation 1.2, the time-evolution of the eigen-frequency may be written,

$$f_2(t) = f_2^0(t) \left(1 - \frac{R_c^3}{R^3(t)}\right)^{1/2} = f_2^0(t) (1 - X(t)) \quad (1.4)$$

where $f_2^0(t) = \frac{1}{2\pi}\sqrt{\frac{8\gamma}{\rho R^3(t)}}$ is the quadrupole eigen-frequency of a neutral droplet of radius $R(t)$ and $R_c = (Q^2/(64\pi^2\gamma\epsilon_0))^{(1/3)}$ the critical radius. Once the droplet radius approaches the critical radius, $f_2(t)$ decreases quickly and reaches the frequency $f_{AC}$ of the AC field in the Paul trap (see figure 1.6). At that point, the amplitude of the quadrupole mode increases due to the resonance, triggering the spontaneous spheroidal deformation of the droplet (see figure 1.5). This also means that, out of resonance, the amplitude of the forced oscillations of the droplet in the Paul trap is so small that it cannot trigger the spontaneous deformation of a droplet. Only at resonance and only for $X \simeq 1$ can the spontaneous deformation of a droplet be triggered by the AC field of the Paul trap. In agreement with Taflin [22], this study discards thus the assumption that the AC field may trigger the spontaneous deformation of the droplet for fissilities $X$ well below unit.
### 1.1.2 How does the deformation develops ?

While the eigenfrequency analyses of the first multipole modes of surface oscillations tells us that the quadupole mode is the first to become unstable, it does not tell us how the droplet deforms initially. Using the deformation energy for small quadrupole amplitudes [5],

$$W(a_2) = \frac{4}{5}(1-X)a_2^2 + \frac{3}{5}(X - \frac{3}{2})a_2^3 + \dots, \quad (1.5)$$

we can discuss the initial spheroidal deformation at the Rayleigh limit. The deformation energy has been plotted in figure 1.7 as a function of the amplitude $a_2$ and for three different fissilities. For $a_2 > 0$, the droplet is prolate spheroidal and for $a_2 < 0$ the droplet is oblate spheroidal, with the axis symmetry as indicated. For $X < 1$, the energy potential wall is parabolic with positive curvature and the droplet shape oscillates between oblate and prolate. For $X = 1$, the quadratic term in $W$ is zero and the cubic term dominates. As a result, for $a_2 > 0$, the barrier has been suppressed and the droplet deforms spontaneously towards prolate spheroidal shapes. As we will see, prolate shape deformations are indeed observed in the experiments. Interestingly, for $X > 1$, the oblate deformation space is now also accessible and the droplets can become oblate. The latter deformation mode can be enforced by pre-deforming the droplet into an oblate shape when reaching

<figure><img src="image_11.png" /><figcaption>Figure 1.7: Deformation energy of a critically charged droplet as a function of the amplitude of the quadupole mode. For a negative amplitude, the deformation is oblate, while for positive it is prolate. The deformation energy is given for subcritical ($X < 1$), critical ($X = 1$) and super-critical ($X > 1$) fissilities.</figcaption></figure>
the Rayleigh limit. However, in the Paul trap, the evaporating droplet is initially sub-critical $X < 1$ and the fissility approaches the limit from the lower end. As a result, the droplet always deforms spontaneously towards a prolate spheroidal, just before crossing the $X = 1$ limit.

<figure><img src="image_12.png" /><figcaption>Figure 1.8: Equilibrium shape of a super-critically charged droplet as a function of the fissility. Deformation is assumed axisymmetric.</figcaption></figure>

Note that in the super-critical case ($X > 1$), the deformation energy shows a minimum for a negative amplitude. I would like to stress that in this analysis, the drop deformation is assumed axisymmetric. The existence of a minimum in the deformation energy means that one can expect the super-critical droplet to stabilize in an oblate shape and that the equilibrium shape depends only on the fissility. Simulations performed with our numerical code, based on the approach given in section 1.4.4, indeed confirmed that the droplet stabilizes in an oblate-like shape. The higher the fissility the flatter the droplet (see figure 1.8). For fissilities above 2, the shape tends towards a disk, with the charge located at the border, where the curvature is largest. The insets of figure 1.8 shows the disk-like droplet for two different fissilities, namely $X = 1.4$ and $X = 2.7$. For larger fissilities, the simulations were not accurate enough to be predictive. Similar results, based on axisymmetric Navier-Stokes calculations, were presented in [58]. The trend suggests that, with increasing fissility, the droplet ends up in a doughnut shape,
as it would reduce the surface energy with respect to the disk shape. While oblate deformations can be calculated numerically for super-critical droplets by enforcing an initial oblate shape, it is much harder to observe them experimentally as the droplet is usually initially sub-critical and deforms spontaneously into the prolate domain when reaching the critical radius due to evaporation. One could however image the case where the charge of a slowly oscillating drop in a Paul trap is suddenly increased above the Rayleigh limit by an intense ion beam so that it becomes super-critically while being oblate deformed<fnref n="4" />. Another reason why such oblate shape may never be observed is that the droplet can deform by breaking the axial symmetry of the drop and elongate along an arbitrary direction orthogonal to the symmetry axis. As our code imposes axial symmetry, such calculations were not performed. Nevertheless, it would be interesting to calculate the stability of the equilibrium shape against small surface perturbations in a future work.

## Observation of the prolate shape evolution

While the work of Duft et al. (2002) [38] confirmed the Rayleigh prediction of the instability of charged droplet close to the Rayleigh limit, it was still unclear how exactly the drope deforms and how exactly the charge was emitted from the droplets. In a pioneering work, Duft and coworkers [39] published in Nature (2003) the observed shape evolution of a critically charged ethylene glycol droplet, by taking snapshots of a fast deforming droplet at the Rayleigh limit (Fig. 1.9). The charged droplet was trapped in a Paul trap and the radius of the evaporating droplet was monitored. From the vertical position of the droplet in the trap, the ratio charge over mass was deduced. From the radius, the initial charge was calculated. A scheme of the setup is given in figure 1.11 and details about the experimental techniques are found in [38, 39, 42, 49]. A similar setup is further discussed in section 1.2.

The injected charged drop is stabilized in the middle of the Paul trap, while it's fissility increases with ongoing evaporation. Once the eigenfrequency of the quadrupole moment of the surface oscillation of the charge drop reached the frequency of the AC field of the Paul trap, the droplet deforms spontaneously into a prolate form. The onset of the deformation is tracked by the photo-multiplier which triggers, with a predefined delay, the flash lamp. This way the snapshots are taken at different deformation times. Putting the time origin at the moment the jet appears, the first snapshot was taken 150 µs before the time origin, when the droplet had already a clear prolate shape. The droplet deforms fast into a rugby-ball shape and then pointed ends appear at the tip. The charge concentrates at

<footnote n="4">⁴Note that the author has not evaluated if such a fast charge accumulation is possible with ion beams.</footnote>
<figure><img src="image_13.png" /><figcaption>Figure 1.9: Snapshots showing the shape evolution of a glycol droplet with a radius of 24 µm after the onset of the instability. Droplet's temperature when injected was (85° C). During evaporation in an atmosphere at room temperature, the liquid cools down so that the temperature, at the moment of instability, is unknown (but estimated close to room temperature). Jet emission is chosen as origin of time. The second jet at the bottom of the droplet is not shown. Extracted from Giglio et al. (2008) [49]</figcaption></figure>

the tips (because of the point effect) increasing locally the coulomb pressure. A jet is formed and charged matter is ejected from both tips. The fine jets of matter together with the concentrated charge at the tip explain intuitively the ~35 % charge loss and only ~2 % mass loss during the process. The total charge of the droplet has decreased by a certain amount and the fissility becomes again sub-critical. The jet breaks up after several microseconds and the droplet returns to its initial spherical shape. With ongoing evaporation the droplet reaches again the Rayleigh limit, triggering a new deformation and jet emission, until the drop has evaporated completely. The series of snapshots allow to realize a "film" of the charge break-up. Apart from Glycol drops, the group of T. Leisner also investigated water droplets and in particular supercooled (−5° C) pure water droplets, which are found in thunderclouds. Figure 1.10 shows snapshots of the shape evolution of such a supercooled water droplet. Note that while the shape evolution is
<figure><img src="image_14.png" /><figcaption>Figure 1.10: Snapshots showing the shape evolution of a supercooled water droplet (-5° C) with a radius of 35 µm after the onset of the instability. Extracted from Giglio *et al.* (2008) [49]</figcaption></figure>

in general similar to the one for Glycol, it differs in detail, motivating thus further investigations.

### 1.1.3 Motivation

The aim of this work is to describe numerically and experimentally the shape evolution of a critically charged droplet. The topic arose from a collaboration of a group of CIMAP supervised by B. Huber and a group of the university of Karlsruhe supervised by T. Leisner. It was found experimentally that the deformation pathway as well as the jet emission of glycol droplets differs in detail from supercooled water droplets [49]. For example, the elongation of the glycol droplet just before jet formation was larger than for supercooled water droplets, which in turn was larger than water droplet at room temperature. In particular, for water droplets at room temperature, no jet emission was observed but charge emission of more than 20 % of the initial charge was still measured, hinting to an extremely thin jet and/or super-fast build up and relaxation of the jet during charge emission, so that the probability of getting a snapshot of the jet was too small. Also, for low viscous water droplets, the return to the initial spherical shape went through strongly deformed intermediate droplet shapes, which were not observed for glycol drops.

On the theoretical front, no reliable simulations of droplet deformations existed so far except for the work of Belelú *et al.* who presented a preliminary study of the viscous droplet deformation in the creeping flow approximation (Stokes flow). However this approximation may break down for the ultrafast low viscous droplet deformation, for which the advective inertial forces may dominate the viscous
forces. We propose in this work to model and simulate the droplet break-up of a
critically charged droplet, from the initial spherical shape to the jet formation and
back to the spherical shape. This allows to highlight the role of the different liquid
parameters like density, viscosity, surface tension and conductivity of the liquid on
the shape deformation and jet formation. We want a model that is able to predict
with high accuracy the different shapes of the deforming droplet and explains in
how far the conductivity or viscosity of the liquid affect the details of the droplet
deformation and jet formation. Analysing the tip shape just before jet formation,
will confirm if the jet is emitted from a Taylor cone as for electrified liquids [8].
Further it gives the intensity of the electric field at the tip at charge break-up
and. Finally details about the jet formation and it's size can be studied. It ought
to be emphasized that the interest in not purely academic as the production of
micro-sized or nano-sized jets is of present considerable importance in modern
technologies that require a fine control of the size of the jets and their decay
droplets. For this purpose, we developed numerical codes able to simulate the
axisymmetric deformation of a charged drop in the presence of an external field
(for example the AC field of the quadrupole trap). The model is based on various
approximations of the Navier-Stokes equation and will be discussed in section 1.4.

On the experimental front, we propose to build our own setup that allowed
us to inject droplets of various liquids, trap them and observe their spontaneous
deformation once it reached the critical fissility. The motivation is to produce new
experimental results which provide guide lines (constraints) to theoretical models.
A particular case represent mercury droplets which have high conductivity and
high density. They represent thus model systems for theoretical study as the
droplet surface is an equipotential and the dynamic viscosity of 1.5 ctP is rather low
(comparable to water) so that simplified numerical models can be used to simulate
the deformation of mercury droplets. Other interesting liquids are water droplets
with additional salt or acid to increase their conductivity. Or a mix of glycol and
water so as to modify finely the viscosity of the liquid. Finally we want to change
the polarity of the charged droplets to study the effect of the conductivity on the
droplet deformations and jet emission as the mobility of positive and negative
charge carriers in the same liquid may differ significantly.

The topic was founded by the ANR 2007 program (ANR-07-JCJC-0062-01,
label "Dynamic") for 3 years for an amount of 135 000 €, including a post-doctoral
position of 12 month that was attributed to D. Allister Radcliffe. The program
founded a new experimental setup for injecting, trapping and observing droplet
deformations of different liquids via highly resolved snapshots.
1.2 Experimental setup

The experimental setup used for trapping and taking snapshots of the droplet is similar to the one presented in [39, 42, 49]. In a nutshell, 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 camera. Delaying the trigger allows taking snapshots at different times of the shape deformation. We give below some details

<figure><img src="image_15.png" /><figcaption>Figure 1.11: Photo of the setup in CIMAP. Vertically above the quadrupole trap is the injector for micro-drops. Flash-lamp, CCD camera photo-multiplier and optics for the laser beam are shown. Credit to Jimmy Rangama (CIMAP).</figcaption></figure>

about the experimental setup.

Individual droplets are generated on demand with a piezo-driven nozzle, charged by influence from a high voltage electrode and injected into a electrodynamic balance operating at ambient pressure. The electromagnetic levitator has been generously borrowed by the group of T. Leisner of Ilmenau [38, 39, 42]. The levitator follows a well proven design [17, 21] and consists of two ring electrodes (inner diameter, 10 mm; vertical distance, 2 mm) which are encapsulated in a chamber. The chamber is machined from massive brass in a shape of an octagon and carries the necessary ports for various optical and electrical probes. In order to levitate the charged droplets, an AC voltage of about 2 kV amplitude and a frequency of
200 Hz is applied to both ring electrodes. The resultant electrical quadrupole field is rotationally symmetric around the vertical axis. For a charged particle with a large charge over mass ratio, the AC field would be sufficient to trap the particle in the center of the trap. For charged micro-droplets however, the ratio charge over mass is too low and a static vertical field for the compensation of the gravitational acceleration on the droplet is needed. The latter is created by superimposing a DC voltage of opposite sign to the top and bottom electrode. The droplet is illuminated by a NeHe laser beam in the horizontal plane and imaged onto a vertically oriented linear CCD array in order to determine its vertical position in the trap. The CCD array then pilots the intensity of the DC field to keep the droplet centred in the trap, via a feed-back regulator. From this voltage, the mass to ratio the of the droplets can be determined with a time resolution of 100 ms. This value serves to adjust the frequency of the AC voltage to assure stable trapping as the droplet evaporates. Under typical conditions, this frequency varies between 200 Hz and 500 Hz. The Laser light scattered from the droplet under 90° is collected

<figure><img src="image_16.png" /><figcaption>Figure 1.12: Schematic experimental setup (cross section through the levitator in the horizontal plane. Unlike the setup used by Duft et al. (2003) the injector is placed vertically above the quadrupole trap. Extracted from Giglio et al. (2008) [49]</figcaption></figure>

by a second microscope objective and imaged onto a photo-multiplier tube. Its signal serves to control the amplitude of the forced quadrupole oscillation of the droplets in the field of the trap. It also serves to track the onset of the instability. Indeed, when the droplet becomes prolate spheroidal, the intensity of the reflected light increases. The fast increase of the signal provides the origin of the time at the onset of the instability which serves to trigger the flash lamp and the CCD camera after a predefined delay for the ultrafast photography of the droplet shape.
<figure><img src="image_17.png" /><figcaption>Figure 1.13: Time resolved light scattering intensity of a levitated droplet around the break-up. Red dashed line gives the trigger level for fixing the time origin. Extracted from Giglio et al. (2008) [49] and modified to include the trigger level.</figcaption></figure>

## Injector

The shape deformation taking place in typically less than 300 µs, the setup can take only one single snapshot of the deforming droplet. The sequence of snapshots given in figure 1.9 and 1.10 represents thus each time a different droplet. In order to have a coherent sequence, the size and charge of the injected droplet as well as its injection velocity must be perfectly reproducible. We chose therefore a piezoelectric injector like the one presented in figure 1.14. The integrated piezoelectric actuator induces a shock-wave into the fluid contained in the head, which causes a droplet to be emitted from the nozzle. The experimental results presented in section 1.3 where obtained using a piezoelectric microdispensers from GeSim. The latter was limited to liquids having a viscosity of less than 5 mPa.s and was not heatable. The liquid was grounded and the injector positioned 2 mm on top of a centred ring electrode, which is polarized so as to create an electric field of about 1 kV/mm between the electrode and the outlet of the dispenser. Under the influence of the electric field, the charges in the jet produced by the injector are separated. When the drop separates from the jet, the drop has accumulated an excess charge. With the computer controlled inject-unit, droplets of same size and
charge are produced with high reproducibility. Once produced, the droplets fall into the trap 15 mm below where they are immobilized by the electric field. It should be emphasized that compared to the setup of our colleagues in Ilmenau, we lacked the possibility to heat the injector as well as heat or cool the Paul trap.

<figure><img src="image_18.png" /><figcaption>Figure 1.14: Scheme of a piezo-electric injector. (1) piezo-driven plates, putting a pressure on the liquid, injecting liquid from the nozzle; (2) Grounded electrode; (3) Glass tube, (4) liquid.</figcaption></figure>

## Taking Snapshots

The trapped droplets have a diameter between 20 to 80 µm. The instability lasts between 20 and 400 µs. For image acquisition, we chose therefore to use a high resolution CCD camera (PCO 2000) which has an ultra-fast shutter (~ 5 µs) equipped with an objective ×20. The light source for illuminating the droplet is a flash lamp NANOLIT (High-Speed Photo-System). The lamp has the particularity of emitting a light flash which lasts about 25 ns for an energy of 9 mJ. The set, "CCD/flash" is triggered during the instability. By varying the delay time after the photo-diode has detected the onset of the instability, a series of time ordered snapshots are obtained. The size of the droplet can be determined from the microscope images within an error of ± 0.5 µm from the first prolate shaped drop or alternatively from the final spherical drop.

### 1.2.1 Deformation pathway

In order to compare easily the deformation pathways of different droplets, the shapes have been fitted with the mathematical expression of a superellipsoid of the form

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

The shapes have been fitted with the constraint of constant volume so that the fit had only two independent parameters, namely a/b and n. Superellipsoid can either fit capsule or spindle-like forms, depending on n. The aspect ratio a/b gives the elongation of the droplet. Adjusting (a/b, n) for the droplet shape for each
<figure><img src="image_19.png" /></figure>

Figure 1.15: Superellipsoids have a mathematically simple expression (eq. 1.6) and can fit most of the observed droplet shapes with just two degrees of freedom, $(a/b, n)$, from pointed rugby-ball to spheroid to capsule with increasing $n$ (generated with the Superellipsoid function in Mathematica (Wolfram))

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

Figure 1.16: Deformation pathway for glycol (black) and supercooled water (violet) droplets. The aspect ratio $a/b$ and exponent $n$ correspond to a shape fitted with a superellipsoid. The arrows give the direction of the time evolution of the path. Extracted from Giglio et al. (2008) [49]
snapshot yields a number of points that can be plotted on an graph with a/b as the abscissa and n as the ordinate, so that the whole droplet deformation path can be represented reasonably well on a 2D graph. In figure 1.16, we give the deformation pathway of glycol and water droplets, where the snapshots have been shown in figure 1.9 and 1.10 respectively. Fitting the shapes by a superellipsoid allows to compare the deformation pathway of droplets having different liquid properties. First, we note that the general behavior is similar in both cases: the instability starts with increasing a/b, while n remains nearly constant at ≈ 2, meaning that the initial deformation remains mainly prolate spheroidal. Then n decreases and the form becomes spindle-like. The charge emission however occurs at a/b ≈ 3.2 and n ≈ 1.4 for water drops and at a/b ≈ 3.5 and n ≈ 1.6 for glycol. After emission, the tip flattens, while n increases quickly well above 2 to fit an elongated capsule, before both parameters come back to their initial values. The return can be monotone or oscillatory, depending on the viscosity of the fluid. Indeed, the water droplet exhibits a double oscillatory behaviour between an oblate and prolate ellipsoid and between sharp and flattened ends. For glycol, the amplitudes of this motion are much smaller than for water, as dissipation is much stronger. The (a/b, n) representation of the droplet shapes is thus a handy tool to analyse in detail the spontaneous shape evolutions of droplets at the Rayleigh limit of liquids with different physical properties. 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 using a software, written by Pierre Salou (at that time PhD student in CIMAP), and then the obtained 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.

1.3 Experimental Results

We wanted to explore the charge break-up of droplets for a large range of viscos-
ity, density, surface tension and conductivity. We wanted to see the influence of
those parameters on the deformation pathway and jet formation as well as charge
emission. The idea was also to gather experimental data that are constraints for
the numerical model presented below.

Mercury droplets

One interesting class of liquids are metallic liquid like mercury, which are electrically conducting and have a low viscosity of about 1.5 mPa.s (1.5 cP), see Table
<table><caption>Table 1.1: Physical properties of various liquids (1 cP = 1 mPa.s).</caption><thead><tr><th></th><th>density (kg/m³)</th><th>surface tension (N/m)</th><th>viscosity (cP)</th></tr></thead><tbody><tr><td>water (20° C)</td><td>1000</td><td>0.072</td><td>1</td></tr><tr><td>water (−5° C)</td><td>1000</td><td>0.076</td><td>2.2*</td></tr><tr><td>glycol (25° C)</td><td>1102</td><td>0.048</td><td>16</td></tr><tr><td>mercury (20° C)</td><td>13500</td><td>0.48</td><td>1.5</td></tr></tbody></table>

*Viscosity of pure water at −5° C is obtained from "Viscosity of deeply supercooled water and its coupling to molecular diffusion", A. Dehaoui, B. Issenmann, and F. Caupin, PNAS September 29, (2015) 112 (39) 12020-12025; https://doi.org/10.1073/pnas.1508996112.

1.1. Mercury droplets should behave closely to perfectly conducting inviscid liquid drops, a class of ideal liquids. The emitted charge carriers in the jet of positively and negatively charged mercury droplets during charge breakup is also fundamentally different so that metallic droplets are a class of their own. We tried thus to inject and trap mercury micro-droplets in our Paul trap. Unfortunately, after several months of testing we did not succeed in trapping mercury droplets. The problem lay mainly within the injector. The meniscus of the mercury liquid at the outlet of the injector tent to oxidize quickly, blocking further injection. We

<figure><img src="image_21.png" /><figcaption>Figure 1.17: Snapshots of the deforming, droplets (ø = 52µm) after the onset of instability. Droplet composition was 60% water and 40% glycol, injected at room temperature.</figcaption></figure>
tried two different injectors: A home-maid $\emptyset = 80\mu m$ injector graciously bor-
rowed by the group of T. Leisner and a microdrop dispenser (MD-K-130) with
a non-heated $\emptyset = 70\mu m$ nozzle head from microdrop. At the end, both injec-
tors were quickly made unusable by the mercury liquid. Also, the weight of the
droplets was such that they fall simply through the trap, even when a DC field
was present to balance the weight of the drops. Finally, the stringent regulations
and constraints associated with the manipulation of toxic mercury made the drop
injection and trapping a difficult task. We had to move the whole setup under a
chemical hood, adding some delay and increasing the manipulation difficulties. As
we never succeeded to repetitively trap mercury droplets, we give up the study of
critically charged mercury micro-droplets. Note that another non-negligible factor
that makes the observation of mercury droplet a difficult task is its low vapour
pressure of 0.002 mbar, reducing the evaporation rate by a factor 1000 with respect
to water droplets at room temperature, meaning that a trapping time of several
hours may be necessary before the droplet reaches the Rayleigh limit.

Glycol-water droplets

We then switched to glycol-water droplets. By varying the percentage of water in
a glycol-water mixture, we were able to alter finely the viscosity of the droplets
between 1 and 20 cP. However, here we were limited by our piezo-electric injector,
which is able to work only with liquids having a viscosity lower than 5 cP. We
succeeded to inject and trap droplets at room temperature, having a mixture of
60% of water and 40 % of glycol. From the data sheet Ethylene Glycol-MEGlobal
<fnref n="5" /> we get the viscosity and surface tension of the mix at 23° C, namely $\mu_{mix} \approx 2.3$
cP and $\gamma_{mix} \approx 0.058$ N/m. We show in figure 1.17 snapshots of the deforming
droplet at different times after the onset of the instability. When injected, the
trapped droplets have an initial radius of 34 $\mu m$. We estimate the density $\rho$ of
the mix to be well approximated by the weighted mean value of water and glycol,
$\rho_{mix} = 1060$ kg/m³. At the onset of the instability, we found the radius to be about
$R = 26 \mu m$, meaning that the droplet carried a charge of about $Q = 4\pi\sqrt{\gamma_{mix}R^3\varepsilon_0}$
$\approx 1.2$ pC. We see that 137 $\mu s$ after the onset of the instability, the droplet emits
a jet of charged matter. We remind that the times below the snapshots give time
intervals relative to a trigger rather than absolute times. Indeed, they are taken
relative to a threshold deformation that triggers the flash lamp. As the threshold
value is put manually to an arbitrary value between 1.2 < a/b < 1.5, the times
may vary from one experiment to another. At the moment of jet emission, the
elongation ratio a/b is about 3.2, similar to the one found for supercooled water
(Fig.1.10 and Fig.1.18). After jet emission, the droplet returns to the spherical

<footnote n="5">⁵https://www.meglobal.biz/products-and-applications/product-literature/</footnote>
<figure><img src="image_22.png" /><figcaption>Figure 1.18: Deformation pathway for droplet of mixed liquids water 60% and glycol 40 % The aspect ratio *a/b* and exponent *n* correspond to a shape fitted with a superspheroid.</figcaption></figure>

shape, with an overshooting at 200 µs into the oblate spheroidal domain. The deformation pathway of the water-glycol droplets is illustrated in figure 1.18, where is compared to the pure glycol droplets. While the initial pathway ($a/b < 2$) are similar, they differ quickly with increasing elongation. In particular, the pointed tips form earlier, as is shown by the parameter *n*, which decreases faster with $a/b$ for the mixed liquid. Also, the mixed liquid droplet becomes oblate shaped, $a/b < 1$, on its "return" path, unlike the pure glycol droplet.

## Pure water droplets

We also studied pure water droplets at room temperature (22° C). De-ionized water was used and graciously made available by Yvette N'Gono from CIMAP laboratory. Inverting the polarity of the polarization ring of the injector, we can inject negatively charged droplets. Of course, the DC field that balances the weight
of the drop needed to be inverted accordingly. In figure 1.32 we compare the
deformation pathway of negatively and positively charged droplets. More than
200 snapshots were taken and fitted by a superellipsoid. 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 [39, 49]. 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 µm. These observations are corroborated by the experimental
results presented in the Ph.D work of R. Müller [52]. Full black squares stand

<figure><img src="image_23.png" /><figcaption>Figure 1.19: Deformation pathway for positively and negatively charged water droplets.</figcaption></figure>
for positively charged while red squared stand for negatively charged pure water droplets. Both pathways are similar but differ in detail. We found that positively charged drops are more elongated (have a larger aspect ratio) at charge emission. Also, while both go through strongly deformed oblate shapes before returning the equilibrium (spherical) shape, the return path of positively and negatively charged droplets seems shifted by $\Delta a/b \simeq 0.4$. At this point we mere note that what differs positively from negatively charged droplets is the mobility of their respective excess charge carriers. For positively charged drops, the excess charge carrier is the proton or more precisely the hydronium cation $H_3O^+$ with a mobility of $36 \times 10^{-8} \text{ m}^2\text{s}^{-1}\text{V}^{-1}$, while for negatively charged water drops it is the hydroxyl anion $OH^-$, with a mobility of $20 \times 10^{-8} \text{ m}^2\text{s}^{-1}\text{V}^{-1}$, thus almost two times slower [53, 32]. The influence of the charge mobility will be studied in detail in section 1.5. Below, in figure 1.20, we give snapshots of positively charged pure water droplets at room temperature are shown. Having a low viscosity, water droplets show strongly deformed shapes after jet emission, with shapes ranging from "diamond" to "peanut" shapes and from oblate to prolate shapes.

<figure><img src="image_24.png" /><figcaption>Figure 1.20: Snapshots of the instability of positively charged water droplets, injected at room temperature. Note the absence of jets and daughter droplets.</figcaption></figure>
1.4 Theoretical models and simulations

Using simulations, we want to understand in more detail how the viscosity and elec-
trical conductivity influence the shape deformation of critically charged droplets.
Some simulations of critically charged droplet deformations have already be per-
formed by several authors. We may cite the work of Belelú et al. [46] and Gawande
et al. [62] who presented a study of the shape deformation of a charged viscous
perfectly conducting drop suspended in air. Using the boundary-element method
(BEM) 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 [38, 39]. However the Stokes flow approximation may
break down for the ultra-fast low-viscous droplet deformation, for which the ad-
 vective inertial forces may dominate the viscous forces. The case of an inviscid
perfectly conducting liquid droplets was discussed by Giglio et al. in [49] 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 defor-
mations. Burton and Taborek [55] 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.
Alastair Radcliffe, who occupied a post-doc position in 2010 in our group SIMUL,
used the finite element methods (FEM) to solve the Navier-Stokes (NS) equation
and to simulate the deformation of perfectly conducting viscous droplets. His
preliminary results [58] will be discussed in section 1.5.5. He later improved the
algorithms and succeded to simulate the drop deformation up to the jet formation
[59]. Related to our study are also the work of Nganguia et al. [60], who studied
the inertia effects on the electrode formation of a viscous drop under a DC electric
field. Collins et al. [50] used FEM with elliptic mesh generation and adaptive im-
plicit time integration to solve the NS equations that described the jet formation
of electrified liquids. They discussed the profile of the jet as a function of the liquid
properties and found that conical (Taylor) tips can occur with perfectly conduct-
ing fluids because charge transport by ohmic conduction is instantaneous so that
the electrostatic pressure continually balances the increasing capillary pressure.
However, for a leaky-dielectric (low conductivity) fluid, jets develops without go-
ing through a conical cusp. In a later work, their simulations unequivocally show
that electrospray daughter droplets, that is droplets that where formed by the
jet of an electrified parent droplet, are coulombically stable at the instant they
are created and that there exists a universal scaling law for the daughter droplets
charge [57]. Alternatively, Gañán-Calvo et al. used the promising Volume of Fluid
(VoF) method for tracking and locating the free surface (fluid-gas interface) [61].
They succeeded to simulate 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, that were validated both numerically and experimentally.

The main aim of this work is to study the shape deformation of a charged droplet as a function of the liquid properties. We propose thus to develop a numerical code, with the requirements that it simulates the whole deformation path of a critically charged droplet, from the onset of the Coulomb instability up to charge breakup and back to the initial sphere. This includes also the forced oscillations of a sub-critically ($X < 1$) charged droplet in the AC quadrupole electric field of the Paul trap (see Fig. 1.5) as well as the slow oblate deformation of super-critically ($X > 1$) charged droplets, (see Fig. 1.8). Below, we detail the theoretical approaches on which our numerical code is based.

### 1.4.1 Navier-Stokes equation for Newtonian incompressible fluids

We define the following quantities that will be used throughout the manuscript. Let $\Omega$ be the volume of droplet bounded by the interface $\partial\Omega$, and $\vec{s}$ a vector pointing at the surface. Let $\vec{n}(\vec{s})$ be the unit vector normal to the surface pointing outwards, and $\vec{\tau}(\vec{s})$ a vector in the corresponding tangent plane. Let $\phi(\vec{s})$ be an arbitrary scalar function defined on the closed surface $\partial\Omega$. We introduce the following surface operators, (where the dependence on $\vec{s}$ has been omitted) [43]

* $\vec{\nabla}_s\phi = \vec{\nabla}\phi - \vec{n}(\vec{n} \cdot \vec{\nabla})\phi$ : tangential gradient (1.7)

* $\vec{\nabla}_s(\phi\vec{\tau}) = \vec{n} \cdot (\vec{\nabla} \wedge (\vec{n} \wedge (\phi\vec{\tau})))$ : surface divergence tangent vect. (1.8)

* $\vec{\nabla}_s(\phi\vec{n}) = -2H\phi$ : surface divergence of normal vector. (1.9)

The latter defines the local mean curvature $H$ of the surface at the surface point $\vec{s}$. Further, we define the following quantities characterizing the liquid droplet:

* $\vec{u}$ : velocity vector field [m/s] of the fluid.

* $\rho$ : density of the fluid [kg/m$^3$], assumed constant ($\vec{\nabla}\rho = 0$).

* $\epsilon_r$ : dielectric constant of the fluid ($\vec{\nabla}\epsilon_r = 0$).

* $\epsilon_0$ : permittivity of the vacuum.

* $\gamma$ : surface tension [N/m] of the fluid, assumed constant ($\vec{\nabla}_s\gamma = 0$).

* $\mu$ : dynamic viscosity [Pa.s] of the fluid.

* $R$ : radius [m] of the spherical droplet.
• $\vec{E}$ : electric field [V/m].

• $V$ : electric potential [V].

• $\sigma$ : surface charge density [C/m²].

• $\kappa$ : electric bulk conductivity [S/m] of the liquid.

• $\lambda$ : charge mobility [m²/s/V] of a charge carrier at the interface.

• $X$ : fissility of the droplet.

We want to simulate the shape-deformation of an initially spherical charged drop in a gas g atmosphere. The liquid is supposed electrically conductive. As a result, the charges accumulate at the droplet surface and the electric field is divergence-free inside and outside the droplet volume,

$$ \vec{\nabla} \cdot \vec{E} = 0 \quad \text{in } \mathbb{R}^3 \qquad (1.10) $$

This corresponds to the Taylor-Melcher leaky dielectric model [9]. In our particular case, the deformation is assumed axisymmetric and the liquid of the droplet incompressible. The density $\rho$ of the liquid is thus constant. The latter condition implies that the velocity vector field of the fluid is divergence-free (solenoidal),

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

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

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

where $T$ is the stress tensor, defined by

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

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

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

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 (1.15) $$

Note that $T$ is symmetric, $T = T^t$. As there are no volume charges, the electric force per unit volume is zero, $\nabla \cdot M = 0$. Thus, in the leaky dielectric model, the Maxwell tensor is divergence-free.
## Boundary conditions for viscous charged drops

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_s$,

$$T_{nn} = \vec{n} \cdot (T^g - T) \cdot \vec{n} = T_{nn}^g - T_{nn} = P_s / \partial\Omega \quad (1.16)$$

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

$$P_s = \gamma(\nabla \cdot \vec{n}) = \gamma(\nabla_s \cdot \vec{n}) = -2\gamma H \quad (1.17)$$

The right hand identity relates the mean curvature $H$ to the surface divergence of the normal vector $\vec{n}$ at a surface point $\vec{s}$, according to (1.8). Having defined the normal vector at the drop interface pointing outward, the mean curvature of a sphere of radius $R$ is negative, $H = -1/R$, so that $P_s > 0$. This definition ensures that the pressure of the liquid is higher than the pressure of the atmosphere! At room temperature, the dynamic viscosity of air is about 50 times lower than that of water. We suppose here 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$^6$), so that at the drop interface $\partial\Omega$ one has:

$$P_s = T_{nn}^g - T_{nn} \quad (1.18)$$

$$= -(p^g - p) + 2(\mu^g - \mu)D_{nn} + M_{nn}^g - M_{nn} \quad (1.19)$$

$$\simeq p - 2\mu D_{nn} + M_{nn}^g - M_{nn} \quad (1.20)$$

We introduced the Coulomb pressure $P_c$, defined as the jump of the Maxwell tensor across the interface

$$P_c = M_{nn}^g - M_{nn} = \frac{\varepsilon_0}{2} \left( (E_n^g)^2 - \varepsilon_r E_n^2 + (\varepsilon_r - 1) E_\tau^2 \right) \quad (1.21)$$

Injecting the surface charge density $\sigma = \varepsilon_0(E_n^g - \varepsilon_r E_n)$ into (1.21) allows expressing the Coulomb pressure as a function of $\sigma$ and the electric field $(E_n, E_\tau)$, 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)$$

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

$$-p = -P_s + P_c - \mu 2 D_{nn} \quad (1.22)$$

<footnote>$^6$As only the gradient of the pressure enters the dynamics, the pressure $p^g$ can be chosen arbitrarily.</footnote>
and the jump condition of the normal component of the stress tensor (1.16) may be expressed as

$$T_{nn} - M_{nn} = -P_s + P_c . \qquad (1.23)$$

In contrast, the shear stress $T_{\tau n}$ is continuous across the free interface,

$$T_{\tau n} \parallel = \vec{\tau} \cdot (T^g - T) \cdot \vec{n} = 0 \qquad (1.24)$$

$$= 2\mu^g D_{\tau n}^g + M_{\tau n}^g - (2\mu D_{\tau n} + M_{\tau n}) = 0 \quad /\partial\Omega . \qquad (1.25)$$

We deduce from (1.25), by neglecting the contribution of the viscous rate-of-strain tensor in the gas atmosphere, that the shear (tangent to the surface) component of the viscous rate-of-strain tensor in the liquid must be equal to the shear component of electric field $E_{\tau}$ at the surface times the surface charge density $\sigma$,

$$\mu 2 D_{\tau n} = (M_{\tau n}^g - M_{\tau n}) = \sigma E_{\tau} / \partial \Omega . \qquad (1.26)$$

Alternatively, boundary condition (1.25) may be expressed in the useful form,

$$T_{\tau n} - M_{\tau n} = \sigma E_{\tau} . \qquad (1.27)$$

where the latter equality is deduced from (1.15). Note that for a perfectly conducting liquid, $E_{\tau} = 0$, yielding the more usual condition that the shear strain rate must be zero at the liquid-gas interface. Equations (1.23) and (1.27) represent boundary conditions that must be fulfilled by a charged viscous drop in atmosphere. Finally the prohibition of mass transfer across the interface $\partial\Omega$ is guaranteed by the kinematic condition,

$$\vec{n} \cdot \left( \frac{d\vec{s}}{dt} - \vec{u} \right) = 0 \quad \text{on} \quad \partial\Omega . \qquad (1.28)$$

Equation (1.28) describes the time evolution of the shape deformation of the droplet$^7$!

$$\frac{\partial s}{\partial t} = u_r - \frac{\partial s}{\partial \theta} \frac{u_{\theta}}{s} .$$

<footnote>$^7$In our code, the shape function $s(\theta, t)$ is defined on the collocation points $\theta_i$, which are the zeros of a high order Legendre polynomial and thus given once for all. We propagate in time the shape function $s(\theta_i, t)$ for fixed angle $\theta_i$ using</footnote>
1.4.2 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 because of the non-zero normal component of the fluid velocity evaluated at the surface (1.28). The charge dynamics at the interface is governed by the surface charge transport (continuity) equation on a closed deforming surface [60, 61].

$$
\frac{D\sigma}{Dt} + u_n \vec{\nabla}_s \cdot \vec{n} + \vec{\nabla}_s \cdot \vec{j}_s = \kappa E_n \quad (1.29)
$$

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 term $u_n \vec{\nabla}_s \cdot \vec{n}$ accounts for the deforming surface. The term $\kappa E_n$ plays the role of a source term, that brings charges from the bulk to the surface or moves charges from the surface into the bulk. The surface charges are dragged by the liquid flow and pushed by the electric field. The surface current $\vec{j}_s$ has thus two contributions,

$$
\vec{j}_s = \sigma(u_\tau + \lambda E_\tau)\vec{\tau} \quad (1.30)
$$

The first is the velocity field of the fluid tangent to the surface $u_\tau$ and the second is the tangent component of the electric field $E_\tau$ times the charge mobility $\lambda$. We neglect the diffusion term which is usually 2 orders of magnitude smaller than the two former terms. Noting that $\vec{\nabla}_s \cdot \vec{n} = -2H\sigma$, (1.9), the surface charge evolution eventually reads

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

The surface divergence operator $\vec{\nabla}_s \cdot$ may not be well known to the reader, so we give a computationally useful expression of the second term in (1.31). Because the droplet deformation and charge dynamics have axis symmetry, $\partial_\phi \sigma = \partial_\phi u_\tau = \partial_\phi E_\tau = 0$, the surface divergence simplifies, using 1.8, to

$$
\vec{\nabla}_s \cdot (\vec{\tau} \sigma (u_\tau + \lambda E_\tau)) = \vec{n} \cdot (\vec{\nabla} \wedge \sigma (u_\tau + \lambda E_\tau) \vec{e}_\phi) \quad (1.32)
$$

In the general case of curvilinear orthogonal coordinates (ξ, η, φ) and for an axis symmetric function F = σ(uτ + λEτ) defined on the surface (ξs(η), η), the expression (1.32) yields,

$$
\vec{n} \cdot \left( \vec{\nabla} \wedge (F(\xi(\eta), \eta) \vec{e}_{\phi}) \right) \Big|_{\xi=\xi_s} = \vec{n} \cdot \vec{e}_{\xi} \frac{1}{h_{\eta}h_{\phi}} \frac{d}{d\eta} (h_{\phi} F(\xi_s(\eta), \eta)) \quad (1.33)
$$
The latter expression can be used for the computation of the surface divergence of an axis symmetric tangent vector defined on the surface.

In spherical coordinates $(r, \theta, \phi)$, with $h_r = 1$, $h_\theta = r$ and $h_\phi = r \sin\theta$, one has at the surface $(R(\theta), \theta)$

$$ \vec{n} \cdot (\vec{\nabla} \wedge (F[r, \theta]e_{\phi})) = \vec{n} \cdot \vec{e}_r \frac{1}{h_\theta h_\phi} \frac{d}{d\theta} (h_\phi F[R(\theta), \theta]) \quad (1.34) $$

$$ = \frac{1}{\sqrt{1 + k^2} R^2(\theta) \sin\theta} \frac{1}{d\theta} (\sin\theta R(\theta)F[R(\theta), \theta]) \quad (1.35) $$

with $k = \frac{1}{R(\theta)} \frac{dR}{d\theta}$.

### 1.4.3 Kinetic energy dissipation

It is interesting to look at the dissipation of the kinetic energy of a viscous fluid. Multiplying (1.12) by $\vec{u}$ and integrating over the droplet volume $\Omega$ yields the time evolution of the mechanical energy of the droplet. Using the identity,

$$ \vec{\nabla} \cdot (\vec{u} \cdot T) = \vec{u}(\vec{\nabla} \cdot T) + T : \frac{1}{2} (\vec{\nabla} \vec{u} + (\vec{\nabla} \vec{u})^t) \quad , \quad (1.36) $$

which holds for a symmetric tensor $T$, the dissipation of the kinetic energy is given by

$$ \frac{d}{dt} E_{\text{kin}} = \frac{d}{dt} \int_{\Omega} \frac{\rho u^2}{2} d\Omega \quad (1.37) $$

$$ = \int_{\Omega} \vec{u}(\vec{\nabla} \cdot T) d\Omega \quad (1.38) $$

$$ = \underbrace{\int_{\partial\Omega} \vec{u} \cdot (T - M) \cdot \vec{n} ds^2}_{\mathcal{T}} - 2\mu \underbrace{\int_{\Omega} D : D d\Omega}_{\mathcal{D}} \quad (1.39) $$

where $\mathcal{T}$ stands for the power of traction (driving power) and $\mathcal{D}$ for the dissipation power. The operator $:$ stands for the double dot product of tensors<fnref n="8" />. We used here the property that, in the leaky dielectric model, the symmetric Maxwell tensor $M$ is divergence-free so that from the tensor identity 1.36 and divergence theorem one has

$$ \int_{\Omega} D : M d\Omega = \int_{\Omega} \vec{\nabla} \cdot (\vec{u} \cdot M) d\Omega = \int_{\partial\Omega} \vec{u} \cdot M \cdot \vec{n} ds^2 \quad (1.40) $$

<footnote n="8">⁸Let $h_{ij}$ and $g_{ij}$ be two square matrices, then $h_{ij} : g_{ij} = \sum_{ij} m_{ij} g_{ij}$ is a scalar</footnote>
We also note that for the identity matrix $I$ times the pressure $p$ one has

$$D : (-pI) = p \vec{\nabla} \cdot \vec{u} = 0 \qquad (1.41)$$

Injecting the boundary conditions (1.23) and (1.27) into (1.39) yields

$$\frac{d}{dt} E_{\text{kin}} = \underbrace{\int_{\partial \Omega} (u_n(-P_s + P_c) + u_{\tau} \sigma E_{\tau}) ds^2}_{\mathcal{T}} - 2\mu \underbrace{\int_{\Omega} D : D d\Omega}_{\mathcal{D}} \qquad (1.42)$$

The dissipation power $\mathcal{D}$ is proportional to the viscosity and involves only the viscous rate of strain tensor $D$. For rotational fluids the strain rate tensor is **not** divergence-free and dissipative forces exist inside the droplet as well as at the interface<fnref n="9" />. The traction power includes the capillary pressure, the Coulomb pressure but also the shear stress of the Maxwell tensor, all evaluated at the interface. Equation (1.42) is valid for all Newtonian viscous fluids satisfying the boundary condition ((1.23) and (1.27) at a free charged surface.

### 1.4.4 Potential Flow approximation: Case of a perfectly conducting inviscid liquid drop

Even by assuming an axisymmetric shape deformation, the numerical propagation of the fluid via the Navier-Stokes (NS) equation with moving boundaries is quite CPU expensive. Several authors [50, 57, 58, 59] used the finite element method (FEM) to solve the NS equation and propagate the droplet interface in time. While they successfully describe the drop deformation, the high CPU cost makes systematic studies challenging. We simplify the description of the fluid dynamics by assuming that the liquid is inviscid (zero viscosity) and the flow irrotational $\vec{\nabla} \wedge \vec{u} = 0$. The latter implies that the velocity field $\vec{u}$ of the fluid can be obtained as the gradient of the scalar velocity potential $\Psi$

$$\vec{u} = \vec{\nabla} \Psi \quad \text{in } \Omega, \qquad (1.43)$$

leading thus to the potential flow (PF) approximation<fnref n="10" />. Within these assumptions, the axis-symmetric Navier-Stokes equation simplifies to the 1D non-linear

<footnote n="9">⁹For irrotational fluids the strain rate tensor is divergence-free so that the last term in (1.42) involving a volume integration can be replaced by a surface integration.</footnote>

<footnote n="10">¹⁰In flow regions where vorticity is known to be important, such as wakes and boundary layers, potential flow theory is not able to provide reasonable predictions of the flow. Fortunately, there are often large regions of a flow where the assumption of irrotationality is valid which is why potential flow is used for various applications. For instance in: flow around aircraft, groundwater flow, acoustics, water waves and macro-sized (cm range) water droplets.</footnote>
time-dependent Bernoulli equation, giving the time-evolution of the velocity potential $\Psi$ on the interface as a function of the interface pressure $p$.

$$ \frac{\partial \Psi}{\partial t} + \frac{\vec{u}^2}{2} + \frac{p}{\rho} = 0 \quad \text{on} \quad \partial \Omega \qquad (1.44) $$

Note that, unlike Betelú et al. [46] who studied the shape deformation using the Stokes equation, the inviscid liquid model ignores viscous properties of the fluid but includes the convection term $\rho \frac{\vec{u}^2}{2}$, accounting for time-independent acceleration of the flow with respect to space and thus accounting for some inertial effects.

In the limit of an inviscid fluid, $\eta = 0$, the boundary condition (1.27) can only be fulfilled for a perfect conductor, $-\nabla_s V = E_\tau = 0$. The liquid is thus considered to be perfectly conducting so that the electric field is zero in the drop volume and the electric potential at the surface and equipotential. We obtain from (1.23) and (1.27) that the pressure at the liquid interface given by the difference of the surface tension pressure and the Coulomb pressure,

$$ p = P_s - P_c = -2\gamma H - \frac{\sigma^2}{2\varepsilon_0} \quad \text{on } \partial\Omega \qquad (1.45) $$

while the jump condition of the tangential component of the stress tensor (1.27) is automatically satisfied. As the velocity field is divergence-free, the velocity potential $\Phi$ satisfies the Laplace equation

$$ \vec{\nabla}^2 \Psi = 0 \quad \text{in} \quad \Omega, \qquad (1.46) $$

The boundary condition necessary to solve the above Laplace equation is given by (1.44). The velocity field of the fluid at the droplet interface is finally computed using (1.43) and the interface propagated using (1.28),

## Surface charge distribution for a perfectly conducting liquid

The liquid being a perfect conductor, the charge distribution is always at equilibrium so that the electric potential is constant all over the surface. We use thus the latter condition to compute the surface charge density for a given droplet shape. The electric field inside a conductor is zero, it is sufficient to compute the electric field outside the droplet volume. As the electric field is divergence-free, the electric potential $V$ satisfies the Laplace equation, with the boundary condition that $V$ is constant at the interface,

$$ \vec{\nabla}^2 V = 0 \quad \text{in} \quad \mathbb{R}^3 \setminus \Omega \qquad (1.47) $$

$$ V = V_0 \quad \text{on} \quad \partial \Omega, \qquad (1.48) $$
The potential V can now be readily computed for an arbitrary equipotential V₀ using, for example, the technique where V is projected onto harmonic functions, see appendix A.1. The surface charge density σ = −ε₀∂V/∂n is deduced from the normal derivative of the electric potential V, with V₀ being determined by the condition that the integration of the surface charge density over the droplet surface equals the total charge Q = 8π√Xγε₀R³ of the droplet,

$$
Q = \int_{\partial\Omega} \sigma(t, \vec{s}) \, ds^2 . \qquad (1.49)
$$

**Dimensionless equations:**

In order to find the quantities that control the dynamics we put equation 1.44 into
dimensionless form. We introduce the characteristic deformation time scale

$$
\tau_0 = \sqrt{\frac{\rho R^3}{\gamma}} \qquad (1.50)
$$

and define the following dimensionless quantities

$$
\tilde{t} = t/\tau_0 \tag{1.51}
$$

$$
\tilde{\Psi} = \Psi \frac{\tau_0}{R^2} \qquad (1.52)
$$

$$
\sqrt{X} \tilde{\sigma} = \sigma \left( \frac{4\gamma\epsilon_0}{R} \right)^{-1/2} \qquad (1.53)
$$

$$
\tilde{H} = H R \tag{1.54}
$$

$$
\vec{u} = \vec{u} \frac{\tau_0}{R} \qquad (1.55)
$$

The dimensionless form of (1.44) redas then,

$$
\frac{\partial \tilde{\Psi}}{\partial \tilde{t}} + \frac{\tilde{u}^2}{2} - 2 \tilde{H} - 2X \tilde{\sigma}^2 = 0 \quad \text{on} \quad \partial \Omega \qquad (1.56)
$$

We remind that ˜H < 0 is a negative quantity. As we consider here a perfectly
conducting, inviscid liquid, the fissility parameter X is the only free parameter
of the model that influences the pathway. The quantities R, γ, ρ define merely
the time scale τ₀ of the droplet deformation. The latter is therefore quite general
within the given hypotheses. Details about the numerical resolution of the PF
model are given in appendix A.1.
<figure><img src="image_25.png" /><figcaption>Figure 1.21: Bottom panel:Calculated shape sequence of an initially critically charged droplet $X = 1$, experiencing charge emission $X = 0.5$ at the moment a bulb appears at the tip. Physical times are obtained by multiplying dimensionless times t, shown under each shape, by the characteristic time $\tau$ of a given droplet. Fluid is nonviscous and perfectly con- ducting. Charge emission is chosen as origin of time. The red lemon-fit has been slightly shifted for better comparison. Top panel: Snapshots of a deforming glycol droplet for comparison.</figcaption></figure>

## Simulated results using the PF approach

In Fig. 1.21, we show snapshots of the simulation before ($X = 1$) and after charge emission ($X = 0.5$). The starting point is a critically charged prolate spheroid of small eccentricity (almost spherical drop)

$$e = \sqrt{1 - a^2/b^2} = 0.2 \ ,$$

where a and b are defined in section 1.2.1. We see that the model is able to reproduce qualitatively the whole initial sequence of the non-linear deformation
dynamics of a critically charged droplet. In particular, for $X = 1$, we end up with the same spindle-like form as the glycol droplet. At the moment, when the droplet starts to form a jet at the tips (snapshot labeled by the dimensionless time 0), we reduce the charge by -30% by changing suddenly the fissility to $X = 0.5$. This allows us to emulate the charge emission and to describe the shape deformation after it, without stopping the calculation. After emission, the droplet retracts similarly to experimental data, adopting a capsule-like form, which, unlike the glycol droplet has a bulge in the middle. The model also shows some weaknesses. The lack of viscosity hinders the dissipation of the excess of potential energy obtained during charge emission. The model is therefore not able to describe the deformation back to equilibrium shape (sphere). Further, as viscosity slows down the flow, we expect that the deformation time we get from the simulation underestimates the experimentally observed one. Indeed, glycol droplets such as the one shown in Fig. 1.21 with $(\rho, \gamma, R) = (1.1 \text{ g/cm}^3, 47 \text{ mN/m}, 24 \text{ µm})$, (see Table 1.1), have a characteristic time scale of $\tau_0 = 18 \text{ µs}$. From Fig.1.21 we note that the glycol droplet needs $22 \text{ µs}$ to attain the spindle-like shape, starting from a prolate spheroid with a ratio $a/b = 2$. For the same deformation, the simulation needs $0.6 \times \tau_0 \approx 11 \text{ µs}$. Compared to the observed $22 \text{ µs}$, the inviscid drop deforms about 2 times faster.

For the supercooled water droplets shown in Fig. 1.10, the droplet needs $30 \text{ µs}$ to deform from the prolate spheroid with a ratio $a/b = 1.6$ ($15 \text{ µs}$) to the lemon-shaped drop ($45 \text{ µs}$). Our water drop has a characteristic time of $\tau_0 = 24 \text{ µs}$ (using $(\rho, \gamma, R) = (1.0 \text{ g/cm}^3, 76 \text{ mN/m}, 35 \text{ µm}$) The simulation gives for the some deformation $0.9 \times \tau_0 \approx 22 \text{ µm}$, meaning the inviscid liquid deforms 1.4 times faster than the supercooled water drop. We assign the trend to the lower viscosity of supercooled water with respect to glycol.

The simulations describe the shape of the tips of the droplet just before surface breakup. It was an open question whether the jets are emitted from a Taylor cone [8] as assumed by de la Mora [48], likely to electrosprays [23], or rather from a hyperboloidal tip as was pointed out by Yarin et al. [37] in the case of sessile and pendant droplets in an electric field. When zooming on the tip, we note that the jets are not emitted from a conical tip, but at best from a hyperboloid which is tangent to a cone of 33°. Yarin et al. showed that the maximal electric potential at which a stationary infinite hyperboloid can exist, corresponds to a hyperboloid tangent to a cone of 33.5° [37]. This critical hyperboloid was the starting point from which jets are formed from pendant droplets for increasing electric fields. Interestingly, it is close to the hyperboloid that best fits the tip, hinting that the last one is close to break-up. Nevertheless, in our case, the hyperboloid is not a satisfying fit either, in the sense that it only fits the tip, instead of the whole droplet. We looked for a geometrical description of the whole shape, just before
charge emission, and found that the droplet is remarkably well fitted by a "lemon" of equation

$$z^2 = A^2 - (B + r)^2, \quad A = 3.02, \quad B = 2.25 \qquad (1.57)$$

with $r = \sqrt{x^2 + y^2}$ being the distance from the symmetry axis. Note that the fit has only one free parameter as A and B are related by the volume constraint. The lemon is shown in the upper inset of Fig. 1.22 as a red contour. When zooming on the tip, we see only a small discrepancy close to the apex (lower inset of Fig. 1.22). We deduce that the tip is tangent to a cone of 39°, thus significantly narrower than the Taylor cone (49.3°) but larger than the one obtained by Betelú et al. (cone of 25° for a viscous droplet [46]). The experimentally observed tip’s opening half-angle of 30° lies between both predictions. The lemon may be considered as the asymptotic shape of a critically charged droplet, if breaking up could be hindered.

For further analysis of the tip behavior, we introduce the local fissility $x(\theta)$ defined as the Coulomb pressure divided by the pressure due to the surface tension, $x(\theta) = X\sigma^2(\theta)/H(\theta)$. In order to be able to break up the surface, the Coulomb pressure must outbalance the stress due to the surface tension $x(\theta) > x_c \ge 1$, where the critical local fissility $x_c$ still needs to be determined. For a spheroid, the Coulomb pressure can be expressed as a function of the Gaussian curvature $K$ of the surface [25], yielding

$$x(\theta) = X\sqrt{K(\theta)}/H(\theta) .$$

As for spheroids, one has always $H^2 \ge K$, independently of the eccentricity. Hence, the local fissility of critically charged spheroids ($X = 1$) cannot exceed the unit $x(\theta) \le 1$ and jets cannot appear at their ends, even for large eccentricities. Critically charged lemons, however, behave differently as the local fissility is no longer bound and can exceed one at the tips. Indeed, for the lemon tip shape depicted in Fig. 1.22, the surface charge density diverges close to the apex as $\sigma(r \to 0) \propto r^{-0.59}$, while the mean curvature diverges as $H(r \to 0) \propto r^{-1}$ (not shown). Hence, the local fissility diverges at the apex as $x(r \to 0) \propto r^{-0.18}$. Such critically charged lemon shapes lead therefore inevitably towards surface break up of the tip. In Fig.1.22, we give the local fissility (upper panel) and charge distribution (lower panel) of the droplet at different times (shapes) before charge emission and compare them to the one calculated for the lemon shape. As soon as the shape deviates from a spheroid to become rugbyball-like (around $t_1 = -0.43$), the local fissility exceeds 1 at the tips $x(0) = 1.08$. This shows that the local fissility can reach values larger than 1 at nonsingular surface points. However, this also shows that $x > x_c = 1$ is not a sufficient criterion for surface break-up but indicates rather the less stable region. At time $t_0 = -0.05$, the charge distribution compares well to the charge distribution on the lemon surface except that it stays bounded at the apex. At that moment, the electric field at the tip is about 6 times larger
<figure><img src="image_26.png" /><figcaption>Figure 1.22: Local fissility and tip formation</figcaption></figure>

than initially. Finally, a bulb appears at the tip, which indicates the formation of a jet. The simulated sequence is quite general for inviscid perfectly conducting liquid droplets. However, properties of the fluid like viscosity and charge mobility may change the details of the deformation dynamics. For polar liquids like water and glycol, the Coulomb pressure on the droplet surface may significantly deviated such as the one of a perfectly conducting droplet.

We compare in Fig. 1.23 the simulated pathway to the experimental data. The simulation should well describe the deformation of a critically charged mercury droplet, which has low viscosity and is an excellent conductor. It predicts that the pointed ends are formed for a ratio *a/b* = 2.6, which is significantly lower than the values obtained for supercooled water (*a/b* = 3.2) and glycol (*a/b* = 3.5).
<figure><img src="image_27.png" /><figcaption>Figure 1.23: Pathway (red curve) of an inviscid perfectly conducting, compared to experimental data of a $\emptyset = 48\mu m$ glycol and supercooled $\emptyset = 70\mu m$ water droplet.</figcaption></figure>

We expect that the difference is due to the viscosity and finite charge mobility of the polar liquids. Indeed, for even lower viscous drops, namely negatively and positively charged water droplets at room temperature, a lower ratio *a/b* between 2.7 and 2.8 was found (see Fig. 1.32). After charge emission, the deformation is mainly driven by the surface tension. The latter tends first to flatten the tips to form a capsule. The slope of the simulated curve after charge emission compares well to the experimental low-viscous water curve. As mentioned before, the model cannot describe the way back to the equilibrium shape and we stop the simulation when *n* = 3.2. However, the comparison with the experimental curve for glycol fails, highlighting the role of viscosity and charge mobility. We deduce that after charge emission, the details of the pathway are dominated by viscous properties of the fluid and, thus, poorly described by the present model. Nevertheless, the inviscid fluid simulations give the asymptotic ratio *a/b* and slope of the curve, after charge emission, that cannot be overcome by low-viscous fluids. The influence of the viscosity and charge mobility of real fluids on the deformation pathway will be studied in the following sections.
### 1.4.5 Viscous Potential Flow (VPF) approximation

The inviscid liquid model suffers clearly from the absence of a dissipation term in the fluid dynamics. While keeping the potential flow approximation, which has a clear CPU cost advantage over to the full Navier-Stokes equation, we relax the condition that the drop liquid is inviscid. This yields the viscous potential flow model (VPF) and the pressure at the liquid interface reads (1.22)

$$p = P_s - P_c + 2\mu D_{nn} \partial\Omega . \qquad (1.58)$$

Compared to the inviscid liquid model, the pressure has now an additional term which is proportional to the viscosity $\mu$. The strain rate tensor $D$ at the interface is computed using the irrotational velocity field $\vec{u}$. Note that, while for an irrotational flow the viscous force per unit volume is zero $\nabla \cdot D = 0$ inside the drop, its contribution to the surface pressure is non-zero

$$D_{nn} \neq 0 \quad \text{and} \quad D_{\tau n} \neq 0 \quad \text{on} \quad \partial\Omega \qquad (1.59)$$

As a result, the jump of the shear (tangential) 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} \quad \partial\Omega . \qquad (1.60)$$

Compared to (1.27) the last equality is not satisfied in VPF. Indeed, VPF has not enough degrees of freedom to enforce that the shear strain rate $D_{\tau n}$ equals the shear stress of the Maxwell tensor. 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 (1.59). The rotational part of the velocity field is missing here to fulfill the condition (1.60). Compared to a rotational flow where condition (1.60) can be enforces, VPF results in an incorrect dispersion relationship, i.e. the damping of a given surface mode is underestimated for low viscosities and overestimated for very large viscosities [51].

While the VPF model is certainly not perfect, it includes after all viscous forces at the free boundary and thus dissipation of the kinetic energy $E_k$. The power equation (1.42) is valid for rotational as well as for irrotational fluids. Using (1.58) and (1.60) and reminding that in the VPF approach the strain rate tensor $D$ is also divergence-free, equation (1.42) leads to

$$\frac{d}{dt} E_{\text{kin}} = \underbrace{\int_{\partial\Omega} u_n(-P_s + P_c)ds^2}_{\mathcal{T}^{\text{VPF}}} - 2\mu \underbrace{\int_{\partial\Omega} u_n D_{nn} ds^2}_{\mathcal{D}^{\text{VPF}}} \qquad (1.61)$$

Compared to the inviscid liquid model, VCP adds the dissipation term $\mathcal{D}^{\text{VPF}}$ proportional to the viscosity $\mu$. The added dissipation allows the fluid to relax to its equilibrium shape and thus to model the shape dynamics of the droplet after jet emission.
### 1.4.6 Viscous Corrected Viscous Potential Flow (VCVPF) approach

By ignoring the zero shear stress boundary condition at the free surface in VPF, we ignored also dissipative forces originating at the surface due to the zero shear stress boundary condition. In order to compensate for the lack of those viscous forces, Joseph and Wang [40, 45] derived a viscous correction formulation for the irrotational pressure. It expresses the idea that a viscous correction of the pressure at the gas-liquid interface is required to compensate for the non-vanishing irrotational shear stress. The accepted 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 (VC)VPF. They proposed thus to add a fictive viscous pressure $p_μ$ to the normal stress,

$$T_{nn} - M_{nn} = -P_s + P_c - p_μ \quad (1.62)$$

The viscous pressure $p_μ$ is assumed to satisfies the Laplace equation in the drop volume<fnref n="11" /> and is defined such as the power of the viscous pressure equals the dissipation power of the unbalanced shear stress at the boundary,

$$\int_{\partial\Omega} \vec{n} \cdot \vec{u}(-p_μ)ds^2 = \int_{\partial\Omega} \vec{\tau} \cdot \vec{u} \underbrace{(-2μD_{τn} + σE_τ)}_{T_{τn}^g - T_{τn}} ds^2 \quad (1.63)$$

Expression (1.63) however does not uniquely define $p_μ$. We will further ask that the power of the viscous pressure balances the power of the jump of the shear stress across the boundary for each surface element $ds^2$,

$$u_n(-p_μ)ds^2 = u_τ(-2μD_{τn} + σE_τ)ds^2 \quad (1.64)$$

Loosely said, the elementary force $p_μ ds^2$ accounts for the missing viscous forces at the free surface element that would originate if the zero shear stress boundary conditions was satisfied. A numerical method for the evaluation of $p_μ$ is given in Appendix A.4. The power of traction in VCVPF becomes using (1.62) and (1.64),

$$\begin{aligned} \mathcal{T}^{\text{VCVPF}} &= \int_{\partial\Omega} \vec{u} \cdot (T - M) \cdot \vec{n} \, ds^2 \\ &= \int_{\partial\Omega} u_n(-P_s + P_c - p_μ) \, ds^2 + 2μ \int_{\partial\Omega} u_τ D_{τn} \, ds^2 \end{aligned} \quad (1.65)$$

$$= \int_{\partial\Omega} (u_n(-P_s + P_c) + u_\tau \sigma E_\tau) ds^2 \approx \mathcal{T}^{\text{VPF}} \quad (1.66)$$

<footnote n="11">¹¹ensuring that the correction is divergence-free</footnote>
which is the same expression than for rotational fluids (1.42), apart that here $\vec{u}$ is purely irrotational. The power of traction $\mathcal{T}^{\text{VCVPF}}$ differs from the one in VPF by the term $\int u_{\tau} \sigma E_{\tau} ds$, which is usually negligible for conducting liquids. So the traction power in VPF and VCVPF can be considered similar. Reminding that the strain rate tensor $D$ is divergence-free, the dissipation power $\mathcal{D}$ reads

$$ \mathcal{D}^{\text{VCVPF}} = 2\mu \int_{\partial\Omega} (u_n D_{nn} + u_{\tau} D_{\tau n}) ds^2 \neq D^{\text{VPF}} \quad (1.67) $$

The dissipation power $\mathcal{D}^{\text{VCVPF}}$ has a contribution from the shear strain rate $\int u_{\tau} D_{\tau n} ds$, not found in VPF (1.61). However, while VCVPF allows to recover the correct expression of the power of traction, it only approaches the power of dissipation of the Navier-Stokes equation, mainly because the rotational part of the velocity field is still missing. Besides, the rotational part of the velocity field is also necessary to satisfy explicitly the jump condition of the shear stress at the free surface.

## Dimensionless time evolution of the velocity potential $\Psi$ and surface charge $\sigma$.

Following the lines of section 1.4.4, the time evolution of the dimensionless velocity potential $\tilde{\Psi}$ at the interface is given in VCVPF

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

where we introduced the dimensionless Ohnesorge number,

$$ \text{Oh} = \frac{\mu}{\sqrt{\rho\gamma R}} \left( \frac{\text{viscosity}}{\sqrt{\text{inertia} \times \text{surface-tension}}} \right) \quad (1.69) $$

that relates the viscous forces to inertial and surface tension forces. Larger Ohnesorge numbers indicate a greater influence of the viscosity. If we omit the viscous pressure term, we recover indeed the time evolution of $\Psi$ in the VPF model, where only the normal rate stress contributed to the dissipation. In VCVPF, the viscous pressure $p_{\mu}$ adds (viscous) forces at the surface that balances the non-zero power of the shear stress at the interface. If, because of finite conductivity, the tangential component of the electric field is non-zero, $p_{\mu}$ also adds 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 appendix A.4.

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

$$ \tau_s^{-1} = \lambda \sqrt{\frac{\gamma}{R^3 \epsilon_0}} \quad \tau_b^{-1} = \frac{\kappa}{\epsilon_0 \epsilon_r} \quad (1.70) $$
allows to define the dimensionless charge carrier mobility $\tilde{\lambda}$ and dimensionless bulk conductivity $\tilde{\kappa}$

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

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

where the characteristic deformation time $\tau_0$ is given by (1.50). Note that $\tilde{\lambda}$ depends on $X$. Further, defining the characteristic surface charge density

$$ \sigma_0 = \sqrt{4\gamma\varepsilon_0/R} \ , $$

introducing the dimensionless quantities

$$ \sigma = \sqrt{X\tilde{\sigma}} \sigma_0 \qquad (1.73) $$

$$ E_{\tau} = \sqrt{X \tilde{E}_{\tau} \sigma_{0}/\varepsilon_{0}} \qquad (1.74) $$

$$ E_n = \sqrt{X \tilde{E}_n \sigma_0 / \varepsilon_0} \qquad (1.75) $$

and using the definitions (1.51,1.54,1.55), allow us to put the surface charge transport with moving boundaries (1.31) into a dimensionless form,

$$ \frac{D\tilde{\sigma}}{D\tilde{t}} + \vec{\nabla}_{\tilde{s}} \cdot (\vec{\tau}\tilde{\sigma}(\tilde{u}_{\tau} + \tilde{\lambda}\tilde{E}_{\tau})) - 2\tilde{H}\tilde{\sigma}\tilde{u}_{n} = \tilde{\kappa}\varepsilon_{r}\tilde{E}_{n} \quad \text{on } \partial\Omega \qquad (1.76) $$

where all the quantities are now dimensionless. Note that the factor $\varepsilon_r$, found in the last term, simply comes from the factor $\varepsilon_r^{-1}$ used in the definition<fnref n="12" /> of the bulk relaxation rate (1.70). The charge dynamics at the interface is controlled by two parameters namely $\tilde{\lambda}$ and $\tilde{\kappa}$. If one of the latter has a value much larger than unit, 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 the perfect conductor. For $\tilde{\lambda} \simeq 1$ and $\tilde{\kappa} \simeq 1$, the charge at the surface may not be at equilibrium during the ultra-fast deformation (in the sense that the potential at the surface is no longer an equipotential), because the charges at the surface are now also driven by the flow of the liquid at the surface. As a result, the electric field inside the

<footnote n="12">¹²The inner field $E_n$, obtained via (A.26), includes the polarization charges that appear at the interface separating two dielectric media (here the liquid of dielectric constant $\varepsilon_r$ and vacuum). As a result, the inner field is scaled by the factor $\varepsilon_r^{-1}$. That is why the bulk relaxation rate also scales like $\varepsilon_r^{-1}$. We show this explicitly, by using the factor $\varepsilon_r^{-1}$ in the definition of the bulk relaxation rate.</footnote>
drop is non-zero and the dimensionless Coulomb pressure $\tilde{P}_s$ found in (1.68) and defined by

$$ \tilde{P}_c = P_c \left( \frac{\sigma_0^2}{2\epsilon_0} \right)^{-1} \qquad (1.77) $$

$$ \tilde{P}_c = \left( \tilde{\sigma} + \epsilon_r \tilde{E}_n \right)^2 - \epsilon_r \tilde{E}_n^2 + (\epsilon_r - 1) \tilde{E}_\tau \qquad (1.78) $$

deviates from the one of the perfect conductor, $\tilde{P}_c = \tilde{\sigma}^2$. Hence, the shape deformation of the drop depends now also sensibly on the parameters $\tilde{\lambda}$ and $\tilde{\kappa}$. Note that, while $\tilde{\lambda}$ does not depend on the radius $R$ of the drop, $\tilde{\lambda}$ scales linearly with $R$. For a given bulk conductivity $\kappa$ of the liquid, larger droplets have thus a larger dimensionless bulk conductivity.

Decay rate: VPF vs VCVPF

<figure><img src="image_28.png" /><figcaption>Figure 1.24: Dimensionless decay rate $\delta$ of the quadrupole amplitude as a function of the dimensionless viscosity $\text{Oh}$ for the VPF and VCVPF models and for the exact calculation. The values of $\text{Oh}$ for which the decay is critical are shown by thin vertical lines and labeled accordingly, where $\text{Oh}_c^{\text{VCVPF}} \simeq 0.56$, $\text{Oh}_c^{\text{VPF}} \simeq 1.41$ and $\text{Oh}_c^{\text{Exact}} \simeq 0.76$.</figcaption></figure>

In the following section we will discuss briefly the differences between VPF and VCVPF regrading the dissipation. Using the small amplitude oscillation approximation, the dispersion relationship of a viscous (uncharged) drop was analysed by Padrino et al. [51] in the case of PF, VPF and VCVPF. The dispersion relationship links the complex frequency to the surface modes (quadrupole, hexapole,
octupole...) of the surface. The dimensionless complex oscillation frequency $\tau_0\Omega$ for the different approaches becomes

$$ \tau_0\Omega_m^{\text{PF}} = \pm i\sqrt{(m+2)m(m-1)} = \pm i(\tau_0\omega_m) \quad (1.79) $$

$$ \tau_0\Omega_m^{\text{VPF}} = m(m-1)\text{Oh} \pm i\sqrt{(\tau_0\omega_m)^2 - (m(m-1)\text{Oh})^2} $$

$$ \tau_0\Omega_m^{\text{VCVPF}} = (2m+1)(m-1)\text{Oh} \pm i\sqrt{(\tau_0\omega_m)^2 - ((2m+1)(m-1)\text{Oh})^2} $$

where $m$ stands for the $m^{th}$ deformation mode of the droplet surface. We note that for the quadrupole mode, the dimensionless decay rate $\delta_2 = \text{Real}(\tau\Omega_2)$ in VCVPF is 5/2 time the one in VPF, making thus VCVPF noticeably different from VPF. For larger $m \gg 1$ the ratio tends toward 2. The decay rate of the quadrupole moment of the surface oscillations as a function of Oh are shown in figure 1.24 and compared to the exact calculation form [51]. Remarkably, it was found that in the low viscosities range, the VCVPF dispersion relationship tends to the exact one as was noted in [51]. Indeed as shown in Fig. 1.24, for $\text{Oh}< 0.1$, the VCVPF model gives a decay rate $\delta_{m=2}$ much closer to the exact one than the VPF model.

We label $\text{Oh}_c$ the critical value of the dimensionless viscosity, for which the imaginary part of the surface mode becomes zero. It is different for VPF and VCVPF:

$$ \text{Oh}_c^{\text{VPF}} = \sqrt{\frac{(m+2)}{m(m-1)}} \quad (1.80) $$

$$ \text{Oh}_c^{\text{VCVPF}} = \sqrt{\frac{(m+2)m}{(2m+1)^2(m-1)}} = \frac{m}{2m+1}\text{Oh}_c^{\text{VPF}} \quad (1.81) $$

In VCVPF, the critical value is shifted towards lower viscosities, compared to VPF. For the quadrupole moment $m = 2$, we get $\text{Oh}_c^{\text{VPF}} \simeq 1.41$ and $\text{Oh}_c^{\text{VCVPF}} \simeq 0.56$, where the latter is much closer to $\text{Oh}_c^{\text{Exact}} \simeq 0.76$ of the exact calculation. Now, our glycol micro-droplets have typically a ($\text{Oh} \simeq 0.4$) and are thus close to critically damping in VCVPF but not in VPF. So glycol droplets are expected to give sensibly different results in the VPF or VCVPF approach. While for low viscous water droplets, ($\text{Oh} \simeq 0.025$), VCVPF and VPF should give similar results, VCVPF is expected to come closer to the exact solutions. The VCVPF model seems thus preferable and we used for our simulations mainly the VCVPF model, keeping in mind that VCVPF calculations should be limited to $\text{Oh} < \text{Oh}_c^{\text{VCVPF}} \simeq 0.56$.

## Charge break up

The ultra fast jet formation and emission are not simulated in this work as the underlying numerical method, based on harmonic function expansion in prolate
coordinate system, is not able to solve the Laplace equation for shapes with almost singular tips<fnref n="13" />. Successful simulations of the jet formation based on FEM approaches and are given in [57, 61, 59]. Now, the jet not only removes charges from the drop's surface but removes also fast flowing matter. As a result, the jet reduces in several µs the charge and the kinetic energy of the droplet. We use two different methods to remove the charge and kinetic energy of the droplet.

1. In the first and simplest method, we simulate the jet emission by reducing suddenly the total charge and kinetic energy of the droplet, by applying a constant factor to the surface charge density $\sigma$ and to the velocity potential $\Psi$.

$$ \sigma(\vec{s}) \rightarrow \sqrt{1 - \frac{\Delta Q}{Q}} \sigma(\vec{s}) \qquad (1.82) $$

$$ \Psi(\vec{s}) \rightarrow \sqrt{1 - \frac{\Delta E_{\text{kin}}}{E_{\text{kin}}}} \Psi(\vec{s}) \qquad (1.83) $$

where $Q$ and $E_{\text{kin}}$ are the charge and the kinetic energy of the droplet just before charge break up and $\Delta Q$ the emitted charge and $\Delta E_{\text{kin}}$ the loss of kinetic energy due to jet emission<fnref n="14" />. The deformation dynamics of the subcritically charged drop $X < 1$ continues then using the same propagator but with reduced charge and kinetic energy. The pieces of information about the amount of emitted charge and lost kinetic energy are usually taken from experimental data.

2. In the second method, we simulate the charge emission by adding for a limited (dimensionless) time $\Delta t$ a source term $\tilde{\sigma}s$ to equation (1.76) that simulates the charge loss at the tip at the rate $\tilde{\tau}_{\text{jet}}$,

$$ \frac{D\tilde{\sigma}}{D\tilde{t}} + \vec{\nabla}_{\tilde{s}} \cdot (\vec{\tau}\tilde{\sigma}(\tilde{u}_{\tau} + \tilde{\lambda}\tilde{E}_{\tau})) - 2\tilde{H}\tilde{\sigma}\tilde{u}_{n} = \tilde{\kappa}\epsilon_{r}\tilde{E}_{n} - s(\tilde{t}, \tilde{z})\frac{\tilde{\sigma}}{\tilde{\tau}_{\text{jet}}} \qquad (1.84) $$

where

$$ s(z, t) = \frac{1}{\tau_{\text{jet}}} \Theta(|z| - z_0) \Theta(t - t_0) \Theta(t_0 + \Delta t - t) \qquad (1.85) $$

<footnote n="13">¹³We will present in section 1.5.4 a numerical method that does allow simulating the jet formation</footnote>

<footnote n="14">¹⁴While the choice of $\Delta Q$ is suggested by direct measurements, the choice of $\Delta E_{\text{kin}}$ may be deduced indirectly from the thickness of the jet and from the deformation pathway of the subcritical drop. We expect $\Delta E_{\text{kin}}$ to decrease with thinner jets. Where the jet is clearly visible on the experimental snapshots, we remove usually as much as $\Delta E_{\text{kin}}/E_{\text{kin}} = 50\%$ of the kinetic energy. In the case where no observable jet is formed, like for negatively charged droplets, we expect the loss of kinetic energy during charge emission to be negligible so that $\Delta E_{\text{kin}}/E_{\text{kin}} = 0$.</footnote>
is a window function that limits the source term in time and space, namely to the tip surface $|z| > z_0$ and to the time interval $\Delta t$, with $\Theta()$ being the step function.

The kinetic energy loss is simulated by adding a source term to the Cauchy momentum equation. The mass loss due to jet emission may be expressed by the continuity equation with the source term $\rho s > 0$

$$ \partial_t \rho + \vec{\nabla} \cdot (\rho \vec{u}) + \rho s = 0 \qquad (1.86) $$

The Cauchy momentum equation describes the non-relativistic momentum transport in the droplet volume. In convective (or Lagrangian) form it is written:

$$ \rho \left( \partial_t \vec{u} + \vec{u} \cdot \vec{\nabla} \vec{u} \right) = \vec{\nabla} \cdot T \qquad (1.87) $$

where the stress tensor $T$ accounts for conservative and non-conservatives (viscous) contributions. We can always add a zero contribution $0\vec{u}$ to the momentum equation,

$$ \rho \left( \partial_t \vec{u} + \vec{u} \cdot \vec{\nabla} \vec{u} \right) + \vec{u} \left( \underbrace{\partial_t \rho + \vec{\nabla} \cdot (\rho \vec{u}) + \rho s}_{0} \right) = \vec{\nabla} \cdot T \qquad (1.88) $$

which yields the Cauchy momentum equation (conservative form) with source term

$$ \partial_t (\rho \vec{u}) + \vec{\nabla} \cdot (\rho \vec{u} \otimes \vec{u}) = \vec{\nabla} \cdot T - \rho s \vec{u} \qquad (1.89) $$

Usually $\Delta t$ is chosen such that the loss of charge $\Delta Q$ due to equation (1.84) equals the requested charge loss. So with this method, the information of the charge loss $\Delta Q$ is still required and cannot be deduced from simulations because the charge loss is not done with the correct velocity field that would emerge when simulating explicitly the jet. In other words, the source term is active until the requested charge loss $\Delta Q$ is attained.

Simulations of charge loss and kinetic energy loss with the second method are underway, but all the presented results where achieved using the first method.

# 1.5 Simulated Results with VCVPF

## 1.5.1 Influence of the dimensionless viscosity Oh on the deformation pathway

We first consider the case where the conductivity of the fluid is assumed sufficiently high so that, during the drop deformation, the electric potential remains
<figure><img src="image_29.png" /><figcaption>Figure 1.25: Deformation pathway of a critically charged droplet as a function of the Ohnesorge number Oh. At the tip appearance ($n = 1.5$), 20% of the initial charge is removed $X \to 0.64$ and kinetic energy reduced by 80%.</figcaption></figure>

an equipotential $V(\vec{s}, t) = V(t)$ at the surface. As a result, at the boundary we have $E_τ = 0$ and $\vec{P}_c = \tilde{\sigma}^2$. We simulated for different Ohnersorge numbers Oh (eq. 1.69), ranging from 0.001 to 0.8, the deformation pathway on critically charged droplets ($X = 1$). The results are shown in Fig. 1.25. The initial shape was a spheroid with a small eccentricity of $e = 0.1$ (almost spherical). We found that with increasing Oh and thus with increasing viscosity of the liquid, the droplet adopted a more and more elongated shape before pointed tips appeared at the ends, typically when $n \approx 1.5$. For Oh=0.001, the maximal elongation was characterized by the ratio $a/b = 2.66$, close to the value of 2.6 found for inviscid fluids (see figure 1.23). At the other extremity, for Oh=0.8, the aspect ratio $a/b$ was 3.82, close to the factor 3.86 found by Betelù et al. [46] in the case of the Stokes flow approximation (creeping flow with vanishing inertia forces). The latter is an approximation of the incompressible Navier-Stokes equation is the limit of large Oh $\gg 1$. The aspect ratio found for $n = 1.5$ as a function of Oh $\in [0, 1.4]$ is summarized in Fig. 1.26. The latter is compared to the VPF model, which shows
a different trend. In particular, VPF exceeds the asymptotic aspect ratio of 3.86 found in the Stokes flow approximation for larger Oh, confirming that VCVPF may be better suited for a larger range of viscosities.

<figure><img src="image_30.png" /><figcaption>Figure 1.26: Aspect ratio *a*/*b* of the droplet at charge break-up, (taken at for *n* = 1.5) as a function of the dimensionless viscosity *Oh*. Pink empty boxes 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 maximal elongation of the drop in the case of Stokes flow approximation [46].</figcaption></figure>

For the simulated pathways in Fig. 1.25, once the tip has formed, 20% of the initial charge are suddenly removed and the fissility parameter drops instantaneously to $X = 0.64$. At the same time, we reduced the kinetic energy by 80 % $^{15}$. After charge emission, the pressure due to the surface tension dominates the Coulomb pressure at the surface and the tips relax. This is represented by the increase of the parameter $n$ from 1.5 to up to 3.0 depending on the Oh number. Because of viscous forces, the accumulated kinetic energy is dissipated and the drops returns finally to the spherical shape, unlike the simulations done with the inviscid liquid (1.23). The shape deformation of viscous drops shows a similar behaviour than the damped forced oscillator. $Oh \simeq 0.2$ corresponds to the value where the droplet is critically damped, meaning that it does no longer overshoot into the

<footnote>$^{15}$This is an arbitrary choice and motivated here so as to allow simulating also very low viscous drops $Oh= 0.001$ for which the low energy dissipation could otherwise prevent the computation of the return path of the low-viscous drop.</footnote>
oblate domain when returning to the asymptotic spherical shape. For lower Oh, the dynamics is underdamped 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 and rounded tips. For Oh> 0.2, the motion is overdamped and the droplet returns slowly to steady state spherical shape. In the particular case of Oh=0.8, the return path is even characterized by $n \le 2$, meaning that the tips stay pointed even after charge emission. As in this case, Oh= 0.8 > Oh$_c$ = 0.56, the prediction may not be meaningful.

## 1.5.2 Comparison of VCVPF calculations with experimental results

### Glycol and supercooled water droplets

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

Figure 1.27: Simulated deformation pathway (dotted pink and dashed black lines) compared to experimental data (blue "x" and red "+" standing for supercooled water and glycol droplets respectively. At charge break-up, ($n = 1.5$ for Oh=0.06 and $n = 1.55$ for Oh=0.25), the charge was reduced by 20% ($X = 0.64$) while the kinetic energy was reduced by 20%.

We superpose in Fig.1.27 the VCVPF calculations with Oh=0.06 and Oh=0.2 to the experimental data already presented in the previous sections, namely glycol
droplets and supercooled pure water droplets. We found a remarkable agreement between Glycol and the Oh=0.25 VCVPF curve. For the supercooled water droplets, the VCVPF curve with Oh=0.06 seems to slightly overestimate the damping on the "way back" but gives an overall excellent agreement. Knowing the radius, surface tension and density of the glycol<fnref n="16" /> and water droplets<fnref n="17" />, we can deduce a viscosity of 9 cP for Glycol and 2.5 cP for supercooled water. In the case of glycol, the deduced viscosity is below the expected one by about 50%. The reason is that VCVPF tends to clearly overestimate the damping for Oh> 0.2. Also, in the case of the Glycol drop, the deformation may not have been observed at room temperature, but slightly above. Indeed, Glycol droplets were injected at about 80° C, then cooled down due to evaporation in the trap for about 3 minutes before reaching the critical fissility. The final temperature of the droplet for which the critical value is reached may thus be slightly above 25° C.

<figure><img src="image_32.png" /><figcaption>Figure 1.28: Snapshots of critically charged droplet shapes at different times before and after charge emission simulated with the VCVPF model. We used here Oh=0.25 and the perfect conductor approximation. The same red lemon shape (same aspect ratio) was used for fitting the experimental snapshot and the simulated one.</figcaption></figure>

We took snapshots of the simulated droplet deformation for Oh=0.25. The snapshots are compared in Fig.1.28 to the experimentally observed shapes of a critically charged glycol drop of same diameter, namely $\emptyset = 48\mu m$. We note the

<footnote n="16">¹⁶With $(\rho, \gamma, R) = (1.1 \text{ g/cm}^3, 48 \text{ mN/m}, 24 \text{ µm})$ for the glycol droplet, Oh=0.25 yields a viscosity of $Oh\sqrt{\gamma\rho R} \approx 9 \text{ cP}$, which is 1/2 of the expected value 18 cP of Glycol at room temperature.</footnote>

<footnote n="17">¹⁷With $(\rho, \gamma, R) = (1.0 \text{ g/cm}^3, 78 \text{ mN/m}, 34 \text{ µm})$ for the water droplet, Oh=0.06 yields a viscosity of 2.5 cP, which is close to the expected viscosity of supercooled water droplets of 2.2 cP at -5° C.</footnote>
excellent agreement. In particular the shapes after charge emission are now well
reproduced by the VCVPF model, and the droplet returns to its initial spherical
shape without overshooting. Both, the simulated and observed shape at t = 0 are
compared to the same lemon profile 1.57 with A = 4.97 and B = 4.27 and shown by
the red contour shapes in Fig.1.28. This indicates that the observed and simulated
shapes, at the moment of charge emission, agree and that the VCVPF model is
sufficient to reproduce the observed data. The "return path" to the sub-critical
spherical shape is also well reproduced.

Glycol-Water mixture

<figure><img src="image_33.png" /><figcaption>Figure 1.29: Simulated deformation pathway for two different viscosities. Black squares are correspond to observed glycol-water droplets and fitted by a superellipsoid. At charge break-up (n = 1.52), the kinetic energy was reduced by 30% and 25% of the initial charge was emitted.</figcaption></figure>

We showed in Fig.1.18 and Fig.1.17 the experimental data of critically charged
droplets made out of a mixture of 60% water and 40% of glycol. The viscosity of
those droplets at room temperature is expected to be around 3 cP. We simulated the deformation pathways of a droplet for various viscosities. We found that best agreement with the experimental data was obtained for droplets having an Ohnesorge number between 0.075 and 0.1, see Fig.1.29. Taking the same initial radius of the droplets than the experimental one (ø = 26µm), and taking the viscosity and surface tension of the glycol-water mixture, γ = 58 N/m and ρ = 1060 kg/m³, we get that Oh=0.075 and Oh =0.1 corresponds to fluids having a viscosity of 3 and 4 cP, respectively. Remarkably, the simulations corroborate the estimated viscosity of the mixture!

Pure water droplets

<figure><img src="image_34.png" /><figcaption>Figure 1.30: Comparison between simulated (top) and observed (bottom) snapshots of
ø = 52µm deionized water droplets at room temperature. Initially, the red droplets are
charged at the Rayleigh limit X = 1. At charge break up, 20 % of the charge is removed
and the fissility for the green drop is X = 0.64. Shape labels: (a) lemon, (b) peanut, (c)
diamond.</figcaption></figure>

Finally, we took snapshots of the ultrafast shape deformation of pure water
droplets at the Rayleigh limit. At the onset of the instability, the trapped water
droplets have a $\emptyset$ = 52µm. Assuming that water at room temperature has a viscosity of 1 cP, a density of 1000 kg/cm³ and a surface tension of 0.072 N/m, the droplets are characterized by a low dimensionless viscosity of Oh=0.023 (see table 1.2), and are thus particularly good candidates for the simulation with our VCVPF code. As mentioned before, no jets were detected on the snapshots, neither for positively nor for negatively charged water droplets at room temperature. 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 showed in Fig. 1.30 the deionized water droplets were positively charged. The measurements in [52] indicate that positively charged pure water droplets emit at room temperature about 20 % 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}$ S/m and a charge mobility of $\lambda = 37 \times 10^{-8}$ m²/s/V (see table 1.2). We compare in Fig. 1.30 the simulated shapes to the observed ones. The red simulated shapes in Fig.1.30 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 [52], 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.

### 1.5.3 Influence of the charge carrier mobility

In this section, we focus on the influence of the charge mobility of excess charge carrier on the deformation pathway. The mobility of the charge carriers in the liquid may be such that the surface $\tau_s$ and bulk $\tau_b$ charge relaxation times (eq.1.70) are comparable to the characteristic deformation time $\tau_0$. In that case, the liquid cannot be considered as a perfect conductor and the charges at the surface may not be at equilibrium for a given shape (the Coulomb potential is not an equipotential at the surface). We consider the particular case where the dimensionless bulk
<figure><img src="image_35.png" /><figcaption>Figure 1.31: Deformation pathway of a droplet with Oh=0.023 (typically a 50 µ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.). Except for the latter, we used the bulk conductivity of deionized water of 5.5 µS, leading to a dimensionless bulk conductivity of $\tilde{\kappa}$ = 0.1. At charge break-up ($n$ = 1.5) the charge was reduced by 20% ($X$ = 0.64) while the kinetic energy was reduced by 80%.</figcaption></figure>

conductivity $\tilde{\kappa} \ll 1$ is sufficiently small so that the charge dynamics is dominated by the surface current $j_s$ and depends mainly on the dimensionless charge mobility $\tilde{\lambda}$ (1.71). The latter gives the ratio between the dimensionless velocity of the charges dragged by the liquid and the velocity of the charges pushed by the electric field at the interface. In figure 1.31, we show the pathway of a critically charged drop as a function of the 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µm droplet of deionized water at room temperature, see table 1.2. 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 with is indistinguishable from the perfect conductor (p.c) case. At $n$ = 1.5 the charge was reduced by 20% ($X$ = 0.64). The effect of the charge mobility on the return path was found negligible and will not be further discussed here.

The reason why smaller charge mobilities $\tilde{\lambda}$ tend to accelerate the tip formation is explained by the competition between the charged 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, $Oh< 0.03$, the velocity at the surface may be so high that the flow drags a large amount of charges to the tips. 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 large concentration of charges at the tips tends to accelerate the formation of pointed ends$^{18}$. 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.

## Positive vs negative drops

In the last example of this section, we compare the simulated shape deformation of positively and negatively charged droplets to our experimental data. 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. 1.32 the deformation pathway of positively (green squares) and negatively (black circles) charged drops. Surprisingly, the pathways differ and the positively charged droplets show at jet emission a more elongated shape with an aspect ratio of $a/b = 2.9$ than the negatively 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 behaviour comes from the different mobility of the charge carriers in water. Indeed, protons $H^+$ and hydroxides $OH^-$ in aqueous solution have both an anomalously large diffusion coefficients, but with the mobility of $H^+$ being almost twice as large as that of $OH^-$ anions [54, 44].

We simulated the pathway with the VCVPF model, with the surface charge dynamics depending on the mobility of the charge carriers via the dimensionless factor $\tilde{\lambda}$ defined by (1.71) and on the bulk conductivity via the dimensionless factor $\tilde{\kappa}$. In table 1.2, we give the relevant parameters of a droplet of deionized water of radius $R = 26\mu m$ at room temperature used in the simulations. We chose for the charge mobility of positive and negative drops respectively $\lambda^+ = 37 \times 10^{-8}$

<footnote>$^{18}$Pointed tips have a higher electric capacity $C$ than rounded tips so that the Coulomb energy $\frac{1}{2}Q^2/C$ diminishes for pointed tips</footnote>
and $\lambda^{-} = 20 \times 10^{-8}$ (m$^2$/s/V), corresponding to the charge mobility of $H^{+}$ and $OH^{-}$ [44] at room temperature. Positively and negatively charged droplets are distinguished only by the charge carrier mobility and, in both cases, the dimensionless charge mobility dominates the dimensionless bulk conductivity $\tilde{\lambda} \gg \tilde{\kappa}$. As a result, the deformation pathway is expected to depend sensibly on $\tilde{\lambda}$, with the bulk conductivity having merely a negligible effect on the deformation dynamics$^{19}$. The two simulated paths are superposed in Fig.1.32 to the experimental data. We found a remarkable agreement between the simulated and observed pathway. In particular, the larger aspect ratio $a/b$ at jet emission of positively charged drops compared to negatively charged drops is well reproduced by the simulation.

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

<figure><img src="image_36.png" /><figcaption>Figure 1.32: Deformation pathway for positively and negatively charged water droplets having a radius of 26 µm. Full lines are simulations with Oh=0.023. The green line used $\tilde{\lambda} = 4$ and the black used $\tilde{\lambda} = 2$. The emitted charge is given in the legend. Black circles and green squares are experimental results of negatively and positively charged deionized water droplets, respectively. As no jets were found on the snapshots, no kinetic energy was removed from the drop during charge emission.</figcaption></figure>

<footnote>$^{19}$for see water, having a bulk conductivity of 5 S/m, one has $\tilde{\kappa} \gg \tilde{\lambda}$. As a result, positively and negatively charge droplets should yield the same deformation path, unlike deionized water</footnote>
<table><caption>Table 1.2: Parameters of positively and negatively critically charged droplet of deionized (pure) water. Upper table: Parameters, common to positive and negative drops; Lower part: electrical properties depending on the sign of the critical charge $Q$.</caption><thead><tr><th>R (μm)</th><th>μ (cP)</th><th>ρ (kg/m³)</th><th>γ (mN/m)</th><th>ε<sub>r</sub></th><th>Oh</th><th>τ<sub>0</sub> (μs)</th></tr></thead><tbody><tr><td>26</td><td>1</td><td>1000</td><td>72</td><td>80</td><td>0.023</td><td>16</td></tr></tbody></table>

<table><thead><tr><th>Sign(Q)</th><th>κ (μS/m)</th><th>λ (m²/s/V)</th><th>κ̃</th><th>λ̃</th></tr></thead><tbody><tr><td>pos</td><td>5.5</td><td>37 × 10<sup>-8</sup></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>0.1</td><td>2.1</td></tr></tbody></table>

<figure><img src="image_37.png" /><figcaption>Figure 1.33: Upper panel: 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. Lower panel: 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 "+" symbol line gives the equilibrium surface charge density $\tilde{\sigma}_{eq}$ in the case of a perfectly conducting liquid (p.c.).</figcaption></figure>
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.1.33 for three different values of $\tilde{\lambda}$. The quantities in Fig.1.33 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/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.

We found that, for negatively charged drops, removing only 20% of the charge, as we did for the positively charge drop, was not sufficient to reproduce the experimental *way back* pathway. We had to remove 40% in order to reproduce the amplitude of the *way back* path. This is corroborated by the observations reported by R. Müller in his thesis<fnref n="20" />. They found that negatively charged droplets with $R = 31$ µm emit about 40% of the initial charge at charge break up. The reason for the latter may be found in the higher accumulated surface charge at the tips with decreasing charge mobility.

What if instead of pure water we had used sea water? Sea water, having a bulk conductivity of 5 S/m yields a dimensionless conductivity $\tilde{\kappa} \simeq 10^4$ which is much larger than the typical dimensionless charge mobility $\tilde{\lambda} < 10$ of excess charge carrier. The charge dynamics would then 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 charged 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.

**Concluding remarks for VCVPF:** We conclude this section by highlighting the remarkable results.

1. The VCVPF model predicts that the aspect ratio of the shape just before charge break-up is increases with increasing dimensionless viscosity Oh. This is because the tip formation is retarded due to higher viscosity, whereas for lower dimensionless viscosities, the tips are formed earlier. We found that the ratio $a/b$ lies between 2.5 for inviscid liquids to 3.86 for Stokes flows, $a/b \in [2.6, 3.86]$.

<footnote n="20">²⁰Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen, Dissertation by René Müller (2010), Fakultät für Mathematik und Naturwissenschaften der Technischen Universität Ilmenau. https://www.db-thueringen.de/receive/dbt_mods_00015683</footnote>
2. The model can be used to predict the viscosity of a liquid (here mix of water and glycol) by fitting the experimental pathway in the $(a/b, n)$ graph. The predictions are more reliable for lower viscosities, $Oh < 0.2$.

3. In the particular case of critically charged pure (deionized) water drops at room temperature, 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. The excellent agreement gives much confidence in the underlying model and numerical algorithms of the simulations.

4. 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. In the case of pure water droplets, the polarity of the charge can eventually be determined by comparing the aspect ratio $a/b$ of positively and negatively droplets just before charge breakup.

### 1.5.4 Boundary element method

In VCVPF, the Laplace solver uses harmonic functions in the prolate coordinate system. While that algorithm is adapted for simulating the time evolution of a critically charged drop, from the initial spherical shape to the pointed lemon shape and back to the sphere, it breaks down at the jet formation. Indeed numerical errors increase exponentially when the ends become pointed and the simulation usually stops just before the Taylor cone and jet could be formed. The culprit is the loss of orthogonality of the harmonic basis functions $\phi_k$ when the shape of the drop deviates strongly from a prolate spheroid. Some eigenvalues of the overlap matrix $S_{k'k} = \int \phi_{k'} \frac{\partial \phi_k}{\partial n} ds$ tend then to zero which makes the matrix almost singular. Also, an increasingly large number of harmonic functions are necessary to represent finely the tip of the pointed drop shape. As a result, the numerical method based on the expansion of harmonic functions in prolate spheroidal coordinate system is no longer appropriate to solve the Laplace equation 1.46 inside the volume of a droplet near charge breakup. We proposed thus to solve the Laplace equation using the Boundary Element Method (BEM$^{21}$). The advantage of BEM is that it can solve the Laplace equation for wherever regular closed surface enclosing a volume and is thus well adapted to describe for example the cone formation at the droplet tip.

<footnote>$^{21}$The Boundary Element Method is well documented on http://www.boundarielements.com) and http://www.boundary-element-method.com by Stephen Kirkup.</footnote>
<figure><img src="image_38.png" /><figcaption>Figure 1.34: Representation of the sphere by truncated conical panels.</figcaption></figure>

The disadvantage of BEM is that it converges slowly ($\propto N^{-1}$) to the exact solution with increasing number $N$ of elements, increasing thus the CPU cost. To optimize the simulation, the numerical code switches thus from the harmonic expansion method to BEM once the overlap matrix $S_{k'k}$ is close to singular. In practice, we switched to BEM algorithm once the $n$ parameter of the superellipsoid fit falls below 1.54. We discretized the axisymmetric droplet surface $\xi(\eta_i)$ using 400 points that clustered close to the apex, so that the density of points close to the apex was higher than close to the equator. Because of axisymmetrie, each BEM element is thus a truncated cone defined by two successive points at the surface (see Fig. 1.34 for example). We used cubic splines[63] to interpolate the surface at any arbitrary point$^{22}$ and to calculate the mean curvature at the grid points.

## Simulated Results

The BEM algorithm allowed us to simulate the tip formation and even the onset of the jet emission of a droplet. In Fig.1.35 we show the time evolution (black lines taken at successive time steps) of the tip of a critically charged droplet with Oh=0. Interestingly, we found that the half-opening angle of the tip just before break-up, approaches the Taylor angle of 49° [8], with an estimated angle of 48.3° near the tip, confirming Taylor's prediction of a conical tip. This result is a refinement of the one shown in Fig.1.22, where an estimated tip angle of 39° was found. Please note the dimensionless x-scale of Fig.1.35; The Taylor angle is measured within a distance of 0.03 from the symmetry axis and the angle drops quickly beyond that region. BEM succeeded also to simulate the onset of the jet formation. The latter is characterized by a kink around 0.01 from the symmetry axis and an exponentially fast increasing of the tip height. Unfortunately, the BEM methods also breaks down at this point as the curvature becomes near singular. The finite thinness of the jet seems to be spurious and limited by the numerical resolution

<footnote>$^{22}$We used zero first derivative at the end points</footnote>
<figure><img src="image_39.png" /><figcaption>Figure 1.35: Simulated time-evolution of the tip shape. *x* and *z* are dimensionless, normalized to the initial droplet radius.</figcaption></figure>

of BEM (the 400 points). In the case of Oh> 0, the jets are expected to have a
finite thinness as was observed in the experiments. Simulations with Oh> 0 have
not been done yet, as some numerical adjustments needed to be done.

**1.5.5 Finite element method**

BEM is a promising technique for solving the Laplace equation for arbitrary
droplet shapes in irrotational flow approaches like VCVPF or Stokes flows like
[46]. However, fluid dynamics of viscous droplets is best described by the full
Navier-Stokes (NS) equation. The latter is usually solved using a Finite Element
Method (FEM<fnref n="23" />), where the boundary conditions are imposed weakly. With Dr.
Radcliff Allister, who occupied the 1-year post-doctoral position founded by the

<footnote n="23">²³The Finite Element Method is introduced by:</footnote>

<footnote>Daryl L. Logan (2011). *A first course in the finite element method*. Cengage Learning. ISBN 978-0495668251.</footnote>

<footnote>Reddy, J. N. (2006). *An Introduction to the Finite Element Method* (Third ed.). McGraw-Hill.
ISBN 9780071267618.</footnote>
ANR DYNAMIC 2007, we proposed to describe the deformation dynamics and Coulomb explosions using a non-conforming, discontinuous Galerkin finite element procedure with moving meshes. Assuming xOy plane symmetry for the droplet deformation, the droplet is discretized in FEM by a mesh of triangles Fig.1.36. The simple equations that model these finite elements are then assembled into a larger system of equations that models the entire problem. FEM then uses variational methods from the calculus of variations to approximate a solution by minimizing an associated error function. Unlike the (VC)VPF models, FEM calculates explicitly the fluid dynamics in the volume. As a result, FEM is clearly more CPU expensive than the (VC)VPF models. However, it is not restricted to potential fluids and is indeed well appropriate for resolving the Navier-Stokes equation for a real (viscous Newtonian fluid). The full Navier-Stokes equations are discretized

<figure><img src="image_40.png" /><figcaption>Figure 1.36: Discretization of the surface by triangles for FEM calculations. Initial mesh of 1908 finite elements used for the prolate deformations. Note the structured "zig-zag diagonal" nature of the mesh elements abutting the axis of rotational symmetry. Extracted from [58].</figcaption></figure>

and solved in an Eulerian manner with a simple predictor-corrector Lagrangian updating of the free boundary location coincident with the droplet surface at each solution iteration. We used linear Crouzeix-Raviart basis functions for the velocity and piece-wise constant pressures. The use of the original and ever popular $P_1^{NC} - P_0$ Crouzeix-Raviart non-conforming velocity-pressure finite element pair [12] was made in anticipation of future developments for flexible mesh interpolation polynomial ("h-p") refinement about the droplet tip during jet formation
and to avoid introducing the additional high-order singularities that appear in the axisymmetric stabilization terms required for otherwise unstable $P_1 - P_0$ or equal-order $P_1 - P_1$ conforming velocity finite element schemes [31]. To avoid spurious mass build-up along inter-element edges and to impose the zero tangential traction applied by two different elements either side of such shared edges equal and opposite (see eq. 1.25), it is necessary to introduce extra terms to penalize any velocity and traction jumps across these internal edges. This is a usual requirement of non-conforming methods and is typically parameterized by a scalar penalty parameter $\beta$. The greater this parameters value, the more severe such jumps are punished, but also the more ill-conditioned the final system matrix becomes. A judicious choice of $\beta = 0.07$ was found to be appropriate for the given problem and its value is discussed in [58].

One of the principle difficulties encountered in the simulation of the droplet lay in the calculation of the mean curvature $H$ (eq. A.14) of its surface. The fourth-order polynomial stencil fitting was in most instances found to be adequate. However, it should be noted that unless using very coarse meshes, the nodes on such stencils could sometimes become approximately collinear and the calculation of mean curvature become very ill-conditioned<sup>24</sup>. The numerical code is validated with the results, presented both for the simple (oscillatory) relaxation of elongated electrically (un)charged droplets to a sphere and for the deformation of supercritically charged initially slightly oblate spheroidal droplets.

## Simulated Results

In a first step, we confirm that the underlying dynamics of the charged droplet are correct, validating the finite element numerical model used. This is performed by considering the relatively small viscous damped oscillations of a prolate elongated droplet after a deformation induced by a small initial irrotational quadrupole velocity field of dimensionless amplitude of $u_0 = 0.07$, $\vec{u} = u_0 \nabla(r^2 P_2(\cos(\theta)))$. Figure 1.37 shows the computed results for just this initial perturbation using the non-conforming scheme described in [58], with blue and red vertical gridlines marking the expected semi-periods for each case respectively. The expected periods of the droplet oscillations with $X = 0.5$ and $X = 0$ were calculated using 1.2. The horizontal displacement for the gridlines was chosen so that the leftmost would correspond with the first peak in each case. We distinguish 3 subsets in Fig.1.37, one with an "s" suffix, one with an asterisk, and one bare. For the "asterisk" subset, the zero shear stress boundary condition was weakly imposed by the penalty parameter $\beta = 0.007$ while the said boundary condition was not imposed

<footnote n="24"><sup>24</sup>This problem was avoided in the VCVPF approach, at least as long no singular tips are formed, as we used a spectral method based on Gauss-Legendre collocation points to interpolate the droplet surface allowing to compute accurately the mean curvature $H$.</footnote>
<figure><img src="image_41.png" /><figcaption>Figure 1.37: Oscillatory relaxation to a sphere after an initial launch with an initial quadrupole velocity of (dimensionless) amplitude of 0.07 for both a charged ($X = 0.5$) and uncharged ($X = 0$) liquid droplet with $\Omega_h = 0.04$.</figcaption></figure>

for the "bare" set. We see that in the low amplitude oscillatory regime, the zero shear stress boundary condition has little to no influence on the damping, as it is not severely violated when not imposed. This is corroborated by the "s" subset with stand for the VCVPF calculations (irrotational flow) and for which the zero shear stress boundary conditions is also not satisfied. The comparison between the VCVPF and FEM approach highlights rather the spurious damping of the FEM method at very low amplitudes due to the numerical errors in the evaluation of the mean curvature at the piece-wise constant discretized surface. We note here that a possible numerical improvement would be to use elliptic mesh generation<fnref n="25" />, allowing a continuous boundary second derivative to be defined easing problems with the calculation of mean curvature and helping suppress surface irregularities.

<footnote n="25">²⁵"Elliptic grid generation" signifies generation of generalized curvilinear coordinate grids through solution of elliptic partial differential equations (PDEs). Usually, such grids are fitted to bounding bodies. Such a grid is smooth and has continuous first and second derivatives, grid lines are appropriately stretched or clustered, and grid lines are orthogonal or nearly so over most of the grid domain.</footnote>
FEM as presented here nevertheless reproduces nicely the damped oscillations of
the viscous drop.

<figure><img src="image_42.png" /><figcaption>Figure 1.38: Simulated full deformation cycles for critically charged ($X = 1$) prolate droplets of various viscosities (given by the Ohnesorge values) and compared against experimental results for Glycol (Oh$\approx$ 0.2) and super-cooled water (Oh$\approx$ 0.09). For the return, $X = 0.4$, meaning about 36% of the cahge was removed.</figcaption></figure>

Confident of our FEM approach, we simulated the deformation pathway of
droplets of various dimensionless viscosities (Oh). Similar to VCVPF, the elon-
gation ($a/b$) of the droplet before pointed ends were formed increased with Oh.
Unfortunately, our attempts to simulate the jet formation, were not successful as
we run into numerical problems, mainly for the evaluation of the mean curvature
close to the tip, owing to significant displacements of the boundary nodes and to
the appearance of slight surface kinks. Also the extreme sensitivity of the charge
distribution to any surface irregularities worsened the numerical stability of the
code. To cure partially the problem, a simple volume conserving re-positioning of
boundary nodes was occasionally used to equalize local boundary node separation
distances when these were significantly different and to remove surface indenta-
tions. The simulated deformation pathways (Fig.1.38) differ significantly from the
VCVPF calculations (Fig.1.25), showing a much too strong influence of Oh on the
maximal elongation in FEM than in VCVPF. For example, the path for Oh=0.1
(blue curve) and Oh=0.09 (green curve) are quite different for only a 10 % variation
of the dimensionless viscosity Oh. However, the simulations appear to capture the
expected general deformation path. At this point, the simulations with the full NS
approach solved via FEM must be considered as preliminary and the effect of the
penalty parameter $\beta$ on the pathway must be studied. Also the numerical errors
arising in the evaluation of the mean curvature and surface charge distributions
may have a non-negligible influence on the pathway.

1.6 Summary and Conclusion

In conclusion, we have studied the deformation dynamics of critically charged mi-
crodroplets from an experimental and theoretical point of view. The project was
founded by the ANR 2007 and allowed us to acquire the necessary material to
inject, trap and take snapshots of fast deforming droplets at the Rayleigh limit.
We inherited a Paul trap of the first generation graciously given by the group of
Prof. Thomas Leisner in Ilmenau and profited from their expertise to to set up
quickly the whole experimental apparatus. Within 6 months (taking into account
the incompressible time for purchasing the needed material), we had a function-
ing setup similar to the one used by the group of T. Leisner. After spending
some unsuccessful time to inject mercury droplets, we succeed in injecting and
trapping mixed water-glycol droplets. Our first snapshots taken in CIMAP and
fitted to a superellipsoid allowed to compare in detail the effect of the viscosity
of the deformation pathway. Compared to pure glycol droplets, our water-glycol
droplets showed mainly a shorter elongation of the drop at charge emission, while
the amount of emitted charge is comparable. We further injected positively and
negatively charged pure water droplets. Their low viscosity confirmed the trend
that with decreasing viscosity the droplets are less elongated at charge break-up.
Moreover the study put into evidence the effect of the influence of the charge car-
riers mobility in the liquid on the deformation pathway. Remarkably, we observed
no jet emission from pure water droplets at room temperature.

At the same time we started the development of our numerical codes based on
a potential flow approximation. In a first step, we developed a numerical model
based on the PF model, valid for inviscid liquids. This allowed us to validate (i)
the time propagator of the surface vector and (ii) the evaluation of the Coulomb
pressure and mean curvature at the surface. We checked for instance the conser-
vation of the total energy and volume of the drop. It showed however quickly its
limitations as the lack of viscosity did not allow to compute the deformation after
charge emission. Also the numerical approach did not allow to compute the jet
formation. In a second step, we improved our model by adding a dissipative term
proportional to the viscosity of the liquid. We based our code on the VCVPF approach, which is discussed in section 1.4.6. We presented an original numerical technique to include the effect of the shear stress in the pressure term, that is an improvement of the technique proposed by Joseph et al. [40] for small surface oscillations. Our VCVPF simulations corroborated the main results of the observed data and highlighted in particular a clear trend of the ratio *a/b* of the drop at charge emission as a function of the viscosity. We found that the ratio varies from *a/b* = 2.6 for an inviscid liquid to *a/b* = 3.9 for a high viscous liquid. The high viscous limit is in agreement with the calculation done for a Stokes flow (creeping liquid), giving confidence in our VCVPF model. The latter keeps the CPU advantage of the PF approach, where the velocity field is obtained from a scalar function defined on the axisymmetric surface (1D scalar function)<fnref n="26" />

Our study on pure water droplets as well as the PhD work by R. Müller [52], showed that, while no jet emission was observed for charged water droplets at room temperature, jet emission and daughter droplets were observed for cooled water droplets, but then only for positively charged drops, not for negatively charged ones. These observations motivated us to refine our numerical models. We lifted the perfect conductor assumption and added the dynamics of the surface charge distribution to the time evolution of the droplet surface. The charge dynamics was assumed to depend on the mobility of the charge carriers at the surface and on the bulk conductivity of the liquid. We demonstrated that the deformation pathway depends on the mobility of the charge carriers at the surface, but only in the cases where the charge relaxation rate of the bulk conductivity can be neglected (which is the case for pure water droplets). Our model succeeded thus to explains clearly why positively charged pure water droplets are more elongated than negatively charged droplets at charge emission.

While our VCVPF model was able to describe the whole deformation pathway for droplets with different viscosities and charge mobilities, we were still not able to simulate the jet emission. Even the tip formation, just before jet emission was badly described be our numerical technique. We added thus to our code a module that solves the Laplace equation with the BEM technique, where the surface is discretized by 400 points distributed non-uniformly, with a higher density of points close to the ends. BEM has a rather high CPU cost and was activated only when the shape approached the tip formation, that is, when the parameter *n* of the superellipsoid was below 1.54. We found the remarkable result that the jet is indeed emitted from a Taylor cone. The half-opening angle of the BEM simulated pointed tip, just before jet formation, was 48.3°, close to the predicted Taylor angle of 49°. Unfortunately we did not succeed to simulate also the jet, as we ran

<footnote n="26">²⁶The Navier Stokes approach needs to propagate a 2D vector field, defined in the volume of the axisymmetric drop.</footnote>
into numerical instabilities in the evaluation of the pressure terms, which became
near singular near the tip. The development of an improved numerical algorithm
is underway.

A. Radcliffe, who occupied a post-doc position in CIMAP for 15 months, de-
veloped a numerical code based on non-conforming finite elements to solve the
2D axisymmetric Navier-Stokes (NS) equations. Finite Element Method (FEM)
have an even higher cost than BEM. For the FEM, we chose triangles which are
implemented in standard libraries. First results showed that the chosen algorithm-
mic succeeded to simulated grossly the surface oscillations on viscous drops. The
code was particularly accurate for high viscosities, but less for lower viscosities.
As such the NS-FEM approach is complementary with VCVPF, which is partic-
ularly accurate for low viscosities. We tried to simulate the jet emission with the
NS-FEM approach but ran into the same problem than with VCVPF, that is,
the pressure terms are evaluated with increasing numerical error close to the tip
formation. A possible solution to this issue is to use the algorithm presented by
Collins et al. [50], which uses FEM with elliptic mesh generation and adaptive
implicit time integration to solve the NS equations. It allows unparalleled accu-
racy in the evaluation of the mean curvature as the technique is able to resolve in
a single simulation length scales that differ by four orders of magnitude. Alter-
natively, Gañán-Calvo et al. used the promising Volume of Fluid (VoF) method
for tracking and locating the free surface (fluid-gas interface) [61]. The VoF be-
longs to the class of Eulerian methods which are characterized by a mesh that is
either stationary or is moving in a certain prescribed manner to accommodate the
evolving shape of the interface.

This 5 year study, involving a small group in CIMAP, namely R. Rangama in
charge of the experimental part, A Radcliffe in charge of the NS-FEM approach
and myself, tackled the Coulomb instability of charged droplets and their amazing
shape deformations. Using various numerical approaches, we succeeded in repro-
ducing with our simulations the observed shape deformations of various liquids
and highlighted trends of shape deformation as a function of the viscosity, charge
mobility and amount of emitted charge. In the case of negatively charged pure
water droplets, where no jets are observed, our simulations are in excellent agree-
ment, giving high confidence in our numerical approach. In a future step, a new
algorithm may allow to study also the jet emission as well as the amount of emitted
charge, mass and lost kinetic energy.
# Perspective

## 1.7 Rotationally Corrected Viscous Potential Flow (RCVPF)

While, in VCVPF, we recover the correct power of traction, the dissipation is still calculated using a purely irrotational velocity vector field. The viscous forces in the volume are still absent and the shear stress boundary condition is still not satisfied. We propose in this section to further refine the VCVPF approach by adding a rotational field to the potential velocity field that accommodates adiabatically so as to satisfy the zero shear stress boundary condition at each time step. We call this original approach the RCVPF approach. In the following we give some preliminary development of the RCVPF model.

In order to overcome the limitations of the previous potential flow models, we add a rotational axisymmetric velocity field to the irrotational velocity field $\vec{\nabla}\Psi(\xi, \eta)$. In a prolate spheroidal coordinate system $(\xi, \eta, \varphi)$ it reads

$$ \vec{u}(\xi, \eta) = \vec{v}(\xi, \eta) + \vec{w}(\xi, \eta) = \underbrace{\vec{\nabla} \Psi(\xi, \eta)}_{\vec{v}} + \underbrace{\vec{\nabla} \wedge (F(\xi, \eta)\hat{e}_{\varphi})}_{\vec{w}} \quad (1.90) $$

where $\Psi(\xi, \eta)$ and $F(\xi, \eta)$ is a scalar function. The decomposition (1.90) is known as the Helmholtz decomposition of an axisymmetric vector field. Because of the xOy plane symmetry, $\Psi$ is an even function and $F$ an odd function with respect to $\eta$. We remind that the incompressible liquid condition requires that $\vec{\nabla} \cdot \vec{u} = 0$, so that that $\vec{\nabla}^2 \Psi = 0$. As was already shown in the previous sections, $\Psi$ can be expanded on harmonic basis functions $\psi_k$. Limiting the expansion to the first L basis elements, the irrotational part of the vector field expands as

$$ \vec{v} = \vec{\nabla}\Psi = \sum_{k=1}^{L} a_k \vec{\nabla}\psi_k = \sum_{k=1}^{L} a_k \vec{k} \quad , \quad (1.91) $$

where we introduced $\vec{k} = \vec{\nabla}\psi_k$ for simplicity. For the rotational vector field $\vec{w} = \vec{\nabla} \wedge F\hat{e}_{\varphi}$, we expand the scalar function $F(\xi, \eta)$ onto the L basis functions $f_{\lambda}(\xi, \eta)$. The latter must satisfy $\Delta f_{\lambda}\hat{e}_{\varphi} \neq 0$ in order to contribute to the rotational field. The rotational part of the axisymmetric vector field is thus expressed as

$$ \vec{w} = \vec{\nabla} \wedge F \hat{e}_{\varphi} = \sum_{\lambda=1}^{L} b_{\lambda} \vec{\nabla} \wedge (f_{\lambda} \hat{e}_{\varphi}) = \sum_{\lambda=1}^{L} b_{\lambda} \vec{\lambda} \quad (1.92) $$

where we introduced $\vec{\lambda} = \vec{\nabla} \wedge (f_{\lambda} \hat{e}_{\varphi})$ to simplify the notation. We can further require the rotational basis elements $\vec{\lambda}$ to be orthogonal to the irrotational vector
basis $\vec{k}$. The orthogonality condition is satisfied if the basis function $f_λ$ are zero on the droplet surface$^{27}$. We define the symmetric tensors $D$ and $W$, standing respectively for the irrotational and rotational contribution of the velocity vector field to the strain-rate tensor :

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

$$W = \frac{1}{2} (\vec{\nabla} \vec{w} + (\vec{\nabla} \vec{w})^t) \qquad (1.94)$$

We remind that $D$ is a divergence-free tensors, unlike $W$. The boundary conditions of the normal and tangential stress now write:

$$T_{nn} - M_{nn} = -p + 2\mu(D_{nn} + W_{nn}) = -P_s + P_c \qquad (1.95)$$

$$T_{\tau n} - M_{\tau n} = 2\mu(D_{\tau n} + W_{\tau n}) = \sigma E_{\tau} \qquad (1.96)$$

With the additional rotational field we have this times enough degrees of freedom to satisfy explicitly the zero shear stress boundary condition (1.96)! The ad hoc viscous pressure $p_μ$ found in VCVPF is replaced here by the strain rate tensor $W$ based on the rotational field $\vec{w}$. The power equation in RCVPF is obtained by projecting the Navier-Stokes equation 1.12 onto the (irrotational) vector $\vec{v}$.

$$\frac{dE_{\text{kin}}}{dt} = \rho \left( \vec{v} \cdot \frac{d\vec{u}}{dt} \right)_{\Omega} = \rho \left( \vec{v} \cdot (\vec{\nabla} \cdot T) \right)_{\Omega} \qquad (1.97)$$

$$= -2\mu (V : (D + W))_{\Omega} + \int_{\partial\Omega} \vec{v} \cdot (T - M) \cdot \vec{n} \, ds^2 \qquad (1.98)$$

$$\qquad (1.99)$$

where $()_Ω$ means integration over the drop volume $Ω$. Injecting the boundary conditions (1.95) and (1.96) yields the power of traction and the power of dissipation

$$\begin{align*}
\int_{\Omega} \vec{\lambda} \cdot \vec{k} \, d\Omega &= \int_{\Omega} \vec{\nabla} \cdot (f_{\lambda} \hat{e}_{\varphi} \wedge \vec{k}) \, d\Omega + \int_{\Omega} f_{\lambda} \hat{e}_{\varphi} \cdot \underbrace{(\vec{\nabla} \wedge \vec{k})}_{0} \, d\Omega \\
&= \int_{\partial\Omega} (f_{\lambda} \hat{e}_{\varphi} \wedge \vec{k}) \cdot \vec{n} \, ds \\
&= - \int_{\partial\Omega} \underbrace{f_{\lambda}}_{0} k_{\tau} \, ds = 0
\end{align*}$$

for arbitrary $k_τ = \vec{k} \cdot \vec{\tau}$. (Note that, as $f_λ$ and $k_τ$ are odd functions at the surface, the product $f_λ k_τ$ is an even function on the surface.)

<footnote>$^{27}$The vectors $\vec{k}$ are orthogonal to $\vec{\lambda}$ in the droplet volume $Ω$ if the functions $f_λ = 0$ on the droplet surface $\partialΩ$,</footnote>
proportional to the viscosity $\mu$,

$$ \frac{dE_{\text{kin}}}{dt} = \underbrace{\int_{\partial\Omega} (v_n(-P_s + P_c) + v_{\tau}(\sigma E_{\tau})) ds^2}_{\mathcal{T}^{\text{RCVPF}} = \mathcal{T}^{\text{VCVPF}}} -2\mu \underbrace{\int_{\partial\Omega} (v_n D_{nn} + v_{\tau} D_{\tau n} + w_n D_{nn} + w_{\tau} D_{\tau n}) ds^2}_{\mathcal{D}^{\text{RCVPF}}} \quad (1.100) $$

A more detailed derivation of the power equation (A.43) is given in the appendix A.5. We recognize the traction integral of the VCVPF model. However, compared to VCVPF, the RCVPF has additional terms in the power of dissipation integral, that involves the rotational vector field $\vec{w}$ evaluated at the boundary. We may further chose the basis function $f_{\lambda}(\vec{s})$ to be zero on the surface. As a result, the basis vectors $\vec{k}$ and $\vec{\lambda}$ are orthogonal in the drop volume and $\vec{n} \cdot \vec{\lambda}(\vec{s}) = 0$ at the surface, so that $w_n(\vec{s}) = \sum_{\lambda} \vec{n} \cdot \vec{\lambda}(\vec{s}) = 0$. In this case, the power of dissipation $\mathcal{D}^{\text{RCVPF}}$ reads,

$$ \mathcal{D}^{\text{RCVPF}} = \mathcal{D}^{\text{VCVPF}} - 2\mu \int_{\partial\Omega} (w_{\tau} D_{\tau n}) ds \quad . \quad (1.101) $$

Within the orthogonality assumption, RCVPF adds, with respect to VCVPF, a single term that involves only the tangential component of the rotational vector field $w_{\tau}$. We will see in the next section that RCVPF allows to recover the correct decay rate of the surface oscillations of a viscous drop, even in the higher viscosity regime, where VCVPF overestimates the rate. RCVPF keeps the advantageous numerical CPU cost of the VCVPF model, where all the quantities are calculated on the interface.

The vector field $\vec{w}$ is activated by the boundary condition (1.96), which can now be satisfied using the additional degrees of freedom given by the set of expansion coefficients $b_{\lambda}$ of $\vec{w}$. Defining the tensor $\Lambda = \frac{1}{2}(\nabla\vec{\lambda} + (\nabla\vec{\lambda})^t)$ so that

$$ W = \sum_{\lambda}^{L} b_{\lambda} \Lambda \quad (1.102) $$

and projecting (1.96) onto the irrotational basis $\{\vec{k}\}$, yields a set of $L$ equations that must satisfy the $L$ coefficients $b_{\lambda}$ of the rotational basis $\{\vec{\lambda}\}$,

$$ \int k_{\tau} D_{\tau n} ds^2 + \sum_{\lambda=1}^{L} b_{\lambda} \underbrace{\int k_{\tau} \Lambda_{\tau n} ds^2}_{R_{\lambda k}} = \int k_{\tau} \sigma E_{\tau} ds^2, \quad 1 \le k \le L \quad (1.103) $$
Solving the latter system of equations implies that the matrix $R_{λk}$ must be non-singular so that $R_{λk}$ can be inverted. As the rotational field $\vec{w}$ that enters the RCVPF model is defined only by the shear stress boundary condition (1.103), we have to restrict $\vec{w}$ to a subspace which yields a non-singular matrix in (1.103). The RCVPF model depends thus on the choice of the subspace of the rotational field that spans $\vec{w}$.

We also want to emphasize that if the rotational part of the velocity field vanishes, $\nabla \wedge \vec{w} \to 0$, the RCVPF model tends toward the VPF model and not toward the VCVPF. This becomes clear, when looking at equation A.45 in the appendix. Indeed, if $W \to 0$, the shear stress condition of RCVPF cannot be satisfied any longer, similar to VPF, and is thus ignored. On the other hand, if the diffusion length of vorticity layer is negligible compared to the droplet radius, so that the vorticity layer is limited to the surface, RCVPF tends to the VCVPF model. This behaviour will be illustrated in section 1.7.

To summarize, the RCVPF model allows to add a dissipation term originating from the rotational part of the field, that is activated by enforcing the shear boundary condition at the free boundary. The numerical approach of the RCVPF model is detailed in the appendix A.6.

## Choice of the rotational basis set

We suppose that the viscous force originating from the tangential shear stress acts inside a thin vorticity layer at the interface. In the case of small quadrupole surface oscillations, we can define the vorticity diffusion lengths $l_μ$ [51]

$$l_μ \approx \sqrt{\frac{\mu}{\rho\omega(\lambda)}} = R\sqrt{\text{Oh}}(\lambda(\lambda + 2)(\lambda - 1))^{-1/4} \quad (1.104)$$

where $\omega(\lambda) = \tau_0^{-1}\sqrt{\lambda(\lambda + 2)(\lambda - 1)}$ is the dispersion relation ship in the PF model$^{28}$. The diffusion lengths corresponds to the thickness of the vorticity layer. For an axisymmetric problem expressed in an orthogonal coordinate system, the rotational part of the velocity field is expressed by (1.92). In prolate coordinate system we suggest the following expansion$^{29}$

$$F(\xi, \eta) = \sqrt{(1 - \eta^2)} \sum_{\lambda} b_{\lambda} \sqrt{(\xi^2 - 1)} P'_{\lambda+N}(\xi) P'_{\lambda}(\eta) \quad (1.107)$$

$$\omega(\lambda) = \text{Im}(\Omega^{\text{VPF}}) = \tau_0^{-1} \sqrt{\lambda(\lambda + 2)(\lambda - 1) - \lambda^2(\lambda - 2)^2 \text{Oh}^2} \quad (1.105)$$

<footnote>$^{28}$Alternatively, we could use the dispersion relationship of the VPF model:</footnote>

<footnote>$^{29}$The rotational field can be made orthogonal to the potential velocity field $\vec{v} = \vec{\nabla}\Phi$ by modifying the the basis functions such that they are zero at the axisymmetric surface defined by</footnote>
with $N \ge 1$ being a positive integer. The choice of $F$ is motivated by the fact that with larger $N$, the radial part of $F$ is non-zero mainly close to the interfaces, meaning that the vorticity layer near the interface is narrower, or that the vorticity layer has a smaller diffusion length at the interface. We are looking thus for a relationship for $N$ of the form

$$N(\text{Oh}, \lambda) \approx \frac{R}{l_{\mu}} = \frac{(\lambda(\lambda + 2)(\lambda - 1))^{1/4}}{\sqrt{\text{Oh}}}. \quad (1.108)$$

We simulated small amplitude surface oscillation as a function of time with the RCVPF model and extracted the damping of the quadrupole mode, $\lambda = 2$. We tuned the free parameter $N$ so as to match the decay time obtained from the exact dispersion relationship of the quadrupole of a viscous drop [51]. Figure 1.39 shows the matching $N$ (black points) as a function of $\sqrt{\text{Oh}}$. We found that the points are well fitted by the a law $N(\text{Oh}) + 1 = 8^{1/4}\text{Oh}^{-1/2}$, in agreement with (1.108)

<figure><img src="image_43.png" /><figcaption>Figure 1.39: Integer $N$ versus the dimensionless viscosity $\text{Oh}$</figcaption></figure>

## First simulations with RCVPF: Damped surface oscillations

We illustrate the RCVPF model in the case of small surface oscillations (mainly quadrupole moment) of a droplet. We compare the damped oscillations given by the RCVPF approach to the VCVPF one and to the exact damping given by

the function $\xi_s(\eta) > 1$,

$$F(\xi, \eta) = \sqrt{(1 - \eta^2)} \sum_{\lambda} b_{\lambda} \left( \sqrt{(\xi^2 - 1)} P'_{\lambda+N}(\xi) - \frac{(\xi_s^2(\eta) - 1)}{\sqrt{(\xi^2 - 1)}} P'_{\lambda+N}(\xi_s(\eta)) \right) P'_{\lambda}(\eta) \quad (1.106)$$
A. Prosperetti [15] and discussed in E. Becker [27] and O.A. Basaran [28] and A.N. Zharov et al. [41]. We show in Fig.1.40 the damping of quadrupole surface oscillation amplitude computed with the VCVPF and RCVPF model and compare it to the exact damping [15].

<figure><img src="image_44.png" /><figcaption>Figure 1.40: For all panels, the x-axis is the dimensionless time $t/\tau_0$ and the y-axis is the dimensionless quadipole amplitude $a_2/R$. Panel (a): Quadupole surface oscillations with $Oh=0.001$. Pink full line uses the RCVPF model with $N=53$ while the blue dotted line uses RCVPF with $N=6$. Dashed-dotted orange line gives the damping of the VPF and the long-dashed blue line the damping in VCVPF approach. Green dashed line corresponds to the exact damping. Panel (b): Damping of quadupole oscillations with $Oh=0.01$. The RCVPF model with $N=18$ (dashed blue) is in agreement with the 'exact' damping (dotted pink). For subcritical $Oh$ values, the VCVPF overestimates the damping. Panel (c): Damping of quadupole moment with $Oh=0.1$. The RCVPF model with $N=4$ (full red line). Again, the VCVPF overestimates the damping. Panel (d): Damping of quadupole moment with the supercritical value $Oh=1$. The RCVPF model with $N=8$ (full red line) is in agreement with the 'exact' damping (dotted pink). VCVPF (green dashed line) underestimates now the damping.</figcaption></figure>

<footnote n="1">1. In panel (a) of Fig.1.40, we consider the low viscous limit with $Oh=0.001$.</footnote>
First we note that the damping of the (VCVPF) and (EXACT) methods are indistinguishable on the presented time range, confirming the validity of VCVPF in this regime. Second, we see that the RCVPF model with $N = 53 \simeq 1.68/\sqrt{Oh}$ (pink full line) is in agreement with the *exact* damping. Comparing the latter to RCPVF with $N = 6$ shows the sensitivity of the damping with respect to $N$ and thus with respect to the thickness of the boundary layer. Indeed, for $N = 6$, the damping is underestimated, as shown by the blue dotted oscillations, and closer to the damping of VPF model.

2. In panel (b) we considered $Oh=0.01$, which is still on the low viscous side and corresponds roughly to our water droplets at room temperature. The VCVPF model slightly overestimates the damping, as can be expected for sub-critical $Oh<Oh_c = 0.56$ values. On the other hand, the RCVPF model with $N = 18 \simeq 1.68/\sqrt{Oh}$ is in total agreement with the 'exact' damping.

3. In panel (c), we used $Oh=0.1$, which is still a sub-critical value. Damping of the quadupole moment in the RCVPF model with $N = 4 \sim 1.68/\sqrt{Oh}$ is in full agreement with the 'exact' damping. Again, the VCVPF overestimates the damping and the error is now significant.

4. Panel (d) shows the damping of the quadupole moment with $Oh=1$. The RCVPF model with $N = 8$ is in agreement with the 'exact' damping implying that the $"N \sim 1.68/\sqrt{Oh}"$ law may not apply for supercritical values of Oh. For **supercritical** values of Oh, the VCVPF now **underestimates** the damping.

We conclude this section by highlighting the excellent agreement of RCVPF with the exact damping in the case of small amplitude surface oscillations. Guided by the physical considerations, we proposed a simple function (1.106) for the rotational part of the velocity field which is parametrized by the integer $N$, accounting for the diffusion length of the vorticity layer of a viscous oscillating drop. At this time, the RCVPF model has not been used to calculate the whole deformation pathway of a critically charged drop but we are confident that RCVPF allows to go a step further than the VCVPF and simulate with even higher accuracy the flow of viscous drops with $Oh$ close to, or above the critical $Oh_c$. RCVPF has a comparable CPU cost advantage than VCVPF, but satisfies in addition the zero shear-stress boundary condition *adiabatically* by adding the necessary degrees of freedom in the form of a rotational velocity field. The latter adds a correction to the VCVPF dissipation term resulting in a damping rate that is indistinguishable from the exact one.
# Appendices
# Appendix A

## A.1 Numerical method of the Potential Flow (PF) model

A Galikin approach has been used by Patzek et al. for the numerical resolution of Bernoulli's equation in the case of *Nonlinear oscillations of inviscid free drops* [26]. We propose here a different numerical approach, making use of the fact the velocity potential and Coulomb potential can both be expanded on harmonic functions in a prolate spheroidal coordinate system.

As the velocity potential satisfies the Laplace equation, $\Psi$ can be written as a linear combination of harmonic functions, having the correct behaviour close to the origin. We choose to solve the Laplace equation in the prolate spheroidal coordinate system<fnref n="1" /> ($\xi, \eta, \phi$), which is well adapted to our boundary problem, while the harmonic functions $f_{nm}$ are still separable

$$f_{nm}(\xi, \eta, \phi) = P_n^m(\xi)P_n^m(\eta)e^{im\phi} \quad \text{in } \Omega \quad (A.1)$$

where $P_n^m$ is a Legendre polynomial<fnref n="2" /> of order $n$ and degree $m$. As $\xi > 1$ we use here the The angular coordinate $-1 < \eta_i < 1$, is discretised over 400 points, corresponding to the zeros of Legendre polynomial of order 400, $P_{400}(\eta_i) = 0$. Because of axial symmetry, only $m = 0$ need to be considered. The expansion coefficients $\{a_n\}$ are chosen such as to satisfy the boundary condition

$$\Psi(\eta_i) = \sum_{n} a_n P_n(\xi(\eta_i)) P_n(\eta_i) \quad \forall i \le 600. \quad (A.2)$$

$$= \sum_{n} a_n f_{n0}(\vec{s}_i) \quad (A.3)$$

<footnote n="1">¹https://en.wikipedia.org/wiki/Prolate_spheroidal_coordinates</footnote>

<footnote n="2">²https://en.wikipedia.org/wiki/Legendre_polynomials</footnote>
<figure><img src="image_45.png" /><figcaption>Figure A.1: Prolate coordinate system $(\xi, \eta, \phi)$. Domain of definition is ; $\xi \ge 1$, $-1 \le \eta \le 1$, $0 \le \phi \le 2\pi$. Shapes with constant $\xi$ are prolate spheroids of eccentricity $e = (\xi^2 - 1)^{-1/2}$. The 2 black points represent the focii and their distance $2a$ defines the length scale of coordinate system.</figcaption></figure>

The vector $\vec{s}_i = (\xi(\eta_i), \eta_i)$ describes the discretized droplet surface in the prolate coordinate system. As the boundary condition (1.68) evolves in time, the set of linear equations (A.2) must be solved at each time step. We use a singular value decomposition (SVD) of the matrix $f_n(\eta_i)$ to calculate the coefficients $\{a_n\}$. The velocity at the surface is finally obtained by taking the gradient of the velocity potential function $\Psi$ (1.43);

$$u_{\xi}(\eta) = \frac{1}{h_{\xi}} \sum_{n} a_{n} \frac{\partial P_{n}(\xi)}{\partial \xi} \Big|_{\xi(\eta)} P_{n}(\eta) \quad (A.4)$$

$$u_{\eta}(\eta) = \frac{1}{h_{\eta}} \sum_{n} a_{n} P_{n}(\xi(\eta)) \frac{\partial P_{n}(\eta)}{\partial \eta} \quad (A.5)$$

where $h_\xi$, $h_\eta$ and $h_\phi$ are scale functions of the prolate coordinate system [4]. A similar technique is used to calculate the Coulomb potential $V$ outside the droplet. For the expansion of $V$, we chose harmonic functions having the correct behaviour in the outer domain of the droplet,

$$g_{nm}(\xi, \eta, \phi) = Q_{n}^{m}(\xi) P_{n}^{m}(\eta) e^{im\phi} \quad \text{in } \mathbb{R}^{3}\setminus\Omega \quad (A.6)$$

where $Q_n^m$ is a Legendre polynomial of second type. For a perfect conductor, even in the presence of an external field $V_{\text{ext}}$, the potential at the surface is an equipo-
tential $V_0$. Assuming an axisymmetric external potential, $m = 0$, the expansion coefficients {$b_n$} are obtained by satisfying the boundary condition

$$V_0 - V_{\text{ext}}(\eta_i) = \sum_n b_n Q_n(\xi(\eta_i)) P_n(\eta_i) \quad \forall i \le 600 \quad (A.7)$$

$$= \sum_n b_n g_{n0}(\vec{s}_i) \quad (A.8)$$

Knowing the initial charge of the droplet $Q_0$, the value of $V_0$ is readily found by the equation

$$Q_0 = \varepsilon_0 \int_{\partial\Omega} \vec{E}^g \cdot \vec{n} \, ds \quad (A.9)$$

where the electric field outside the droplet is computed using

$$\vec{E}^g = - \sum_n b_n \left( \frac{1}{h_\xi} \frac{\partial g_{n0}(\xi, \eta)}{\partial \xi} \vec{e}_\xi - \frac{1}{h_\eta} \frac{\partial g_{n0}(\xi, \eta)}{\partial \eta} \vec{e}_\eta \right) \quad (A.10)$$

In the case where $V_{\text{ext}}$ is zero, the surface charge distribution at the equilibrium $\sigma = \varepsilon_0 \vec{n} \cdot \vec{E}^g$ depends thus only on the droplet shape $\vec{s} = (\xi(\eta), \eta)$. A schematic matrix representation of how to determine $V_0$ as a function of $Q_0$ is given in the footnote<fnref n="3" /> It allows to compute the capacity $C = Q_0/V_0$ of the droplet. We used a maximum of $n \le 200$ harmonic functions. Careful testing of the results convinced us that we were using a sufficient number of terms.

$$\sum_n G_{in} b_n = V_i^0 + V_i^{\text{ext}} \quad (A.11)$$

$$\sum_{j,n} M_{jn} b_n = Q_0 \quad (A.12)$$

The $b_n$ is eliminated in equation (A.12) by inverting equation (A.11).

$$\underbrace{\sum_{j,n,i} M_{jn} G_{ni}^{-1} V_i^0}_{C} = - \underbrace{\sum_{j,n,i} M_{jn} G_{ni}^{-1} V_i^{\text{ext}}}_{0} = Q_0 \quad (A.13)$$

The left hand term in equation (A.13) computes the capacity $C$ and gives the potential $V_0 = Q_0/C$.

<footnote n="3">³$V_0$ is the unknown and is computed knowing $Q_0$. We define $g_{n0}(\vec{s}_i) = G_{ni}$ and $\varepsilon_0 \left( \frac{1}{h_\xi} \frac{\partial g_{n0}(\vec{s}_j)}{\partial \xi} \vec{e}_\xi - \frac{1}{h_\eta} \frac{\partial g_{n0}(\vec{s}_j)}{\partial \eta} \vec{e}_\eta \right) \vec{n} ds = M_{jn}$. The sum over $j$ corresponds to the surface integration of (A.9). Relations (A.8) and (A.9) become respectively</footnote>