Growing discharge trees with self-consistent charge transport: the collective dynamics of streamers

When a high electric voltage is suddenly applied to ionizable matter, electric breakdown frequently takes the form of growing filaments, and these filaments can form a complex tree structure. Discharge trees are observed in streamer coronas around tip or wire electrodes, in the streamer coronas ahead of propagating lightning leaders [1] and in the (hot) leaders themselves. Streamer discharge trees also appear in transient luminous events such as jets [2], gigantic jets [3] and sprites [4] between thunderclouds and the ionosphere. Streamer and leader trees are a generic response to high voltage pulses; they appear in various gases, liquids and solids in plasma and high voltage technology.

Our understanding of such non-thermal, filamentary electrical discharges is remarkably unbalanced. On the one hand, we are now reaching a very detailed knowledge of their microphysics; this includes models of electron energy distributions [5], and of transport coefficients and cross-sections of the main reactions, at least for air and other common gas compositions. This knowledge translates into sophisticated and reasonably accurate models of single streamers [6–11], the initiation of the streamer branching [5, 11–13], the merging of two nearby streamers [14, 15] and the influence of surrounding mesoscopic inhomogeneities [16, 17]. On the other hand, we barely understand most macroscopic processes in a fully developed corona or streamer tree involving hundreds or thousands of mutually interacting plasma filaments. The large scale transport of charge, the internal electric fields and the influence of the many surrounding streamers on a single streamer have rarely been discussed in the literature. However, these mechanisms are relevant for the propagation of long sparks [18–20] and the approach of lightning leaders toward protecting rods. The overall tree structure also determines which volume fraction of the medium is 'treated' by the discharge, creating radicals, ions and subsequent chemical products relevant for plasma technology and for the production of greenhouse gases during a thunderstorm.

Furthermore, as our results will show, the collective dynamics of a streamer tree exhibit some counter-intuitive features that cannot be directly derived from microscopic models but are nevertheless required to explain the observed behavior of streamers. For example, streamer reconnection [21–26] is naturally explained by the opposing charge polarity between inner and outer branches of a tree.

Most studies on the growth of electrical discharge trees descend from the dielectric breakdown model (DBM) that Niemeyer et al [27] proposed in 1984 to explain the fractal properties of some electrical discharges such as Lichtenberg figures that propagate over a dielectric surface. In their model, a discharge tree expands in discrete time-steps by the stochastic addition of new segments with a probability that depends on the local electric field.

We are not aware of many models of fully three-dimensional (3D) streamer trees not based on the DBM. Only Akyuz et al [28] modeled streamers as a tree of connected, perfectly conducting cylinders that propagate according to simple rules based on the value of the electric field surrounding the tips. The computations required to solve the electrostatic problem limited their simulations to small trees with less than 10 branches.

The original DBM as well as [28] assume that the channels in the tree are perfectly conducting, but there is strong experimental evidence that the electric potential decreases along a discharge channel.

The common approach to introduce a potential decay along a streamer channel and inside the streamer corona is to assume that the electric field inside a streamer has a fixed value, the so-called stability field. For example, in air at standard temperature and pressure the stability field of positive streamers is thought to be 4–5 kV cm⁻¹. A fixed stability field is used to model the streamer corona that precedes a leader in a long spark discharge [29–31] or the enormous streamer trees in sprite discharges high above thunderstorms [32].

However, the concept of a fixed field inside streamer channels lacks any theoretical support. Rather, it is based on a phenomenological interpretation of experiments that nevertheless have not measured the internal streamer fields. Originally, the concept of stability field referred to the minimum average applied field for sustained streamer propagation in a gap between parallel plates [33, 34] where the discharge was initiated from a protruding electrode. The existence of such a minimum field around 4–5 kV cm⁻¹ was interpreted [35, 36] in terms of a now discarded model of streamers as isolated propagating patches of charge. Later it was found that the relation between the applied potential at the originating electrode U and the longest streamer length L is roughly linear with U/L ≈ (4.5–5) kV cm⁻¹ in air [37]. By relating this observation to the existing concept of a stability field, the results were interpreted as indicating that the stability field was the electric field inside the streamer channel. However, even the earliest numerical simulations of two-dimensional (2D) streamers [38] already showed a clearly non-constant electric field in the channel. As we will see, this variation is enhanced by the collective dynamics of a streamer tree. Indeed, our results will show that the assumption of a constant electric field in all streamers is in contradiction with a consistent charge transport model, as long as conductivity stays finite.

Recent simulations of density models resolving the inner structure of streamers already have established the relevance of a self-consistent charge transport model for the dynamics of streamer channels and, in particular, for the dynamics of the electric field in the channel. For upper-atmospheric streamers, Liu [9] and Luque and Ebert [39] independently showed that the re-brightening of sprite streamer trails is due to a second wave associated with a significant increase of the electric field in the sprite channel; Luque and Gordillo-Vázquez [40] postulated later that sprite beads are also caused by persisting and localized electric fields. These electric fields may only persist due to a finite conductivity in the streamer channel [41], which also sets their decay times.

To our knowledge, the only DBM-inspired models that treat the charge transport self-consistently appear in the context of discharge trees in dielectrics [42], generated when a solid insulator is subjected to an intense, repetitive electrical stress [43].

In the present paper, we first outline the general structure of a model for growing discharge trees that consistently incorporates charge conservation. Then we introduce the simplest model for a streamer corona as a tree structure of linear channel segments with a finite fixed diameter and with a finite fixed conductivity. The streamer channel tips advance and branch according to simple, phenomenologically motivated rules. We analyze the internal electric fields and the transport of charge in fully branched, extensive streamer coronas. This is a stepping stone toward more realistic and detailed models and, although many improvements of our approach are straightforward, we have often kept complexity to a minimum in order to focus on the overall qualitative behavior of streamer trees with realistic conductivities and consistent charge transport, which appears to be largely unexplored in the existing literature.

The paper is organized as follows: in section 2 we give general prescriptions for discharge tree models with self-consistent charge transport, which are then particularized into the simplest streamer tree model, which we have implemented. We present the most relevant results of the model in section 3. Finally, section 4 concludes with a short summary and discussion.

Numerous experimental observations of discharge streamers and leaders show a structure of branching filaments [24, 44, 45]. The understanding that has evolved over decades since Raether's seminal work in the 1930s [46] is that streamers are able to penetrate into areas where the background field would be too low to maintain an ionization reaction; as they are conducting filaments they enhance the electric field at their tips to values above the breakdown value which allows them to grow there. While photography with nanosecond resolution shows these active streamer heads in air as glowing dots [47], simulations of single streamer channels (that are typically performed with 2D fluid models) [6, 9, 38, 48, 49] reveal the inner structure sketched above: long ionized filaments develop a thin surface charge layer around their whole ionized body and maintain in this manner a low field in their interior while enhancing it at their tips. The tip propagates with a velocity comparable to the electron drift velocity in the enhanced field or even faster, while the lateral surface charge layer for positive streamer filaments consists of much heavier positive ions and is depleted of electrons (while negative streamers also can develop some lateral dynamics as their surface charge consists of an electron overshoot, this can weaken the field enhancement at their tips during their evolution) [49, 50]. Therefore the growth at least of positive streamers can be modeled as the growth of a conducting channel at its tip only. While streamer propagation has a long history of experimental and theoretical investigations, streamer branching and streamer interaction are now being investigated as well by experiments [45, 51–53] and theory [12, 13, 54]. Leader dynamics, though investigated in less detail up to now, is believed to evolve in a similar manner through field enhancement at the channel tip—however, its conductivity is maintained over longer times by Ohmic heating and its trajectory is paved by a streamer corona.

As discussed in the introduction, to model the discharge tree as a growing network of conductors was already suggested by Niemeyer et al [27]. In this work, we concentrate on elaborating the charge conservation and the charge transport within the discharge tree. We remark that charge content and electric field distribution are typically experimentally not accessible, except when a streamer discharge propagates over a dielectric surface [55]. Therefore these features have to be derived theoretically or observed indirectly, e.g. through streamer reconnections. The geometric structure of the network with its charge content and the external electric field determine the actual electric field distribution; this field distribution together with the conductivity distribution within the network determines the consecutive charge transport in the tree, and the local field distribution at the tip determines growth and branching of the tree tips. The tip dynamics determines the diameter, conductivity and tree structure of the newly grown parts of the network.

Let us now discuss the general structure of such a model with reasonable approximations, before introducing the simplest manifestation of such a model in section 2.2, the numerical implementation in section 2.3, and the particular choice of model parameters for positive streamers in ambient air in section 2.4. The goal of the work is to overcome the limitation of current fluid models to model only a few filaments, and to implement our current microscopic understanding into a coarse grained tree model.

2.1.1. Linear channel parts: radius R, line charge density q, line conductivity σ and electric current I.

A schematic diagram of a linear channel part is provided in figure 1. We parameterize the channel length with a longitudinal or arc length coordinate s, and we assume these parts to be cylindrically symmetric with a radius R(s,t). In general, we can assume that the radius varies slowly over the arc length s. The electric charge typically resides in the surface of the channel. It can be assumed to be cylindrically symmetric as long as other charges stay at a distance much larger than the channel radius. According to standard electrodynamics, the electric field created by the charge of the channel is determined only by the line charge density and not by the channel radius at distances much larger than the channel radius.

Zoom In Zoom Out Reset image size

The conductivity of the channel is provided by the densities n_e,± and the mobilities μ_e,± of the electrons and of the positive and negative ions inside the channel; these densities of charged species have been created at the streamer tip. The current flowing through a cross-section of the channel is

$$\begin{equation} I(s, t) = \int 2\pi r\,{\rm d}r \;E_s(r, t)\left(\mu_{\mathrm{e}}n_{\mathrm{e}} + \mu_+n_+ + \mu_-n_-\right)(r,s,t), \end{equation} \tag{ 1 }$$

where E_s is the longitudinal component of the electric field along s. Inside the space charge layer the electric field does not essentially change in the radial direction, and it is oriented along the channel [56]—otherwise the current would flow into or out of the channel walls and would change the charge content very rapidly; hence as long as charges change slowly, the field is directed along the axis and E_s = E.

Therefore we can rewrite (1) as Ohm's law,

$$\begin{equation} I(s,t) = \sigma(s,t)\;E(s,t), \end{equation} \tag{ 2 }$$

where we have defined a line conductivity σ(s,t) as the integral of the conductivity over the channel cross-section,

$$\begin{equation} \sigma(s,t) = \int 2\pi\;r\,{\rm d}r\;\left(\mu_{\mathrm{e}}n_{\mathrm{e}} + \mu_+n_+ + \mu_-n_-\right)(r,s,t). \end{equation} \tag{ 3 }$$

We define the line charge density q(s,t) by the integral of the charge density over the channel cross-section

$$\begin{equation} q(s,t) = \int 2\pi\;r\,{\rm d}r\;{e}\;\left(n_+ - n_{\mathrm{e}} - n_-\right)(r,s,t), \end{equation} \tag{ 4 }$$

where e is the elementary charge. The line conductivity is the inverse of the resistance per length, and the line charge density is the charge per length.

The conservation of electrical charge implies that

$$\begin{equation} \partial_t q(s,t) + \partial_s I(s,t)=0. \end{equation} \tag{ 5 }$$

For radius R(s,t) or line conductivity σ(s,t) particular dynamical equations could be implemented that incorporate a physical understanding of the channel dynamics. Alternatively they can be considered as fixed after they have been generated by the motion of the channel head.

2.1.2. Head radius, charge, velocity and branching.

The charge distribution in the discharge head and channel together with the external field determine the electric field distribution at the head. The head velocity in general depends not only on the electric field in some particular spot but also on the electric field E_enh and electron density distribution in the whole ionization region at the discharge head; and the shape of this region is strongly determined by the head radius R. The velocity of the head or tip can therefore be considered as a function of radius R, electric field E_enh, polarity ± and of the gas type and conditions,

$$\begin{equation} {\bf v}_{\rm tip}^\pm = {\bf v}^\pm(\mathbf{E}_{\rm enh},~R,\mbox{ gas type and pre-ionization}). \end{equation} \tag{ 6 }$$

For the velocity of streamers in air, Naidis [57] has suggested a particular analytic approximation.

For branching of the channel tip, an appropriate distribution as a function of the head parameters has to be found. For positive streamers in air, both experimental [45, 51–53] and theoretical [12] studies have been presented; they constitute the start of quantitative investigations.

The channel conductivity is also created at the channel tip. Particular results for ionization degrees for streamers in air will be discussed later. For leaders, also a reduced medium density due to thermal expansion contributes to increasing the electrical conductivity of the channel.

2.1.3. Electric field.

The electric field is given by the external field plus contributions due to the charges in the tree. In the density approximation, the electric potential is given by the classical equation

$$\begin{equation} \phi({\bf r}) = \phi_{\rm ext}({\bf r}) + \frac1{4\pi\epsilon_0}\int {\rm d}{\bf r}'\; \frac{{e}\,(n_+-n_{\mathrm{e}}-n_-)({\bf r}')}{|{\bf r}-{\bf r}'|}. \end{equation} \tag{ 7 }$$

We recall that the electrical charge density e (n₊ − n_e − n₋) is non-vanishing essentially only in the walls of the channels, at the radius R. When approximating the channel by a line as above, the kernel in (7) has to be modified by a regularization to avoid unphysical singularities for | r −  r'| → 0. We use

$$\begin{equation} \phi({\bf r}) = \phi_{\rm ext}({\bf r}) + \frac1{4\pi\epsilon_0}\int {\rm d}s\; \frac{q(s)}{|{\bf r}-{\bf r}(s)|+R}. \end{equation} \tag{ 8 }$$

We selected this regularization after tests with some other kernels that lead to oscillatory behavior of charges. Further investigations are under way.

2.1.4. The general setup of this model.

In the current paper, we will make a number of assumptions to make the model as simple as possible. This will allow us to identify the key new features induced by consistent charge transport, without having to wonder whether properties are due to certain other model features.

In this simplest model, we assume that all channel parts and tips have the same time-independent radius R and line conductivity σ. This amounts to considering the ion density as fixed after the initial ionization wave, which is justified by the low ion mobility and by the relatively long time scales of chemical processes such as attachment that would otherwise affect the ion density. We also assume that the streamer head velocity is proportional to the local electric field and that branching is a Poisson process depending on the length of the streamer segment.

Together with the electric potential being fixed at the boundary of the simulation domain, and with the location of the electrode that supplies the electric current, these assumptions characterize the physical model.

We shall describe now the numerical implementation of the model described above. This numerical implementation, along with all the input files used in this article, is freely available⁵.

As sketched in figure 2, we replace the continuous arc lengths s of the different linear channel parts by the set i = 1,...,N of N charged nodes at positions r_i, each containing a time-dependent charge q_i(t), and a time dependent electric potential ϕ_i(t) is attributed to each node. The tree evolves through two coupled mechanisms. Firstly, due to the electric field, charge is transported along the edges. Secondly, each channel grows or branches at its tip according to the local conditions. In our model, we alternate between these two evolutions: to evolve our system from time t to time t + Δt we first calculate the electric field and transport the charge in the tree for an interval Δt, and then we add new nodes at the tips of existing channels, allowing some channels to branch eventually. The choice of the numerical time step Δt is discussed in appendix A. We describe now the steps of the simulation.

Zoom In Zoom Out Reset image size

2.3.1. Electric field with boundary conditions.

We assume that the stem of the discharge tree is connected to an upper planar electrode located at z = 0 that creates a constant background electric field E₀. This electrode together with the set of charges q_i with i = 1...N within the discharge tree creates an electric potential

$$\begin{equation} \phi_j = \frac{1}{4\pi\epsilon_0}\sum_{i=-N}^N{\frac{q_i}{\ell_{ij} + R}} + \phi_{\rm ext}(\mathbf{r}_j),\quad \ell_{ij}= |{\bf r}_i-{\bf r}_j|, \quad \phi_{\rm ext}(\mathbf{r}) = - \mathbf{E}_0 \cdot \mathbf{r} \end{equation} \tag{ 9 }$$

at node j, according to equation (4). Here i = −1... − N parameterizes the mirror charges introduced to keep the electrode at potential zero: for each charge q_i located at r_i = (x,y,z) a mirror charge q_−i = −q_i is located at r_−i = (x,y,−z). The node i = 0 is taken as the root of the tree; it is located at the origin, and it is discharged by contact with the electrode. Therefore q₀ = 0.

For details of the numerical solution of the electrostatic problem (9), refer to appendix A.

2.3.2. Charge transport within the tree.

During the relaxation phase, electric currents flow along the conductor links according to Ohm's law, where the current through each link is calculated from the potential difference between its two endpoints as

$$\begin{equation} I_{ij} = \sigma E_{ij},\quad E_{ij} = -\frac{\phi_i - \phi_j}{\ell_{ij}}. \end{equation} \tag{ 10 }$$

Due to these currents, the charge at node i changes as

$$\begin{equation} \frac{{\rm d}q_i}{{\rm d}t} = \sum_{j\in{\mathrm{neigh}(i)}} I_{ij}, \end{equation} \tag{ 11 }$$

where neigh(i) stands for the set of nodes connected to i. For the root node i = 0, q₀ = 0 is maintained because the current I₀₁ is exactly balanced by the current drawn from the electrode.

At each time step, we integrate the set of ordinary differential equations and (11), coupled with (9), from t to t + Δt. In our implementation, we used the real-valued variable-coefficient ordinary differential equation solver [58].

2.3.3. Growth of tree tips.

Each streamer in the tree grows at its tip, and we model this growth by adding a new node T' ahead of the old terminal node T after time Δt at the location

$$\begin{equation} \mathbf{r}_{T'} = \mathbf{r}_T + \mathbf{v}_T \Delta t, \end{equation} \tag{ 12 }$$

see figure 2(c).

The tip velocity v_T depends on the electric field distribution around the terminal node T. We approximate this distribution by the electric field in the node T generated by the background field and the charges of all other nodes plus the term F_T

$$\begin{equation} \mathbf{E}_T = \mathbf{E}_0 + \frac{1}{4\pi\epsilon_0}\sum_{\substack{j=-N, j\neq T}}^N{\frac{q_j \mathbf{e}_{jT}} {(|\mathbf{r}_j - \mathbf{r}_T|+ R)^2}} + \mathbf{F}_T, \end{equation} \tag{ 13 }$$

where e_jT is a unit vector pointing from r_j to r_T .

The term F_T accounts for the contribution of the terminal node T. In the limit Δt → 0, as the separation between nodes decreases, the charge contained in the terminal node becomes negligible compared with the many charges in the channel at distances shorter than R. For finite Δt, the term F_T accounts for the contribution of these many charges that are now summed up into the terminal charge q_T :

$$\begin{equation} \mathbf{F}_T = \frac{q_T \mathbf{e}_{\mathcal{P}(T)T}}{4 \pi \epsilon_0 R^2}, \end{equation} \tag{ 14 }$$

where $\mathbf {e}_{\mathcal {P}(T)T}$ is the unit vector that points toward T from its predecessor $\mathcal {P}(T)$.

Note that the terminal charge q_T is calculated in the same manner as all other charges q_i in the tree. The terminal node is added to the tree with q_T  = 0 and then it is charged through the newly created conducting link according to (10) and (11).

We will assume that the tip velocity is proportional to E_T through a model parameter that we name head mobility, μ_H. Since the charges that enter into equations (13) and (14) change continuously during the time interval Δt, we advance the streamer tips with a velocity that is linearly interpolated from its values at t and at t + Δt:

$$\begin{equation} \mathbf{v}_T = \frac{1}{2} \mu_{\rm H} \left[\mathbf{E}_T(t) + \mathbf{E}_T\left(t + \Delta t\right) \right]. \end{equation} \tag{ 15 }$$

Assuming a linear dependence of the tip velocity on the electric field is a strong simplification that nevertheless can be easily removed to incorporate more realistic dependences. In appendix B we study one of them, where we impose a minimum electric field for the streamer propagation.

2.3.4. Branching.

We model streamer branching only phenomenologically. The overwhelming majority of streamer observations show branching into two descendant branches—for rare examples of branching into three channels, see [52]. Therefore in our implementation we only considered binary trees.

Currently, microscopic models shed some light on the mechanisms of branching but are not mature enough to provide the quantitative predictions that our model requires. Therefore we adopted the simplifying assumption that branching is a Poisson process characterized by the length ℓ_branch along the streamer channel. Hence the probability that the streamer tip at T branches during a time step Δt is

$$\begin{equation} p = v_T\; \Delta t/\ell_{\mathrm{branch}}. \end{equation} \tag{ 16 }$$

We always ensure that the time step Δt is such that p ≪ 1. Once the algorithm has decided that a tip branches, the locations of its two descendant nodes are calculated as shown in figure 2(d); the locations of the two new tips r_T'± are symmetrical with respect to the location of the straight path (12):

$$\begin{equation} \mathbf{r}_{T'\pm} = \mathbf{r}_T + \mathbf{v}_T \Delta t \pm \mathbf{\delta r}, \end{equation} \tag{ 17 }$$

where δr is a random vector in the plane perpendicular to v_T with a bi-dimensional Gaussian distribution with standard deviation ℓ_sib.

Our model contains five dimensional parameters, the radius R of the discharge channel, the mobility μ_H of the channel head, the line conductivity σ, the average channel length ℓ_branch between two branching points and the initial separation ℓ_sib between two new branches. These parameters have to be chosen appropriately for the system under consideration, like streamers or leaders in different gases and at different pressures and temperatures.

For positive streamers in air at standard temperature and pressure we now estimate their values from phenomenological observations. These values are listed in table 1.

Streamer radius R. Depending on the applied voltage, visible streamer diameters in air at standard temperature and pressure vary between a minimum of  ≈ 0.12 mm [59] and 3 mm in the experiments of Briels et al [51] for sharply pulsed voltages of up to 100 kV, and increase up to the order of 1 cm in the experiments of Kochkin et al [20] with a Marx generator delivering pulses with voltages in the range of 1 MV. Due to the projection of the radiation into the 2D image plane and the non-homogeneous excitation of emitting species in the streamer head, the radiative or visible diameter is about half the electrodynamic diameter that parameterizes the extension of the space charge layer around the streamer tip, i.e. the visible diameter approximates the electrodynamic radius. Numerical simulations [7, 49] show radii in the range of 0.1–1 mm, similarly to the measurements of [51]. As streamers of minimal diameter generically do not branch, we have here chosen an electrodynamic radius of R ≈ 1 mm.

Head mobility μ_H. It was found in experiments [51] as well as in simulations [49] that the velocity of a positive streamer strongly depends on its radius. The analysis of Naidis [57] showed that the velocity of a uniformly translating streamer also depends on the peak electric field. This is because the peak field together with the radius determine the size of the region around the streamer head where the electric field is above the breakdown value and where the ionization grows. Naidis' numerical data for a fixed radiative diameter of 1 mm suggest a roughly linear approximation v ≈ μ_HE_p, μ_H ≈ 900  cm² V⁻¹ s⁻¹, where E_p is the peak electric field at the streamer ionization front.

Line conductivity σ. The electrical conductivity inside a streamer channel is dominated by free electrons. Most numerical simulations [7, 38, 60–63] agree on a value of about n₀ ≈ 10¹⁴ cm⁻³ electrons on the streamer axis, and a further analysis of the relation between peak field E_p and ionization density n₀ behind the front can be found in [64]. If we assume a quadratic decay of the density away from the axis up to a radius R, we obtain

$$\begin{equation} \sigma = 2 \pi {e} \mu n_0 \int_0^R r \left(1 - \frac{r^2}{R^2}\right) \, {\rm d}r = \frac{\pi}{2}e\mu n_0 R^2, \end{equation} \tag{ 18 }$$

where e is the elementary charge and μ ≈ 380  cm² V⁻¹ s⁻¹ is the electron mobility [65]. The expression (18) yields σ = 9.6 × 10⁻⁷ cm Ω⁻¹.

Branching ratio ℓ_branch/R. Briels et al [66] measured an approximately linear relationship between the average branching distance and the streamer radius for positive streamers in air. We use their value ℓ_branch/R ≈ 10, where R is the electrodynamic streamer radius.

Initial separation ℓ_sib between sibling branches. Finally, we used the arbitrary value 0.1R for ℓ_sib. The only constraints on this value are that it is much smaller than ℓ_branch and that it is of the order of vΔt, where v is a typical streamer velocity. Below, we will find that the effect of the value of ℓ_sib on the simulations is quite weak.

3.1.1. Simulation and overall structure.

Figure 3 shows a streamer tree simulated with the parameters of table 1 and an external electric field E₀ = 15 kV cm⁻¹ pointing downwards. This field corresponds to about half of the classical breakdown field. We colored the edge between two connected nodes i and j according to the mean electric field in the link, defined as

$$\begin{equation} E_{ij} = \frac{\phi_i - \phi_j}{\ell_{ij}}. \end{equation} \tag{ 19 }$$

We chose the order of the labels i and j such that the electric field is positive in the direction of streamer propagation.

Zoom In Zoom Out Reset image size

The xz projection of the streamer tree in figure 3 (upper left) has an approximately diamond shape; in the upper part the tree becomes wider at lower altitude due to the repulsion between the heads whereas in the lower part the tree gets thinner because the branches close to the center propagate faster. The diamond shape is typical in sprites [23] and in laboratory streamers [20, 47] captured before they contact the lower electrode. In needle-plane discharges, the strong divergence of the electric field around the needle electrode produces a sharper widening of the tree during the initial stages of evolution, hence in the upper part of the discharge.

We name the discharge structure in figure 3 a collective streamer front; it can be interpreted as a 'streamer of streamers'. The many positive charges at the tips of the lower channels have a role akin to the continuous space charge layer in a single streamer. Below them, they enhance the field around the center axis; above, the field is screened. In a single streamer, the charge is transported to the boundary due to the enhanced conductivity of the streamer channel; in a streamer tree, there is a coarse-grained conductivity arising from the many conductive filaments inside the tree. Figure 4 illustrates this phenomenon by plotting the electrostatic potential and electric field in a region around the streamer tree. The equipotential lines are further apart inside the tree, indicating a lower electric field, whereas they are compressed in the volume directly in front of the tree, where the electric field is significantly enhanced. Note that the electric field plotted in figure 4 reaches higher values than the internal electric fields in figure 3; this reveals the enhancement of the electric field close to the streamer heads but outside the channels.

Zoom In Zoom Out Reset image size

3.1.2. Non-constant electric fields inside the streamer.

3.1.3. Field reversal.

A salient feature of the fields shown in figure 3 is that in some channels the fields have opposite sign, transporting charge backwards. Although seemingly paradoxical, this results from some streamers outrunning others, as outlined in figure 5. The charges in a streamer create a field E_ch that opposes the external field E₀. Normally E_ch and E₀ add up to an internal field weaker than the external field but with the same orientation. Suppose, however, that the streamer is overrun by a few neighboring streamers carrying charges that screen E₀ inside the original streamer. Then only E_ch remains inside the channel, which thus starts to discharge. In that case the streamer halts, leaving a 'dead' channel behind.

Zoom In Zoom Out Reset image size

However, our algorithm, as described in section 2.3.3. adds new nodes to the tree tips even for very small values of the velocity defined in (15). The resulting slow growth of these dead channels is most often irrelevant for the overall dynamics of the streamer tree but may result in unphysical behavior, such as streamer channels slowly turning backwards.

This problem is solved by a field–velocity relation more realistic than the linear one in (15). In appendix B we discuss the inclusion of a realistic threshold electric field for streamer propagation.

The distribution of charges in the same simulation as in figure 3 appears in figure 6. To focus on the charge density inside the streamer channels, we have truncated the color scale, which would be otherwise dominated by the charges at the streamer heads.

Zoom In Zoom Out Reset image size

Figure 6 shows that while the lower part of the tree closer to the streamer tips is charged positively, the innermost segments are negatively charged. This resembles the negative charging of the upper regions of sprite streamers [39] and arises from an analogous mechanism. The many channels in the external branches transport a large amount of charge. The fewer channels in the inner sections collect this charge that then gets stuck due to the lower collective conductivity. Hence it brings about a negatively charged inner core in the tree.

In figure 7 we plot the current entering the streamer tree of figure 3 as well as the total net charge content of the tree as functions of time. The simulation reaches only up to 90 ns, the time of the first reconnection between streamer channels (see section 3.6.) and up to that point the current is increasing at an accelerating rate. For longer times the entering current is limited by the lost of conductivity in the upper channels due to electron attachment, presently not included in our model. It is therefore difficult to establish a direct comparison with empirical estimations of total charge, such as those by Ortega et al [67]. However, our estimation of a few amperes is similar to the peak current reported in [67]. Our simulation would require some hundreds of nanoseconds to reach the thousands of nanocoulombs measured in that experiment.

Zoom In Zoom Out Reset image size

We turn now to the influence of the line conductivity σ of the streamer channel on the propagation and shape of the streamer tree. We focus on this parameter because a straightforward dimensional analysis (see appendix C) shows that changing the line conductivity while keeping a fixed applied electric field is equivalent, after rescaling time, to a change in the external electric field with a fixed line conductivity. Therefore the analysis described here translates directly into a study of the influence of the applied field.

3.4.1. Branching angles.

At this point, it is helpful to suppress the randomness of the model and focus on an even simpler system. We run simulations where we impose a single branching point at z = −1 cm. In each of these simulations, we multiplied by a factor from 10⁻² to 10² the line conductivity discussed above and listed in table 1, here denoted σ₀. Figure 8 shows the results.

Zoom In Zoom Out Reset image size

The left panel of figure 8 shows the influence of the line conductivity on branching angles. Channels with a higher conductivity lead to wider branching. The reason is that charge moves more easily along the channel and then accumulates faster at the streamer tips. Equivalently, one can say that more electrons move upwards, leaving a higher positive charge in the tips. The electrostatic repulsion between both heads is thus stronger and they diverge more sharply.

However, figure 8 shows that this mechanism is quite weak. Although it is theoretically possible to infer the channel conductivities from branching angle measurements, such as those by Nijdam et al [45], the dependence seems too weak to be useful, given the natural variation and the measurement uncertainties of branching angles. In figure 8 we mark with arrows the branch-to-branch angles 30° and 50° from the branching point to underline that all conductivities agree with the branching angles of (39.7 ± 13.2)° reported in [45] for positive streamers in air at atmospheric pressure.

3.4.2. Velocity.

As we mentioned above, ℓ_sib does not substantially influence the simulations as long as it stays within reasonable physical bounds. To investigate this, we run simulations where we changed ℓ_sib from one tenth to twice the value in table 1. As in the previous section, in these simulations we forced the streamers to branch uniquely at a prescribed location z = −1 cm. The outcome appears in figure 9.

Zoom In Zoom Out Reset image size

Simulations with very different ℓ_sib behave similarly. After a short transient, the electrostatic repulsion between the two sibling branches strongly dominates their propagation. About 1 cm below the branching point, the trajectories of simulations with different ℓ_sib are barely separated. We conclude that ℓ_sib, which was introduced as a numerical parameter, does not influence the results much.

Let us now use our model to investigate the reconnection of streamer channels inside a tree. In a reconnection event, a streamer head is attracted toward a pre-existing channel. This should not be confused with streamer merging, where two streamer heads expand to form a single channel [14, 15].

Streamer reconnection has been observed both in laboratory discharges [21, 22] and in high-speed sprite observations [23–26]. Nijdam et al [22] reviewed the recorded examples of reconnection and extended them with new experimental data. Using stereoscopy, they were able to discriminate between actual reconnection and ambiguous observations resulting from projecting the 3D streamers into the camera plane. They concluded that reconnection of positive streamers in laboratory experiments is indeed frequent but consists in a thinner, slower streamer moving toward the channel of a thicker, faster streamer that had already contacted the cathode. After this contact, the ionized streamer channel charges negatively and attracts the streamer heads surrounding it, still positively charged. Although commonplace in the laboratory, this mechanism does not explain the observations of streamer reconnection in sprites, where a lower electrode does not exist. Here we will limit ourselves to the study of the latter kind of reconnection, where a lower electrode does not exist or is not essential. We henceforward restrict the meaning of reconnection to this type of event only. In this restricted sense, reconnection has not been unambiguously observed in laboratory experiments.

We frequently observe reconnection events in our model. Figure 10 shows an example where, for clarity, we searched for the earliest reconnection event in a set of 25 simulations. This selection biases the sample toward a higher amount of branching, as can be seen by comparing with figure 3. However, the pattern shown in figure 10 is generic to all the reconnection events that we found in our simulations. It consist of a lagging streamer being attracted to the stem of a sub-tree that has propagated much farther. The picture shows that the reason is that, as explained in section 3.2., the inner branches of the tree acquire a negative charge; usually, most of the channels in that volume are similarly negatively charged but if a lagging streamer propagates through the inner sections of the tree, its positive charge is attracted and reconnects to a negative, inner branch. To put it concisely, the extremal branches are attracted toward the internal ones.

Zoom In Zoom Out Reset image size

In figure 11 we zoom into the reconnection of figure 10 and plot two snapshots of the charge distribution. We see that as the head approaches the channel, it induces a significant, additional negative charge in the pre-existing channel. The relevance of these induced charges in a conductive channel was pointed out by Cummer et al [23]. Nevertheless, our simulations suggest that the initial attraction of a head toward a channel is possible only in cases where that channel has the opposite charge. The induced charges dominate only when the head is already very close to the channel.

Zoom In Zoom Out Reset image size

We speculate that reconnection (in our restricted sense) has not been observed in laboratory discharges because their innermost branches do not charge negatively or do not do it strongly enough. We offer two possible reasons for this. (a) The needle-electrode geometry most often employed in the laboratory, by imposing higher and divergent electric fields around the anode, discharges the negative charges in that region faster and reduces streamer interaction. (b) The reduced propagation length imposed by the cathode does not allow the tree enough time to reconnect. Most likely, there is a combination of both (a) and (b) at play; and finally, in the laboratory experiments [22], only for sparse trees with less than about 50 streamers can the full 3D structure be reconstructed, which gives a bias in the observations.

To investigate further whether we should expect to see streamer reconnection in laboratory experiments, we can tune the parameters in our model and make reconnection more or less likely. In particular, we may force the streamers to branch more or less frequently by varying the parameter ℓ_branch. We used values from 0.35 cm to 5.5 cm⁻¹ and for each value we run ten simulations up to the time of the first reconnection. The results are plotted in figure 12.

Zoom In Zoom Out Reset image size

For the standard value ℓ_branch = 2 cm the plot indicates that we need a gap of about 7 cm between electrodes to have a significant chance of observing reconnections; if ℓ_branch would increase to 2.85 cm, one would need a gap of more than 12 cm. Given the uncertainties and approximations in our model and point (a) discussed above we believe that laboratory discharges would also reconnect if given enough space.

A necessary condition for the numerical calculation of the model is that it converges for decreasing time step Δt. To check this, we run deterministic simulations (with ℓ_branch = 0) with an external electric field E₀ = 15 kV cm⁻¹ and various Δt. Figure A.1 shows the length of the streamer channel as a function of time; the simulations converge to a solution once the time steps are shorter than about 0.25 ns.

Zoom In Zoom Out Reset image size

Therefore in all simulations in this paper we use Δt = 0.25 ns.

We are calculating all interactions between pairs of charged nodes and therefore our computation time scales as O(N²). This is the main limitation on the size of trees that we can efficiently simulate. To overcome this limitation we also implemented the fast multipolar method (FMM) which is able to solve the electrostatic problem with O(N) computations up to an arbitrarily good approximation. However, the kernel in (9) is not the Poisson kernel for R ≠ 0 and although we restricted the FMM only for distant interactions with r_ij ≫ R, we run into problems around the cutoff. Besides, we found that due to the overhead of the FMM, it was advantageous only for N larger than a few thousands and all the simulations reported here are below that threshold. Each of the simulations that we show took a few hours on a modern desktop computer.