Simulation model for the propagation of second mode streamers in dielectric liquids using the Townsend--Meek criterion

A simulation model for second mode positive streamers in dielectric liquids is presented. Initiation and propagation is modeled by an electron-avalanche mechanism and the Townsend--Meek criterion. The electric breakdown is simulated in a point-plane gap, using cyclohexane as a model liquid. Electrons move in a Laplacian electric field arising from the electrodes and streamer structure, and turn into electron avalanches in high-field regions. The Townsend--Meek criterion determines when an avalanche is regarded as a part of the streamer structure. The results show that an avalanche-driven breakdown is possible, however, the inception voltage is relatively high. Parameter variations are included to investigate how the parameter values affect the model.


INTRODUCTION TO STREAMERS
Dielectric liquids are widely used for insulation of high power equipment, such as transformers, since liquid insulation has good cooling properties, high electrical withstand strength, and recovers from an electrical discharge within short time [1]. Electric breakdown in liquids is preceded by the formation of a prebreakdown channel called a streamer [2]. A partial discharge, a local electric breakdown, changes the electric field distribution, which could cause another local breakdown, and in this way, a streamer may propagate through a liquid. A streamer bridging the gap between two electrodes, for instance an energized part and a grounded part, lowers the electrical withstand strength and may cause a complete electric breakdown, possibly destroying the equipment [1].
A streamer consists of a gaseous and partly ionized structure, originating in one location and branching out in filaments as it propagates through the liquid. This structure may be observed through shadowgraphic or schlieren photography since its refractive index di ers from the surrounding liquid [3]. Streamers are classified as positive or negative, depending on the polarity of the initiation site. Streamer experiments are often carried out in needle-plane gaps since a strongly divergent field allows control of where the streamer initiates, the polarity of the streamer, and also enables the study of streamers that initiate, propagate, and then stops without causing an electric breakdown [2,3]. Conversely, in a gap with a uniform field, inception governs the breakdown probability, since an initiated streamer is always able to propagate the gap due to the high background field.
The nature of streamers has been investigated for decades [1][2][3][4][5][6][7][8][9], but is still not well understood since the complete picture involves numerous mechanisms at many scales. For streamer initiation, the importance of field-ionization [10], space-charge limited current [11], Joule heating [11], bubble nucleation [12], micro-bubbles [13], micro-cracking [13], and electron avalanches [11,14,15], have been discussed [2], and these mechanisms could be relevant for streamer propagation as well. However, the situation with a propagating streamer is even more complex. There is a gaseous channel and a liquid medium, with an interface in between. The channel is a partly ionized, low-temperature plasma, having a varying conductance [6,16]. For positive streamers, it is common to define four distinct modes of propagation, mainly characterized by their speed [2]. The 1st mode is only seen for very sharp points and propagates in a bubbly or bushy fashion with a speed of the order of 100 m s 1 , the 2nd mode is faster, of the order of 1 km s 1 , and has a branched or tree-like structure. The even faster 3rd and 4th modes propagates at speeds of the order of 10 km s 1 and 100 km s 1 , respectively. A streamer can change mode during propagation.
The development of models is important for improving electrical equipment as well as the prevention of equipment failure. An early simulation model for liquid breakdown uses a lattice to investigate the fractal nature of the streamer structure as a function of the electric field E [17], and has been expanded to incorporate needle-plane geometry [18], a 3D-lattice [19], statistical time [20], availability of seed electrons [21], and varying conductance of the streamer channels [22]. Charge generation and transport in an electric field have also been solved by a finite element method (FEM) approach, to simulate streamer propagation in 2D and 3D, adding impurities to generate streamer branching [23,24]. A major di erence between breakdown in gases and liquids is that a phase change is involved when making the streamer channel in liquids. The phase change is di cult to model, but it is possible to make approximations [25], or to focus on the plasma within the channel [26].
Both lattice and FEM simulations require considerable computational power, and therefore, the simulations are often done for either very short timescales or very simplified models. The work presented here is based on [27], which chooses a di erent approach. It is a computational model for 2nd mode positive streamers in non-polar liquids, driven by electron avalanches in the liquid phase. A point-plane geometry is modeled, with the point being a positively charged hyperbolic needle. Cyclohexane is used as a model liquid, since it is a well defined system used extensively in experiments [5,8,10,14,28].
The model and the theoretical background is presented in section 2, as well as the parameters and the algorithm used for the simulations. In section 3, the results are given and discussed. First a baseline is established, then parameter variations and alternative parameter values are investigated. A general discussion, outlining the weaknesses and strengths of the model, is given in section 4. Finally, the main conclusions are summarized in section 5. Appendix A contains additional details on the coordinate system used in the model.

FIGURE 1
The hyperbolic needle and a streamer head, with relevant variables shown. The distance to the plane is usually far greater than illustrated here.

SIMULATION MODEL AND THEORY
The model is built on the assumption that electron avalanches occur in the liquid phase, and that these govern the propagation of 2nd mode, positive streamers [27]. Applying a potential to the needle in a needle-plane geometry gives rise to an electric field. A number of anions and electrons, assumed to be already present in the liquid, are accelerated by the electric field. Subsequently, electron multiplication occurs in areas where the electric field is su ciently strong, turning electrons into electron avalanches. An avalanche is assumed to be "critical" if it reaches a magnitude given by the Townsend-Meek criterion [29], and the position of such an avalanche is regarded as a part of the streamer. Then the electric field is reevaluated, accounting for the potential of both the needle and the streamer. This work investigates liquid cyclohexane as the insulating liquid, with the option to add dimethylaniline (DMA) as an additive, but the model can be used for other base liquids and additives as well, if the parameter values are available.

Geometrical and electrical properties
A hyperbolic needle electrode with a tip radius r p is placed at a distance d g from a planar electrode, as illustrated in Fig. 1 where all important geometric variables are shown. In prolate spheroid coordinates ( , ⌫, ; a), a hyperboloid is represented by a single coordinate ⌫, see Appendix A for details and definitions. With these coordinates, the potential is (cf. (A15)) and the electric field is (cf. (A17)) where C i is a constant. The subscript i refers to a given hyperboloid (the needle or a streamer head), hence, the subscript in ⌫ i implies a transformation to a coordinate system centered at hyperboloid i, The constant C i (cf. (A16)) is given by the potential at the surface ⌫ i (r i ), which is valid for a sharp needle, r p~zi . For the needle, which is the applied potential. Calculating the electric field in (2) is the most expensive part of the computer simulation, although explicit calculation of the trigonometric functions can be avoided (cf. Appendix A).

Electrons and ions in dielectric liquids
Naturally occurring radiation is of the order of D r = 1 mSv per year [30] and may produce electron-cation pairs by ionizing neutral molecules. The production rate is [31] R e = D r ⇢ G , where the density ⇢ is 0.78 kg_l for cyclohexane. The yield G is usually given in events per 100 eV. Hydrocarbons typically have an ion yield G ion of about 4 [32], and for cyclohexane it is 4.3 [33]. However, the free electron yield G free is much lower, about 0.15 [33,34], which implies that most electrons recombines geminately. This gives a production of R e = 2.3 ù 10 8 m *3 s *1 . The recombination process is rapid, and the electron lifetime is [31] where ✏ 0 is the vacuum permittivity, ✏ r = 2.0 is the typical relative permittivity for hydrocarbons, r 0 is the recombination distance, e is the electron mobility, and e is the elementary charge. Inserting the thermalization distance (the most likely distance) r 0 = 5.9 nm [33] and a mobility e = 45 mm 2 V *1 s *1 [34,35], yields ⌧ r = 1.7 ps. The average drift velocity v d of an electron or ion is given by its mobility and the local electric field E, In liquids where the electron mobility is low ( e < 10 2 mm 2 V *1 s *1 ), the electron is regarded as localized, and electron transport is explained either through a hopping or a trapping mechanism [36,37]. The drift velocity is proportional to the electric field when the electric strength is low, however, for low-mobility liquids, it becomes superlinear in high fields [31,36]. The lifetimes of free electrons and ions can be related to the reaction rates. The reaction rate constants k r are found by the Debye relation [31,38], where ± is the mobility of the respective reacting species. This relation holds as long as the mobilities are low (< 10 4 mm 2 V 1 s 1 ) [31]. In cyclohexane, the ion mobility is of the order of 10 2 mm 2 V 1 s 1 to 10 1 mm 2 V 1 s 1 [10,11,33,[39][40][41] and the electron mobility is of the order of 10 mm 2 V 1 s 1 [33,34,42,43]. Using e = 45 mm 2 V *1 s *1 and ion = 0.1 mm 2 V *1 s *1 , yields k r = 4.1 ù 10 *13 m 3 _s for electron-ion recombination and k r = 1.8 ù 10 *15 m 3 _s for ionion recombination according to (8). This implies that there is a far greater number of anions than electrons. However, small impurities, such as O 2 , have higher mobilities [31]. The low-field conductivity for the liquid is given by the number density of charge carriers n i for species i and their mobilities, By assuming that the measured conductivity is due to ions only and that the ions are similar in number and mobility, the number density of the anions is which yields n ion = 6.2 ù 10 12 m *3 for = 0.2 pS_m [41,44]. A similar result is obtained by considering a steady-state condition, dn e dt = R e * k r n e n p * n e ⌧ a = 0 , where n e is the electron density, n p is the cation density, and t is the time. If the electron attachment time ⌧ a is large [45], which yields n e = 2.4 ù 10 10 m *3 . However, ⌧ a is assumed small, about 200 ns [24], which implies that n ion˘np . Using (12) with the ion-ion recombination rate yields n ion = 3.6 ù 10 11 m *3 , about an order of magnitude lower than what obtained from (10). With rapid attachment, (11) is and yields n e = 46 m *3 , which shows that the assumption n ion˘np holds.

Electron avalanches
The main concept the model is that electrical breakdown is driven by electron avalanches occurring in the liquid phase [8,14,27]. A number of anions, calculated by (10), is considered as the source of electrons by an electron-detachment mechanism. These electrons initiates the avalanches. As shown in section 2.2, the number of anions is far greater than the number of electrons, and it is also far greater than the number of electrons produced within a simulation (a volume less than 1 cm 3 and a time less than 1 s). The needle electrode and the streamer creates an electric field E. It is assumed that an electron will detach from its anion for field strengths above E d = 1 MV_m. The movement s of each electron or anion i is calculated by The simulation time step t is chosen low enough, typically 1 ps to 10 ps, to ensure that s is less than 0.1 µm. For a positive streamer, the negative charged species move towards higher field strengths. Increasing the electric field strength, increases the kintic energy an electron gains between colliding with molecules as well as lowering the ionization potential (IP) of the molecules [9], which increases the probability of impact ionization. As electron attachment processes dominate at low field strengths, an electric field exceeding E a = 0.2 GV_m is required for electron multiplication to be observed in cyclohexane [14]. An electron avalanche occurs when electron multiplication is dominant and the number of electrons N e grows rapidly. The growth of such an avalanche is modeled as [29] where ↵ is the average number of electrons generated per unit length. For discharges in gases, ↵ is assumed to be dependent on the type of molecules, the density, and the electric field strength [46]. Assuming that the same holds for a liquid, considering a constant liquid density [14,15], yields The maximum avalanche growth ↵ m and the inelastic scattering constant E ↵ are dependent on the liquid and are found from experimental data [14,47]. Equation (15) leads to an exponential growth of electrons in an avalanche, where N 0 is the initial number of electrons, and Q e is introduced as a measure of the avalanche size. At each simulation step, Q e for each avalanche is increased by For discharges in gases it is assumed that an electron avalanche becomes unstable when the electron number N e exceeds some threshold N c , which is known as the Townsend-Meek avalanche-to-streamer criterion [29]. In the model, an avalanche obtaining this criterion is removed and its position is considered as a part of the streamer channel. Assuming that an avalanche starts from a single electron, the criterion N e > N c is rewritten as The Meek constant Q c is typically 18 in gases [29,48]. A recent study on liquids found values in the range 5 to 20 when evaluating a number of experiments [47]. Another study found Q c = 23 by considering the field required for propagation [8], in contrast to the field required for initiation, which is more common.

Additives
Additives with low IP have proven to facilitate the propagation of 2nd mode streamers, since such additives lower the voltage required for propagation and for breakdown, whereas they increase the voltage required for 4th mode streamers [2]. This is likely a consequence of an increased number of branches, which may increase the electrostatic shielding and thereby reducing the electric field at the streamer heads [7,28]. To account for the e ect of low-IP additives on electron avalanche growth, the mole fraction c n of the additive and the IP di erence between the base liquid I b and the additive I a , is used to modify the expression for ↵ in (16) as [8] where the parameter k ↵ = 2.8 eV *1 is estimated from experiments. [8] For example, an additive with an IP di erence of 3.1 eV from the base liquid, in a concentration c n of 0.1 %, yields ↵ ® = 6.9↵. Equation (20) is derived assuming that ionization is caused by electrons in the exponentially decaying, high-energy tail of a Maxwellian distribution, and that the introduction of an additive does not significantly change the energy distribution [8].

Streamer representation
The model focuses on the processes occurring in front of the streamer. The streamer is represented by a collection of hyperboloids, approximating the electric field in front of the streamer. The streamer channel, and in particular its dynamics, is not included in the model. The streamer hyperboloids are referred to as "streamer heads", and the initial streamer consists of only one streamer head: the needle. The needle, one other streamer head, and relevant variables, are shown in Fig. 1.
The potential V at position r is given by a superposition of the potential V i in (1) of each streamer head, where the coe cients k i are introduced to account for electrostatic shielding between the heads. The electric field is found FIGURE 2 For given a streamer head i (shown), other positions are considered to be within, behind, in front, and/or within join distance.
in a similar manner, where E i in (2) is the electric field arising from streamer head i. The electric field arising from a streamer head is strongly dependent on its tip radius r p . Experiments have shown that there exists a critical tip radius for the inception of 2nd mode streamers, which is r p = 6 µm for cyclohexane [5,49]. Whenever an electron avalanche meets the Townsend-Meek criterion in (19), a new streamer head is added at the position of the avalanche. The potential at the tip of streamer head i is given by where V 0 is the potential at the needle, E s is the electric field within the streamer channel, and l i is the distance from the tip of the needle to the tip of streamer head i, again see Fig. 1 for definitions. Equation (23) is used to find C i through (4). The shielding coe cients k i ensure that the combined potential of all the streamer heads equals the potential at the tip of each streamer head, and are obtained by a non-negative least squares (NNLS) routine [50]. The problem actually solved numerically is stated in a slightly di erent form. Defining which only depend on the geometry and not on the potentials, (25) is rewritten as which is computationally more convenient to solve.
It is desirable to keep the number of streamer heads to a minimum since it is expensive to calculate the electric field from a head. Also optimization of the potential becomes more dicult and unstable as it tend to become a more overdetermined problem with more heads present, especially when the heads are close or "within" each other. Streamer heads located within another streamer head are removed, that is, if then streamer head j is removed, which is the same as being above the ⌫ 0 -line in Fig. 2. In addition, if the tip of one streamer head is within a certain distance d m of the tip of another streamer head, the heads are merged and only the streamer head closest to the plane is kept (see Fig. 2). Physically, this is motivated as charge transferred from one streamer head to another located closer to the grounded plane. Finally, since fewer heads implies less calculation and faster simulations, streamer heads with a shielding coe cient below a given threshold, are also removed. When k c is chosen su ciently low, only streamer heads that are to a large degree shielded by other heads are removed, and removing them have thus little e ect on the simulation results.

Region of interest
Anions, electrons, and avalanches are here referred to as "seeds". The seeds are placed as anions, but can become electrons or avalanches, depending on the local electric field strength, which is illustrated in Fig. 3. To save computational cost, especially for simulations in large gaps, seeds are limited to a region of interest (ROI) surrounding the leading tip, see Fig. 3. The ROI is a cylinder defined by a radius from the centerline (x 2 + y 2 = r 2 ), a distance in front of the leading streamer head, and a distance behind the leading head. Seed avalanches that obtain a critical size, seeds that collide with a streamer head, and seeds that fall behind the ROI, are removed and replaced by new seeds. A new seed is placed one ROI length from the old seed in the z-direction, with random placement within the ROI radius for the x-and y-coordinates. The seed density is thus kept constant as the ROI moves together with the leading streamer head. For a given configuration, it is possible to calculate the time t i for an electron to travel from a given point to the needle. This is achieved by numeric integration of v *1 d dl along an electric field line (constant ), using h ⌫ = dl_d⌫ (cf. (A12)), Similarly, the maximum avalanche size Q i , is computed by An illustration of (31) and (32) is found in Fig. 4. Both v d in (7) and ↵ in (16) are functions of the electric field E in (2), which makes numeric integration straightforward in prolate spheroid coordinates. The time, t i , provides an indication of how large the ROI should be. Given that a slow streamer may propagate at 1 km s 1 = 1 mm µs 1 , the ROI should be chosen so wide that seeds on the sides does not have enough time to collide with the passing streamer. According to Fig. 4, a width of V n * Needle tip curvature r n 6.0 µm Streamer tip curvature [5] r s 6.0 µm Field in streamer [6,51] E s 2.0 kV mm 1 Electron detachment threshold E d 1.0 MV m 1 Avalanche threshold [14] E a 0.2 GV m 1 Scattering constant [14] E ↵ 3.0 GV m 1 Max avalanche growth [14] ↵ m 200 µm 1 Meek constant [8] Q c 23 Electron mobility [42,43] e 45 mm 2 /Vs Anion mobility [11] ion 0.30 mm 2 /Vs Ion conductivity [41] ion 0.20 pS m 1 Base liquid IP [52] I b 10.2 eV Additive IP [53] I a 7.1 eV Additive IP di . factor [8] k ↵ 2.8 eV 1 Additive number density c a,n 0.0 1.5 mm gives about 1 µs before collision, both from the sides and from below. As the streamer should propagate about the same length, or more, in this time, is a reasonable value. However, a somewhat wider ROI should be used to account for a streamer propagating o -center, and for branched propagation.
Further, Fig. 4 shows that Q i is large in the front, but quickly declines for seeds behind the streamer head. This gives an indication on how far behind the streamer head an avalanche may obtain critical size, which is how far behind the streamer head it is interesting to extend the ROI. However, the ROI should also extend far enough behind the leading streamer head to enable the propagation of secondary branches. Even though t i and Q i give good indications of how big the ROI should be, it is important to verify the settings after the simulation, or vary the ROI to verify that the results are not a ected.

Parameters
The model parameters may be divided in two groups: physical parameters and parameters for the numerical algorithm. The values of the physical parameters summarized in Tab. 1 are given by the properties of the simulated experiment or based on values available in the literature for the base liquid (cyclohexane) and the additive (dimethylaniline). Since not all the parameter values are available and some are uncertain, a sensitivity analysis is carried out in this work to investigate the influence of individual parameters. Parameter values needed by the simulation algorithm, which are not based on physical properties, are given in Tab. 2 and include the size of the ROI and certain criteria for stopping a simulation. The initial setup is given by V n , d g , and r n . Then the number fraction of seeds n ion is calculated using ion and ion , according to (10), and whether a seed is considered as an anion, an electron, or an avalanche is given by E d and E a . The electron multiplication probability is given by (16), using E ↵ and ↵ m . If an additive is present, then (20) is also applied, where I b , I a , c a,n , and k ↵ are used. Equation (18) gives the growth of an avalanche, using t and e . Finally, the Townsend-Meek criterion, stated in (19), uses Q c to evaluate whether the avalanche has obtained a critical size. The streamer branching is regulated by d m and k c , by (29) and (30), while the streamer head potential, and thus also the electric field at the tip, is dependent on E s and r s through (23).

Algorithm
A simulation begins by reading an input file that is used to initialize the various data classes used by the program, including random placement of seeds within the ROI, thereafter, a loop is executed until the simulation is complete. These main steps are shown in Fig. 5. The first and most expensive step of the algorithm is the update of the seeds, which is detailed in Fig. 6. First, the electric field is calculated for all seeds (each anion, electron, and avalanche). All the avalanches are treated separately in a loop, where they are moved, the electrons are multiplied, and the field is calculated for their new positions. This loop, in Fig. 6, is performed until either N MSN steps are done, an avalanche becomes critical (obtaining the Townsend-Meek criterion), or an avalanche collides with the streamer. Then, all other seeds (anions and electrons) are moved, using a time step equal to the total time used by the avalanches. The next step in Fig. 5 is to update the streamer structure. Any critical avalanches are added to the streamer, and the streamer structure is optimized by removing heads using (28) and (29) and correcting the scaling using (27) to set k i for each streamer head. Thereafter, if there is a new leading streamer head, the ROI is updated. In the "clean-up" part, seeds behind the ROI, critical seeds, and seeds that have collided with the streamer, are removed and replaced by new seeds. A number of criteria can be set to determine when the simulation loop in Fig. 5 should end. For instance, total simulation time, total CPU time, and number of iterations. However, simulations presented in this work ended for one of three reasons: the leading head reached the planar electrode (z i < z min ), low propagation speed (< v min ), or long time between critical avalanches (> t ava max ). The final step of the loop is saving data, and finalizing a simulation ensures that all temporary data is properly saved to files.
The implementation has been done in Python [54] using NumPy [55] extensively. During initialization, the seed for random numbers is set in NumPy to ensure reproducible results. The input parameters are given in a JSON-formatted file, which is used for initiation of the simulation. Simulation results are saved with Pickle and illustrated using Matplotlib [56].

SIMULATION RESULTS AND DISCUSSION
The model involves numerous parameters, some of which is given by the experimental setup (e.g. gap distance), others by properties of the liquids (e.g. mobilities), and some are purely for the simulation procedure (e.g. time step). In the first part, the default parameters given by Tabs. 1 and 2 show the basic behavior of the model. Thereafter, a sensitivity analysis is presented, indicating the influence of various parameters. Mainly the propagation speed is used to indicate the di erences, but the number of streamer heads, their scaling k i , the propagation length, and the degree of branching are also investigated. Ten simulations are carried out at each voltage, using the numbers 1 to 10 in the random number generator generating di erent initial configurations of the seed distribution.

Simulation baseline
Simulations have been performed for a range of voltages, using the parameters in Tabs. 1 and 2. These simulations are used as a baseline in the sensitivity analysis. As seen from the streak plots in Fig. 7, a voltage exceeding 60 kV is needed for a breakdown. For lower voltages, the streamer propagates less than 100 µm before the simulation is terminated, either because of waiting too long for an avalanche or because of very slow propagation speed. Above the breakdown voltage, the time to breakdown is reduced as the voltage is increased, and the

FIGURE 5
The main steps of the simulation algorithm. The algorithm for the seeds is detailed in Fig. 6   streamers tend to accelerate towards the end of their propagation. The average propagation speed, shown in Fig. 8 tells a   FIGURE 8 Streamer average speed versus leading head position, that is, the average gradient of the "streaks" shown in Fig. 7.
similar story, but it also indicates that the propagation speed slows down a bit after the first few steps. The speed reduction is possibly due to branching, however, by looking at the

FIGURE 9
Streamer trails, xz-and yz-projection for a range of voltages 60 kV to 120 kV, using the same legend as in Fig. 7. Each dot represents the position of a streamer head at some point of the propagation. The streamers are plotted with an o set to improve the readability.

FIGURE 10
Actual streamer head scale k i (left) and total number of streamer heads (right). Data are taken every 1 % of the gap. The dashed lines are moving averages calculated by loess-regression [57].
streamer in Fig. 9, it is clear that the degree of branching is very low, but the streamer gets thicker with increasing voltage. This implies that even though branching is not apparent, there are several streamer heads present. The number of streamer heads may increase when the electric field strength increases (at higher voltages or closer to the plane) as seen in Fig. 10. Values of k i lower than one implies that the streamer heads shield each other to some degree (cf. (21)), as seen in Fig. 10, but not enough to stop a propagating streamer. It is of interest to investigate how the leading head is a ected by shielding, and the average scaling indicates this. The propagation speed can be described by the time it takes to get a critical avalanche

FIGURE 11
The leading streamer head is moved in a sequence of discrete "jumps" in the z-direction. The average jump length and the standard deviation of the jumps are found for each individual simulation. The dotted lines are interpolated to the average, and the bars covers the minimum and maximum values for ten simulations at the same voltage.
in front of the leading streamer head combined with the distance the leading head is moved, where the latter is presented in Fig. 11. Increased voltage increases both the maximum and the average propagation "jumps", especially when the streamer is in the final part of the gap. The propagation speeds in Fig. 8 are somewhat low for 2nd mode streamers, which should be 1 km s 1 to 10 km s 1 [2]. Many, if not most, of the simulation parameters a ects the propagation speed. In the case of the electron mobility e , it is easy to see that the propagation speed is directly proportional to e , since it only a ects the movement of the electrons (cf. (14)). For most other parameters, it is not that simple.

E ect of avalanche parameters
The avalanche mechanism is the most important part of the model. For this reason, parameters relevant to the avalanche growth, given in (16) and (19), are especially important. To get an avalanche, however, a seed electron is needed. A doubling of the concentration of seeds n ion , gives about a doubling in the propagation speed, as seen in Fig. 12. The figure shows the average speed for the mid 50 % of the gap, that is for a position from 0.75 mm to 2.25 mm. Since streamers terminated in the first quarter of the gap are not shown, the figure also indicates that the breakdown voltage is dependent on n ion , as increasing n ion allows propagation at lower voltages. The streamer is represented by one or more heads, and propagates as new heads are added in front of current heads. As such, the leading head moves in a series of discrete "jumps". The average streamer  head jump length seems independent of n ion , indicating that the linear increase in propagation speed is caused by a reduction in the time required for an electron to become a critical avalanche. At n ion = 2 ù 10 12 m *3 , the average distance between seeds is 79 µm, while the average jump length is about 6 µm, so n ion would have to be increased by some orders of magnitude to a ect the streamer jump distance. Inhomogeneities on the order of 10 11 m 3 was introduced by [24] to explain branching, but this e ect is not found here. An upper estimate on the ions available can be calculated from (12) by using G ion instead of G free when calculating R e in (5) and using a low estimate of k r = 10 *3 mm 2 V *1 s *1 [24,40], yielding n ion = 1.8 ù 10 13 m *3 and an average distance of 38 µm between seeds. As such, the simulations in Fig. 12 cover the most interesting range.
The baseline results in section 3.1 do not show any stopping of streamer propagation mid-gap. The streamers either  stop within the first 100 µm or cause a breakdown. This occurs when the supply of electrons is constant and E s is too low to create a high voltage drop along the streamer. Increasing the electron detachment threshold E d reduces the number of electrons available, which in turn reduces the density of electrons as electrons are swept out, see Fig. 13. This results in a negative feedback loop where a lower density of electrons decreases the speed (Fig. 12) and the decreased speed results in a lower rate of ions turning into electrons. The propagation length is shown as a function of the needle potential and E d in Fig. 14.
By considering E d = 15 MV m *1 , three di erent regimes is identified. Up to 70 kV, a few avalanches may occur, but then the propagation stops. Above 90 kV, the propagation is fast enough to provide a stable rate of new electrons, enabling the propagation to continue. In between, the initial electrons allow the streamer to propagate, but the electron density is decreasing and the streamer eventually stops.
The electric field is important for electron movement and multiplication, and E ↵ in (16) is therefore an important parameter. The strong influence of E ↵ is seen in Fig. 15, where the propagation speed may increase by an order of magnitude when E ↵ is reduced by 50 %. This makes sense as E ↵ enters exponentially in (16). The propagation speed of 2nd mode streamers is weakly dependent on the applied voltage [2], however, for E ↵ = 1 GV m *1 in Fig. 15, the dependence is much stronger than for the other values. Reducing E ↵ facilitates streamer propagation and the breakdown voltage is thus strongly influenced. Both E ↵ and ↵ m are based on experimental results, and are very important to the model. Instead of varying ↵ m , however, the Meek-constant Q c is varied. From (16), (18) and (19), it is clear that the avalanche size Q e is linearly dependent on ↵ m , which implies that doubling Q c has the same e ect as halving ↵ m . The speed is not as a ected by Q c as intuitively expected, see Fig. 16, and changing Q c by a factor of 4 only changes the speed by a factor of 2. However, Q c cannot change much before the simulation becomes unphysical. For instance, consider a conducting sphere of r = 6 µm with a charge q = exp(Q c ). The electric field at the surface is where e is the electron charge and ✏ is the permittivity. For Q c equal 15, 20, and 25, the electric field becomes

FIGURE 17
Streamer propagation length as a function of needle potential and electric field in streamer channel E s . Each marker is a simulation and the dotted lines are interpolated to the average. Note that up to 8 kV m 1 , the results overlap to a high degree.
6.5 ù 10 7 V m 1 , 9.7 ù 10 9 V m 1 , and 1.4 ù 10 12 V m 1 , respectively. Increasing Q c by a little gives too high fields, and a decrease results in low fields. This can, however, be "fixed" by changing the radius. For instance, Q c = 15 and r = 1 µm, results in 2.4 ù 10 9 V m 1 , which is more reasonable. To consider the electron avalanche as a charged sphere is of course a simplification, but the majority of the charge does build up over a length of some µm, and this is also the size used for the streamer heads, which makes the analogy reasonable. While it would seem like increasing Q c does not make sense, one should remember that it actually has the same e ect on the model as decreasing ↵ m , and the value of that parameter is not certain. For instance, according to [14], ↵ m = 200 µm *1 , but [47] finds ↵ m = 130 µm *1 , however, the latter study also finds E ↵ = 1.9 GV m *1 , and changing this parameter has a big impact on the model, as discussed above.

E ect of streamer parameters
The streamer structure is responsible for propagating the electric field from the needle into the gap. The electric field in the streamer channel E s gives a voltage drop from the needle to the streamer head. The electric field in front of a streamer head is also dependent on the tip radius of curvature r s and the potential scaling of the streamer head k i . The scaling depends on the potential and position of all the streamer heads, that is, the entire "streamer". Both the streamer head merge distance d m and the potential shielding threshold k c may be important for the streamer configuration.  Fig. 14 demonstrates streamers stopping as a result of a reduction in the seed electron density, however, it is common to explain stopping as a result of an electric field E s in the streamer channel resulting in a lower field strength at the streamer head [16]. A high E s is needed to a ect the results (see Fig. 17), conversely, when E s is low, the streamer either stops quickly or causes a breakdown. When E s is high, the propagation speed is reduced throughout the gap and the propagation may stop somewhere in the gap, see Fig. 18 for E s = 16 kV mm *1 , which is in contrast to Fig. 7    the propagation, but becomes important when a streamer has reached some length. When E s = 8 kV mm *1 , the potential is reduced by 24 kV across the gap, but this e ect is barely seen (Fig. 17), since only a few streamers stop mid-gap. However, at 16 kV mm *1 the e ect is clearly present as many of the streamers stop mid-gap. Notice that at 75 kV to 85 kV, in Fig. 17 the average propagation length is increased from about 1.7 mm to 2.6 mm, giving an apparent electric field of only 11 kV mm *1  and not 16 kV mm 1 . This is perhaps an e ect of the field increasing as the gap is getting smaller. Also, actual experiments show stopping lengths that are increasing linearly with voltage in the first part of the gap, followed by more scatter and superlinear behavior towards the end of the gap [8,51,58]. This behavior is not seen in Fig. 17, possibly because E s is kept constant in the simulations, while it has been found to vary with applied voltage [6]. Streamers are subject to re-illuminations, associated with current pulses, which could change the electric field in the streamer channel, however, the propagation of the streamer head seems to be una ected by these e ects [6].
The curvature radius r s of a streamer head is an interesting parameter since a sharper tip gives a higher field and a larger volume where electron multiplication may occur. Changing r s from 1.5 µm to 12 µm only changes the speed by a factor of 2, see Fig. 19. Further increase to 24 µm decreases the speed, and increases the breakdown voltage. Simulations with smaller r s tend to have more streamer heads, scaled to a lower potential, than the simulations with a larger r s , indicated in Fig. 20, although the e ect is not visible for the smallest r s in that figure. The increased number of streamer heads seems to act as a regulating mechanism, however, the number of branches is not increased, but there are more streamer heads present simultaneously in the same branch. This is similar to the situation in Figs. 9 and 10, where an increased voltage does not increase the number of branches, but instead increases the streamer thickness.
An increase in voltage increases the speed (Fig. 8) as well as the number of streamer heads, while decreasing the scaling of the heads as demonstrated in Fig. 10. The parameters k c and d m are used to remove streamer heads, and therefore they could have a big impact on the model, since the scaling, which the electric field depends on, is strongly dependent on the number of streamer heads as well as their configuration. Also, these parameters are purely a consequence of the model, and do not have an origin in a physical property. Simulation results for varying k c are found in Fig. 21 and show that the propagation speed is not that a ected, except for k c = 40 %. This figure also indicates that the breakdown voltage is una ected, since all the values of k c are present for all the voltages. Setting k c = 40 % restricts the streamer to one head in most situations, and keeping two heads in rare occasions, which gives an upper bound to the propagation speed for each voltage. From a computational point of view, it is preferable to set k c high as fewer streamer heads implies less calculation. From a physical point of view, however, it does not make sense to just remove charges from the system, so k c should be reasonably low. According to Fig. 21, k c can be as high as 10 % without any particular impact on the results.
The influence of the streamer head merge distance d m is shown in Fig. 22. For the lower values, many streamer heads are present at the same time, which in turn lowers the potential scaling of each head, increases the breakdown voltage, and moderates the propagation speed. Increasing d m increases propagation speeds, up to the limit where there is mainly just a single active streamer head. Fig. 22 also indicates that at low voltages, the streamers propagate with a single head, but when the voltage is increased and more heads are possible, the propagation speed is moderated. As d m is increased, the voltage needed to have several heads is also increased, and the propagation speed is thus higher. The set of streamers presented in Fig. 23 shows that the thickness of the streamers is dependent on k c and d m , which is an indication of the number of streamer heads present during propagation. However, the figure does not indicate a change in the number of major branches.

E ect of additives
Adding small amounts of an additive increases the electron multiplication according to (20). The e ect should be similar to an increase of ↵ m , or a decrease in Q c , as discussed above and shown in Fig. 16. This is indeed the case, the propagation speed increases and the breakdown voltage decreases with increasing content of an additive with low ionization potential, see Fig. 24. When the liquid consists of c a,n = 10 % additive (mole fraction) it cannot be argued to be a "small amount" of additive. Even as little as 1 % could be too much. As mentioned in section 2.4, an addition of just 0.1 % increases the avalanche growth by a factor of 6.9, when using (20) and the parameters in Tab. 1. A decrease in breakdown voltage and an increase in propagation speed is also found in experiments with low-IP additives [3,8,28], however, increased branching is also seen in the experiments in contrast to the simulation results here.

Increased speed and branching
The above sections illustrate how the model behaves and how it is a ected by the various parameters. In order to reduce the initiation voltage and increase the propagation speed, the avalanche parameters are changed to E ↵ = 1.9 GV m *1 and ↵ m = 130 µm *1 , and the number of seeds is increased to c s = 8 ù 10 12 m *3 . In addition, the merge distance is changed to d m = 12.5 µm and the streamer head tip radius to r s = 3 µm in order to facilitate branching. Also, using E s = 8 kV mm *1 should be enough for some of the streamers to stop mid-gap. Most of the predicted results are found: the speed in Fig. 25 is clearly increased compared to Fig. 8, the amount of small branches is larger in Fig. 26 than in Fig. 9, and the decrease  in streamer head scaling and increase in streamer head number is seen by comparing Fig. 27 to Fig. 10. The propagation voltage is somewhat lower than the base case, around 60 kV. The streamer propagation begins at high speed, then slows down towards the middle of the gap, before the speed increases towards the end of the gap, see Fig. 25. This change does not seem to be correlated to the number of streamer heads, which is fairly constant for most of the propagation (Fig. 27). Branching may have an e ect, and streamer branching is illustrated in Fig. 28, showing 6 snapshots of a single simulation. As the streamer splits into two major branches, the number of electron avalanches surrounding the streamer heads decreases. The branches propagate at di erent speeds, and the faster one gains a higher potential and thus creates more electron avalanches. As the two branches approaches the end of the gap, one gains speed, while the other one stops.

DISCUSSION OF THE MODEL
Using a Laplace field is of course a simplification compared to a Poisson field. In fact, neither positive nor negative charges are accounted for in the model. The potential is simply calculated by assuming a constant field in the streamer channel, and then superimposing the streamer heads. Including the charge of the avalanches and the ions left behind could improve the model.

FIGURE 29
Maximum avalanche size at a distance from the needle tip. The unmarked lines use the same values as the baseline simulations from [14], the < indicates parameter values from [47] as used in section 3.5, and the † indicates formulation and parameter values from [15].

FIGURE 30
Total computational time for the simulations shown in Fig. 21 Streamers that terminate in the beginning of the gap require little time, while streamers that slowly bridges the gap requires the most computational time.
For the needle and the streamer heads, using a space charge limited field (SCLF) [59,60] would provide a more physically correct field distribution, but would also increase the computational requirements drastically. Using an SCLF rather than a Laplace field, gives a reduction of the electric field where the field is the strongest, since the maximum field is limited [59], with a corresponding increase everywhere else. The SCLF is time-dependent [60], and the e ect increases with time until an steady-state is obtained. The overall e ect on the model would be an increase in average jump length, as most jumps would be longer and the shortest ones would not occur. While an SCLF can give more accurate results for slow streamers, a Laplace field could be good enough for fast streamers, since the SCLFregion expands at some finite speed. However, the avalanche parameters in (16) were estimated using a Laplace field, so the current model is internally consistent.
The inception of 2nd mode streamers has been estimated to somewhat less than 15 kV for cyclohexane [5], however, for a propagating 2nd mode streamer, 33 kV was found for a 10 mm gap [8]. Since the model uses this as a criterion for inception (getting a critical avalanche, but no movement), a high propagation voltage is actually to be expected. This is well illustrated by the maximum avalanche size in Fig. 29, obtained by integration of ↵. Streamer propagation is possible when Q f > Q c (cf. (19)). The baseline simulations are performed inserting parameters from [14] in (16) to calculate ↵. At 33 kV, the maximum possible streamer jump is less than a µm, however, at 60 kV (the breakdown voltage), the value is increased to about 6 µm, possibly indicating that a strong field is needed over some distance. Changing to parameters from [47] decreases the propagation voltage by increasing the possible jump length, however, the decrease is not enough to enable for inception of 2nd mode streamers at 15 kV. As such, Fig. 29 indicates that streamer inception at 15 kV is not possible with this model when considering a Laplace field, calculating the electron multiplication with (16), and using the Townsend-Meek criterion for inception of 2nd mode streamers. Using the parameters of either [14] or [47] gives too low avalanche size. According to [15], the correct way of calculating electron multiplication in a dense medium is where I is the ionization potential, e is the electron charge, and E ⌫ is given by properties of the liquid. With this formulation, electron multiplication is more dependent on the electric field, implying that the electron avalanches become shorter, are closer to the streamer heads, and grow faster where the field is strong, which is illustrated in Fig. 29 using values for n-hexane [15]. The propagation velocity is somewhat low, which is to be expected since the inception voltage is too high. Changing parameters to values that lowers the inception voltage also increases the speed at a given voltage. As mentioned, the speed is proportional to the electron mobility, and it is the low-field mobility that has been used. For low-mobility liquids, such as cyclohexane, the mobility is expected to have a superlinear dependence on the electric field [36,43]. For this reason, one study multiplies the mobility by 2.5, to make it similar to the gas phase mobility [25], which would increase the streamer propagation speed by the same factor. Conversely, limitations to the maximum speed of electrons have been introduced [61], which would e ectively control the maximum speed of a streamer branch. The speed is also proportional to the concentration of seeds (see Fig. 12), which was calculated from the low-field conductivity of the liquid (see (10)). However, for breakdown in non-polar liquids, the conductivity is not important [2], and hence, it seems unreasonable for this parameter to be as important as demonstrated here. The equilibrium density of ions can also be calculated based on cosmic radiation (17), but obtaining > 10 11 m *3 ions, when the production is Ì 10 8 m *3 s *1 , implies that a long time is needed. It is therefore an approximation to simulate a situation where this density is kept constant. By changing the simulation conditions such that all the gap is included in the ROI and such that seeds are not replaced, it can be verified that the seeds present at the beginning of the experiment is not enough. They are swept out very fast if they are electrons and not ions. Increasing E d so that most seeds remain as anions changes this by allowing the low-mobility anions to live longer before entering the high-field area and ionize into molecules and electrons. Even so, it seems clear that some mechanism for generation of new seeds is warranted. New seeds could be generated in the high-field areas, and near the electrodes. The Zener model [62] for breakdown in solids has been used also for charge generation in liquids [24]. However, photoionization could also have an important role in the generation of new charges [2,7]. In addition, when ionizing neutral molecules, the field-dependent ionization potential [9] should also be taken into account. This kind of additions add complexity to the model, but Monte Carlo (MC) [63] methods can aid in keeping the added computational cost low. There are also some parts of the current model where MC could be reasonable to use. For instance, for electron detachment from an anion and for avalanche growth from a single electron, since a large number of electrons is needed to model an avalanche through the average growth ↵.
The degree of branching is lower than desired, with more or less only one major branch, and thus the simulations resemble more the 3rd mode or the start of the 4th mode than the 2nd mode of a streamer. It is worthwhile noting that streamers branch far less in cyclohexane than in mineral oil, but the addition of low-IP additives increases the branching [28]. The shapes of the simulated streamers do resemble the shape of streamers in longer gaps [28], however, while including additives in the model increases the propagation speed, the degree of branching is not increased. Although branching is thought of as a mechanism for regulating the propagation speed, it could be the other way around. With nothing to hold it back, the foremost head should have the strongest electric field and the fastest propagation. If something is regulating the speed or field of the foremost head, however, then other heads are given a better chance of propagation, increasing the number of branches, which in turn may regulate the electric field of all the branches. In the present model, there is nothing holding the foremost head back, since the only time scale included is that of the electron avalanches. If, for instance, the time required for bubble nucleation or the time for charges to move through the streamer structure (streamer dynamics) is important, it may result in a disadvantage for the foremost head. This is, however, not included, and the potential of each streamer head is instantly updated each simulation step. The shape chosen for the streamer heads could also be a major reason for the low degree of branching. For a hyperboloid, the electric field declines as r *1 in front, and the high-field region extends much further in the front than on the sides. Conversely, the field from a monopole declines like r *2 in all directions, and could as such facilitate branching. In such a model, however, the high field would be in a region closer to the streamer heads, making an SCLF approach even more relevant.
The simplicity of the presented model comes with several limitations, as discussed above, however, a simple model is also a good place to start. It makes it possible to identify whether a certain mechanism is important or not at a relatively low computational cost. Consider Fig. 30, which shows that the computational time for breakdown streamers averages to about one hour, using a single core on a regular desktop computer. The simulation time is of course strongly dependent on the number of seeds, streamer heads, and simulation steps, but with such a low base case, it is possible to perform a lot of simulations to gather statistics on a normal desktop computer. Contrary to lattice models, the presented model is based on physical processes, and the results are thus easier to evaluate. FEM models may be better in the end, but for now, such models cannot model a complete breakdown. They are also simplified, for example in the sense that phase changes are not accounted for [61]. Both lattice and FEM models demands much computational power and the mesh size becomes an important parameter, however, this is avoided in the model presented. Instead of dealing with processes at discrete point or in discretized elements, the model deals with discrete points that move. This approach makes sense when considering charge generated by electron avalanches at some distance from the streamer structure, or a streamer moving in discrete steps. For details on processes inside or very close to the streamer, however, a FEM approach seems more reasonable, and could provide valuable input to models on a larger scale.

CONCLUSION
A simple simulation model for streamer propagation has been presented. The streamer is represented by a collection of hyperbolic streamer heads, and is responsible for propagating the electric field from the needle electrode. In high-field areas, electrons detach from ions present in the liquid, and may turn into avalanches. If an avalanche meets the Townsend-Meek criterion, a new streamer head is added at its position, causing the streamer to propagate. As demonstrated, the model has some limitations, the inception voltage is too high while the degree of branching is low. These issues are discussed and explained, and directions for a systematic way of further developments are described. The main feature missing in the model is a proper representation of the dynamics of the streamer channel, however, the charge generation and the electric field calculation can be improved as well. The approach to streamer propagation applied here is di erent from that used by other models. The principle behind the model is simple, it is founded on physical mechanisms, and provides interesting information about how an avalanche-driven breakdown may occur. The simple model has its advantages in that it can be used to identify important mechanisms, without demanding excessive computational power.
Prolate spheroid coordinates involves a set of hyperbolas and ellipsoids revolved around the center axis, forming hyperboloids and prolate spheroids. The two focal points, of the hyperbolas as well as ellipsoids, are located at a distance a from the plane. The hyperbolic coordinate is À [0, ÿÎ, the elliptic coordinate is ⌫ À [0, ⇡], and rotation about the center is given by À [0, 2⇡]. The definition used here is and a constant ⌫ gives a hyperbola, z 2 a 2 cos 2 ⌫ * x 2 + y 2 a 2 sin 2 ⌫ = 1 . (A5)

FIGURE A1
In prolate spheroid coordinates, spheroids (blue) are given by a constant , and hyperboloids (red) have a constant ⌫. Here, ⌫ is given in units of ⇡.

FIGURE A2
Electric potential (left) and electric field strength (right), for a region close to a needle (center, gray) placed 10 mm from a grounded plane. The contour lines give a qualitative impression of how the respective magnitudes change as a function of position. A linear scale is used for both sides, and the magnitudes are linearly dependent on the potential of the needle.
Transformation from Cartesian to prolate spheroid coordinates is obtained through 2a cos ⌫ = p * m , (A7) where p =˘x 2 + y 2 + (z + a) 2 , (A9) m =˘x 2 + y 2 + (z * a) 2 , and are the distances between a given point and the two focal points. Prolate spheroid coordinates exists in many forms. In some cases, it is easier to work with substitutions such as ⇠ = sin ⌫, however, starting with trigonometric functions allows for greater flexibility through relations such as sin 2 + cos 2 = 1.
Scale factors h are useful to define when transforming between coordinate systems. The scale factor for ⌫, for instance, is found from Solving this, and the similar expressions for the other coordinates, yields h = a sinh sin ⌫ .
These are useful when defining the spatial derivative, The electric potential V and the electric field E are found by solving the Laplace equation, ( 2 V = 0. For a system where the hyperboloids represent equipotential surfaces, V = V (⌫), the Laplace equation is satisfied for [64] V (⌫) = A + C ln tan where the constants A and C are defined by boundary conditions. Given V (⌫ = ⇡_2) = 0 at the xz-plane and V (⌫ = ⌫ 0 ) = V 0 at the ⌫ 0 -hyperboloid, yields A = 0 and C = V 0 ln tan(⌫ 0 _2) .
Consequently, the electric field E = *(V becomes where Ç ⌫ is unit length in the direction of ⌫, Equations (A15) and (A17) are both illustrated in Fig. A2. The figure shows the di erences in behavior between the electric potential and the electric field, the latter increases rapidly close to the tip of the hyperboloid. Explicit transformation between Cartesian and prolate spheroid coordinates requires trigonometric and hyperbolic functions, which are costly when it comes to computations. There is, however, no need to calculate , ⌫, and explicitly, as both the potential (A15) and the electric field (A17) may be obtained by using (A6) and (A7), and trigonometric relations such as 2a sin ⌫ = 2a˘1 * cos 2 ⌫ =˘4a 2 * (p * m) 2 . (A19)