domain: info math

text:
Amaury Curiel

Antoine Genitrini

Lexicographic unranking algorithms for the Twelvefold Way *

Keywords: 2012 ACM Subject Classification Mathematics of computing → Combinatorial algorithms; Theory of computation → Generating random combinatorial structures Twelvefold Way, Set partitions, Unranking, Lexicographic order Digital Object Identifier 10.4230/LIPIcs

The Twelvefold Way represents Rota's classification, addressing the most fundamental enumeration problems and their associated combinatorial counting formulas. These distinct problems are connected to enumerating functions defined from a set of elements denoted by N into another one K. The counting solutions for the twelve problems are well known. We are interested in unranking algorithms. Such an algorithm is based on an underlying total order on the set of structures we aim at constructing. By taking the rank of an object, i.e. its number according to the total order, the algorithm outputs the structure itself after having built it. One famous total order is the lexicographic order: it is probably the one that is the most used by people when one wants to order things. While the counting solutions for Rota's classification have been known for years it is interesting to note that three among the problems have yet no lexicographic unranking algorithm. In this paper we aim at providing algorithms for the last three cases that remain without such algorithms. After presenting in detail the solution for set partitions associated with the famous Stirling numbers of the second kind, we explicitly explain how to adapt the algorithm for the two remaining cases. Additionally, we propose a detailed and fine-grained complexity analysis based on the number of bitwise arithmetic operations.

Introduction

The Twelvefold Way, a classification from the 1960s by Rota, was introduced to address the most fundamental enumeration problems associated with their combinatorial counting formulas. It has been extensively discussed in Stanley's book [START_REF] Stanley | Enumerative Combinatorics[END_REF]Section 1.9]. The distinct problems are related to the enumeration of functions defined from a set of elements denoted by N into another set denoted by K. The respective cardinalities of these sets are denoted as n and k. Each set may consist of either distinguishable or indistinguishable elements, resulting in consideration of four pairs of sets. Additional constraints pertain to the properties of the functions, whether they are injective, surjective, or arbitrary. Consequently, we encounter twelve cases when enumerating these functions. The counting solutions are well-known, as presented in Stanley's book [START_REF] Stanley | Enumerative Combinatorics[END_REF]Section 1.9]. In Table 1, we illustrate the classical combinatorial object enumerating each set of functions, in contrast to Stanley, who directly presents the counting solution.

In this paper, our focus lies in the generation of these classical combinatorial objects. To initiate our exploration, we arrange each object within a given class in lexicographic order. Subsequently, given the rank of an object, our goal is to construct it directly. This process is referred to as a lexicographic unranking algorithm. For instance, among the six permutations of {1, 2, 3}, the first one (with rank 0 in lexicographic order) is [START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF][START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF][START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF] • stand respectively for the Stirling numbers of the second kind and binomial coefficients; and p k (n) is the number of integer partitions of n into k positive integers.

by the second one (with rank 1), which is [START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF][START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF], and so forth, culminating with the last one (rank 5) being [START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF][START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF]. Consequently, the lexicographic unranking algorithm for the function with rank 4 returns [START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF][START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF]. In Table 1, we provide references to such algorithms for 9 out of the 12 cases. However, for cases 3, 7, and 9, no knowledge about lexicographic unranking algorithms seems available in the literature. This paper introduces an approach to unranking in lexicographic order for the set partitions of an n-set into k blocks (case 9). Furthermore, we present extensions of this approach to address cases 3 and 7.

The problem of unranking objects emerges as one of the most fundamental challenge in combinatorial generation, as seen in [START_REF] Shablya | Method for Developing Combinatorial Generation Algorithms Based on AND/OR Trees and Its Application[END_REF], and is applicable in various domains such as software testing [START_REF] Myers | k-out-of-n:G System Reliability With Imperfect Fault Coverage[END_REF], optimization [START_REF] Grebennik | Developing software for solving some combinatorial generation and optimization problems[END_REF], or scheduling [START_REF] Tamada | Parallel Algorithm for Learning Optimal Bayesian Network Structure[END_REF]. In different contexts, it serves as the core problem for generating complex structures, as observed in phylogenetics [START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF] and bioinformatics [START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF]. As mentioned earlier, to unrank, one must first establish a total order over the objects in question. The often-utilized order is the lexicographic order due to its ease of handling, leading to extensive study in the literature. However, Ruskey notes in [19, p. 59] that lexicographic generation is typically not the most efficient, thus requiring particular care in lexicographic unranking. Knuth dedicates a section to the lexicographic generation of combinatorial objects in [START_REF] Knuth | The Art of Computer Programming[END_REF], relating it to the special case of Gray codes.

Other combinatorial objects are also explored in Ruskey's and Kreher and Stinson's books on combinatorial generation [START_REF] Ruskey | Combinatorial Generation[END_REF][START_REF] Kreher | Combinatorial Algorithms: generation, enumeration, and search[END_REF]. Skiena focuses on the practical implementation of such algorithms [START_REF] Skiena | The Algorithm Design Manual, Third Edition[END_REF].

Usually, the approach for constructing structures using a recursive decomposition schema involves leveraging this decomposition to build a larger object from smaller ones. This method is extensively detailed in the well-known book by Nijenhuis and Wilf [START_REF] Nijenhuis | Combinatorial algorithms[END_REF]. The approach has been systematically applied to decomposable objects in the context of analytic combinatorics, initially for recursive generation [START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF], and later for unranking methods [START_REF] Martínez | A generic approach for the unranking of labeled combinatorial classes[END_REF]. Related work. Let us first quickly detail the classical unranking methods for the Twelvefold Way. As indicated in Table 1, cases 1, 8, and 11 are straightforward. In fact, an n-sequence in K consists of a word of length n over the finite alphabet K, making lexicographic unranking direct. Cases 8 and 11 are extreme situations, both corresponding to the Iverson bracket [n ≤ k]. As a result, the enumeration problems contain either one function (only when n ≤ k) or none. The unranking method is trivial.

Cases 4, 5, and 6 are all associated with the enumeration of subsets and are directly related to combination enumerations. Various algorithms to solve such lexicographic unranking problems are relevant in the literature. In [START_REF] Genitrini | Lexicographic Unranking of Combinations Revisited[END_REF], we present a survey of the most efficient methods with a modern algorithm complexity analysis. Moreover, we introduce a new algorithm based on the factoradics number system, which is at least as efficient as the others.

Case 12 is associated with integer partition enumerations, and [START_REF] Kreher | Combinatorial Algorithms: generation, enumeration, and search[END_REF] presents an efficient recursive algorithm. This algorithm follows lexicographic order but for the reverse standard form of printing a partition. In standard form, partitions print the components from the largest to the smallest, whereas this algorithm is based on the reverse printing (from the smallest component to the largest one). It appears that, currently, there is no existing lexicographic unranking method specifically designed for the standard form of printing. case 10 can be considered an extension of case 12, much like case 7 is an extension of case 9.

The last three cases pertain to set partition problems. Various combinatorial objects, such as permutations with a specific pattern [START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF], graph coloring [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF], walks in graphs [START_REF] Dzhumadil | Walks, partitions, and normal ordering[END_REF], or trees for phylogenetics [START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF], are enumerated by set partitions. In a recent paper [START_REF] Lutz | Active clustering for labeling training data[END_REF] the uniform random generating for set partitions for given n and k is studied, in the context of clustering algorithms. However, as far as we know, there is no lexicographic algorithm that takes arguments n, k and the rank r, returning the r-th partition in lexicographic order. Instead, there exists another classical object called a restricted growth sequence that is in bijection with set partitions (see [START_REF] Mansour | Combinatorics of Set Partitions[END_REF][START_REF] Ruskey | Combinatorial Generation[END_REF]). The unranking approaches presented in these works return such restricted growth sequences in lexicographic order. However, the natural bijection from restricted growth sequences to partitions does not preserve the lexicographic order. Main results. To develop an efficient unranking generator for set partitions, we first introduce the lexicographic order over set partitions. Some care must be taken since we are dealing with sets of integers. Therefore, we use a standard printing of a set partition to obtain a canonical representation. We then introduce an ad hoc combinatorial algorithm to unrank set partitions in lexicographic order. Due to the very large integers manipulated in the algorithms, of order of n ln n bits, our algorithm computes the necessary ones on-the-fly in a lazy paradigm. The correctness and complexity of the algorithm are managed based on specific combinatorial properties derived throughout the paper. Finally, we present some experiments using a Go 1 implementation for our algorithm. We leverage the simple and efficient parallelism mechanism provided by this language to significantly reduce space consumption without degrading time consumption for large values of n. Organization of the paper. Following the introduction of the paper, Section 2 highlights the combinatorial aspects of set partitions and presents some preliminary properties. In Section 3, we introduce our method for unranking set partitions, providing key insights into proving the correctness and complexity of our approach. Additionally, we present ideas for running calculations in parallel and share experiments that validate our parallel approach. Finally, Section 4 presents extensions of our unranking algorithm to address cases 3 and 7 from the Twelvefold Way. As an example let S be {1, α, 2, 3, 4, β, 6, 12}. The collection {2, 3}, {4, 6, 12}, {β, 1, α} is a partition of S in 3 blocks. In the rest of for paper, the set of positive integers from 1 to n is denoted by n . We can identify a set S of n elements with n , thus from now we will only be interested in partitions for n . ▶ Definition 2. Let 1 ≤ k ≤ n be two positive integers and S be n . The set of k-partitions of S is denoted by P n k . The sequential form of a partition of P n k (i.e. a k-partitions of S) is such that for all i ∈ k , the block B i contains the smallest integer from n not present in ∪ j<i B j . Furthermore for each block, it is represented in the increasing order of its elements.

Preliminaries

Context of set partitions

For example {1}, {2, 3, 5}, {4, 6} is a 3-partition of 6 represented in its sequential form. The sequential form is a canonical representation of the partition. As a shortcut, we will from now represent a partition simply as 1/235/46. In the paper we have chosen to use the terminology and the notations from Mansour [START_REF] Mansour | Combinatorics of Set Partitions[END_REF].

▶ Fact. Let 1 ≤ k ≤ n be two positive integers. The number of partitions in P n k is the Stirling number of the second kind denoted by n k . It satisfies the following recurrence:

n k = n-1 k-1 + k • n-1 k if 1 < k < n; 1 otherwise, i.e. when k = 1 or k = n. ( 1
)
This sequence is stored in OEIS A0082772 . We now introduce a natural order over k-partitions.

▶ Fact. Let A and B two subsets of positive integers. We say that A ≤ B iff either

A = B, or A ⊂ B and max(A) < min(B \ A), or B ⊂ A and min(A \ B) < max(B), or min(A \ B) < min(B \ A).

XX:5

The relation ≤ is a total order over subsets of n .

For example {1, 3} ≤ {1, 3, 4} and {1, 3} ≤ {1, 4}. But we also have {1, 3, 4} ≤ {1, 4}.

▶ Fact. Let 1 ≤ k ≤ n be two positive integers. The lexicographic order3 over partitions from P n k , in sequential form, is well defined using the latter order to compare two blocks: in fact a partition in k blocks is a Cartesian product of k subsets of positive integers.

There is another classical representation for partitions called canonical form in [START_REF] Mansour | Combinatorics of Set Partitions[END_REF]. A partition in k blocks is represented as a word over a k-letters alphabet. For example the partition 1/235/46 is represented by the word 122323. The ith letter is the index of the block containing the integer i. Using this representation we can also define a lexicographic order over partitions, but here we compare partitions that do not necessarily contain the same numbers of blocks. The lexicographic order over the sequential form is not compatible with the lexicographic order used for the sequentical form we are interested in. This can be noted in the Table 2.

▶ Definition 3. Using the lexicographic order over the sequential form for partitions in P n k , we define a ranking function assigning to each partition its rank corresponding to the number of k-partitions smaller than it in the lexicographic order.

Rank Partition Canonical form [START_REF] Mansour | Combinatorics of Set Partitions[END_REF] ▶ Definition 4. Let 1 ≤ k ≤ n be two positive integers. Let P be a partition from P n k , represented in the sequential form as

B 1 /B 2 / . . . /B k . An integer subset p is called prefix of P if p ⊂ B 1 and p ≤ B 1 .
For the partition 12/35/4, there are three possible prefixes ∅, 1 and 12. We can further extend the definition of prefixes of a partition by letting S being any subset of n . Thus removing the first block of the latter partition gives 35/4, we define prefixes of the 2-partition (of {3, 4, 5}) to be ∅, 3 and 35. Here we formalize this extension.

▶ Definition 5. The definition of a prefix p of a partition is extended to any set S partitioned in a sequence of blocks (with the first one being denoted by B 1 ) such that p ≤ B 1 .

Combinatorial properties

We are now interested in counting results for partitions sharing the same prefix. These are the core results for our unranking algorithm.

▶ Proposition 6. Let 1 ≤ k ≤ n be two positive integers. Let ℓ and d be two integers such that either ℓ = d = 1 or 1 < ℓ ≤ d. Let S n k (ℓ, d) be the number of partitions in P n k accepting as prefix of length ℓ: 1 α 2 α 3 . . . α ℓ-1 d. We have

S n k (ℓ, d) = min(n-k-ℓ+1,n-d) u=0 n -ℓ -u k -1 n -d u .
Proof. First if ℓ = d = 1, then in the sequential form the first block necessarily contains 1.

Thus S n k (1, 1) = |P n k | = n k = n-k u=0 n-1-u k-1 n-1
u . The latter equality is given e.g. in [7, p. 251, Table 251].

In the second case when 1 < ℓ ≤ d, we aim at counting the number of partitions in P n k accepting 1 α 2 α 3 . . . α ℓ-1 d as a prefix. In order to exhibit a combinatorial interpretation, we rewrite S n k (ℓ, d) as

S n k (ℓ, d) = min((n-ℓ)-(k-1),n-d) u=0 n -(ℓ + u) k -1 n -d u .
Once the prefix is given, it remains to complete the first block B 1 from the partition, and then to calculate how we can further partition the other elements in the next blocks. The variable u in the sum corresponds to the number of elements that are appended to the prefix to complete B 1 . Its value ranges from 0 up to the maximal number of elements that we can append i.e. (n -ℓ) -(k -1) because at least k -1 among the remaining n -ℓ elements must be assigned to the other k -1 blocks. Obviously the number of possible elements u is also upper bounded by the number of remaining elements, i.e. n -d. Once the number u of elements for the completion of B 1 is given, we choose u elements greater than d: the number of possibilities is given by the binomial coefficient. Finally it remains to build the other blocks of the partition: we partition n -(ℓ + u) elements into k -1 blocks. Hence the formula is proved. ◀

We introduce an example using Table 2 for P 5 3 . If we are interested in the prefix 13, then there are 3 partitions without completing block B 1 , in the sum, when u = 0 we get 3 2 = 3. The other possible value is u = 1 with the general term being 2 In order to get a formula that is more efficient to calculate, we observe that the latter numbers S n k (ℓ, d) depend essentially in three variables instead of four. The proof is direct with some variable renaming.

▶ Proposition 7. Let n, k, d be integers with 0 ≤ k ≤ n and 0 ≤ d ≤ n. By defining Sn k (d) = min(n-k,n-d) u=0 n -u k n -d u , we get S n k (ℓ, d) = Sn-ℓ k-1 (d -ℓ).
We note that S n k (1, 1) = Sn-1 k-1 (0) = n k . Note that the 3-dimension sequence S seems not to be stored in OEIS. There exits several generalizations of Stirling numbers, but none of them apparently corresponds to our sequence S.

▶ Corollary 8. The numbers Sn

k (d) satisfy the following recurrence:

Sn k (d) =        Sn-1 k-1 (d -1) + k • Sn-1 k (d -1) if 1 ≤ k ≤ n and 1 ≤ d ≤ n; n+1 k+1 if d = 0 and 0 ≤ k ≤ n; 0 otherwise. (2)
Note the later recurrence is similar to the one satisfied by Stirling numbers of the second kind (but with here a third variable d giving some kind of level of numbers). The proof is provided in Appendix A.

▶ Proposition 9. Let n, k, d be integers with 0 ≤ k ≤ n and 0 ≤ d ≤ n. The function Sn k (d) can be represented as a binomial transform:

Sn k (d) = min(n-k,d) u=0 (-1) u n + 1 -u k + 1 d u .
The main idea of the proof consists in proving that the two expressions given in Propositions 7 and 9 are satisfying the same recurrence and thus are equal.

Proof. In order to prove this new expression for S, we just have to prove that this expression satisfy the recurrence stated in Corollary 8. Substituting d by 0 we get the base case. We now consider the case where the three integers n, k, d satisfy 0 ≤ k ≤ n and 1 ≤ d ≤ n. Using Proposition 9 in the case where 0 < k < n (the cases k = 0 or k = n are obvious) we have

Sn-1 k-1 (d -1) + k • Sn-1 k (d -1) = min(n-k,d-1) u=0 (-1) u n -u k d -1 u + k • min(n-1-k,d-1) u=0 (-1) u n -u k + 1 d -1 u
By using factorization and Stirling numbers of the second kind recurrence, we obtain:

Sn-1 k-1 (d -1) + k • Sn-1 k (d -1) = min(n-k,d-1) u=0 (-1) u n + 1 -u k + 1 - n -u k + 1 d -1 u .
After having telescoped the two sums we get the stated result. ◀

Finally, given two prefixes, one being smaller than the second one, the next proposition allows to compute how many partitions are in-between the two prefixes. More formally:

▶ Proposition 10. Let 1 ≤ k ≤ n be two positive integers. Let d 1 ∈ n \ {1}, d 0 ∈ [d 1 -1]
and ℓ > 1 be integers. The number of elements of P n k that admit a length-ℓ prefix satisfying 1 α 2 . . . α ℓ-2 d 0 d1 (for all d1 ranging from d 0 + 1 to d 1 ) is given by

R n k (ℓ, d 0 , d 1 ) = Sn-ℓ k-1 (d 0 -ℓ) -Sn-ℓ k-1 (d 1 + 1 -ℓ).
The proof of the proposition is given in Appendix A.

Methods for unranking set partitions

Merging the combinatorial properties stated in the previous section, we are now ready to design algorithms to unrank set partitions in the lexicographic order.

Unranking algorithm design

Our aim consists in constructing the r-th partition related to a pair n, k in sequential form for the lexicographic order. The constructions follows the next main lines. The global idea consists in building the normalization of the partition. So we build together its block pattern and its reversed factoradics (seing the partition as a size-n permutation).

The building of the blocks is going from left to right;

The construction of a block is also from left to right, component by component using a binary search approach; Finally once the block pattern and the reversed factoradics are set, a slight adaptation of the lexicographic permutation unranking algorithm gives the result.

The details for the correctness of our approach lies on the ranking function associated to the set partitions. Details about this function are presented in Appendix B.

We first present in detail the main function Unranking of Algorithm 1. Using a loop, at each turn it defines the next block of the partition and then refine the value of the rank related to the remaining part of the partition. The result B returned by next_block contains the indices of the components of the block that has been calculated and acc allows to update the rank so that it is related to the remaining part of the partition that must still be computed. With our previous definition, B is the normalization of the corresponding partition block. At the end of the function a dynamic extraction is executed in an array containing elements from 1 to n according to the indices in Res. n ′ := n 3:

Res := [] 4:

while k > 1 do 5:

(B, acc) := next_block(n, k, r) 6:

Append(Res, B) 7:

r := r -acc 8:

n := n -len(B) 9:

k := k -1 10:

Append(Res, [0, 0, . . . , 0]) 11:

Res :=Extract(n ′ , Res) 12:

return Res 1: function Extract(n, R) 2:

L := [1, 2, . . . , n] 3:

P := [] 4:
for r in R do 5:

p := [] 6:
for i in r do 7:

Append(p, L[i]) 8:

Remove(L, i) 9:

Append(P, p) 10:

return P Remove(L, i) removes element with index i in L.

1: function next_block(n, k, r) 2:

Block := [0]; acc := n-1 k-1 3: if r < acc then 4:
return (Block, 0) 5: d0 := 1; index := 2; inf := 2; sup := n 6:

complete := F alse 7:

while not complete do 8:

while inf < sup do 9: The function next_block takes parameters n, k, r and returns essentially the first block of the r-th partition in P n k . In fact, using Table 2 the call next_block(5, 3, 16) returns 0 1 1 (instead of 1 3 4), the latter block being obtained through a dynamic extraction of the element 0 in [START_REF] Ali | Peptide Combination Generator: a Tool for Generating Peptide Combinations[END_REF][START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF][START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Dzhumadil | Walks, partitions, and normal ordering[END_REF][START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF] then the element 1 is extracted in the remaining part [START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF][START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Dzhumadil | Walks, partitions, and normal ordering[END_REF][START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF] and finally the element 1 in [START_REF] Bodini | Strict monotonic trees arising from evolutionary processes: Combinatorial and probabilistic study[END_REF][START_REF] Dzhumadil | Walks, partitions, and normal ordering[END_REF][START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF]. Constructing the blocks of indices instead of the blocks of values allows to neglect about the remaining elements for the further blocks construction. Note that obviously the last block of the partition contains only the indices 0 (Line 10 from Unranking Algorithm) and the first element of a block is always index 0, both due to the sequential form. Finally while calling Unranking [START_REF] Flajolet | A calculus for the random generation of labelled combinatorial structures[END_REF][START_REF] Defant | Highly sorted permutations and Bell numbers[END_REF][START_REF] Myers | k-out-of-n:G System Reliability With Imperfect Fault Coverage[END_REF], at the end of Line 10, Res contains

mid := ⌊(inf + sup)/2⌋ 10: if r >= acc + R n k (index -1, d0,
[[0 1 1][0][0]].
Reading the components from right to left we get the factoradics 0 0 1 1 0 of the number 8 corresponding to the lexicographic rank of the permutation [1, 3, 4, 2, 5] (cf. [START_REF] Genitrini | Lexicographic Unranking of Combinations Revisited[END_REF] for details).

▶ Theorem 11. Unranking(n, k, •) is a lexicographic unranking algorithm for set partitions from P n k .

Proof key-ideas. To derive the proof, we first rely on the factoradics of the rank number, and the Ranking function, both detailed in Appendix B. Hence the core property of our algorithm relies on the function next_block. Its call with parameters n, k, r returns the first block of the r-th partition of P n k (in sequential form). After k -1 calls to it we get the pattern of the partition and also the factoradics for the underlying permutation. ◀

Complexity analysis and experiments for unranking

In our implementation in Go4 , we offer two approaches for the necessary Stirling numbers calculations: either a precomputation of them or a computation on the fly of those that are needed at each step. We never precompute the 3 dimension table Sn k (d). In fact, in many bad cases these numbers are of order of n!, thus precomputing would be too expensive while only few of the numbers are needed. We compute the necessary numbers Sn k (d) on the fly. First let us recall the behavior of the sequence of Stirling numbers of the second kind when k is ranging from 1 to n.

▶ Fact. Let 1 ≤ k ≤ n be two positive integers. The sequence n k is unimodal and its maximum is reached when k n ∼ n/ ln n. Around this value, we have log n kn = Θ(n log n) 5 . Furthermore, we have an upper bound valid for all 1 ≤ k ≤ n:

log n k ≤ (n -k) log k + log n k ≤ (n -k) log k + k log n • e k .
See the fundamental paper of Rennie and Dobson [START_REF] Rennie | On stirling numbers of the second kind[END_REF] to get a proof for these results.

In the following we propose six distinct implementations of the function S presented in Proposition 7 and underlying the function R from Proposition 10. S_v1: direct implementation of the formula stated in Proposition 7; S_v2: implementation of the formula from Proposition 7 taking into account the symmetry of binomial coefficients, thus the sum contains half of the terms in comparison to S_v1 (and thus half of the multiplications); S_v3: direct implementation of the formula stated in Proposition 9; S_v4: implementation of the formula from Proposition 9 taking into account the symmetry of binomial coefficients, thus the sum contains half of the terms in comparison to S_v3; S_v5: is our most efficient algorithm without precomputations. The calculation way consists in deciding according whether a call to S_v2 or to S_v4 should be the most efficient, according to the number of terms in the sums interacting with Propositions 7 and 9; S_v6: all necessary Stirling numbers of the second kind are precomputed and this precomputation is not taken into account in the further time complexity experiments.

XX:10 Lexicographic unranking for the Twelvefold Way

The integers computed during the unranking algorithm are very large, thus a classical complexity in the number of arithmetical operations is not precise. We hence are interested in the bit-complexity, corresponding the the number of atomic operations on digits.

▶ Theorem 12. For the time complexity, the algorithm Unranking based on each of the function S_v• has a bit-complexity bounded by

O (n -k) 3 M (n) n ln n log k + k(n -k) 2 M (n) n ln n log n • e k ,
where M (n) is the bit complexity for the multiplication of two numbers, each one containing n bits.

The naïve multiplication algorithm satisfies M (n) = Θ(n 2 ). But using, for example, Karatsuba algorithm, we obtain M (n) = Θ(n log 3 ) for the time complexity. In Go6 , as soon as the integers are greater than 2 40 , Karatsuba multiplication algorithm is used. In our context, almost all cases are thus based on the latter algorithm.

Proof. We are interested in a worst case complexity analysis when n is large and for k ranging in n . We are using the same kind of analysis in bit complexity as the one presented in [START_REF] Genitrini | Lexicographic Unranking of Combinations Revisited[END_REF]Section 4.3]. We compute an upper bound of the complexity in the central range of the Stirling numbers of the second kind. The central range, when n tends to infinity, is observed when k = Θ(n/ ln n). A detailed similar analysis is presented in the paper [START_REF] Louchard | Asymptotics of the Stirling numbers of the second kind revisited[END_REF]. In our context each Stirling number necessitates log n k bits to be stored. They are multiplied by binomial coefficients containing at most n bits. Thus Stirling numbers are separated in blocks of n bits in order to use a multiplication of similar sizes numbers, inducing a time complexity bounded by log n k /n • M (n). Furthermore the number of calls the the function S is O((n -k) ln n) induced by the repetitive calls to the binary search algorithm. Compiling all these upper bounds gives the stated bit-complexity. ◀

For approach S_v6, the following result establishes that the precomputation is negligible in terms of time complexity compared to the unranking itself. However note that the memory complexity is quadratic instead of linear (in n) by using the precomputation step.

▶ Proposition 13. The bit-complexity for the Stirling numbers precomputation is bounded by

O k(n -k) 2 M (log k) + k 2 (n -k) log n log k M (log k) .
In order to get the Stirling numbers on the fly, we use parallel compuatations. In fact, for each block determination, we observe that only two neighbors columns from the triangle of numbers are needed. Thus during the determination of a block, we compute in parallel the next two columns that will be necessary for the next block. Thanks to this approach, we present an algorithm which uses a reasonable quantity O(n) of memory and with essentially the same time complexity than the algorithm where Stirling numbers are exhaustively precomputed, thus needing O(n 2 ) memory size.

In Figure 1, we run experiments7 by fixing n = 1000 and k ranging from 2 to 992 with steps of 15 units. For each value of k, 500 uniform samples are computed and the average time for the building of the partition is drawn for each Algorithm S_v1 up to S_v6. Obviously Algorithms S_v2 and S_v4 are better than their naïve versions respectively S_v1 and S_v2. It is interesting to note that the optimization S_v5, obtained by computing the most efficient formula between Propositions 7 and 9. Finally we remark that the Algorithm S_v5 is almost as efficient as S_v6 where all precomputation of Stirling numbers have been stored before the computation of the partition. Strangely, for the smallest values of k, we note that S_v5 is even faster than S_v6. This is probably due to the RAM accesses: in fact in some preliminary experiments with computers equipped with DDR5 RAM Algorithm S_v6 is always faster than S_v5, and this is what is expected.

Extension and conclusion

As we observe in Table 1, both enumeration cases 3 and 7 are some extended version of the enumeration case 9. An adaptation for the Ranking function allows to rank the families counted by cases 3 and 7; then adapting the unranking algorithm solves these cases.

Ordered set partitions

Recall Stirling numbers of the second kind are counting the numbers of surjective functions from set N to set K, where the elements of N are distinguishable and also for those in K.

We can represent these functions as set partitions. Now, what happens when elements of K are distinguished? These functions are counted by ordered Stirling numbers of the second kind. In addition, they can be represented as ordered set partitions, which are similar to set partitions except that the order of the subsets matters. For instance, while in the world of unordered set partitions, elements 14/25/3; 14/3/25; 25/14/3; 25/3/14; 3/14/25 and 3/25/14 are equivalent and represented by the partition 14/25/3 in sequential form, in the world of ordered set partition, the 6 elements are all different.

▶ Proposition 14. Let 1 ≤ k ≤ n, be two integers with n being the cardinality of set N . The number of ordered set partition of N in k (non empty disjoint) subsets is k! • n k . The family of these partition is denoted by O n k .

The proof is direct: the blocks in the sequential form of a set partition are distinguishable, thus permuting them gives the associated ordered set partitions.

▶ Fact. Let 1 ≤ k ≤ n, be two integers, the lexicographic order on set partitions P n k is easily extended to get the lexicographic order for the ordered set partitions from O n k .

We can now derive the enumeration core result in this new context. T

n k (ℓ, d) = min(n-k-ℓ+1,n-d) u=0 k! n -ℓ -u k -1 n -d u .
This formula is the analogous to S n k (ℓ, d). Using the same variable changes, we also get a three variable function, like Sn k (d). Then we can deduce an adaptation of our first algorithm by replacing Stirling numbers of the second kind by ordered Stirling numbers of the second kind and using the latter formula.

Bell's set partitions

In the twelvefold way, the remaining case where no lexicographic unranking algorithm is known is when f is an arbitrary function from a set N whose elements are distinguishable to a set K where elements are indistinguishable. We denote by F the family of these functions. Such functions can be represented as unordered set partitions with at most k blocks where k is the numbers of elements in K. Let K i ⊂ K be a subset of i distinguishable elements and B i the functions that are surjective from N into K i . We have

B n k = k i=1 B i and for a given i ∈ k , |B i | = n i . Obviously |B n k | = k i=0 |B i | = k i=0 n i .
U n k (ℓ, d) = min(n-k-ℓ+1,n-d) u=0 k i=1 n -ℓ -u i -1 n -d u .
Again, his formula is the analogous to S n k (ℓ, d). Using the same variable changes, we also get a three variable function, like Sn k (d). Then we can deduce an adaptation of our first algorithm by replacing Stirling numbers of the second kind by Bell's number and using the latter formula.

As a final remark, the correctness of both previous algorithms is directly hanging to the one for the set partition algorithm. What is remaining is their complexity analysis: it is not difficult, and it will be written in a long version of this paper.

B Appendix related to the ranking algorithm

In order to introduce the way to rank a set partition, we first recall the definition of the factorial number system, or factoradics. It is a mixed radix numeral system in which the representation of integers relies on the use of factorial numbers. The factoradics are deeply related to the lexicographic (un)ranking of permutations. Details and references are presented in our paper [START_REF] Genitrini | Lexicographic Unranking of Combinations Revisited[END_REF].

▶ Definition 17. Let u be a positive integer and let n be the unique integer satisfying (n-1)! ≤ u < n!. Then there exists a unique sequence of integers (f ℓ ) ℓ∈{0,...,n-1} , with 0 ≤ f ℓ ≤ ℓ for all ℓ, such that:

u = f 0 • 0! + f 1 • 1! + • • • + f n-2 • (n -2)! + f n-1 • (n -1)!
The finite sequence f 0 f 1 . . . f n-1 is called the factoradic decomposition (or factoradics) of u (note that f 0 = 0 for all values for u).

Take the number u = 41 196 as an example, we obtain the following decomposition: 41 196 = 0

• 0! + 0 • 1! + 0 • 2! + 2 • 3! + 1 • 4! + 1 • 5! + 1 • 6! + 0 • 7! + 1 • 8!
, thus its factoradics is 0 0 0 2 1 1 1 0 1. Obviously 0s can be appended to the right of the factoradics without changing the integer.

▶ Fact. Let n be a positive integer. The factoradics with n components relates the lexicographic rank of a permutation and the permutation itself.

This deep link is, for example, described in our previous paper [START_REF] Genitrini | Lexicographic Unranking of Combinations Revisited[END_REF].

Let us further introduce two new notations. Let 1 ≤ k ≤ n be two positive integers. Let P be a partition from P n k , represented in its sequential form B 1 /B 2 / . . . /B k . For i ∈ k , we denote by t i the number of elements of the block B i and for j ∈ [t i ], the notation B i,j refers to the j-th element of B i . We further denote by t <i the accumulated number of elements in B 1 / . . . /B i-1 , i. e. t <i = i-1 j=1 t j .

▶ Definition 18. Let 1 ≤ k ≤ n be two positive integers. Let P ∈ P n k be a partition written in sequential form. The normalization of P , denoted by P consists in the collection B1 / B2 / . . . / Bk , such that each element Bi,j is equal to B i,j translated by -|{ℓ | ∃i 0 ≤ i and B i0,ℓ ≤ B i,j }|.

For example, for P the set partition from P 10 3 written in sequential form : P = 13/25679/4810, we get P = 01/01112/000. Using the normalization of a set partition is such that, once the first blocks have been computed, the rest of the partition can be constructed without considering the elements appearing in the previous blocks. We gain some kind of independence between the blocks. Obviously this is only a relative independence. Once we want to get the set partition itself, we take care of the remaining elements when the first blocks have been determined.

▶ Theorem 19. Let 1 ≤ k ≤ n be two positive integers. Let P ∈ P n k be a partition written in sequential form and P = B1 / B2 / . . . / Bk its normalization. The rank of P for the lexicographic order is given by

Ranking(P ) = k-1 i=1 ti ℓ=2 n -t <i -ℓ + 1 k -i + R n-t<i k-i+1
(ℓ, Bi,ℓ-1 + ℓ -1, Bi,ℓ + ℓ -1) .

XX:16 Lexicographic unranking for the Twelvefold Way

Proof key-ideas. Let 1 ≤ k ≤ n be two positive integers. Let P ∈ P n k be a partition and P its normalization written as B1 / B2 / . . . / Bk . The Ranking function takes into account the block pattern of P and the rank of of P seen as a size n permutation (i.e. by omitting the block pattern of P ) to produce the lexicographic rank of P seen à a partition from P n k . First of all, from P (or P ) we deduce the values for n and k. Then, in order to compute the rank associated of P , we iterate the following process from left to right, block by block. We do not need to consider the remaining values once a block has been built due to the normalization P .

So let us suppose we are interested in the first block. We iterate over all prefixes p smaller than B1 , and accumulate the numbers of partitions accepting p as a prefix, because all of them are smaller than P in the lexicographic order. These numbers are calculated using functions S n k and R n k . Once done we get the first possible value for the rank of P , and also its maximal value again by using S n k (the upper bound is not necessary for the rest of the process). We the start again with the next block, keeping in memory the accumulated value. ◀

Let us present on our running example the process. So we are dealing with P whose normalization is 01/01112/000. We thus deduce n = 10 and k = 3. The accumulate rank is 0 for the moment. Starting with the first block 01 we increase the accumulate rank by 255 corresponding to the number of set partition whose first normalized block is reduced to 0, i. e. 9 2 . We further increase the rank by 3025 that corresponds to all set partitions whose first block is of length at least 2 and that accept 0 as a prefix, but not 01. At this point, we have used all information from the first normalized block and memorized it into the accumulated rank whose value is now 3280. The second block starts with 0 and the number of set partitions with a block reduced to 0 is now 7 1 = 1, i. e. it remains 7 elements to fill the last block. Then the number of partitions whose (first) normalized block is of length at least 2 and that accept 0 as a prefix, but not 01 is 63. Thus the accumulated rank is now 3344. At this point we suppose the the second block starts with 01 and first it is reduced to 01, there is only 1 in this situation (all other elements fill the last block). We now are interested in the number of partitions whose (first) normalized block is of length at least 3 and that accept 01 as a prefix, but not 011. Their count is 0 (010 is not an acceptable prefix). Thus the accumulated rank is now 3345. It remains two steps as the previous ones to finally obtain the rank 3351 for the set partition P .

Finally, as a conclusion, taking the lexicographic rank of a partition, the values n and k, and following an analogous process as the previous one (but not the reversed process because it is more convenient to do the exhibition of the normalized partition from left to right) we obtain an algorithm for lexicographic unranking.

2 21

 2 = 2 as it appears in the table.

Algorithm 1 1 :

 11 Lexicographic unranking of the partition with rank r in P n k function Unranking(n, k, r) 2:

Figure 1

 1 Figure 1 Time (in seconds) for unranking a partition in P 1000 k when k is ranging in 1000 .

▶ Proposition 15 .

 15 Let 1 ≤ k ≤ n be two integers. Let ℓ and d be two integers such that either ℓ = d = 1 or 1 < ℓ ≤ d. Let T n k (ℓ, d) be the number of ordered set partitions in O n k accepting of the prefix of length ℓ: α 1 . . . α ℓ-1 d.

, followed XX:2 Lexicographic unranking for the Twelvefold Way

	elts of N	elts of K	f is arbitrary	f is injective	f is surjective
			1.		2.	3.
	dist.	dist.	n-sequence in K	n-permutation of K	composition of N
						with k subsets
	enumeration	k n		k n	k! • n k
	lex. unranking	easy		[6, Section 5]	Section 4.1
			4.		5.	6.
	indist.	dist.	n-multisubset of K	n-subset of K	composition of n
						with k terms
	enumeration	k+n-1 n	k n	n-1 n-k
	lex. unranking		see survey [6] and references therein
			7.		8.	9.
	dist.	indist.	partition of N	partition of N	partition of N
			into ≤ k subsets	into ≤ k elements	into k subsets
	enumeration	k i=0	n i	[n ≤ k]	n k
	lex. unranking	Section 4.2	easy	Theorem 11
			10.		11.	12.
	indist.	dist.	partition of n	partition of n	partition of n
			into ≤ k parts	into ≤ k parts {1}	into k parts
	enumeration	p k (n + k)	[n ≤ k]	p k (n)
	lex. unranking	[19, Section 4.8]	easy	[19, Section 4.8]

Table 1

 1 The Twelvefold Way *

* The notation k n corresponds to the product k • (k -1) • • • (k -n + 1); [n ≤ k] is the Iverson bracket returning 1 when n ≤ k and 0 otherwise;

  Let S be a set of n distinguishable elements. A partition π of S in k blocks is a collection B 1 , B 2 , . . . , B k of disjoint non-empty subsets of S such that for every element from S belongs to exactly one B i , for i ∈ {1, . . . , k}.

▶ Definition 1.

Table 2

 2 Ranking of the 3-partitions of 5

	Rank Partition Canonical form [14]

  The cardinality of B n k is counted by the k-restricted Bell numbers and finally, when k = n, we get the Bell numbers. Let 1 ≤ k ≤ n be two integers. Let ℓ and d be two integers such that either ℓ = d = 1 or 1 < ℓ ≤ d. Let U n k (ℓ, d) be the number of Bell's set partitions in B n

	▶ Fact. Let 1 ≤ k ≤ n, be two integers, the lexicographic order on set partitions P n k is also
	lexicographic for B n k .

▶ Proposition 16.

k accepting of the prefix of length ℓ: 1 α 2 . . . α ℓ-1 d.

The Go language offers routines to manage concurrency.

OEIS stands for the On-line Encyclopedia of Integer Sequences.

The lexicographic order of partitions from P n k in sequential form is a total order.

Go implementation and the material used for repeating the experiments are all available here.

In this paper we use the notation log for the logarithm in basis 2.

Go documentation for big integers manipulations.

The experiments provided in this paper are driven by using a PC equipped with an Intel Xeon X5677 processor, 32GB of DDR4-SDRAM and running Debian GNU/Linux 12.

This work was supported by the anr-fwf project PAnDAG ANR-23-CE48-0014-01.

A Appendix related to the formulas for S

Here we provide the proof of the recursive equation satisfied by Sn k (d).

Proof of Corollary 8. We aim at proving that this recurrence is satisfied by the formula Sn k (d) given in Proposition 7. Let 0 ≤ k ≤ n be two integers. We get the value

k+1 according to [7, p. 251]. We now consider the case where the three integers n, k, d satisfy 0 ≤ k ≤ n and 1 ≤ d ≤ n. Using Proposition 7 in the case where 0 < k < n (the cases k = 0 or k = n are obvious) we have

Using the classical recurrence for Stirling numbers of the second kind (cf. Equation ( 1)) we obtain the stated result. ◀ Proof of Proposition 10. Using the notations from the proposition statement, we are interested in the number of set partitions from P n k accepting as prefix 1 α 2 . . . α ℓ-2 d 0 d1 , when d1 ranges in {d 0 + 1, . . . , d 1 }. So using the approach exhibited in the proof of Proposition 6, we get

In fact the first sum computes the number of partition accepting prefix 1 α 2 . . . α ℓ-2 d 0 , but the first block contains at least one more element (since u starts at 1 ). To those partitions we remove the ones whose first block is completed only with elements greater than d 1 . Then if we add the terms u = 0 in both sums, this terms are opposite, thus we can rewrite

as stated in Proposition 10. ◀

C Appendix related to the unranking algorithm

Let us introduce a second experiments in Figure 2 for larger set partitions. We do not observe here anymore Algorithm S_v5 that is more efficient than S_v6. When we compute the ratio of time necessary for both algorithms, we note when n is increasing the ratio is tending to 1 until when k is almost n (then there is almost nothing to compute for obtaining a set partition while the computing of a linear number of Stirling numbers of the second is costly and present using S_v5 algorithm).


E Giglio

J Rangama

S Guillous

T Le Cornu

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

Keywords: numbers: 47, 20, Dr, 47, 20, Ma

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

I. INTRODUCTION

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

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

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

II. DYNAMICS OF THE DROPLET DEFORMATION

We want to simulate the nonlinear time-evolution of the droplet shape that develops after the onset of the Coulomb instability of an initially spherical charged drop of radius R in a gas g atmosphere. A liquid of dynamic viscosity µ and density ρ is supposed electrically conductive. As a result, the charges accumulate at the droplet surface and the electric field is divergence-free ∇ • E = 0 inside and outside the droplet volume. This corresponds to the Taylor-Melcher leaky dielectric model [18,[START_REF] Saville | The Taylor-Melcher Leaky Dielectric Model[END_REF]. The droplet deformation is assumed to be axisymmetric and the liquid of the droplet incompressible. The latter condition implies that the velocity vector field of the fluid is divergence-free (solenoidal),

∇ • u = 0 (1)
ensuring that the time evolution of the droplet occurs at constant volume. The Navier-Stokes equation for an incompressible (Newtonian) fluid in the droplet volume Ω is expressed by the equation

d u dt = ∂ u ∂t + u • ∇ u = 1 ρ ∇ • T ( 2
)
where T is the stress tensor, defined by

T = -pI + 2µD + M (3)
with I being the identity matrix, p the pressure in the liquid, D the strain-rate tensor,

D = 1 2 ∇ u + ( ∇ u) t (4)
and M the Maxwell stress tensor

M i,j = ε 0 ε r E i E j - ε 0 ε r 2 E 2 δ i,j , (5)
with ε r being the dielectric constant of the fluid. Note that T is symmetric, T = T t . As there are no volume charges, the electric force per unit volume is zero, ∇•M = 0. Thus, in the leaky dielectric model, the Maxwell tensor is divergence-free.

Boundary conditions for viscous charged droplets

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

T nn || = n • (T g -T ) • n = T g nn -T nn = P γ ( 6
)
The latter is due to the phenomenon of surface tension γ, and is given by the Young-Laplace equation

P γ = γ( ∇ • n) = -2γH , (7)
where H is the mean curvature at a surface point and γ is considered constant. Having defined the normal vector n at the drop interface pointing outward, the mean curvature of a sphere of radius R is negative, H = -1/R, so that P γ > 0. This definition ensures that the pressure in the liquid is higher than the pressure of the atmosphere. We suppose that the viscosity of the outer gas atmosphere is negligible µ g µ and that the pressure of the outer atmosphere p g is constant (and set to zero), so at the drop interface ∂Ω one has

P γ = T g nn -T nn = -(p g -p) + 2(µ g -µ)D nn + M g nn -M nn p -2µD nn + M g nn -M nn (8)
We introduce the Coulomb pressure P c , defined as the jump of the normal component of the Maxwell tensor across the interface

P c = M g nn -M nn = ε 0 2 (E g n ) 2 -ε r E 2 n + (ε r -1)E 2 τ ( 9
)
where τ stands for the tangential vector at the surface along the polar angle. Injecting the surface charge density σ = ε 0 (E g n -ε r E n ) into (9) allows expressing the Coulomb pressure as a function of σ and the electric field (E n , E τ ), evaluated at the liquid side of the liquid-gas interface of the droplet.

P c = ε 0 2 σ ε 0 + ε r E n 2 -ε r E 2 n + (ε r -1)E 2 τ ( 10
)
Finally rearranging (8) we obtain an equation for the liquid pressure at the interface,

p = P γ -P c + 2µD nn . (11)
The shear stress T τ n , in contrast, is continuous across the free interface,

T τ n || = τ • (T g -T ) • n = 0 (12) = 2µ g D g τ n + M g τ n -(2µD τ n + M τ n ) = 0 . ( 13
)
We deduce from [START_REF] Nganguia | Electrohydrodynamics of a viscous drop with inertia[END_REF], by considering the viscous shear stress at the gas side of the interface to be negligible, that the viscous shear stress in the liquid must be equal to the Maxwell shear stress across the interface,

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

d s dt -u( s) • n = 0 . (15)
where u( s) is the fluid velocity at the boundary represented by s and d s/dt the velocity of the boundary. Equation ( 15) describes thus the shape deformation of the droplet.

Charge dynamics at the interface

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

j s = σ(u τ + λE τ ) τ (16)
has thus two contributions. The first is the velocity field of the fluid tangent to the surface u τ and the second is the tangent component of the electric field E τ times the charge mobility λ. We may safely neglect the diffusion term which is usually 2 orders of magnitude smaller than the two former terms. The surface charge evolution eventually reads [19][START_REF] Saville | The Taylor-Melcher Leaky Dielectric Model[END_REF][START_REF] Castellanos | Electrohydrodynamics[END_REF],

Dσ Dt -2Hσu n + ∇ s • ( τ σ(u τ + λE τ )) = κE n (17)
where the convective time derivative means the derivative when we follow the surface along a direction normal to itself, that is Dσ Dt = ∂σ ∂t + u n ∂σ ∂n . The second left hand term (LHT), 2Hσu n , accounts for the deforming surface.

The third LHT involves the surface divergence operator ∇ s • which is defined in the supplementary material (36) and accounts for the charge transport along the surface.

The last term, proportional to the bulk conductivity κ, plays the role of a source term, which moves charges from the bulk to the surface or moves charges from the surface into the bulk. Eventually, the divergence-free electric field inside, E, and outside , E (g) , the droplet are deduced by requiring that τ • ( E (g) -E) = 0 and n • ( E (g) -ε r E) = σ/ε 0 at the interface. Note that ( 17) differs from the one used in [12], where the drag of charges by the fluid velocity u τ was omitted and where a constant surface conductivity was used instead of λσ.

A. Viscous Corrected Viscous Potential Flow approximation (VCVPF)

The VPF approach

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

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

∂Ψ ∂t + u 2 2 + p ρ = 0 , (19)
where the liquid pressure at the interface reads (11) ,

p = P γ -P c + 2µD nn , (20)
with the difference that here the strain rate tensor D is computed using only the irrotational velocity field ∇Ψ. Equations ( 19) and ( 20) yield the viscous potential flow model (VPF). Compared to the inviscid liquid model used in [4], the pressure has now an additional term which is proportional to the viscosity µ. Note that, while for an irrotational flow the viscous force per unit volume is zero inside the drop, ∇ • D = 0, its contribution to the surface pressure is generally non-zero

D nn = 0 and D τ n = 0 . (21)
As a result, the jump of the shear stress T τ n across the interface violates the zero shear stress boundary condition at a free surface,

T τ n -M τ n = µ2D τ n = σE τ . (22)
Indeed, VPF does not have enough degrees of freedom to enforce that the shear strain rate D τ n equals the shear stress of the Maxwell tensor. The rotational part of the velocity field is missing here to fulfill the condition [START_REF] Padrino | [END_REF].

In the particular case of a perfectly conducting liquid, where the shear stress of the Maxwell tensor is zero, the shear strain rate D τ n should vanish at the interface, which is in contradiction with [START_REF] Castellanos | Electrohydrodynamics[END_REF]. Compared to a rotational flow where the condition ( 22) can be enforced, VPF results in an incorrect dispersion relationship for small surface oscillations, i.e. the damping of surface modes is underestimated for low viscosities [START_REF] Padrino | [END_REF].

The VCVPF model

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

p = P γ -P c + 2µD nn -p µ . (23)
The viscous pressure p µ satisfies the Laplace equation in the drop volume Ω and is defined such as the power of the viscous pressure balances the power of the non-zero shear stress across the boundary for each surface element ds 2 ,

u n (-p µ )ds 2 + u τ (-2µD τ n + σE τ )ds 2 = 0 . ( 24
)
The viscous pressure being assumed harmonic, Eq. ( 24) defines uniquely the viscous potential p µ in the droplet volume. Combining ( 19), ( 23) and ( 24) yields the time evolution of the velocity potential in the viscous correction for VPF (VCVPF) model. The viscous pressure p µ adds dissipation at the interface, that originates from the shear stress boundary condition. However, energy dissipation within the bulk is still absent. In the case of finite conductivity, the tangential component of the electric field may be non-zero, and p µ accounts also for the traction originating from the non-vanishing shear stress of the Maxwell tensor σE τ at the interface. The technicalities for the evaluation of the viscous pressure p µ are detailed in the Appendix Sec. 2. We want to emphasize that the VCVPF model can also be derived from the dissipation model, where no ad hoc viscous pressure is introduced, but where the zero shear stress boundary condition is implicitly satisfied [START_REF] Padrino | [END_REF].

Numerical approach

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

Validity of VCVPF

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

B. Dimensionless quantities

We introduce the characteristic time τ 0 , giving the timescale of the deformation, and the characteristic capillary pressure P • γ of the spherical drop,

τ 0 = ρR 3 γ , P • γ = 2γ/R . ( 25
)
Defining the dimensionless time t = t/τ 0 , velocity potential Ψ = Ψ τ 0 R 2 , capillary pressure -H = P γ /P • γ and viscous pressure pµ = p µ /P • γ , allows putting the time evolution of the velocity potential Ψ in VCVPF into a dimensionless from

∂ Ψ ∂ t + ũ2 2 + 2 -H -X Pc + Oh Dnn -pµ = 0 , ( 26
)
where the quantities with a tilde are dimensionless. We introduced in (26) the dimensionless viscosity Oh (Ohnesorge number),

Oh = µ √ ργR (27)
which relates the viscous forces to inertial and surface tension forces. Larger Oh indicate a greater influence of the viscosity. Note that pµ depends on Oh via (24). We further introduced the fissility parameter

X = P • c P • γ = Q 2 64π 2 ε 0 γR 3 , (28)
defined as the ratio of the Coulomb pressure over the surface pressure for the spherical droplet shape. If X < 1 the spherical shape of droplet is stable. If X = 1, the spherical droplet is unstable and deforms spontaneously. The dimensionless Coulomb pressure Pc found in ( 26) is defined by

P c = Pc P • c = X Pc P • γ , yielding Pc = σ + ε r Ẽn 2 -ε r Ẽ2 n + (ε r -1) Ẽτ , ( 29
)
where σ and Ẽ are given by the relations σ = σ 4Xε 0 γ/R and E = Ẽ 4Xγ/(ε 0 R).

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

τ -1 s = λ Xγ R 3 ε 0 , τ -1 b = κ ε 0 ε r ( 30
)
allows to define the dimensionless charge carrier mobility λ and bulk conductivity κ

λ = τ 0 τ s = λ 4Xρ ε 0 (31) κ = τ 0 τ b = κ ρR 3 γε 2 0 ε 2 r (32)
Note that λ is proportional to √ X and thus to the charge |Q| carried by the droplet. Equation ( 17) can now be put into the dimensionless form

Dσ D t -2 H σũ n + ∇ s • ( τ σ(ũ τ + λ Ẽτ )) = κε r Ẽn . ( 33
)
Note that the factor ε r , found in the last term, simply originates from the definition of the bulk relaxation time τ b [START_REF] Weisstein | Superellipsoid." From MathWorld -A Wolfram Web Resource[END_REF]. The charge dynamics at the interface is controlled by two parameters namely λ and κ. If λ 1 or κ 1, then the charge relaxation dominates the rate of motion of the fluid and the surface charge distribution is close to equilibrium, meaning that the electric potential at the surface is an equipotential. We recover thus the case of a perfectly conducting liquid, where the surface charge distribution is always at equilibrium during the droplet deformation, σ = σeq , so that Pc = σ2 eq depends merely on the droplet shape s. From equation (26) we deduce that, in the particular case of a perfectly conducting liquid, the shape deformation depends only on the dimensionless viscosity Oh and on the fissility X. If the droplet is critically charged, X = 1, the shape deformation depends only on the dimensionless viscosity Oh.

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

III. SHAPE DEFORMATION OF CRITICALLY CHARGED DROPLETS

We took snapshots of the ultrafast shape deformation of pure water droplets at the Rayleigh limit. The experimental setup used for trapping and taking snapshots of the droplet is similar to the one presented in [3][4][5]. A charged droplet is injected into a Paul trap, where it evaporates until it reaches the critical fissiliy 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 of the time evolution of the shape deformation. At the onset of the instability, the trapped water droplets have a = 52µm. Assuming that water at room temperature has a viscosity of 1 cP, a density of 1000 kg/m 3 and a surface tension of 0.072 N/m, the droplets are characterized by a low dimensionless viscosity of Oh=0.023 (see table I), and are thus particularly good candidates for the simulation with our VCVPF code. Surprisingly no jets were detected on the snapshots, neither for positively nor for negatively charged water droplets at room temperature. We took more than 30 snapshots of the moment where the elongated droplet shapes have pointed ends, but found no trace of small daughter droplets or jets on the snapshots, unlike for glycol and supercooled water [3,4]. We concluded that the jet was too thin to be observed and estimated that the jet radius, if any, must be smaller than the resolution of the snapshots of 0.36 µm. These observations are corroborated by the experimental results presented in the Ph.D work of R. Müller [START_REF] René | Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen[END_REF]. 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 the deionized water droplets were positively charged. The measurements in [START_REF] René | Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen[END_REF] 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 κ b = 5.5 × 10 -6 S/m and a charge mobility of λ = 37 × 10 -8 m 2 /s/V (see table II). We compare in Fig. 1 the simulated shapes to the observed ones. The red simulated shapes in Fig. 1 are critically charged droplets (X = 1) and show the ultra-fast shape deformation after the onset of the instability. Once the pointed tips are formed, we removed, as suggested by [START_REF] René | Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen[END_REF], 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 subcritical shape deformation (green shapes). Especially, the peanut and diamond shape have been found by the simulation. The general agreement validates our approach for describing low-viscous liquid droplets.

Superellipsoid:

In order to compare easily the deformation pathways of different droplets, the shapes are fitted with the mathematical expression of a superellipsoid of the form [START_REF] Weisstein | Superellipsoid." From MathWorld -A Wolfram Web Resource[END_REF]

z a n + x 2 + y 2 b n = 1 . (34)
The shapes are fitted with the constraint of constant volume so that the fit had only two independent parameters, namely a/b and n. Expression (34) can fit either capsules, spheroids or pointed shapes, depending on n.

The aspect ratio a/b gives the elongation of the droplet. Adjusting (a/b, n) for the droplet shape for each snapshot yields a number of points that can be plotted on a graph with a/b as the abscissa and n as the ordinate, so that the complete droplet deformation path can be represented reasonably well on a 2D graph. Fitting the shapes by a superellipsoid allows us to compare the deformation pathway of droplets having different liquid properties. Such a graphical representation of the deformation pathway has already been used in [4] to highlight the differences in the shape deformation of glycol and supercooled water droplets or in [10,11] to show the influence of the viscosity on the deformation pathway. Below, we give examples of the fit parameters (a/b, n) found for various observed shapes. We took ∼ 270 snapshots of the fast deforming water droplet and fitted the shapes on the snapshots by the superellipsoid function. The resulting set of (2.9 , 1.53) (2.1 , 2.0) (0.65, 1.9) (1.1, 1.5) (2.7 , 2.9) FIG. 2: Fit parameters (a/b, n) found for various observed shapes using the superellipsoid equation (34). Red contours show the corresponding superellipsoid shape.

(a/b, n) was reported on a graph, yielding the pathway of the droplet deformation, see Fig.

A. Influence of the viscosity

In a first step, we focus on the influence of the dimensionless viscosity Oh on the deformation pathway. We assume here that the fluid is perfectly conducting, so that Oh is the only free parameter controlling the deformation pathway. We simulated for different Ohnesorge numbers Oh, ranging from 0.001 to 0.8, the shape deformation of critically charged droplets (X = 1). The pathways for several values of Oh are shown in Fig. 4. In all simulations, the initial shape was a spheroid n = 2 with a ratio a/b = 1.03 (almost spherical). First, we note that the general deformation pathway is similar for all cases: the instability starts with a/b = 1.03 and n = 2. Then a/b increases, while n remains nearly constant at 2, meaning that the initial deformation remains mainly prolate spheroidal. With ongoing elongation of the drop, n starts to decrease and the shape becomes lemon-like with pointed tips. When looking in more detail, we find that with in-creasing Oh and thus with increasing viscosity of the liquid, the droplet adopts a more and more elongated shape before pointed tips appear at the ends, typically when n ∼ 1.5. Indeed, the viscous forces at the interface are stronger for fast moving/deforming boundaries, or alternatively, higher surface modes are more damped than lower modes. As a result, for larger Oh, the tip formation is retarded during the droplet elongation and consequently, tips appear for shapes with higher aspect ratios. For Oh=10 -3 , the maximal elongation is characterized by the ratio a/b = 2.66, close to the value of 2.65 found for inviscid fluids [4]. At the other extremity, for Oh=0.4, the aspect ratio a/b is 3.72, close to the factor 3.86 found by [7,8] in the case of Stokes flow (creeping flow with vanishing inertia forces) corresponding to large Oh 1. The trend is shown in Fig. 5, where we give the aspect ratio a/b of the droplet at the moment n = 1.5 as a function of Oh ∈ [0, 1.4]. For comparison, we also add the results obtained for the VPF model. In the low viscosity domain Oh< 0.01, when the dimensionless viscous pressures Oh Dnn and pµ are negligible compared to the difference between the capillary and the Coulomb pressure, VPF and VCPPF give similar pathways and final aspect ratios. However because VCVPF recovers in that domain the correct dispersion relationship, as discussed in a previous section, VCVPF is preferable over VPF also for low viscous liquids. For larger dimensionless viscosities Oh> 0.8, the aspect ratio in VCVPF tends to the one obtained in the Stokes flow approximation, while VPF clearly overestimates the elongation of high viscous drops confirming that VCVPF may be better suited for a larger range of viscosities. For water droplets, characterized by Oh=0.023, the differences in the pathways between VPF and VCVPF are indeed small but nevertheless notable. Typically, a VCVPF calculation with Oh=0.023 is similar to a VPF calculation with Oh=0.046.

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

The shape deformation of viscous drops shows a similar behaviour than the damped forced oscillator. For Oh< 0.2, the dynamics is under-damped and the return to the sphere goes through oblate shapes as indicated by the ratio a/b < 1. In the (a/b, n) graph, the droplet oscillates with decreasing amplitude, between the prolate and oblate shapes and between pointed (diamonds) and flattened (capsule) tip shapes. We note that Oh 0.2 corresponds to the value were the droplet is critically damped, meaning that it no longer overshoots when returning to the asymptotic spherical shape. For Oh> 0.2, the motion is overdamped and the droplet returns slowly to the stable spherical shape. In the particular case of Oh≥ 0.4, the return path is even characterized by n ≤ 2, meaning that the tips stay pointed for a while even after charge emission. This is a surprising result and the prediction may not be meaningful as VCVPF calculations with Oh≥ 0.4 are expected to deviated significantly from full NS calculations, as was shown in [START_REF] Padrino | [END_REF] in their study of small drop surface oscillations. FIG. 6: Deformation pathway of a droplet with Oh=0.023 (typically a = 50µm water droplet at room temperature) for various dimensionless charge carrier mobility λ, ranging from 0.1 to 10 and for the perfect conductor (p.c.).

B. Influence of the charge carrier mobility

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

The reason why smaller charge mobilities λ tend to accelerate the tip formation is explained by the competition between the charges being dragged by the fluid and the charges being pushed by the tangent component of the electric field at the surface. For low-viscous drops, Oh< 0.03, the velocity at the surface may be so high that the flow drags a large amount of charges to the tips. The large concentration of charges at the tips tends to accelerate the formation of pointed ends. If the accumulated charge density at the tips exceeds the equilibrium charge density σ 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 τ -1 s and the closer the charge distribution stays to the equilibrium charge distribution σ eq . That is why the deformation pathway tends to the p.c. case for λ 1. On contrary, for λ 1 the surface charges move mainly with the fluid to the ends, where they accumulate. The finite charge mobility has thus a non-negligible effect on the shape deformation and should be taken into account when looking at the pathway in detail. We will see below that the charge mobility allows to discriminate between positively and negatively charged pure water droplets.

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

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

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

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

It seems thus that the reason to why positively charged pure water droplets have a aspect ratio a/b larger than negatively charged ones is because H + has a higher charge mobility than OH -. For such low-viscous droplets, Oh= 0.023, the drag of the charges by the fluid velocity is in competition with the push by the electric field. It happens that for OH -the latter effect is less efficient so that more charges accumulate at the tips, accelerating the tip formation. The charge velocity tangent to the surface, given as the sum of the drag velocity ũτ and "electric field" velocity λ Ẽτ is illustrated in the upper panel of Fig. 8 for three different values of λ. The quantities in Fig. 8 are given as a function of the dimensionless distance r/R, corresponding to the distance from the symmetry axis to a surface point, such that r → 0 corresponds to the tips (see inset). We note that close to the tips (r/R 0.06), the velocity of the charges is reduced by a factor 2 when increasing λ = 1 to λ = 4. The accumulated charge density at the tips (lower panel) for λ = 1 was also found to be larger about by a factor 2 than for λ = 4, accelerating thus the tip formation. The higher charge density may also explain why negatively charged pure water droplets emit more charge than positively charged pure water droplets, as was observed in [START_REF] René | Stabilität und Zerfallsdynamik hochgeladener Flüssigkeitströpfchen[END_REF]. If instead of pure water we had used sea water, having a bulk conductivity of 5 S/m so that κ 10 4 λα for all excess charge carrier α, the charge dynamics would be dominated by the bulk current κ Ẽn and the liquid well approached by a perfect conductor. We expect thus positively and negatively charge droplets of sea water to yield the same deformation path and emit the same amount of charge at charge break-up. The latter prediction may be verified in future experiments.

IV. CONCLUSION

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

SUPPLEMENTARY MATERIAL 1. Surface divergence

Let τ be a vector tangent to the surface defined by the vector s. The surface divergence operator

∇ s • (Ψ( s) τ ( s)) (35)
for arbitrary function Ψ( s) defined on the surface, may not be well known to the reader, so we give a computationally useful expression of the surface divergence found in (17). Let (ξ > 1, -1 ≤ η ≤ 1, 0 ≤ ϕ ≤ 2π) be curvilinear orthogonal coordinates and r(ξ, η, ϕ) an arbitrary vector. Because the droplet deformation and charge dynamics are axisymmetric, the function Ψ(ξ, η) = σ(u τ + λE τ ) does not depend on the the azimuthal angle ϕ. The surface vector s(η) depends thus only on η. The general expression simplifies in that case to where τ = e ϕ ∧ n is the tangent vector and h η and h ϕ the scale factors associated to the curvilinear coordinate system and evaluated at the surface. The term d()/dη in (36) can be accurately computed using for example a Gauss-Lobatto interpolation scheme, which is appropriate if the endpoints are know (here h ϕ ( s)Ψ( s) = 0 for η = ±1).

Numerical method for computing the viscous potential

The VCVPF model is only useful if one can compute the viscous pressure p µ , which enters the normal stress term. We propose here a numerical method to compute the viscous pressure in VCVPF. Making use of the fact that the Laplace operator is separable in prolate spheroid coordinate system (ξ ≥ 1, -1 ≤ η ≤ 1, 0 ≤ ϕ ≤ 2π), we get a simple expression for the harmonic functions f m that satisfy the Laplace equation ∇ 2 f m = 0. As we want f m (ξ, η) = f m (ξ, -η) to be axisymmetric and even and defined at the origin f (1, 0) < ∞, we have f m (ξ, η) = P 2m (ξ)P 2m (η) (37

)
where P m () are Legendre polynomials of order m. Let {f m } be a basis of harmonic function. We may now expand the viscous pressure on the harmonic basis f m ,

p µ (ξ, η) = M m=1 b m f m (ξ, η) (38)
The expansion coefficients {b m } are deduced using the boundary condition (24),

M m=1 b m f m ( s i )u n ( s i )ds i Fm,i = u τ ( s i )(2µD τ n ( s i )ds i di -u τ ( s i )σ i ( s i )E τ ( s i ))ds i ei ( 39
)
where s i≤N is an ensemble of N points at the surface. Equation (39) is a system of N equations and M unknowns of the form

M m=1 b m F m,i = d i -e i , 1 ≤ i ≤ N (40)
If N ≥ M the system is overdetermined and can be solved using a singular value decomposition (SVD), which finds the coefficient {b m } as least squares solutions.

FIG. 1 :

 1 FIG. 1: 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 breakup, 20 % of the charge is removed and the fissility for the green droplet is X = 0.64. Shape labels: (a) lemon, (b) peanut, (c) diamond.

3 .FIG. 3 :

 33 FIG. 3: Graph (a/b, n) showing the deformation pathway of a critically charged deionized water droplet. Starting with the spherical shape (1, 2), arrows give the time evolution of the path.

4 FIG. 4 :

 44 FIG.4: Deformation pathway of a critically charged droplet as a function of the Ohnesorge number Oh. At the pointed tip appearance (n = 1.5), the fissility is reduced from X = 1 to X = 0.64.

FIG. 5 :

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

FIG. 7 :

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

FIG. 8 :

 8 FIG. 8: (a) Dimensionless velocity of the charges tangent to the surface, ũτ + λ Ẽτ , as a function of the dimensionless distance r/R defined in the inset. Positive values indicate velocities directed towards the tips. (b) Dimensionless surface charge density σ as a function of r/R. The quantities are calculated for three different charge mobilities λ. The blue line with pluses gives the equilibrium surface charge density σeq in the case of a perfectly conducting liquid (p.c.).

∇

  s • ( τ ( s) Ψ( s)) = n • ∇ ∧ (Ψ( r) e ϕ ) r= s = n • e ξ 1 h η h ϕ d dη (h ϕ ( s)Ψ( s)) (36)

TABLE I :

 I Parameters of deionized (pure) water droplets at room temperature.

	26	1	1000	72	0.023	16

R (µm) µ (cP) ρ (kg/m 3 ) γ (mN/m) Oh τ0 (µs)

Acknowledgments

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