Nicolas Monrolin
email: nicolas.monrolin@enac.fr

Franck Plouraboué
email: franck.plouraboue@imft.fr

Multi-scale two-domain numerical modeling of stationary positive DC corona discharge/drift-region coupling

Keywords:

Ionic wind

 result in flux continuity conditions at the boundary of the two domains. These coupling conditions are enforced by Lagrange multipliers, within a variational formulation, leading to a hierarchy of non-linear coupled problems. The proposed approach is both monolithic and two-domains: two asymptotic regions, an inner-one associated with corona discharge, and an outer-one, the ion drift region. Numerical convergence and validations of the finite element implementation is provided. A comparison with various experimental results convincingly demonstrate the applicability of the method, which avoids tuning parameters dedicated to each specific configuration, but, on the contrary, exclusively relies on known and measurable physical quantities (e.g., ion mobilities, photo-ionization coefficient, ionization electric field, Townsend discharge coefficient, etc...).

Introduction and context

DC-corona discharge is a complex phenomenon arising within a gas when the electric field reaches a threshold for which electron collisions cascade and produce positive and negative ion charged molecules in some confined regions. These confined regions are called 'corona' or 'glowing regions' where a cold plasma is set-up and ejects unipolar charges in a second region called the 'drift region' in the gas where electrons die-away. Since these unipolar charges can further collide with neutral gas molecules in the 'drift region', under the action of an applied electric field, they can then generate net momen-tum and produce ionic wind there. A back-coupling between these two regions comes from the action of photo-ionization. Light is indeed emitted from the 'glowing region' into the 'drift region' and produces a small amount of secondary electrons in a thin zone of the 'drift region', of crucial importance to sustain the cold plasma creation. This very brief and synthetic description of DC-corona discharge depicts its complexity, so that its modeling raises challenges. If one adds the fact that the time-scales associated with charge creation and electro-drift can be very different, one realizes that the dynamics of corona, (e.g. associated with so-called streamers), is even more challenging [START_REF] Pancheshnyi | Numerical simulation of filamentary discharges with parallel adaptive mesh refinement[END_REF][START_REF] Soria-Hoyo | A PIC based procedure for the integration of multiple time scale problems in gas discharge physics[END_REF][START_REF] Papadakis | New high quality adaptive mesh generator utilized in modelling plasma streamer propagation at atmospheric pressures[END_REF][START_REF] Li | Spatially hybrid computations for streamer discharges: II. Fully 3D simulations[END_REF][START_REF] Villa | An asymptotic preserving scheme for the streamer simulation[END_REF][START_REF] Villa | Stability of the discretization of the electron avalanche phenomenon[END_REF][START_REF] Villa | An efficient algorithm for corona simulation with complex chemical models[END_REF]. Furthermore the detailed physics of the modeling associated with the various non-stationary aspects of corona render its comparison with experimental results (e.g. the so-called Trichel pulses) delicate, either using commercial codes [START_REF] Tran | Numerical modelling of negative discharges in air with experimental validation[END_REF] or more elaborated ones [START_REF] Villa | Simulation of the AC corona phenomenon with experimental validation[END_REF], albeit feasible in 2D [START_REF] Zhu | Nanosecond surface dielectric barrier discharge in atmospheric pressure air: I. Measurements and 2D modeling of morphology, propagation and hydrodynamic perturbations[END_REF]. Nevertheless, at intermediate voltages, above the inception voltage, a steady-state can be sustained, the modeling of which is still difficult when coupled with drift-region. Here, we focus our interest on the numerical computation of steady-state DC-corona discharge which is already a difficult issue, as, for example, studied in [START_REF] Feng | Application of Galerkin finite-element method with Newton iterations in computing steady-state solutions of unipolar charge currents in corona devices[END_REF] for the drift region or in [START_REF] Hasan | DC negative corona discharge in atmospheric pressure helium: transition from the corona to the 'normal' glow regime[END_REF][START_REF] Benilov | Multiple solutions in the theory of dc glow discharges and cathodic part of arc discharges. Application of these solutions to the modeling of cathode spots and patterns: a review[END_REF][START_REF] Bieniek | Modelling cathode spots in glow discharges in the cathode boundary layer geometry[END_REF][START_REF] Ferreira | Simulation of pre-breakdown discharges in high-pressure air I: the model and its application to corona inception[END_REF] for the corona region. From the applicative view-point, the modeling of steady corona discharge is relevant in many applications such as electrostatic precipitors [START_REF] Davis | Wire-duct precipitator field and charge computation using finite element and characteristics methods[END_REF], EHD (Electro-Hydro-Dynamic) gas pump [START_REF] Chang | Narrow-flow-channel-driven EHD gas pump for an advanced thermal management of microelectronics[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], 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] and xerography, i.e. electrophotography. In these applications, many configurations involve corona discharges generated from wires into a cavity, the wall of which are placed at reference potential. In these cases the Kaptzov assumption (which is correct for a wire in a infinite domain, or centered into an axi-symmetric cavity) might oversimplify the real electric field at emitters, so that a more elaborated approach taking care of the corona discharge physics is necessary.

Historically, many approaches have tried to avoid the modeling of the complete coupling between glowing region and drift region. Most of these approaches relied on experimental measurements, providing some approximate expression of the electric field and the charge density at the edge of the glowing region. More precisely, these approaches are generally calibrated for air at atmospheric pressure, and provide the current-potential law Iφ needed to set the charge distribution and the electric field at the frontier between glowing and drift regions. For a single cylindrical electrode (called the emitter), inside a finite co-cylindrical geometry Townsend's law has been successfully used [START_REF] Townsend | Electricity in Gases[END_REF][START_REF] Roth | Industrial Plasma Engineering, Principles[END_REF][START_REF] Shrimpton | Charge Injection Systems[END_REF]. Considering non axisymmetric drift region problems whilst using axi-symmetric charge injections and/or electric potential (such as Peek's law) has also been used (e.g. in point/plane configuration [START_REF] Adamiak | Simulation of corona discharge in point-plane configuration[END_REF], cylinder/cylinder configurations [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], etc...) which might be a fair approximation is some cases. Nevertheless, in general non axi-symmetric configurations not only the parameters of currentpotential law (and/or charge injection-electric field law) have to be adapted, but also the hypothesis of axi-symmetrical emitted charges has to be reconsidered. For example, based upon experimental measurements [START_REF] Chapman | Corona point current in wind[END_REF][START_REF] Vogel | Experimental evaluation of discharge characteristics in inhomogeneous fields under air flow[END_REF] have shown that the current-potential law is modified in the presence of external air flow in the drift region in a tip/plane configuration. More recently, the modification of the charge injection boundary conditions has also permitted to reproduce experimental measurements in a point-to-ring configuration [START_REF] Guan | Experimental and numerical investigation of electrohydrodynamic flow in a point-to-ring corona discharge[END_REF]. Other experimental evidences calls for non axi-symmetric charges injections, such as the observations of light intensity variations (in Dielectric Barrier Discharge, i.e DBD, configurations) resulting from the gas flow effects [START_REF] Pereira | Effect of external flow velocity on momentum transfer of dielectric barrier discharge plasma actuators[END_REF], as recently confirmed by [START_REF] Moreau | Ionic wind produced by positive and negative corona discharges in air[END_REF]. In this context [START_REF] Nguyen | Computational study of glow corona discharge in wind: biased conductor[END_REF] has recently proposed to use a Robin boundary condition for the charge density n injection at the drift region edge, , n(x)| = β(E(x)| -E p ).

The boundary condition associated with the drift region is clearly resulting from the interaction between the various fields (electric potential, ions, electrons) between the glowing and drift regions. This is why many modeling approaches have considered a coupled multi-domain or 'hybrid' approaches in order to model the physics of DC corona [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF][START_REF] Seimandi | A two scale model of air corona discharges[END_REF][START_REF] Wettervik | A domain decomposition method for three species modeling of multi-electrode negative corona discharge -with applications to electrostatic precipitators[END_REF][START_REF] Chuan | A hybrid approach for corona discharge in needle electrode configuration: in a large-scale space[END_REF]. A major issue in this area is to foresee a relevant modeling using physical parameters only, (kinetically based parameters available from open data-bases) but avoiding the need of dedicated phenomenological parameters. Efforts toward this direction have been addressed using multi-domain approaches within a partitioned strategy, iteratively seeking for the solution in each sub-domain with a fixed point method. Nevertheless, in many problems a similar partitioned strategy is known to be less stable than a monolithic one. Monolithic fully-coupled approaches have also indeed been pursued to numerically compute the non-linear elliptic/hyperbolic problem associated with electric potential, electrons, ions charge creations, electro-drift and secondary photo-ionization (Cf [START_REF] Montijn | An adaptive grid refinement strategy for the simulation of negative streamers[END_REF][START_REF] Villa | An efficient algorithm for corona simulation with complex chemical models[END_REF] among others). These monolithic fully-coupled approaches might be interesting in order to get physically detailed, chemical composition of corona [START_REF] Villa | An efficient algorithm for corona simulation with complex chemical models[END_REF]. They have been mainly applied to very simple corona geometries, since the numerical complexity of the complete physics is difficult to address in complex domains.

In this paper we propose an alternative method both monolithic and two-domain, derived from the asymptotic analysis of the fully-coupled problem, producing two asymptotic regions, an inner-one associated with corona discharge, and an outer-one, the drift region. This method generalizes the analytical axi-symmetrical analysis performed in [START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF] to domains having any regular shapes for which no analytical solution is available. The approach combines the advantage of being stable and efficient so as to be able to address potentially complex domains in the drift region. The numerical approach is also inspired by domain decomposition techniques [START_REF] Babuška | The finite element method with Lagrangian multipliers[END_REF][START_REF] Magoulès | Lagrangian formulation of domain decomposition methods: a unified theory[END_REF] using Lagrange multipliers defined at the interface between two domains to match suitable boundary conditions between the various fields involved. The idea behind our approach is to gain understanding on the corona discharge mechanisms so as to set-up an asymptotic hierarchy of main coupled effects, whilst retrieving irrelevant ones.

The paper is organized as follows. Section 2.1 describes the constitutive model, its underlying physics, the geometrical setting and context as well as its dimensionless formulation. Section 3 discusses its asymptotic formulation and develops on the resulting multi-scale/two-domain strategy. Section 4 provides the numerical details of the implementation, the convergence study, and validation test-cases combining previous analytic, numerical and experimental results. Finally section 5 showcases some illustrations and comparison with previously published experimental results.

Corona discharge model

Constitutive equations

We consider a positive DC-corona discharge arising into an infinite two-dimensional configuration sketched in Fig. 1a. Even though, the general ideas and method proposed in this paper might be generalized to 3D, there are hereby distinctly derived in 2D for notations and methodological simplifications.

As mentioned in the introduction, the effective fluid model of the positive DC corona is considered. The production of positive ions, electrons and negative ions (respective density n p , n e and n n ) is governed by the impact ionization coefficient α and the attachment coefficient η. The ionization coefficient dependency with electric field follows the standard Townsend form α = β exp(-E i /E), [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] where β and E i are two physical parameters which depends on the gas composition, thermodynamic conditions and they are supposed to be known. E i is the ionization electric field, i.e., the field beyond which the corona discharge lightens. The impact ionization coefficient α is assumed to vanish at low electric field intensity E = ∇ϕ . The complete set of equations describing the electric potential ϕ, electron density n e , positive and negative ion charges densities n p and n n is

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ 2 ϕ = e 0 (n e + n n -n p ), ∇ • j p = α j e + S, ∇ • j e = (α -η) j e + S, ∇ • j n = η j e , ( 2
)
where e is the elementary charge, j e = μ e n e ∇ϕ, j p = -μ p n p ∇ϕ, and j n = μ n n n ∇ϕ are the local fluxes of the electron, positive and negative ion charges, j e = |j e |, associated with their respective mobility (i.e. μ e for the electrons, μ p , μ n for the positive and negative ion charges). α is the impact ionization Townsend coefficient (1) and η is the attachment coefficient. In the following we introduce notation

α ef = α -η, (3)
and consider that α ef is a known smooth function of E. Furthermore, Appendix C shows that both coefficients α and η (and thus α ef ) have an exponential dependence with the inverse of the local electric field, similar to (1) that we will be subsequently used. Finally, S is the source term associated with secondary ionization which is one complex aspect of corona discharge modeling. Note that, in this formulation, photo-ionization provides an equally balanced source term for electrons and positive charges, since it both generates an electron and a positive charge out of a neutral molecule. Such a balance is not always taken into account, but this point will be discussed further in the next section. Even if secondary ionization is very small compared to the impact ionization, it is necessary to explain the onset and to sustain the discharge. Photoionization is the source of secondary electron and results from a non-local creation coming from a convolution of the charge flux with a radiative kernel. In 3D, using position vector R = r + ze z built from horizontal position r and vertical distance along z,

S(R) = γ g(R, R ) (α(R ) -η(R )) j e (R )d 3 R , (4)
where, again, j e (r ) = |j e (r )| and S(R) is the number of photo-ionizing events at position R per unit time and volume.

The coefficient γ is the secondary electron efficiency, identical to the one introduced by Zheng [START_REF] Zheng | Self-sustained criterion with photoionization for positive dc corona plasmas between coaxial cylinders[END_REF], a dimensionless small quantity, i.e. γ 1, to account for the photo-ionization cross-section and probability as in [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF][START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF]. The photon radiative kernel may have different forms [START_REF] Naidis | Conditions for inception of positive corona discharges in air[END_REF][START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF][START_REF] Penney | Photoionization measurements in air, oxygen, and nitrogen[END_REF][START_REF] Janalizadeh | A framework for efficient calculation of photoionization and photodetachment rates with application to the lower ionosphere[END_REF]. We hereby derive a general theory which can be adapted to any (regular) form of kernel. In this paper, we restrict our attention to 2D problems being translationaly invariant along z. In this context, we derive in Appendix B a specific 2D kernel from a well-established 3D one [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF][START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF]. Hence, in the hereby considered context cylindrical coordinates are used and (4) reduces to S(r) = γ G(r, r ) (α(r )η(r )) j e (r )d 2 r , [START_REF] Papadakis | New high quality adaptive mesh generator utilized in modelling plasma streamer propagation at atmospheric pressures[END_REF] where the photo-ionization source S(r) also being invariant along z, because it only depends on r. For the sake of simplicity, in the following, each time we will specify the 2D domain of integration, we will omit the differential increment in all integrals, i.e. (5) will be denoted S(r) = γ G(r, r ) (α(r )η(r )) j e (r ). [START_REF] Li | Spatially hybrid computations for streamer discharges: II. Fully 3D simulations[END_REF] The boundary conditions associated with problem (2) are based upon notations of Fig. 1a. The electric potential ϕ fulfills Dirichlet boundary conditions on the electrodes, with a high tension ϕ a applied at emitter and a reference zero potential at collector, i.e.

ϕ| ∂ e = ϕ a , ϕ| ∂ c = 0. (7)
Both n p and n e fulfill a purely hyperbolic problem so that one upstream boundary condition for each field is needed. In a positive corona discharge, for positive charges traveling along the electric field lines from the emitter surface ∂ e toward the collector one ∂ c , zero positive charges flux is set at the emitter

j p • n| ∂ e = 0. (8)
Symmetrically, for the electrons and negative charges traveling against the electric field, a zero flux inlet boundary condition is set at the collector j e • n| ∂ c = 0, [START_REF] Villa | An efficient algorithm for corona simulation with complex chemical models[END_REF] so that it is assumed that no electrons are injected at collector, which might result in no electrons at all. But a few secondary electrons are created by photo-ionization nearby the emitter that will feed the corona discharge. This simplified framework is meaning full since we assume that the electric field at the emitter is much larger than the one nearby collector so that the generated photo-ionization source term is much smaller there, thus of negligible effect.

Dimensionless formulation

The first main physical parameters associated with the corona discharge are the applied electric potential difference ϕ a , between the emitter and the collector, being at distance L apart, with a resulting applied electric field magnitude of ϕ a /L. It is interesting to compare this applied field to the "internal" one defined by the electric ionization field E i used in Townsend relation [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF]. From this comparison a small asymptotic parameter ε is defined as in [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF][START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF]

ε = ϕ a L E i . ( 10
)
This ratio being small indicates that the applied electric field is small compared to the ionization field of the discharge. Dimensionless variables are chosen from the external (outer or drift region) length reference L by

r = r L , φ = ϕ ϕ a , nk = n k n k , â = a L , ( 11
)
with k ≡ e, p, n for electrons, positive ions and negative ions respectively and a the emitter radius. The reference number

density n k is n k = 0 ϕ a eL 2 μ p μ k . ( 12
)
Note that contrary to Durbin & Turyn [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] we differentiate the adimensionalization for ions and electrons so that ne ∼ O [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] in the corona region 1 and np ∼ O [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] in the drift region 2 . This is why the small parameter δ μ = n e /n p = μ p /μ e lateron appears in [START_REF] Shrimpton | Charge Injection Systems[END_REF]. δ μ typically takes values smaller than 10 -2 in air. In the following, we also use inner (corona region)

variable scaling R = r L ≡ 1 r, (13)
so that r = R. Defining outer non-dimensional gradient ∇ ≡ ∂ r, and inner ones as ∇R ≡ ∂ R also lead to ∇ = 1 ∇R . It is interesting to mention that, here, the chosen reference number density n k differs from previous contributions [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF][START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF] since it does not contain the electric current I . This choice is justified because the total current I is a priori unknown, but was taken as control parameter for easier theoretical derivations in [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF][START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF]. Since the purpose of this contribution is to provide a numerical formulation based on known imposed parameters, the current being one result of the computation, we built n k on known parameters. Doing so, the non-dimensional equation for the electric potential, will not contain unknown parameter (such as dimensionless current denoted J in [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF][START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF]). Using dimensionless electric field in (1), as in [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF], the reaction coefficients scale as follows

α ε = Lα = β ε exp - 1 ε Ê , ( 14
) η ε = Lη, ( 15
) αef ε = Lα ef ≡ L(α -η), ( 16
)
with β = β Lε, and Ê = | ∇ φ|. Both α and η dimension being the inverse of a reference length-scale, ( 14)-( 16) state that, this length-scale is the inner one L. For the sake of brevity, in the following we use αef = αη as the effective ionization coefficient.

Let us now consider the non-dimensionalization of the photo-ionization term [START_REF] Papadakis | New high quality adaptive mesh generator utilized in modelling plasma streamer propagation at atmospheric pressures[END_REF]. First, it is important to mention that, since the convolution integral arises over = 1 ∪ 2 , it can be decomposed into two distinct contributions from the corona discharge domain 1 and the drift domain 2 . In these contributions, since the reference length-scale is L in 1 (resp. L in 2 ), the electric field respectively scales as E = ϕ a L Ê in 1 (resp. E = ϕ a L Ê in 2 ). Thus, using previously defined non-dimensionalization and particularly ( 14)-( 16) in ( 5) leads to

S(r) = μ e n e ϕ a L 2 γ 2 ⎛ ⎜ ⎝ 1 G(r, r ) αef (r ) ˆje (r )d 2 r + 2 G(r, r ) αef (r ) ˆje (r )d 2 r ⎞ ⎟ ⎠ . ( 17
)
Then, one needs to consider the non-dimensionalization of the hereby considered 2D photo-ionization kernel G. In most contributions, photo-ionization kernels g(R) are discussed and defined in 3D, with R 2 = |rr | 2 + z 2 the 3D Cartesian distance, z being the direction orthogonal to the hereby considered plane. As detailed in Appendix B, the relation between

g(R) and G(r, r ) ≡ G(|r -r |) being G(|r -r |) = R g(R) 4π R 2 dz. ( 18
)
Then, non-dimensionalization of kernel g(R) leads to g(R) = ĝ(R)/L (Cf Appendix B for more details), and from ( 18)

G(|r -r |) = 1 L 2 R ĝ 4π R2 dẑ = 1 L 2 Ĝ(|r -r |). (19)
From using [START_REF] Chang | Narrow-flow-channel-driven EHD gas pump for an advanced thermal management of microelectronics[END_REF] in ( 17) leads to

S(r) = μ e n e ϕ a L 2 γ ⎛ ⎜ ⎝ ˆ 1 Ĝ(r, ε R ) αef (R ) ˆje (R )d 2 R + 1 ˆ 2 Ĝ(r, r ) αef (r ) ˆje (r )d 2 r ⎞ ⎟ ⎠ , ( 20
)
where we have now re-scaled coordinates in the corona using inner variable R (13), and defining ˆ 1 being a dimensionless (order O (1)) domain 1 . Now, realizing that the second term of (20)'s r.h.s. is small because both the Townsend coefficient α and the attachment term η decay as exp(-1/ ) in region ˆ 2 , so does αef from ( 3) and ( 14)-( 16), dominating over any algebraic power in , one gets,

S(r) = μ e n e ϕ a L 2 γ ⎡ ⎢ ⎣ ˆ 1 Ĝ(r, ε R ) αef ( R ) ˆje ( R )d 2 R + O exp(-1/ ) ⎤ ⎥ ⎦ . ( 21
)
So that, one can then define the non-dimensional photo-ionization kernel Ŝ from S =

μ e n e ϕ a L 2 γ Ŝ, i.e. Ŝ(r) = ˆ 1 Ĝ(r, ε R ) αef ( R ) ˆje ( R ) ≡ ˆ 1 G(r, ε R ) αef ˆje R . ( 22
)
Then, a multipole asymptotic expansion of ( 22), together with the form of ( 19) reads,

Ŝ(r) = G(r) ˆ 1 αef ˆje R + ε∇G(r) • ˆ 1 αef ˆje R R +O ( 2 ), (23)
Ŝ(r) = Ŝ0 (r)+ ε Ŝ1 (r) +O ( 2 ), ( 24
)
neglecting quadrupolar O ( 2 ) corrections. Using reference charge density [START_REF] Zhu | Nanosecond surface dielectric barrier discharge in atmospheric pressure air: I. Measurements and 2D modeling of morphology, propagation and hydrodynamic perturbations[END_REF], outer dimensionless variable r (13) in ( 2)

whilst using non-dimensionalization ( 21), leads to the following dimensionless drift region formulation It is interesting to note that the non-dimensionalization leading to (25) produces a smaller contribution of electron density compared to positive charge in the drift region. The main reason is based on flux considerations: the electron current density at the emitter should balance the ion current density at collector. The ratio between the maximum number density of unipolar positive ions n p and the maximum number density of electrons n e is then given by the mobility ratio δ μ . One might question this hierarchy in the corona region 1 since the ion number density decreases drastically near the emitter surface: np ∼ O [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] in drift region but np = 0 at the emitter whilst ne ∼ O [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] in the corona region and ne = 0 at the collector. In practice this is not a concern since in the corona region, a re-scaling of the coordinates produces O ( 2 ) small term in front of (25)'s r.h.s., leading to negligible charge effect at leading order in the electrostatic problem [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF]. In a nutshell, the space charge plays an important role only in the drift region and is strongly dominated by the positive ions charge, there. Last but not least, it is important to realize that the negative charges concentration do not play an active role in the problem. First, in the corona region, negative charges does not contribute to the potential (as any other charges), for the aforementioned reason of having a negligible impact on electrostatic problem [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF]. Built into the corona region by attachment coefficient η from electron flux, negative charges only migrate to the emitter so as to produce, together with the electrons, the necessary (negative) charge flux balance to the positive charges drifting away from it. Secondly, in the drift region, the only source term for negative charges in [START_REF] Chapman | Corona point current in wind[END_REF] is the product of attachment coefficient η with electron flux. As discussed just after [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], η decay as exp(-1/ ) in the drift region, leading to negligible production of negative charges flux, thus leading to negligible negative ion number density there. This is why, in the sequel, negative charges are not considered.

∇2 φ = -(n p -δ μ ne -nn ), ( 25
) ∇ • ĵp = α ˆje + γ Ŝ(r), (26)
∇ • ĵe = α -η ˆje + γ Ŝ(r), (27)
∇ • ĵn = η ˆje , (28)
Dimensionless problem ( 25)-( 28) is complemented with dimensionless boundary conditions

φ| ∂ ˆ e = 1, ϕ| ∂ ˆ c = 0, ( 29
) ĵp • n| ∂ ˆ e = 0, (30)
and

ĵe • n| ∂ ˆ c = 0. (31)
Thus ( 25)-( 28) associated with boundary conditions ( 29)-( 31) and source term [START_REF] Townsend | Electricity in Gases[END_REF] represents a coupled non-linear nonlocal system of equations. In the following we show how a multi-scale approach can be used to transform it into two coupled local problems, with notations provided in Fig. 1b.

Multi-scale asymptotic expansion

We now seek for a regular asymptotic expansion with respect to parameter of the problem, neglecting O (δ μ ), O ( exp(-1/ ) ) as well as O ( 2 ), but keeping O (γ ) and O ( ) terms, i.e.

( φ, np , ne , α, αef ) = ( φ0 , n0 p , n0 e , α0 , α0 ef ) + ε( φ1 , n1 p , n1 e , α1 , α1 ef ) + O 2 , δ μ , exp(-1/ ) . ( 32
)
We also subsequently define α0 ef ≡ αef ( Ê0 ), α0 ≡ α( Ê0 ), whilst, obviously, Ên = | ∇ φn | for n = 0, 1. Furthermore, from Taylor expanding the electric field expansion Ê = Ê0 + Ê1 + O ( 2 ) in ( 14), leads to α = α0 + α1

+ O ( 2 ) with α1 = α0 Ê1 ( Ê0 ) 2 , α1 ef = ∂ α0 ef ( Ê0 ) ∂ Ê = α1 - ∂ η0 ∂ Ê . ( 33
)
Since from (3), αef = αη, whilst also using notation η0 ≡ η( Ê0 ). Some explicit relation for η(E) and its derivative are given in (C.2) and (C.3). In the following, we will index the fields φ, Ê, ne , np by j , j = 1, 2 for specifying into which domain they fall under.

Corona domain 1 problem

At leading order, the corona problem reads

∇2 R φ0 1 = 0, ( 34
) ∇R • (n 0 p 1 ∇R φ0 1 ) = -α0 n0 e 1 Ê0 1 , ( 35
) ∇R • (n 0 e 1 ∇R φ0 1 ) = α0 ef n0 e 1 Ê0 1 . ( 36
)
Note that, surprisingly, there is no more source term on the right-hand-side of [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF], as opposed to many other two-region modeling for corona models already proposed in the literature, (e.g. [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF][START_REF] Seimandi | A two scale model of air corona discharges[END_REF]), some of them not derived from asymptotic considerations [START_REF] Wettervik | A domain decomposition method for three species modeling of multi-electrode negative corona discharge -with applications to electrostatic precipitators[END_REF]. This issue is much more benign than what could be though at first sight. As a matter of fact, since [START_REF] Seimandi | An asymptotic model for steady wire-to-wire corona discharges[END_REF] is expressed in internal variable R which is stretched upon the external one, r, R = r/ , the resulting re-scaling of the Laplacian applied on the right-hand-side of (25) multiplies it by an O ( 2 ) term. This means that the charge effect on the corona region only adds a very small correction to the potential. Furthermore, taking into account this correction when discarding other O ( 2 ) terms associated with the coupling between 1 and 2 is not asymptotically consistent. At order O ( ), we have

∇2 R φ1 1 = 0, ( 37
) ∇R • (n 0 p 1 ∇R φ1 1 + n1 p 1 ∇R φ0 1 ) = -α0 n0 e 1 Ê1 1 -α0 n1 e 1 Ê0 1 -α1 n0 e 1 Ê0 1 , ( 38
) ∇R • (n 0 e 1 ∇R φ1 1 + n1 e 1 ∇R φ0 1 ) = α0 ef n0 e 1 Ê1 1 + α0 ef n1 e 1 Ê0 1 + α1 ef n0 e 1 Ê0
1 .

(39)

Drift domain 2 problem

For the potential and positive charges in the drift domain, at leading order, the electrostatic (25) and positive charges conservation problem [START_REF] Adamiak | Simulation of corona discharge in point-plane configuration[END_REF] reads

∇2 φ0 2 = -n 0 p 2 , ( 40
) ∇ • (n 0 p 2 ∇ φ0 2 ) = γ Ŝ0 , ( 41
) because the α term is O ( exp(-1/ ) ) in 2 . At order O ( ), we have ∇2 φ1 2 = -n 1 p 2 , ( 42
) ∇ • (n 0 p 2 ∇ φ1 2 + n1 p 2 ∇ φ0 2 ) = γ Ŝ1 , . (43)
Finally, in the following, we will not solve the electron problem in the drift domain 2 , but, for now, we leave it as in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF], but for neglecting the contribution of the αef term which is O ( exp(-1/ ) ), without expanding it in , i.e.

∇ • (n

e 2 ∇ φ 2 ) = γ Ŝ(r). (44)
Since the photo-ionization term Ŝ(r) is evanescent, i.e. exponentially decaying along r from [START_REF] Townsend | Electricity in Gases[END_REF], so does the electron density in the drift region. Hence, except for a small evanescent region of width λ, i.e., a very thin layer λ/L in dimensionless

In the following, we will use ∂ n ≡ ∇ • n for the projection of gradient operator to the outward normal of a boundary.

This leads to

n0 p ∂ n φ0 | ∂ e 1 = 0, (51)
and,

n1 p ∂ n φ0 | ∂ e 1 + n0 p ∂ n φ1 | ∂ e 1 = 0. ( 52
)
In 1 , no further condition is needed. To enforce ion flux continuity, 2 must be fed with the ion flux coming from 1 . This leads to a "one-way" coupling, i.e. np 2 directly depends on np 1 but not reciprocally

ĵp 2 • n 2 | = -ĵp 1 • n 1 | , ( 53
)
with again a minus sign because of the normal. So the inlet condition of boundary condition of 2 is given by 1 and no outlet condition is required.

• Symmetrically, for the electrons and negative charges traveling against the electric field, the inlet boundary condition is set at the collector:

ĵe • n 2 | ∂ c 2 = 0, (54)
leading to,

n0 e ∂ n φ0 | ∂ c 2 = 0, (55)
and

n1 e ∂ n φ0 | ∂ c 2 + n0 e ∂ n φ1 | ∂ c 2 = 0. (56)
We assume that no electrons are injected in at the collector, which should result in no electrons at all. But a few secondary electrons are created by photo-ionization in 2 that will feed 1 through the interface

ĵe 1 • n 1 | = -ĵe 2 • n 2 | . ( 57
)
This is again a "one-way" coupling, since ĵe 1 directly depends on ĵe 2 and not vice-versa. Furthermore, photo-ionization in drift domain 2 depends on the ionization rate in 1 , in a rather complex way. Thus, given [START_REF] Townsend | Electricity in Gases[END_REF] in [START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF] in 2 domain, leads to

∇ • ĵe = γ G(r)M 0 + ε∇G(r) • M 1 + O ( 2 ) , ( 58
)
with the multi-polar expansion associated with hereby defined mono-polar scalar M 0 and dipolar vector M 1

M 0 = 1 αef ˆje R , (59) M 1 = 1 αef ˆje R R , ( 60
)
whilst, again, omitting the differential increment in the integrals. Inserting expansion [START_REF] Moreau | Ionic wind produced by positive and negative corona discharges in air[END_REF] in (59), one finds

M 0 = M 0 0 + εM 1 0 + O ( 2 ), with (61) M 0 0 = 1 α 0 ef ˆj0 e R , ( 62
) M 1 0 = 1 α1 ef n0 e 1 Ê0 1 + α0 ef n1 e 1 Ê0 1 + α0 ef n0 e 1 Ê1 1 R . ( 63
)
And, similarly, inserting expansion [START_REF] Moreau | Ionic wind produced by positive and negative corona discharges in air[END_REF] in (60) keeping only the leading order contribution to the dipolar correction,

M 0 1 = 1 α0 ef ˆj0 e R R .
(64)

Then, (58) reads,

∇ • ĵe = γ G(r)M 0 0 + ε G(r)M 1 0 + ∂ r G(r)M 0 1 • e r + O ( 2 ) . ( 65
)
Realizing from (32) that the electron flux ĵe follows the same regular asymptotic expansion ĵe = ĵ0 e + ε ĵ1 e + O (ε 2 ).

(66)

We seek to solve, at each order, the electron flux coming from photo-ionization only. At leading order in , the forcing term displays an axi-symmetrical radial dependence,

∇ • ĵ0 e = γ G(r)M 0 0 . ( 67
)
The solution of (67) ĵ0 e can be decomposed into a general (conservative, i.e. divergence-free) contribution ĵ0 eG and a particular solution ĵ0 e P whose divergence equals the right-hand-side photo-ionization term of (67), i.e. ĵ0 e = ĵ0 eG + ĵ0

e P , and

∇ • ĵ0 eG =0, (68)
∇ • ĵ0 e P =γ G(r)M 0 0 . ( 69
)
Since we consider no-incoming electron from any other source, the general conservative contribution, being unique, is zero, ĵ0 eG = 0. Hence, we are left with finding the particular solution ĵ0

e P . From the axi-symmetry of both the source term and the boundary , we can assume that ĵ0 e P = j 0 e (r)e r and thus develop the divergence operator in cylindrical coordinates, only keeping the radial part. Integrating between r (dimensionless radius of boundary ) and infinity leads to

j 0 e (r ) = M 0 0 γ 1 r ∞ r G(r)rdr = M 0 0 γ 0 (r ),
∇ • ĵ1 e = γ G(r)M 1 0 + ∂ r G(r)M 0 1 • e r . ( 72
)
The same consideration applies, at this order and the general conservative contribution to ĵ1

e is thus zero. The electron flux (72) thus results from two contributions. An axi-symetric one, provided by the M 1 0 term, and a dipolar one resulting from M 0 1 . The first one, is similar to the leading-order in , having an amplitude M 1 0 instead of M 0 0 . Seeking for a particular solution, from the axi-symmetry of the boundary , as well as the radial dependence of the photo-ionization term on the right-hand-side of (72), one finds that,

ˆj1 e (r ) = M 1 0 γ 0 + M 0 1 • e r γ 1 , ( 73
)
γ 1 (r ) = γ r ∞ r ∂ r G(r)rdr. ( 74
)
Again, one can find an explicit expression for γ 1 (r ) using

γ 1 (r ) = γ G( r ) - γ r ∞ r G(r)dr. ( 75
)
Hence, provided flux at interface (70) and (73) we found the electron flux in domain 2 to be

ĵe 2 • n 2 | = ˆj0 e | + ε ˆj1 e | (θ) + O (ε 2 ), ( 76
)
with,

j 0 e | = γ 0 M 0 0 , ( 77
)
j 1 e | (θ) = γ 0 M 1 0 + γ 1 M 0 1 • n| = γ 0 M 1 0 + γ 1 M 0 1 • e x cos θ + M 0 1 • e y sin θ . ( 78
)
Hence, it is interesting to realize that the only relevant quantities associated with photo-ionization for the corona discharge modeling are some integrals of the kernel associated with γ 0 (r ) in (71) and γ 1 (r ) in (75) functions. The prescribed functions γ 0 (r ) and γ 1 (r ) are explicitly computed for the derived 2D kernel in (B.10) and (B.12). Relations (77) and (78) provide the electron flux, at the interface between 1 and 2 , given the electron flux inside 1 .

In the forthcoming section we detail how the missing boundary conditions associated with continuity of electrical potential φ, positive charge flux ˆjp and electron flux ˆje at are taken care of by Lagrange multipliers.

Two-domain variational formulation

The two sub-domains are defined on Fig. 2a. A conformal mesh is used at the interface between the sub-domains, as depicted on Fig. 2b. In this section, all numerical fields will be supposed regular and derivable so as they should pertain to functional space H 1 . This choice is a first simple framework, but additional complexity might lead to discontinuous solutions. For example, if some forced convection from an external flow is super-imposed to the electro-drift convection, it is expected that, at large Péclet number, some sharp variations of the positive charge density field develop an up-stream front. If the Péclet number is very large, this front will become a shock, and the regularity of the solutions might be lost. Nevertheless, at finite Péclet number, a physically reasonable hypothesis, regularity will be preserved. Hence, in a more general context, preserving the solution regularity necessitates the addition of diffusive terms in all density fluxes. For the sake of simplicity, we did not consider diffusion here, but it can easily be added to the formulation, so as to preserve regularity.

Corona domain 1 variational formulation

The corona problem ( 34)-( 36) in 1 is associated with variational formulation involving test functions (u

1 , v 1 , w 1 ) ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - 1 ∇ φ0 1 • ∇u 1 + ∂ e 1 u 1 ∂ n φ0 1 - ∂ e 1 n0 p 1 ∂ n φ0 1 v 1 + 1 n0 p 1 ∇ φ0 1 • ∇ v 1 - 1 α0 n0 e 1 Ê0 1 v 1 ∂ e 1 n0 e 1 ∂ n φ0 1 w 1 - 1 n0 e 1 ∇ φ0 1 • ∇ w 1 - 1 α0 ef n0 e 1 Ê0 1 w 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎝ 0 0 0 ⎞ ⎠ , ( 79
)
with, again, Ê0 1 = | ∇ φ0 1 |, is the modulus of the dimensionless electric field. Using boundary conditions ( 46) and ( 51), whilst leaving the contribution of boundary to Lagrange multipliers to be defined later-on, we specialize variational formulation (79) so as to define the following bi-linear functional F 0 1 , for all test functions, i.e. ∀(u

1 , v 1 , w 1 ) ∈ U 1 ×V 1 ×W 1 F 0 1 ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ φ0 1 n0 p 1 n0 e 1 ⎞ ⎟ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ ⎤ ⎥ ⎦ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - 1 ∇ φ0 1 • ∇u 1 + ∂ e 1 u 1 ∂ n φ0 1 1 n0 p 1 ∇ φ0 1 • ∇ v 1 - 1 α0 n0 e 1 Ê0 1 v 1 ∂ e 1 n0 e 1 ∂ n φ0 1 w 1 - 1 n0 e 1 [∇ φ0 1 • ∇ w 1 + α0 ef Ê0 1 w 1 ] ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 80
) with U 1 ⊂ H 1 ( 1 ), V 1 ⊂ H 1 ( 1 ) and W 1 ⊂ H 1 ( 1 )
being functional spaces in 1 . Note that on the second line of (80), the positive charge flux at ∂ 1 = ∂ e 1 ∪ found in (79) is now omitted since the contribution of ∂ e 1 is zero from boundary condition [START_REF] Zheng | Current-voltage characteristics of dc corona discharges in air between coaxial cylinders[END_REF] and the contribution of is left to Lagrange multiplier λ p . The same consideration applies for the third line, but for the fact that the contribution of ∂ e 1 on the electron flux is not zero, whereas the contribution of will be left to Lagrange multiplier λ e . The linearized form of (80) is found from computing its Fréchet derivative ∀(u

1 , v 1 , w 1 ) ∈ U 1 × V 1 × W 1 δF 0 1 ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ φ0 1 n0 p 1 n0 e 1 ⎞ ⎟ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ , ⎛ ⎜ ⎝ δ φ0 1 δ n0 p 1 δ n0 e 1 ⎞ ⎟ ⎠ ⎤ ⎥ ⎦ = (81) ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - 1 ∇δ φ0 1 • ∇u 1 + ∂ e 1 u 1 ∂ n δ φ0 1 1 (δ n0 p 1 ∇ φ0 1 • ∇ + n0 p 1 ∇δ φ0 1 • ∇)v 1 - 1 δ[ α0 n0 e 1 Ê0 1 ]v 1 ∂ e 1 (δ n0 e 1 ∂ n φ0 1 + n0 e 1 ∂ n δ φ0 1 )w 1 - 1 (δ n0 e 1 ∇ φ0 1 • ∇ + n0 e 1 ∇δ φ0 1 • ∇ + δ[ α0 ef n0 e 1 Ê0 1 ])w 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 82
)
with, this time, test functions in

U 1 ⊂ H 1 0 ( 1 ), V 1 ⊂ H 1 ( 1 )
and W 1 ⊂ H 1 ( 1 ) since homogeneous Dirichlet boundary condition is applied to δ φ0 on ∂ 1 , where we have used notation for chain rule δ[ α0 n0

e 1 Ê0 1 ] = δ α0 n0 e 1 Ê0 1 + α0 δ n0 e 1 Ê0 1 + α0 n0 e 1 δ Ê0 1 , (83) δ[ α0 ef n0 e 1 Ê0 1 ] = δ α0 ef n0 e 1 Ê0 1 + α0 ef δ n0 e 1 Ê0 1 + α0 ef n0 e 1 δ Ê0 1 .
(84)

Here, δ Ê0

1 = Ê0 1 ( φ0 1 + δ φ0 1 ) -Ê0 1 ( φ0 1 ) = | ∇( φ0 1 + δ φ0 1 )| -| ∇ φ0 1 | is found from linearizing Ê0 1 . More precisely denot- ing Ê2 1 = ∇ φ 1 • ∇ φ 1 and δ Ê 1 = Ê 1 ( φ 1 + δ φ 1 ) -Ê 1 ( φ 1 ), since Ê2 1 ( φ 1 + δ φ 1 ) = ∇ φ 1 • ∇ φ 1 + 2∇ φ 1 • ∇ δϕ 1 + O (δϕ 1 ) 2 , and since Ê2 1 ( φ 1 + δ φ 1 ) -Ê2 1 ( φ 1 ) = δ Ê 1 (2 Ê 1 + O (δϕ 1 )), one gets δ Ê0 1 = ∇ φ0 1 • ∇δ φ0 1 Ê0 1 , ( 85
) so that δ α0 ef reads δ α0 ef = ∂ α0 ef ( Ê0 1 ) ∂ E 0 1 δ Ê0 1 = ∂ α0 ∂ E 0 1 - ∂ η0 ∂ E 0 1 δ Ê0 1 , = α0
( Ê0 1 ) 2 - ∂ η0 ∂ E 0 1 ∇ φ0 1 • ∇δ φ0 1 Ê0 1 . ( 86
)
The derivative of η is provided in Appendix C. At order O ( ), the variational formulation of ( 37)-( 39) associated with boundary conditions (47), [START_REF] Tirumala | Comparative study of corona discharge simulation techniques for electrode configurations inducing non-uniform electric fields[END_REF], and (78), again, leaving the contribution of to Langrange multipliers, reads ∀(u

1 , v 1 , w 1 ) ∈ U 1 × V 1 × W 1 F 1 1 ⎡ ⎣ ⎛ ⎝ φ1 1 n1 p 1 n1 e 1 ⎞ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ ⎤ ⎦ = (87) ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - 1 ∇ φ1 1 • ∇u 1 + ∂ e 1 u 1 ∂ n φ1 1 1 (n 1 p 1 ∇ φ0 1 + n0 p 1 ∇ φ1 1 ) • ∇ v 1 - 1 ( α1 n0 e 1 Ê0 1 + α0 n1 e 1 Ê0 1 + α0 n0 e 1 Ê1 1 )v 1 ∂ e 1 (n 1 e 1 ∂ n φ0 1 + n0 e 1 ∂ n φ1 1 )w 1 - 1 (n 1 e 1 ∇ φ0 1 + n0 e 1 ∇ φ1 1 ) • ∇ w 1 - 1 ( α1 ef n0 e 1 Ê0 1 + α0 ef n1 e 1 Ê0 1 + α0 ef n0 e 1 Ê1 1 )w 1 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 88
) with again test functions in U 1 ⊂ H 1 0 ( 1 ), V 1 ⊂ H 1 ( 1 ) and W 1 ⊂ H 1 ( 1 )
. It is interesting to note, that, as expected

F 1 1 ⎡ ⎣ ⎛ ⎝ φ1 1 n1 p 1 n1 e 1 ⎞ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ ⎤ ⎦ = δF 0 1 ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ φ0 1 n0 p 1 n0 e 1 ⎞ ⎟ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ , ⎛ ⎝ φ1 1 n1 p 1 n1 e 1 ⎞ ⎠ ⎤ ⎥ ⎦ , ( 89
)
whilst implicitly using α1 and α1 ef given in [START_REF] Nguyen | Computational study of glow corona discharge in wind: biased conductor[END_REF] on the left-hand-side of (89), rather than δ α0 and δ α0

ef given in (86).

Drift domain 2 variational formulation

Again leaving the contribution of boundary to Lagrange multipliers to be defined later-on, we specialize variational formulation of ( 40)-( 41) with homogeneous Dirichlet boundary condition [START_REF] Naidis | Conditions for inception of positive corona discharges in air[END_REF] at ∂ c 2 to define the following bi-linear

functional ∀(u 2 , v 2 ) ∈ U 2 × V 2 F 0 2 φ0 2 n0 p 2 , u 2 v 2 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 2 -∇ φ0 2 • ∇u 2 + 2 n0 p 2 u 2 + ∂ c 2 ∂ n φ0 2 u 2 - ∂ c 2 n0 p 2 ∂ n φ0 2 v 2 + 2 n0 p 2 ∇ φ0 2 • ∇ v 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , ( 90
)
where ∂ n denote the outward derivative at 2 edges. Test functions u 2 and v 2 are chosen beyond functional space U 2 and

V 2 with U 2 ⊂ H 1 0 ( 2 ) and V 2 ⊂ H 1 ( 2 ). The Fréchet derivative of (90) reads ∀(u, v) ∈ U 2 × V 2 δF 0 2 φ0 2 n0 p 2 , u 2 v 2 , δ φ0 2 δ n0 p 2 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ - 2 ∇(δ φ0 2 ) • ∇u 2 + 2 δ n0 p 2 u 2 + ∂ c 2 ∂ n δ φ0 2 u 2 - ∂ c 2 δ n0 p 2 ∂ n φ0 2 v 2 + n0 p 2 ∂ n δ φ0 2 v 2 + 2 δ n0 p 2 ∇ φ0 2 • ∇ v 2 + n0 p 2 ∇δ φ0 2 • ∇ v 2 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 91
)
At order O ( ), similarly with previous section one gets

F 1 2 φ1 2 n1 p 2 , u 2 v 2 = δF 0 2 φ0 2 n0 p 2 , u 2 v 2 , φ1 2 n1 p 2 .
(92)

Coupled formulation

Finally building together variational forms (80) and (90), with the leading order O ( 0 ) boundary conditions at (48)-( 49)

and (77), we arrive at the following monolithic non-linear formulation

F 0 1 ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ φ0 1 n0 p 1 n0 e 1 ⎞ ⎟ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ ⎤ ⎥ ⎦ + ⎛ ⎜ ⎝ -λ ϕ u 1 0 -λ e w 1 ⎞ ⎟ ⎠ = 0, ∀(u 1 , v 1 , w 1 ) ∈ U 1 × V 1 × W 1 , ( 93
)
F 0 2 φ0 2 n0 p 2 , u 2 v 2 + λ ϕ u 2 λ p v 2 = 0, ∀(u 2 , v 2 ) ∈ U 2 × V 2 , ( 94
)
μ ϕ ( φ0 2 -φ0 1 ) = 0, ∀μ φ ∈ R, (95)
μ p (λ p -n0

p 1 ∂ n φ0 1 ) = 0, ∀μ p ∈ R, ( 96
)
μ e (λ e -ˆj0 e | ) = 0, ∀μ e , ∈ R, ( 97
)
thanks to Lagrange multipliers associated with potential (λ ϕ ), positive charge flux continuity (λ p ), and imposed electron flux ˆj0

e | (77) (λ e ) at , whose test functions are respectively denoted μ ϕ , μ p and μ e . The functional space associated with Lagrange multiplier adjoin test-functions on is R ⊂ H -1/2 ( ) [START_REF] Magoulès | Lagrangian formulation of domain decomposition methods: a unified theory[END_REF].

It is interesting to note that the Lagrange multiplier associated to the charged particles only comes one-side in 2 in (96) and not in 1 since the last term in the second line of (93) is zero. This one-side coupling results from the purely hyperbolic nature of the positive charge problems, both in 1 and 2 . At order O ( ) the coupled linear system reads

F 1 1 ⎡ ⎣ ⎛ ⎝ φ1 1 n1 p 1 n1 e 1 ⎞ ⎠ , ⎛ ⎝ u 1 v 1 w 1 ⎞ ⎠ ⎤ ⎦ + ⎛ ⎜ ⎝ -λ ϕ u 1 0 -λ e w 1 ⎞ ⎟ ⎠ = 0 ∀(u 1 , v 1 , w 1 ) ∈ U 1 × V 1 × W 1 , ( 98
)
F 1 2 φ1 2 n1 p 2 , u 2 v 2 + λ ϕ u 2 λ p v 2 = 0 ∀(u 2 , v 2 ) ∈ U 2 × V 2 , (99) μ ϕ ( φ1 2 -φ1 1 ) = 0 ∀μ φ ∈ R, ( 100
)
μ p (λ p -n1

p 1 ∂ n φ0 1 -n0 p 1 ∂ n φ1 1 ) = 0 ∀μ p ∈ R, ( 101
)
μ e (λ e -ˆj1 e | ) = 0 ∀μ e ∈ R, (102)
with imposed electron flux ˆj1 e | (78).

Numerical method

Now considering the numerical discretization of the problem and the resulting meshing issues, we chose to re-express the Corona problem in 1 with external variable r rather than internal one R. This issue renders the numerical implementation much easier, since, upon the convention of common spatial variable scaling of domains 1 and 2 , their interface lies at the mesh conformal interface. Choosing a distinct scaling for the inner region discretization 1 and the outer region 2 would have resulted into two distinct interfaces of two distinct meshes, to map one-another, obviously a less convenient choice of discretization.

Newton solution for leading order

Non-linear problem (93)-( 97) is solved using a Newton method. In this monolithic formulation, the unknown fields

( φ0 1 , n0 p 1 , n0 e 1 , φ0 2 , n0 p 2 , λ 0 ϕ , λ 0 p , λ 0 e ),
δF 0 1 ⎡ ⎢ ⎣ ⎛ ⎜ ⎝ φ0k 1 n0k p 1 n0k e 1 ⎞ ⎟ ⎠ , ⎛ ⎝ u v w ⎞ ⎠ , ⎛ ⎝ δ φ0 δ n0 p δ n0 e ⎞ ⎠ ⎤ ⎥ ⎦ ≡ δF 0 1 ⎡ ⎣ δ φ0 δ n0 p δ n0 e ⎤ ⎦ , ( 103
)
in (81) and

δF 0 2 φ0k 2 n0k p 2 , u v , δ φ0 δ n0 p ≡ δF 0 2 δ φ0 δ n0 p , ( 104
)
in (91). Linearizing (93)-(97) leads to the formal Jacobian operator for the Newton iteration (

δF 0 1 ⎡ ⎣ δ φ 1 δ np 1 δ ne 1 ⎤ ⎦ + ⎛ ⎜ ⎝ -δλ 0 ϕ u 1 0 -δλ 0 e w 1 ⎞ ⎟ ⎠ = -F 1 ⎡ ⎢ ⎣ φ0k 1 n0k p 1 n0k e 1 ⎤ ⎥ ⎦ - ⎛ ⎜ ⎜ ⎝ -λ 0k ϕ u 1 0 -λ 0k e w 1 ⎞ ⎟ ⎟ ⎠ , (105) δF 0 2 δ φ 2 δ np 2 + δλ 0 ϕ u 2 δλ 0 p v 2 = -F 2 φ0k 2 n0k p 2 - λ 0k ϕ u 2 λ 0k p v 2 , ( 106
) μ ϕ (δ φ 2 -δ φ 1 ) = -μ ϕ ( φ0k 2 -φ0k 1 ), ( 107
) μ p (δλ 0 p -δ np 1 ∂ n φ0k 1 -n0k p 1 δ φ 1 ) = -μ p (λ 0k p -n0k p 1 ∂ n φ0k 1 ), ( 108
) 109
Where, from ( 62) and ( 77)

, δ ˆj0 e | in (109) is δ ˆj0 e | = γ 0 1 δ αef n0k e 1 Ê0k 1 + α0k ef δ ne 1 Ê0k 1 + α0k ef n0k e 1 δ Ê 1 .
(110)

Which from using (85) and (86) reduces to

δ ˆj0 e | = γ 0 1 α0k ef δ ne 1 Ê0k 1 + (n 0k e 1 + 1 ( Ê0k 1 ) 2 ) ∇ φ0k 1 Ê0k 1 • ∇δ φ 1 , (111) δ ˆj0 e | = δ J 0 en [δ ne 1 ] + δ J 0 eϕ [δ φ 1 ], (112)
defining the two linear (functional) forms

δ ˆJ 0 en [δ ne 1 ] = γ 0 1 α0k ef Ê0k 1 δ ne 1 , ( 113
) δ ˆJ 0 eϕ [δ φ 1 ] = γ 0 1 (n 0k e 1 + 1 ( Ê0k 1 ) 2 ) ∇ φ0k 1 Ê0k 1 • ∇δ φ 1 . ( 114
)
Implementation details of this leading order O ( 0 ) are provided in D.1. At order O ( ), the Jacobian of the Newton is very similar to the leading order, and again details are provided in D.2. It is shown that, as expected, the O ( ) correction is zero if dipolar corrections are zero; i.e., if M 0 1 = 0. Finally, the inversion of the Jacobian matrix is performed using a sparse LU factorization using the UMFPACK library [START_REF] Davis | Algorithm 832: Umfpack v4.3-an unsymmetric-pattern multifrontal method[END_REF].

Validation

Mesh convergence study

A conformal mesh is generated in the sub-domain 1 while the unstructured mesh generator of FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF], is used for the sub-domain 2 . A mesh convergence study has been carried out. The numerical value of each variable at each DOF is stored in the solution vector u and the quadratic norm of the numerical error, uu ref 2 is estimated by comparing the solution u to a reference solution u ref computed with a fine mesh: 34000 triangles, i.e. 221000 DOFs1 in 1 and 69000 DOFs in 2 . In sub-domain 1 , the three variables are φ 1 , n p 1 and n e 1 , so that the error reads

u 1 -u ref 1 2 = N ϕ1 i=1 ( φi -φref ,i ) 2 + N p1 i=1 (n p,i -np,ref ,i ) 2 + N e1 i=1 (n e,i -ne,ref ,i ) 2 N D O F 1 , ( 115
) with N ϕ1 + N p1 + N e1 = N D O F 1 the number of DOFs in 1 .
Similarly in 2 the variables are φ 2 and np 2 and the error then reads

u 2 -u ref 2 2 = N ϕ2 i=1 ( φi -φref ,i ) 2 + N p2 i=1 (n p,i -np,ref ,i ) 2 , N D O F 2 (116) with N ϕ2 + N p2 = N D O F 2 the number of DOFs in 2 .
A first convergence study, Fig. 3a, shows that the error scales as

N -1 D O F 1 in 1 and N -1 D O F 2 in 2 .
A similar study shows that when N (number of mesh nodes along ) increases at fixed N rad1 (number of nodes in the inner radial direction, Cf Fig. 2b), the error decreases rapidly in 1 (∼ N -2 D O F 1 ) while the convergence rate in 2 remains unchanged (because the unstructured mesh in 2 remains identical to the previous case). The last study, Fig. 3b, shows that the error in 1 decreases very rapidly (∼ N -4 D O F 1 ) when the radial resolution N rad1 increases at fixed azimuthal resolution N . However, increasing N rad1 provides a limited improvement since the error saturates for a given N . The saturation threshold depends on N . It is interesting to note that the error in sub-domains 2 depends on N rad1 despite the mesh T h2 is rigorously the same, whatever the value of N rad1 . This is because the boundary condition on is given by the solution in 1 .

The last investigation deals with r , the position of the boundary . The impact on the solution is assessed by the net ion flux at the collecting electrode Îco . In the range of r considered, the current varies by 3%, showing the robustness of the solution regarding this parameter. Fig. 3c shows that the ratio r /a should not be smaller than 3, otherwise Îco is underestimated, because the ionization region is "cropped". On the other hand, when r is too large, the net current slightly increases, because the space charge term is neglected in the potential equation ( 34) in 1 . Indeed, the asymptotic smallness of the space charge term compared to the potential derivatives is preserved only when r → 0. The yellow-star curve shows that this effect disappears when the space charge term is added into the potential equation. the emitter has a known fixed value (determined by Peek's law as any similar corona onset criterion), (ii) the prescribed electron flux is iteratively determined so as to match the experimental current value. This illustrates the above-mentioned parametric dependence of most previously published DC-corona numerical modeling. The obvious limitation from using phenomenological parameter's is they can only apply to a specific configuration at hand, i.e. a specific gas, at given pressure and temperature for a cylindrical emitter. On the contrary our approach can provide generic predictions for any gas, in any thermodynamic conditions, and possibly non-cylindrical electrode. Finally it is interesting to mention that the Kaptzov's model predictions depicted in Fig. 6a, 6b and 6c are not performed using the usual Peek's law at the emitter but the asymptotic prediction provided in [START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF]. This choice gives a much better prediction for the Kaptzov approximation model compared with the two-domain model (for example, Peek's law gives a total current off by more than a factor 3 compared with the two-domain approach for the 0.33 dimensionless shift in the configuration of [START_REF] Tirumala | Comparative study of corona discharge simulation techniques for electrode configurations inducing non-uniform electric fields[END_REF] explored in Fig. 6c).

Other external shapes

We now consider a square shape for the external collector as depicted in Fig. 5b. A systematic current-voltage characteristic curves exploration is performed in Figs. 7a and 7b as an illustration. One can observe in Fig. 7a that, for an applied voltage of 4.2kV, the current level-off by more than a factor two for an emitter offset being one-half of the domain size. The offset effect on the current is highly affected by the applied voltage, as illustrated in Fig. 7a.

Finally we illustrate in Fig. 8 the versatility of our formulation associated with an unstructured finite-element grid which is able to handle complex shapes. In principle Lipschitz-continuous boundary is needed, for the formulation to apply, and provide a solution. From the numerically performed tests, we experience a weak sensitivity of the numerical convergence to cusps and bumps as the one illustrated in Fig. 8a, with rather sparse grid meshing of ∂ 2 . On the contrary, one still needs a

Conclusion

We present an asymptotically-based two-domain approach of corona discharge within a monolithic formulation. Focusing in 2D configurations (translationally invariant along z), we present a Lagrange multiplier approach taking care of the coupling between each domain within a weak-formulation and a finite element method. The proposed formulation allows to take into account both the detailed physics of the corona discharge as well as the associated Townsend discharge coupled with the drift domain. We believe three main benefits can result from the presented approach • It avoids using any dedicated parameter modeling but instead uses intrinsic kinetically based physical parameters only (but for the ionization parameters γ , the experimental measurement of which is known to be difficult).

• The finite element discretization is flexible and robust: in principle, it can easily adapt to various geometries, with an automatic mesh generator [START_REF] Hecht | New development in freefem+[END_REF].

• There is no prescribed surface electric field value and no assumption on the curvature of the emitter surface. This approach can handle non axi-symmetric emitting surface.

As perspectives, most of the presented analysis could be adapted to negative discharge condition, with a distinct treatment of electron flux at the emitter.

Even for positive DC-corona discharges, at high cathode electric field, one might also consider electron's emission resulting from ion-cathode collisions (i.e collector). In this case, it is possible to show that this effect can be handled from considering an additional contribution in the secondary electron flux j e | , as detailed in Appendix F.

With a suitable adaptation of the secondary photo-ionization kernels, and their multi-polar expansion, three dimensional configurations could also be addressed by a similar asymptotically-based two-domain approach. To be more specific, the presented approach could be adapted to point-plane source configurations. For a single point-plane tip being φ-invariant, the problem has to be re-expressed in the rz plane rather than in the rφ plane considered here. In this case, three additional issues have to be faced. First, the photo-ionization effective kernel G will no longer be z-invariant but will depend on both r and z, G(r, z). Second, the electron density in the drift region should be numerically computed from solving [START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF] from defining an additional variational problem. Third, the numerical solution of the electron density in the drift region should be coupled with the non-linear coupled formulation (93)-( 97) with an additional Lagrange multiplier for the electron flux at . These developments encompass the scope of the present paper, but should deserve proper attention for future developments and application of the presented method.

CRediT authorship contribution statement

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement

This work was supported by the CNES (research contract 5100015475) and the French Occitanie region.

Appendix

G(r, r ) ≡ G(ρ ≡ |r -r |) = ẑ ∈R g( ρ 2 + ẑ 2 ) 4π (ρ 2 + ẑ 2 ) (B.5) G(ρ) = 1 4π ln( λ 2 λ 1 ) ẑ ∈R exp -λ 1 ρ 2 +ẑ 2 -exp -λ 2 ρ 2 +ẑ 2 (ρ 2 + ẑ 2 ) 3/2 (B.6)
This integration along ẑ can be performed so as to obtain an explicit 2D Kernel G G(ρ) = ρ 4π ln( λ 2 λ 1 )

λ 3 1 G 3 0 0 3 0 (-1 2 ,-1,-3 2 ) ( λ 1 ρ 2 ) 2 -λ 3 2 G 3 0 0 3 0 (-1 2 ,-1,-3 2 ) ( λ 2 ρ 2 ) 2 (B.7)
where G m n p q a 1 ,...,a p b 1 ,...,b q z is the Meijer G-function [54], and λ 3 j with j = 1, 2 the cube of parameters λ j of 3D kernel (B.2).

Using known relation between Meijer G-function, i.e. with σ 1 = 1, σ 2 = -1. Note we have used dimensionless λi ≡ λ i L but the product λ i r is dimensionless, so that λi r = λ i r .

Furthermore, since the typical physical range of r is μm, since at atmospheric pressure 1/λ 1 = 1.9 mm and 1/λ 2 = 33 μm, only λ 1 r 1, whilst λ 2 r > 1 (see Figs. B.9 and B.10). Now γ 1 (r ) can be find an explicitly using (75) and relation

∞ z G 3 0 0 3 0 (-1 2 ,-1,-3 2 ) z dz = -G 4 0 2 4 (1,1) ( 1 2 ,0,0,-1 2 ) z , (B .11)
The attachment coefficient in pure dry air at N = 2.5e25 m -3 ( P = 1013 hPa and

T = 293 K) is η = β η exp - E η E (E/E η ) D η , (C.2) with β η = B η N = 2391 m -1 and E η = C η N = 7.428 MV.m -1 . Differentiation gives ∂η ∂ E = E η E -D η η E . (C.3)
And the non-dimensional version η = ηL derivative writes

∂ η ∂ Ê = 1 η Ê -D η η Ê (C.4) with η = ϕ a L E η .

Appendix D. Finite element Newton Jacobian discretization

D.1. Leading order Jacobian discretization

We now consider two triangulation T h1 of 1 and T h2 of 2 being conformal at . The left-hand-side operator of (105)-( 109) is associated with the following Jacobian discrete matrix

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ δF 0 1 0 0 -•, u 1 0 0 0 0 0 0 0 0 0 0 0 -•, w 1 0 0 0 δF 0 2 •, u 2 0 0 0 0 0 0 •, v 2 0 -•, μ ϕ 0 0 •, μ ϕ 0 0 0 0 -np 1 ∂ n (•), μ p -(•)∂ n φ 1 , μ p 0 0 0 0 •, μ p 0 -δ Ĵ0 eϕ 0 -δ Ĵ0 en 0 0 0 0 •, μ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ • ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ δ φ 1 δ np 1 δ ne 1 δ φ 2 δ np 2 δλ ϕ δλ p δλ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = B 0k . (D.1)
Using the space P n (T h ) of Lagrange-P n finite elements (n = 2, 3), the matrix δF 0 1 is built from bi-linear form (81) on triangulation T h1 in 1 \ , and, similarly, δF 0 2 built from bi-linear form (91) on triangulation T h2 in 2 \ . Furthermore, (transposed) vectors δ Ĵ0 en and δ Ĵ0 eϕ are build into triangulation T h1 from linear form (113) and (114).

The non-diagonal terms result from the coupling between Lagrange multipliers and unknowns fields at in (93)-( 97) are matrix resulting from variational formulation and built upon linear forms acting on field φ in L, with test function μ with notation A(•), μ defined as

A(•), μ : L → R φ → A(φ)μ, (D.2)
where operator A is either scalar multiplication by a field, or a scalar multiplication combined with normal derivative ∂ n , and L being either H 1 for potential or charge density fields or H -1/2 ( ) for Lagrange multipliers. To be more precise, in (D.1), the linear operator A(φ) is either identity, i.e.

A(φ) = φ, or A(φ) = φ∂ n φ 1 , or A(φ) = np 1 ∂ n (φ).
B 0k = - ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F 0 1 ⎛ ⎜ ⎝ φ0k 1 n0k p 1 n0k e 1 ⎞ ⎟ ⎠ F 0 2 φ0k 2 n0k p 2 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ - ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ -λ k ϕ , u 1 0 λ k e , w 1 λ k ϕ , u 2 λ k p , v 2 ϕ 0k 2 , μ ϕ -ϕ 0k 1 , μ ϕ λ k p , μ p -n 0k p 1 ∂ n ϕ 0k 1 , μ p λ k e , μ e -ˆj0k e | , μ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (D .3)
the matrix F 0 1 is built from bi-linear form (80) on triangulation T h1 in 1 \ , and, similarly, F 0 2 built from bi-linear form (90), whereas, the right-hand-side vector of (D.

D.2. Linear system for dipolar correction

From the electron flux dipolar correction (78) and (63) one defines the following (functional) linear form,

ˆJ 1 en [n 1 e 1 ] = γ 0 1 α0 ef Ê0 1 n1 e 1 , (D.4)
and, from (78), ( 63), ( 64) and [START_REF] Nguyen | Computational study of glow corona discharge in wind: biased conductor[END_REF],

ˆJ 1 eϕ [ φ1 1 ] = γ 0 1 (n 0 e 1 + 1 ( Ê0 1 ) 2 ) ∇ φ0 1 Ê0 1 • ∇ φ1 1 .
(D.5)

The linear system (98)-(102) admits the following matrix formulation

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ F 1 1 0 0 -•, u 1 0 0 0 0 0 0 0 0 0 0 0 -•, w 1 0 0 0 F 1 2 •, u 2 0 0 0 0 0 0 •, v 2 0 -•, μ ϕ 0 0 •, μ ϕ 0 0 0 0 -np 1 ∂ n (•), μ p -(•)∂ n φ 1 , μ p 0 0 0 0 •, μ p 0 -Ĵ1 eϕ 0 -Ĵ1 en 0 0 0 0 •, μ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ • ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ φ1 1 n1 p 1 n1 e 1 φ1 2 n1 p 2 λ ϕ λ p λ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ = B 1 . (D.6)
Where F 1 1 is built from bi-linear form (87) on triangulation T h1 in 1 \ , and, F 1 2 built from bi-linear form (92) on triangulation T h2 in 2 \ . Furthermore, (transposed) vectors Ĵ1 en and Ĵ1 eϕ are built into triangulation T h1 from linear form (D.4) and (D.5). Right-hand-side vector B 1 in (D.6) is given by

B 1 = - ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ 0 0 0 0 0 0 0 -γ 0 γ 1 M 0 1 • n| , μ e ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ , (D.7)
showing that the dipolar corrections are exclusively due to the (dipolar) feeding of the electron flux from drift domain associated with dipolar vector M 0 1 given by the integral (64) over the leading order fields. Hence, as expected, in the case of axi-symmetrical drift domains, since M 0 1 = 0, there is no dipolar corrections.

Appendix E. Variational formulation of Kaptzov assumption

The assumption, known as Kaptzov's assumption, that the surface electric field at the emitter E a remains constant is often used in simplified corona modeling. It turns out to be a good approximation [START_REF] Zheng | Current-voltage characteristics of dc corona discharges in air between coaxial cylinders[END_REF] as long as the electric field value is correctly chosen. A theoretical justification was recently suggested for axisymmetric configurations [START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF], based on Durbin's work [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF]. In a nutshell, it relies on the idea that the corona discharge reaches a self-sustained regime only for one fixed value of the emitter surface electric field E ac : if E a < E ac , the discharge does not generate enough secondary electrons to be self-sustained, while the case E a > E ac leads to an excess of secondary electron. In the absence of any limiting mechanism this should cause an exponential increase of current. But the natural current amplification due to successive Townsend's avalanches is counterbalanced by a decrease of the electric field due the surrounding positive space charge, leading to an "equilibrium" electric field E ac . This value depends on gas properties and the emitter radius, see [START_REF] Monrolin | Revisiting the positive DC corona discharge theory: Beyond Peek's and Townsend's law[END_REF][START_REF] Naidis | Conditions for inception of positive corona discharges in air[END_REF][START_REF] Zheng | Self-sustained criterion with photoionization for positive dc corona plasmas between coaxial cylinders[END_REF] for example, and probably on its geometry even if this last point hast not been widely investigated yet. The competition of the two aforementioned processes can also lead to oscillations [START_REF] Morrow | Streamer propagation in air[END_REF].

In most simplified numerical simulations, the ionization layer shrinks and is replaced by an effective boundary condition: at fixed emitter potential φ a (resp. electric field E a ), the injected charge density at the corona discharge frontier is adjusted until the surface electric field E a (resp. emitter potential φ a ) reaches the desired value (in fact in this approach the emitter frontier is approximated by the corona surface). This approach necessitates to run the corona simulation several times. The following approach is fully coupled: the constitutive equations and the constraint on space charge are solved simultaneously in a few Newton-Raphson steps. We solve the outer equations in the domain = 2 ∀μ ϕ ∈ R.

λ ϕ v + ∂ c np ∂ n φ v -np ∇ φ • ∇ v = 0 ∀v ∈ V (E.
(E.9)

This approach is very close to the one described by Feng [START_REF] Feng | Application of Galerkin finite-element method with Newton iterations in computing steady-state solutions of unipolar charge currents in corona devices[END_REF] as well as more recently in [START_REF] Coseru | Numerical study of ElectroAeroDynamic force and current resulting from ionic wind in emitter/collector systems[END_REF] except for an additional proposed regularization. There is however a difference: the boundary condition for the ion transport equation is specified as a flux instead of a prescribed Dirichlet value.

Appendix F. Cathode electron emission

In this section we consider the possible additional contribution of Cathode electron emission in our approach. Let us compute the electron flux in the drift region. We start with the electron flux conservation equation at leading order ∇ • j 0 e = S(r).

Integrating over the drift region volume, whilst denoting J the total contribution of boundaries to flux leads to

Fig. 1 .

 1 Fig. 1. (a) Schematic representation of the positive corona discharge problem: collector size L is much larger than emitter diameter a (in black). The origin of the position vector r in (xy) plane is the emitter center. (b) The two-domain approach of section 3: = 1 ∪ 2 , boundary is the interface between

  where ˆje = | ĵe | = ne Ê, ˆjp = | ĵp | = np Ê, and ˆjn = | ĵn | = nn Ê.

  details concerning G and the explicit computation of γ 0 . Now considering the order O ( ), inserting (66) in (65), and[START_REF] Moreau | Ionic wind produced by positive and negative corona discharges in air[END_REF] in (59) one finds

) μ e δλ 0 e

 0 δ ˆj0 e | =μ e (λ 0k e -ˆj0k e | ).

Fig. 6 .

 6 Fig.6. Numerical results for current ratio compared to previous modeling and experiments[START_REF] Tirumala | Comparative study of corona discharge simulation techniques for electrode configurations inducing non-uniform electric fields[END_REF] for a non-axisymetric electrode arrangement, for an emitter of 50 μm (a) and 100 μm (b). (c) Non-dimensional positive ion density for a=50 μ, L=3 mm, at V = 3800 V: half-left obtained using Kaptzov assumption (E a = 1.44 × 10 7 V/m), on the half-right, obtained with the two-scale/two-domain ionization method.

Nicolas Monrolin :

 Monrolin Data curation, Methodology, Software, Validation, Visualization, Writing -review & editing. Franck Plouraboué: Conceptualization, Formal analysis, Funding acquisition, Methodology, Project administration, Validation, Writing -original draft, Writing -review & editing.

Fig. B. 9 .R = ρ 2 + ẑ 2 (B. 3 )

 923 Fig. B.9. γ 0 r /γ versus dimensional r plotted between [10 μm, 400 μm] using 1/λ 1 = 1.9 mm, 1/λ 2 = 33 μm.

  (B.8) with variable change z = (λ i ρ/2) 2 with i = 1, 2, one gets γ 0 (r )

  Finally it is important to stress that each Lagrange multipliers λ, as well as their adjoint test-functions μ (not to be confused with mobility) are discretized upon Lagrange-P 1 finite elements over the conformal sets of points of T h1 ∩ = T h2 ∩ . Finally the right-hand-side of (105)-(109) described by vector B 0k in (D.1) is given by

  3) uses notation (D.2). In (D.3), the electron flux at the edge of the corona discharge is ˆj0k e | = γ 0 1 α0k ef ˆj0k e from (62) and (77).

  φ = -n p ∂ n φ| = E ac and φ| ∂ c = 0 ( E . 1 ) ∇ • (n p Ê) = 0 np Ê • n| = λ ϕ (E.2) φa | -1 = 0, (E.3) with an additional constraint on potential enforced by equation (E.3). At the first glance the problem defined by the elliptic equation (E.1) and constraint (E.3) seems ill-posed because at boundary , the potential has a Neumann condition ∂ n φ| = E ac and a Dirichlet one φa = 1. This over-determination is balanced by the need for an additional upstream condition for the ion flux of the hyperbolic electro-convection problem. As in [27], this additional flux condition is estimated thanks to the introduction of a new unknown, the Lagrangian multiplier λ ϕ , which is the flux of positive ions entering the domain at the emitter, estimated so that the constraint on the potential at the emitter is fulfilled. The variational formulation of this set of equation writes -∇ φ • ∇u + np u + E ac u = 0 ∀u ∈ U (E.4)

5 )

 5 μ ϕ ( φ -1) = 0 ∀μ ϕ ∈ R, (E.6)and can be solved iteratively. The linearized system at each step writes:-∇δ φ • ∇u + δ np u = -⎛ ⎝ -∇ φ • ∇u + np u + E ac u ∂ n φ + np ∂ n δ φ)v -(δ np ∇ φ + np ∇δ φ) • ∇ v = -⎛ ⎜ ⎝ λ ϕ v + ∂ c np ∂ n φ v -np ∇ φ • ∇ v φ =μ ϕ ( φ -1)

J 0 e 2 S 2 S(r) = γ 2 G 1 (α -η) j e , which leads to J 0 e | = γ 2 G 1 ( 0 e | ∂ c = γ c 1 2 G(r)rdr 1 (α -η) j 0 e + γ c 1 α j 0 e . 2 ) 2 G

 0222102101211e22 | -J e | ∂ c = (r).In this paper, we have considered a zero electron-flux contribution at the cathode J e | ∂ c = 0. Furthermore we have shown that the net photo-ionization rate, at leading order, boils down to (r, 0)rdr• (r, 0)rdr • αη) j 0 e .This first initial footstep shows that the secondary electron flux is proportional to the integral of the ionization rate inside the corona region. If, on the contrary, we take into account the contribution of additional ion-cathode secondary electrons,then J e | ∂ c = γ c J p | ∂ cwith γ c a secondary electron emission coefficient for ion-cathode collisions. To evaluate the net ion flux we integrate the ion conservation equation ∇ • j 0 p = α j 0 e , over the corona volume 1 , delimited by the boundary , to obtain:The total flux of positive ions at the edge of the corona is (at the leading order) the same as the one hitting the collector surface because of the flux conservation in the drift region J Since the ionization rate α is orders of magnitude greater than the attachment rate η the integral can be approximated by 1 (αη)This approximation shows that ion-cathode collision can be taken into account from choosing an adequate value of the secondary ionization coefficient γ ef f = γ (r)rdr + γ c .

  with embedded Lagrange multipliers λ's, are simultaneously computed. At each k + 1

	step, their increments
	(δ φ0	1 , δ n0 p 1 , δ n0 e 1 , δ φ0	2 , δ n0 p 2 , δλ 0 ϕ , δλ 0 p , δλ 0 e ) are found, from the inversion of the Jacobian operator (105)-(109) knowing
	step k fields ( φ0	1 , n0 p 1 , n0 e 1 , φ0	2 , n0 p 2 , λ 0 ϕ , λ 0 p , λ 0

e ) k . For notation simplification we will now denote

A. Nomenclature

  , v k , w k Test functions for variational formulation in k , k = 1, 2 U k , V k , W k Functional spaces of test functions in k , k = 1, 2μ ϕ , μ p , μ e Test functions for Lagrangian mutlipliers λ ϕ , λ p , λ e

	g(R, R ) 3D photo-ionization kernel
	G(r, r ) 2D photo-ionization kernel
	j e	Electron charge flux
	j p	Positive ion charge flux
	j n	Negative ion charge flux
	L	Typical distance between Emitter and collector
	λ	Typical distance decay of photo-ionization kernel
	λ e	Lagrange multiplier for electron flux at
	λ p	Lagrange multiplier for positive charge flux at
	λ ϕ	Lagrange multiplier for electric potential at
	μ e	Electron mobility
	μ p	Positive ion mobility
	μ n	Negative ion mobility
	M k 0	monopolar moment of the photo-ionization kernel k = 0, 1
	M 1	Dipolar moment (2D vector) of the photo-ionization kernel
	n	Outward normal to any domain edge
	N	Neutral gas number density
	n e	Electron number density
	n p	Monopolar positive ion number density
	n n n k	Monopolar negative ion number density ≡ e 0 ϕ a L 2 μ p reference number density, k = e, p, n μ k ≡ 1 ∪ 2 2D domain between emitter and collector
	1	Corona discharge 2D domain
	2	Drift region 2D domain
	S	Photo-ionization source term
	∂ e 2 ∂ c 2 ϕ	Edge of emitter Edge of collector Electric potential
	ϕ a R ≡ r L R	Applied electric potential (at the emitter) Rescaled coordinate in 1 ≡ r + ze z 3D coordinate
	r	2D coordinate in
	H 1	Sobolev functional space
	u k R	Functional space of Lagrangian multipliers in
	T hj	Triangular mesh for finite element P j , j = 1, 2
	a α	Emitter diameter ≡ βe -E i E Townsend discharge ionization coefficient
	α ef f δ μ	≡ α -η Effective ionization coefficient ≡ μ p μ e Small parameter: positive ion to electron mobility ratio
	η	Attachment coefficient
	0 ≡ ϕ a L E i	Electric constant/vacuum permittivity Small parameter comparing the applied electric field to the ionization field
	e	Electron elementary charge
	E	Electric field
	E i	Ionization electric field
	γ	Photo-ionization coefficient
		Frontier between corona domain 1 and drift domain 2

DOFs: Degrees Of Freedom.

Appendix B. Derivation of the 2D photo-ionization kernel and its relevant integrals γ 0 and γ [START_REF] Durbin | Analysis of the positive DC corona between coaxial cylinders[END_REF] As discussed in [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF][START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF] 3D photo-ionization kernels are given by

where the relative 3D distance between points R and R is denoted R = R -R , and where the radiative intensity I( R ) is here taken proportional to the electron flux in [START_REF] Ferreira | Simulation of pre-breakdown discharges in high-pressure air I: the model and its application to corona inception[END_REF], i.e. I( R ) = γ α ef f ( R ) j e ( R ). The Kernel g is given by [START_REF] Bourdon | Efficient models for photoionization produced by non-thermal gas discharges in air based on radiative transfer and the Helmholtz equations[END_REF][START_REF] Zheleznyak | Photoionization of nitrogen and oxygen mixtures by radiation from a gas discharge[END_REF]

where λ 1 = χ min P O 2 , λ 2 = χ max P O 2 , and χ min = 0.035 Torr -1 cm -1 , χ max = 2 Torr -1 cm -1 and where P O 2 is the partial pressure of molecular oxygen. At atmospheric pressure, P O 2 = 150 Torr, so that 1/λ 1 = 1.9 mm and 1/λ 2 = 33 μm. In this appendix we derive an effective 2D kernel for a cylindrical geometry based upon the integration of 3D kernel (B.2) in (B.1).

Doing so, we decompose each 3D vector R into a cylindrical radial term and a longitudinal one along e z , i.e: and R = r + ẑ , with, obviously, ẑ = ẑ e z . Now, for any given plane transverse to the cylinder, one can take this plane as the origin for the cylindrical coordinate system, so that R = r. Furthermore, one has to recognize that, because the source term is translationaly invariant along the z direction, then I( R ) = I(r ). Defining notation ρ = r-r for the relative 2D (cylindrical) radial position vector, having amplitude ρ = |ρ|, one gets