Laplacian instability of planar streamer ionization fronts - an example of pulled front analysis
aa r X i v : . [ n li n . PS ] J un Laplacian instability of planar streamer ionizationfronts — an example of pulled front analysis ∗ Gianne Derks † Ute Ebert ‡ Bernard Meulenbroek § October 24, 2018
Abstract
Streamer ionization fronts are pulled fronts propagating into a linearly unstable state;the spatial decay of the initial condition of a planar front selects dynamically one specificlong time attractor out of a continuous family. A stability analysis for perturbations in thetransverse direction has to take these features into account. In this paper we show how toapply the Evans function in a weighted space for this stability analysis. Zeros of the Evansfunction indicate the intersection of the stable and unstable manifolds; they are used todetermine the eigenvalues. Within this Evans function framework, a numerical dynamicalsystems method for the calculation of the dispersion relation as an eigenvalue problem isdefined and dispersion curves for different values of the electron diffusion constant and ofthe electric field ahead of the front are derived. Numerical solutions of the initial valueproblem confirm the eigenvalue calculations. The numerical work is complemented with ananalysis of the Evans function leading to analytical expressions for the dispersion relationin the limit of small and large wave numbers. The paper concludes with a fit formula forintermediate wave numbers. This empirical fit supports the conjecture that the smallestunstable wave length of the Laplacian instability is proportional to the diffusion lengththat characterizes the leading edge of the pulled ionization front.
Keywords : Pulled front, stability analysis, streamer ionization front, dispersion relation,wave selection of Laplacian instability.
AMS subject classifications : 37L15, 34L16, 35Q99. ∗ G. Derks acknowledges a travel grant of the Royal Society, which initiated this research, and a visitorgrant of the Dutch funding agency NWO and the NWO-mathematics cluster NDNS + to finish the work. Thework was also supported by a CWI PhD grant for B. Meulenbroek. † Department of Mathematics, University of Surrey, Guildford, Surrey, GU2 7XH, UK([email protected]). ‡ Cluster ‘Modelling, Analysis and Simulation’, Center for Mathematics and Computer Science (CWI), P.O.Box 94079, 1090 GB Amsterdam, and also at Eindhoven Univ. Techn., The Netherlands ([email protected]). § While contributing to this paper, at CWI Amsterdam, now at Faculty of Electrical Engineering, Math-ematics and Computer Science, Delft Univ. Techn., P.O. Box 5031, 2600 GA Delft, The Netherlands([email protected]). aplacian instability of pulled streamer ionization fronts A streamer is the first stage of electric breakdown in large volumes, it paves the way of sparksand lightning, but also occurs without successive breakdown in phenomena like sprite dis-charges above thunderclouds or in corona discharges used in numerous technical applications.Recent reviews of relevant phenomena can be found in [20, 45]. Considered as a nonlinearphenomenon, the streamer is a finger shaped ionized region that propagates by self generatedfield enhancement at its tip into nonionized media. It has multiple scales as described in [20];as a consequence one can investigate a hierarchy of models on different levels of refinementthat are reductions of each other, starting from the reduction from a particle to a continuummodel [32] to the reduction from a continuum model to a moving boundary model [9] up tothe formulation of effective models for complete multiple branched streamer trees withoutinner structure that are known as “dielectric breakdown models” [38, 39, 40, 7]. All these re-ductions are the subject of current research; the present paper analyzes the stability of frontsin the continuum model; the resulting dispersion relation provides a test case for movingboundary approximations.Specifically, simulations of the simplest continuum model for negative streamers [15, 16, 47]have established the formation of a thin boundary layer around the streamer head. This layeris an ionization front that also carries a net negative electric charge. (Positive streamerswith positive net charge occur as well, but are not the subject of the present study.) Theconfiguration of the charge in a thin layer leads to the above mentioned field enhancement atthe streamer head that creates high ionization rates and electron drift velocities and hence letsthe streamer rapidly penetrate the non-ionized region. More recent numerical investigationsshow that the boundary layer or front can undergo a Laplacian instability that generates thestreamer branch [4, 42, 36, 37]. (We remark that an additional interaction mechanism incomposite gases like air somewhat modifies this picture [33] while the present analysis appliesto negative streamers in simple gases like pure nitrogen or argon.)
The streamer can be considered as a phenomenon where an ionized phase is separated froma non-ionized phase by a moving thin front. This concept [24, 4] implies that streamersshow similar features as moving boundary problems like viscous fingers, solidification frontspropagating into undercooled liquids, growth of bacterial colonies or corals in a diffusive fieldof food etc. Quantitative predictions within such models require a proper understanding ofthe front dynamics, in particular, of their stability against perturbations in the transversaldirection. This stability determines whether perturbations of the front position will grow orshrink, and on the long term whether the streamer will branch or not. As a first insight, onewould therefore like to analyze the stability of planar fronts against transversal perturbations,more specifically, the growth or shrinking rate s ( k ) of a linear perturbation with transversalwave length 2 π/k .The ionization front in the model for a negative streamer in a pure gas as treated in[15, 16, 47, 24, 25, 4, 42, 36, 37], including electron diffusion, creates a so-called pulledfront that has a number of peculiar mathematical properties: ( i ) for each velocity v ≥ v ∗ ,there is a dynamically stable front solution where the stability is conditional on the spatial aplacian instability of pulled streamer ionization fronts z → ∞ (where z is the spatial variable along the front); ( ii ) theconvergence towards this front is algebraically slow in time [21, 22]; ( iii ) this slow dynamicsis determined in the leading edge of the front that in principle extends up to z → ∞ and in thedynamically relevant space it will cause Fredholm integrals in the linear stability analysis todiverge, therefore curvature corrections cannot be calculated in the established manner [23],( iv ) the unconventional location of the dynamically relevant region ahead of the front alsorequires particular care in numerical solutions with adaptive grid refinement [37]. For thecalculation of the dispersion relation, which can be phrased as an eigenvalue problem for s ( k ) ,these features pose two challenges: first, the condition on the one-dimensional dynamicalstability and algebraic convergence properties, which are typical for pulled fronts, will leadto an apparently degenerate eigenvalue problem. Second, in a neighborhood of the origin,the dispersion curve s ( k ) is near the continuous spectrum. Hence numerical calculations ofthe eigenvalue problem with finite difference, collocation or spectral methods often lead tospurious eigenvalues. A dynamical systems method involving stable and unstable manifoldsavoids this problem. The stable and unstable manifolds are at least two-dimensional and anexterior algebra approach is employed to calculate the manifolds accurately.In [17, 4, 3], the treatment of pulled fronts and more-dimensional stable/unstable mani-folds was circumvented by neglecting the electron diffusion that acts as a singular perturba-tion. In this way, the leading edge of the front together with its mathematical challenges isremoved and the eigenvalue problem can be solved using shooting on the one-dimensional sta-ble/unstable manifolds. The resulting problem is characterized by two length scales, namelythe length scale 2 π/k of the transversal perturbation, and the longitudinal length scale ofelectric screening through the front that will be denoted by ℓ α . The dispersion relation inthis case shows a quite unconventional behavior, namely a short wave length instability whoseconsequences are further investigated in [35, 19]. In the present paper, we analyze the dis-persion relation including diffusion, mastering the above challenges and deriving quantitativeresults through a combination of analytical and numerical methods. The Evans function is an analytic function whose zeros correspond to the eigenvalues of aspectral problem, usually a linearization about a coherent structure like a front or solitarywave. It was first introduced in [26] and generalized in [1]. In the last decade, the Evansfunction has been applied in the context of many problems and various extensions and gener-alizations have been found, see the review papers [30, 44] and references in there. One of thefirst uses of the Evans function in the analysis of a planar front can be found in [46], whichanalyzes the stability of a planar wave in a reaction diffusion system arising in a combustionmodel. In the current paper we will show how pulled fronts can be analyzed with the Evansfunction by using weighted spaces in its definition.To define the Evans function, one writes the eigenvalue problem as a linear, first order,dynamical system with respect to the spatial variable z . Along the dispersion curve s ( k ) , thedynamical system has a solution which is bounded for all values of z . This can be phrasedin a more dynamical way as: the manifold of solutions which are exponentially decaying for z → + ∞ (stable manifold) and the manifold of solutions which are exponentially decayingfor z → −∞ (unstable manifold) have a non-trivial intersection along the dispersion curve. aplacian instability of pulled streamer ionization fronts s and k , which vanishes if thestable and unstable manifolds have a non-trivial intersection. Hence the Evans function canbe viewed as a Melnikov function or a Wronskian determinant, see also [29] or references inthere.In case of a pulled front, the definition of the stable manifold, and hence the Evansfunction, is not straightforward. The temporal stability of the asymptotic state of the pulledfront at + ∞ is conditional on the spatial decay of the perturbation. So this decay conditionshould be included in the definition of the stable manifold, otherwise the dimension of thismanifold might be too large. We will show that this condition can be built in the definitionof the stable manifold by considering the stable manifold in a weighted space. The Evansfunction is defined by using the weighted space for the stable manifold. Hence the dispersioncurve s ( k ) can be found as a curve of zeros of this Evans function. In section 2, we recall the model equations and the construction and properties of planarfronts. In particular, we summarize the multiplicity, stability, dynamical selection and con-vergence rate of these pulled fronts. In section 3, the stability of these fronts is investigatedas an eigenvalue problem for the dispersion relation s ( k ) of a linear perturbation with wavenumber k . The dispersion relation depends on the far electric field E ∞ and the electrondiffusion D as external parameters. In the stability analysis of the pulled ionization fronts,a constraint is imposed on the asymptotic spatial decay rate of the perturbations. Thisconstraint corresponds to the decay condition for the one-dimensional stability, but has tobe chosen quite subtly to avoid problems with the algebraic decay of the front solution. Aconsequence of the decay condition is that the eigenvalue problem (dispersion relation) issolved in a weighted space. In this weighted space, the apparent degeneracies have disap-peared, the stable and unstable manifolds of the ODE related to the eigenvalue problem arewell-defined and intersections of those manifolds are determined by using the Evans function.In section 3.4, dispersion relations for positive s are derived numerically for a number ofpairs of external parameters ( E ∞ , D ) . The numerical implementation of the Evans functionuses exterior algebra to reliably solve for the higher dimensional stable and unstable man-ifolds. In section 4, the numerical dispersion relation is tested thoroughly and confirmedwith numerical simulations of the initial value problem for the complete PDE model forthe particular values ( E ∞ , D ) = ( − , .
1) where D = 0 . E ∞ = − E ∞ , D ) analytically or a larger range of ( E ∞ , D )numerically.In section 5, explicit analytical asymptotic relations for the dispersion relation s ( k ) arederived for the limits of small and large wave numbers k . For k = 0 , two explicit eigenfunc-tions are known (which are related to the translation and gauge symmetry in the problem).These explicit solutions lead to expressions for the solutions on the stable manifold for smallwave numbers. The interaction between the slow and fast behavior on this manifold leadsto an asymptotic dispersion relation for small k . For large wave numbers, the eigenvalueproblem for the dispersion relation is dominated by a constant coefficient eigenvalue problem.An eigenvalue exists only if this constant coefficient system has no spectral gap. Using expo-nential dichotomies and the roughness theorem, the asymptotics of the dispersion relation is aplacian instability of pulled streamer ionization fronts k fits thedata very well, while the asymptotic limit for large k is not yet applicable in the range where s ( k ) is positive. After a discussion of relevant physical scales, we suggest a fit formula joiningthe analytical small k asymptotic limit with our physically motivated guess. This formulafits the numerical data well for practical purposes and strongly supports the conjecture thatthe smallest unstable wave length is proportional to the diffusion length that determines theleading edge of the pulled front. In this section we describe the streamer model and summarize the features of planar ionizationfronts as solutions of the purely one-dimensional model as a preparation for the stabilityanalysis in the dimensions transversal to the front. In particular, we recall the multiplicity ofthe front solutions that penetrate a linearly dynamically unstable state, and the dynamicalselection of the pulled front.
We investigate negative fronts within the minimal streamer model, i.e., within a “fluid ap-proximation” with local field-dependent impact ionization reaction in a non-attaching gas likeargon or nitrogen [24, 25, 17, 4, 42]. The equations for this model in dimensionless quantitiesare ∂ t σ − D ∇ σ − ∇ · ( σ E ) = σ f ( | E | ) , (2.1) ∂ t ρ = σ f ( | E | ) , (2.2) ∇ · E = ρ − σ , E = −∇ φ , (2.3)where σ is the electron and ρ the ion density, E is the electric field and φ is the electrostaticpotential. For physical parameters and dimensional analysis, we refer to discussions in [24,25, 17, 4, 42]. The electron current is approximated by diffusion and advection − D ∇ σ − σ E .The ion current is neglected, because the front dynamics takes place on the fast time scaleof the electrons and the ion mobility is much smaller. Electron–ion pairs are assumed tobe generated with rate σf ( | E | ) = σ | E | α ( | E | ) where σ | E | is the absolute value of electrondrift current and α ( | E | ) the effective impact ionization cross section within a field E . Hence f ( | E | ) is f ( | E | ) = | E | α ( | E | ) . (2.4)For numerical calculations, we use the Townsend approximation α ( | E | ) = e − / | E | [24, 25, 17,4, 42]. For analytical calculations, an arbitrary function α ( | E | ) ≥ α (0) = 0 and therefore f (0) = 0 = f ′ (0) . We will furthermore assume that α ( | E | ) is monotonically increasing in | E | , this is a sufficient criterion for the front to be apulled one [22]. The electric field can be calculated in electrostatic approximation E = −∇ φ .Mathematically, the model (2.1)-(2.3) describes the dynamics of the three scalar fields σ , ρ and φ . It is a set of reaction-advection-diffusion equations for the charged species σ and ρ coupled nonlinearly to the Poisson equation of electrostatics. aplacian instability of pulled streamer ionization fronts It follows immediately from (2.1)-(2.3) that there can be two types of stationary states of thesystem, one characterized by σ ≡ E ≡ f ( | E | ) = 0 implies | E | = 0 .).The stationary state with σ ≡ σ , there is no temporal evolution for σ ≡ ρ hasan arbitrary spatial distribution. The electric field E = −∇ φ then is determined by thesolution of the Poisson equation −∇ φ = ρ and by the boundary conditions on φ . In certainionization fronts in semiconductor devices [43], it is essential that the equivalent of ρ does notvanish in the non-ionized region. In the gas discharges considered here, on the other hand, itis reasonable to assume that the non-ionized initial state with σ ≡ ρ ≡ E ≡ E ≡ ∇ · E = 0 follows immediately, and therefore electron and iondensities have to be equal σ = ρ . In the absence of a field, the electrons diffuse ∂ t σ = D ∇ σ while the ions stay put ∂ t ρ = 0 . Therefore, these densities only can stay equal if ∇ ρ = 0 .Simulations [15, 16, 47, 24, 25, 4, 42, 36, 37] show that this occurs typically only if ρ ishomogeneous (though counter examples can be constructed). An ionization front separates such different outer regions: an electron-free and non-conductingstate with an arbitrary electric field E ∞ ahead of the front from an ionized and electricallyscreened state with arbitrary, but equal density σ − = ρ − of electrons and ions. In particular,we are interested in almost planar fronts propagating into a particle-free region ρ = σ = 0(where therefore ∇ φ = 0 ), and we study negative fronts, i.e., fronts with an electron surplusthat propagate into the electron drift direction towards an asymptotic electric field E ∞ < ∇ φ = −∇ · E = 0 that the electric field ahead of the frontis homogeneous.We assume that the front propagates into the positive z direction; the electric field aheadof a negative front then is E → E ∞ ˆ z , E ∞ < z → ∞ . (Here ˆ z is the unit vector in the z -direction.) It is convenient to introduce the coordinate system ( x, y, ξ = z − vt ) movingwith the front velocity v = v ˆ z . A planar, uniformly translating front is a stationary solutionin this co-moving frame, hence it depends only on the co-moving coordinate ξ , and will bedenoted by a lower index 0. A front satisfies D∂ ξ σ + ( v − ∂ ξ φ ) ∂ ξ σ + σ ( ρ − σ ) + σ f = 0 ,v∂ ξ ρ + σ f = 0 ,∂ ξ φ + ρ − σ = 0 , (2.5)where f = f ( | E | ) . This system can be reduced to 3 first order ordinary differential equa-tions. First, due to electric gauge invariance, the system does not depend on φ explicitly, butonly on E = − ∂ ξ φ . Using the variable E instead of φ reduces the number of derivativesby one. Second, electric charge conservation ∂ t q + ∇ · j = 0 can be rewritten in co-moving aplacian instability of pulled streamer ionization fronts − v∂ ξ q + ∂ ξ j = 0 . Therefore it can be inte-grated once − vq + j = c , ∂ ξ c = 0 . In the present problem, the space charge is q = ρ − σ and the electric current is j = − D∂ ξ σ − σ E . Furthermore, as there is a region withvanishing densities σ = 0 = ρ ahead of the front, the integration constant c vanishes inthis region, and therefore everywhere. Thus the planar front equations (2.5) can be writtenas D ∂ ξ σ = v ( ρ − σ ) − E σ ,v ∂ ξ ρ = − σ f ( | E | ) , (2.6) ∂ ξ E = ρ − σ , where ∂ ξ φ = − E decouples from the other equations. The planar front equations implythat E ( ξ ) < ξ when E ∞ < σ ρ E ξ → + ∞ → E ∞ and σ ρ E ξ →−∞ → σ − σ − , (2.7)and the electrostatic potential φ connects φ − (for ξ → −∞ ) with − E ∞ ξ + φ + (for ξ → + ∞ ). The ionization density σ − behind the front and the electrostatic potential difference φ + − φ − have to be determined for arbitrarily chosen electric field E ∞ ahead of the frontand for arbitrary, but sufficiently large, front velocity v . (We remark that only the potentialdifference φ + − φ − matters due to the gauge invariance of the electrostatic potential as oneeasily verifies on the equations.) The fronts can be constructed as heteroclinic orbits in athree-dimensional space as demonstrated in [25].The diffusion constant D is obviously a singular perturbation. For D = 0 , the frontequations can be solved analytically [25, 3], i.e., one can find explicit expressions for theparticle densities σ [ E ] , ρ [ E ] and for the front coordinate ξ [ E ] as a function of theelectric field E . For the negative fronts treated here, the front is continuous as function of D and the limit D → D = 0 , while for positivefronts ( E ∞ > The non-ionized state ( σ, ρ, E ) = (0 , , E ∞ ) with a nonvanishing electric field E ∞ is linearlyunstable under the temporal dynamics of the PDE (2.1)-(2.3). In fact, this spatial regionahead of the front dominates the dynamics, cf. the discussion in [25, 22]. Therefore, forfixed E ∞ , there is a continuous family of uniformly translating solutions, parametrized bythe velocity v ≥ v ∗ [24, 25, 21, 22], where v ∗ ( E ∞ ) = | E ∞ | + 2 p D f ( | E ∞ | ) . (2.8)The dynamics of uniformly translating fronts with velocity v > v ∗ are dominated by a flatspatial profile in the leading edge of the front σ v ( ξ ) ξ →∞ ∼ e − λξ with λ < Λ ∗ = r f ( | E ∞ | ) D , (2.9) aplacian instability of pulled streamer ionization fronts v and decay rate λ are related through v ( E ∞ , λ ) = | E ∞ | + Dλ + f ( E ∞ ) λ , (2.10)and therefore v ( E ∞ , λ ) > v ∗ ( E ∞ ) ≡ v ( E ∞ , Λ ∗ ) for λ = Λ ∗ . The spatial profile (2.9) with λ < Λ ∗ cannot build up dynamically from some initial condition with larger λ ; and it willdestabilize if perturbed with an initial condition with smaller λ , therefore such flat and fastfronts can be approached dynamically only by initial conditions with exactly the same profile(2.9) in the leading edge. For a thorough discussion of this dynamics, we refer to [22].In practice, the continuum approximation for the electron density breaks down for verysmall densities in the leading edge and the initial electron distribution satisfies a decay con-dition of the form lim ξ →∞ σ ( x, y, ξ, t = 0) e λξ = 0 for all λ < Λ ∗ , (2.11)if the penetrated state is really non-ionized. Therefore the velocity v ∗ is called the “selected”one, because it is the generic attractor for most physical initial conditions. Mathematicallyspeaking, the profile with velocity v ∗ (the selected front) is the only profile that can build updynamically from steeper initial conditions.Therefore the condition (2.11) on the spatial decay of the initial electron distributionexcludes all front solutions with velocity v > v ∗ as long time attractors of the dynamics. Ifthe criterion (2.11) is satisfied, then the selected front with speed v ∗ is dynamically stableand is approached with the universal algebraic convergence rate in time [21, 22] v ( t ) = v ∗ − ∗ t + O (cid:18) t / (cid:19) . (2.12)However, without the spatial decay condition on the initial condition, the selected front isformally not stable (although this is physically irrelevant). This will lead to specific problemsand solutions in the transverse stability analysis presented in the next section.The spatial profile of the electron distribution in the selected front is σ v ∗ ( ξ ) ξ →∞ ∼ ( aξ + b ) e − Λ ∗ ξ , a > . (2.13)To summarize, if the analysis is restricted to initial conditions with a sufficiently rapid spatialdecay in the electron densities (2.11), then the fronts have only one free external parameter,namely the field E ∞ ; it determines the asymptotic front velocity (2.8) and profile (2.13) aftersufficiently long times. Furthermore, the equivariance in the system gives that the positionof the front and its electrostatic potential are free internal parameters. The spatial decay behind the front will be important in the analysis as well, therefore werecall the basic behavior. For ξ → −∞ , the electron density approaches σ v ∗ ( ξ ) ξ →−∞ = σ − + c e λ − ξ , c > , (2.14)and the electric field decays with the same exponent E ( ξ ) = − ( c/λ − ) e λ − ξ . For D = 0 , σ − ( E ∞ , D = 0) = Z | E ∞ | α ( x ) dx (2.15) aplacian instability of pulled streamer ionization fronts D > σ − decreases by a correction of order of D , more precisely, σ − ( E ∞ , D ) = σ − ( E ∞ ,
0) + O ( D ) , σ − ( E ∞ , D > < σ − ( E ∞ ,
0) (2.16)was proved in the appendix of [32]. The eigenvalue λ − is given by λ − = √ v ∗ + 4 Dσ − − v ∗ D , (2.17)where both v ∗ and σ − depend on E ∞ and D . For small D , λ − can be approximated as λ − = σ − v ∗ + O ( D ) = Z | E ∞ | α ( x ) dx | E ∞ | + O ( √ D ) . (2.18)As an illustration, the spatial profiles of electron and ion density and the electric field ofthe selected front solution for a range of fields E ∞ and diffusion constants D are plotted inFigure 1. σ ρ E=-10 σ ρ
E=-5 σ ρ
E=-1 0 0.5 1 1.5 2 2.5-50 -40 -30 -20 -10 0 10 20 30 40 σ ρ
D=0.1 σ ρ
D=0.01 σρ D=0-1-0.9-0.8-0.7-0.6-0.5-0.4-0.3-0.2-0.1 0-50 -40 -30 -20 -10 0 10 20 30 40E=-10 E=-5 E=-1 -1-0.8-0.6-0.4-0.2 0-50 -40 -30 -20 -10 0 10 20 30 40D=0.1 D=0.01 D=0
Figure 1: The pulled planar front solutions on the left for varying E ∞ = − − −
10 andfixed D = 0 . E ∞ = − D = 0 . .
01 and 0 . Theupper panels show scaled electron and ion densities σ ( ξ ) /σ − ( E ∞ , D ) and ρ ( ξ ) /σ − ( E ∞ , D )and the lower panels the corresponding scaled electric fields E ( ξ ) / | E ∞ | as a function ofthe spatial coordinate ξ . The fronts are displayed in a staggered way. The normalizationfactors σ − ( E ∞ , D ) in the upper panels are σ − ( − ,
0) = 0 .
149 , σ − ( − , .
01) = 0 .
148 , σ − ( − , .
1) = 0 .
144 , σ − ( − , .
1) = 2 .
832 , σ − ( − , .
1) = 7 .
169 . aplacian instability of pulled streamer ionization fronts First we will introduce the transversal perturbation setting and discuss an apparent degener-acy of the dispersion relation. However, it turns out that the constraint on the spatial decay ofthe electron density “selects” a single dispersion relation for every far field E ∞ . This relationthen is calculated numerically based on dynamical systems techniques involving intersectionsof stable and unstable manifolds. Results for different fields E ∞ and diffusion constants D are presented. Suppose that there is a linear transversal perturbation of the uniformly translating front σ ( x, y, ξ, t ) = σ ( ξ ) + δ σ ( x, y, ξ, t ) + O ( δ ) , ξ = z − v ∗ t, (3.1)and similarly for ρ and φ . The linearized equation for σ , ρ and φ follows from Eqs. (2.1)-(2.3). By decomposing the perturbations into Fourier modes in the transversal directions x and y , by using isotropy in the transversal ( x, y ) plane and by using a Laplace transformationin t , the ansatz ( σ , ρ , φ ) = e ikx + st ( σ k , ρ k , φ k )( ξ ) (3.2)can be used for each Fourier component. The resulting equation can be written as a linearfirst order system of ODEs, using the extra variables τ k = ∂ ξ σ k and E k = − ∂ ξ φ k . Introduce w = ( τ k , σ k , ρ k , E k , φ k ) and the linear system is given by ∂ ξ w = M ( ξ ; E ∞ , k, s ) w , (3.3)with M = − E + v ∗ D σ − ρ − f + s + Dk D − σ D − ∂ ξ σ − σ f ′ D
01 0 0 0 00 − f v ∗ sv ∗ σ f ′ v ∗ − − k − . In the matrix M , the abbreviated notations f = f ( | E | ) and f ′ = ∂ η f ( η ) (cid:12)(cid:12) η = | E | are used.For the terms with f ′ , we have used that E < E | E | = − M depends on k , but not on k itself, the matrix is invariant under thetransformation k → − k . Thus if s ( k ) = s ∗ , then also s ( − k ) = s ∗ and vice versa. Thereforeit is sufficient to determine the dispersion relation for k > k < k > | k | . As will be shown later, thedispersion relation is not an analytic function of k near k = 0 and its expansion near k = 0is linear in | k | .For future use, we remark that the linearization matrix M does not involve any ξ -dependent terms in the fourth and fifth row and implies that E k and φ k are related by E ′ k = − φ k . Thus the E k -component of any solution of the linearized system (3.3) can be aplacian instability of pulled streamer ionization fronts E k ( ξ ) = c e kξ + c e − kξ + 12 Z ξξ h e k ( ξ − η ) + e − k ( ξ − η ) i [ ρ k ( η ) − σ k ( η )] dη, (3.4)where the constants c and c are determined by the value of E k and φ k at ξ = ξ . The linearized problem (3.3) is a spectral problem with the spectral parameters s and k . Ifthe asymptotic matrices M ± ( E ∞ , k, s ) = lim ξ →±∞ M ( ξ ; E ∞ , k, s ) exist and are hyperbolic(i.e., no eigenvalues on the imaginary axis), then the system (3.3) has a bounded solutionif and only if the unstable manifold from ξ = −∞ and the stable manifold from ξ = ∞ have a non-trivial intersection. So we will focus in this section on determining the stable andunstable manifolds.The behavior of the unstable manifold at the back of the front is given by the asymptoticmatrix M − ( E ∞ , k, s ) = lim ξ →−∞ M ( ξ ; E ∞ , k, s ) = − v ∗ D σ − + s + Dk D − σ − D sv ∗ − − k − . For s > k = 0 , this matrix has two negative and three positive eigenvalues: ± k, sv ∗ , µ −± = − v ∗ D ± p v ∗ + 4 D ( σ − + s + Dk )2 D . (3.5)Thus the unstable manifold is three dimensional. We remark that µ − + ( s = 0 = k ) is identicalto the spatial decay rate λ − (2.17) behind the unperturbed front.Finding the stable manifold ahead of the front is less straightforward. Normally the stablemanifold ahead of the front would be characterized by the matrix lim ξ → + ∞ M ( ξ ; E ∞ , k, s ) .For s > s + Dk < f ( E ∞ ) this matrix exists and has two positive and three negativeeigenvalues: ± k, sv ∗ , − Λ ∗ ± r s + Dk D = − p f ( E ∞ ) ± √ s + Dk √ D , (3.6)Thus the stable manifold is three dimensional and a dimension count gives that the inter-section of stable and unstable manifold is generically one dimensional. So for small valuesof s and k , a continuous family of eigenvalues seems to exist. This feature is related to theinstability of the asymptotic state at + ∞ , to the continuous family of uniformly translat-ing solutions for all v ≥ v ∗ ( E ∞ ) , and to the instability of fronts against perturbations withsmaller spatial decay rate λ , as discussed in the previous section. The continuous family ofeigenvalues s for fixed wave number k is eliminated by applying the analysis only to frontswith a sufficiently rapid spatial decay (2.11). This condition will be imposed in the definitionof the stable manifold; it will make the spectrum discrete. aplacian instability of pulled streamer ionization fronts e w = D w , D = diag( e (Λ ∗ − β ) ξ , e (Λ ∗ − β ) ξ , , ,
1) (3.7)where β ∈ (0 , Λ ∗ ) will be fixed later and depend on k and Λ ∗ . The freedom in the choiceof β is reminiscent of the fact that the decay condition holds for any λ < Λ ∗ , but not for λ = Λ ∗ . The system for e w is e w ξ = f M ( ξ ; E ∞ , k, s ) e w , (3.8)with f M = D · M · D − + ( ∂ ξ D ) · D − = − E + v ∗ D + Λ ∗ − β σ − ρ − f + s + Dk D − σ D e (Λ ∗ − β ) ξ − ∂ ξ σ − σ f ′ D e (Λ ∗ − β ) ξ
01 Λ ∗ − β − f v ∗ e − (Λ ∗ − β ) ξ sv ∗ − − σ f ′ v ∗ − e − (Λ ∗ − β ) ξ − k − Note that if β = 0 , then the asymptotic matrix ahead of the front (at ξ = + ∞ ) does notexist because e Λ ∗ ξ σ ( ξ ) grows linearly in ξ according to (2.13). To get an asymptotic matrixahead of the front, it is necessary that 0 < β < Λ ∗ . In this case, the asymptotic matrix is f M + ( E ∞ , k, s ) = lim ξ →∞ f M ( ξ ; E ∞ , k, s ) = − Λ ∗ − β − f ∞ + s + Dk D ∗ − β − f ∞ v ∗ sv ∗ − k − where f ∞ = f ( | E ∞ | ) . The matrix f M + has the eigenvalues ± k, sv ∗ , and µ + ± = − β ± r s + Dk D . (3.9)Hence for s > < β < min(Λ ∗ , k p s/ ( Dk )) , there are two negative and threepositive eigenvalues. Thus the stable manifold of (3.8) is two dimensional.For the original unscaled system (3.3) this means that only the two-dimensional sub-manifold given by D − acting on the stable manifold of (3.8) is relevant for the transverseinstability. This submanifold will be called the stable manifold of (3.3) from now on. Withthis convention, the dispersion relation is a well-defined curve s ( k ) and the curve is such thatat s = s ( k ) , the linearized system (3.3) has a bounded solution which satisfies the spatialdecay condition (2.11). Note that for both asymptotic matrices f M + and M − , the eigenvalues ± k become a degenerate eigenvalue 0 at k = 0 . This leads to square root singularities andit can be expected that the dispersion relation s ( k ) will be a function of √ k = | k | for k issmall. This will be confirmed in section 5. aplacian instability of pulled streamer ionization fronts The occurrence of an intersection of the stable and unstable manifolds will be measuredwith the Evans function. Our numerical method to determine the dispersion curve as aneigenvalue problem is based on a definition of the Evans function in an exterior algebraframework and uses similar ideas as in [2, 11, 12, 8, 10, 14]. The approach of following thestable/unstable manifolds at ξ = ±∞ with a standard shooting method and checking theirintersection using the Evans function, works only if these manifolds are one-dimensional orhave co-dimension one; in the present model, this is the case in the singular limit D = 0 anda shooting method was used in [3]. Otherwise, any integration scheme will inevitably justbe attracted by the eigendirection corresponding to the most unstable (stable) eigenvalue.Exterior algebra can be used to avoid this problem for higher dimensional manifolds andto preserve the analytic properties of the Evans function. Recently, a different method tocalculate the Evans function for higher dimensional manifolds has been proposed in [27].This method uses a polar coordinate approach and looks like a more suitable method for veryhigh dimensional problems.To calculate the evolution of the two dimensional stable and three dimensional unstablemanifold in a reliable numerical way, we will use the exterior algebra spaces V ( C ) and V ( C ) , respectively. The advantage of these spaces is that in V l ( C n ) , an l -dimensionallinear subspace of C n can be described as a one-dimensional object, being the l -wedge productof a basis of this space. Also, the differential equation on R (or C ) induces a differentialequation on the spaces V l ( C ) : W ξ = M ( l ) ( ξ ; E ∞ , k, s ) W , W ∈ V l ( C ) . (3.10)Here the linear operator (matrix) M ( l ) is defined on a decomposable l -form w ∧ . . . ∧ w l , w i ∈ C , as M ( l ) ( w ∧ . . . ∧ w l ) := ( Mw ) ∧ . . . ∧ w l + . . . + w ∧ . . . ∧ ( Mw l ) (3.11)and it extends by linearity to the non-decomposable elements in V l ( C ) . General aspectsof the numerical implementation of this theory can be found in [2]. The general form of thematrices M (2) and M (3) can be found in the appendix.To determine the three-dimensional unstable manifold for ξ ∈ ( −∞ ,
0] , we will use (3.10)with l = 3 . Since the induced matrix M (3) ( ξ ; E ∞ , k, s ) inherits the differentiability andanalyticity of M ( ξ ; E ∞ , k, s ) , the following limiting matrix exists: M (3) − ( E ∞ , k, s ) = lim ξ →−∞ M (3) ( ξ ; E ∞ , k, s ) . The set of eigenvalues of the matrix M (3) ± ( E ∞ , k, s ) consists of all possible sums of threeeigenvalues of M ± ( E ∞ , k, s ) (see Marcus [34]). Therefore, for s > k = 0 , there is aneigenvalue ν − of M (3) − , which is the sum of the 3 positive eigenvalues of M − , i.e., ν − = k + sv ∗ − v ∗ D + p v ∗ + 4 D ( σ − + s + Dk )2 D (note that the subscript “ − ” in ν − refers to exponentially decaying behavior at −∞ , not tothe sign of ν − , which is obviously positive). The eigenvalue ν − is simple and has real part aplacian instability of pulled streamer ionization fronts M (3) − (as M − is hyperbolic). We denote theeigenvector associated with ν − as W − e , i.e., M (3) − W − e = ν − W − e . This vector can always beconstructed in an analytic way (see [31, pp. 99-101], [10, 12, 28]). In this case it is easy todetermine an explicit analytical expression for the eigenvector as M − is quite sparse. Theunstable manifold corresponds to the solution W − ( ξ ) of the linearized system (3.10) (with l = 3 ) which satisfies lim ξ →−∞ e − ν − ξ W − ( ξ ) = W − e .The stable manifold can be determined in a similar way. As indicated in the previoussection, the scaled system (3.8) will be used to determine the stable manifold. For the stablemanifold with ξ ∈ [0 , ∞ ) , we will use (3.10) with l = 2 and the scaled matrix f M . As before,the asymptotic matrix M (2)+ ( E ∞ , k, s ) = lim ξ →∞ f M (2) ( ξ ; E ∞ , k, s ) . exists. Now the eigenvalues of M (2)+ ( E ∞ , k, s ) consists of all possible sums of two eigenvaluesof f M ± ( E ∞ , k, s ) . Therefore, for s > k = 0 , M (2)+ has an eigenvalue ν + , which is the sumof the 2 negative eigenvalues of f M + , i.e., ν + = − r s + Dk D + k − β ! As before, this eigenvalue is simple and has real part strictly less than any other eigenvalue of M (2)+ . The eigenvector associated with ν + will be denoted by W + e , i.e., M (2)+ W + e = ν + W + e The stable manifold of the scaled system (3.8) corresponds to the solution W + ( ξ ) of thelinearized system (3.10) (with l = 2 and M = f M ) which satisfies lim ξ →∞ e − ν + ξ W + ( ξ ) = W + e .To get the stable manifold of the original unscaled system, the inverse scalings matrix D − ( ξ )has to be used. For arbitrary ξ ≥ V ( C ) is quitecomplicated, but we will only need the original stable manifold at ξ = 0 . And at ξ = 0 , thescalings matrix is the identity matrix. Hence at ξ = 0 , the scaled stable manifold and theoriginal stable manifold are the same and W + e (0) describes the stable manifold of (3.3) at ξ = 0 .With the stable and unstable manifold as found above, the Evans function can be definedas ∆( E ∞ , k, s ) = W − (0; E ∞ , k, s ) ∧ W + (0; E ∞ , k, s ) , s > , k = 0 . (3.12)Thus the Evans function ∆ is more or less the determinant of the matrix formed by a basisof the unstable manifold at ξ = 0 and a basis of the stable manifold at ξ = 0 . If this functionis zero, then the bases are linearly dependent, hence the two manifolds have a non-trivialintersection.We have focused on the case s > − Dk < s < To calculate the Evans function numerically, first the front solution has to be determinednumerically as it appears explicitly in the linearization matrix M ( ξ ; E ∞ , k, s ) . The front is an aplacian instability of pulled streamer ionization fronts β = min(Λ ∗ , k ) . For the stable manifold, the linearized equation on V ( C ) c W + ξ = hf M (2) ( ξ ; E ∞ , k, s ) − ν + ( E ∞ , k, s ) I i c W + , c W + ( ξ ) (cid:12)(cid:12) ξ = L ∞ = W + e ( E ∞ , k, s ) , is integrated from x = L ∞ to ξ = 0 , using the second order Gauss-Legendre Runge-Kutta(GLRK) method, i.e. the implicit midpoint rule. Here the scaling c W + ( ξ ) = e − ν + ξ W + ( ξ ) en-sures that any numerical errors due to the exponential growth are removed and c W + ( ξ ) (cid:12)(cid:12) ξ =0 = W + ( ξ ) (cid:12)(cid:12) ξ =0 is bounded. The eigenvector W + e ( E ∞ , k, s ) can be determined explicitly as wedgeproduct of the relevant eigenvectors of M + ( E ∞ , s, k ) thanks to the sparse nature of this ma-trix.For the unstable manifold, the linearized equation on V ( C ) c W − ξ = h M (3) ( ξ ; E ∞ , k, s ) − ν − ( E ∞ , s, k ) I i c W − , c W − ( ξ ) (cid:12)(cid:12) ξ = L ∞ = W − e ( E ∞ , s, k ) , is integrated from x = − L ∞ to ξ = 0 , also using the implicit midpoint rule and introducingthe rescaling c W − ( ξ ) = e − ν − ξ W − ( ξ ) to remove potential exponential growth. Again, theeigenvector W − e ( E ∞ , k, s ) can be determined explicitly as wedge product of the relevanteigenvectors of M − ( E ∞ , s, k ) .At ξ = 0 , the computed Evans function is (see (3.12))∆( E ∞ , k, s ) = W − (0) ∧ W + (0) = c W − (0) ∧ c W + (0) . (3.13)For s = 0 = k , the center-stable and the center-unstable manifold have a two-dimensionalintersection, due to the translation and gauge invariance, see section 5.1 for details. In orderto determine the dispersion curve, we start near k = 0 and s = 0 and then slowly increase k and determine for which s ( k ) the Evans function ∆( E ∞ , k, s ( k )) vanishes.The numerical errors in the calculation of the Evans function are mainly influenced bythe step size used in the numerical integration with the GLRK method and errors in thenumerically determined front. The numerical integration uses the step size δx = 0 .
01 . Wehave performed various checks with a decreased step size and this shows that the error in thevalue of s for fixed k is largest (order 10 − ) if k is small and decreases for larger k (order10 − ). The accuracy of the front has been checked and is such that the error in the frontgives a negligible error (compared to the error due to the error in the step size) in the valueof s ( k ) . It turns out that the scheme is not very sensitive to errors in the front (at least forthe E ∞ and D values considered).In the following sections, we will present data for the dispersion curve for varying electricfield E ∞ and diffusion coefficient D . A more detailed discussion of the data, relation withanalytical asymptotics and some empirical fitting can be found in section 6. First we consider how the dispersion curve depends on the electric field E ∞ ahead of the front,while the diffusion coefficient is fixed to D = 0 . aplacian instability of pulled streamer ionization fronts E ∞ = − E ∞ = − E ∞ = −
10 . The figure shows that the shape of the dispersioncurve stays similar, but the scales of s and k increase when E ∞ increases. The dispersioncurves can be characterized by the maximal growth rate s max and the corresponding wavenumber k max where s ( k max ) = s max as well as by the wave number k > s ( k ) = 0that limits the band 0 < k < k of wave numbers with positive growth rates. s ( k ) k E=-10E=-5E=-1 (a) E ∞ = − − −
10 and fixed D = 0 . s ( k ) kD=0D=0.01D=0.1 (b) Fixed E ∞ = − D = 0 . .
01 and 0 .
Figure 2: Dispersion curves s ( k ) : (a) for varying E ∞ and fixed D = 0 . E ∞ = − D . The pairs ( E ∞ , D ) shown are the same as in Fig. 1. The data forthe singular limit D = 0 is taken from [3]. Next we consider the effect of varying the diffusion coefficient D , while keeping the electricfield ahead of the front fixed at E ∞ = − D = 0 ), the dispersion curve stays positive and is monotonically increasing to the saturationvalue s ( k ) = f ( | E ∞ | ) / k → ∞ . Our numerics shows that if diffusion is present, thisis not the case anymore. This is not surprising as diffusion is a singular perturbation. InFigure 2(b), the dispersion curve is shown for D = 0 . D = 0 .
01 and D = 0 ; the data for D = 0 is taken from [3]. It shows that the growth rate s ( k ) has a maximum s max if diffusionis present and becomes negative for k larger than some k . Furthermore for decreasingdiffusion D , the maximal growth rate moves upward towards the saturation value f ( | E ∞ | ) / D = 0 . This suggests that some features of the dispersion curve behave regularly in D ,in spite of the fact that D is a singular perturbation. For example, for a finite wave numberinterval, the limit of the dispersion curves for D → D = 0 .However, the asymptotic profile for large values of the wave number is obviously singularin D . This duality can also be found in the front itself: the velocity and the profile of theionization density and the electric field of the uniformly translating negative front dependregularly on D = 0 , while the profile of the ionization density is singular, as discussed insection 2.4 and shown in Fig. 1. aplacian instability of pulled streamer ionization fronts In the previous section, we have determined the dispersion relation s ( k ) for transversal pertur-bations of ionization fronts as a temporal eigenvalue problem of the PDE system linearizedabout the uniformly translating planar front. Since we are dealing with pulled fronts (cf.sections 1 and 2.4), the problem is unconventional: both the velocity v ∗ of the uniformlytranslating planar front and the dispersion relation s ( k ) of its transversal perturbations areunique only if the spatial decay constraint (2.11) is imposed. Furthermore a longitudinallyperturbed planar front approaches its asymptotic profile and velocity algebraically slowly intime (2.12). Therefore it is worthwhile to test the predicted dispersion relation on directnumerical simulations of the corresponding initial value problem.In this section, we will therefore simulate the temporal evolution of a perturbed planarfront by numerically solving the full nonlinear PDEs (2.1)-(2.3), and we will determine the dis-persion curve from a number of simulations with perturbations with different wave vectors k .This is done for far field E ∞ = − D = 0 . k as U ( x, z, t ) ≈ U ( ξ ) + δ U ( ξ, t ) e ikx + st , ξ = z − v ∗ t, U = ( σ, ρ, φ ) . (4.1)If δ e st is small enough, the solution is in the linear regime, and s ( k ) can be determinedfrom the evolution of the perturbation after U ( ξ, t ) has relaxed to some time independentfunction. Therefore in the numerical simulations, we choose δ for each wave number k insuch a manner that t is large enough to extract meaningful growth rates and that δ e st issmall enough that the dynamics at the final time is still well approximated by the linearizationabout the planar front.Furthermore, an appropriate choice of the initial condition reduces the initial transienttime during which U ( ξ, t ) in the co-moving frame still explicitly depends on time t . Ideally,such an initial condition is of the form U ( x, z,
0) = U ( ξ ) + δ U ( ξ ) cos kx etc., where U is a solution of the linearized system (3.3). To find an approximation for U ( ξ ) , we usethat the instability acts on the position of the front, i.e., we write the perturbed front as U ( ξ + δe ikx + st ) ≈ U ( ξ ) + δ e ikx + st ∂ ξ U ( ξ ) . Therefore we choose U ( ξ ) = ∂ ξ U ( ξ ) and theinitial condition as U ( x, z,
0) = U ( z ) + δ ∂ z U ( z ) cos kx. (4.2)As ∂ ξ U ( ξ ) is a solution of the linearized system for k = 0 = s , this choice will be veryefficient for small values of k and require longer transient times for larger k .To solve the full 2D PDE, the algorithm as described in [4, 42] is used, while adaptivegrid refinement as introduced in [37] was not required. For fixed k , the PDE with initialcondition (4.2) is solved on the spatial rectangle ( x, z ) ∈ [0 , L x ] × [0 , L z ] . The length of thedomain in the transversal x -direction, L x , is such that exactly 5 wave lengths fit into thedomain, i.e., L x = πk , and periodic boundary conditions are imposed in this direction byidentifying x = 0 with x = L x . On the boundaries in the longitudinal z -direction, Neumannconditions for the electron density are imposed. The potential is constant far behind thefront and the electric field is constant far ahead of the front; therefore for the potential φ ,the Dirichlet condition φ = 0 is imposed at z = 0 , and the Neumann condition ∂ z φ = − E ∞ at z = L z accounts for the far field ahead of the front. aplacian instability of pulled streamer ionization fronts σ max ( x, t ) = max z ∈ [0 ,L z ] σ ( x, z, t ) (4.3)evaluated across the front. The reason is as follows. First, Figure 1 shows the spatial profilesof planar fronts for different electric fields E ∞ and illustrates that for fixed D , the maximumof the electron density σ max as well as the asymptotic density σ − behind the front stronglydepend on the field E + immediately ahead of the front, where for a planar front the closeand the far field are identical: E + = E ∞ . Second, the modulation of the front positionleads to a modulation of the electric field E + immediately before the front (cf. discussion insection 5.2); therefore σ max as a function of E + is modulated as well. A (a) The maximal value of the electron density σ max ( x, t ) for t = 50 as a function of the transversalcoordinate x . The perturbation has wave number k = 0 .
45 , the transversal length L x = 10 π/k leavesspace for 5 wave lengths that are clearly visible.
20 40 60 80 100t-10-5 l og ( A ) (b) The logarithm of the amplitude of the front mod-ulation log A as a function of time t for the same k . Figure 3: Examples of data of the initial value simulation from which the growth rate s ( k )shown in Fig. 4 are determined.An example of σ max ( x, t ) as a function of the transversal coordinate x for a fixed time t is plotted in Fig. 3(a). The amplitude of the wave modulation is determined by the Fourierintegral A ( t, k ) = k π Z πk σ max ( x, t ) cos kx dx. In Figure 3(b) we plot log A against time t for k = 0 .
45 . Note that k = 0 .
45 is close to k = 0 .
482 (see Figure 2(a) and Table 1) where the growth rate vanishes, s ( k ) = 0 , thereforethe growth rate in the present example is small and particularly sensitive to numerical errors.Figure 3(b) shows an initial temporal transient before steady exponential growth is reached(where exponential growth manifests itself as a straight line in the logarithmic plot). This istypically observed for the larger k -values ( k > . U ( z ) in the initial condition (4.2) is not optimal. For k < . U ( z ) ≈ ∂ z U ( z ) for small values of k . aplacian instability of pulled streamer ionization fronts s ( k ) , a least squares algorithm is used to fit the best linethrough the data points ( t, log A ) , and the initial transient time is ignored for larger valuesof k . For each value of k , the growth rate is determined with various choices of δ . Theresulting growth rate s ( k ) is indicated in Figure 4 with crosses X and the error bars arerelated to the various choices of δ . dynamical systemssimulation Figure 4: The dispersion curve s ( k ) for E ∞ = − D = 0 . × with errorbars indicate results of simulations of the full initial value problem as discussed in section 4and demonstrated in Fig. 3(a). For comparison, the results of the dynamical systems methodfrom section 3.4 are indicated with + symbols.Fig. 4 also shows the dispersion relation for ( E ∞ , D ) = ( − , .
1) determined with thedynamical systems method in the previous section 3.4; these numerical results are denotedwith + and can now be compared with the results of the initial value problem from thepresent section. Around the maximum of the curve, the agreement between the numericalresults of the two very distinct methods is convincing. For larger values of k , the differencesincrease, but the error bars of the initial value problem results increase as well. Furthermore,the plotted error bars are an underestimation as they only account for the errors discussedabove that emerge from the choice of the initial condition and from the time interval ofevaluation and therefore from possible initial transients and from a possible transition tononlinear behavior. Additional errors can be due to the numerical discretization and timestepping of the s themselves. We therefore conclude that the two results agree within thenumerical error range of the initial value simulations over the whole curve. aplacian instability of pulled streamer ionization fronts k ≪ and k ≫ Having determined the dispersion relation numerically for different values of electric field E ∞ and diffusion constant D in section 3, and having tested the correctness of the eigenvaluecalculation against numerical solutions of the initial value problem in section 4, we now willanalytically derive asymptotic expressions for the dispersion relation for small and large valuesof the wave modes k . It will be shown that these asymptotic limits are s ( k ) = k E ∞ dv ∗ dE ∞ , k ≪ − Dk , k ≫ k asymp-totic limit that was presented in [3] for the singular limit D = 0 , and we correct the resultproposed in [5]; and we also rigorously derive the large k asymptotic limit, in agreement withthe form proposed in [5]. k ≪ Translation invariance and electrostatic gauge invariance give two explicitly known boundedsolutions of the linearized system (3.3) at k = 0 and s = 0 . These are u ′ ( ξ ) = ( σ ′′ ( ξ ) , σ ′ ( ξ ) , ρ ′ ( ξ ) , E ′ ( ξ ) , − E ( ξ )) T and e = (0 , , , , T . Note that e is a solution of the linearized system (3.3) for k = 0 and s arbitrary.From the asymptotics of (3.3) for k = 0 = s at ξ = −∞ , we see that the only exponentiallydecaying solution at ξ = −∞ is given by u ′ ( ξ ) . This solution is related to the only positiveeigenvalue µ − + (see (3.5)). For ξ → + ∞ , the solution u ′ ( ξ ) → − E ∞ e , hence this solution isnot exponentially decaying for ξ = + ∞ . However, it is easy to obtain an explicit exponentiallydecaying solution at ξ = + ∞ , this is u ′ ( ξ ) + E ∞ e .From the eigenvalues in (3.5) it follows that for 0 < k ≪ < s ≪ ξ → −∞ involves one eigenfunction with a fast exponentialdecay (related to the eigenvalue µ − + ) and two eigenfunctions with a slow exponential decay(related to the eigenvalues k and sv ∗ ). Similarly, from the eigenvalues in (3.9), it followsthat the two-dimensional stable manifold at ξ → + ∞ involves one eigenfunction with a fastexponential decay (related to the eigenvalue − Λ ∗ + β + µ + − ) and one eigenfunction with a slowexponential decay (related to the eigenvalue − k ). Recall that the stable manifold is definedas a subset of the full stable manifold to account for the spatial decay condition (2.11).We focus on approximating an exponentially decaying solution on the stable manifold. Aswe have seen above, in lowest order, this solution is w s ( ξ ; E ∞ , ,
0) = u ′ ( ξ ) + E ∞ e + O ( k + s ) . To determine the next order, we will use the slow behavior of the asymptotic system andwrite w s ( ξ ; E ∞ , k, s ) = u ′ ( ξ ) + E ∞ (0 , , , k, e − kξ + kU s ,k ( ξ ) + sU s ,s ( ξ ) + O (( s + k ) ) . aplacian instability of pulled streamer ionization fronts E ∞ (0 , , , k, e − kξ ) is an approximation of theslow behavior on the stable manifold, while the other three terms are related to the fast decay.Because of the slow decay, the expansion is only valid on a ξ -interval with kξ = o (1) , hence ξ shouldn’t be too large.Substitution of these expressions into the linearized system (3.3) gives that U s ,s and U s ,k have to satisfy( D ξ − M ( ξ ; E ∞ , , U s ,s = (cid:16) σ ′ D , , ρ ′ v ∗ , , (cid:17) ;( D ξ − M ( ξ ; E ∞ , , U s ,k = − E ∞ (cid:16) − τ − σ f ′ D e − kξ , , σ f ′ v ∗ e − kξ , , − e − kξ (cid:17) . By analyzing the unperturbed system, we can find a particular solution of the first equation.The front solution u = ( τ , σ , ρ , E , φ ) satisfies τ ′ = − v ∗ + E D τ + ( σ − ρ ) σ D − σ f D σ ′ = τ ρ ′ = − σ f v ∗ E ′ = − ( σ − ρ ) φ ′ = − E Differentiating this system with respect to E ∞ gives( D ξ − M ( ξ ; E ∞ , , ∂ u ∂E ∞ = dv ∗ dE ∞ (cid:18) − τ D , , σ f ( v ∗ ) , , (cid:19) = − dv ∗ dE ∞ (cid:18) σ ′ D , , ρ ′ v ∗ , , (cid:19) . Hence U s ,s = − (cid:18) dv ∗ dE ∞ (cid:19) − ∂ u ∂E ∞ + a homogeneous solution.The asymptotic behavior of ∂ u ∂E ∞ is ∂ u ∂E ∞ ∼ (0 , , , , − ξ ) , ξ → ∞ . (5.1)The polynomial growth in the φ k -component will need to be canceled by the behavior of theother terms which involve k and hence will give a relation between s and k .In fact, the E k -components of the linearized system (3.3) for any s or k can be expressedby the integral equation (3.4). For all solutions on the stable manifold, the limit E k ( ∞ ) =lim ξ →∞ E k ( ξ ) is well-defined, so we can write the E k -components on the stable manifold as E sk ( ξ ) = c e − kξ − Z ∞ ξ h e k ( ξ − η ) + e − k ( ξ − η ) i [ ρ k ( η ) − σ k ( η )] dη, In the integral, σ k must satisfy the decay condition (2.11) on the stable manifold and hencewill have fast exponential decay. Furthermore, from (3.3), it can be seen that ρ k ( ξ ) = c e sξv ∗ + v ∗ Z ξξ e − s ( ξ − η ) v ∗ [ σ ( η ) f ′ ( η ) E k ( η ) − f ( η ) σ k ( η )] dη . As the term inside the integral hasfast exponential decay, we get that on the stable manifold c = 0 and ρ k has fast exponential aplacian instability of pulled streamer ionization fronts E sk has fast exponential decay for ξ large.Since φ k = − E ′ k , we get on the stable manifold for k small and ξ -values not too large, say ξ ∼ k − (hence kξ ∼ k )( E sk , φ sk )( ξ ) = ( k, e − kξ + O ( e − Λ ∗ ξ )= ( k, (cid:16) − kξ + O ( k √ k ) (cid:17) = ( k, − kξ ) + O ( k √ k ) . The exponentially decaying solution on the stable manifold is given by w s ( ξ ; E ∞ , k, s ) = u ′ ( ξ ) − s (cid:18) dv ∗ dE ∞ (cid:19) − ∂ u ∂E ∞ ( ξ )+ E ∞ (0 , , , k, e − kξ + kU s ,k ( ξ )+ O (( s + k ) ) , and the arguments above show that the order k contribution in the ( E k , φ k ) -componentsis given fully by E ∞ (0 , , , k, e − kξ and that kU s ,k ( ξ ) does not contribute to those com-ponents at this order. So it follows that the polynomial growth in the φ k -component of ∂ u ∂E ∞ ( ξ ) as given by (5.1) has to be canceled by the φ k -component in E ∞ (0 , , , k, e − kξ ,i.e., s (cid:16) dv ∗ dE ∞ (cid:17) − = kE ∞ or s = c ∗ k + O ( k ) , with c ∗ = E ∞ dv ∗ dE ∞ , (5.2)and v ∗ given in Eq. (2.8). Equation (5.2) establishes the small k limit of the dispersionrelation s ( k ) . k ≪ asymptotic limit There is also a physical argument for the asymptotic limit (5.2) that generalizes the calculationin section IV.C of reference [3] to nonvanishing
D > k ≪ π/k is the largest length scale of the problem. It is much largerthan the inner longitudinal structure of the ionization front. On the length scale 2 π/k , thefront can therefore be approximated by a moving boundary between ionized and non-ionizedregion at the position z f ( x, t ) = z + v ∗ ( E ∞ ) t + δ e ikx + st , (5.3)and the local velocity of this perturbed front is v ( x, t ) = ∂ t z f ( x, t ) = v ∗ ( E ∞ ) + s δ e ikx + st . (5.4)The electric field in the non-ionized region is determined by E = −∇ φ , where φ is thesolution of the Laplace equation ∇ φ = 0 together with the boundary conditions; these are E → E ∞ ˆ z for z → ∞ fixing the field far ahead of the front and φ ( ... ) = O ( k ) ≈ φ ( x, z, t ) = − E ∞ ( z − z − v ∗ t ) + E ∞ e − k ( z − z f ) δ e ikx + st + O ( δ ) , for z ≥ z f ,E + ( x, t ) = E ∞ + k E ∞ δ e ikx + st + O ( δ ) , for z = z f , (5.5) aplacian instability of pulled streamer ionization fronts E + ( x, t ) = lim ǫ ↓ E ( x, z f + ǫ, t ) is the electric field extrapolated onto the boundary fromthe non-ionized side.As the perturbation is linear, the front is almost planar δ ≪ π/k . Therefore it willpropagate with the velocity v ∗ ( E + ) = | E + | + 2 p Df ( | E + | ) (2.8) of the planar front in thelocal field E + . Inserting E + from (5.5) and expanding about E ∞ , we get v ( x, t ) = v ∗ ( E + ( x, t )) = v ∗ ( E ∞ ) + ∂ E v ∗ (cid:12)(cid:12)(cid:12) E ∞ k E ∞ δ e ikx + st + O ( δ ) . (5.6)Comparison of (5.6) and (5.4) immediately gives the dispersion relation s = c ∗ k + O ( k ) (5.2),that generalizes the result s ( k ) = | E ∞ | k + O ( k ) that was derived in [3] for the singular limit D = 0 . k ≫ The asymptotic limit for k ≫ k is large and that s + Dk is positive, but not small, i.e., k ≫ , s + Dk > , and s + Dk = o (1) (5.7)and show that this does not allow for bounded solutions. With the assumptions above on s and k , the dominant contributions in the matrix M on the whole axis ξ are M ∞ = s + Dk D sv ∗ − k − + O (1) . (5.8)Here the three entries − sv ∗ are necessary for a nonvanishing determinant.We want to use the Roughness Theorem [13] for exponential dichotomies to show thatfor k large and s not close to − Dk , the exponential dichotomy of the constant coefficientODE is close to the exponential dichotomy of the full system. So first we recall the definitionof an exponential dichotomy, which gives projections on stable or unstable manifolds. Definition 1 ([13])
Let A be a matrix in R n × n , u ∈ R n , and J = R − , R + , or R . Let Φ ( y ) be a solution matrix of the linear system dudy = A ( y ) u, y ∈ J. (5.9) The linear system (5.9) is said to possess an exponential dichotomy on the interval J if thereexist a projection P and constants K and κ s < < κ u with the following properties: | Φ ( y ) PΦ − ( y ) | ≤ K e κ s ( y − y ) , for y ≥ y , y, y ∈ J | Φ ( y )( I − P ) Φ − ( y ) | ≤ Ke κ u ( y − y ) , for y ≥ y, y, y ∈ J An extension for PDEs of this definition can be found in [41].The Roughness Theorem for exponential dichotomies states the following. aplacian instability of pulled streamer ionization fronts Theorem 2 (Roughness Theorem [13])
Consider the system dudy = [ A + A ( y )] u, (5.10) with A ∈ R n × n a hyperbolic matrix and u ∈ R n . Then for all δ > there exists a δ > such that for all matrix functions A : R → R n × n with k A k L ∞ ( R + , R n × n ) < δ , thesystem (5.10) has an exponential dichotomy on R + (and R − ) with its dichotomy exponentsand projections δ -close to those of dudy = A u (in the L ∞ ( R + , R n × n ) norm). A constant coefficient linear system does not have bounded solutions. So if the exponentialdichotomy of the linearized system (3.3) is close to the one of the constant coefficient systemwith the matrix M ∞ as in (5.8), the linearized system (3.3) does not have bounded solutionseither. We will show that this is the case if s and k satisfy the assumptions (5.7).First we introduce some scaling and coordinate transformations. Define the small pa-rameter ε = k and the scaled spatial variable, the transformation matrix and transformedvector η = kξ = ξε , T ( ε ) = diag( ε, , ε, ε,
1) and b w ( η ) = T ( ε ) w ( εη ) . Now (3.3) can be written as ∂ η b w = hc M ( ε, s ) + ε c M ( η ; E ∞ , ε ) i b w , (5.11)with c M ( ε, s ) = ε sD sεv −
10 0 0 − and c M ( η ; E ∞ , ε ) = − E + v ∗ D ε σ − ρ − f D − σ D − ∂ ξ σ − σ f ′ D
00 0 0 0 00 − ε f v ∗ − σ f ′ v ∗ − ε , where ξ = εη. The eigenvalues and eigenvectors of the constant coefficient matrix c M are ± , with eigenvectors w ± = (0 , , , ∓ , ± q ε sD , with eigenvectors w ± = (cid:18) ± q ε sD , , , , (cid:19) ; εsv ∗ , with eigenvector w = (0 , , , , . If | s | ≪ ε , then the matrix c M ( ε, s ) is not hyperbolic at ε = 0 . However, this problem isnot fundamental as it is known that there is a hyperbolic splitting in the full problem (see aplacian instability of pulled streamer ionization fronts ε is close to zero. The spectral gap disappears if s ≈ − Dk (or ε s ≈ − D ) and in this casethe following arguments will not work. The spectral gap allows us to define a weight functionwhich moves the spectrum away from the imaginary axis. To be specific, define e w ( η ) = e νη b w ( η ) , with ν = sgn( s ) , if | εs | < v ∗ , if | εs | ≥ v ∗ Then e w ( η ) satisfies the ODE ∂ η e w = (cid:16)hc M ( ε, s ) + ν I i + ε c M ( η ; E ∞ , ε ) (cid:17) e w (5.12)and the spectrum of c M ( ε, s ) + ν I is bounded away from zero for all ε small (as long as s + Dk is not small). The system ∂ η e w = hc M ( ε, s ) + ν I i e w (5.13)has an exponential dichotomy with projection P such that the range of P is the span of alleigenvectors of the negative eigenvalues and the kernel of P is the span of all eigenvectorsof the positive eigenvalues. Then P is the projection on the stable subspace of the linearsystem (5.13) and I − P is the projection on the unstable subspace.Clearly the matrix c M ( η ; E ∞ , ε ) is uniformly bounded for all η and ε small. Thusapplying the Roughness Theorem 2 gives that for ε small there is an exponential dichotomyfor the system (5.12) on R + with projection P sε ( η ) onto the stable subspace such that P sε ( η )is ε -close to P for all η ≥ R − with projection P uε ( η ) onto the unstable subspace such that P uε ( η ) is ε -close to I − P for all η ≤ P sε (0) is ε -close to the range of P andthe range of P uε (0) is ε -close to the range of I − P , so the range of P sε (0) and the range of P uε (0) have only a trivial intersection for ε small.As the weight function e νη has been chosen such that no eigenvalue crosses the imaginaryaxis, it affects only the value of the dichotomy exponentials, not their sign, nor the stableand unstable manifolds. Hence the stable and unstable manifolds of (5.11) have a trivialintersection only. And the same holds for the stable and unstable manifolds of (3.3), as theonly difference between the systems (5.11) and (3.3) is a scaling. So it can be concluded thatthe linear system (3.3) does not have any bounded solutions for ε small ( k large) and s notclose to − Dk .If s is close to − Dk , i.e., s = − Dε (1 + o (1)) , then the matrix b A ( ε, s ) has a positiveand negative eigenvalue of order o (1) and the spectral gap will disappear in the limit ε → s = − Dk (1 + o (1)) , there is a possibility forbounded solutions to exist. As the dispersion curve indicates a bounded solution of theODE (3.3), this implies that for ε near zero, hence k large, the dispersion curve is s ( k ) = − Dk (1 + o (1)) , k → ∞ . (5.14)So far we have used that s + Dk > s = − Dk + f ( | E ∞ | ) . However, one should aplacian instability of pulled streamer ionization fronts < β < Λ ∗ the edgebecomes s = − Dk + β . By taking the limit for β → s = − Dk is an edge of the continuous spectrum. Thus with the decay condition, either the dispersioncurve satisfies s ( k ) ≥ − Dk or it ends at the continuous spectrum. In section 3 we have derived dispersion relations for a number of fields E ∞ and diffusionconstants D by numerically solving an eigenvalue problem, we have confirmed these calcula-tions by a numerical solutions of the initial value problem in section 4, and we have derivedanalytical asymptotic limits to these dispersion relations in section 5. This sets the stage forcomparing the numerical results to the analytical asymptotic limits and for deriving physicallyguided empirical fits to the numerical dispersion curves where the analytical asymptotic limitsare not applicable. The small k -data derived in section 3.4 and the analysis of section 5.1are shown to be consistent in section 6.1. After showing in section 6.2 that a simple cross-over formula joining the asymptotic behavior of the small wave numbers with the asymptoticbehavior of the large wave numbers is not satisfactory, we give a data collapse, empirical fitsand arguments on relevant scales in section 6.3. k asymptotic limits First the asymptotic relation (5.2) for small k is tested on the numerical results. Beyondthe results visible in the plots of figure 2, the numerical dispersion relation for E ∞ = − D = 0 . k , and the result is shown in adouble logarithmic plot in Figure 5 that zooms in on the small k behavior. Also plotted isthe analytical asymptotic limit (5.2). The comparison between numerical data and analyticalasymptotic limit is convincing in the range of small k .For the other values of E ∞ and D presented in figure 2, the dispersion relation s ( k ) againis fitted very well by the asymptotic limit (5.2) for small values of k . This will be illustratedin Figure 6 below. It is quite suggestive to join the small k asymptotic limit (5.2) with the large k asymptoticlimit (5.14) into one cross-over formula s ( k ) = c ∗ k − Dk , c ∗ ( E ∞ , D ) = E dv ∗ dE (cid:12)(cid:12)(cid:12)(cid:12) E ∞ = | E ∞ | f ′ ( | E ∞ | ) s Df ( | E ∞ | ! . (6.1)A formula similar to s ( k ) = c ∗ k − Dk was suggested in [5], but with a different prefactorinstead of c ∗ . However, we now will confirm once more the correctness of the prefactor c ∗ ,and we will show that the large k asymptotic limit is not yet applicable in the range ofpositive growth rates s ( k ) .If (6.1) holds, then the dispersion relation for the rescaled variable S = Ds/c ∗ as afunction of the rescaled wave number κ = Dk/c ∗ becomes S = κ − κ . Therefore the formula(6.1) can easily be tested on the numerical data from figure 2 by plotting them in rescaledvariables S and κ with appropriate values for D and c ∗ ( E ∞ , D ) for each curve. The result is aplacian instability of pulled streamer ionization fronts -12-10-8-6-4-2 0-12 -10 -8 -6 -4 -2 0 l og ( s ) log(k) E=-1, D=0.1 Figure 5: A log - log plot of the dispersion curve s ( k ) for E = − D = 0 . k . Red squares and solid line: Data from the numerical evaluation ofthe eigenvalue problem. Blue dashed line: Analytical asymptotic limit log s = log c ∗ + log k according to (5.2).shown in Fig. 6, together with the parabola κ − κ . The plot illustrates that the asymptoticlimit (5.2) indeed is a very good fit to all data for small k .For larger k , the curves differ quantitatively. In particular, S vanishes for κ between0.014 and 0.035 for the numerical dispersion curves while the formula (6.1) predicts this tohappen for κ = 1 . Also the maximum of the dispersion curve S max is never higher than 0.0027for the numerical data while formula (6.1) predicts 0.25. Of course, this is not in contradictionwith the analytical results in Section 5.3 for large k . Rather it says that the positive part ofthe dispersion curve lies completely in the range of small k , where the asymptotic limit forlarge k is not applicable. We conclude that cross-over formula (6.1) is not an appropriate fitfor the numerically derived dispersion relations. We finish this section with a data collapse and arguments on relevant scales that guide em-pirical fits.
First, we investigate whether the numerical data for different E ∞ and D can be collapsedonto one curve. This is done by determining the maximum of the dispersion curve s max andthe wave number k where the growth rate vanishes, s ( k ) = 0 , from the numerical data foreach pair ( E ∞ , D ) . In figure 7 all curves are plotted as s/s max and k/k with their respective s max and k . The plot shows that the curve shapes are very similar, but they do not coincidecompletely. For example, there seems to be a small drift in the position of the maximum. aplacian instability of pulled streamer ionization fronts S κ E=-1, D=0.1E=-5, D=0.1E=-10, D=0.1E=-1, D=0.01
Figure 6: Labeled curves: The numerical dispersion curves from Figure 2 plotted as S = Ds/c ∗ over κ = Dk/c ∗ . Dotted line on the left: the parabola κ − κ that would bepredicted by (6.1) as far as it fits into the plotted region. D = 0 case In a second step, we investigate which physical or mathematical mechanisms can suppressthe growth rate s ( k ) for much smaller values of k than suggested by the large k asymptoticlimit (5.14). In a first overview, there are three length scales in the problem. The transversalperturbation is characterized by its wave length 2 π/k . In the longitudinal direction, the frontis characterized by two length scales, the electric screening length ℓ α and the diffusion length ℓ D , cf. Fig. 1, ℓ α = 1 α ( E ∞ ) and ℓ D = 1Λ ∗ = s Df ( E ∞ ) . (6.2)For vanishing diffusion D , the diffusion length vanishes, and the screening length ℓ α has tobe compared to the wave length of the perturbation; in [3] it was shown that it determinesthe cross-over from the small to the large k asymptotic limit of the dispersion relation: s ( k ) = (cid:26) | E ∞ | | k | , for | k | ≪ k α , | E ∞ | k α , for | k | ≫ k α ; for D = 0 , where k α = 12 ℓ α . (6.3)The actual curve for D = 0 and E ∞ = − s ( k ) = | E ∞ | k h (cid:16) kk α (cid:17) p i /p (6.4)for positive real p reproduces the asymptotic limits (6.3), but does not fit the full numericalcurve for E ∞ = − p . The functional form of (6.4) will servebelow as an inspiration for our empirical fits for D > aplacian instability of pulled streamer ionization fronts s / s m a x k/k E=-1, D=0.1E=-5, D=0.1E=-10, D=0.01E=-1, D=0.01
Figure 7: The numerical dispersion curves from Figure 2 plotted as s/s max over k/k ; here s max = max k s ( k ) and k > s ( k ) = 0 are determined from the respective curve.Searching for why (6.4) does not properly fit the data, one realizes that there are actuallytwo different definitions of the screening length possible:1 ℓ + α = α ( E ∞ ) and 1 ℓ − α = λ − = Z | E ∞ | α ( x ) dx | E ∞ | + O ( √ D ) , (6.5)with λ − from (2.17). The dimensions of both quantities are the same, and they approach eachother if α is constant in a large part of the integration interval [0 , | E ∞ | ] ; this is the case withthe Townsend approximation α ( x ) = e − /x for | E ∞ | ≫ ℓ + α characterizes theslopes of the fields near the discontinuity of σ [3], while ℓ − α characterizes the decay (2.14) ofthe fields far behind the front for ξ → −∞ . The analysis in [3] shows that linear perturbationswith wave numbers k ≫ ℓ + α , while smaller k could couple to the larger spatial structure characterized by ℓ − α , thisconjecture will be tested on the numerical data below and asks for future analysis. D > ℓ D emerges as another length scale in thefront. As illustrated in figure 1, instead of the discontinuous electron density in the front for D = 0 , a diffusive layer of width ℓ D = 1 / Λ ∗ (2.13) builds up in the leading edge. While D increases, the dispersion relation decreases as shown in figure 2(b). As the diffusive layeris the main new feature of the front for D > s is created within this boundary layer. The physical mechanism is that diffusion cansmear perturbations of short wave length out, hence suppressing their growth. This processmainly takes place in the diffusive layer because gradients are largest in this region. Thisidea has inspired an attempt in [5] to calculate s ( k ) by local analysis within the diffusionlayer. In principle, such an approach combined with proper matched asymptotic expansionscould work. However, the calculation in [5] was intrinsically inconsistent [18], disagrees with aplacian instability of pulled streamer ionization fronts k and therefore fits the numerical results even worse thanformula (6.1), cf. Fig. 6.We have tested whether the diffusion length ℓ D = 1 / Λ ∗ plays a role in the dispersionrelation by plotting the numerical data from Fig. 2 this time for the rescaled variables s = s/ ( c ∗ Λ ∗ ) over K = k/ Λ ∗ . The result is shown in Figure 8. It shows that the numericaldispersion curves are well approximated by s ( k ) ≈ c ∗ k + O ( k ) , for k → , (6.6) k ≈ Λ ∗ / , where s ( k ) = 0 , (6.7)The numerical evidence from Fig. 8 summarized in (6.7) together with the physical explana-tion above suggest the following conjecture. Conjecture 1
The largest unstable wave number of the Laplacian instability is proportionalto the inverse diffusion length.
We remark that the data gives k ≈ / (4 ℓ D ) , while the cross-over formula (6.1) would suggestthat k ≈ ℓ α /ℓ D , highlighting again its inadequacy for intermediate k values. Figure 8 alsoshows that the value of the wave number for which the maximum of the dispersion curve isattained, lies in the range of k max = 0 . k to 0 . k . s / Λ * c * k/ Λ *E=-1, D=0.1E=-5, D=0.1E=-10, D=0.1E=-1, D=0.01 (a) Using a = 3 /α ( | E ∞ | ) in formula (6.8). s / Λ * c * k/ Λ *E=-1, D=0.1E=-5, D=0.1E=-10, D=0.1E=-1, D=0.01 (b) Using a = 3 c ∗ /f ( | E ∞ | ) in formula (6.8). Figure 8: The data of Figure 2 plotted as s = s/ (Λ ∗ c ∗ ) over K = k/ Λ ∗ . The lines are givenby the empirical formula (6.8).The data in Figure 8 suggest an empirical formula of the form (for s ≥ s ( k ) = c ∗ k a k (cid:18) − k Λ ∗ (cid:19) or s = K a Λ ∗ K (1 − K ) , (6.8)where the parameter a will depend on the external parameters D and E ∞ . The factor c ∗ k creates the correct asymptotic limit (6.6) for k ≪ − k/ Λ ∗ ) creates thenon-trivial zero of the dispersion relation at k (6.7). The form of the numerator is inspiredby (6.4), and the proper asymptotic limit (6.3) for large k and D = 0 would be reached for a = 2 ℓ + α + O ( √D ) . Obviously, the empirical formula (6.8) is not valid in the asymptotic range k ≫ s < s ≈ − Dk . aplacian instability of pulled streamer ionization fronts K max , s max ) of (6.8), it follows that K / s max = 1 / a . (The number 1/4 directly stems from the factor 4 in (1 − K ) .) Thisrelation indeed fits the numerical curves quite well, therefore the factor 4 is supported twiceindependently. Relevant numerical data for this and other fits is collected in Table 1.The value for a is less obvious. The empirical formula (6.8) gives the following relationbetween a Λ ∗ and the maximum of curve1 − K max K = a Λ ∗ = 1 − √ s max s max . The empirical values for those quotients are given in Table 1.The limit of D = 0 and k ≫ a = 2 ℓ + α + O ( √D ) , but the fit is unconvincing(cf. also the discussion for D = 0 above). However, we found that a = 3 ℓ + α fits the datareasonably well. Formula (6.8) with this value of a together with the numerical data arevisualized in Figure 8(a). The fit is quite good for the lower two curves, the upper twodisplay some discrepancies. The main problem is the value of s max with a relative errorbetween 2% and 24% , while the position of K max has a relative error as low as 0 .
5% to 7% .As f ( | E ∞ | ) /c ∗ = α ( | E ∞ | ) + O ( √ D ) , another possible fit is a = 3 c ∗ /f ( | E ∞ | ) , it is displayedin Figure 8(b). The fit is quite good for the upper two curves, but now the fit has somediscrepancies for the lower two curves. And the value of K max has a larger error between 2%and 10% while the position of s max has a much smaller error of only 1% to 10% . Obviously,these observations ask for further analytical investigation. Note finally the striking relationbetween λ − = 1 /ℓ − α and the value of k max for larger values of the electric field in Table 1.As a basis for future work, all characteristic numerical data is collected in this table. In this paper, we have found dispersion curves for negative streamer ionization fronts bynumerically solving an eigenvalue problem, we have verified this prediction on the numericalsolution of an initial value problem, we have derived analytical expressions for the asymptoticsof the curve for large and small wave lengths, and we have presented a physically motivatedfit formula to the numerical curves for intermediate wave lengths. The investigation is ofinterest for two reasons: because pulled fronts like these ones are mathematically challengingto investigate, and because explicit predictions on the linear stability of ionization fronts helpto interpret numerical and experimental observations of propagating and branching streamerdischarges.The ionization front is a pulled front, i.e., the front is part of a family of traveling waves,which propagate into a temporally unstable steady state. For the dynamics with one spatialvariable, most traveling waves in this family are attractors only for waves with exactly thesame asymptotic decay profile. The exception is the pulled front, which has the steepestdecay of all waves in the family and is an attractor for waves with a sufficiently fast decay(therefore excluding the slower decay rates for the other traveling waves in the family). Theinstability of the state ahead of the front and the related spatial decay condition imply thatonly a submanifold of the stable manifold in the transverse instability problem is relevant forthe transverse instability analysis. This submanifold is identified by introducing a weightedsolution space that excludes solutions with a too slow decay rate. We have integrated the aplacian instability of pulled streamer ionization fronts E ∞ , D ) ( − , .
01) ( − , .
1) ( − , .
1) ( − , .
1) ( − , . s max k max k k /k max v ∗ c ∗ ∗ = p f ( | E ∞ | ) /D α ( | E ∞ | ) = 1 /ℓ + α σ − λ − = 1 /ℓ − α ∗ /α ( | E ∞ | ) 49.5 15.6 17.2 23.4 31.53Λ ∗ c ∗ /f ( | E ∞ | ) 55.5 21.6 21.7 27.0 34.8 K = k / Λ ∗ s max = s max / ( c ∗ Λ ∗ ) 0.0118(2) 0.0196(1) 0.0193(2) 0.0175(1) 0.0155(1) k max /k α K / s max − K max ) / (4 K ) 40(6) 17.7(1) 21(2) 22(2) 31(2)(1 − √ s max ) / s max ∗ that characterizes the leading edge of the pulled ionization front;though the evidence up to now is only numerical. As such behavior is physically reasonable,the next step would be to derive it analytically, e.g., by a local analysis in the diffusive layerand matched asymptotic expansions. Such an expansion could be based on the limiting casewhere the diffusion length ℓ D = 1 / Λ ∗ is much smaller than the screening length ℓ α .The calculated dispersion curves also contribute to understanding the stability of actualstreamers. Two- and three-dimensional time dependent simulations [4, 42, 36, 37, 33] of the aplacian instability of pulled streamer ionization fronts A Matrices in exterior algebra spaces
In this appendix, we give explicit expressions for the matrices M ( l ) acting on the exterioralgebra space V l ( C ) for l = 2 , e , . . . , e be the standard basis for C . Then aninduced basis on V ( C ) is given by a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e , a = e ∧ e . The matrix M (2) : V ( C ) → V ( C ) can be associated with a complex 10 ×
10 matrixwith entries such that M (2) a i = X j =1 M (2) ij a j , i, j = 1 , . . . , , (A.1)where, for any decomposable x = x ∧ x ∈ V ( C ) , M (2) x := Mx ∧ x + x ∧ Mx . Let M be an arbitrary 5 × M = m m m m m m m m m m m m m m m m m m m m m m m m m , (A.2)then M (2) takes the explicit form M (2) = d m m m − m − m − m m d m m m − m − m m m d m m m − m m m m d m m m − m m d m m − m − m − m m m d m m − m − m m m m d m m − m m − m m d m − m − m m − m m m d m − m m − m m − m m d where d ij = m ii + m jj . aplacian instability of pulled streamer ionization fronts M (3) : V ( C ) → V ( C ) can be associated with a complex10 ×
10 matrix. First we define an induced basis on V ( C ) by b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e , b = e ∧ e ∧ e . The matrix for M (3) ∈ C × has entries such that M (3) b i = X j =1 M (3) ij b j , i, j = 1 , . . . , , (A.3)where, for any decomposable x = x ∧ x ∧ x ∈ V ( C ) , M (3) x := Mx ∧ x ∧ x + x ∧ Mx ∧ x + x ∧ x ∧ Mx . If M is given by (A.2), then M (3) takes the explicit form M (3) = d m m − m − m m m m d m m − m − m m m m d m m − m − m − m m d m − m m m − m m m d m m − m − m m − m m d m m m − m m d m − m m m − m m m d m − m m − m m − m m d m m − m m m − m m d where d ijk = m ii + m jj + m kk . Acknowledgments
We thank Bj¨orn Sandstede for helpful discussions on the roughness theorem. We thankWillem Hundsdorfer and Ren´e Reimer for advice and help with the simulations in section 4.
References [1]
J. Alexander, R. Gardner & C.K.R.T. Jones . A topological invariant arising inthe stability analysis of traveling waves , J. Reine Angew. Math. (1990), pp. 167–212.[2]
L. Allen & T.J. Bridges . Numerical exterior algebra and the compound matrixmethod , Numerische Mathematik (2002), pp. 197–232.[3] M. Array´as and U. Ebert . Stability of negative ionization fronts: Regularization byelectric screening? , Phys. Rev. E (2004), article no. 036214 (10 pages).[4] M. Array´as, U. Ebert and W. Hundsdorfer . Spontaneous branching of anode-directed streamers between planar electrodes , Phys. Rev. Lett. (2002), articleno. 174502 (4 pages). aplacian instability of pulled streamer ionization fronts M. Array´as, M.A. Fontelos, and J.L. Trueba . Mechanism of branching in negativeionization fronts , Phys. Rev. Lett. (2005), article no. 165001 (4 pages).[6] A. Back, J. Guckenheimer, M.R. Myers, F.J. Wicklin & P.A. Worfolk . DsTool: Computer assisted exploration of dynamical systems , Notices Amer. Math. Soc. (1992), pp. 303-309.[7] A.D.O. Bawagan . A stochastic model of gaseous dielectric breakdown , Chem. Phys.Lett. (1997), pp. 325–331.[8]
K.B. Blyuss, T.J. Bridges and G. Derks . Transverse instability and its long-termdevelopment for solitary waves of the (2+1)-dimensional Boussinesq equation , Phys. Rev.E (2003), article no. 056626 (9 pages).[9] F. Brau, A. Luque, B. Meulenbroek, U. Ebert, L. Sch¨afer . Construc-tion and test of a moving boundary model for negative streamer discharges , preprint:arXiv:0707.1402 .[10]
T. Bridges, G. Derks & G. A. Gottwald . Stability and instability of solitary wavesof the fifth-order KdV equation: a numerical framework , Physica D (2003), pp. 190–216.[11]
L.Q. Brin . Numerical testing of the stability of viscous shock waves , Math. Comp. (2001), pp. 1071–1088.[12] L.Q. Brin & K. Zumbrun . Analytically varying eigenvectors and the stability of viscousshock waves , Seventh Workshop on Partial Differential Equations, Part I (Rio de Janeiro,2001). Mat. Contemp. (2002), pp. 19–32.[13] W. Coppel . Dichotomies in stability theory . Lecture Notes in Mathematics. Springer,1978.[14]
G. Derks & G. A. Gottwald . A robust numerical method to study oscillatory insta-bility of gap solitary waves , SIAM J. Appl. Dyn. Sys. (2005), pp. 140–158.[15] S.K. Dhali, P.F. Williams . Numerical simulation of streamer propagation in nitrogenat atmospheric pressure , Phys. Rev. A (1985), pp. 1219–1221.[16] S.K. Dhali, P.F. Williams . Two-dimensional studies of streamers in gases , J. Appl.Phys. (1987), pp. 4696–4707.[17] U. Ebert and M. Array´as , Pattern formation in electric discharges , pp. 270-282 in:
Coherent Structures in Complex Systems (eds.: Reguera, D. et al ), Lecture Notes inPhysics (Springer, Berlin 2001).[18]
U. Ebert and G. Derks , Comment on [5] , 1 page submitted to Phys. Rev. Lett.[19]
U. Ebert, B. Meulenbroek, L. Sch¨afer . Rigorous stability results for a Laplacianmoving boundary problem with kinetic undercooling , SIAM J. Appl. Math. (2007),pp. 292–310. aplacian instability of pulled streamer ionization fronts U. Ebert, C. Montijn, T.M.P. Briels, W. Hundsdorfer, B. Meulenbroek, A.Rocco, E.M. van Veldhuizen . The multiscale nature of streamers , Plasma SourcesScience and Technology (2006), pp. S118–S129.[21] U. Ebert, W. van Saarloos , Universal algebraic relaxation of fronts propagating intoan unstable state and implications for moving boundary approximations , Phys. Rev. Lett. (1998), pp. 1650–1653.[22] U. Ebert, W. van Saarloos . Front propagation into unstable states: Universal al-gebraic convergence towards uniformly translating pulled fronts , Physica D (2000),pp. 1–99.[23]
U. Ebert, W. van Saarloos . Breakdown of the standard Perturbation Theory andMoving Boundary Approximation for “Pulled” Fronts , Phys. Rep. (2000), pp. 139–156.[24]
U. Ebert, W. van Saarloos and C. Caroli , Streamer Propagation as a PatternFormation Problem: Planar Fronts , Phys. Rev. Lett. (1996), pp. 4178–4181.[25] U. Ebert, W. van Saarloos and C. Caroli , Propagation and Structure of PlanarStreamer Fronts , Phys. Rev. E (1997), pp. 1530–1594.[26] J.W. Evans . Nerve axon equations IV. The stable and unstable impulse , Indiana Univ.Math. J. (1975), pp. 1169–1190.[27] J. Humpherys and K. Zumbrun . An efficient shooting algorithm for Evans functioncalculations in large systems . Physica D (2006), pp. 116–126.[28]
J. Humpherys, B. Sandstede and K. Zumbrun . Efficient computation of analyticbases in Evans function analysis of large systems . Numerische Mathematik (2006),pp. 631-642.[29]
T. Kapitula , The Evans function and generalized Melnikov integrals . SIAM J. Math.Anal. (1999), pp. 273-297.[30] T. Kapitula and B. Sandstede . Eigenvalues and resonances using the Evans func-tion . Discrete and Continuous Dynamical Systems (2004), pp. 857-869.[31] T. Kato . Perturbation Theory for Linear Operators . Second Edition, Springer Verlag:Heidelberg (1984).[32]
C. Li, W.J.M. Brok, U. Ebert, J.J.A.M. van der Mullen . Deviations from thelocal field approximation in negative streamer heads , J. Appl. Phys. (2007), articleno. 123305 (14 pages).[33]
A. Luque, U. Ebert, C. Montijn, W. Hundsdorfer . Photoionisation in negativestreamers: fast computations and two propagation modes , Appl. Phys. Lett. (2007),article no. 081501 (3 pages).[34] M. Marcus . Finite Dimensional Multilinear Algebra, Part II , Marcel Dekker: New York(1975). aplacian instability of pulled streamer ionization fronts
B. Meulenbroek, U. Ebert, L. Sch¨afer . Regularization of moving boundaries in aLaplacian field by a mixed Dirichlet-Neumann boundary condition: exact results , Phys.Rev. Lett. (2005), article no. 195004 (4 pages).[36] C. Montijn, U. Ebert, W. Hundsdorfer . Numerical convergence of the branchingtime of negative streamers , Phys. Rev. E (2006), article no. 065401 (4 pages).[37] C. Montijn, W. Hundsdorfer, U. Ebert . An adaptive grid refinement strategy forthe simulation of negative streamers , J. Comp. Phys. (2006), pp. 801–835.[38]
L. Niemeyer, L. Pietronero, H.J. Wiesman . Fractal Dimension of Dielectric Break-down , Phys. Rev. Lett. (1984), pp. 1033–1036.[39] L. Niemeyer, L. Ullrich, N. Wiegart . The mechanism of leader breakdown in elec-tronegative gases , IEEE Trans. Electr. Insul. (1989), pp.309–324.[40] V.P. Pasko, U.S. Inan, T.F. Bell . Mesosphere-troposphere coupling due to sprites ,Geophys. Res. Lett. (2001, pp. 3821–3824.[41] D. Peterhof, B. Sandstede and A. Scheel . Exponential dichotomies for solitary-wave solutions of semilinear elliptic equations on infinite cylinders . Journal of DifferentialEquations, (1997) pp. 266–308.[42]
A. Rocco, U. Ebert and W. Hundsdorfer . Branching of negative streamers in freeflight , Phys. Rev. E, (2002), article no. 035102(R) (4 pages).[43] P. Rodin, U. Ebert, W. Hundsdorfer, I.V. Grekhov . Superfast fronts of im-pact ionization in initially unbiased layered semiconductor structures , J. Appl. Phys. (2002), pp. 1971–1980.[44] B. Sandstede . Stability of travelling waves . In: Handbook of Dynamical Systems II (BFiedler, ed.). North-Holland (2002), pp 983-1055.[45]
S.M. Starikovskaia . Plasma assisted ignition and combustion , J. Phys. D: Appl. Phys. (2006), pp. R265–R299.[46] D. Terman . Stability of Planar Wave Solutions to a Combustion Model , SIAM J. Math.Anal. (1990), pp 1139-1171.[47] P.A. Vitello, B.M. Penetrante, J.N. Bardsley . Simulation of negative streamerdynamics in nitrogen , Phys. Rev. E49