A moving boundary model motivated by electric breakdown: II. Initial value problem

A moving boundary model motivated by electric breakdown: II. Initial value problem

C.-Y. Kao, F. Brau, U. Ebert, L. Schäfer, S. Tanveer
Department of Mathematics,
The Ohio State University, OH 43210, USA,
Centrum Wiskunde & Informatica (CWI), P.O.Box 94079,
1090GB Amsterdam, The Netherlands,
Groupe de Physique Nucleaire Théorique,
Université de Mons-Hainaut,
Académie universitaire Wallonie-Bruxelles,
Place du Parc 20, 7000 Mons, Belgium
Department of Physics, Eindhoven University of Technology,
Eindhoven, The Netherlands,
Fachbereich Physik, Universität Duisburg-Essen,
Lotharstr. 1, 47048 Duisburg, Germany.
July 12, 2019

An interfacial approximation of the streamer stage in the evolution of sparks and lightning can be formulated as a Laplacian growth model regularized by a ’kinetic undercooling’ boundary condition. Using this model we study both the linearized and the full nonlinear evolution of small perturbations of a uniformly translating circle. Within the linear approximation analytical and numerical results show that perturbations are advected to the back of the circle, where they decay. An initially analytic interface stays analytic for all finite times, but singularities from outside the physical region approach the interface for , which results in some anomalous relaxation at the back of the circle. For the nonlinear evolution numerical results indicate that the circle is the asymptotic attractor for small perturbations, but larger perturbations may lead to branching. We also present results for more general initial shapes, which demonstrate that regularization by kinetic undercooling cannot guarantee smooth interfaces globally in time.
PACS: 47.54.-r
Keywords: moving boundary, kinetic undercooling regularization, initial value problem, Laplacian instability, electric breakdown

1 Introduction

Propagating fronts in Laplacian growth occur naturally in quite a number of physical problems including viscous fingering [1, 2, 3, 4, 5], electro-chemical growth, dendritic crystal growth for small undercooling [6, 7, 8], and void migration in a conductor [9, 10, 11]. More recently, it has been shown that this class of problems includes the ’streamer’ stage of electric breakdown [12, 13, 14, 15, 16, 17, 18, 19], which will be described below. A central issue in these problems is the stability of curved fronts. In a limiting case, most of these models reduce to the classic Saffman-Taylor problem [1], which is known to be ill-posed [20, 21]. Numerical as well as formal asymptotic results [4, 8, 5, 22] suggest that one branch of steadily propagating finger or bubble solutions in a Hele-Shaw cell is stabilized by surface tension regularization, though only recently some mathematically rigorous results [23, 24] are available to justify nonlinear stability to small disturbances in the special case of a nearly circular bubble. Besides surface tension, other regularizations [9, 10, 11, 25] have also been analyzed. In the present paper we study both the linear and the nonlinear initial value problem for one such regularization, in particular, the stability of a steadily propagating circular shape. This regularization is called kinetic undercooling in the crystal growth context111In crystal growth, kinetic undercooling is likely to be the more important regularization compared to the Gibbs-Thompson effect for large undercooling where the Laplacian growth model becomes questionable. This happens when the time scale over which the interface evolves becomes comparable to the time scale of heat diffusion., but has a different physical interpretation for streamers.

During the streamer stage of electric breakdown the discharge paves its way through a nonconducting medium, leaving behind a weakly ionized conducting channel. The basic growth mechanism is impact ionization due to electrons strongly accelerated in the local electric field. In a sufficiently strong field, a thin space charge layer forms around the head of the streamer. This layer screens the field in the inner ionized region to a very low level, and the growth of the streamer is driven by the electrons moving and multiplying in the strong self-enhanced field ahead of the curved ionization front.

For sufficiently strong external fields, the thickness of the electron layer is small compared to the radius of the streamer head. Therefore Lozansky and Firsov suggested (mainly in the Russian literature, but also in [26]) that this layer can be modeled as an interface separating the ionized from the non-ionized region. Probably, the idea is even older, since a similar concept was already proposed by Sämmer in the German literature in 1933 [27]. However, a deeper study of the implications of this concept started only later [12, 13, 14, 28, 29, 30] where the problem is placed in the context of other Laplacian moving boundary problems. The validity of the moving boundary approximation for negative streamers is discussed in [15] for simple gases like pure nitrogen or argon, and in [31] for air. The dimensional analysis and the proposed regularization mechanism of the moving boundary problem for negative streamer ionization fronts are discussed in detail in our previous papers [16, 17, 18, 19].

In dimensionless form, the model is defined as follows. The normal velocity of the interface is given by the drift velocity of the electrons, which is proportional to the local electro-static field . In appropriate units it takes the dimensionless form:


where the super-script denotes the limiting value as the interface is approached from the exterior (the non-ionized region) and is the outward normal on the interface. Outside the streamer the electric potential obeys the Laplace equation:


An analytical and numerical analysis of the underlying physical model formulated in terms of partial differential equations for charge densities and field suggests the interfacial condition




Far from the streamer, the electric field tends to a constant222 A correction of order to the electric field E can occur only if the streamer carries a net electric charge. We here concentrate on the analysis of streamers that are globally electrically neutral.


where is the unit vector in -direction. Eqs. (1)–(5) define our model.

In two dimensions, a simple solution to the free boundary problem posed by this model takes the form of a uniformly translating circle. Our previous work in [17, 18, 19] and the present paper are primarily concerned with the linear and nonlinear stability of this solution to small perturbations. It is to be noted that the circular shape differs from an actual streamer shape. However, the front half of a circle roughly resembles the shape of the front part of a streamer. Since growth of disturbances is found to be most pronounced in this advancing part of the interface, we expect stability features found here to be qualitatively relevant for an actual streamer and more generally for curved fronts.

In the special case , the linearized evolution of small perturbations can be determined exactly in our model [17, 18]. The case of general is treated in part I [19] of this series of papers and in the present manuscript. In [19], we discussed the spectrum of the linear operator which results from the linear stability analysis of the circular solution. Restricting ourselves to an appropriate space of analytic perturbations we found a pure point spectrum. Asymptotically in time, except for the trivial translation mode, all eigenmodes were found to decay exponentially in time. These eigenmodes are singular at the back of the bubble; nonetheless, as evidenced in the present paper, this singularity is not reflected by the actual linear evolution near the bubble back. The usual asymptotic form of the solution for large time: , where is the eigenfunction corresponding to the eigenvalue , fails in a neighborhood of the rear of the bubble, though it holds elsewhere.

In the present paper, we consider the initial value problem. For the linearized evolution, analytical results are obtained in the limit of strongly localized disturbances of the circle. Also the large time behavior of general perturbations can be studied analytically. Numerical calculations confirm these results. Together with the eigenvalue analysis of the first paper [19], clear evidence of linear stability is presented. The full nonlinear evolution of a perturbed interface is calculated numerically. Our results suggest that, similar to linear evolution, small enough perturbations of a circular bubble grow in the front part of the bubble, but eventually decay as interfacial distortions advect to the bubble rear. Nonetheless, when is small but nonzero, the large transients in the linear regime make nonlinearity important even when the initial perturbation is exponentially small in . Furthermore, when the perturbations are larger, the circle is no longer an attractor of the dynamics and the propagating structure branches. For general initial shapes, we give some numerical evidence that the undercooling regularization condition can not guarantee a smooth interface globally in time. For some initial conditions, the interface tends to develop a sharp corner in the back. Other initial conditions lead to the separation of the moving body into two parts.

This paper is organized as follows. In section 2 we present equations derived earlier in a conformal map setting. Section 3 is devoted to the linear evolution of perturbations of the circle. Subsection 3A recalls previous results, and in subsection 3B we present rigorous results on the growth of a strongly localized perturbation. We continue the discussion of localized perturbations in subsection 3C and explain at an intuitive level how strongly localized perturbations are generically advected to the rear of the circle, increasing in amplitude in the front-half before decreasing in the back half. Mathematically, the advection is described by a one-parameter family of conformal maps which is a subgroup of the automorphisms of the unit disk. The important role of this subgroup has been previously established for the exactly solvable case [17, 18]. In subsection 3D we discuss the anomalous behavior found at the back of the circle in the large time limit. In subsection 3E we give arguments indicating that an initially analytic interface stays analytic for all finite times, but singularities initially outside the physical region of interest approach the back of the circle for . Provided the perturbation for stays analytic in the closed unit disk, except for the point , we in subsection 3F prove that it asymptotically reduces to a constant. This implies that the perturbation just leads to a shift in space with respect to the unperturbed propagating circle. In subsection 3G, we present numerical solutions of the linear evolution equations. These calculations support the asymptotic results derived in the previous subsections. For disturbances, not necessarily localized, we present evidence that on any part of the interface not containing a neighborhood of the bubble rear, the decay rate of the disturbance matches what is expected from the prior spectral analysis [19].

Section 4 presents a numerical study of the nonlinear evolution for different perturbations. We first consider perturbations of a circular bubble. It is shown that the circular bubble can be nonlinearly stable if the perturbation is small. However, when the perturbation is large enough, the front may start to branch. Furthermore, we study the nonlinear evolution for more general initial configurations. It is shown that the formation of a cusp precisely on the back side of the moving body can not be excluded. We also observe that the body might split into two parts.

2 Equations resulting from conformal mapping

As already explained repeatedly [16, 17, 18, 19], we assume the streamer to be a simply connected compact domain in the -plane. The area of is conserved under the dynamics and equals in dimensionless units. Identifying the -plane with the closed complex plane , we introduce a conformal map that maps the unit disk in the -plane to the complement of in the -plane


The Laplace equation (2) and the boundary condition (5) are incorporated in the definition of a complex potential .


Both functions and are analytic for . The physical potential is related to as


The remaining boundary conditions (1), (3) take the form


The problem reduces to solving these two equations, respecting the analyticity properties of and .

A simple solution corresponding to a steadily translating circle is given by


In physical space it describes a unit circle moving with constant velocity in -direction. For small and smooth distortions of this circle, it is appropriate to look for solutions of the form


where and are analytic in and is a small parameter. Since the area is conserved, it can be shown that the residue of the pole in (11) remains unchanged to first order in . Substituting (12) into equations (9), (10) we in first order in find a system of two partial differential equations, from which can be eliminated. The final equation for takes the form


where we introduced the rescaled time variable


Eqs. (13), (14) determine the linearized evolution that will be discussed in Section 3. We will assume that the initial interface is analytic, i.e., that all singularities of are outside the closed unit disk , though much of the analysis is valid for a sufficiently smooth interface as well333Analyticity is not crucial, except in §3E, 3F..

3 Analysis of infinitesimal perturbations

3a Summary of previous results

In part I [19] we have analyzed the eigenvalue problem, resulting from Eqs. (13), (14) via the ansatz . We have shown that the spectrum is purely discrete and that the real part of all eigenvalues is negative, except for the trivial value , which corresponds to a simple shift of the circle. An infinite set of real negative eigenvalues was found. All eigenfunctions, except for , are singular at at the back of the circle. Thus the expansion of a regular initial condition in terms of eigenfunctions has to break down in the neighborhood of , which indicates that in that neighborhood some anomalous relaxation shows up. Furthermore, we found that as , any eigenvalue tends to zero and the corresponding eigenvector tends to a constant. A similar behavior of the spectrum was found for a steadily moving circle in a Hele-Shaw cell with surface tension regularization [22] and this degeneracy is not unexpected since the unregularized problem is mathematically ill-posed [20, 21].

Here, we consider the initial value problem defined by Eqs. (13), (14). Our analysis is guided by previous results [17, 18] on the special case where the general time dependent solution is known analytically; it is


where the function is given by the initial condition,


and is defined as


The properties of these solutions are discussed and visualized in detail in [17, 18]. Here we in particular note that the essential time dependence of is contained in the transformation


, , defines a one-parameter family of automorphisms of the unit disk, with fixed points . The point is stable, whereas is unstable in the following sense: as , i.e. , all the complex -plane, except for , is mapped into a neighborhood of . This results in an advective dynamics. Any perturbation not centered precisely at is advected towards , where it vanishes asymptotically. As , only a shift of the circle is left:


However, it is to be noted that the limit is not uniform, and no matter how large is, there is a neighborhood of , where may change dramatically. We note that advection of distortions from the front to the sides has been observed in viscous fingering and crystal growth models with surface tension and has been derived from somewhat heuristically simplified models [4, 32]. We further note that in the limit a purely advective dynamics results [18]:

Expecting the automorphism and the resulting advective dynamics to play an important role also for we transform the PDE (13), (14) from variables to variables , introducing the notation


This results in the normal form of a hyperbolic PDE:




3b Localized perturbations; rigorous results

Consider for general an initial perturbation that is centered at for and has ‘width’ in the sense that decays rapidly with when . The decay rate will be specified more precisely below Eq. (30). To study this problem, we first write (22) as an integral equation:


where and denote derivatives of . Then integration by parts in s replaces by which is written as

Here is a reference point in the tail of the perturbation chosen such that . We assume that is so large that is negligible.

Then, after some algebraic manipulation, we are able to rewrite (24) as the following equation for


where, with the understanding that , ,




With considered known444Since cannot be determined without considering the full non-local problem on , part of the expression (29) for is not known. Nonetheless, if a disturbance is localized, the contribution to from will be relatively small. In any case, in order to study the evolution in the -scale, we are not prevented from considering as known. we determine the solution to the integral equation (26) for , , where and are some suitably chosen positive values independent of . Now it is clear from the expression for and that they are uniformly small in the norm when is sufficiently small. We now choose the norm


where the positive weight function obeys

For example, for and 1 for would suffice for our analysis. We define to be localized if is finite, and , can be chosen such that is negligibly small for . Now it is clear from (26) that the linear operator has the contractive property


It follows that there exists a unique solution to the integral equation (26) if is small enough and that for

provided and are in the above specified range. For a perturbation localized in the sense given above, our result reduces to


We note that in general will not vanish for . This is the reason for restricting to the interval given above and indicates that for the localized perturbation will sit on top of a dynamically generated delocalized background of amplitude .

A detailed discussion of the result (32) will be presented in the next subsection.

3c Localized perturbations; formal intuitive arguments

It is useful to obtain the result (32) through a more formal, yet intuitive, reasoning. This will also be helpful in our subsequent treatment of the long-time asymptotics in the anomalous region near the back of the bubble. We again restrict the analysis to the unit circle , or correspondingly to . According to Eq. (19), the two angular coordinates and are related through


Initially, (at ), and obviously are identical. In terms of , the PDE (22) takes the form




We now search for a solution that during its evolution stays localized near a fixed angle , with an angular width . We use the ansatz


where again


and is assumed to vanish rapidly for . With this ansatz, Eq. (34) takes the form


For we neglect the term and the -dependence in the argument of to find an approximate solution of the form


which is the same as (32).

Before we evaluate this result we briefly discuss its limitations, as resulting from the present derivation. In view of the assumptions and , the use of the zero order result is justified provided


This is valid for all times provided , i.e., for initial conditions which essentially vanish in the forward direction . For the condition (40) is violated if becomes of the order , and therefore the approximation becomes invalid in the large-time limit . This special role of perturbations in the forward direction is not unexpected since for such perturbations advection is ineffective.

Neglecting the term has more serious consequences. Substituting into Eq. (3C) an ansatz of the form

one finds that the result for violates the condition for . A localized initial condition dynamically generates a delocalized contribution, with an amplitude proportional to , in full accord with the rigorous discussion of the previous subsection. Again this result is not unexpected since the eigenfunctions of the operator , Eq. (14), are delocalized. Assuming that we can expand an initially localized perturbation in terms of eigenfunctions we must expect that the balance of the expansion coefficients , which for leads to localization, is destroyed by the time evolution. With these limitations in mind, we now discuss the result (39).

According to Eq. (39), if expressed in the variable the evolution of the perturbation is most simple. Neither its position nor its shape change. Only the overall amplitude varies with time. For , i.e., if is at the front half of the circle, increases up to a time given by


and then decreases again. For decreases monotonically. For any , we find the asymptotic behavior


For a perturbation centered precisely at the back of the circle , exponential relaxation


holds for all . We recall that the localized approximation must break down it becomes of the order of the amplitude of the delocalized background. Nevertheless we will argue in subsection 3D that a contribution with asymptotic time behavior generally shows up.

Using Eq. (33) to transform back to we see that the center is convected along the circle, reaching for . A little calculation yields the velocity of this advection


This result has a simple interpretation. Recalling that we are working in a frame moving with the velocity of the unperturbed circle, we identify the velocity (44) as the projection of onto the tangent to the circle at the instantaneous location of the perturbation.

In terms of the overall amplitude of the perturbation takes the simple form


It increases as long as the perturbation is on the front half of the circle and decreases on the backside. The maximum, reached for , strongly depends on the initial position .

Defining the scale factor of the width of the perturbation as


we find


Thus the width behaves similarly to the amplitude, except that for it varies much less. For it vanishes like .

So far we considered perturbations localized away from the tip of the circle. For , Eq. (39) still holds for times such that

cf. Eq. (40). It describes the initial increase and broadening of the perturbation. Advection, of course, is absent. For the width becomes of order and the local approximation clearly becomes invalid.

On the qualitative level these results are most similar to the exact results found for [17, 18] and resemble the dynamics of a localized perturbation found in the context of viscous fingering [32].

Figure 1: Evolution of a strongly localized perturbation for . Curves 1, 2, 3 correspond to times , 1.84 and 4.59, respectively. a) Re  as function of the comoving angular coordinate . Broken lines: local approximation. Full lines: exact result with the shift of the circle subtracted. b) Perturbed interface in the physical plane in the system of local tangential and normal coordinates as explained in the text.

The quantitative performance of the local approximation is illustrated in Fig. 1, where for the exact evolution of a localized perturbation is compared to our approximation. From the exact result (16) the contribution representing a simple shift of the circle, has been subtracted. The initial condition is chosen as


Fig. 1a shows as function of for three different times. Curve 1 shows the initial condition, where by construction the exact form and the approximation coincide. Curve 2 shows the perturbation when it is largest, in -space being located near . Curve 3 is taken at some later time. Evidently in this example the local approximation, (broken lines), is quite accurate. Very similar results are found for , which therefore is not shown. Fig. 1b shows the effect of this perturbation in physical space. In evaluating , Eq. (12), we choose the amplitude . To combine the three curves into one plot, we introduced a time-dependent rotation of the coordinate system such that or are measured along the normal or the tangent to the unperturbed circle at the center of the perturbation, (i.e., at angle , since the inversion contained in the conformal map induces a sign change of the angles). In this representation the exact solution and the approximation cannot be distinguished within the resolution of the plot. We note that in physical space the shape of the perturbation varies due to interference with the unperturbed circle.

3d Asymptotic relaxation near

In discussing the asymptotic relaxation we prefer to rewrite (22) in terms of , using . Inserting the explicit form (23) of and multiplying by we find


Here stands for

cf. Eq. (18). Keeping only the leading -dependence in the coefficients of the derivatives, we reduce Eq. (3D) to


For we neglect the term to find


where and depend on the initial condition and of course cannot be fixed by this asymptotic argument.

Since the derivative in Eq. (49) is multiplied by , the neglect of the term involving can be justified only for . In terms of


this implies that we deal with a neighborhood of that is contracted to this point like . This range of is complementary to the region where an expansion in terms of eigenfunctions can be expected to be valid asymptotically.

In the result (50) the -dependence is suppressed by a factor , which for , , vanishes much faster than . Thus -dependent corrections of order will dominate the asymptotic relaxation at the back of the circle. Noting the presence of in the coefficients of the differential equation, it is natural to determine the structure of these terms with the ansatz




From Eq. (3D) with replaced by we find


where the are integration constants. Generally is found to be a polynomial in of degree . In this analysis we assumed . For the ansatz (53) has to be modified. In particular a term proportional to has to be included. We note that the exact result for shows such a contribution [18].

To transform our result back to -space we introduce


Eq. (19) yields


Thus has an expansion of the form


where the again depend on the initial condition. For the terms of order dominate over the contribution . For we therefore in a region of size near expect to see a very smooth asymptotic relaxation of the interface, with only a few coefficients depending on the initial condition. In contrast, for the asymptotic relaxation is determined by the term , which will depend on the initial condition in a complicated way. For this is illustrated in Fig. 5.2 of Ref. [18]. In the next subsection we will argue that the function picks up contributions due to singularities of the initial condition, which for are driven towards . We finally note that the results discussed here resemble the behavior of the low order eigenfunctions . As shown in part I [19] of this series, these functions near develop a singularity of the form , implying that the derivatives at exist for all orders .

To illustrate our results we consider a perturbation centered at . As initial condition we choose


and we calculate the function


We expect to find the limiting behavior




respectively. Whereas depends on the initial condition, the limit (61) is universal. The results shown in Figure 2 conform to these expectations.

Figure 2: Time evolution of the initial condition (59). a) as function of for times as given and for . b) Phase of for the same values of and . c) Again as function of , but now for .

Fig. 2a shows for several values of and for . It illustrates the approach to the limiting form , which within the accuracy of the plot is in fact reached for . Fig. 2b shows the corresponding phase of . Here the approach to the limit is slower, but is definitely visible. Fig. 2c shows results for , . Here seems to approach a limiting curve which clearly shows remainders of the initial peak. (We should note that is symmetric: , and that the peak at , of course, is rounded, which however is not visible on the scale of the plot). We finally recall that the -range shown here in terms of corresponds to a small region near . Specifically for it corresponds to .

For the asymptotic relaxation our results predict


This prediction is tested in Figure 3 by plotting results for as function of for several values of . The expected behavior is reasonably well observed.

Figure 3: as a function of for several values of . The initial condition is given in (59). The lines indicate the expected slope for or for , respectively.

3e Analyticity of the interface

If we assume the initial interface to be analytic, all singularities of the initial perturbation have to be outside the closed unit disk, . We here argue that under the linearized dynamics the singularities stay outside for all finite times . For they approach , and contribute to the anomalous behavior found in subsection 3D.

This argument is based on the recurrence relation for the coefficients in the Taylor expansion


The evolution equation (13), (14) yields


The singularities of are determined by the behavior of the in the limit . For simplicity, we consider an initial condition whose singularity closest to is a branch point at with behavior for nonintegral or a pole with a non-positive integer. Then for