Accurate and efficient calculation of photoionization in streamer discharges using fast multipole method
aa r X i v : . [ phy s i c s . c o m p - ph ] J un Accurate and efficient calculation of photoionizationin streamer discharges using the fast multipolemethod
Bo Lin
Department of Mathematics, National University of Singapore, 10 Lower Kent RidgeRoad, Singapore 119076E-mail: [email protected]
Chijie Zhuang
State Key Lab of Power Systems and Department of Electrical Engineering,Tsinghua University, Beijing 100084, ChinaE-mail: [email protected]
Zhenning Cai
Department of Mathematics, National University of Singapore, 10 Lower Kent RidgeRoad, Singapore 119076E-mail: [email protected]
Rong Zeng
State Key Lab of Power Systems and Department of Electrical Engineering,Tsinghua University, Beijing 100084, ChinaE-mail: [email protected]
Weizhu Bao
Department of Mathematics, National University of Singapore, 10 Lower Kent RidgeRoad, Singapore 119076E-mail: [email protected]
Abstract.
This paper focuses on fast and accurate three-dimensional simulation ofthe photoionization in streamer discharges based on the classical integral model derivedby Zheleznyak et al . The simulation is based on the integral form of the photoionizationrate with the kernel-independent fast multipole method. The accuracy of this methodis studied quantitatively for different domains and various pressures in comparison withother existing models based on partial differential equations (PDEs). The comparisonindicates the numerical error of the fast multipole method is much smaller than thoseof other PDE-based methods, with the reference solution given by direct numericalintegration. Such accuracy can be achieved with affordable computational cost, andits performance in both efficiency and accuracy is quite stable for different domains and alculation of photoionization in streamer discharges by FMM pressures. Meanwhile, the simulation accelerated by the fast multipole method exhibitsgood scalability using up to 1280 cores, which shows its capability of three-dimensionalsimulations using parallel (distributed) computing. The difference of the proposedmethod and other efficient approximations are also studied in a three-dimensionaldynamic problem where two streamers interact.PACS numbers: 02.60.Cb, 02.70.-c, 52.80.-s Keywords : streamer discharge, parallel computing, fast multipole method (FMM),photoionization
1. Introduction
As a natural phenomenon of non-thermal filamentary discharge with a large amountof applications, streamer discharge happens when an insulating medium such as air isexposed to a sufficiently strong electric field, where electron avalanche occurs and formsfilamentary streamers. The filamentary streamer discharges are pivotal for many gasdischarges in nature [14], including the frequently seen lightning [36] and the spritedischarges in high altitude [25, 30]. It has mature industrial applications [5, 2, 37]including dust precipitator, ozone production, and water purification [40, 19]. A reviewof streamer discharge and its role in these industrial applications can be found in [13].Streamers can be classified into positive ones and negative ones.The photoionizationplays an important role in the propagation of streamers in air, especially for positiveones. In particular, the photoionization provides seed electrons ahead of the tips, whichare required by the propagation of positive streamers [43, 20, 34, 48, 49]. Besides, thestochastic photoionization is shown to have an impact on the branching of streamer[42, 3, 29].The photoionization is important in streamer discharges; its modelling andsimulation have attracted continuous attention. The classical model for oxygen-nitrogenmixture derived by Zheleznyak et al. in [46] is widely utilized in the simulation of positivestreamers [32, 39], and was improved in [33, 18] to gain better accuracy and has beenextended to a stochastic version in [9].Direct calculation of the classical integral model requires a large amount ofcomputation, especially in three dimensions (3D) where streamer discharges inherentlyhappen. To ease the numerical difficulty and reduce the computational cost, someapproximation methods are proposed in [6, 24, 27] based on the kernel expansion andconversion to Helmholtz equations. Moreover, modeling of photoionization based on theradiative transfer equation (RTE) also provides good results [7].Less than two decades ago, the kernel-independent fast multipole method (FMM)was proposed to compute particle interactions efficiently and accurately [45, 44]. It canbe easily applied to the convolutional integrals [28], and its computational complexity iscomparable to the method of fast Fourier transform (FFT). Compared with FFT, FMM alculation of photoionization in streamer discharges by FMM
2. Model formulation
To make the contents self-contained, we briefly review commonly used approaches forphotoionization calculations.
The widely used photoionization model derived by Zheleznyak et al. [46] describes thephotoionization rate by S ph ( ~x ) = Z Z Z V ′ I ( ~y ) g ( | ~x − ~y | )4 π | ~x − ~y | d ~y, ∀ ~x ∈ V, (1)where ~x = ( x, y, z ) T , V ′ is the source chamber in which the photons are emitted, and V is the collector chamber where the photons are absorbed, I ( ~y ) is proportional to theintensity of the source radiation: I ( ~y ) = ξ p q p + p q ωα S i ( ~y ) , (2)where ξ is the photoionization efficiency, p q is the quenching pressure, p is the gaspressure, ω and α are the excitation coefficient of emitting states without quenchingprocesses and the effective Townsend ionization coefficient, respectively, with ωα being acoefficient to be determined by experiments, and S i is the effective ionization rate. Thefunction g ( r ) = g ( | ~x − ~y | ) in (1) is given by g ( r ) p O = exp( − χ min p O r ) − exp( − χ max p O r ) p O r ln( χ max /χ min ) , (3)where r = | ~x − ~y | , p O is the partial pressure of oxygen, χ max = 2 cm − Torr − and χ min = 0 .
035 cm − Torr − are the maximum and minimum absorption coefficients of O alculation of photoionization in streamer discharges by FMM g ( r ) in (3), we follow [46, 6, 33] to write g ( r ) /p O on the left-hand side so that theright-hand side is dependent directly on the product p O r . Interested readers may referto [46] and [33] for more details.Clearly, Eq. (1) is a convolution in three dimensions. A naive numericalimplementation of (1) requires a whole domain quadrature for every point ~x ∈ V ,which requires a time complexity O ( N ) with N being the total number of degrees offreedom. One idea to reduce the computational cost is to use a coarse grid in the weakfield at the price of possibly losing some accuracy, as [20] did in the three dimensionalcases with cylindrical symmetry. Instead of a straightforward computation of the integral, the efficiency can besignificantly enhanced by converting it into a problem of differential equations at theexpense of losing some accuracy. One important and pioneer work was done in [27],which approximates the photoionization kernel as the sum of the fundamental solutionsof a number of partial differential equations. The function g ( r ) defined by (3) isapproximated as follows: g ( r ) p O ≈ p O r N E X j =1 C j exp( − λ j p O r ) , (4)where λ j and C j (1 ≤ j ≤ N E ) are constants that can be fit numerically [27, 6].Consequently, it suffices to take the linear combination of S ph ,j to approximate theintegral (1) S ph ( ~x ) ≈ N E X j =1 C j S ph ,j ( ~x ) , (5)where S ph ,j ( ~x ) is the solution of the following modified Helmholtz equation( − ∆ + ( λ j p O ) ) S ph ,j ( ~x ) = ( p O ) I ( ~x ) . (6)The modified Helmholtz equation (6) can be solved efficiently by numerous fast ellipticsolvers like multigrid-preconditioned FGMRES method [23]. N E = 2 was used in [27], and the constants λ j and C j were chosen to fit the low-pressure experimental data from [35] (the misprint of these constants is corrected in [12]). N E = 3 was suggested in [6] for a better fitting for the range 1 < p O r <
150 Torr · cm,and the constants λ j and C j are chosen to fit the function in Eq. (4) since it agrees wellwith both the low-pressure experimental data from [35] and the experimental data inatmospheric air from [1], as indicated in [31]. While zero boundary conditions were usedin [27], it is suggested in [6] that the boundary condition for Eq. (6) can be providedby computing the integral (1). In this paper, we take the three-term exponentialapproximation and adopt the coefficients in [6] listed in Table 1. alculation of photoionization in streamer discharges by FMM Table 1.
Coefficients of three-exponential ( N E = 3) approximation in (4) [6]. j C j (cm − Torr − ) λ j (cm − Torr − )1 1 . × − . . . . . Another type of differential equations that can facilitate the computation of thephotoionization rate is the radiative transfer equation. In [38, 6, 7], the following multi-group approximation of the steady-state radiative transfer equation is chosen to describethe intensity of radiation Ψ j for the j -th group of spectral frequency: ~ω · ∇ Ψ j ( ~x, ~ω ) + κ j Ψ j ( ~x, ~ω ) = n u ( ~x )4 πc τ u , j = 0 , , · · · , N ν , (7)where ~ω ∈ S is the solid angle defined on the unit sphere, κ j is the absorption coefficient, n u is the density of the species with the excited state u , c is the speed of light and τ u is the radiative relaxation time for the state u . Here the scattering and the change infrequency of the photons during collisions with molecules have been neglected [38, 7].For photoionization in air, κ j = λ j p O , and for simplicity, only one excited state isconsidered n u ( ~x ) τ u = I ( ~x ) ξ , (8)with λ j to be determined by data fitting [46, 38, 6]. The photoionization rate is thenproportional to the weighted sum of the integral of Ψ j over ~ω ∈ S : S ph ( ~x ) = N ν X j =1 A j ξ p O c Z S Ψ j ( ~x, ~ω )d ~ω, = N ν X j =1 A j ξ p O c Z Z Z V n u ( ~y ) c τ u exp( − λ j p O | ~x − ~y | )4 π | ~x − ~y | d ~y, (9)where A j are also parameters which can be fit according to the experimental data, andwe have applied the analytical solution of the radiative transfer equation to express theright-hand side. To determine the parameters, it is noticed that (9) is identical to (1) if N ν X j =1 A j p O exp( − λ j p O r ) = g ( r ) , r = | ~x − ~y | , (10)where g ( r ) is given in (3), and the coefficient A j and λ j (1 ≤ j ≤ N ν ) are determinedby fitting the left hand side of (10) with g ( r ) in the range 0 . < p O r <
150 Torr · cm [6].The results for three-group ( N ν = 3) approximation are shown in Table 2. alculation of photoionization in streamer discharges by FMM Table 2.
Coefficients of three-group ( N ν = 3) approximation in (10) [6]. j A j (cm − Torr − ) λ j (cm − Torr − )1 0.0067 0.04472 0.0346 0.11213 0.3059 0.5994Instead of computing the integral in (9), a more efficient way to get the intensityfunction Ψ j is to solve (7) as an differential equation. For example, in [7], adirect solver of (7) is employed for two-dimensional axisymmetric discharges usingthe finite volume method for both space and angular variables. However, theradiative transfer equation (7) is still a five-dimensional partial differential equation.Further reduction of dimensionality can be realized by the improved Eddington or SP approximation [38, 6, 47]. In [21], the simplified P N ( SP N ) approximations ofoptically thick radiative heat transfer equations are theoretically derived by asymptoticanalysis. SP N approximations are introduced in [38] to obtain a fast numericalsimulation for the photoionization source term mainly with monochromatic (one-group)approximation. The SP N approximations for photoionization are further improved in[6], and extended to multi-group approximation, including the three-group SP methodwhich approximates the isotropic part of the solution by [38, 6] Z S Ψ j ( ~x, ~ω )d ~ω = γ φ j, ( ~x ) − γ φ j, ( ~x ) γ − γ , (11)where γ n = h − n q i with n = 1 ,
2, and φ j, ( ~x ) and φ j, ( ~x ) are solutions ofthe following two Helmholtz equations − ∆ + ( λ j p O ) µ ! φ j, ( ~x ) = λ j p O µ n u ( ~x ) c τ u , (12) − ∆ + ( λ j p O ) µ ! φ j, ( ~x ) = λ j p O µ n u ( ~x ) c τ u , (13)with the coefficients µ n = r + ( − n q ( n = 1 , ∇ φ j, · ~n + α ( λ j p O ) φ j, = − β ( λ j p O ) φ j, , (14) ∇ φ j, · ~n + α ( λ j p O ) φ j, = − β ( λ j p O ) φ j, , (15) alculation of photoionization in streamer discharges by FMM ~n is the outward unit normal vector, α n = (cid:16)
34 + ( − n − q (cid:17) and β n = (cid:16) − n q (cid:17) ( n = 1 ,
3. Fast multipole method for accurate and efficient evaluation of integral
As can be seen from the Sections 2.2 and 2.3, different methods based on differentialequations have been proposed to approximate the integral (1) or (9), leading to muchhigher numerical efficiency. However, these methods may suffer numerical issues, i.e.the approximation errors might be significant in some cases. On the other hand, despitethe high computational cost[12], the results calculated from the integral form are freeof further approximations, therefore, these results are often used as reference solutions[6, 38, 8]. Moreover, the integral form can be easily extended to stochastic versions[9, 41]. The importance of the integral form inspires us to tackle the original integrationproblem (1) directly using fast algorithms. The exponential decay of the kernel withrespect to the distance (see (3)) reminds us to adopt the fast and accurate fast multipolemethod [44, 45], which utilizes the low-rank structure of far-away interactions to gainsignificant speed-up.The fast multipole method used in this paper is established based on the fastevaluation of the numerical quadrature of (1). For convenience, we discretize S ph and n e on the same mesh. In general, the integral (1) can be discretized as S ph ( ~x i ) = N pt X j =1 G ( ~x i , ~y j ) I ( ~y j ) , i = 1 , · · · , N pt , (16)where G ( · , · ) is the discrete kernel function calculated from the corresponding functionin (1) and the numerical quadrature weights. In this paper, we apply the midpointquadrature rule on each grid cell unless ~x i and ~y j are in the same grid where the second-order Gauss-Legendre quadrature rule is alternatively applied. We remark that this isnot essential and other numerical quadratures can also be used.The kernel-independent adaptive fast multipole method [44] does not require theimplementation of multipole expansions [15, 10] of the kernel function. Based on ahierarchical tree, it uses a continuous equivalent density on a surface enclosing a box torepresent the potential generated by sources inside the box. Given a set of N pt pointsin three dimensions, a hierarchical octree is constructed adaptively such that each leafcube of the tree contains no more than m points, where m is a selected constant. Thisoctree can be built from a sufficiently large root cube to contain all N pt points, and thensubdivided to equal-sized sub-cubes recursively if the current cube contains more than m points. For illustrative purpose, an example of the hierarchical tree in two dimensions(2D), i.e. quadtree, is shown in Figure 1.To sketch the idea, we consider the simple case where the source points are uniformlydistributed. This corresponds to the case when the uniform mesh is applied in the alculation of photoionization in streamer discharges by FMM + ++++ +++++ + ++++ +++++ + ++++ +++++ + ++++ +++++ Figure 1.
An example of a hierarchical tree in 2D, with N pt = 10, m = 2. The arrowsshow the construction procedure, and the circles with “+” denote the N pt points. B F ( B ) N ( B ) Figure 2.
Cross section of near range N ( B ) and far range F ( B ) of a box B in 3D.The blue thick side is the boundary of B , green part is N ( B ) and red part is F ( B ). discretization of I ( ~y ) in (2). In this case, for each target point ~x i in a cube or box B , fastmultipole method splits the summation (16) into two parts, namely, near interactionsand far interactions: S ph ( ~x i ) = X ~y j ∈N ( B ) G ( ~x i , ~y j ) I ( ~y j ) + X ~y j ∈F ( B ) G ( ~x i , ~y j ) I ( ~y j ) , (17)where N ( B ) and F ( B ) are the near range and far range of B , respectively. For ~y j ∈ N ( B ), the interactions with all ~x i ∈ B are calculated directly. For the points ~y j ∈ F ( B ), the interactions can be approximated with controlled accuracy due to thelow-rankness of G ( ~x i , ~y j ). If a box is centered at ~c with side length 2 r , then N ( B ) isdefined as a box centered at ~c with side length 6 r , and F ( B ) is the domain outside N ( B ) (See Figure 2).In (17), the summation for points ~y j ∈ F ( B ) can be approximated using thehierarchy tree. The idea is composed of two parts: 1) represent the potential generatedfrom source points inside any box B by some equivalent source points enclosing B ; 2)represent the potential generated from source points in F ( B ) by other equivalent sourcepoints enclosing B , which gives an approximation to the summation for points ~y j ∈ F ( B )in (17). The first part is implemented by post-order traversal of the hierarchical tree. If B is a leaf box, the potential generated from the source points inside B is represented byseveral equivalent points surrounding the box, as is called the multipole expansion to bedefined in (18). If B is not a leaf box, its multipole expansion can be accumulated from alculation of photoionization in streamer discharges by FMM B from original source points in F ( B ) by a small number of equivalentsource points in F ( B ) calculated from the first part. This is the idea of the second part,and we similarly represent the potential generated from source points in F ( B ) by someequivalent source points surrounding B , as is called the local expansion to be defined in(19). The second part is implemented by pre-order traversal of the hierarchical tree. Ifa non-root box B is embedded in its parent box P ( B ), its local expansion is calculatedfrom: the accumulation of the local expansion of P ( B ), which is called “L2L translation”to be defined in (22); and the multipole expansion of the boxes in N ( P ( B )) but notadjacent to B , as it is implemented by the operation called “M2L translation” to bedefined (21).We now show more details about the kernel-independent FMM: firstly introducemultipole expansion and local expansion in the FMM, and then show three translationsamong them: M2M (multipole to multipole), M2L (multipole to local) and L2L (localto local). For simplicity, we would like to neglect the vector symbol on ~x and ~y whenintroducing FMM. Multipole expansion
Multipole expansion of a box B is used to represent the potentialin F ( B ), generated by the source inside B . Two surfaces of the cube are introduced forthe approximation, upward equivalent surface y B,u and upward check surface x B,u . Theequivalent surface y B,u should be taken to enclose B , and check surface x B,u enclosesequivalent surface y B,u . Moreover, both y B,u and x B,u should locate inside N ( B ). Seethese two box surfaces in Figure 3. B +++ x B,u y B,u (1) (2) B +++ y B,d x B,d (1)(2)
Figure 3.
Cross section of equivalent surfaces and check surfaces in multipoleexpansion (left subfigure) and local expansion (right subfigure) of box B . Dashed lineswith red dots denote equivalent surfaces, where red dots can be viewed as equivalentsources. Dotted lines with blue dots denote check surface, where blue dots can beviewed as check points. Green shadow is the near range of B . Circles with “+” denotesource points. Step (1) in blue arrow is the evaluation of potential on check surface,and step (2) in red arrow is the calculation of equivalent density on equivalent surface. An upward density function φ B,u ( y ), or the density φ B,uk = φ B,u ( y B,uk ) on several alculation of photoionization in streamer discharges by FMM y B,uk ∈ y B,u , is introduced to represent the potentialin F ( B ), generated by the source inside B . If the upward check potential q B,u ( x ) atthe check surface x B,u , evaluated from the source in B , is equal to the potential q B,u ( x )evaluated from the equivalent source φ B,uk , then these source density points φ B,uk canbe used to represent the potential outside the check surface x B,u including F ( B ). Thisis because of the uniqueness of the Dirichlet boundary value problem (similar to themethod of image charges in electrostatics). The equality is written as X k ∈ Id ( y B,u ) G ( x B,uj , y
B,uk ) φ B,uk = q B,u ( x B,uj ) = X i ∈ I Bs G ( x B,uj , y i ) I ( y i ) , ∀ j ∈ Id ( x B,u ) , (18)where I Bs is the index set of the source points inside B , Id ( y B,u ) is the index set ofdiscrete source points on y B,u and Id ( x B,u ) is the index set of discrete check points on x B,u . A prescribed number m is used to denote the number of discrete equivalent sourcepoints at each side of y B,u , and this number is identical to the number of check pointsat each side of x B,u . Local expansion
Local expansion is used to represent the potential inside a box B ,generated by the source in F ( B ). Similar to the multipole expansion, a downwardequivalent surface y B,d with downward equivalent density φ B,d on it, is introduced. Atthe same time, downward check surface x B,d with downward check potential q B,d is usedto check the equality of potential generated by the source in F ( B ) and the one generatedby φ B,d . Different from the multipole expansion, y B,d should enclose x B,d , since in thelocal expansion we want to approximate the potential inside B . Again both y B,d and x B,d should locate between B and F ( B ). Two surfaces are shown in Figure 3 as anexample, with evaluation procedure.The downward equivalent density satisfies: X k ∈ Id ( y B,d ) G ( x B,dj , y
B,dk ) φ B,dk = q B,d ( x B,dj ) = X i ∈ I F ( B ) s G ( x B,dj , y i ) I ( y i ) , ∀ j ∈ Id ( x B,d ) , (19)where I F ( B ) s is the index set of the source points in F ( B ), Id ( y B,d ) is the index set ofsource points on y B,d and Id ( x B,d ) is the index set of check points on x B,d . Again aprescribed finite number (related to m ) of index is chosen in Id ( · ). M2M translation
M2M translation translates the upward equivalent density φ B,u ofa box, to the upward equivalent density φ P ( B ) ,u of its parent box P ( B ). The idea issimilar to (18), with an upward check surface of P ( B ) as x P ( B ) ,u , and the correspondingupward check potential q P ( B ) ,u . The equality is given as q P ( B ) ,u ( x P ( B ) ,uj ) = X k ∈ Id ( y P ( B ) ,u ) G ( x P ( B ) ,uj , y P ( B ) ,uk ) φ P ( B ) ,uk = X i ∈ Id ( y B,u ) G ( x P ( B ) ,uj , y B,ui ) φ B,ui , ∀ j ∈ Id ( x P ( B ) ,u ) . (20) alculation of photoionization in streamer discharges by FMM q P ( B ) ,u , we evaluate the upward equivalent density φ P ( B ) ,u which is marked asred arrow in the same subfigure. This implementation, which is adding potential to thecheck surface and then calculating the equivalent density from check potential, is alsoapplied to the calculation of downward equivalent density. Therefore, we also indicatethe implementation by blue and red arrows in other subfigures related to M2L and L2Ltranslations in Figure 4. B P ( B ) x P ( B ) ,u y P ( B ) ,u y B,u (1)(2)
A B (1) (2) x B,d y A,u y B,d B P ( B ) y P ( B ) ,d y B,d x B,d (1) (2)
Figure 4.
Cross section of M2M translation (left subfigure), M2L translation (middlesubfigure) and L2L translation (right subfigure). Dashed lines with red dots denoteequivalent surfaces. Dotted lines with blue dots denote check surface. Step (1) in bluearrow is the evaluation of potential on check surface, and step (2) in red arrow is thecalculation of equivalent density on equivalent surface. P ( B ) denotes parent box of B . M2L translation
Two boxes A and B are well-separated if A ⊂ F ( B ) and B ⊂ F ( A ).If two boxes A and B are in same size and well-separated, M2L translation can beused to translate the multipole expansion of A to local expansion of B . In other words,M2L translation calculates the downward equivalent density of B from the upwardequivalent density of A , which accumulates the potential in B from the source in A . Seethis procedure in Figure 4, which satisfies X k ∈ Id ( y B,d ) G ( x B,dj , y
B,dk ) φ B,dk = q B,d ( x B,dj ) = X i ∈ Id ( y A,u ) G ( x B,dj , y
A,ui ) φ A,ui , ∀ j ∈ Id ( x B,d ) . (21)Since these two boxes A and B are in same size, fast Fourier transform can be used tospeed up the calculation in (21), as indicated in [44]. L2L translation
For a box B , F ( P ( B )) ⊂ F ( B ), therefore L2L translation is used tocalculation the local expansion of B from the local expansion of P ( B ). This specifies alculation of photoionization in streamer discharges by FMM B from the source in F ( P ( B )). Similar as (20), the equation is q B,d ( x B,dj ) = X k ∈ Id ( y P ( B ) ,d ) G ( x B,dj , y P ( B ) ,dk ) φ P ( B ) ,dk = X i ∈ Id ( y B,d ) G ( x B,dj , y
B,di ) φ B,di , ∀ j ∈ Id ( x B,d ) . (22) Outline of the algorithm
The outline steps of the kernel-independent FMM is presentedas follows:1. Tree construction: to construct a hierarchical tree in pre-order traversal, such thateach leaf box contains no more than m source points.2. Upward pass: to calculate the multipole expansion for leaf boxes, and use M2Mtranslation for multipole expansion of all non-leaf boxes in a post-order traversal ofthe hierarchical tree.3. Downward pass: for non-root boxes, use local expansion, M2L and L2L translationsto accumulate the potential from far range in a pre-order traversal of the tree.4. Target potential: for each leaf box in pre-order traversal of the tree, sum up the nearinteractions with the potential calculated in the last step, get the target potential.We would like to remark that we tested the backward stable pseudo-inverse trickindicated in [28] for (1), and find few difference with results given by the original pseudo-inverse in [44] for our problems. Therefore, we implement the kernel-independent FMMby the original pseudo-inverse with prescribed number m = 6 (number of equivalentsource or check points at each side of the enclosing box) in this paper. We find that m = 6 gives a good balance between accuracy and computational cost. Readers mayconsider increasing m to obtain a more accurate result, or decreasing m for a fastercomputation.It should be emphasized that in order to capture the multiscale structure ofstreamers, a non-uniform mesh may be adopted in the simulation like in [13, 4, 29]. Theaforementioned fast multipole framework still works for non-uniform and unstructuredmeshes.
4. Results and comparison for computing photoionization
In this section, we compare the performance, in terms of accuracy and efficiency, ofdifferent methods for the evaluation of the photoionization S ph defined in (1) with(2). We take V = V ′ = [0 , x d ] × [0 , y d ] × [0 , z d ] cm and denote its center as ~x = ( x , y , z ) T = ( x d / , y d / , z d / T cm. The box V is partitioned uniformly by n x × n y × n z cells, with n x , n y and n z the number of cells along x , y , z directions,respectively.For simplicity, different numerical methods to be compared are summarized inTable 3. The numerical simulations were performed on the Tianhe2-JK cluster alculation of photoionization in streamer discharges by FMM Table 3.
Notations of several methods introduced in this paper.Notation of method Brief descriptionClassical Int Direct calculation on (1), with (2) and (3).Helmholtz zero BC Three-term summation on (5),by solving (6) with zero boundary condition.Helmholtz Int BC Three-term summation on (5),by solving (6) with integral boundary condition from (1). SP Larsen BC Three-group summation on (9),by solving (12)–(13) with boundary conditions (14)–(15). SP Int BC Three-group summation on (9),by solving (12)–(13) with integral boundary condition (1).FMM classical Int Fast multipole method based on (1), with (2) and (3). located at Beijing Computational Science Research Center. More details can befound at . In ourcomputations via the MPI parallelism, excepted stated otherwise, we always use 4 nodeswith 20 cores in each node in the simulation.The accuracy of different numerical methods is quantified by the following relativeerrors: E V := k S numph ( ~x ) − S refph ( ~x ) k k S refph ( ~x ) k × , E δ ( ~x ) := 1 N tot X | ~x − ~x |≤ δ | S numph ( ~x ) − S refph ( ~x ) | S refph ( ~x ) × , (23)where k · k is the standard discrete l -norm on V , ~x ∈ V , δ > S refph ( ~x ) is the reference result calculated by the (discrete) Classical Intmethod, S numph ( ~x ) is the numerical approximation by a numerical method, and N tot is thenumber of grid points located within a δ radius of ~x . In fact, here E V and E δ ( ~x ) canbe regarded as the global relative error over the whole domain V and the local relativeerror over a ball centered at ~x with radius δ , respectively. The first example is to compute the photoionization rate S ph ( ~x ) in (1) generated froma single Gaussian emission source, which is taken from [6]. The Gaussian ionizationsource S i ( ~x ) in (2) is given as S i ( ~x ) = 1 . × exp (cid:0) − (cid:0) ( x − x ) + ( y − y ) + ( z − z ) (cid:1) /σ (cid:1) cm − s − , (24)where σ > p q = 30 Torr, p = 760 Torr, ξ = 0 . ω/α = 0 . p O = 150 Torr. Wetake δ = 5 σ in (23). alculation of photoionization in streamer discharges by FMM n x = n y = 320 and n z = 160 because directcomputation of the classical integral (1) is too time-consuming even if parallel computingis utilized.Similar to [6], we demonstrate the influence of different ranges of p O r by consideringtwo different sizes of the domain V :(i) x d = y d = 0 . z d = 0 . σ = 0 .
01 cm;(ii) x d = y d = 0 .
04 cm, z d = 0 .
02 cm, σ = 0 .
001 cm.The numerical results are shown in Figures 5 and 6, and the relative errors are thenshown in Tables 4 and 5. z (cm) S ph ( c m - s - ) Classical IntFMM classical IntHelmholtz Int BCHelmholtz zero BC x (cm) y ( c m ) z (cm) S ph ( c m - s - ) Classical IntFMM classical IntSP Int BCSP Larsen BC x (cm) y ( c m ) Figure 5.
Photoionization rate S ph calculated from one Gaussian source in (24). x d = y d = 0 . z d = 0 . σ = 0 .
01 cm. The figures in the left column are S ph along line x = y = 0 . S ph on the plane z = 0 . × , 2 × ,2 × , 2 × cm − s − . As it is clearly shown in Figure 5 and Figure 6, the FMM classical Int methodalways gives the most accurate results, especially for the smaller domain. In all thefigures, the lines for “FMM classical Int” almost coincide with the lines for “ClassicalInt”. The deviations of the solutions of the other four methods from the reference resultsare clearly observable, especially in the central area where the peak locates. Near theboundaries, the methods based on modified Helmholtz equations and SP equations areaccurate only when the boundary values are computed from direct integration. Tables alculation of photoionization in streamer discharges by FMM Table 4.
Time usage and relative error of methods indicated in Table 3, for the caseof single Gaussian source x d = y d = 0 . z d = 0 . σ = 0 .
01 cm.Method Time usage (s) E V E δ ( ~x )Classical Int 184248 — —FMM classical Int 27.1897 0.21% 1.30%Helmholtz zero BC 3.66044 25.33% 16.37%Helmholtz Int BC 3.76133+4606.20 a SP Larsen BC 12.9268 12.05% 8.49% SP Int BC 7.30268+4606.20 a a time usage to compute the boundary values, which is estimated from Classical Intmethod, with multiplication to a factor 2( n x × n y + n x × n z + n y × n z ) / ( n x × n y × n z ). Table 5.
Time usage and relative error of methods indicated in Table 3, for the caseof one Gaussian source x d = y d = 0 .
04 cm, z d = 0 .
02 cm, σ = 0 .
001 cm.Method Time usage (s) E V E δ ( ~x )Classical Int 183761 — —FMM classical Int 27.6406 0.22% 0.53%Helmholtz zero BC 3.87327 83.26% 48.03%Helmholtz Int BC 4.09442+4594.03 a SP Larsen BC 17.3571 68.92% 16.67% SP Int BC 7.88867+4594.03 a a Estimated from the time usage of Classical Int method, with multiplication to afactor 2( n x × n y + n x × n z + n y × n z ) / ( n x × n y × n z ). SP Int BC method take the boundary values from the Classical Int method, and thetime to compute the boundary conditions is also included in Table 4 and Table 5 for a faircomparison. Both tables show that the FMM classical Int method is significantly fasterthan the Classical Int method. In fact, the integration only for the boundary nodesis already much more expensive than the FMM classical Int method. For the threeefficient methods, including FMM classical Int, Helmholtz zero BC, and SP LarsenBC, their computational times have similar magnitudes, and the speed-accuracy trade-off can be observed, meaning that higher computational cost yields better accuracy.Nevertheless, the remarkably lower numerical error and the mildly higher computationalcost of the FMM classical Int method indicate its outstanding competitiveness among allthe approaches for computing the photoionization rates. Additionally, the time usagesof FMM classical Int method are stable for different problem settings, while that of SP Larsen BC method varies significantly (see Tables 4 and 5). alculation of photoionization in streamer discharges by FMM z (cm) S ph ( c m - s - ) Classical IntFMM classical IntHelmholtz Int BCHelmholtz zero BC x (cm) y ( c m ) z (cm) S ph ( c m - s - ) Classical IntFMM classical IntSP Int BCSP Larsen BC x (cm) y ( c m ) Figure 6.
Photoionization rate S ph calculated from one Gaussian source in (24). x d = y d = 0 .
04 cm, z d = 0 .
02 cm, σ = 0 .
001 cm. Left-hand size subfigures are S ph along line x = y = 2 cm, while right-hand size subfigures are contour line of S ph on z = 0 .
01 cm plane, with contour values 2 × , 2 × , 2 × cm − s − . The linecolor and format in the right-hand size subfigures are same as the one in the left-handsize in same row. The second example is to compute the photoionization rate S ph ( ~x ) in (1) generated froma single Gaussian radiation source, which is taken from [25, 26], in order to test the effectof the partial pressure of oxygen p O in the kernel g given in (3). The Gaussian sourceof radiation I in (2) is taken as [7] I ( ~x ) = 4 πξc exp (cid:20) − ( x − x ) + ( y − y ) + ( z − z ) σ (cid:21) cm − s − , (25)where c is the speed of light. Similar to [7], we take ξ = 0 . σ = 0 .
01 cm, c = 2 . × cm · s − , δ = 5 σ and V = [0 , . × [0 , . × [0 , .
4] cm . Wefix the ratio of the partial pressure of oxygen and the air pressure p O /p = 0 . V is uniformly partitioned by 256 × ×
320 cells.In this example, when the partial pressure of oxygen p O in (3) is lower, thephotoionization rate decays slower as the distance from the emission source increases.Therefore, we would like to test the performance of different methods under differentpressures p O . The robustness of the methods under different pressures is of great alculation of photoionization in streamer discharges by FMM Table 6.
Time usage and relative error of methods indicated in Table 3, for the caseof one Gaussian in (25) with p O = 160 Torr.Method Time usage (s) E V E δ ( ~x )Classical Int 292196 — —FMM classical Int 24.7371 0.15% 1.24%Helmholtz zero BC 5.59405 31.93% 16.49%Helmholtz Int BC 5.60994+6391.79 a SP Larsen BC 17.6286 21.31% 8.74% SP Int BC 11.2337+6391.79 a a Estimated from the time usage of Classical Int method, with multiplication to afactor 2( n x × n y + n x × n z + n y × n z ) / ( n x × n y × n z ). significance in practical applications such as sprite discharges.For comparison purpose, two different pressures are considered: (i) p O = 160 Torr;(ii) p O = 10 Torr. The photoionization rate along the central vertical line is plotted inFigures 7 and 8, and the time usage and the numerical error are shown in Tables 6 and7. z (cm) S ph ( c m - s - ) Classical IntFMM classical IntHelmholtz Int BCHelmholtz zero BC z (cm) S ph ( c m - s - ) Classical IntFMM classical IntSP Int BCSP Larsen BC
Figure 7.
Photoionization rate S ph along line x = y = 0 .
125 cm, calculated from oneGaussian in (25) with p O = 160 Torr. Figure 7 shows that in the high-pressure case, all methods give similar resultsdespite obvious mismatch of the peak values. As the pressure decreases, the discrepancybetween different methods becomes more obvious, as shown in Figure 8. The curvesgiven by the Helmholtz and SP methods (with both boundary conditions) deviatesignificantly from the curves of Classical Int method in the low-pressure case. On thecontrary, the results of the FMM classical Int method and the reference Classical Intmethod are in good agreement regardless of the air pressure. The values of the errorsprovided in Tables 6 and 7 again show the advantage of the FMM classical Int methodin accuracy. In fact, for the case of low air pressure, the time used by the FMM classicalInt is quite close to the method of SP Larsen BC.In order to see the relationship between the error and computational time with alculation of photoionization in streamer discharges by FMM z (cm) S ph ( c m - s - ) Classical IntFMM classical IntHelmholtz Int BCHelmholtz zero BC z (cm) S ph ( c m - s - ) Classical IntFMM classical IntSP Int BCSP Larsen BC
Figure 8.
Photoionization rate S ph along line x = y = 0 .
125 cm, calculated from oneGaussian in (25) with p O = 10 Torr. Table 7.
Time usage and relative error of methods indicated in Table 3, for the caseof one Gaussian in (25) with p O = 10 Torr.Method Time usage (s) E V E δ ( ~x )Classical Int 292426 — —FMM classical Int 24.8619 0.19% 0.44%Helmholtz zero BC 6.17337 88.73% 65.87%Helmholtz Int BC 6.01811+6396.82 a SP Larsen BC 23.5598 77.74% 35.11% SP Int BC 12.2226+6396.82 a a Estimated from the time usage of Classical Int method, with multiplication to afactor 2( n x × n y + n x × n z + n y × n z ) / ( n x × n y × n z ). respect to different pressures, we compute more numerical examples under the samesettings with different partial pressures of oxygen ranging from 10 Torr to 160 Torr.Results for the three most efficient method, i.e., FMM classical Int, Helmholtz zero BCand SP Larsen BC methods, are plotted in Figure 9. The FMM classical Int methodalways provides the most accurate results for all pressures, and the error is basicallystable as the pressure varies. For the other two methods, the error increases as pressuredecreases. Moreover, the global relative error E V and the local relative error E δ ( ~x ) nearthe center ~x of the box V of the FMM classical Int method are, in general, at leasttwo order of magnitudes less than those of the other two methods. The time usage ofthe FMM classical Int method is also independent of pressure while the computationtimes of the other two methods increase slightly as the pressure becomes lower. As thepressure decreases, the time cost between the FMM classical Int method and the SP Larsen BC method trends to be the same. alculation of photoionization in streamer discharges by FMM partial pressure of oxygen (Torr) -3 -2 -1 r e l a t i v e - no r m e rr o r t i m e u s age ( s ) FMM classical IntHelmholtz zero BCSP Larsen BC
Figure 9.
Error (red color, solid line) and time usages (blue color, dotted line) ofFMM classical Int, Helmholtz zero BC and SP Larsen BC methods with different airpressures.
5. Results and comparison for computing streamer discharges
To further compare their performances of different methods for treating thephotoionization S ph ( ~x ) in (1), we study the dynamics of streamers with photoionization,where S ph appears as the source term of the transport of charged particles. Thegoverning equations for streamer discharges are given as [4, 23]: ∂n e ∂t − ∇ · ( µ e ~En e ) − ∇ · ( D e ∇ n e ) = S i + S ph ,∂n p ∂t + ∇ · ( µ p ~En p ) = S i + S ph , − ∆ φ = eε ( n p − n e ) , ~E = −∇ φ, (26)where e and ε are the elementary charge and the vacuum dielectric permittivity,respectively; n e := n e ( ~x, t ) and n p := n p ( ~x, t ) are the densities of electrons and positiveions, respectively; µ e and µ p are the mobility coefficients for electrons and positive ions,respectively; D e = diag( D e,x , D e,y , D e,z ), where D e,x , D e,y and D e,z are the diffusioncoefficients in x , y , z directions, respectively; φ and ~E denote the electric potential andelectric field, respectively. Here the photoionization rate S ph is given in (1)-(3) with S i defined as S i = µ e n e α | ~E | , with α := α ( | ~E | ) = 5 . p exp( − p/ | ~E | ) cm − , (27)while α is taken from [11], p is the air pressure, n e and ~E are the solution of (26).The streamer discharge between two parallel plates are used for comparision. Thecomputational domain V is set to be a three-dimensional axis-aligned hyper-rectangle.For the Poisson equation, the Dirichlet boundary conditions are applied on the two facesperpendicular to the z axis, i.e., φ = φ on the upper face and φ = 0 on the bottomface; and the homogeneous Neumann boundary conditions are applied on the other four alculation of photoionization in streamer discharges by FMM n e and the inflow boundaries for n p . The initial value is set as n e ( ~x, t = 0) = n p ( ~x, t = 0) = ˜ n ( ~x ) . (28)Coefficients in (26) are taken from [11, 23]: µ e = 2 . × /p cm / (V · s), µ p = 2 . × /p cm / (V · s), p = 760 Torr, φ = 52 kV and D e = diag(2190 , , / s. Twophysics constant are taken as: e = 1 . × − C and ε = 8 . × − F · cm − . The other physics parameters in (1)-(3) are chosen as [38, 6]: p q = 30 Torr, ξ = 0 . ω/α = 0 . − . The other elliptic equations (6), (12) and (13) are solved by the samealgebraic elliptic solver, with weaker stopping condition that the relative residual is lessthan 10 − . In this subsection, we consider the interaction of two double-headed streamers andcompare the numerical results of the three most efficient methods: FMM classical Intmethod, Helmholtz zero BC method and SP Larsen BC method.The initial value ˜ n in (28) is taken as˜ n ( ~x ) = 10 (cid:16) exp (cid:0) − (cid:0) ( x − . + ( y − . + ( z − . (cid:1) / (0 . (cid:1) + exp (cid:0) − (cid:0) ( x − . + ( y − . + ( z − . (cid:1) / (0 . (cid:1) (cid:17) cm − . The computational domain is fixed as V = [0 , . × [0 , . × [0 ,
1] cm , which ispartitioned by a uniform grid of 512 × × t = 2 . × − ns, where 1 ns = 1 × − s. In order to see the interaction with respectto different p O , we pick two values as p O = 1 Torr and p O = 150 Torr, respectively, inour simulations.We first compare the three methods by observing the electron density. The contoursof the electron densities at 1.5 ns are shown in Figure 10, where the curves of differentmethods are plotted as different line styles and colours. Generally, the results of thethree different methods are in good agreement in most part of the domain for bothpartial pressures of oxygen, while some differences can be observed at the heads ofstreamers. The differences are particularly obvious at the head of positive streamer,which is zoomed in the same figure. The generally good agreement can be attributedto a stronger influence of the impact ionization comparing to the photoionization in the alculation of photoionization in streamer discharges by FMM x (cm) z ( c m ) FMM classical IntHelmholtz zero BCSP Larsen BC x (cm) z ( c m ) FMM classical IntHelmholtz zero BCSP Larsen BC
Figure 10.
Contours of different electron density values at n e = 1 × , 5 × ,9 × , 1 . × cm − , on plane y = 0 .
25 cm at 1 . p O = 150 Torr (left) and p O = 1 Torr (right). Additionally, Figure 10 also displays the difference among three methods withrespect to different p O . As expected from Section 4.2, the difference between threemethods are smaller in higher partial pressure of oxygen (150 Torr), and more observablewhen p O is lower (1 Torr). This implies the validity of using the Helmholtz zero BCmethod and SP Larsen BC for the photoionization in higher p O and also the necessityof using the FMM classical Int method in lower p O .Besides observing the electron density at a fixed time 1.5 ns in Figure 10, the thirdcomponent E z of the electric field ~E = ( E x , E, y, E z ) T along the line x = y = 0 .
25 cmat 0.5 ns, 1.0 ns and 1.5 ns is also shown in Figure 11. As expected from Figure 10, thedifferences of the three methods are generally small, while the difference are easier to beobserved near the heads of streamers (the leftmost and rightmost minimum points). Thedifference is larger when p O is small as 1 Torr, and the deviation increases over time,which is consistent with the results in [7]. One can see that, in the result of the FMMclassical int method, the head of the streamer propagates slightly faster than the othertwo, which is possibly due to the underestimation of photoionization using the Helmholtzzero BC method and SP Larsen BC method. As a summary, these results indicate thatthe accurate approximation of the photoionization could be significant in simulationswith long-time propagation of streamers, especially when the partial pressure of oxygenis low.
As demonstrated previously, one advantage of the FMM method is the scalability inparallel computing with distributed memory, which means the ability to reduce the alculation of photoionization in streamer discharges by FMM z (cm) -150-100-500 E z ( k V / c m ) FMM classical IntHelmholtz zero BCSP Larsen BC z (cm) -150-100-500 E z ( k V / c m ) FMM classical IntHelmholtz zero BCSP Larsen BC
Figure 11.
The third component E z of the electric field ~E along line x = y = 0 .
25 cm,at 0 .
5, 1 . . p O = 150 Torr (left) and p O = 1 Torr (right). Table 8.
Time usage (s) using different nodes over two meshes. 20 cores are used ineach node.Number of nodes 1 2 4 8 16 32 64Mesh size: 256 × ×
160 3843.02 2030.53 1057.19 569.106 304.626 157.651 90.8405Mesh size: 512 × ×
320 31198.5 16329.4 7998.79 4093.52 2116.02 1140.97 621.827 execution time as the number of processes increases. In this subsection, we study thescalability of the FMM classical Int method, which is quantified by the relative speed-up, defined by the ratio of the execution time using the smallest number of cores overthe execution time of the parallel program.In this test, the governing equation is again (26), and the initial value ˜ n in (28) isset as one Gaussian,˜ n ( ~x ) = 10 exp (cid:0) − (cid:0) ( x − . + ( y − . + ( z − . (cid:1) / (0 . (cid:1) cm − . All physics parameters are similar as those in previous subsection except statedotherwise. We set the computational domain as [0 , . × [0 , . × [0 , .
2] cm , andadopt two uniform meshes with 256 × ×
160 and 512 × ×
320 grid cells. Thesimulation is run until 5 × − ns with a fixed time step 1 × − ns. It should benoted that S ph is evaluated twice in each time step, and therefore the FMM classicalInt method is applied 100 times in one simulation.The time usage for the FMM classical Int method in whole simulation (100evaluations) using different numbers of CPU cores is given in Table 8 and plotted inFigure 12, where a satisfactory scalability can be observed.
6. Conclusion
This paper focuses on the accurate and efficient calculation of the photoionization,and proposes the kernel-independent fast multipole method to directly compute the alculation of photoionization in streamer discharges by FMM number of nodes r e l a t i v e s peedup IdealFMM classical Int 10 number of nodes r e l a t i v e s peedup IdealFMM classical Int
Figure 12.
Relative speedups over two meshes: 256 × ×
160 (left) and512 × ×
320 (right). 20 cores are used in each node. photoionization rate efficiently.Quantified accuracy and time usage of the fast multipole method are studied incomparison of the classical integral model and existing approximation models basedon conversion to differential equations. The comparison shows when calculating thephotoionization, the fast multipole method outperforms previous approximations in thefollowing senses: (i) it is significantly efficient (or computationally cheaper) comparedwith the direct calculation by the classical integral under similar accuracy; (ii) it is muchmore accurate compared with those PDE-based approximations (with simple or efficientboundary conditions) under similar computational cost, and it is highly accurate whenthe pressure is not too high; (iii) there are no fitting parameters and it is easy to beextended to unstructured meshes; and (iv) it is more robust with respect to domainsizes and pressures.In summary, in terms of efficiency and accuracy as well as applicability toarbitrary domain with unstructured mesh, the fast multipole method demonstratesbetter performance than those existing numerical methods for the calculation ofphotoionization in streamer discharges in the literature.Finally, we remark here that it is straightforward to apply the kernel-independentfast multipole method for computing photoionization models with different integralforms, and thus we provide a practical tool for the validation of newly proposed integralmodels of photoionization in streamer discharges. Future works include applying themethod to other integral models and taking the stochastic effect into consideration.
Acknowledgments
This work was partially supported by the National Science Foundation of China underproject 51921005 (C. Zhuang and R. Zeng), the Academic Research Fund of the Ministryof Education of Singapore under grants R-146-000-305-114 (Z. Cai) and R-146-000-290-114 (B. Lin and W. Bao). Some computations were done on the Tianhe2-JK cluster atthe Beijing Computational Science Research Center under the kind support by Professor alculation of photoionization in streamer discharges by FMM
References [1] M Aints, A Haljaste, and L Roots. Photoionizing radiation of positive corona in moist air. In
Proc. 8th Int. Symp. on High Pressure Low Temperature Plasma Chemistry (Estonia, 21–25July 2002) , 2002.[2] Natalia Yu Babaeva, Ananth N Bhoj, and Mark J Kushner. Streamer dynamics in gases containingdust particles.
Plasma Sources Science and Technology , 15(4):591, 2006.[3] B Bagheri and J Teunissen. The effect of the stochasticity of photoionization on 3d streamersimulations.
Plasma Sources Science and Technology , 28(4):045013, 2019.[4] Delphine Bessi`eres, Jean Paillol, Anne Bourdon, Pierre S´egur, and Emmanuel Marode. A newone-dimensional moving mesh method applied to the simulation of streamer discharges.
Journalof Physics D: Applied Physics , 40(21):6559, 2007.[5] Annemie Bogaerts, Erik Neyts, Renaat Gijbels, and Joost Van der Mullen. Gas discharge plasmasand their applications.
Spectrochimica Acta Part B: Atomic Spectroscopy , 57(4):609–658, 2002.[6] A Bourdon, V P Pasko, N Y Liu, S C´elestin, P S´egur, and E Marode. Efficient models forphotoionization produced by non-thermal gas discharges in air based on radiative transfer andthe helmholtz equations.
Plasma Sources Science and Technology , 16(3):656–678, aug 2007.[7] Julien Capeill`ere, Pierre S´egur, Anne Bourdon, S´ebastien C´elestin, and Sergey Pancheshnyi. Thefinite volume method solution of the radiative transfer equation for photon transport in non-thermal gas discharges: application to the calculation of photoionization in streamer discharges.
Journal of Physics D: Applied Physics , 41(23):234018, nov 2008.[8] S´ebastien Celestin.
Study of the dynamics of streamers in air at atmospheric pressure . Theses,Ecole Centrale Paris, December 2008.[9] O. Chanrion and T. Neubert. A pic-mcc code for simulation of streamer propagation in air.
Journal of Computational Physics , 227(15):7222 – 7245, 2008.[10] Hongwei Cheng, Leslie Greengard, and Vladimir Rokhlin. A fast adaptive multipole algorithm inthree dimensions.
Journal of computational physics , 155(2):468–498, 1999.[11] SK Dhali and PF Williams. Two-dimensional studies of streamers in gases.
Journal of AppliedPhysics , 62(12):4696–4707, 1987.[12] A Dubinova, D Trienekens, U Ebert, S Nijdam, and T Christen. Pulsed positive discharges in air atmoderate pressures near a dielectric rod.
Plasma Sources Science and Technology , 25(5):055021,sep 2016.[13] Ute Ebert, Carolynne Montijn, Tanja MP Briels, Willem Hundsdorfer, Bernard Meulenbroek,Andrea Rocco, and Eddie M van Veldhuizen. The multiscale nature of streamers.
PlasmaSources Science and Technology , 15(2):S118, 2006.[14] Ute Ebert and Davis D Sentman. Streamers, sprites, leaders, lightning: from micro-to macroscales.
Journal of Physics D: Applied Physics , 41(23):230301, 2008.[15] Leslie Greengard and Vladimir Rokhlin. A fast algorithm for particle simulations.
Journal ofcomputational physics , 73(2):325–348, 1987.[16] Leslie F Greengard and Jingfang Huang. A new version of the fast multipole method for screenedcoulomb interactions in three dimensions.
Journal of Computational Physics , 180(2):642–658,2002.[17] Nail A Gumerov and Ramani Duraiswami. Fast multipole method for the biharmonic equation inthree dimensions.
Journal of Computational Physics , 215(1):363–383, 2006.[18] Ming Jiang, Yongdong Li, Hongguang Wang, Pengfeng Zhong, and Chunliang Liu. Aphotoionization model considering lifetime of high excited states of n2 for pic-mcc simulationsof positive streamers in air.
Physics of Plasmas , 25(1):012127, 2018.[19] Ravindra P Joshi and Selma Mededovic Thagard. Streamer-like electrical discharges in water: alculation of photoionization in streamer discharges by FMM Part ii. environmental applications.
Plasma Chemistry and Plasma Processing , 33(1):17–49,2013.[20] A A Kulikovsky. The role of photoionization in positive streamer dynamics.
Journal of PhysicsD: Applied Physics , 33(12):1514–1524, may 2000.[21] Edward W Larsen, Guido Th¨ommes, Axel Klar, Mohammed Seaıd, and Thomas G¨otz. SimplifiedPN approximations to the equations of radiative heat transfer and applications.
Journal ofComputational Physics , 183(2):652 – 675, 2002.[22] Zhi Liang, Zydrunas Gimbutas, Leslie Greengard, Jingfang Huang, and Shidong Jiang. A fastmultipole method for the rotne–prager–yamakawa tensor and its applications.
Journal ofComputational Physics , 234:133–139, 2013.[23] Bo Lin, Chijie Zhuang, Zhenning Cai, Rong Zeng, and Weizhu Bao. An efficient and accurate MPI-based parallel simulator for streamer discharges in three dimensions.
Journal of ComputationalPhysics , 401:109026, 2020.[24] Ningyu Liu, S´ebastien C´elestin, Anne Bourdon, Victor P Pasko, Pierre S´egur, and EmmanuelMarode. Application of photoionization models based on radiative transfer and the helmholtzequations to studies of streamers in weak electric fields.
Applied Physics Letters , 91(21):211501,2007.[25] Ningyu Liu, Joseph R Dwyer, Hans C Stenbaek-Nielsen, and Matthew G McHarg. Sprite streamerinitiation from natural mesospheric structures.
Nature communications , 6(1):1–9, 2015.[26] Ningyu Liu and Victor P Pasko. Effects of photoionization on propagation and branching ofpositive and negative streamers in sprites.
Journal of Geophysical Research: Space Physics ,109(A4), 2004.[27] Alejandro Luque, Ute Ebert, Carolynne Montijn, and Willem Hundsdorfer. Photoionization innegative streamers: Fast computations and two propagation modes.
Applied Physics Letters ,90(8):081501, 2007.[28] Dhairya Malhotra and George Biros. PVFMM: A parallel kernel independent FMM for particleand volume potentials.
Communications in Computational Physics , 18(3):808830, 2015.[29] Robert Marskar. 3d fluid modeling of positive streamer discharges in air with stochasticphotoionization.
Plasma Sources Science and Technology , 2020.[30] Dana Moudry, Hans Stenbaek-Nielsen, Davis Sentman, and Eugene Wescott. Imaging of elves,halos and sprite initiation at 1ms time resolution.
Journal of Atmospheric and Solar-TerrestrialPhysics , 65(5):509–518, 2003.[31] G V Naidis. On photoionization produced by discharges in air.
Plasma Sources Science andTechnology , 15(2):253–255, mar 2006.[32] Sergey Pancheshnyi. Role of electronegative gas admixtures in streamer start, propagation andbranching phenomena.
Plasma Sources Science and Technology , 14(4):645, 2005.[33] Sergey Pancheshnyi. Photoionization produced by low-current discharges in o2, air, n2and CO2.
Plasma Sources Science and Technology , 24(1):015023, dec 2014.[34] SV Pancheshnyi, SM Starikovskaia, and A Yu Starikovskii. Role of photoionization processes inpropagation of cathode-directed streamer.
Journal of Physics D: Applied Physics , 34(1):105,2001.[35] G. W. Penney and G. T. Hummert. Photoionization measurements in air, oxygen, and nitrogen.
Journal of Applied Physics , 41(2):572–577, 1970.[36] Samaneh Sadighi, Ningyu Liu, Joseph R Dwyer, and Hamid K Rassoul. Streamer formation andbranching from model hydrometeors in subbreakdown conditions inside thunderclouds.
Journalof Geophysical Research: Atmospheres , 120(9):3660–3678, 2015.[37] WJM Samaranayake, Y Miyahara, T Namihira, S Katsuki, T Sakugawa, R Hackam, andH Akiyama. Pulsed streamer discharge characteristics of ozone production in dry air.
IEEETransactions on Dielectrics and Electrical Insulation , 7(2):254–260, 2000.[38] P S´egur, A Bourdon, E Marode, D Bessieres, and J H Paillol. The use of an improved Eddingtonapproximation to facilitate the calculation of photoionization in streamer discharges.
Plasma alculation of photoionization in streamer discharges by FMM Sources Science and Technology , 15(4):648–660, jul 2006.[39] J Stephens, M Abide, A Fierro, and A Neuber. Practical considerations for modelingstreamer discharges in air with radiation transport.
Plasma Sources Science and Technology ,27(7):075007, 2018.[40] Pavel ˇSunka. Pulse electrical discharges in water and their applications.
Physics of plasmas ,8(5):2587–2594, 2001.[41] Jannis Teunissen and Ute Ebert. 3d PIC-MCC simulations of discharge inception around a sharpanode in nitrogen/oxygen mixtures.
Plasma Sources Science and Technology , 25(4):044005, jun2016.[42] Zhongmin Xiong and Mark J Kushner. Branching and path-deviation of positive streamersresulting from statistical photon transport.
Plasma Sources Science and Technology ,23(6):065041, 2014.[43] Won J Yi and P F Williams. Experimental study of streamers in pure n2and n2/o2mixtures anda ≈
13 cm gap.
Journal of Physics D: Applied Physics , 35(3):205–218, jan 2002.[44] Lexing Ying, George Biros, and Denis Zorin. A kernel-independent adaptive fast multipolealgorithm in two and three dimensions.
Journal of Computational Physics , 196(2):591 – 626,2004.[45] Lexing Ying, George Biros, Denis Zorin, and Harper Langston. A new parallel kernel-independentfast multipole method. In
SC’03: Proceedings of the 2003 ACM/IEEE conference onSupercomputing , pages 14–14. IEEE, 2003.[46] M. B. Zheleznyak, A. Kh. Mnatsakanian, and Sergei Vasil’evich Sizykh. Photoionization ofnitrogen and oxygen mixtures by radiation from a gas discharge.
High Temperature Science ,20:423–428, 1982.[47] C. Zhuang, B. Lin, R. Zeng, L. Liu, and M. Li. 3-D parallel simulations of streamer discharges inair considering photoionization.
IEEE Transactions on Magnetics , 56(3):7513804, 2020.[48] C. Zhuang and R. Zeng. A local discontinuous Galerkin method for 1.5-dimensional streamerdischarge simulations.
Applied Mathematics and Computation , 219:9925–9934, 2013.[49] C. Zhuang and R. Zeng. A WENO scheme for simulating streamer discharge with photoionizations.