{'fraction_non_alphanumeric': 0.046493259374943904, 'fraction_numerical': 0.020159046439406178, 'mean_word_length': 4.279878684484883, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Moving boundary approximation for curved streamer ionization fronts: Solvability analysis

The minimal density model for negative streamer ionization fronts is investigated. An earlier moving boundary approximation for this model consisted of a"kinetic undercooling"type boundary condition in a Laplacian growth problem of Hele-Shaw type. Here we derive a curvature correction to the moving boundary approximation that resembles surface tension. The calculation is based on solvability analysis with unconventional features, namely, there are three relevant zero modes of the adjoint operator, one of them diverging; furthermore, the inner/outer matching ahead of the front has to be performed on a line rather than on an extended region; and the whole calculation can be performed analytically. The analysis reveals a relation between the fields ahead and behind a slowly evolving curved front, the curvature and the generated conductivity. This relation forces us to give up the ideal conductivity approximation, and we suggest to replace it by a constant conductivity approximation. This implies that the electric potential in the streamer interior is no longer constant but solves a Laplace equation; this leads to a Muskat-type problem.


I. INTRODUCTION
Streamers are growing plasma channels that in sufficiently strong electrostatic field expand into ionizable matter like air or other gases, liquids or solids. The field accelerates the electrons sufficiently to ionize gas molecules upon collisions. The self-focussing of the field at the tip of the propagating streamer strongly supports this process. The richness of spatio-temporal structures that form as streamers undergo successive branching events during their propagation make them an important example of nonequilibrium pattern formation. They are of fundamental interest not only for this reason, but also because they determine the initial stages of electric breakdown in sparks and lightning, in technological applications as well in natural processes [1,2,3,4]. In this paper we focus on negative (anode-directed) streamers in simple gases like pure nitrogen or argon. Their dynamics can be described by a minimal model of nonlinear reaction-advection-diffusion equations, as described, e.g., in [5,6]. The model describes the evolution of electron and ion densities and their nonlinear coupling to the local electric field.
Many simulations [7,8,9,10,11,12,13,14] since the early 1980'ies have shown that the ionized interior of the streamer finger is preceded and surrounded by a narrow curved space charge layer that enhances the electric field in the non-ionized region ahead of the front and screens it in the ionized interior. In particular, in strong background fields after a sufficiently long evolution, the width of the ionization front can be much smaller than its radius of curvature [13,14,15]. This separation of scales enables one to consider the front as an infinitesimally thin, sharp moving interface Γ. Many inhomogeneous systems involve domains of well defined phases separated by thin interfaces. These include non-equilibrium systems like solidification occurring by dendritic growth  (30)) of the ionization front. On the right: the moving boundary approximation (MBA) with the ionized region, Ω − , the non-ionized region, Ω + , and the sharp interface, Γ. [16,17,18,19,20,21,22], bacterial growth [23,24] or many diffusive-reactive systems [25]. In Fig. 1, we present a generic solution for the net charge density of the minimal PDE model showing the separation of scale together with the resulting moving boundary approximation. The original nonlinear dynamics is then replaced by a set of linear field equations (typically Laplace) on both sides of Γ, Ω + and Ω − , with appropriate boundary conditions at the interface, Γ + and Γ − respectively, and further away from it [21,26,27,28] (and references therein); the nonlinearity enters through the motion of the boundary. The interface dynamics is then typically related to gradients of the Laplacian fields at its vicinity, Γ + and Γ − .
In the context of streamer dynamics, the concept of an interfacial approximation was proposed already in the early 1970'ies by Lozansky and Firsov before any numerical simulations [29]. At the time, their model was in competition with other streamer models and concepts. Lozansky and Firsov proposed to consider the streamer interior, Ω − , as ideally conducting, i.e., the electric potential φ to be constant in the interior. The exterior, Ω + , is non-ionized and therefore does not contain space charges, solving in electrostatic approximation (that is justified in [6]) The interface was assumed to move with the local electron drift velocity v = ∇φ + .
From here, superscripts ± attached to fields, potential and densities indicate their limit value as they approach the interface from Ω + and Ω − , respectively. In particular, we denote φ + ≡ φ| Γ + and φ − ≡ φ| Γ − . However, the solutions of this interfacial model were hardly explored. Only 30 years later, the interfacial model was taken up again in the streamer context [30], and Meulenbroek et al. showed that the model actually allowed for spontaneous streamer branching. However, the model was also found to be mathematically ill posed; in the context of fluid dynamics, this is explain for example in Ref. [31] and references therein. To resolve this problem, the boundary condition φ + − φ − = 0 was replaced by the regularizing boundary condition which was proposed in [32] and derived in planar front approximation in [33]. The boundary condition accounts for the finite width of the charged layer (the ionization front) that leads to a finite variation of the electric potential across the front. The boundary condition in the limit of large electric fields actually turns out to be identical to the "kinetic undercooling" boundary condition that was applied to crystal growth under certain conditions [34,35]. The solutions of the interfacial model with this boundary condition (3) implemented were studied in [32,33,36,37], and the results of these papers show that the interfacial model employed there is the simplest one that regularizes the motion -therefore we call this model the minimal regularized model. However, it also was shown recently that this analysis applies only to boundaries that at time t = 0 are sufficiently many times differentiable [38].
In this paper we compute curvature corrections to the boundary condition (3). A systematic expansion of slightly curved fronts about planar fronts for streamers was first suggested in [5], and [6] contains the analysis of planar fronts as first step in this research program. In analogy with other weakly curved fronts, a perturbative expansion of the curved front about the planar front with subsequent solvability analysis was the method at hand. Solvability analysis means that the inner front structure is integrated out and replaced by conditions for interfacial motion that are matched to the dynamics on the outer scale.
However, such an analysis could not be performed on the streamer model of [5,6] as the fronts are pulled; for pulled fronts, we refer to [39,40], and for the impossibility of a solvability analysis for pulled fronts to [41]. The formal reason for the non-applicability of the method are diverging integrals in the leading edge of the ionization front; the physical reason is the algebraically slow front relaxation of the leading edge that actually stretches out through the whole non-ionized region and does not make part of the inner front region that is to be integrated out. The pulled nature of the front also requires special care with grid refinement in numerical solutions [14].
However, the leading edge that pulls the front, is diffusive, and it is a physically and mathematically meaningful approximation to neglect electron diffusion in strong fields, where electron motion due to drift dominates over the diffusive motion [42,43]. We have employed this approximation in the interfacial model and we have checked that its results approximate simulations with electron diffusion well [33,44]. We here employ the same approximation which now makes it possible to perform a solvability analysis.
In fact, neglecting electron diffusion and assuming that the state ahead of the front is completely non-ionized, converts the inner-outer matching problem from one nontypical situation to the other. With electron diffusion, the inner front region was not separated from the outer one. Without electron diffusion, however, the inner region finishes precisely where the electron density discontinuously jumps to 0, and the matching between inner front region and the non-ionized outer region takes place precisely on this line rather than on an extended spatial region.
The "moving boundary condition" (MBC) for a slightly curved front dynamics now can be systematically derived from the original nonlinear field equations: A perturbation of a planar front is assumed whose curvature in the direction transverse to the front motion is much smaller than the front width, and solvability analysis is used to connect the perturbed values of the fields ahead and behind the curved front. Recently, such approach has been successfully applied to derive MBC equations for the dynamics of discontinuity curves which appear in nonlinear diffusion problems [24].
The solvability analysis poses technical challenges. We will show that the electrostatic field behind a slowly evolving curved front does not vanish but is rather constant and proportional to the front curvature. Consequently, the electric potential diverges linearly behind the front. Another divergence appears in one of the zero modes of the solvability analysis. This divergence is how-ever necessary to cancel the one of the electric potential behind the front in order to get a meaningful boundary condition for the jump in the electric potential across the interface Γ. To our knowledge, this nontrivial computational aspect is new, and we are not aware of similar types of solvability analysis in implementations of sharp interface approach to front dynamics in other physical systems.
Furthermore, since the potential is no longer constant behind the front, the ideal conductivity approximation (φ = 0 in Ω − ) must be relaxed and is replaced by a constant conductivity approximation for the ionized streamer interior (∇ 2 φ = 0 in Ω − ).
Explicitly, our front dynamics consists of the following equations: where the coefficients Q i depend on the electrostatic field ahead the front, and are given by analytic formulas derived from the planar front solution, see Eqs. (129)-(131). The paper is organized as follows: In Sec. II we summarize previous results on the minimal model and its planar front solutions essential for our analysis. In Sec. III A, we start our formal derivation by assuming a slightly curved front with ǫ = κℓ − α ≪ 1, where ℓ − α and κ are the width and the curvature of the front respectively. We construct a perturbative expansion of the solution to the minimal model dynamics in powers of ǫ around the zeroth order planar front solution. We focus on the leading order behavior in ǫ and describe it as a solution to a set of four inhomogeneous linear ODEs. In Sec. III B, we explain and develop the solvability analysis formalism: Using zero modes of the adjoint linear differential operator we derive relations between deviations of the electrostatic field and ion/electron densities ahead and behind the curved front, Γ + and Γ − , from their values at the planar front solution. These relations turn out to be the required MBC equations for the curved fronts dynamics. The required zero modes are calculated analytically in Sec. III C. In Sec. III D we conclude our inner analysis and derive the MBC equations. The final form of the MBC equations (124)-(128), obtained by matching inner and outer regions, is derived in Sec. IV where we give also the large field limit of the interfacial model together with an elementary algorithm to solve it. In Sec. V, we discuss the constant conductivity approximation for the ionized region. In Sec. VI we conclude by highlighting the major results of this paper and point to future directions.

II. COLLECTION OF SOME PREVIOUS RESULTS
In this section we summarize previous results required for the present analysis. We briefly describe the minimal model, we discuss the front velocity and the coupling between the front and the non-ionized exterior region, and we cite several properties of planar uniformly translating fronts required for our further analysis.

A. The minimal streamer model
The model for negative streamers in simple nonattaching gases like nitrogen and argon as used, e.g., in [5,9,10,11,12,13,14,33] consists of a set of three coupled partial differential equations for the electron density σ, the ion density ρ and the electric field E. In dimensionless units the model reads A general discussion about dimensions for this model can be found, e.g., in [4,5,6,15]. The first two equations are the continuity equations for the electrons and the ions while the last is the Coulomb equation for the electric field generated by the space charge ρ − σ of electrons and ions. Here the electrostatic approximation for the electric potential φ was introduced [6]. Notice that the electron particle current is here taken as the drift current only −σE, neglecting electron diffusion. For the effect of a diffusive contribution, see the discussion in the introduction as well as [43]. The ion particle current is neglected because due to their larger mass the ions are much less mobile. The term σ|E|α(|E|) is the generation rate of additional electron ion-pairs; it is the product of the absolute value of the drift current times the effective cross section α taken as field dependent; a commonly used form is the Townsend approximation However, our analysis holds for the more general case of

B. Ionization fronts and their velocity
We consider planar ionization fronts propagating in the positive x-direction into a medium that is completely non-ionized beyond a certain point x f (t) and we work in the comoving coordinate As is easily verified by rewriting Eq. (9) in the comoving coordinate ξ with the help of (11) as wherex is the unit vector in the x-direction, the velocity of the front is determined purely by the value of the field E + on the discontinuity of the electron density Notice that, in order to simplify notations later, E + is a positive quantity. Obviously, this relation between front velocity and local field is not restricted to planar fronts only, but holds generally: As the space charge ρ − σ never diverges, the field E is continuous across r f .

C. Coupling the front to the outer non-ionized region
In the purely one-dimensional planar setting, the field does not vary ahead of the front: ∂ x E = 0 for x > x f (t). On the other hand, the curved ionization front around the tip of a streamer finger creates the characteristic field enhancement ahead of the streamer tip. However, the region x > x f (t) can be left purely to the outer scale analysis of the Laplace equation and has no further influence on the front solution except for determining the field E + at r f (t). The matching between inner and outer region here is completely concentrated on the line r f (t). As also remarked already in the introduction, these two spatial regions only decouple if the electron diffusion D (that is included in numerical streamer models) is neglected. Putting D = 0 eliminates the leading edge and the pulled dynamics of the front [6,42].

D. Uniformly translating planar ionization fronts
If E + is time independent, ∂ t E + = 0, the front propagates uniformly with velocity v = E + . Uniformly translating front solutions of the minimal model (9)-(11) in a field E + depend only on the comoving coordinate ξ = x − vt. They have been discussed in many previous papers [5,6,42], for a recent thorough discussion, we refer, in particular, to section II of [33] and for the comparison of fronts with or without electron diffusion D to section 2 of [43].
Analytically the front solution in the ionized region ξ < 0 is given implicitly by the equations Fig. 2 shows the shape of the solution for E + = 1. At the shock front, the electron density jumps from zero to while ion density and electric field are continuous across the interface For ξ → −∞ the densities approach the value which means that the space charge density ρ−σ vanishes, and the electric field vanishes as well Far behind the front, where E is so small that α(|E|) ≃ 0, the profiles of E and σ decay exponentially: For more details on the asymptotic behavior at ξ = 0 and ξ → −∞, we refer to section II of [42]. The exponential decay in (28) suggests a natural definition of the front width as The length scales within the front are further elaborated in section 6.3.2 of [43] where next to ℓ − α also another scale ℓ + α = 1/α(E + ) was defined that characterizes the inner front region close to ξ = 0 for D = 0. (Nonvanishing electron diffusion D introduces another scale ℓ D .) The length scale ℓ + α is obviously of no relevance for matching of the outer scale, but ℓ − α is the appropriate quantity. From the uniformly translating planar front solution, the boundary condition (3) was derived in [33] as follows. Integration of the electric potential, ∂ ζ φ(ζ) = −E(ζ), from −∞ to 0 leads to Substituting the variable E(ζ) = −x with the help of (22), and using the limit (27) for the lower limit of integration, we obtain where the function ρ(x) is defined in (21). As discussed in detail in [33], the function Q 0 depends nonlinearly on E + (see Fig. 5), but asymptotically Finally, let us also recall an important property of the planar front solution which is a useful relation used later in the paper This relation is obtained by taking the derivative of Eq. (22) and using Eq. (20).

E. The minimal regularized boundary model
The moving boundary model derived and evaluated in [32,33,36,37] is given by Eqs. (1)-(3) together with the ideal conductivity approximation φ = const. in Ω − that will further be discussed in Sec. V. This model can be considered as a minimal regularized boundary model. While the model with the boundary condition φ + − φ − = 0 treated in [30] lacks regularization and for generic initial conditions breaks down within infinitesimal time, the model of [32,33,36,37] according to our analysis does stay regularized as long as the initial contour is infinitely many times differentiable. It consists of the exact description of the non-ionized region (19), of the exact relation between velocity and local field (18) and furthermore of two approximations for the full dynamics of the PDE's (9)- (11). First, the jump of the electric potential across the front is modeled using the approximation (33) of a planar uniformly translating front in the limit of large field E + ; corrections to this behavior will be studied through the solvability analysis from section III on. Second, the streamer interior is modeled as ideally conducting: φ = const. in Ω − . This approximation will be discussed and improved in Sec. V.

III. BOUNDARY APPROXIMATION FOR SLIGHTLY CURVED FRONTS
To derive curvature corrections to the minimal regularized boundary model (section II E) from the underlying PDE's described in section II A, we follow the general approach as described, e.g., in [21,24,26,27,28,41,45] for solidification fronts in supercooled melts, nonlinear diffusion fronts or interface between two inmiscible fluids, and adapt it to our problem as necessary.

A. Perturbative expansion about planar fronts
In a first step, we expand slightly curved fronts about the uniformly translating planar fronts that were recalled in section II D. Weak curvature means that the width of the front ℓ − α (30) is much smaller than the mean radius of curvature R = 1/κ of the front. The quantity is then a small parameter. The inner region of the curved front can be expanded about the planar front by introducing a curvilinear orthogonal coordinate system, with coordinates (η, ζ) locally tangential and normal to the moving front. We expand the inner region as where [φ 0 (ζ), σ 0 (ζ), ρ 0 (ζ)] and v 0 = E + are the planar uniformly translating front solutions (20)-(22) within a fixed field E + . Assuming that the variations in the transversal direction η are of the order of the radius of curvature R, it is useful to rescale the transversal coordinate asη = ǫη. A divergence ∇ · J then can be expanded in ǫ as and the Laplacian ∇ 2 φ as well: Therefore, the equations for [φ 1 , σ 1 , ρ 1 ] and v 1 depend only on ζ as common.
Substituting the expansion (36) in the minimal model equations (9)-(11) we find to O(ǫ 0 ) the planar front Eqs. (20)- (22). We obtain to O(ǫ) the equations: where we used and where These equations are linear for (σ 1 , ρ 1 , E 1 , φ 1 ) with inhomogeneous terms v 1 ∂ ζ σ 0 , v 1 ∂ ζ ρ 0 and E 0 . The structure of this system of ODEs reads B. The structure of the solvability analysis for streamer fronts In this subsection, using the zero modes of the adjoint operator L * of L, we formally derive corrections to the boundary relations derived from the planar front solution. The explicit expressions of the coefficients in these relations are calculated in later subsections.

Formal procedure of solvability analysis
The boundary conditions are evaluated through a solvability analysis. We first introduce its general formalism.

Consider a system of inhomogeneous linear ODEs
where L is a linear differential operator that can be written as and where M and D are matrices whose entries depend on ζ. The adjoint operator of L is given by where t stands for the transpose of the matrix. Assume that a zero mode, U * , of L * is known Multiplying Eq. (50) with (U * ) t , integrating over [a, b] and integrating by part, one easily verifies that Here the arbitrary integration boundaries a and b were introduced to control possible divergences of the integral and the boundary term. The relation (54) is the basis of the solvability analysis performed in this work. As we will show in the next sections, taking the limits a → −∞ and b → 0, relation (54) yields linear relations between the perturbed electrostatic field and potential on both ends of the inner region, thus providing boundary conditions for the Laplace equations (4) and (5) in the outer regions. It should be remarked that the different coupling to the outer regions on both sides of the front is reflected by the asymmetric integration boundaries. At ζ = 0, there is an abrupt transition from the inner front region to the outer non-ionized region as described in section II C; here the matching region shrinks to a line. For ζ → −∞, there is an exponential decay to charge neutrality and the outer ionized region as discussed later in section V.

Zero modes of L * and MBC equations
Denoting the components of the zero modes as relation (54), with a = ζ c and b = 0, leads to the relation where In order to obtain the MBC equations we have to consider the limit ζ c → −∞. This limit has to be analyzed with caution due to the mutual cancellation of diverging terms. In order to extract from Eq. (56) appropriate MBC equations, we use some properties of the perturbed field and densities. First, from Eq. (18) and the expan- Second, asymptotic analysis of Eq. (39) for ζ → −∞ leads to which ensures charge neutrality at ζ = −∞ to order ǫ. Third, we show below that ψ σ (0) does not diverge. Since from Eq. (18), v 0 + E 0 vanishes at ζ = 0, the finite value of σ 1 (0) does not enter the MBC equations. Finally, since ρ is a continuous function and ρ 0 vanishes for ζ > 0 (26), ρ 1 (0) must vanish as well: With these relations, Eq. (56) takes a simpler form where we assumed |ζ c | ≫ ℓ − α , such that E 0 (ζ c ), σ 0 (ζ c ) and ρ 0 (ζ c ) reach their plateau values (28).
In order to extract the actual MBC from Eq. (62), we must find the zero modes of L * . The number of independent MBC equations is related to the number of zero modes that give rise to converging integrals in the expressions for the coefficients A, B (see Eqs. (57) and (58)). The calculation of these zero modes and the consequent evaluation of the required integrals are the subject of this section and the next one. 3. The asymptotic behavior of zero modes at ζ → −∞ The construction of the adjoint operator L * from Eq. (52) is straightforward. The system of equations for the zero modes, Eq. (53), is: Equation (66) immediately gives As we explained above, the zero modes which are relevant for the derivation of MBC equations are those that give rise to converging integrals in Eqs. (57) and (58). It is thus essential to analyze first their asymptotic behavior in the limit ζ → −∞. In this limit, Equation (70) immediately gives: where c E is another constant. Substituting Eq. (71) in (68) and integrating we obtain where c σ is a third constant of integration. Finally, the asymptotic expression for ψ ρ is obtained by substituting Eqs. (71) and (72) into (69) and integrating: where the fourth constant of integration c ρ was introduced.
The set of four independent constants c σ , c ρ , c E , c φ corresponds to four independent solutions of the set of four ODE's (63)-(66).
The exponential divergence of the coefficient of c σ in equation (72) and the exponentially decaying profile of σ 0 (ζ) with the same length scale (28) imply that the integral (57) converges only if c σ = 0. This condition reduces the number of relevant zero modes, and hence the number of MBC equations derived from Eq. (62), to three.

Initial conditions for zero modes at ζ → −∞
In order to extract the three relevant zero modes from Eqs. (63)-(66) in the whole interval ζ ∈ (−∞, 0), it is natural to specify initial conditions α , such that the zero modes are properly described by their asymptotic behavior, Eqs. (71), (72) and (73). Obviously, our final result must not depend on our choice of ζ ini . Using Eqs. (67), (71) and (72), and imposing the condition c σ = 0 we obtain the relation: (74) From Eqs. (71)-(73) and the constraint (74) we obtain is nonvanishing, the zero modes diverge linearly with ζ as ζ → −∞. Since at least one of the three required independent relevant modes must have ψ φ (ζ ini ) = 0, we must include the computation of such a diverging mode in our analysis. For computational purposes, it is convenient to choose the three independent orthogonal sets of initial conditions at ζ = ζ ini so that only one of them contains a nonvanishing value of ψ φ (ζ ini ). It can easily be shown that choosing another set of three independent initial conditions, lead to identical results for the MBC. With our convenient choice we can immediately extract from Eqs. (75)-(77) the initial conditions at ζ = ζ ini (taking Eq. (74) into account) for the two nondiverging modes [ψ and  where the linearity of Eqs. (63)-(66) was used to scale one of the coefficients to unity in each initial condition (78) and (79). The last set of initial conditions, for the linearly diverging zero mode [ψ φ ], is determined by requiring it to be orthogonal to the two vectors (78) and (79) (and to satisfy the constraint (74)). This condition gives, up to a normalization factor: Thus, in order to obtain the three zero modes relevant for the MBC equations, from Eq. (62), we have to solve Eqs. (63)-(66) in the regime ζ ∈ (−∞, 0) under the initial conditions, Eqs. (78)-(80). Computation of the linearly diverging mode [ψ φ ] and extracting nondiverging MBC equations from Eq. (62) corresponding to this mode involves nontrivial technical complications, which for the sake of clarity we defer to the following section. We first discuss the two MBC equations corresponding to the two converging zero modes [ψ

The nondiverging zero modes
In this subsection we derive a general form of two MBC equations. The exact expressions for the coefficients are deferred to section III C. As we show below, the general form of the two MBC equations we derive here is sufficient to deduce an important property: The perturbed electrostatic field E 1 (−∞) does not vanish in the inner region behind the front. This observation emphasizes the necessity to derive the third MBC equation with the aid of the linearly diverging zero mode.
For the two nondiverging modes we can take ζ ini = ζ c → −∞. Substituting Eqs. (78) and (79) in Eq. (62) we obtain the two MBC equations: and (82) Let us reconsider now our perturbative expansion (36). As we emphasized in previous sections, the only physical source of perturbations to the planar front solution in our analysis are the transverse perturbations which correspond to nonvanishing κ = λǫ. Thus, we formally write the perturbation to front velocity as ǫ v 1 , requiring it to vanish for ǫ = 0. The calculations presented so far, however, do not involve any assumption on the actual value of v 1 , except ǫ v 1 ≪ v 0 . By keeping κ = 0 in our former analysis, we can thus consider the formal limit process: ǫ → 0 with ǫ v 1 small (compared to v 0 ), but finite. This simply means that we consider a perturbation of the electric field for a pure planar front which implies which implies This important identity will be confirmed in Sec. III D after finding the exact solution for the converging zero modes, and evaluating the actual values of A (i) + ψ (i) Furthermore, substituting the identity (84) in Eqs. (81) and (82) we obtain the equation: This relation implies that the electric field behind our curved fronts does not vanish if B (1) = B (2) . Finally, let us use the relation (85) to derive a correction to the boundary condition (3). Integrating the equation ∂ ζ φ 1 = −E 1 , we get: Since Eq. (85) implies that E 1 (ζ) reaches a constant value E 1 (−∞) as ζ → −∞, substituting this in Eq. (86) we obtain the asymptotic form of φ 1 (ζ): where W 1 is a constant. This relation will be used in Sec. IV A.

C. Exact solutions for the zero modes
In this section we find exact solutions of Eqs. (63)-(66) under the three sets of initial conditions, Eqs. (78)-(80), thus obtaining the three zero modes relevant for the MBC equations derived from Eq. (62). The calculations are lengthy and many technical details are delegated to appendix A. Since the final expressions we obtain are considerably complicated, we support our analytic calculations by comparison to numerical solutions of the system of ODEs (63)-(66). To avoid cumbersome notations, we omit the superscript distinguishing between different zero modes in general expressions that apply for all modes.

Analytic solution
Our starting point is Eqs. (63)-(66). With some algebraic manipulations (see appendix A), these equations can be brought to the form: Notice the structure of this set of equations: Eq. (88) is an inhomogeneous first order ODE for ∂ ζ (ψ σ + ψ ρ ) and can thus be solved analytically. Once this solution is substituted in Eq. (89), this last equation becomes also an inhomogeneous first order ODE for ψ σ , and is thus also analytically solvable. Finally, substitution of the solutions to Eqs. (88) and (89) yields an inhomogeneous first order ODE for ψ E , which is again analytically solvable. In the rest of this section we follow this procedure to obtain general analytic expressions with appropriate constants of integrations, which are then determined by imposing the initial conditions, Eqs. (78)-(80). The general solution of Eq. (88) has the form: where the expression of C(ζ) is given by (A12). Integrating from ζ = −∞ to ζ (the sum ψ σ + ψ ρ does not diverge at ζ = −∞, see Eqs. (75) and (76)) and using the change of variables y = E 0 (x) to calculate the resulting integrals, we obtain the final form where The explicit steps leading to this expression are given in appendix A 1. Substituting Eq. (92) in Eq. (89), we obtain after integration and using the same variable transformation y = E 0 (x), the complicated expression (A29) given in appendix A 2. This expression can be written in the following compact form where a σ (ζ) is given by Eq. (A32) and b σ (ζ) is a regular function of ζ whose exact form does not affect our analysis since we will consider later the limit |ζ ini | ≫ 1 and the last term of Eq. (94) will be as small as we want. Notice that if c φ = 0, as is the case for the linearly diverging mode [ψ ρ , ψ φ ], the first term on the RHS of Eq. (94) diverges as the auxiliary parameter ζ ini → −∞. As we show below, this singular dependence on ζ ini cancels out and does not appear in the final form of the MBC equation extracted from Eq. (62) corresponding to this mode. Since in deriving Eq. (62) in Sec. III B 2 we assumed, however, that ψ σ (0) is finite, we must carry the ζ ini -dependent terms that appear in Eq. (94) and other related expressions throughout our analysis.
Substituting Eqs. (92) and (94) (whose full form is given by Eq. (A29)) in Eq. (90), integrating and using again the variable transformation y = E 0 (x) we obtain the cumbersome expression (A41) given in appendix A 3 (notice that to simplify this expression, we used for the first time the explicit form (12) of α(x)). This expression can also be written in the following compact form similarly to Eq. (94) where a E (ζ) is given by Eq. (A44) and b E (ζ) is a regular function of ζ whose exact form does not affect our analysis.

Comparison with numerical results
To check the analytic formulas (92)

General expression for the relevant coefficients
In this section, we use the analytic solutions derived in Sec. III C 1 to provide compact expressions for the coefficients A, A + ψ E (0) and B, required for deriving the MBC equations from Eq.

(A57) (see apeendix A 4). This expression can be written in the compact form
where a A is given by Eq. (A59) and the exact form of b A is not required for our analysis. The exact form of the coefficient A + ψ E (0) is given in Eq. (A60) which is obtained using the expressions (96) and (95) for A and ψ E respectively (see Eqs. (A41) and (A57) for their complete expressions). The compact form for A + ψ E (0) is given by Notice that this last expression does not diverge when ζ ini → −∞. The term a A + a E (0) is given by Eq. (A62). Using the definition of B, Eq. (58), and the expression (95) for ψ E (ζ), we obtain the complicated form of B, Eq. (A68). The compact form for B is given by where a B is given by Eq. (A70) and the exact form of b B is again not required our analysis.

D. Completing inner region analysis
We here complete the inner region analysis by returning to Eq. (62) and using the results obtained in Sec. III C to derive MBC equations for the curved front.
First, let us rewrite Eq. (62) in the form where we used Eqs. (67), (77) and (93) (taking already the limit ζ c → −∞). To further facilitate our analysis, let us first write the expressions of the coefficients A+ ψ E (0) and B in the following form The coefficients γ j and β j , j = 1, · · · , 4, depend only on E + and can be easily extracted from Eqs. (A62) and (A70): (105) The contributions of last the term on the RHS of Eqs. (100) and (101) can be made arbitrarily small by choosing |ζ ini | sufficiently large and can thus be neglected in the rest of our analysis. The first MBC equation is obtained by using the values of the nondiverging zero mode corresponding to the initial condition (78) in Eqs. (99), (100) and (101). We obtain where we took the limit ζ c → −∞. The second MBC equation is obtained in a similar way by using the values of the second nondiverging zero mode corresponding to the initial condition (79) in Eqs. (99), (100) and (101). Using Eq. (110), we obtain As was pointed out already in Eq. (85), this last equation implies that the electrostatic field does not vanish behind the front. Notice also that, as mentioned in the discussion following Eq. (84), it is now easy to verify explicitly the identity Using the nondiverging zero modes to extract MBC equations (110) and (111) from (99) was relatively straightforward since all terms are finite at both end points ζ = 0 and ζ = ζ c → −∞. Naively, one may anticipate that since ψ E (ζ c ) ∼ ζ c it is impossible to extract physically meaningful MBC equation corresponding to this mode from (99). We show however below that taking appropriately the asymptotic limits ζ c → −∞ and the auxiliary parameter |ζ ini | ≫ 1 leads to cancellation of the diverging terms and to an additional nontrivial MBC equation.
In order to derive this third MBC equation, using the initial condition (80), we first substitute (87) in Eq. (99) and use Eqs. (100) and (101) to obtain where c φ is given by (80) (using Eq. (67)). Notice that the divergence of ψ E and φ 1 with ζ c cancels out provided |ζ c | is large enough such as E 1 (ζ c ) has reached its plateau value. We can now take the limit ζ c → −∞ such as the limits in Eqs. (110) and (111) are reached and can be used in Eq. (113). We then also notice that the two terms proportional to ζ ini in Eq. (113) cancel. After some manipulations, our third MBC equation can be written as where In Fig. 5, we present the dependence of the coefficients in Eqs. (32), (111)    section. In Fig. 6, we present the dependence of the coefficients in Eq. (110) on the variable E + . Notice the structure of the MBC equations (110), (111) and (114) derived above. As we will show below, the last two equations give rise to the MBC equations (126) and (127), whereas the first one determines the excess density of ionized matter behind the curved front which is thus "slaved" by the dynamics specified by Eqs. (124)-(128).
Let us discuss Eqs. (110), (111) and (114), in view of results obtained in Ref. [42] for planar front solution. For purely planar front, the only possible perturbation of the quantities (σ 0 , ρ 0 and E 0 ) are due to modification of the electric field ahead of the front E 0 (0) = E ∞ . As already discussed in Sec. III B 5, in our formalism this is equivalent to keeping in Eqs. (110) and (114) only the terms proportional to v 1 . In this case we also have E 1 (−∞) = 0, which implies that Eq. (87) reduces to A first test concerns our Eq. (110). In Ref. [42] the following relation was derived with E 0 (0) = −E + . This relation actually follows directly from the definition of σ − (25). Since for infinitesimal perturbations dσ 0 (−∞) = σ 1 (−∞) and dE 0 (0) = E 1 (0) we obtain from Eq. (118) Keeping on the RHS of Eq. (110) only terms proportional to v 1 and using (103) we recover Eq. (119). The second test stems from comparing our MBC Eq. (114) to the planar front relation Eq. (32), which for infinitesimal perturbations, φ 1 and E 1 , recast the form A simple computation leads to (121) Keeping again only terms proportional to v 1 on the RHS of Eq. (114), using Eq. (117) and comparing with (120) implies that the following equality must be satisfied. To check this we use the definition (115) together with (102) to obtain Using Eq. (21), integrating by parts the second integral (using d x (1/ρ(x)) = α(x)/ρ(x) 2 ) and taking the appropriate limit of boundary terms, Eq. (123) becomes identical, as required, to Eq. (121). These last two tests demonstrate consistency of our calculation with the previous results of [42].

A. The new moving boundary conditions
In this section, we derive a refined version of the boundary model which includes curvature correction. We first write the equations of the model and we discuss it below. It consists of Laplace equations for both outer regions and the moving boundary conditionŝ where the coefficients Q i depend on the electrostatic field ahead the front, and are given by analytic formulas derived from the planar front solution as follows and where ρ, σ − and ω 2 are given by (21), (25) and (116) respectively. From inner analysis performed in the previous sections, and as already discussed in Sec. II C, we know that the electron and ion densities vanish precisely at the discontinuity line of the front part of the ionization front where the matching with the outer region Ω + has to be done. Since there is no charge in Ω + the only consistent equation for the electric potential in this region is given by (124). On the other hand, we know also that the space charge density vanishes at the back part of the inner region (see Eqs. (26) and (60)). Consequently Eq. (125) for the outer region Ω − is also consistent with the inner analysis. Another possible equation for the electric potential in Ω − is given by ∇ 2 φ = σ − ρ with the necessary boundary condition σ − ρ = 0 on Γ − in order to be consistent with the inner analysis. However this equation introduces new unknowns, the electron and ion densities in Ω − , for which additional evolution equations would be needed. In Sec. V, we discuss the constant conductivity approximation introduced by Eq. (125).
Eq. (126) has required us to take the finite conductivity of the ionized region into account (Eq. (125)); it expresses the field behind the ionization front in terms of the field ahead of it, of the curvature and of the conductivity in the ionized region.

B. The field inside a propagating streamer
Eq. (126) has an interesting direct physical consequence for streamers whose radius and field enhancement evolve slowly during propagation within a comoving coordinate frame. Indeed, simulations of single streamers show that the field inside the streamer is constant in very good approximation (see for example Fig. 4 in [33]). Such a constant interior field is a common ingredient of phenomenological streamer modeling. Here we can calculate this field: on the streamer axis we find That the field inside the streamer might approach some fixed constant value, is a common interpretation of experimental results. Here we have derived an explicit relation for this field; more precisely, it depends on the conductivity inside the streamer, on the front velocity and on the curvature. This explicit prediction will be discussed further and tested on simulations elsewhere [46].

C. Large field limit
Since the width of the front ℓ − α (v 0 ), see Eqs. (29) and (30), is a monotonic decreasing function of v 0 , we expect that the sharp interface will be a better approximation for large fields. Moreover, as seen in Fig. 5, the coefficients Q i involved in the boundary conditions (126) and (127) become linear in this regime which simplifies the model. In this section, we compute the asymptotic behaviors of the coefficients Q i .
From the definition (129) of Q 0 (v 0 ) we have Consequently, we simply have Using the same method, the limit of the constant conductivity behind the front (25) takes the form since We then deduce that the asymptotic behavior of Q 2 (131) is given by The computation of the asymptotic behavior of Q 1 is a bit more lengthly. The expression of this coefficient is given by five terms. Let us compute the limit for each one and use the notation q i with i = 1, 2, · · · , 5 with Q 1 = 5 i=1 q i . Using the results obtained above we easily find that The third term, actually contained in the expression of β 1 (106), is given by The factor in front of the integral tends to 1 for large v 0 . The asymptotic behavior of the integralq 3 is obtained as above: (146) We then find The fourth term is given by To obtain the asymptotic behavior we also use the derivative: We thus obtain Which means that q 1 + q 2 + q 3 + q 4 ∼ O((v 0 ) 0 ) for large v 0 . Since we know from numerics that the asymptotic behavior of Q 1 is linear, this behavior is encoded in the last term: With the help of the following changes of variables x = v 0 q and y = v 0 p we get From the definition of ρ(x) (21) we have and from the definition of α(x) (12) we know that where θ(z) is the Heaviside's function. We then get Finally the asymptotic behavior of Q 1 (v 0 ) is given by For large electric fields (E + > ∼ 4), the boundary conditions (126)-(128) reducê

5.
Advance the curve (x f (y, t), y) in a normal directionn (towards Ω + ) at a raten · ∇φ + found in step 4. Eqs. (25), (26) show that the degree of ionization dµ α(µ) behind a uniformly translating front depends on the electric field strength E + . In the appendix of [47], it is shown that precisely the same degree of ionization σ − = σ − (E + ) from Eq. (25) is reached behind a planar front even if it does not propagate uniformly, if only the electric field E + is the same at the moment when the electron discontinuity passes the point of observation. Behind a slightly curved front, the degree of ionization will be approximately the same.
As a consequence, the conductivity further behind the front varies as a function of evolution history, as the immobile ions with their positive charge act as a memory. E.g., when a streamer emerges from an avalanche in a homogeneous background field E back , the conductivity in the evolving channel will increase from σ − (|E back |) immediately after the avalanche phase to σ − (|E enh (t)|) where E enh (t) is the enhanced field at the tip of the streamer at time t. This increasing electron and ion density along the axis of the streamer finger is typically seen in simulations, e.g., in Fig. 1 of [48] for an emerging streamer or in Fig. 3(a) of [13] for a fully developed long streamer.
When an electric field is applied to a body with spatially varying conductivity, in general, space charges can form in the interior, since charge conservation ∂ t q+∇·j = 0 with q = ρ − σ and j = σE can be rewritten in the form in agreement with (9)-(11), i.e., the space charge ρ − σ will decay on the Maxwell time scale 1/σ only if there is no conductivity gradient ∇σ in the direction parallel to the electric field E. However, in practice such space charge effects are typically small if the evolution is slow and consequently |∇σ| is small.

B. The approximation of ideal conductivity
Lozansky and Firsov in their Russian textbook and in their short English article [29] suggested to neglect the weak field in the ionized interior completely together with possible weak space charge effects, and to approximate the streamer interior as ideally conducting, i.e., to assume infinite conductivity and Inspection of simulations show that this is a valid approximation for the weak interior fields in comparison to the strong fields in the exterior. Therefore the approximation was adapted in the model [32,33,36,37], and its validity was further tested on simulations in section IV of [33].

C. The approximation of constant conductivity
A shortcoming of the ideal conductivity approximation is that it does not allow for a potential drop along the length of the streamer which is an important observable in experiments. In fact, the electric field along the axis of a long streamer is frequently assumed to be constant; it is an important point of scientific discussion and an ingredient of coarse grained models for complete streamer trees.
The approximation of constant conductivity in the interior is justified if the width of the streamer and the field enhancement at its head vary relatively slowly during propagation as it is the case in many simulations. We therefore here adapt the approximation of constant conductivity σ − in the ionized streamer interior. According to Eq. (160), space charges ρ − σ then decay on the time scale 1/σ − . For a front propagating with velocity v, this corresponds to a spatial decay on the length ℓ − α = v/σ − (30), generalizing the explicit asymptotic behavior (28) derived for uniformly translating planar fronts. Behind this layer, there are no space charges and the Laplace equation approximates the ionized interior as well: From a mathematical point of view, a finite conductivity in the interior might lead to an additional regularization mechanism. We have a model with two Laplacian fields on both side of the interface. This type of model is called Muskat, or Muskat-Leibenzon, problem in hydrodynamics and describes the evolution of the interface between two immiscible fluids in a porous medium or Hele-Shaw cell under applied pressure gradients or fluid injection/extraction. The so-called Hele-Shaw problem is thus the one-phase version of the Muskat problem, see for example [49,50,51,52].

VI. CONCLUSION AND OUTLOOK
We have derived the curvature corrections to a moving boundary approximation of streamer ionization fronts within the minimal model. We also have introduced a constant conductivity approximation in the ionized streamer interior. The matching procedure between inner and outer region followed the standard procedure at the back end of the front. However, ahead of the front, the matching region shrunk to a line. The solvability analysis was based on integrating from −∞ up to this line. It involved three zero modes of the adjoint operator, one of them diverging. Nevertheless, the boundary approximation could be carried through, and even fully analytically.
The analysis has two important merits. First, the curvature correction to the potential jump φ + − φ − across the interface might remove non-regular solutions [38] found in the "minimal regularized boundary model" [32,33,36,37]. Second, it provides a very practical relation between the front curvature and the electric field behind the front.
The derived boundary approximation accounts for propagating streamers in fairly homogeneous fields and far from electrodes. On the other hand, ionization fronts driven by monopole fields (e.g., close to sharp needle electrodes) or during branching undergo fast temporal changes of the electric field E + on the ionization front. They cannot be approximated by uniformly translating fronts anymore; and time dependent corrections need to be taken into account. In a forthcoming paper, we will explore how far charge conservation principles will carry us towards a full derivation of a moving boundary approximation for the minimal streamer model.

Acknowledgments
We acknowledge inspiring discussions with W. van Saarloos. The work of F.B. was supported by The Netherlands Organisation for Scientific Research (NWO) through Contract No. 633.000.401 within the program "Dynamics of Patterns," and the guest visit of B.D. to CWI Amsterdam by NDNS+ cluster.

APPENDIX A: SOLVING THE ZERO MODES EQUATIONS
In this appendix, we construct the exact expression for the zero modes. The system of ODE we need to solve for the zero modes is, according to Eqs. (63)-(66) the following Summing Eq. (A1) with Eq. (A2) leads to We will derive below an equation for ψ σ + ψ ρ . Once we get it, the complete system is solved with the help of (A4) and (A3). For this purpose, we take the derivative of Eq. (A2) and replace ∂ ζ ψ E by its expression given by (A3) to get Taking the derivative of (A4) and using the fact that (see Eqs. (34) and (44)) we obtain We finally multiply (A5) by E 0 /v 0 and sum it with (A7) to get the equation for This is formally a first order ODE and can be solved exactly. Once ψ σ + ψ ρ is known, ψ σ is computed with (A4) and then ψ E is obtained from (A3).

Expression of ψσ + ψρ
To solve (A8), we use the ansatz ∂ ζ (ψ σ + ψ ρ )(ζ) = C(ζ)e g(ζ) , with (A9) . (A10) The equation for C reads After integration we obtain, However from the asymptotic analysis of the system of ODE for the zero mode (see (75) and (76)), we know that The solution we search takes thus the form The explicit form of g(ζ) is computed using (A10) Using the change of variable we obtain where we use the notation ρ(|E 0 (ζ)|) = ρ(E 0 (ζ)) ≡ ρ 0 (ζ) in this Appendix. Substituting Eq. (A19) in Eq. (A14) we obtain (A20) With the change of variable (A17) we obtain where Again, with the change of variable (A17) one can compute that Consequently, the final form for the solution of Eq. (A8) is given by

Expression of ψσ
We have now to solve Eq. (A4). We rewrite it using the exact solution (A24). In particular we have Remark that it is easy to verify with the help of this relation that we can recover the asymptotic behavior (75).
After an integration of (A27) we obtain From the asymptotic analysis of (A27), and as we show explicitly below, one can verify that only the first integral of this last expression diverges when ζ ini → −∞. Using the change of variable (A17) we finally get In this last expression only the first integral diverges when ζ ini → −∞. The singularity can be explicitly ex-tract thanks to the relation where we used Eq. (22). The integral of this last relation is always finite. Expression (A29) can then be written as in the main text where The expression b σ (ζ) is not useful for our analyze since the contribution of this term in Eq. (A31) can be made as small as we want since |ζ ini | can be as large as we want (but finite).

Expression of ψE
Equation (A3) can be written as where we used the expression (A27) for ∂ ζ ψ σ . Using the relation (45) for f ′ 0 , this last expression reduces to where we used also the relation (20) between σ 0 and ρ 0 . It is also easy to see that in the limit ζ → −∞ we recover the asymptotic form (77). The expression of ψ E is obtained after an integration of (A34) Using the change of variable (A17) we get The two last terms diverge when ζ ini → −∞ but their sum is finite. With the help of the definition (22) we can write this sum as a finite quantity. Indeed, we know that and we substitute it in (A36) to get x 2 x 0 dy (v 0 + y) ρ(y) 2 − σ 0 (ζ)ψ σ (ζ) The only part of this relation that diverge when ζ ini → −∞ is contained in the expression of ψ σ (ζ). We can still simplify this last relation by using the following formulas Since we know that E 0 , σ 0 and ρ 0 decay exponentially for ζ → −∞, see (28), and since the zero modes diverge at most linearly in this limit, we can already replace ζ c by −∞ to compute A and B.
We know the exact expression of the zero modes and we start to compute the relevant quantities for the relation (62) which will give us the boundary conditions we search for. We use in this section the notation (ψ σ + ψ ρ )(ζ) = c σρ + c φ h(E 0 (ζ)), h(E 0 (ζ)) = 1 v 0 From the definition (57) of A and Eq. (A46), we have Using (34) we can compute the derivation of the ion density of the planar front solution We substitute it in Eq. (A48) to obtain The two last integral can be simplified. Indeed, using an integration by parts we have where we have still to replace ∂ ζ ψ σ and ∂ ζ E 0 by their expression (see (A27) and (34)). The last integral of (A50) can also be simplified by an integration by parts Notice that the only term that diverge when ζ ini → −∞ is contained into ψ σ (0). Now this formula can be simplified further since actually the two last integrals cancel. Indeed, with an integration by parts we get where we used the notation It is easy to verify that H(0) = 0, consequently we have The final expression of A is then rather simple and reads This expression can be written as in the main text where a A = σ 0 (0)a σ (0) + c φ

Expression of A + ψE(0)
We just need to sum the expression (A57) of A with the expression (A41) of ψ E evaluated at ζ = 0. We get with c ζini defined by (A42). This expression can be written as in the main text where a A + a E (0) =c ζini + c σρ v 0 α(E + ) − c φ α(E + ) and wherec ζini is defined by (A45).

Expression of B
Using the change of variables (A17) in the definition (58) of B we get Replacing ψ E by its exact form (A41) we get Now we need to compute The repeated integrals can be computed by parts to get For the first integral in Eq. (A66), we use Eq. (A30) in order to extract the part which diverges as ζ ini → −∞. Now we substitute Eq. (A66) in Eq. (A64), and in the final result, in order to eliminate a maximum of double integral, we use the following relation