J Lemetayer

C Marion

D Fabre

F Plouraboué

Multi-inception patterns of emitter array/collector systems in DC corona discharge

Keywords: ionic wind, corona discharge, ElectroHydroDynamic (EHD), ElectroAero-Dynamic (EAD), electric propulsion, multi-inception, emitter/collector systems, SDMC, finite element, drift-region model, Kaptzov hypothesis

Multiple emitters systems have been previously used so as to increase charge density in the drift region, many times without producing sensible increment neither in total current nor ionic wind. This contribution focuses on analyzing the detailed physics behind this failure, that is named "multiple emitters unscalability". It is established that multiple emitters un-scalability is related to the inability of multiple corona discharge inceptions when increasing the emitter number and/or density. This confirms recent findings that corona discharge inception is shielded by electro-static interactions between emitters. This contribution demonstrates that this shielding can be balanced by emitter/collector electrostatic interactions depending on the considered configuration.

For sufficiently close collector-emitter distances, ignition starts at the array center, whereas, on the contrary, when the collector is distant, the ignition not only starts at the array's periphery but might also be limited there.It is also demonstrated that emitter/emitter electrostatic interactions can be balanced by emitter/collector ones, depending of their chosen configuration. This lead to a variety of multi-inception patterns, the condition of which are analyzed. Intermediate configurations for which the collector is neither sufficiently close nor distant from the emitter array center provide a variety of multi-inception patterns that are hereby analyzed. Combining finite element computations of multiinception drift-diffusion modeling with experimental measurements, provides a coherent picture explaining why multiple emitters sources systems do not lead to full ignition, and also exhibit conditions for which it does, leading to multiple emitters scalable systems.

Introduction

Ionic wind is of interest in distinct applications such as electrostatic precipitors [START_REF] Davis | Wire-duct precipitator field and charge computation using finite element and characteristics methods[END_REF][START_REF] Lu | Numerical simulation on the fine particle charging and transport behaviors in a wire-plate electrostatic precipitator[END_REF], Aero-Electro-Dynamic (AED) [START_REF] Wilson | An investigation of ionic wind propulsion[END_REF][START_REF] Masuyama | On the Performance of electrohydrodynamic propulsion[END_REF][START_REF] Konstantinos | Electro-Hydrodynamic (EHD) Thrust Analysis in Wire-Cylinder Electrode Arrangement[END_REF][START_REF] Moreau | Electrohydrodynamic force produced by a wire-to-cylinder dc corona discharge in air at atmospheric pressure[END_REF][START_REF] Gilmore | Electrohydrodynamic thrust density using positive corona-induced ionic winds for in-atmosphere propulsion[END_REF][START_REF] Shibata | Performance prediction of electrohydrodynamic thrusters by the perturbation method[END_REF][START_REF] Christopher | Electroaerodynamic thruster performance as a function of altitude and flight speed[END_REF][START_REF] Monrolin | Electrohydrodynamic Thrust for In-Atmosphere Propulsion[END_REF][START_REF] Vladislav | In-atmosphere electrohydrodynamic propulsion aircraft with wireless supply onboard[END_REF][START_REF] Plouraboue | Flying with ionic wind[END_REF][START_REF] Xu | Flight of an aeroplane with solid-state propulsion[END_REF][START_REF] Xu | A dielectric barrier discharge ion source increases thrust and efficiency of electroaerodynamic propulsion[END_REF][START_REF] Krishna | A laser-microfabricated electrohydrodynamic thruster for centimeter-scale aerial robots[END_REF][START_REF] Xu | Higher thrustto-power with large electrode gap spacing electroaerodynamic devices for aircraft propulsion[END_REF][START_REF] Sankar Vaddi | Analytical model for electrohydrodynamic thrust[END_REF][START_REF] Leng | An empirical formula for electrohydrodynamic thrust in the point-to-grid configuration with small electrode gaps[END_REF][START_REF] Gomez-Vega | Performance of decoupled electroaerodynamic thrusters[END_REF], gas and ionic pumps [START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF][START_REF] Chang | Narrow-flow-channel-driven EHD gas pump for an advanced thermal management of microelectronics[END_REF][START_REF] Zhang | Numerical and experimental investigation on configuration optimization of the large-size ionic wind pump[END_REF][START_REF] Zhang | Experimental and Numerical Investigation on Flow Characteristics of Large Cross-Sectional Ionic Wind Pump With Multiple Needlesto-Mesh Electrode[END_REF][START_REF] Ye | The Effects of Inlet Blockage and Electrical Driving Mode on the Performance of a Needle-Ring Ionic Wind Pump[END_REF][START_REF] Zeng | Analysis of a two-stage ionic wind pump with multiple needle-to-mesh electrodes for cooling electronics[END_REF][START_REF] Kawamoto | Development of Ionic Wind Pump Used in Martian Environment[END_REF], particle analyzer [START_REF] Mazumder | Development of a dust particle analyzer for in-situ and simultaneous measurements of size and charge distributions of martian dust[END_REF][START_REF] Sankar Vaddi | Behavior of ultrafine particles in electrohydrodynamic flow induced by corona discharge[END_REF], miniaturized heat cooler [START_REF] Johnson | Impingement cooling using the ionic wind generated by a low-voltage piezoelectric transformer[END_REF][START_REF] Wang | Review on the recent development of corona wind and its application in heat transfer enhancement[END_REF][START_REF] Ho Shin | Analysis of thermal flow around heat sink with ionic wind for high-power LED[END_REF][START_REF] Qu | Experimental Investigation on Flow and Heat Transfer Characteristics of a Needle-Cylinder Type Ionic Wind Generator for LED Cooling[END_REF][START_REF] Kalman | Enhancement of heat transfer by means of a corona wind created by a wire electrode and confined wings assembly[END_REF] and xerography, i.e electro-photography. In most of these applications, the increase of ionic wind is interesting so as to achieve larger effects in either convection, heat or mass transfer, and/or propulsion.

Ionic wind is mostly generated from the drift of created charges produced in corona discharges localized nearby thin high-electric field regions which, most often, are high-tension small wires, fittingly called emitters. Since the created charges are drifting inside the domain and collected into a reference potential electrode, this collecting electrode is in return called the collector.Increasing the charge creation directly leads to a corresponding level-off of the resulting ionic wind. This is why, many attempts have been pursued to consider multiple emitter sources in order to produce multiple sites for charge emissions [START_REF] Wilson | An investigation of ionic wind propulsion[END_REF][START_REF] Masuyama | On the Performance of electrohydrodynamic propulsion[END_REF][START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF][START_REF] Fylladitakis | Design of a Prototype EHD Air Pump for Electronic Chip Cooling Applications[END_REF][START_REF] Lee | Parallel integration of ionic wind generators on PCBs for enhancing flow rate[END_REF][START_REF] Belan | Compared ionic wind measurements on corona and DBD plasma actuators[END_REF]. This is why charge sources densities increase has been pursued in many applications.

As stated in [START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF], "the challenge in the field is no longer one of understanding the physics but of engineering the device (electrodes, geometry, operation) to optimize the flow". This is why an active area of research in EHD pumps has been devoted to the analysis of multiple electrodes (multiple emitters/multiple collector) in needle to grid or wire to grid arrays [START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF][START_REF] Fylladitakis | Design of a Prototype EHD Air Pump for Electronic Chip Cooling Applications[END_REF][START_REF] Lee | Parallel integration of ionic wind generators on PCBs for enhancing flow rate[END_REF][START_REF] Belan | Compared ionic wind measurements on corona and DBD plasma actuators[END_REF]. When run within parallel emitters systems the resulting ionic wind could, in some cases, produce a smaller flow speeds relative to the discharges operating individually [START_REF] Wilson | An investigation of ionic wind propulsion[END_REF][START_REF] Masuyama | On the Performance of electrohydrodynamic propulsion[END_REF][START_REF] Belan | Compared ionic wind measurements on corona and DBD plasma actuators[END_REF][START_REF]Handbook of Mathematical Functions[END_REF][START_REF] Lee | Parallel integration of ionic wind generators on PCBs for enhancing flow rate[END_REF]. Better results have been obtained in some configurations from introducing a shielding layer in-between emitter and collector electrodes [START_REF] Lee | Parallel integration of ionic wind generators on PCBs for enhancing flow rate[END_REF][START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF]. Shielding effects arising from electrostatic interactions between emitters and collectors have thus been identified as one important effect arising in those applications [START_REF] Arif | The influence of design parameters on the occurrence of shielding in multi-electrode ESPs and its effect on performance[END_REF][START_REF] Sh | A semi-analytical electrostatic dipole approach for systems of parallel metallic cylinders[END_REF][START_REF] Guerra-Garcia | Corona Discharge in Wind for Electrically Isolated Electrodes[END_REF]. Moreover, the possible effect of hydrodynamic interactions between the ionic-wind jets created in these multiple emitters systems has also been suspected as possible explanation, for modified behaviours [START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF]. Islamov [START_REF] Sh | A simple dipole approach for the calculation of electrostatic field and onset corona voltage in wires-to-cylinders electrode configurations[END_REF] recently analyzed the electrostatic shielding effect between an array of emitters. Using a nice semi-analytic method based on combining bi-cylindrical elementary solutions developed in [START_REF] Sh | A semi-analytical electrostatic dipole approach for systems of parallel metallic cylinders[END_REF], he found that the inception criteria significantly differ inside the emitter array as opposed to its far end, leading to preferential inception at the emitter array periphery. A similar preferential peripheral ignition has also been reported in wire-plate electrostatic precipitators configuration [START_REF] Arif | The influence of design parameters on the occurrence of shielding in multi-electrode ESPs and its effect on performance[END_REF]. In this contribution it is shown that this scenario not always happens. Depending on the considered configuration of emitters, various inception patterns can be found, either showing preferentially central or outskirts ignitions. Furthermore, considering secondary and further emitter's ignitions many possible ignition patterns can be found depending on the chosen emitter/collector configuration. These secondary ignitions are not satisfactorily captured by a simple electro-static analysis because they are influenced by charges built in already ignited emitters. Hence, a non-linear driftdiffusion modeling which takes into account the influence of created charges in driftregion is needed for a better prediction of multi-ignited configurations.

The paper is organized as follows. Section 2 details theoretical and experimental methods. Section 2.1 explains how the non-linear drift-diffusion problem, is solved for multi-inception systems. This approach follows [START_REF] Kaptsov | Elektricheskie Yavleniya v Gazakh i Vakuume[END_REF][START_REF] Tirumala | Comparative study of corona discharge simulation techniques for electrode configurations inducing non-uniform electric fields[END_REF][START_REF] Chen | Numerical investigation of an electroaerodynamic driven aeroplane: Electrical properties, ionic wind and flight performance[END_REF][START_REF] Xu | Higher thrustto-power with large electrode gap spacing electroaerodynamic devices for aircraft propulsion[END_REF][START_REF] Kaptsov | Elektricheskie Yavleniya v Gazakh i Vakuume[END_REF][START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF] for which a Kaptzov approximation is used for corona discharge modeling. Nevertheless, a dynamically selected boundary condition is adopted here, where either inception or non-inception is possible for every emitter. The choice setting the correct emitters boundary conditions is provided in Appendix A.2. Section 2.2 details the experimental methods used, with set-up and signal post-processing detailed in 2.2.1 and 2.2.2. The results for multi-inception of multiple emitters/collectors systems are presented in Section 3. Section 3.1 details the multi-inception patterns diagrams of various configuration for which experimental measurements are confronted with numerical predictions. Finally section 3.2 analyses the current-voltage relations as well as the overall ionic wind production from measuring the resulting propulsive thrust.

Methods

Drift-diffusion modeling of multiple ignited emitter : multi-Kaptzov model

As commonly done (see, e.g. [START_REF] Johnson | Recent advances in electrohydrodynamic pumps operated by ionic winds: A review[END_REF][START_REF] Khaddour | Numerical solution of the corona discharge problem based on mesh redefinition and test for a charge injection law[END_REF][START_REF] Adamiak | Corona discharge in the hyperbolic point-plane configuration: direct ionization criterion versus approximate formulations[END_REF]), the ionic wind modeling used here relies on a drift-diffusion equation governing the overall charge density ρ (corresponding for positive DC corona to positive charges only created in the corona discharge region). As in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], a non-dimensional potential φ = ϕ/V and a non-dimensional charge density ρ = ρ

Where ∂Ω e I stands for ignited emitters boundaries, ∂Ω e O for non-ignited ones, Ω c for collector boundaries, and ∂Ω E for lateral external boundaries (Cf Fig 1 and (A.13) for more precise definitions). Furthermore, Êa = E P eek R c /V is the non-dimensional electric field at the ignited emitters as predicted by the Peek law (A.11). Note that only one boundary condition is applied at the collector, as a consequence of the hyperbolic dominated nature of Eq. ( 2).

For pre-ignited states (i.e. for voltages V below the threshold of initial ignition), the problem reduces to the simple electrostatic problem (with ρ = 0) which is fairly easy to solve. A resolution method using finite-elements is presented in Appendix A.1. Since the problem is linear, the electric field at each emitter is directly proportional to the applied voltage, so the Kaptzov condition directly indicates the first ignition threshold as well as the emitter(s) where this initial ignition will occur.

Above the ignition threshold, one has to solve the full nonlinear problem (3). In [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], a resolution method based on finite-element formulation and Newton iteration was introduced to solve this set of equations considering a single emitter or a pair of emitters placed in a symmetric way and considered as both ignited. The method can be directly used for situations with a larger number of emitters considered as all ignited, as done for instance in [START_REF] Zhao | Numerical analysis of forces in an electrostatic levitation unit[END_REF][START_REF] Martins | Modeling of an EHD corona flow in nitrogen gas using an asymmetric capacitor for propulsion[END_REF][START_REF] Colas | Ionic wind generation by a wire-cylinder-plate corona discharge in air at atmospheric pressure[END_REF][START_REF] Martins | On the influence that the ground electrode diameter has in the propulsion efficiency of an asymmetric capacitor in nitrogen gas[END_REF][START_REF] Chen | Numerical investigation of an electroaerodynamic driven aeroplane: Electrical properties, ionic wind and flight performance[END_REF][START_REF] Adamiak | Quasi-stationary modeling of the DBD plasma flow control around airfoil Quasi-stationary modeling of the DBD plasma flow control around airfoil[END_REF][START_REF] Xu | Higher thrustto-power with large electrode gap spacing electroaerodynamic devices for aircraft propulsion[END_REF].

Nevertheless, in real systems, all emitters may not ignite simultaneously, due to the electrostatic shielding between emitters. Consequently, different boundary conditions have to be applied on ignited and non-ignited emitters. A resolution method, extending the one of [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF] but considering an assumed ignition pattern, is presented in appendix Appendix A.2. Having at hand a resolution method considering an assumed ignition pattern is still not sufficient, since for a given geometrical configuration and value of the applied voltage, the ignition pattern is not known a priori. Therefore, an algorithm to identify the correct ignition pattern is needed. Such an algorithm is presented in detail Appendix A.3. In short, an assumed ignition pattern is first selected a priori. The problem is solved using the numerical method of appendix Appendix A.2, and the consistency of the assumed ignition pattern is checked a posteriori by physical arguments. Namely, if an emitter was assumed as "non-ignited" and the computed electric field at this emitter exceeds the value corresponding to the Peek law (A.11), the hypothesis has to be revised as "ignited" ; on the other hand, if an emitter was assumed "ignited" and the computed charge density predicts negative values in its vicinity the hypothesis has to be revised as "non-ignited". The problem is then re-solved with the revised hypothesis and the process is possibly iterated up to convergence.

Finally, as in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], we numerically checked the effect of Neumann boundary conditions applied at the external border ∂Ω E of a finite domain-size. We consider domain size being 15 to 20 times larger than the emitter-collector distance, and found a very weak influence of the domain size (less than 0.5%) on the quantity of interest (total current, thrust). Beside computational methods and validations, quantity of interest such as total current intensity and thrust are evaluated from the numerical computations as in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF]. Concerning the current, it is evaluated from one dimensional integration of the local current intensity -µρ∇ φ • n along the emitter frontier having normal n as well as the collector one. These two distinct estimations are found closed within 0.1%. The thrust is evaluated from the bulk integration (in 2D a twodimensional integration) of the local Coulomb force -ρ∇ φ. An alternative estimate for the thrust can also be obtained from a surface integration (in 2D, a one dimensional integration) of the Maxwell stress projection over the emitter and collector borders as detailed in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF]. These 2D or 1D integration formulations for the thrust provides a very similar estimate up to 0.5%. The 'incepted' emitters, here nearby the center, are surrounded by a red halo, to exemplified an inception pattern. Definitions of the various geometrical parameters are given.

Experimental method

Set-up

A PolyTetraFluoroEthylene (PTFE) frame, previously used in [START_REF] Monrolin | Electrohydrodynamic Thrust for In-Atmosphere Propulsion[END_REF], supported the electrodes of length L e = 39 cm. Emitters were copper wires of radius R e = 50 µm, each of them maintained straight in tension by a spring, and the collector was a steel cylinder of radius R c = 5 mm. Slots into the supports of the emitters enabled to adjust the space d between the emitters by step of 5 mm, while the distance L (face-to-face) between the emitters and the collector was chosen at 45 mm or 100 mm. The frame was hung around 50 cm below the digital scale (Mettler Toledo ME3002) used to measure the thrust generated by the system. The voltage applied to the emitters was monitored on the High Voltage (HV) supply (Iseg HPp 700505, voltage range 0/+70kV maximum intensity 5mA) with an accuracy of ±35 V and the current was measured across a resistance in the line of the collector with an accuracy of ±1.5 µA. For each geometrical configuration, additionnaly to the thrust and current measurements, images of the light emitted by the corona discharge around the emitters were recorded by an UV camera (Princeton PI-MAX 1024 X 1024) mounted with a 200 mm lens, see Figure 2. Its spectral domain lays within 200 and 550 nm. All experiments were carried out in ambiant air at a temperature T = 18 ± 1°C and a relative humidity RH = 45 ± 5%.

Signal post-processing

Labview was used to monitor the devices and to collect all returned data. Data for the thrust, the applied voltage and the current were acquired at 1 Hz, 10 Hz and 100 Hz, respectively, and averaged over 20 s for each working point. Regarding the visualization of the light intensity emitted by the corona discharge, 10 images were summed for each working point. As explained on Figure 3, a background correction based on images with the HV supply off was also applied. Due to the collector being between the mirror and the emitters (Figure 2), some of these ones were hidden. However, the ignition of each emitter of a pair happens at the same voltage and they emit a similar light intensity (6.6 % maximum difference for all configurations except for 5E/1C L=45mm where a maximum 12% difference was found). Thus, for the sake of clarification half of the image has been extracted (red dash in Figure 3) before applying a symmetry based on the central axis and concatenating both parts. A similar method was used to get the complete intensity profile, with an additional step consisting in averaging data along the y-axis in the blue box to improve the signal-to-noise ratio.

Experimental uncertainties

Uncertainties on the thrust, the voltage or the current measurements can be approached by two different ways. The first one considers the systematic error of each material involved in the measurement. Thus, the uncertainty on the thrust depends on the balance accuracy so ∆ thrust = 0.1 mN, and the uncertainty on the voltage is based on the power supply accuracy, ∆ V = 35 V. The uncertainty on the current is estimated from the accuracy of the voltage given by the acquisition box (DAQ, National Instruments), ∆ V I = 2.4 mV where V I is the measured tension used to estimate the current. The second method is to consider the statistical error so the standard deviation of the measurements. Indeed, for each working point (for a given voltage delivered by the power supply), data are acquired during a long enough time to ensure the statistical convergence of the mean and r.m.s values. The systematic error being only an estimation of the measurement error, the statistical error approach is prioritized. However, if the statistical error is smaller than the systematic one, the latter is used. Errors on the thrust and the voltage are directly obtained, but the error on the current is calculated through the method of uncertainty propagation to take into account the uncertainty on the resistors used to measure the tension V I proportional to the current. Error bars on Figures 10 and 11 report the uncertainty on the measurements of thrust and current over voltage ratio. The inception thresholds were determined with a certain accuracy depending on the method used. On the one hand, the uncertainty for the threshold obtained by optical diagnostic corresponds to half of the voltage step used during experiment, i.e. 0.25 or 0.5 kV depending on the configuration. On the other hand, the uncertainty determined by a fit of the current curve is estimated by propagating the systematic errors on the voltage and the current measurements (given above) through the fit equation I/V = a(V -b), with b the inception threshold.

Results

Multi-inception patterns

Figure 4 provides a synoptic view of the inception patterns for all considered configuration. The first column displays the experimentally observed succession of patterns and the corresponding threshold values. The second column displays the same results as predicted using the multi-Kaptzov model of Sec. 2.2. Let us consider, first, the experimental results. The 3E/1C configurations are relatively simple; for proximate configuration for which collector/emitter interaction dominates (i.e N e = 3, L = 45mm, d = 40, 60mm) the first ignition starts in the central position, whereas for distant configurations where emitter/emitter interactions dominates (i.e N e = 3, L = 100mm, d = 40, 60mm) the opposite arises with peripheral ignition. For the 5E/1C configurations, an interesting pattern is observed : ignition happens first for the peripheral emitters, then the central ones and ultimately the intermediate ones. This is qualitatively consistent with the simple view that for large emitter/collector distance, the ignition patterns are determined by a dominant emitter/emitter shielding inhibiting central ignition and favoring peripheral ones. This inhibiting effect also happens for the 9E/1C configurations. For the L = 45mm configuration ignition is first peripheral, and then happens successively on emitter pair E 1 (following the numbering convention of Figure 1(b)), closely after the central one (the corresponding threshold values V c E1 and V c E0 being within a difference of 0.5kV), and finally for emitter pair E 2 . In this case, ignition of emitter pair E 3 is not observed as arc discharge occurs. Finally, for the 9E/1C, L = 100mm configuration, only the peripheral emitters are observed to ignite before arc discharge occurs (around 70 kV)

Consider now the numerical results displayed in the second column of figure 4. It is striking to note that in most cases, the pattern successions are correctly reproduced. The only differences are found for 9E/1C configurations. Here, simulation predicts that ignition happens first on the peripheral emitters, then on the central one, and lately on the intermediate ones starting from the center. For L = 45mm, the observed difference with the experiments lies in inversion of the second and third thresholds (however, as already discussed, in the experiments these two thresholds are very close to each other). The simulation also predicts ignition of all emitters for voltage values above 65kV , a range where arc discharge always happens in experiments (around 45 kV).

Comparing the threshold values it is found that the numerically computed values always over predict the experimentally observed ones. The relative error (last column of the figure) is smallest for the cases where the emitter-to-emitter distance d is largest and grows as d is decreased. Figure 5 characterizes this behavior by plotting the threshold error estimate (up-right column of Figure 4) versus the emitter/emitter distance d. -One observes that a 1/d trend fits well the modeling error. In 2D, such 1/d trend is the signature of electro-static dipolar charges distribution effects which has not been taken into consideration when a local Kaptzov model is used on each emitter. This might be a limitation of the multi-Kaptzov model for which charge distribution around emitters is indeed uniformly imposed, from the imposed uniform Peek's electric field [START_REF] Franck | Dielectric phenomena in high voltage engineering[END_REF] from Gauss theorem. Note that the emitter diameter being 0.1 mm, the experiments still consider a very spaced-out emitter array for which, in the closest configuration, the distance apart is hundred time larger than emitter's diameter.

Hence, Figure 4 demonstrates that a drift-diffusion modeling provides a faithful estimate of ignition pattern for multiple emitter systems, and a reasonable quantitative estimate of successive ignition thresholds. The observation that the threshold prediction precision is not spoiled when increasing the ignition order, suggests that the space-charge effects are well captured by the drift-diffusion/multi-Kaptzov model.

We now characterize the inception patterns in Figures 6,7 & 8 from the direct observation of their glow using a UV camera (Cf. 2.2.1).

Electrons created in the corona discharge gain kinetic energy from the electric field acceleration so as to be able to excite, ionize or dissociate molecules in air such as N 2 , O 2 or H 2 O. Each molecule going under these transitions generates a light emission at a particular wavelength that is particularly intense in UV domain [START_REF] Grum | Spectral emission of corona discharges[END_REF][START_REF] Czech | Optical emission spectroscopy of point-plane corona and back-corona discharges in air[END_REF][START_REF] Mraihi | Electrical and spectroscopic analysis of mono-and multi-tip pulsed corona discharges in air at atmospheric pressure[END_REF]. By measuring the intensity of this light, the density of charge created in the corona discharge can be qualitatively estimated. Figure 6a displays one example of UV camera image of the emitter array for N e = 3, showing a more intense central glowing discharge. The corresponding intensity profile, deduced using the procedure explained in Sec. 2.3.2, is represented in Figure 6b. The upper intensity profile of figure 6b illustrates the peripheral ignition represented in the bottom pattern of the third line of Figure 4. Similarly Figure 7 illustrates the emitter array intensity profile for the three ignition patterns of the fifth line of Figure 4. Interestingly enough, even though the central ignition is less precocious than peripheral one, at higher voltage, its intensity is slightly higher than peripheral peak. This suggests a more intense charge creation in the center for high voltages in this proximal (i.e L = 45mm) configuration.

For configurations with N e = 9 emitters, a rather similar scenario is observed when the first ignition occurs at the peripheral ones, followed by a gradual ignition of the central ones. This can be observed in Figure 8 corresponding to the seventh line of Figure 4. Note that, in this case, the experiment shows that the three central emitters ignite almost simultaneously. Interestingly, the measured light intensities indicate that Figure 9 displays the multi-inception drift-diffusion numerical prediction in various configurations. In these figures, color levels indicate the electric charge density ρ from blue without charge, to red with maximal charge concentration. The configurations are chosen to illustrate various patterns observed experimentally. The values of applied tension V used in the computations were chosen with the same relative difference as in the corresponding experimental ignition threshold of the considered pattern indicated in Figure 4. It is interesting to note that, again, as for experimental ignition profiles, when the emitter density is increased, emitter/collector interaction is not strong enough to maintain a dominating central ignition pattern Figure 8: Same conventions as Figure 7 for N e = 9 emitter system for distance L = 45mm and d = 10mm.

(Cf Figure 9c,d). Moreover, one can still find a clue of this interaction on the developed fields (Cf Figure 9e, f ), i.e when all emitters are ignited, weakening central charge density. This means that emitter/collector interaction competition with emitter/emitter one is important to take into account when considering geometrical settings of emitters to keep charge density as homogeneous as possible. Ignition patterns and profiles provide detailed quantitative view about sources response to the applied voltage but one needs to complement these measurements with more "system oriented" ones for applications.

Current and ionic wind generation in multiple emitters/collectors systems

The overall current-voltage and free stream mass-flow rate are now considered. The free stream mass-flow rate is not directly measured but the thrust, i.e the Electro-Aero-Dynamic (EAD) propulsion force due to momentum transfer into the air associated with ionic wind is considered instead. They are nevertheless closely related. The free stream mass-flow rate, Q, i.e the flow rate of air mass (associated with momentum transfer) given by Q = ρ a S∆V , with air density ρ a , surface of the propulsive system S and ∆V being the difference between inlet and outlet velocity of air resulting from ionic wind. Considering an emitter/collector source system having span along its principal axis, the total thrust is related to the free stream mass-flow rate by relation T = Q 2 /ρ a S. Hence, the thrust is directly proportional to the square of the free stream mass-flow rate, and as such an interesting measurement to perform for applications (this choice is also justified because T is much easier to measure than Q).

Figure 10 provides both the experimentally measured current-voltage curve and the thrust one. In each sub-figures of Figure 10, after the first voltage threshold, the intensity/voltage ratio displays a linear trend associated with a quadratic variation of the intensity with voltage, as predicted by Townsend law [START_REF] Nijdam | The physics of streamer discharge phenomena[END_REF]. At high-voltage, the departure from this linear trend (more visible in Figure 10b than Figure 10a) reveals the streamer-regime leading to current intensity greater than the Townsend regime. A corresponding parallel increase of thrust in the streamer-regime is not clearly observed, however. Since the expected linear relationship between thrust and current suggests a parallel increase to arise, its absence is puzzling. Nevertheless, a similar observation is reported in [START_REF] Xu | Higher thrustto-power with large electrode gap spacing electroaerodynamic devices for aircraft propulsion[END_REF], with the interpretation that the streamer regime arises with reverse corona, the production of negative charges from collector to emitter producing a reverse ionic wind contribution. This might indeed explain the absence of thrust increase associated with the streamer regime. In each sub-figure of Figure 10, the first and second inception threshold are reported with dotted lines. It can be observed that successive inceptions neither lead to visible change in the current intensity slope, nor in the thrust one. Finally comparing Figure 10a and 10b an almost factor three improvement in thrust, can be observed in the configuration L = 45mm which is the closest to collector. Figure 11 provides similar measurements as those reported in Figure 10, but for an increasing number of emitters being closer one another. In each case, the successive ignition thresholds are reported with vertical dotted lines. Again, both current regimes are observed : Townsend's one having a linear I/V behavior versus V , and the streamer one. Albeit one can observe in Figure 11d a concomitance between the transition toward the streamer regime and a specific ignition threshold, this coincidence does not occur in any other configurations. Hence, there is no evidence that the ignition patterns are driving the occurrence of the streamer regime. Each couple of sub-figure (a-b) and(c-d) compare, for the same configuration (same number of emitters and same emitter/emitter distance d), the influence of the emitter/collector distance L. As in Figure 10, L has a strong effect on the current intensity as well as on the thrust. Nevertheless, in figure 11a-b, for which d = 40mm the improvement from L = 100mm to L = 45mm is closer to a factor two than the observed factor three in Figure 10. Hence the emitter closeness levels down the influence of the emitter/collector one, an additional effect of the emitter/emitter electro-static shielding. Furthermore, one should question the influence of the ignition patterns and the number of ignited sources on the observed current and thrust. As previously found in the literature quantitative improvements in the provided intensity current and thrust is not observed for denser emitter sources. Nevertheless, for fully ignited proximal configurations (i.e L = 45mm) the current intensity as well as the thrust performance are also neither degraded for denser emitter array, as for example observed when comparing Figures 11a and11c for V = 35kV giving approximately the same 22mN thrust. Hence, it is found that, when the entire emitter array is fully ignited, having denser emitter sources does not degrade the thrust (i.e the free stream mass-flow). Since previously reported measurements in some multi-emitter/collector configurations have reported performances degradation [START_REF] Wilson | An investigation of ionic wind propulsion[END_REF][START_REF] Masuyama | On the Performance of electrohydrodynamic propulsion[END_REF] in denser emitter arrays, also confirmed by numerical simulations [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], the reported results are noteworthy. The explanation for previous observations might lie in the possible non-ignition of a large fraction of the emitters, conjugated with sources shielding. Although the last effect is also present in fullyignited configurations, it might be compensated by more broadly distributed density charge emission. Figure 11 also displays numerical predictions of the multi-Kaptzov model for intensity and thrust (deduced from electric potential and charge density computations in the same way as described in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF]). Concerning the current/intensity curves since ignition thresholds are within 10 -20% accurate, the current predictions are not precise at low current intensity, although experiments uncertainties are obviously larger there, as expected. Hence, in many cases (e.g. Figures 11a,b,c,d) the slope of the predicted linear regime closely follow the experimental one. Predictions and experiments hardly collapse in this linear regime however, since approximated ignition threshold result in missing its abscissa origin. Furthermore, at higher voltage, since the streamer regime is not included into the modeling, the experiments and numerical predictions differ quite significantly in each case. Surprisingly, the comparison of the thrust predictions obtained without any adjusted parameter with experiments do not suffer much from the lack of streamer-regime consideration in the modeling. Indeed, Figures 11a,b,c,d display a close matching between the observed thrust (black circle) and its prediction (straight line). More precisely the numerical predictions indeed produce a 20% accurate estimate of the observed thrust in the 3E/1C configuration for emitter/collector distance L = 45mm (Fig. 11a), and a 10% accurate one for 3E/1C & L = 100mm, 5E/1C & L = 45mm, 5E/1C &L = 100mm (Fig. 11b,c,d). This observation is consistent with the previously evoked contribution of a reverse corona [START_REF] Xu | Higher thrustto-power with large electrode gap spacing electroaerodynamic devices for aircraft propulsion[END_REF] being responsible for observed additional currents not contributing to bulk charges drift, hence poorly contributing to EAD thrust. Hence, possibly, the streamer regime associated with the deviation of the I/V versus V slope might be also partly attributable to a reverse corona effect. In the 9E/1C configuration at distance L = 45mm (Fig. 11e), the numerical prediction comparison accuracy is degraded to 45% accurate in Figure 11e. This degradation of the prediction is quantitatively consistent with the ignition threshold mismatch already observed in this configuration in the 7th line of Figure 4. This degradation might result from the increased electro-static shielding of emitters in this 9E/1C configuration, resulting in a degraded prediction for charge production, and finally for the produced thrust. As observed for the current intensity, the emitter's successive ignitions do not produce a slope discontinuity of the resulting thrust. This insignificant signature of successive emitter's ignition is consistently observed experimentally and numerically giving a supplementary support to the modeling significance.

We finally consider the relationship between the total emitted UV light integrated over the entire emitter array with the produced current intensity and the resulting thrust. Figure 12 display the total recorded emitted UV light variation with thrust, current-voltage ratio, and current (inset). In this figure, we report the previously discussed regimes, i.e. the linear and the streamer regime with a vertical dotted line reported from the critical value of I/V = 0.63(SI) of Fig. 11d and I/V = 1.74(SI) of Fig. 11e. In the linear regime we indeed found that the current intensity linearly varies (dotted lines) with the current intensity in Figure 12's insets. Similar observations have been found for the other configurations (not shown). Surprisingly, one can observe in Figure 12 a very linear trend between the UV-light intensity and the thrust over the entire range, including the streamer regime (brown points versus brown dotted lines). This observation shows a weak effect of the streamer regime on the UV-light/thrust linear relation. Albeit a positive correlation between emitted UV light and current was expected from the known physical relationship between secondary photo-ionization and unipolar charge ejection, it was not clear how the total thrust will be related to it. Such a simple linear trend clarifies the interpretation of the provided measurements. Furthermore this observation enables to establish the interest of the emitted light as a supplementary relevant diagnosis for ionic wind flow generation and propulsion issues.

Conclusion

This contribution analyses experimentally as well as numerically the multi-inception patterns of emitter array sources in DC-corona discharge. Depending on the considered configurations, i.e collector/emitter interaction dominated (resp. emitter/emitter interaction dominated), central preferential ignition is found (resp. peripheral). A systematic analysis of successive ignition patterns of various configurations has been performed experimentally from direct optical visualization. These ignition patterns are very satisfactorily compared with predictions from a multi-Kaptzov drift-diffusion model solved using a finite-element method. For proximal configurations for which collector/emitter interaction dominates, a full ignition is always found (starting from a central one). On the contrary for large emitter/collector distances, emitter/emitter interaction dominates, resulting in central shielding, peripheral first ignition, and possibly, depending on emitter's density, a lack of full ignition. This contribution thus provides a deeper understanding of complex ignition patterns arising in multi-emitter arrays. It is shown that full ignition might not be prevented by emitter shielding. When full ignition is achieved, emitter shielding however still shrinks unipolar charge emission in central regions. Nevertheless, it is shown that full ignition can be achieved from a proper design of the emitter/collector system, a first condition for multi-emitter sources scalability. this case, the following electrostatic problem is first considered

∇ 2 φ = 0 on Ω, φ = 1 on ∂Ω e , φ = 0 on ∂Ω c , ∇ φ • n = 0 on ∂Ω E , (A.1)
with imposed Dirichlet boundary conditions at emitters and collectors boundaries ∂Ω e & ∂Ω c , and Neumann one (zero electric field in the far field) at external boundaries [START_REF] Hinch | Think Before You Compute[END_REF]. A weak formulation of (A.1) is obtained from multiplying with a test function ϕ † and integrating over domain Ω. Using the divergence theorem so that boundary conditions are taken care of whilst taking into account the zero contribution of ∇ φ • n on Ω E , the weak-formulation of (A.1) leads to

∀ϕ † , ∂Ω e ∪∂Ω c ϕ † ∇ φ • n - Ω ∇ φ • ∇ϕ † = 0, (A.2)
where the integration measure is omitted (a usually done in weak-formulations and in this Appendix with subsequently considered ones). Furthermore, in all weakformulations considered here, a finite-element discrete mesh is considered upon which the discretisation of weak-formulation operators are performed. Automatic mesh raffinement procedures are used nearby emitters so as to obtain the necessary accuracy to capture strong potential gradients. More detailled discussion about mesh raffinement procedures can be found in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF].

We introduce Lagrangian multipliers associated with the local electric fields at emitter and collector, λ e,j = ∇ φ • n on ∂Ω e j , (A. (A.9)

The finite-element dicretization of the linear operator (A.9) involves the construction of a matrix A acting upon the unknown vector q assembling the total degree-of-freedoms of the problem unknowns formulated as the linear problem Aq = b, with

A =   -Ω ∇ϕ † ∇ ∂Ω e ϕ † ∂Ω c ϕ † ∂Ω e λ † e 0 0 ∂Ω c λ † c 0 0   , q =   φ λ e λ c   , b =   0 ∂Ω e λ † 0   ,
The inversion q = A -1 b provides the electric field amplitude at each emitter i from integrating λ e along each emitter boundary ∂Ω e i . The electrostatic problem being linear, all electrostatic properties are proportional to the applied potential V . In particular, considering our non-dimensionalization choices, the electric field E j at the emitter number j can be expressed as E j = Êj V /R c where Êj is the non-dimensional value. The latter can be deduced from the computed non-dimensional electric potential φ by averaging over the corresponding surface, namely: Following the classical Kaptzov assumption, it is expected that ignition occurs as soon as the electric field reaches the critical value given by the Peek law on at least one of the emitters, namely max j∈[1,Ne] (E j ) = E P eek , where .11) This leads to the following prediction for the first ignition threshold : In Ref. [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], Lagrange multipliers representing the ion flux at the ignited emitters were introduced. In the previous section, on the other hand, Lagrange multipliers representing the electric field at non-ignited emitters were used. Both ideas are combined by defining Lagrange multiplier at emitter j as (A.17)

E P eek = 3.1 1 + 0.308 √ R e , ( A
V c = E P eek R c max j∈[1,Ne] ( Êj ) (A.12) Appendix A.
As identified in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], the hyperbolic-dominated nature of the drift-diffusion equation leads to convergence difficulties when associated to Dirichlet conditions at the ignited emitters. A regularization method was introduced in this reference, which consists of replacing the Dirichlet condition by φ = 1 + ε λ ∂ 2 λ e j ∂s 2 (A. [START_REF] Leng | An empirical formula for electrohydrodynamic thrust in the point-to-grid configuration with small electrode gaps[END_REF] where s is the curvilinear coordinate along the boundary and ε λ a small regularization parameter fixed here to 10 -3 . Test functions λ e, † j associated with the generalized Lagrangian multipliers (A. This nonlinear problem is solved iteratively using a Newton algorithm very similar to that used in ref. [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF] which is not detailed here.

Appendix A.3. Ignition pattern selection algorithm

This appendix details the algorithm used to identify the correct ignition pattern for a given configuration, which is also summarized in Table 1.

We introduce the emitter state indicator θ e being a vector in {0, 1} Ne whose components θ e j ∈ {0, 1} gives the ignition state (1 for ignited, 0 for "off") of the j'th emitter with j = 1, ..N e . The computation starts assuming a given ignition pattern, represented by state vector θ e , and solves a set of boundary conditions compatible with this assumption. These boundary conditions being prescribed electric potential at every emitters and collector, and a supplementary Kaptzov assumption applied at ignited emitter only. Kaptzov assumption means that a constant electric field having Peek's value (A.11) is supposed on every ignited emitter. This is why we name this method ""a multi-Kaptzov model".

The validity of the initial emitter state θ e assumption is verified a posteriori. For this, the following criteria is adopted : (a) if an emitter was assumed as "non-ignited" and the computed electric field at this emitter exceeds the value corresponding to the Peek law, then the hypothesis has to be revised as "ignited" ; (b) on the other hand, if an emitter was assumed "ignited" and the computed charge density is negative, then the hypothesis has to be revised as "non-ignited" because this hypothesis leads to non-physical results. Thus, after solving the problem with the Newton method, the validity of the emitter state hypothesis is checked, revising the state vector θ e accordingly. This result in possibly solving again the problem with a new hypothesis, until convergence.

Algorithm 1 Multi-Kaptzov algorithm procedure. Here ε n = 10 -15 . For E P eek we use the usual Peek's formula (A.11) also discussed in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF][START_REF] Franck | Dielectric phenomena in high voltage engineering[END_REF]. Require: State vector of emitters θ e Require: Electric potential value V Require: Emitter configuration : radius R ∈ R Ne , and positions X ∈ R 2Ne vectors Require: (Refined) Mesh T h around emitters 1: Solve harmonic problem (A.9) over emitters configuration 2: while flag= 0 do if E i > E P eek then θ e i = 1, flag= 0

Figure 1 :

 1 Figure 1: Sketch of the numerical domain and the related notations : ∂Ω e I for ignited emitters domain frontier (enlarged for clarity), ∂Ω e O for non-ignited ones (here only one non-ignited emitter), ∂Ω c for collector frontier, ∂Ω E = ∂Ω Eu ∪∂Ω Ei ∪∂Ω E d ∪∂Ω Es for lateral frontier, the domain size being L Ω × H Ω .

R 2 c

 2 ε0V is considered. Dimensionless form of the Poisson equation and charge conservation read ∇ 2 φ = -ρ P e = µV R c /D ρ is the Péclet number based on the ion diffusivity D ρ and positive charge mobility µassumed constant. The corresponding boundary conditions are φ = 1 and ∇ φ • n = Êa on ∂Ω e I , φ = 1 and (ρ∇ φ + 1 P e ∇ρ) • n = 0 on ∂Ω e O , φ = 0 on ∂Ω c , ∇ φ • n = 0 and (ρ∇ φ + 1 P e ∇ρ) • n = 0 on ∂Ω E ,

Figure 2 :

 2 Figure 2: (a) Sketch of the experimental set-up. (b) Focus on the emitter/collector electrode systems under study and definition of the various geometrical parameters. The 'incepted' emitters, here nearby the center, are surrounded by a red halo, to exemplified an inception pattern. Definitions of the various geometrical parameters are given.

Figure 3 :

 3 Figure 3: Processing steps from a raw image to 'incepted' emitters visualization and corona discharge light intensity profile. Example with 3 emitters. (a) Sketch of the setup top view. (b) raw image example of the emitter UV light images, with the averaging frame along y direction in blue, and Region Of Interest (ROI) in dotted red lines (c) ROI image extraction with added symmetrical part (d) y-averaged light intensity profile along x in (c)'s ROI. (This figure has been performed using Inkscape® free software).

Figure 4 :

 4 Figure 4: Sketch of multi-inception patterns for all considered configuration, varying N e from 3 to 9. The charts in the ignition state column are depicting the emitter's ignition pattern, being white for ignited emitters, black for non-ignited ones. Left column : experimental observations, Middle column : numerical predictions using drift-diffusion multi-Kaptzov modeling. Right column : relative error e i = |(V C E,exp -V C E,num )/V C E,exp | for ignition threshold prediction of emitter i (i = 0 for central emitter).

Figure 5 :

 5 Figure 5: Error in the ignition threshold (last column of Figure 4) prediction versus the emitter/emitter distance d. A 1/d trend is represented with dotted lines.

Figure 6 :Figure 7 :

 67 Figure 6: Multi-inception for N e = 3 emitter system for distances L = 45mm and d = 60mm. (a) Experimental multi-inception UV image. (b) Experimental profile of corona discharge light-intensity for various applied voltage at emitters.

Figure 9 :

 9 Figure 9: Electric charge density ρ [C.m -3 ] predicted by the multi-Kaptzov simulation, for configurations (a) 3E/1C L = 45mm d = 60mm, V = 13kV; (b) 3E/1C L = 100mm d = 40mm, V = 17kV; (c) 5E/1C L = 45mm d = 20mm, V = 19kV; (d) 5E/1C L = 100mm d = 20mm, V = 22.5kV; (e) 5E/1C L = 45mm d = 20mm, V = 30kV; (f) 5E/1C L = 45mm d = 20mm, V = 57kV;

  (a) 3E/1C, L = 45mm, d = 60mm (b) 3E/1C, L = 100mm, d = 60mm

Figure 10 :

 10 Figure 10: Thrust (black circle symbols) and current intensity (diamond brown symbols) versus applied potential for configuration 3E/1C, for emitter spacing d = 60mm when varying the distance to L (Cf Fig. 2) between emitter/collector. The blue dashed line provides the quadratic trend I ∼ V (V -V c E0 ) (affine relation between I/V and V ).

  (a) 3E/1C L = 45mm, d = 40mm (b) 3E/1C L = 100mm d = 40mm (c) 5E/1C L = 45mm d = 20mm (d) 5E/1C L = 100mm d = 20mm (e) 9E/1C L = 45mm d = 10mm

Figure 11 :

 11 Figure 11: Experimental Thrust (black circle symbols) and current intensity (diamond brown symbols) versus applied potential for two configurations (i.e L = 45mm & L = 100mm Cf Fig. 2) when densifying emitters sources (i.e diminishing parameter d). For each sub-figure continuous curves represent the prediction of the multi-inception Kaptzov model presented in section 2.1. As in Fig. 10 the blue dashed line provides the quadratic trend I ∼ V (V -V c E0 ). Numerical predictions are represented with continuous curves for the thrust prediction (continuous black line) and the current intensity (continuous brown line)

  (a) 5E/1C L = 100mm d = 20mm (b) 9E/1C L = 45mm d = 10mm

Figure 12 :

 12 Figure 12: Total UV light recorded from emitter's array versus current, thrust and current/voltage ratio in two configurations. (a) Same configuration as Fig. 11d (b) same configuration as Fig. 11e.

3 )

 3 λ c,j = ∇ φ • n on ∂Ω c j . (A.4)Usage of the Lagrange multipliers allows to introduce an associated set of test functions λ † e,j and λ † c,j which are used to enforce the Dirichlet boundary conditions of the potential ∀λ † e,j , the set of Lagrangian multipliers and the set of associated test functions by introducing notation

  multiplier λ I,j and non-ignited Lagrangian multiplier λ O,j defined as λ e I,j = ρ ∇ φ + 1 P e ∇ρ • n on ∂Ω e j , (A.16) λ e,c O,j = ∇ φ • n on ∂Ω e,c j .

  [START_REF] Krishna | A laser-microfabricated electrohydrodynamic thruster for centimeter-scale aerial robots[END_REF] are introduced to enforce the boundary conditions. The complete problem (1)-(2) associated with boundary conditions (3) fulfills the following variational formulation ∀ϕ † ,

4 : 5 : 9 :

 459 Solve non-linear drift problem (A.19) over emitters configurationfor i ∈ {1, N e } do if n i < ε n then θ e i = 0, flag= 0

  2. Drit-diffusion modeling of multi-ignition

	∂Ω e 0 =	∂Ω e j .	(A.14)
	{j∈[1,Ne]/θ e j =0}
	This section now considers the weak-formulation of drift-diffusion problem. Since
	we now consider situations above the first ignition onset (i.e above the first voltage
	threshold associated with the first emitter ignition) multi-ignition can occur. In this
	general context, the set of emitters boundaries ∂Ω e can then be split into two distinct
	sub-groups ∂Ω e I and δΩ e O , such as ∂Ω e = ∂Ω e I ∪ ∂Ω e O in order to distinguish the groups
	of ignited emitters (denoted by index I) from non-ignited ones (denoted by index O
	for "off") such that
	∂Ω e I =	∂Ω e j ,	(A.13)
	{j∈[1,Ne]/θ e j =1}

Acknowledgments

This work as been supported by the RTRA STAE Grant 'IOWIND' number R010-L00-T00/2019-CIF-R-103.

Appendix A. Details about the multi-ignition method Appendix A.1. Harmonic modeling of multi-ignition .

In this first section, the modeling of the first emitter ignition is considered for the electrostatic problem in the abscence of any bulk charge within the air domain. In