A kinetic Monte Carlo study of positive streamer interaction with complex dielectric surfaces

We present a computational study of positive streamers in air propagating over dielectric plates with square channels running orthogonal to the propagation direction. The study uses a newly developed non-kinetic Particle-In-Cell model based on Îto diffusion and kinetic Monte Carlo, which does not introduce artificial smoothing of the plasma density or photo-electron distributions. These capabilities permit the computational study to use high-resolution grids with large time steps, and also incorporates geometric shielding for particle and photon transport processes. We perform Cartesian 2D simulations for channel dimensions ranging from 250 µm to 2 mm, and track streamers over a distance of 4 cm and times ranging up to 300 ns, for various voltages ranging from 15 kV to 30 kV. These baseline simulations are supplemented by additional studies on the effects of transparency to ionizing radiation, streamer reignition, dielectric permittivity, and 3D effects. The computer simulations show: 1) Larger channels restrict streamer propagation more efficiently than narrow channels, and can lead to arrested surface streamers. 2) Geometric shielding of ionizing radiation substantially reduces the number of starting electrons in neighboring channels, and thus also reduces the onset point of streamer reignition. 3) Decreasing the streamer channel separation leads to slower streamers. 4) Increasing the dielectric permittivity increases the discharge velocity. The results are of generic value to fields of research involving streamer-dielectric interactions, and in particular for high-voltage technology where streamer termination is desirable.


Introduction
Streamer discharges [1] are fascinating and complex phenomena that occur in a wide range of environments, from lightning * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. strikes in the atmosphere to plasma technologies in industry. At their core, streamers are narrow channels of ionized gas that propagate through dielectric media, driven by high electric fields at the streamer tips. These structures exhibit a rich variety of dynamics, including branching, propagation, extinction, or bridging into a leader or spark.
Because of their unique properties, streamers have found numerous applications in fields such as materials processing, pollution control [2], and aerodynamics [3]. For example, atmospheric pressure plasma jets based on streamer discharges have been used to promote wound healing [4]. However, despite the significant progress that has been made in these areas, many fundamental questions regarding the underlying mechanisms that govern streamer behavior still remain.
Solid dielectrics can affect the initiation, propagation, and termination of streamers through a variety of mechanisms, such as streamer attraction through dielectric polarization, formation of sheaths [5,6], or geometric restriction of the available ionization and photoionization volumes. The interaction between streamers and dielectrics is an important area of research that has broad implications for both fundamental physics and practical applications, e.g. for applications in plasma medicine or high-voltage technology.
Despite the many advances in our understanding of streamers and their interactions with dielectrics, there are still significant challenges in modeling and simulating these phenomena. For example, accurate modeling of streamer propagation requires a detailed understanding of the electric field distribution, ionization processes, and electron transport. Additionally, the complex geometries and material properties of dielectrics can make it difficult to accurately capture the interactions between streamers and their environment. Addressing these modeling challenges requires a multi-physics approach that captures the behavior of the streamer from the appearance of the first electron to the propagation and termination stages.
Recent research has investigated the complex interaction between streamers and solid dielectrics, with several recent studies highlighting the important role of dielectric surface patterns. For example, Meyer et al [7,8] have presented a combined computational and experimental study for streamers propagating over non-planar surfaces consisting of planar dielectrics with circular and square channels running orthogonal to the streamer propagation direction. Meyer et al [8] used high-speed imaging and performed computer simulations, and an experimental follow-up study [9] demonstrated the relevance of these effects for high-voltage technology (albeit in a different field configuration). Tangentially connected computational studies have been presented by Wang et al [10] and Konina et al [11], the latter focusing on streamer interaction with pores and droplets. It warrants mention that a common theme in all these studies is that they employ a continuum description of photoionization (formally speaking, a first-order truncation of the radiative transfer equation). While this formulation [12] has been highly enabling for the low-temperature plasma community, it has the disadvantage of causing artificial smoothing of the photoionization profile, and being susceptible to radiative leakage around geometric obstructions. The importance of discrete photoionization on overall streamer morphology has been demonstrated by several authors [13,14]. Meyer et al [8] also point out the deficiencies of the continuum model in the context of morphologically complex surfaces, as it can lead to an artificial increase of the plasma density behind objects, which puts into question some computational observations (like streamer reignition).
In this paper we consider a computational study of positive streamers propagating over square channels similar to the ones in [8]. As the study in Meyer et al [8] only considered a single square-channel surface, we focus here on additional main parameters: (i) The size of the square channels, in particular channels that are substantially larger than the streamer dimensions. (ii) The applied voltage. (iii) Geometric shielding of ionizing radiation. (iv) Distance between the channels.
(v) The effects of dielectric permittivity.
To provide insight into these physical dependencies, our work uses a newly developed non-kinetic Particle-In-Cell model based inÎto diffusion and kinetic Monte Carlo (KMC), termed Ito-KMC, which has the advantage of being fully stochastic and discrete [15]. The outline of this paper is as follows: The numerical model is presented in section 2 and various computational results are presented in section 3. Specifically, section 3.1 provides calculations of the mean streamer velocities, and section 3.2 provides an overview of the morphological evolution of the streamers for the various surfaces. Section 3.3 provides a brief study of the role of opaque boundaries which obstruct ionizing radiation. A parameter study on variations of the geometries that facilite so-called streamer reiginition is presented in section 3.4. The effects of dielectric permittivity are given in section 3.5, and a 3D simulation is presented in section 3.6. The paper is then closed with a few concluding remarks in section 4.

Simulation model
We use a three-species model for discharges in air, consisting of electrons, positive ions, and negative ions. Our model uses anÎto-KMC formulation of the drift-diffusion model for streamers, where we replace the conventionally used macroscopic advection-diffusion model for the electrons by an equivalent microscopic advection-diffusion model. The transport equation for the electrons in theÎto-KMC model is where X(t) is the electron position, V is the electron drift velocity, and D is the electron diffusion coefficient. N is a normal distribution with mean value zero and a standard deviation of one. We close the velocity relation in the local field approximation (LFA) and set the electron velocity and diffusion coefficients as where v e = −µ e E is the conventional fluid drift velocity (µ e is the electron mobility), and D e is the conventional fluid diffusion coefficient. Note that this is simply a microscopic representation of the conventionally used macroscopic advectiondiffusion model for the electrons, and that statistical averaging over identical particles yields standard drift-diffusion statistics. In essence, we are simply sampling the macroscopic evolution of the plasma using computational particles that represent average electrons. Deposition of particles on the mesh, and interpolation of the macroscopic quantities v and D to the particle position X are done using a cloud-in-cell scheme. Details regarding the underlying algorithm are given in [15]. A standard drift-diffusion model is used for positive and negative ions, whose densities are indicated by n ± , respectively. The equation of motion for the ions is where v ± , D ± , and S ± are the drift velocities, diffusion coefficients, and source terms for positive (+) and negative (-) ions, respectively. Charge conservation on dielectric surfaces also dictates that the surface charge density σ is locally conserved by where J σ is the charge flux onto the surface. Here, this charge flux is composed of impinging electron particles (which are deposited directly onto σ) and ion drift. Finally, the electric field E = −∇Φ is obtained by solving the Poisson equation (6) for the potential Φ: Here, ρ is the space charge density and ϵ r is the relative permittivity.
Our model differs from a standard streamer model in the usage of anÎto model for the electrons, and a KMC model for the plasma chemistry (discussed in the next section). The primary advantage of using theÎto-KMC formulation vis-avis a fluid formulation lies in the underlying discretization of the electron transport equation. With a fluid model, there is a fundamental requirement on the grid spacing ∆x [15,16] and a Courant-Friedrichs-Lewy (CFL) condition on the time step. However, theÎto-KMC formulation has no such restrictions, permitting us to use large computational time steps even for high-resolution grids.

Plasma chemistry
In addition to the transport and field equations above, we consider stochastic evolution of the number of particles in each grid cell by using the KMC algorithm as proposed by Cao et al [17]. Letting ⃗ X = (X e , X + , X − ) ⊺ denote the physical number of electrons X e , positive ions X + , and negative ions X − in a grid cell, we advance over a time step ∆t as where ⃗ R is a set of reactions and ∆t is a time step. This is only a symbolic representation of the integration algorithm; the full algorithm is substantially more complicated [15,17]. The Table 1. Simplified air plasma chemistry used for the example simulations. The notation ∅ indicates an untracked species (N 2 or O 2 depending on context) incorporated directly into the rate constant for the reaction. As we only sample ionizing photons, the photoionizing reaction γ + ∅ → e + M + is deterministic; each photon absorbed in the gas generates one electron-ion pair.

Reaction
Fluid rate KMC propensity Reference [22,23] γ + ∅ → e + M + -- [22] number of electrons X e (t) in a grid cell is immediately available since they are represented by computational particles. However, the positive and negative ions are represented as densities and so the number of these particles in a grid cell with volume ∆V is taken as where the lower brackets ⌊ and ⌋ round downwards (to a nonnegative integer). After running the KMC algorithm for the new number of particles in each grid cell, the reactive products are then added or subtracted directly from the grid cells. The remainder from equation (8) is, of course, kept during this process.
The plasma kinetics that we use is summarized in table 1. The electron diffusion coefficient D e , mobility µ e , temperature T e , ionization frequency k α , and attachment frequency k η are field-dependent and are computed by using BOLSIG+ [18] and the SIGLO database [19]. The ion mobility is set to 2 × 10 −4 m 2 (V s) −1 , and the electron-ion and ion-ion recombination rates are where T = 300 K is the gas temperature and T e = T e (E) is the electron temperature.
There is a subtle difference in the way the KMC algorithm and deterministic reaction rate equations operate, as the former uses chemical propensities a r ( ⃗ X) rather than deterministic reaction rates k. The propensities are defined such that a r ( ⃗ X)dt is the probability of exactly one reaction of type r occurring in the time interval dt. For example, the propensity for the reaction e + ∅ cα − → e + e + M + is where c α is the KMC rate for the reaction. Observe that we use symbols k for fluid rates and c for KMC rates. These are different rates, but for unipolar reactions one may show that k = c, i.e. the KMC rate and fluid rates are numerically equivalent [20]. However, for bi-particle reactions of the type e + M + kep − → ∅, observe that the reaction rate equation for this reaction would be ∂ t n e = k ep n e n + . (12) For the KMC algorithm, one may also derive the reaction rate equation in the limit X e , X + ≫ 1 as [20] ∂ t X e = a ep , where a ep = c ep X e X + is the propensity for the reaction e + M + cep − → ∅. One may realize that c ep cannot be a constant, since this would imply that particle reaction probabilities are independent of the distance between the particles. Since n e = X e /∆V, we find and thus we can associate the KMC rate c ep = k ep /∆V. While this scaling of the KMC rate by the grid volume may seem unconventional from a fluid point of view, it is microscopically correct as particles in a small grid volume are more likely to react than particles in a large grid volume.

Photoionization
We use the Zheleznyak photoionization model [22] including the corrections by Pancheshnyi [23] for modeling photon transport for the reaction e + ∅ kγ − → e + γ + ∅. This is a lumped representation of impact excitation of molecular nitrogen to an excited state N * 2 , i.e. e + N 2 → e + N * 2 , followed by spontaneous emission N * 2 → N 2 + γ and photoionization of molecular oxygen, i.e. γ + O 2 → O + 2 + e. The lumped rate constant is where ν Z (E) is a lumped function that accounts for excitation efficiencies and photoionization probabilities [23], while the factor p q /(p q + p) accounts for collisional quenching (i.e. non-radiative de-excitation) of N * 2 . The quenching pressure is p q = 40 mbar and the gas pressure is p = 1 bar. When a photon is generated within the reaction step we draw a random absorption coefficient for each photon as where 059 PHz, and f is a random number sampled from a uniform distribution on the interval [f 1 , f 2 ]. Note that this is a lumped photoionization model for photoionization of O 2 in the spectral range 980 Å to 1025 Å. The values of K 1 and K 2 correspond to the minimum and maximum absorption coefficients of O 2 in this range [24]. More elaborate models that also take into account the excited states of N 2 are also possible [14,25], but experiments [26] indicate that the simplified model above is sufficiently accurate for describing streamer discharges in air. After κ f has been determined for a computational photon, its propagation distance is determined by drawing a random number from an exponential distribution with parameter κ f .

Particle and photon interaction with boundaries
Our geometries are described with a signed distance field, given by a function S (x) that provides the distance to the geometry, where the material boundary is given by S (x) = 0. As the particle solvers (both electrons and photons) do not have CFL conditions, all intersection tests between particles/photons and the geometry are done using a ray-tracing algorithm on S (x), see figure 1. For example, for a photon whose emission and absorption positions are Y 0 and Y 1 , respectively, we first compute the distance to the geometry at Y 0 as and we then move the photon a distance d 0 along the photon path, This process is then recursed until either the photon strikes a boundary or reaches its final position Y 1 . The ray-tracing algorithm is conceptually simple and provides a CFL-free algorithm where photons can be propagated independent of the numerical grid. However, it also relies on a fast representation of the signed distance field. Fortunately, the simulations in this paper use very simple geometries that do not require surface tessellations, so these tests are typically very fast. We ignore photon scattering in the gas, but remark that inclusion of scattering only requires shallow modifications to the underlying ray-tracing algorithm.

Geometries
Our computational study is inspired by the experiments in Meyer et al [8], where we considered a pizza-slicer electrode and a dielectric plate with relative permittivity ϵ r = 3. In the experiments, the geometry was chosen this way so that the streamer inception region was well-defined while the streamer simulations could be done in Cartesian 2D coordinates. The full 3D geometry with indicated field stress (in arbitrary units) is shown in figure 2. For boundary conditions, we used a homogeneous Neumann boundary condition on the top y-face and homogeneous Dirichlet boundary (i.e. grounded potential) everywhere else. Full 3D streamer simulations at reasonable mesh resolutions end up using billions of grid cells and particles, and thus put a large strain on computational hardware. For simplicity, we invoke the same simplification as in Meyer et al [7,8] and reduce the simulations to Cartesian 2D coordinates. The fundamental approximation is then that the plasma has infinite extension along the z-axis indicated in figure 2. More precisely, we investigate the x-directed propagation of the streamer(s) that initiate right below the electrode disk, and assume that these streamers do not have substantial curvature along the z-coordinate. Although we have been unwilling to lavish vast computational resources on fully 3D simulations, we provide further justification for this simplification in section 3.6 where we present a fully 3D simulation. Figure 3 shows the full computational domain, including a comparison between the corresponding 2D and 3D field calculations, evaluated at z = 0 (the center of the 3D geometry). For simplicity, the 2D calculations also ignore the wheel holder. We find that the error introduced by the 2D reduction of the field is less than 2% on the center axis, and this error arises largely due to neglection of the curvature of the wheel in the z-direction. Further away from the symmetry plane, i.e. moving along the z-coordinate in figure 2, this error increases substantially. Overall, the larger source of 2D-versus-3D discrepancy is probably due to the assumption of translational invariance of the plasma along the z-axis. Experiments indicate that this is only qualitatively true [7,8].
We consider four basic sizes of the square channels on the dielectric surface. The depths and widths of each channel are identical, and range from w = 250 µm to w = 2 mm. A flat surface is used as a reference surface for these calculations. Except for the investigations in section 3.4, the spatial separation ∆ between the channels is always one channel width, i.e. ∆ = w. Figure 4 show the field distribution near the inception region for the various surface profiles.

Initial conditions
For initial conditions we drew ten computational electrons randomly located in a circle with a 500 µm radius centered on the electrode tip. Electrons whose stochastically sampled positions ended up inside the electrode were removed before the simulation started, so that the simulations started from 7 electrons. We point out that, strictly speaking, the computational electrons are lines in 2D Cartesian coordinates.
In experiments and simulations, streamers only appear if there is a starting electron in a region where there is a sufficiently high field over a large enough distance to cause an avalanche-to-streamer transition. The streamer morphology is in general dependent both on the applied voltage pulse, which varies in time, but also where on the voltage pulse the starting electron appears, i.e. on the so-called statistical time lag. Here, we wish to lock our study into fewer degrees of freedom; the peak voltage and the surface morphology, and we therefore dispense with the extra complications and degrees of freedom associated with transiently varying pulses and statistical time lags. Our simulations therefore used step voltages of U ∈ 15 kV-30 kV in steps of 5 kV.

Computational framework and numerical parameters
We run 20 baseline computer simulations-four different voltages for each of the five surfaces. These baseline simulations are then supplemented by additional simulations where we vary the dielectric permittivity, turn on/off transparent boundaries, change the channel spacing, etc. This is done in order to shed additional light on emerging transformations of the underlying streamer behavior for variations in the geometric and physical models. All calculations are done using the chombo-discharge [27] computer code, which is publicly available.
We use Cartesian adaptive mesh refinement (AMR) in our simulations (see inset in figure 2), refining grid cells if and coarsening if where ∆x is the grid resolution and α and η are the Townsend ionization and attachment coefficients, respectively. The simulations used adaptive grids with a finest resolutions of ∆x ≈ 2.44 µm. We used constant-sized time steps of ∆t = 25 ps for all simulations, and integrated for 300 ns or until the streamer hit the grounded electrode (i.e. side walls). The computer simulations typically ran in 1-6 h on 512 CPU cores. We could probably have run the calculations on fewer resources since there were not always enough work for the various MPI ranks, but 512 CPU cores is the minimum allocation on the supercomputer that we used.

Travel curves
We extracted travel curves for the streamer head, measuring its left-and right-most position along the x-axis for all simulations. We used the position of the maximum value of the electron ionization source term as a proxy for the streamer head, and we ignored contributions from backward traveling streamers inside the channel spaces (which appeared in some of the simulations). Figure 5 shows the streamer head position for the various simulations, focusing on a quantitative velocity comparison between different surfaces for a fixed applied voltage. We also include the surface profiles on the right-hand side of this figure in order to more readily associate trends in the travel curves with local morphology changes in the insulator surface. We observe the following features: • The mean streamer velocity increases with voltage for all surfaces. • The mean streamer velocity decreases as the channel dimensions increase.
From these trends we conclude that the patterned surfaces result in substantially slower streamers, which is in agreement with past work [8]. It is hardly surprising that increasing the voltage also increases the streamer velocity, since an increased voltage is equivalent to increasing the background electric field. However, we again point out that we are not including a study of the statistical time lag in these simulations. In earlier work [8] we showed that the streamer velocity is also dependent on this factor, which is not studied further in this paper.

Morphology
In previous work [7,8] we have shown that the streamer morphology is substantially dependent on the local surface morphology. Figure 6 shows snapshots of the plasma density for the four profiled surfaces at various time instants for the simulation with U = 25 kV, corresponding to the travel curves in figure 5(c). The frames in figure 6 are separated by 2.5 ns, with the exception of the first and second frame for w = 2 mm since the streamer takes a long time climbing out of the channel. As in previous work [8] we find that the propagation is at least partially influenced by reignition in subsequent channels. This is particularly clear for w = 2 mm where a secondary streamer appears in the next channel (see frame with t = 47.5 ns) before the primary streamer has crossed into it. However, for the surface with the smallest channel spacing w = 250 µm, the streamer tends to hop over channels and no streamer reiginition processes could be identified.
It is well-known that positive streamers propagating over dielectric surfaces lead to the formation of a cathode sheath on the dielectric [5,6]. The sheath is completely analogous to the sheath formation on metallic cathodes, and arises because low-energy electrons travel outwards from the surface, leaving a surplus of positive ions hovering over it. The field in the sheath can become quite high, ranging between 30 kV mm −1 -60 kV mm −1 in our simulations for the various surfaces. Unfortunately, we cannot study this sheath in detail because of the inherent limitations in the LFA. Within this approximation, secondary electrons that appear in the sheath either due to photoionization or electron impact ionization are born with energies given by the LFA closure, i.e. the ionization coefficient is given parametrically as a function of E. Since the field in the sheath is quite high, the LFA artificially raises the energies of the secondary electrons in the sheath, and ultimately the LFA-based model predicts an artificially high ionization rate in these regions.

Screening of ionizing radiation
Positive streamers rely on free electrons ahead of them, and in air these are supplied through photoionization of molecular oxygen. As electron avalanches grow toward the streamer head, they extend the positive streamer channel forward and during this process they also generate new seed electrons. Due to the reliance on this source of electrons, it is reasonable to expect that geometric screening of ionizing radiation affects positive streamer propagation. Figure 7 shows the photoionization profile during streamer inception for one of the computer simulations, showing clear-cut shadows in the channel spaces. If we were to use a continuum approximation of the photoionization process, as in Marska [5] and Bourdon et al [12], this feature is much less pronounced even when appropriate boundary conditions [5,28] are used. The underlying issue in the isotropic approximation is that the radiative flux F γ is proportional to the gradient of the radiative intensity Ψ, i.e. F γ ∝ ∇Ψ, which will always lead to artificial diffusion of ionizing radiation around geometric obstacles. Helmholtz reconstructions [12,29] of the photoionization profile are quite popular in the streamer simulation community. However, the underlying equations that underpin the Helmholtz reconstruction are equivalent to the Eddington equations, a first order closure of the radiative transfer equation which is strictly speaking only valid for highly scattering media [28]. Since the Helmholtz and Eddington reconstructions have no directional dependence on the equations of motion, they fail to accurately capture shadows. To the best of our knowledge, higher order closures or usage of discrete ordinates have not been investigated in the context of streamer discharges. Direct solutions to the RTE have been reported by Capeillère et al [30], but the method requires solutions in higher-order space (six computational coordinates for three-dimensional physical space).
To test the importance radiative shielding, we next run computer simulations with geometries that are completely transparent to the ionizing radiation. In practice, we simply turned off the ray-tracing photon interaction with the geometry to achieve this functionality. We select the computer simulation with w = 1 mm and U = 25 kV, and compare the results with the original computer simulation where the dielectric is completely opaque. Figure 8 shows snapshots of the plasma density, comparing the temporal evolution for the opaque and transparent geometries. The panels on the left hand side of this figure show the results for the opaque geometry, and the panels on the right hand side show the corresponding results with the transparent geometry. Sub-figures (a)-(g) show the plasma density on a logarithmic color scale at various time instants indicated in each panel. Observe that in e.g. figure 8(c) there is some plasma in the channel spaces in the right hand side panels at x ∼ 6.5 mm and x ∼ 8 mm. However, the corresponding channels in the left hand side panels are geometrically screened and thus do not contain this artificially raised plasma density. The evolution in the two models differ as to when the reignition in the next channel takes place. In figure 8, reiginition for the opaque geometry is barely discernible on the right-hand side of the channel just before the primary streamer descends into the channel space. This can be seen in the left panels in figures 8(d) and (f) as a small blob of plasma on the right-hand side of the channel walls. For the transparent geometry, streamer re-ignition in the channels is more discernible, and it also occurs earlier. Re-ignition of this kind can be seen in the right-hand side panel in figure 8(d), where almost the entire channel pore space is filled by plasma before the primary streamer has climbed over the profile. A similar process can be see in figure 8(f) (right column). As the only difference between the two simulations is due to geometric transparency of ionizing radiation, we conclude that the premature reignition process is simply due to artificial accumulation of seed electrons in the neighboring channel space.
Running the simulations in figure 8 further, we observe significant deviations in both morphology, mean velocity, and thus ultimately streamer range. In the computer simulations with correct screening of the ionizing radiation, the streamer tends to traverse the channel gap before later filling it with plasma (see section 3.2). With artificial transparency turned on, the channel space is first filled by a reignited streamer, and the primary streamer connects to the tail of secondary streamer. The cumulative outcome is that the premature reignition process leads to a deeper penetration of the primary streamer into the channel, and it will consequently also have a harder time climbing out of it. The travel curve for the transparent geometry is given in figure 9(a), comparing it with the corresponding data from figure 5. For the transparent geometry, the streamer propagates with approximately 2/3 of the velocity of the opaque streamer simulation. Figures 9(b) and (c) show snapshots of the plasma density on a linear scale for this simulation at t = 100 ns, where the morphology changes are clearly discernible.

Role of channel spacing
As observed in Meyer et al [8], streamers may in principle reignite in the neighboring channel if there is at least one starting electron there, and the field is sufficiently high over a long enough distance. The first requirement requires photoionization directly into the channel, and in Meyer et al [8] we showed that this is feasible. The second requirement depends on the distance between the channels: As the distance between the channels increase the field in the next channel space is reduced. Vice versa, decreasing the distance between the channels leads to an increase in the field in the channels. Likewise, since the intensity of ionizing radiation drops off as 1/r in Cartesian 2D, the number of ionizing photons in the channel spaces will be larger for narrowly spaced channels. Consequently, one can expect that reignition processes become increasingly relevant as the distance between the channels is reduced.
Next, we run simulations with w = 2 mm and U = 25 kV, but vary the distance between the channels from ∆ = 0.5 mm to 2 mm in steps of 500 µm. This corresponds to varying the number of channels on the surface from 17 to 27. Figure 10 shows snapshots of the plasma density for the four surfaces after t = 100 ns. The data is shown on a truncated logarithmic color scale for improved visual exposition of the streamer paths. In particular, as we decrease ∆ we find that reignition in the neighboring channel occurs earlier, which is what we expected based on the simple electrostatic and geometric evaluations discussed above. This can be seen in figure 10(d) where the channel centered at x = 12.5 mm contains a streamer, but where the primary streamer in the preceding channel centered at x = 10 mm has not yet climbed out of it. In figure 10(c) the channel spacing is twice that of figure 10(d), but no reignition can be observed despite the fact that the primary streamer has climbed further out of the channel than the corresponding streamer in figure 10(d). Figure 11 shows the streamer travel curves when varying ∆. For reference, we have also included the flat surface. The data for the flat surface and for the surface with ∆ = 2 mm correspond to the data given in figure 5(c). For the surface with the narrowest investigated channel spacing ∆ = 500 µm the discharge propagation slows down substantially, and the discharge terminates after x ∼ 15 mm. We believe that the streamer termination for this surface occurs due to (1) the increased number of channels that the streamer needs to climb out of, and (2) the reignited streamer in the neighboring channel. Figure 12 shows the plasma density and potential distribution during the reignition process. The reignited streamer seen in figure 12 at t = 50 ns is a double-headed streamer. It has a negative streamer tail which propagates backwards towards the anode, and a positive streamer head that propagates toward the grounded side. The primary streamer head is electrically attracted to the secondary streamer tail since they carry opposite charge polarities, and we believe that the slow upwards propagation for this case is due to the fact that the attraction between the primary streamer head and the secondary streamer tail occurs directly through the dielectric. They are also electrically disconnected and there is a comparatively large potential difference between them, approximately 4 kV in figure 12 at t = 60 ns. When the streamers finally connect after t ∼ 90 ns this potential difference is transferred to the secondary streamer head, which then becomes the primary streamer.

The role of dielectric permittivity
Next, we consider three additional simulations with relative dielectric permittivities ϵ r = 1, ϵ r = 6, and ϵ r = 9. As in section 3.4, we select the simulation with w = 2 mm and U = 25 kV as the baseline geometry and applied voltage. Figure 13 shows the plasma density after t = 200 ns for the four simulations, including the baseline case with ϵ r = 3. We find that with ϵ r = 1 the streamer does not propagate over the surface at all, whereas the streamer velocity and ranges otherwise increase with increasing ε r . The result is comparatively straightforward to explain since increasing the permittivity of the dielectric is tantamount to increasing its polarization. Since the polarization decreases the electric field in the dielectric and increases it in the gas phase, it leads to a de-facto increase in the amount of electron impact ionization in the gas, and thus a faster streamer. In previous work, Li et al [6] showed a very low streamer velocity dependence on the relative dielectric permittivity for positive streamers propagating over planar surfaces. The difference with that work, however, is that the dielectric is oriented parallel to the discharge gap, and there is therefore no permittivity-induced modification of the background electric field in the gas. Here, the dielectric is orthogonally oriented with respect to the background electric field, and we therefore obtain a different polarizability effect from the dielectric. Figure 14 shows the travel curves for the simulations in figure 13 when we vary the dielectric permittivity. We find that the velocity of the streamers increase with increasing permittivity, which we attribute to the fact that increasing the polarization of the dielectric increases the electric field in the gas phase.

3D simulation
In previous sections we restricted our computer simulations to Cartesian 2D simulations. While the background electric field on the investigated computational axis had a low error in the 2D calculations (see figure 3), we did not judiciously demonstrate that this error remains sufficiently small also when the streamers start propagating. This is a point worthy of investigation, especially since branching is a fundamental capability of streamers. Filamentation of the surface discharge implies one extra dimension of curvature in the space charge layer of the streamer, and thus additional field amplification at the discharge tip(s). To investigate these effects, we consider a fully 3D simulation with w = 500 µm and U = 25 kV. The computational domain and boundary conditions are the same as in figure 2, and measure 8 cm × 8 cm × 32 cm. We also use the same procedure for sampling the initial electrons, drawing 10 electrons in a sphere with a 500 µm radius centered at the disk electrode edge that is closest to the dielectric. Unlike the 2D simulations where these electrons are idealized line charges, they represent physical electrons in the 3D simulation.
We have restricted the finest AMR grid resolution to ∆x ≈ 20 µm, which is eight times larger than the grid resolution in the corresponding Cartesian 2D simulation. The effective domain is therefore 4096 × 4096 × 16 384 cells. The 3D computer simulation was run on 8192 CPU cores for about 12 h, and used around 500 million grid cells and 4 billion computational particles. Unfortunately, since the 3D computer simulations strain our available computational hardware, we did not perform this analysis for all the corresponding 2D simulation cases. Figure 15 shows an isosurface of the electron density n e = 10 18 m −3 after t ≈ 36 ns. In the simulation we make several observations: (i) There are additional discharges occurring further up on the disk electrode, partially in agreement with experimental observations of a similar setup [8]. The experiments performed by Meyer et al [8] used a similar geometry, but high-speed imaging also revealed initiating discharges on the full circumference of the disk electrode, but in our the 3D simulation we only observe such discharges on the lower part of the electrode. However, the voltage in the experiments [8] used a modified lightning impulse with a peak voltage up to 45 kV, while the voltage in our simulations is a constant step-voltage of 25 kV, which is presumably why these extra discharges are not seen in our computer simulation. (ii) There is some filamentation of the surface discharge, but (1) there are several side branches and (2) the curvature in the z-direction is smaller than the curvature in the xy-plane, see inset in figure 15(a). In other words, the surface filaments do not appear as three-dimensional filaments in the same way that streamers in bulk air do. Instead, they are wider than they are thick. This implies that the self-enhanced field of the discharge is mostly due to xy-curvature in the space charge layer. Since this is also the part of the discharge that we studied in the 2D simulations, we conclude that the 2D approximation is qualitatively valid for this part of the streamer. (iii) Clearly, there are numerous features in the full 3D simulation that can not be predicted by our 2D simulations, which are unable to describe any z-variations in the discharge. For example, the longitudinal spread of the discharge (along z) cannot be predicted by the 2D simulations. We point out that there is also lateral spread inside the channels, as seen in the inset in figure 15(a). To provide additional insight into the quantitative relevance of the 2D simulations, we perform an apples-to-apples comparison between the 2D and 3D simulation. Figure 16 shows the plasma density on a truncated logarithmic scale after t ≈ 36 ns. Here, the left portion of the figure with x ⩽ 0 shows the plasma density as it was calculated in the Cartesian 2D simulations. The right-hand side of the figure shows the corresponding plasma density from the 3D simulation, sliced through z = 0. In general, the morphological features that we observed in the 2D simulations are also observed in the 3D simulations. For example, there is little or no plasma in the channels ahead of the streamer, which is due to geometric screening of ionizing radiation. The 3D streamer simulation is approximately 25% slower than the corresponding 2D simulation. We have not fully investigated the cause of this discrepancy, but remark that the spatial resolution in the 3D simulation is substantially coarser than in the 2D simulation. Numerical convergence is consequently one of the factors contributing to this discrepancy. Other factors are e.g. due to differences in streamer initiation morphologies, filamentation, and different contributions of stochastic noise due to both particle and photoionization processes (which are mostly suppressed in the 2D simulations). In general, the 2D simulations reproduce the pertinent features of the 3D discharge quite well. We conclude that the 2D approximation is valid for the transverse spread of the discharge, which explains the good quantitative agreement observed between 2D computer simulations and experiments in Meyer et al [8].

Summary
We have investigated the propagation of positive streamers in air over dielectric surfaces with a modified surface morphology consisting of rectangular channels running orthogonal to the propagation direction. The calculations were primarily done in Cartesian 2D, using a model that also incorporates geometric screening of ionizing radiation. Four basic surfaces were investigated, and an idealized flat surface was used as a reference surface. The results can be summarized as follows: • Larger channels with narrow spacings are more efficient for restricting the propagation length of positive streamers in air. • Geometric screening of ionizing radiation can not be ignored. By letting the geometry be artificially transparent to radiation, the streamer morphology and mean velocity changed significantly. • Decreasing the channel spacing led to slower streamers.
• Increasing the dielectric permittivity increases the streamer velocity.
We substantiated the two-dimensional simplification by presenting a fully 3D simulation of one of the 2D simulation cases. Although we found that 3D simulations showed some degree of filamentation, a direct comparison with the 2D simulations showed the transverse spread of the discharge was quantitatively well approximated by the two-dimensional approximation.
The numerical results in this paper show that positive streamers propagating over dielectric surfaces can be substantially slowed down by engineering the surface pattern, thus inhibiting surface discharges. Our results have substantial value for high-voltage technology where streamer termination is a desirable outcome. However, there are still some open questions left behind in our study. For example, we have not investigated the role of pre-ionization or chemical composition of the plasma (we use a simplified reaction set). Secondary electron emission from the dielectric, and negative streamers, have not been studied. The present study focus quite a bit on the role of ionizing radiation in the presence of obstructive boundaries. As negative streamers do not require photoionization ahead of them, several major questions are left open regarding the interaction between negative streamers and such complex surfaces.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors. https://github.com/chombodischarge/chombo-discharge/.