Noise and error analysis and optimization in particle-based kinetic plasma simulations
NNoise and error analysis and optimization in particle-based kineticplasma simulations
E. G. Evstatiev ∗1 , J. M. Finn , B. A. Shadwick , and N. Hengartner Sandia National Laboratory, Albuquerque, NM Tibbar Plasma Technologies, Los Alamos, NM University of Nebraska-Lincoln, Lincoln, NE Los Alamos National Laboratory, Los Alamos, NM
16 August 2020
Abstract
In this paper we analyze the noise in macro-particle methods used in plasma physics and fluid dynamics,leading to approaches for minimizing the total error, focusing on electrostatic models in one dimension. Webegin by describing kernel density estimation for continuous values of the spatial variable x , expressing thekernel in a form in which its shape and width are represented separately. The covariance matrix C ( x, y ) of the noise in the density is computed, first for uniform true density. The band width of the covariancematrix is related to the width of the kernel. A feature that stands out is the presence of constant negativeterms in the elements of the covariance matrix both on and off-diagonal. These negative correlations arerelated to the fact that the total number of particles is fixed at each time step; they also lead to the property ´ C ( x, y ) dy = 0 . We investigate the effect of these negative correlations on the electric field computed byGauss’s law, finding that the noise in the electric field is related to a process called the Ornstein-Uhlenbeckbridge , leading to a covariance matrix of the electric field with variance significantly reduced relative to thatof a Brownian process.For non-constant density, ρ ( x ) , still with continuous x , we analyze the total error in the density estimationand discuss it in terms of bias-variance optimization (BVO). For some characteristic length l , determined bythe density and its second derivative, and kernel width h , having too few particles within h leads to too muchvariance; for h that is large relative to l , there is too much smoothing of the density. The optimum betweenthese two limits is found by BVO. For kernels of the same width, it is shown that this optimum (minimum)is weakly sensitive to the kernel shape.We repeat the analysis for x discretized on a grid. In this case the charge deposition rule is determinedby a particle shape . An important property to be respected in the discrete system is the exact preservationof total charge on the grid; this property is necessary to ensure that the electric field is equal at both ends,consistent with periodic boundary conditions. We find that if the particle shapes satisfy a sum rule, theparticle charge deposited on the grid is conserved exactly. Further, if the particle shape is expressed as theconvolution of a kernel with another kernel that satisfies the sum rule, then the particle shape obeys the sumrule. This property holds for kernels of arbitrary width, including widths that are not integer multiples ofthe grid spacing.We show results relaxing the approximations used to do BVO optimization analytically, by doing numericalcomputations of the total error as a function of the kernel width, on a grid in x . The comparison betweennumerical and analytical results shows good agreement over a range of particle shapes.We discuss the practical implications of our results, including the criteria for design and implementationof computationally efficient particles that take advantage of the developed theory. The particle-in-cell (PIC) method has been an indispensable tool of numerical modelers in fluid dynamics andkinetic plasma physics for several decades now [1–8] and the variety of kinetic plasma problems to which ithas been applied keeps increasing. The success of this method has spurred more recent developments [9–14],with emphasis on geometrical aspects; for example, advantage has been taken of the Hamiltonian nature of the ∗ Corresponding author: [email protected] a r X i v : . [ phy s i c s . p l a s m - ph ] F e b lasov-Maxwell system [15,16]. To ensure conservation properties such as momentum, energy, charge, etc., someformulations [4, 11, 12, 17] rely on a variational method, related to the Hamiltonian prescription, while othersdevise specific spatial and temporal discretizations for providing conservation properties [9, 10, 17].The first particle methods applied to plasma simulations [2] showed significant effects of noise due to a fewfactors: first, the deposition of a particle’s charge applied to only one grid node, an approach called the nearestgrid point (NGP) method; second, very few particles were used, due to the limited computational power in theearly 1970s; third, guidelines for noise minimization were quite limited.The introduction of finite size computational particles by Birdsall and Langdon [3] significantly elevated theusefulness of the particle method. Indeed, a hierarchy of shapes, varying in size and smoothness, were proposedto address issues of noise [18, 19] as well as frequency aliasing [8].The recognition that the level of particle noise scales as N − / p , where N p is the number of computationalparticles, together with ever-increasing demands for accuracy and fidelity of simulations has again put the issueof noise in particle methods in the spotlight. Advancement has come with the emergence of hybrid kinetic-fluidmethods and the δf method [20–24]. In spite of this progress, noise is still a limiting factor in particle codes: in δf and hybrid methods particles are used to describe only a subset of the distribution function, however, noiseis still an important factor for the particle part of the computations as well as for the fluid-particle coupling(for time-evolving fluids). In full kinetic treatments as well in hybrid and δf methods, noise in the density is anespecially serious problem for quasineutral plasmas, in which the local net charge density is small.The present work presents a new analysis of the statistics of noise in particle-based methods. Specifically,we analyze the error in the estimation of the particle density in terms of a finite number of computationalparticles of finite size, a special case of kernel density estimation , with a focus on the bias-variance trade-off [25, 26]. We also analyze the error in the electric field computed from the charge density, showing thatcertain negative correlations in the density noise lead to properties of the electric field related to the Ornstein-Uhlenbeck bridge [27–29], a generalization of the
Brownian bridge [30], a Brownian process with boundaryconditions at each end. We concentrate on a 1D (one-dimensional) electrostatic (ES) formulation with periodicboundary conditions, with overall charge neutrality and immobile ions, leaving generalizations such as to higherdimensions and electromagnetic models for future work.In Sec. 2 we establish the framework used in estimating the electron density and its noise properties in a D electrostatic Vlasov-Poisson system. In this system, the only source of noise is the estimated charge density.(Electromagnetic models also involve noise in the estimated current density.) These issues related to densityestimation are introduced with a continuous, i.e. non-discretized, spatial variable x . We also discuss the variouskernels that can be used, show how the kernel width and its shape (smoothness) enter, and summarize someproperties of these kernels that relate to discretization and that will enter in later sections.In Sec. 3, and in the next section, we continue to restrict our attention to continuous x . We introducethe covariance matrix for the noise in the density, focusing in this section on a system in which the “true”density is uniform. We discuss the origin of certain negative terms in the covariance matrix. (These negativeoff-diagonal terms represent negative correlations.) We show that in computing the electric field by Gauss’s law,these negative correlations and the boundary conditions on the electric field lead to properties associated withthe Ornstein-Uhlenbeck bridge. We characterize noise in the computed electric field in terms of its covariancematrix. We illustrate with a kernel involving a delta function. Issues associated with the Brownian bridge arediscussed in more depth in Appendix A and issues related to relaxing the delta function restriction to give anonzero kernel width are discussed in Appendix B.In Section 4 we generalize to non-uniform density and discuss the application of bias-variance optimization [25]to find the optimal kernel width. This optimization in the presence of non-uniform density minimizes the totalerror in the density estimated with a kernel of a specific shape and width. This error consists of a variance term(noise) caused by the finite number of particles and a bias term, a smoothing of the density that occurs becauseof the finite kernel width. For the remainder of this paper we refer to noise as the error due to having a finitenumber of particles and the more general term error as including the bias. Issues relating to the scaling of thekernel width are discussed in Appendix C.In Section 5 we discuss the density and electric field on a discrete grid, where a particle shape for the chargedeposition enters. We discuss the importance of a sum rule; obeying this sum rule is a sufficient condition forthe net charge on the grid to be exactly zero when the ion charge is subtracted. This requirement assures thatthe electric field at the endpoints are equal. We also show that for a general kernel, the sum rule is obeyed if theparticle shape is a convolution of the kernel with another kernel that already satisfies the sum rule. We discussthe covariance matrix of the noise on the discrete grid for various particle shapes.In Sec. 6 we compute the total error (bias plus variance) numerically for various shapes and compare withthe analytic theory of sections 3 and 4. 2n Sec. 7 we summarize and discuss our results. Particle methods are hybrid Lagrangian-Eulerian in nature: computational macro-particles are allowed to movewith continuous positions and velocities, while charge densities and fields are resolved on a fixed computationalgrid. The connection between particles and grid is via a particle shape, which specifies a charge deposition rule.It has been traditional [8] to use particle sizes that are integer number of computational cells wide, althoughsuch restriction is not necessary; a related issue, which we also discuss in Sec. 5.1, is that the width of a particleshape and its degree of smoothness need not be related. The latter point has been emphasized in Ref. [12].Furthermore, grid size is many times determined subjectively by the modeler according to a desired resolution,accuracy, the particular physics problem under consideration, etc. This resolution may be increased a few timesto determine convergence of the numerical results, while also increasing the number of particles; a typical quantitythat is kept constant is the average number of particles per cell. Of course, every time the grid resolution andparticle number are increased, the demand for computational resources increases and for large problems thisstrategy quickly becomes prohibitive. The analysis in this paper aims to provide a systematic way of minimizingnoise and error in the charge density by selecting optimal size and number of particles and, as a consequence, tominimize the computational resources required to achieve a given accuracy.Our discussion will focus on 1D electrostatic models, in periodic geometry. As will be shown in the following,the charge density and its gradients are essential for the analysis of noise and error. Therefore, we considerworking with quantities that are periodic functions of x on the real interval [0 , , but over all real values of theparticle velocity v . This choice is advantageous for the presentation of the ideas, postponing grid discretizationto later sections.We introduce a representation of the electron distribution function in phase space ( x, v ) ∈ [0 , × ( −∞ , ∞ ) in terms of N p number of finite-size computational particles, f e ( x, v, t ) = N p (cid:88) µ =1 q µ K ( x − ξ µ ) δ ( v − ˙ ξ µ ) , (1)where f e is the estimated phase space distribution, q µ is the computational particle charge, ξ µ is the computationalparticle position, and ˙ ξ µ is its velocity. We use q µ > in Eq. (1) and throughout the presentation and the negativesign of the electron charge is added explicitly in places where it is used, e.g., in Gauss’s law (so strictly speaking q µ is the weight and the macro-particle charge is − q µ ). The general form of the kernel is K ( x, ξ ) , however, dueto the periodic boundary conditions it assumes the translationally invariant form K ( x − ξ ) (see below). Thesubscript “ e ” in Eq. (1) and everywhere throughout the paper stands for “estimated.” Also, for the rest of thepaper we use the term particle in place of computational particle or macro-particle. This particle is usuallycomprised of many physical particles.By integrating over velocity space, we obtain the estimated density of the electrons at any spatial point x interms of the positions of all of the particles: ρ e ( x ) = N p (cid:88) µ =1 q µ K ( x − ξ µ ) . (2)We take the special case, in which all the q µ are equal. With periodic boundary conditions on [0 , , no particlesare gained or lost, so ´ ρ e ( x ) dx is conserved. We normalize to ´ ρ e ( x ) dx = 1 and assume immobile ions withuniform, fixed density ρ ( i ) ( x ) . Overall neutrality is assumed, i.e. ´ ρ ( i ) ( x ) dx = 1 as well.The form in Eq. (2) is the usual form of kernel density estimation, used in statistics and machine learning [25].The kernel K ( x ) is usually assumed to satisfy the following conditions, which do not present practical limitations:• Normalized to unity, ˆ K ( x ) dx = 1 ; (3)• Symmetric, K ( x ) = K ( − x ) , x ∈ [0 , (4)• Translationally invariant, K ( x, ξ ) = K ( x − ξ ) , x, ξ ∈ [0 , (5)• Nonnegative, K ( x ) ≥ , x ∈ [0 , (6)• Has compact support. (7)3ondition (3) ensures the density normalization discussed above while conditions (4)–(7) are chosen out ofconvenience but are not essential for the theory development. The normalization condition on the kernel andthe condition ´ ρ e ( x ) dx = 1 imply (cid:80) N p µ =1 q µ = 1 while the assumed equal and constant particle charges lead to q µ = 1 /N p .At this stage, there is no grid, so the kernel width is not related to a grid spacing. We, in fact, express akernel of width h as K ( x ) = 1 h K f (cid:16) xh (cid:17) , (8)where K f is the fundamental kernel with support [ − / , / . Thus, K f contains all the information on theparticle shape, including its smoothness, while its width is independently set by h . (To be specific, we assumethat K f is defined on the real line, and after scaling to form K ( x ) , it is extended to be periodic with periodunity.)Examples of fundamental kernels are given in Table 1 and illustrated in Fig. 1. The boxcar, linear andquadratic kernels are the convolutional particle shapes of Ref. [7, 8], scaled to the unit interval; the trapezoidalkernel is discussed in the following sections. Another important kernel is the Epanechnikov kernel [31]. Sec. 4.1discusses the BVO process, optimal kernel width, etc., where the shape factors ´ K f ( x ) dx and ´ x K f ( x ) dx ,which are of order unity, play a prominent role.In sections 5 and 6 we will apply these results involving the kernel K ( x ) in the presence of a uniform grid x i , i = 1 , , . . . , N g with spacing ∆ = 1 /N g . There we will construct a particle shape S ( x ) , which satisfiesconditions (3)-(7) for a kernel and also satisfies a sum rule N g (cid:88) i =1 ∆ S ( x i − ξ ) = 1 (9)for an arbitrary particle position ξ , the discrete analog of the normalization condition in Eq. (3). We will showthat a sufficient condition for S ( x ) to satisfy the sum rule is that it be a convolution of a kernel of arbitrarywidth δ with either another particle shape or a finite element of width i ∆ , i = 1 , , . . . ; thereby, S ( x ) can havean arbitrary width h = i ∆ + δ .As a last comment in this section, we note that not every kernel when scaled to a grid satisfies the sum rule;among our examples, the Epanechnikov kernel scaled to the grid spacing (width h = i ∆ + δ ) does not obey thesum rule (9) for any δ and i = 1 , , . . . . In contrast, the boxcar, linear, quadratic, and trapezoidal kernels, whenscaled to h = ∆ , , , , correspondingly [cf. Eq. (8)], do satisfy the sum rule. We shall discuss these issuesfurther in Sections 5 and 6. At the same time, since all particle shapes satisfy conditions (3)–(7), any particleshape may be used as a kernel in the density estimation expression (2). Based on the kernel representation from Sec. 2, we analyze in this section what is typically known as “noise” inthe PIC method. We note that—as will become clear in the next section—noise is only one part of the totalerror that we make when estimating the “true” density (or field) with the help of a finite number of particles, theother part being the bias error. When estimating the error in a uniform density, as we do in this section, onlythe noise part appears, the bias part being zero for this case. The focus is on the covariance matrix betweenthe noise in the density at different spatial points. Based on the density covariance matrix, we consider thecovariance matrix between the noise in the electric field at different spatial points.
In this section we introduce the mathematical method of the analyses of noise, and later that of error, whilealso deriving the uniform density correlations with the important negative contributions resulting from the fixednumber of particles in a numerical simulation. Let us denote the true electron density by ρ ( x ) . Because of ourchoice of normalization, the true electron density (or true density) satisfies all the properties of a probabilitydensity function , and the normalization condition (3) guarantees that ρ e ( x ) does also. As discussed above, at We will call the measure of the support the width of the kernel. We restrict our attention to uniform grids strictly for convenience. ernel DefinitionBoxcar (top-hat) K fB ( x ) = (cid:40) , | x | ≤ otherwise .Linear (tent) K fL ( x ) = (cid:40) − | x | ) , | x | ≤ otherwise .Quadratic K fQ ( x ) = 9 − x , | x | ≤ / (cid:0) − | x | (cid:1) , / ≤ | x | ≤ / otherwise .Trapezoidal K fT ( x ) = | x | ≤ / (cid:0) − | x | (cid:1) , / ≤ | x | ≤ / otherwise .Epanechnikov K fE ( x ) = (cid:40) (cid:0) − x (cid:1) , | x | ≤ otherwise , Table 1: Examples of fundamental kernels, with width unity and normalized to have ´ K ( x ) dx = 1 . x F un d a m e n t a l k e r n e l BoxcarLinearQuadraticTrapezoidalEpanechnikov
Figure 1: Illustration of the fundamental kernel examples from Table 1, normalized to unity. The importantissues of bias-variance optimization and the sum rule in Eq. (9) are discussed in the text.this point in our analysis we do not consider grid discretization and hence we do not relate the kernel width h to the grid spacing ∆ . Using the continuous spatial variable x , the average of any quantity f ( x ) is calculated as (cid:104) f (cid:105) = ˆ f ( ξ ) ρ ( ξ ) dξ (10)or (cid:104) f (cid:105) = ˆ f ( ξ, η ) ρ ( ξ ) ρ ( η ) dξdη (11)5or a function of two random variables, etc. In this way, we can calculate the expected (statistical expectation)value of the estimated density over the true density ρ ( x ) as (cid:104) ρ e ( x ) (cid:105) = (cid:32)(cid:88) µ q µ (cid:33) ˆ K ( x − ξ ) ρ ( ξ ) dξ = ˆ x + h/ x − h/ K f (cid:18) x − ξh (cid:19) ρ ( ξ ) dξh = ˆ / − / K f ( η ) ρ ( x + hη ) dη, (12)where (cid:80) µ q µ = 1 has been used, the symmetry of K has been used, and we have defined η = ( ξ − x ) /h . Expandingfor small h , we find (cid:104) ρ e ( x ) (cid:105) = ρ ( x ) ˆ / − / K f ( η ) dη + h ρ (cid:48)(cid:48) ( x ) ˆ / − / K f ( η ) η dη + · · · = ρ ( x ) + h ρ (cid:48)(cid:48) ( x ) ˆ / − / K f ( η ) η dη + · · · . (13)where the normalization ´ / − / K f ( η ) dη = 1 condition has been used and the symmetry of K f has again beenused to conclude ´ / − / K f ( η ) ηdη = 0 . (Henceforth, we omit integral limits to improve readability.) We defer theissues of a non-uniform density ρ ( x ) to a later section.For the uniform density case considered in this section, ρ ( x ) = 1 (hence ρ (cid:48)(cid:48) ( x ) = 0 ), and we obtain (cid:104) ρ e ( x ) (cid:105) = ρ ( x ) = 1 . Now defining the fluctuations ˜ ρ e ( x ) by ρ e ( x ) = (cid:104) ρ e ( x ) (cid:105) + ˜ ρ e ( x ) = 1 + ˜ ρ e ( x ) (note that (cid:104) ˜ ρ e ( x ) (cid:105) = 0 ), we can writethe variance as V ( x ) = (cid:104) ˜ ρ e ( x ) (cid:105) = (cid:104) ρ e ( x ) (cid:105) − , or V ( x ) = V d ( x ) + V o ( x ) − where V d ( x ) and V o ( x ) denote the diagonal ( µ = ν ) and off-diagonal ( µ (cid:54) = ν ) terms in V ( x ) , as defined below.We find V d ( x ) = (cid:88) µ q µ (cid:104) K ( x − ξ µ ) (cid:105) = 1 N p ˆ K ( x − ξ ) dξ = 1 N p ˆ K ( ξ ) dξ (14)and V o ( x ) = (cid:88) µ (cid:54) = ν q µ q ν (cid:104) K ( x − ξ µ ) K ( x − ξ ν ) (cid:105) = N p ( N p − N p ˆ K ( x − ξ ) K ( x − η ) dξdη = N p ( N p − N p (cid:18) ˆ K ( ξ ) dξ (cid:19) , = 1 − N p . (15)At this point we notice the general scaling V d = (1 /N p h ) ´ K f ( η ) dη ∼ /N p h , which follows by using Eq. (8).The quantity N p h is the expected number of particles over the width of the kernel. These results lead to V ( x ) = V d ( x ) − /N p . (16)To relate the estimated density at arbitrary spatial points, x and y , we must compute the covariance matrix C ( x, y ) = (cid:104) ˜ ρ e ( x )˜ ρ e ( y ) (cid:105) = (cid:104) ρ e ( x ) ρ e ( y ) (cid:105) − . (17)We again use a decomposition into diagonal and off-diagonal terms: C ( x, y ) = C d ( x, y ) + C o ( x, y ) − . (18)6he two contributions are C d ( x, y ) = C d ( x − y ) = 1 N p ˆ K ( x − ξ ) K ( y − ξ ) dξ = 1 N p ˆ K ( x − y ) , (19)where ˆ K is the convolution of K with itself (a legitimate kernel according to Eqs. (3) - (7)), and C o ( x, y ) = N p ( N p − N p ˆ K ( x − ξ ) K ( y − η ) dξdη = N p ( N p − N p (cid:18) ˆ K ( ξ ) dξ (cid:19) = 1 − N p . (20)Putting (18), (19), and (20) together, we find C ( x, y ) = 1 N p (cid:104) ˆ K ( x − y ) − (cid:105) . (21)In the special case K ( x ) = δ ( x ) we have ˆ K ( x − y ) = δ ( x − y ) and from Eq. (21) we obtain C ( x − y ) = 1 N p [ δ ( x − y ) − . (22)Notice the translational invariance form C ( x − y ) and the presence of constant negative contributions − /N p in the expressions for the variance ( x = y ) and the off-diagonal (correlation) terms ( x (cid:54) = y ). In particular, wehave ˆ C ( x, y ) dy = ˆ C ( x − y ) dy = 0 . (23)Indeed, ´ C d ( x − y ) dy = (1 /N p ) ´ ˆ K ( x − y ) dy = 1 /N p since ´ ˆ K ( x − y ) dx = 1 . We emphasize that the propertyin Eq. (23) is general, i.e., for any kernel. An alternative proof is given as follows. Recalling that ´ ρ e ( y ) dy = 1 ,it follows that ´ ˜ ρ e ( y ) dy = ´ ( ρ e ( y ) − dy = 0 . Then from the definition of correlations (17) we obtain ˆ C ( x, y ) dy = (cid:28) ˜ ρ e ( x ) ˆ ˜ ρ e ( y ) (cid:29) dy = 0 . The result Eq. (23) implies that the function defined by u ( x ) = 1 is in the null space of the covariance matrix,i.e., is the eigenfunction with zero eigenvalue.The significance of the negative correlations is further discussed in the next section. In particle codes, noise in the density leads to noise in the electric field, which in turn affects particle orbits. Inthis section we quantify the effect of density noise on the electric field. The quantification of errors in particleorbits due to errors in the electric field, leading in turn to density errors, i.e., “closing the loop,” will be thesubject of future work.In the electrostatic model the electric field is computed from Gauss’s law (we use the dimensionless form), dEdx = ρ ( i ) − ρ e = ρ q , (24)where ρ q is the estimated density, ρ e is the estimated electron density, and again, ρ ( i ) = 1 is the fixed backgroundion density. Because of the assumption that the net charge is exactly zero, ´ ρ q ( x ) dx = 0 , we have E (0) = E (1) ,consistent with the assumed periodic boundary conditions. We also specify that there is no applied potentialacross the system, so that ˆ E ( z ) dz = 0 . (25)7o incorporate condition (25) in the solution of (24), we start with the general expression ˆ E ( x ) = ˆ xx ρ q ( z ) dz , (26)where ˆ E satisfies Eq. (24) and x is an arbitrary initial point of integration. We calculate the integral of ˆ E ( x ) over the periodic domain [0 , : R = ˆ ˆ E ( x ) dx = ˆ dx ˆ xx ρ q ( z ) dz = x ˆ xx ρ q ( z ) dz (cid:12)(cid:12)(cid:12)(cid:12) x =1 x =0 − ˆ xρ q ( x ) dx = ˆ x ρ q ( z ) dz − ˆ zρ q ( z ) dz = ˆ ρ q ( z ) dz − ˆ x ρ q ( z ) dz − ˆ z ρ q ( z ) dz = − ˆ x ρ q ( z ) dz − ˆ x ρ q ( x ) dx, where in the last line we have used ´ ρ q ( z ) dz = 0 . The quantity R depends on x . The quantity R needs tobe subtracted from ˆ E ( x ) in order to obtain an expression E ( x ) satisfying Eq. (25): E ( x ) = ˆ E ( x ) − R = ˆ xx ρ q ( z ) dz + ˆ x ρ q ( z ) dz + ˆ x ρ q ( x ) dx = ˆ x ρ q ( z ) dz + ˆ xρ q ( x ) dx ≡ E ( x ) + E . (27)Expression (27) indicates the unsurprising fact that in a periodic system the initial point of integration x canbe chosen arbitrarily regardless of the functional form of E ( x ) .To compute correlations, we use (27), with E ( x ) = (cid:104) E ( x ) (cid:105) + ˜ E ( x ) , to find C E ( x, y ) = (cid:104) ˜ E ( x ) ˜ E ( y ) (cid:105) = (cid:104) ˜ E ˜ E (cid:105) + (cid:104) ˜ E ˜ E ( y ) (cid:105) + (cid:104) ˜ E ( x ) ˜ E (cid:105) + (cid:104) ˜ E ( x ) ˜ E ( y ) (cid:105) . (28)Using ρ q = 1 − ρ e = − ˜ ρ e , we find C E = (cid:104) ˜ E ˜ E (cid:105) = ˆ zdz ˆ w (cid:104) ˜ ρ e ( z )˜ ρ e ( w ) (cid:105) dw ; (29)we also find C E ( x ) = (cid:104) ˜ E ( x ) ˜ E (cid:105) = ˆ wdw ˆ x (cid:104) ˜ ρ e ( w )˜ ρ e ( z ) (cid:105) dz, (30)and similarly for C ( y ) , and (e.g., for x > y ), C E ( x, y ) = ˆ x dz ˆ y (cid:104) ˜ ρ e ( z )˜ ρ e ( w ) (cid:105) dw. (31)Putting these together, we have C E ( x, y ) = ˆ zdz ˆ wC ( z, w ) dw + ˆ wdw ˆ x C ( w, z ) dz + ˆ zdz ˆ y C ( z, w ) dw + ˆ x dz ˆ y C ( z, w ) dw. (32)For the special δ -function case of Eq. (22), a substitution into (32) yields C E ( x, y ) = 1 N p (cid:20) min ( x, y ) − xy + x ( x − y ( y − (cid:21) , (33)which is extended outside < x, y < to be periodic in both arguments. Equation (33) can be cast into theform C E ( x, y ) = 12 N p (cid:20) −| x − y | + ( x − y ) + 16 (cid:21) , (34)8howing explicitly the translational invariance C E ( x, y ) = C E ( x − y ) and thus independence of the initial pointof integration, x , in addition to the symmetry C E ( − x ) = C E ( x ) . The first term in Eq. (33), ∝ min ( x, y ) , isthe Brownian motion result, obtained by assuming a Poisson probability distribution of particle numbers in eachdifferential region. The sum of the first two terms represents the Brownian bridge [32], a random walk withnegative correlations that force the electric field to be equal at both ends E (0) = E (1) = 0 ; the physical originof this condition is the net neutrality of the plasma ´ ρ q ( x ) dx = 0 . The complete result (33) using the zeropotential assumption (25) can be identified as the Ornstein-Uhlenbeck bridge [27]. The three different cases arediscussed in more detail in Appendix A.The electric field correlations C E ( x − y ) in Eq. (34) are plotted in Fig. 2 for a fixed value of y , showinga maximum at a cusp at x = y , with C E being negative over about three parts and positive over about twoparts of the range of x − y . Also shown in Fig. 2 is the covariance matrix for the Poisson case (random walk)and for the Brownian bridge, neither having translational invariance, both with cusps at x = y . Notice thatfor the Poisson case C E (1 , y ) (cid:54) = C E (0 , y ) . For the Brownian bridge case, we have C E (1 , y ) = C E (0 , y ) = 0 . Itis easy to show from Eq. (33) that for periodic boundary conditions C E for the Ornstein-Uhlenbeck bridge ismaximal at the cusp at x = y (and at every periodic image), it has a smooth minimum exactly between themaxima, and has continuous derivatives at the endpoints x = 0 , . Notice that the Brownian bridge correlationsare positive, but less than those of the Poisson case for x > y , whereas correlations for the Ornstein-Uhlenbeckcase are negative over an appreciable region of ( x, y ) and are significantly smaller in magnitude. Finally, notethe aperiodic behavior of C E for the random walk and the cusps at x = 0 and x = 1 for the Brownian bridge.We have also computed C E ( x − y ) for ˆ K equal to the linear tent function kernel, which is the convolution of twoboxcar kernels of finite width h (rather than ˆ K ( x ) = δ ( x ) ). The result for h = 0 . is shown shown in Fig. 2 asthe “smooth” Ornstein-Uhlenbeck bridge.The implications of the lower values of the correlations to particle methods is as follows. Smaller variances( x = y ) are desirable because they correspond to lower noise level. For x (cid:54) = y , the lower level of the magnitudeof the correlations is expected because of the property ´ C E ( x, y ) dx = 0 . Also, in a physical system, i.e., in thelimit of large N p , such correlations vanish; therefore lower C E ( x, y ) is expected to improve the fidelity of thenumerical results. x C E ( x , y = . ) Random WalkBrownian bridgeOU bridge, cuspOU bridge, smooth
Figure 2: Comparison of the electric field covariance matrix C E ( x, y ) , with y = 0 . , for the random walk, theBrownian bridge, and the Ornstein-Uhlenbeck bridge. The black curve is the smooth Ornstein-Uhlenbeck bridgewith a linear (tent) kernel of width h = 0 . . Note that only the two Ornstein-Uhlenbeck bridge cases are periodicand have C E significantly reduced relative to the other two processes.Recall that ´ C ( x, y ) dx = 0 [cf. Eq. (23)]. A similar general result can be derived for the electric fieldcorrelations C E . Indeed, we have ˆ C E ( x, y ) dy = (cid:28) ˜ E ( x ) ˆ ˜ E ( y ) dy (cid:29) = 0 , (35)9hich vanishes because of the relation (25). In particular, (35) can be verified by a direct calculation for thespecial case of the covariance matrix (33) (or (34)). As noted for the density covariance matrix, the covariance C E also has an eigenfunction with eigenvalue zero, namely ´ C E ( x, y ) u ( y ) dy = 0 for u ( y ) = 1 ; we will returnto this point in Sec. 5. In this section we discuss a statistical study of non-uniform density distributions. We now use the more generalterm “error” or “statistical error” instead of “noise,” as we will show that noise is only part of the total error,characterized by the variance , the other important contribution being the bias . Let us evaluate the mean-square difference (error) Q (henceforth simply error; the actual error can be calculatedas √ Q ;) between the estimated density, ρ e ( x ) , and the true density ρ ( x ) , where now ρ ( x ) is not assumed to beconstant. The quantity Q is Q = (cid:68)(cid:0) ρ e ( x ) − ρ ( x ) (cid:1) (cid:69) = (cid:28)(cid:16) ρ e ( x ) − (cid:104) ρ e ( x ) (cid:105) + (cid:104) ρ e ( x ) (cid:105) − ρ ( x ) (cid:17) (cid:29) , (36)where we remind the reader that (cid:104) f (cid:105) is given by Eq. (10) or (11). This quantity equals (omitting the argument x for clarity) Q = (cid:68) ( ρ e − (cid:104) ρ e (cid:105) ) (cid:69) + 2 (cid:10) ( ρ e − (cid:104) ρ e (cid:105) ) ( (cid:104) ρ e (cid:105) − ρ ) (cid:11) + (cid:68) ( (cid:104) ρ e (cid:105) − ρ ) (cid:69) . (37)We recognize that the factor ( (cid:104) ρ e ( x ) (cid:105) − ρ ( x )) in the middle and third terms is not a random variable; since theother factor in the middle term is zero, we find Q = Q + Q , (38)with Q = (cid:10) ρ e (cid:11) − (cid:104) ρ e (cid:105) , (39) Q = (cid:0) (cid:104) ρ e (cid:105) − ρ (cid:1) . (40)To proceed, we go back to the Taylor expansion Eq. (13), noting that in a non-uniform density ρ (cid:48)(cid:48) ( x ) (cid:54) = 0 ,and write it as (cid:104) ρ e ( x ) (cid:105) = ρ ( x ) + B ( x ) + O ( h ) (41)with B ( x ) = h ρ (cid:48)(cid:48) ( x )2 ˆ ζ K f ( ζ ) dζ . (42)The quantity B ( x ) is called the statistical bias . Since ´ (cid:104) ρ e ( x ) (cid:105) dx = 1 , the bias satisfies ´ B ( x ) dx = 0 , consistentwith the periodic boundary conditions on ρ ( x ) .For the first term in Q we have (cid:104) ( ρ e ( x )) (cid:105) = (cid:88) µν q µ q ν (cid:104) K ( x − ξ µ ) K ( x − ξ ν ) (cid:105) = (cid:88) µ q µ ˆ K ( x − ξ ) ρ ( ξ ) dξ + (cid:88) µ (cid:54) = ν q µ q ν (cid:18) ˆ K ( x − ξ ) ρ ( ξ ) dξ (cid:19) , (43)where we have again split the sums into diagonal terms ( µ = ν ) and off-diagonal terms ( µ (cid:54) = ν ). Again, changingvariables and Taylor expanding, we find (cid:104) ( ρ e ( x )) (cid:105) = ρ ( x ) N p h ˆ K f ( ζ ) dζ + O (cid:18) hN p (cid:19) + (cid:18) − N p (cid:19) (cid:104) ρ e ( x ) (cid:105) ≈ ρ ( x ) N p h ˆ K f ( ζ ) dζ + (cid:104) ρ e ( x ) (cid:105) − N p ( ρ ( x ) + B ( x )) = ρ ( x ) N p h ˆ K f ( ζ ) dζ + (cid:104) ρ e ( x ) (cid:105) − N p ρ ( x ) , (44)10sing Eq. (41) and neglecting terms of first and higher orders in h/N p . Note that the − ρ ( x ) /N p term in (44)does not arise from the Taylor expansion (13). For the purpose of the present argument we neglect that termsince ρ ( x ) is of order one and N p is typically a large number in particle simulations. However, recall that this isthe same factor responsible for the negative correlations in Sec. 3, where although small, it had a non-negligiblecumulative effect; we will revisit its importance in Sec. 6.We find that the (cid:104) ρ e ( x ) (cid:105) terms cancel in Eqs. (39) and we are left with Q = ρ ( x ) N p h ˆ K f ( ζ ) dζ. (45)We also have from Eq. (40), (41), and (42) Q = B ( x ) = (cid:18) h ρ (cid:48)(cid:48) ( x )2 ˆ ζ K f ( ζ ) dζ (cid:19) . (46)Eqs. (45), (46) give Q = Q + Q = ρ ( x ) N p h ˆ K f ( ζ ) dζ + (cid:18) h ρ (cid:48)(cid:48) ( x )2 ˆ ζ K f ( ζ ) dζ (cid:19) . (47)The first term, Q = V , is the variance (the diagonal terms in the covariance matrix) and the second term, Q ,is the square of the bias , Q = B . (Note that the addition of the bias to ρ ( x ) in Eq. (41) is analogous to thesmoothing obtained by diffusion of the density over a time interval t , ρ ( x ) → (cid:0) Dt∂ x + · · · (cid:1) ρ ( x ) , where D isa diffusion coefficient and h ´ ζ K f ( ζ ) dζ/ → Dt .) Writing C = ˆ / − / K f ( ζ ) dζ , C = ˆ / − / ζ K f ( ζ ) dζ , (48)we have Q = V + B = ρ ( x ) C N p h + ρ (cid:48)(cid:48) ( x ) C h . (49)Clearly the factors C and C are related to the kernel shape, whereas the kernel width is represented by h .The interpretation of the two contributions in the result Eq. (49) is as follows: the bias is an error caused byestimating the spatially varying density ρ ( x ) using a kernel of width h , i.e., it is a finite size particle effect; thevariance is an error (noise) due to the finite number of particles. The balance between the two effects is reachedwhen ρ ( x ) N p h ∼ ρ (cid:48)(cid:48) ( x ) ρ ( x ) h ∼ (cid:18) hl (cid:19) , (50)where we have defined the density gradient length scale l = (cid:112) ρ ( x ) / | ρ (cid:48)(cid:48) ( x ) | and have assumed the shapecoefficients C and C are of order unity. We see that the bias error term dominates for h large comparedto l (more smoothing of the density); more specifically, when (cid:18) lh (cid:19) (cid:28) ρ ( x ) N p h ≡ N h . (51)One recognizes the product N h = ρ ( x ) N p h as the typical number of particles within the kernel width h . Thevariance error dominates when the opposite inequality holds. The condition (50) will be revisited in Sec. 6 wherenumerical examples are presented.For a more quantitative description, optimizing over h for fixed x , we find a minimum at h = h opt = (cid:18) ρ ( x ) C N p ρ (cid:48)(cid:48) ( x ) C (cid:19) / , (52) Q min = 54 (cid:32) ρ ( x ) | ρ (cid:48)(cid:48) ( x ) | / C C / N p (cid:33) / , (53) Q (cid:48)(cid:48) ( h opt ) = 5 (cid:18) ρ ( x ) | ρ (cid:48)(cid:48) ( x ) | C C N p (cid:19) / . (54)11 ernel C C (cid:16) C C / (cid:17) / (cid:0) C /C (cid:1) / Boxcar /
12 0 .
370 2 . Linear (tent) / /
24 0 .
353 3 . Quadratic /
20 1 /
36 0 .
356 4 . Trapezoidal / /
108 0 .
350 3 . Epanechnikov / /
20 0 .
349 3 . Table 2: The values of the coefficients C and C for the fundamental kernels in Table 1 and Fig. 1, includingthose used in Fig. 3. The quantity in column is the factor appearing in Q min . The quantity in column is thefactor appearing in h opt as well as in the width W Q from Eq. (58). Note that column varies little between thekernels, but column varies by almost a factor of two.This process of minimizing Q is called bias-variance optimization [25, 26]. We see that the very factor that hasmade particle methods so useful—the finite size of computational particles—is not without its drawbacks, leadingto the bias error in the density estimation. However, our result provides a guideline for taking advantage of thisfactor as it varies oppositely to the other error contribution, that of the variance (noise). Thus we arrive at thetrade-off between bias and variance error embodied in the BVO process just described.Eqs. (52) and (53) suggest that the optimal value of h depends on x . A reasonable alternative is to let ρ ( x ) → ´ ρ ( x ) dx = 1 and ρ (cid:48)(cid:48) ( x ) → ´ ρ (cid:48)(cid:48) ( x ) dx , i.e., to integrate Eq. (49), leading to the mean integratedsquare error result h = h opt , av = (cid:32) C N p (cid:0) ´ ρ (cid:48)(cid:48) ( x ) dx (cid:1) C (cid:33) / , (55) Q min , av = 54 (cid:32) (cid:0) ´ dx | ρ (cid:48)(cid:48) ( x ) | (cid:1) / C C / N p (cid:33) / , (56) Q (cid:48)(cid:48) ( h opt , av ) = 5 (cid:32) ( ´ ( ρ (cid:48)(cid:48) ( x )) dx ) / C C N p (cid:33) / . (57)If | ρ (cid:48)(cid:48) ( x ) | does not vary too much, it is possible to take advantage of the fractional power in Eq. (56) to use akernel of width h ≈ h opt , av throughout the whole simulation domain. This is especially important in cases inwhich a choice is made to have a fixed relation between the kernel width h and a uniform grid spacing ∆ .Note the dependence of the quantities in Eqs. (52), (53), and (54) on N p , namely h opt ∝ N − / p , Q min ∝ N − / p and Q (cid:48)(cid:48) ( h opt ) ∝ N − / p . This shows that h opt is quite insensitive to the number of particles. It is interesting tonote that Q min has a slightly weaker scaling that the usual ∝ N − p scaling of the variance alone, implying scaling /N / p vs. / (cid:112) N p for the error ∝ Q / .Kernels with compact support are typically used in particle simulations, for computational efficiency. Amongall kernels with compact support, with width equal to one, and having ´ K f ( ζ ) dζ = 1 , the Epanechnikov kernelminimizes the factor C C / in Q min [31, 33]. However, the factor ( C C / ) / in Q min varies little betweendifferent kernels, so that the kernel shape has little influence on Q min .A plot of the Q vs. h [cf. Eq. (56)] is shown in Fig. 3, using the mean integrated square error approximationand ´ ρ ( x ) dx = ´ ρ (cid:48)(cid:48) ( x ) dx = 1 . The three curves Q ( h ) correspond to the boxcar, quadratic, and Epanechnikovkernels, defined in Table 1. The coefficients C and C , calculated from Eq. (48), are given in the first twocolumns of Table 2. The number of particles is taken to be N p = 10 . The range of h is chosen so that thesections dominated by variance (small h ) and by bias (large h ) are clearly seen, as well as the intermediate valueswhere a minimum is attained. Note the shape dependencies h opt , av ∝ ( C /C ) / , Q min ∝ ( C C / ) / . Thewidth of the minimum of Q ( h ) is proportional to W Q ∝ ( Q min /Q (cid:48)(cid:48) ( h opt )) / ∝ ( C /C ) / , (58)which is the same factor appearing in the expressions for h opt . These combinations are also listed in Table 2.The values in the third column, i.e., ( C C / ) / , confirm that the shape has a minimal effect on Q min and isconsistent with the slightly lower minimum of the Epanechnikov kernel compared to the other two kernels inFig. 3, as discussed above. It is also easy to see that the location of h opt for the three curves in Fig. 3 is consistentwith the values in the last column of that table; for example, note that the ratio of the values of h opt , av for the12 .00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 h Sq u a r e d e rr o r , Q = V + B BoxcarQuadr. splineEpanechnikov
Figure 3: A sketch of the error, Eqs. (49), with ρ ( x ) and ρ (cid:48)(cid:48) ( x ) averaged, as a function of h for the boxcar,Epanechnikov and quadratic spline kernels, showing a slightly lower minimum at h opt for the Epanechnikov.Also, note the higher values of h opt and the broader minima for the Epanechnikov and quadratic spline kernels,in agreement with the results summarized in Table 2.Epanechnikov kernel to the boxcar kernel, equal to . / . ≈ . . Also, note that the last column in Table 2is the factor entering in the width W Q .An important feature of the error curve is the relatively broad minimum, which means Q min does not varysignificantly over a relatively large range of h values around its optimal value. This has the practical consequencethat when ρ ( x ) and ρ (cid:48)(cid:48) ( x ) have a modest variation over the simulation domain, a fixed kernel of width close to h opt , av still provides a near optimal density estimate; hence, we have another justification for using the averagedquantities, Eqs. (55) and (56). Fig. 3 shows that the minimum for the Epanechnikov kernel is broader by lessthan a factor than that for the boxcar kernel, while for the quadratic spline that factor is closer to ; bothof these observations are consistent with the values in the last column of Table 2. Quantitative comparison ofEqs. (52)–(56) against numerical simulations is done in Sec. 6, which relaxes the approximations of this analysis.Next we illustrate how our results can be applied to algorithms of practical importance, i.e., including griddiscretization. So far we have obtained results in continuous spatial variables. For numerical purposes, we need to performgrid discretization. Clearly, there is no universal discretization and different problems may benefit from differentdiscretizations. For illustration purposes, in this section we choose a particular finite difference discretizationthat is commonly used in electrostatic PIC algorithms [8] but repeating the analysis for other discretizations isstraightforward, including the use of finite elements.
An important step in obtaining a complete particle algorithm is the connection of Lagrangian particles with aEulerian grid. This connection is given by a charge deposition rule and traditionally done with so-called splinefunctions [8]. Although spline functions of varying degree of smoothness and width are available, they have thefollowing two limitations: (i) their width is an integer number of cell widths, i ∆ , i = 1 , , . . . ; and (ii) theirwidth and smoothness are strictly related, with smoother particles being wider. The smoothness of a particlebecomes important when force interpolation from the grid to the particle position, especially when a particlecrosses cell boundaries. Since in this work we are not addressing the full PIC cycle, we emphasize the importanceof particle width over smoothness. As well, if one were to use particles arbitrarily related to the grid spacing ∆ to minimize noise and error, per our theoretical developments, the above restrictions may present a drawback.For example, when high grid resolution is desired (small ∆ ) while h opt is (relatively) large, that would require aparticle that spans a large number of cells. If splines are used, they would be of high order and computationallyexpensive because of the larger number of floating point operations associated with high order polynomials. In13 article shape DefinitionBoxcar (NGP) S B ( x ) = (cid:40) , (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ otherwise .Linear spline S L ( x ) = (cid:40) − (cid:12)(cid:12) x ∆ (cid:12)(cid:12) , (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ otherwise .Quadratic spline S Q ( x ) = − (cid:0) x ∆ (cid:1) , (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ / (cid:0) − (cid:12)(cid:12) x ∆ (cid:12)(cid:12)(cid:1) , / ≤ (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ / otherwise .Trapezoidal S T ( x ) = , (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ / (cid:0) − (cid:12)(cid:12) x ∆ (cid:12)(cid:12)(cid:1) , / ≤ (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ / otherwise , Table 3: Examples of particle shapes. These shapes are similar to the kernels of Table 1 and Fig. 1, but haveinteger valued cell width and satisfy the sum rule, Eq. (9).fact, this is why practitioners rarely go beyond fourth order spline functions. The advantage of being able tochoose separately particle smoothness and width becomes obvious.In this section we address the relaxation of the two limitations discussed above, those associated with thesmoothness and the width of particle shapes. The former was discussed in Ref. [12], where the smoothness ofparticle shapes was decoupled from its width. An example of a cubic particle shape depositing charge on threegrid points (same as the quadratic spline) was given therein; by the same method, for example, one could devisea quadratic particle wider than three cells, etc. The key element that allows this generalization is the distinction between the kernel K ( x ) and the particle shape (factor) S ( x ) , alluded to in Sec. 2. Therefore, before proceedingto address the limitation associated with the particle width, we return to a discussion of the difference between K ( x ) and S ( x ) .It is the sum rule property that distinguishes a kernel from a particle shape, see Eq. (9). We will require thata particle shape satisfies the sum rule, whereas we will not impose this requirement on a kernel. That propertystates that the sum of the fractions of the computational particle’s charge deposited on the grid sum exactlyto the charge carried by the particle, and this is true at any (continuously varying) particle position ξ . As aconsequence, the total charge of the system after being deposited on the grid is also conserved. For example, fordensity that integrates to unity on a uniform grid with spacing ∆ , we have ´ ρ e ( x ) dx = (cid:80) N p µ =1 q µ = 1 . Considerthe amount of charge a single particle deposits on the grid point x i . We use a charge deposition rule based on aparticle shape S ( x ) , which gives for the fraction of that charge q µ ∆ S ( x i − ξ µ ) . We find that the density on thegrid point x i due to depositing the charge from all the particles is ρ e ( x i ) ≡ ρ e,i = N p (cid:88) µ =1 q µ S ( x i − ξ µ ) . (59)For particles with equal charges, q µ = 1 /N p , we obtain N g (cid:88) i =1 ∆ ρ e,i = N g (cid:88) i =1 ∆ N p (cid:88) µ =1 q µ S i ( ξ µ ) = N p (cid:88) µ =1 N p N g (cid:88) i =1 ∆ S i ( ξ µ ) = N p (cid:88) µ =1 N p = 1 . (60)We have defined S k ( ξ µ ) ≡ S ( x k − ξ µ ) and invoked the sum rule, Eq. (9). Thus, the sum rule implies that the totalcharge assigned to the grid is preserved. If a particle shape does not satisfy the sum rule (i.e. it obeys conditions(3)–(7) for a kernel but not the sum rule), the lack of exact charge conservation would allow E (1) (cid:54) = E (0) . Ofcourse, the total charge of the system is always conserved as long as no association with a computational grid ismade, as discussed in Sec. 2. Examples of particle shapes (satisfying the sum rule) are given in Table 3. The firstthree are familiar from Ref. [8]; the last (trapezoidal) particle shape is discussed below. The charge depositionon the grid point x i , associated with each of the particles in Table 3, is found by the substitution x = ξ µ − x i ,for i = 1 , , . . . , N g .We now describe the generalization associated with particle width : a particle shape satisfying the sum rule(9) is not required to have a width equal to an integer number of grid cells. In fact, such a shape can have14almost – see below) completely arbitrary width relative to the grid. To see this, consider a known “primary”particle shape that satisfies properties (3)–(7) and the sum rule (9), say S ( x ) , with (cid:80) i ∆ S ( x i − ξ ) = 1 for anyvalue of ξ ∈ [0 , . Then for an arbitrary kernel K ( x ) satisfying properties (3)–(7), we perform the convolution S ( x ) = ˆ K ( y ) S ( x − y ) dy. (61)The so-obtained new particle shape S ( x ) satisfies the normalization (3): ˆ S ( x ) dx = ˆ K ( y ) dy ˆ S ( x − y ) dx = ˆ K ( y ) dy = 1 . (62)The sum rule (9) is also easily verified: N g (cid:88) i =1 ∆ S ( x i ) = ˆ K ( y ) (cid:88) i ∆ S ( x i − y ) dy = ˆ K ( y ) dy = 1 . (63)It is easy to verify that properties Eq. (4)–(7) are inherited by S ( x ) as well. The width of S ( x ) equals to thesum of the widths of S ( x ) and K ( x ) . The procedure just described allows a kernel K ( x ) of arbitrary width;we conclude that this construction allows one to generate arbitrary width particle shapes satisfying the sumrule, including such that are non-integer number of cells wide. The only condition on the width of S ( x ) is thatit cannot be less than the width of S ( x ) , hence the qualifier “almost” above; this is usually not a limitation.We stress that if a particle shape S ( x ) is obtained from a kernel K ( x ) and another particle shape, S ( x ) , theirfunctional forms are different (in addition to their widths being different).We note that the convolution described by Eq. (61) is the easiest way to obtain a particle shape that satisfiesthe sum rule. However, this is a sufficient but not necessary condition. As well, choosing a familiar particleshape obeying the sum rule as a primary, S , is the easiest way to ensure S ( x ) satisfies the sum rule; again,this choice is sufficient but not necessary. In other words, it is possible that other methods of obtaining particleshapes that satisfy the sum rule exist.The examples in Table 3 satisfy the sum rule and can either be used as primary to obtain other particleshapes or directly in a simulation. In fact, all particles from Table 3 can in turn be obtained by the convolutionformula (61) as follows: the boxcar shape is the convolution of itself with a delta-function kernel; the linearshape is the convolution of the boxcar shape with a boxcar kernel of width ∆ ; the quadratic spline shape is theconvolution of the boxcar shape with a linear kernel (tent function) of width ; and finally, the trapezoidalshape is the convolution of the boxcar shape and a boxcar kernel of width .An example of a particle shape that is non-integer number of cells wide is given next. Consider the convolutionof the boxcar particle shape of width ∆ from Table 3 with the boxcar kernel of width δ with ≤ δ ≤ ∆ , K δ ( x ) = 1 δ (cid:40) , (cid:12)(cid:12) xδ (cid:12)(cid:12) ≤ otherwise . (64)The convolution formula (61) gives the following trapezoidal particle shape, of width ∆ + δ : S ∆+ δ ( x ) = ˆ K δ ( y ) S B ( x − y ) dy = 1∆ , ≤ | x | ≤ ∆ − δ δ ∆+ δ − | xδ | , ∆ − δ ≤ | x | ≤ ∆+ δ , otherwise . (65)The particle shape (65) transforms into the usual ∆ -wide boxcar shape in the limit δ → and into the usual -wide linear shape in the limit δ → ∆ . Note that ´ S ∆+ δ ( x ) dx = 1 , as expected. The fractional width particleis illustrated in Fig. 4. The charge deposition rule is found by substitution of ξ − x i into (65), where x i is thenearest grid point to the particle position and | ξ − x i | ≤ ∆ / . The result is given in Table 4 and the constructionby means of a convolution assures that the sum rule is satisfied for arbitrary values of δ .We remark that the dependence of the particle shape (65) on the fractional width δ is not a scalingtransformation via a fundamental kernel such as changing h in (8). Instead, this is a parametric shape transformation (with parameter δ ). Nevertheless, changing the shape via δ also changes the support of S ∆+ δ and provides anothermeans of attaining the optimal width, h opt = ∆ + δ . 15 .6 0.4 0.2 0.0 0.2 0.4 0.6 x S + ( x ) + Figure 4: Illustration of the fractional width particle shape, Eq. (65), with δ/ ∆ = 0 . and h = 1 . . Becausethis quantity is a convolution of two functions, one of which satisfies the sum rule in Eq. (9), it satisfies the sumrule for an arbitrary value of δ/ ∆ , and therefore for an arbitrary value of h/ ∆ < . Charge deposition rule Range S i − ( ξ ) = 0 S i ( ξ ) = ≤ | ξ − x i | ≤ ∆ − δ S i +1 ( ξ ) = 0 S i − ( ξ ) = δ (cid:104) − ∆ − δ − ( ξ − x i ) (cid:105) S i ( ξ ) = δ (cid:104) ∆+ δ + ( ξ − x i ) (cid:105) − ∆2 ≤ ξ − x i ≤ − ∆ − δ S i +1 ( ξ ) = 0 S i − ( ξ ) = 0 S i ( ξ ) = δ (cid:104) ∆+ δ − ( ξ − x i ) (cid:105) ∆ − δ ≤ ξ − x i ≤ ∆2 S i +1 ( ξ ) = δ (cid:104) − ∆ − δ + ( ξ − x i ) (cid:105) Table 4: Charge deposition rule corresponding to the fractional particle shape (65).The existence of the fractional width particle shape (65) was noted in Ref. [18] and has been previously knownin particle hydrodynamics. Its derivation, however, has been based on area weighting arguments and not on themore general convolution method given by (61). Again, the fractional width trapezoidal shape (65) is not thesame as the trapezoidal shape listed in Table 3.Another approach to obtaining particle shapes of arbitrary support is based on the finite element method ofdiscretizing a system of equations [34]. In particular, in particle algorithms based on a variational principle [4,12],the convolution method of obtaining particle shapes emerges naturally [12] and the role of the above primaryparticle shape S , which provides the connection to the grid, is taken by finite element basis functions. Forexample, tent functions of width and height (linear Lagrange finite elements, which are basically the sameas the linear spline shape function in Table 3 except with a different amplitude), Ψ i ( x ) = Ψ( ξ − x i ) , have theproperty that at any x (including at grid points x i ) N g (cid:88) i =1 Ψ i ( x ) = 1 . (66)16xcept for a factor ∆ , this is the sum rule that we require of particle shapes. That is, to obtain a particleshape S ( x ) satisfying the sum rule (9), we perform the convolution S ( x ) = ´ K ( y )Ψ( x − y ) dy/ ∆ with thefinite element. (Note that translation invariance is also satisfied). The unit normalization follows from the finiteelement property ´ (1 / ∆)Ψ i ( x ) dx = 1 , i = 1 , , . . . , N g , and the kernel normalization (3). The shapes (top tobottom) from Table 3 may also be obtained by the finite element method from convolutions.To summarize, relaxing the limitations associated with particle smoothness and width allows one to deviseparticle shapes that are both computationally efficient and suitable to take advantage of the BVO guidelinesdiscussed in the previous section, as well as assuring that the sum rule is satisfied. We consider a uniform grid in x on [0 , with vertices at x i = ∆ i , i = 0 , · · · , N g and grid spacing ∆ = 1 /N g .Returning to ρ ( x ) = 1 for simplicity, we define the estimated density at cell centers , x i +1 / = ( i + 1 / , as ρ e,i +1 / = ρ e ( x i +1 / ) , ≤ i ≤ N g − . Note that given a particle shape, e.g., from Table 3, charge deposition on cell centers amounts to simplysubstituting x = ξ − x i +1 / in S ( x ) . Recall that for uniform density we have (cid:104) ρ e,i +1 / (cid:105) ≡ (cid:10) ρ e ( x i +1 / ) (cid:11) = 1 ;then the discrete approximation to ´ ρ e ( x ) dx = 1 is N g − (cid:88) i =0 ρ e,i +1 / ∆ = (cid:88) i,µ q µ S i +1 / ( ξ µ )∆ = (cid:88) µ q µ (cid:88) i ∆ S ( x i +1 / − ξ µ ) = 1 . (67)For the covariance matrix we have, from Eqs. (19), (21) C i +1 / ,j +1 / ≡ C ( x i +1 / , x j +1 / ) = C d,i +1 / ,j +1 / − N p = 1 N p (cid:20) ˆ S ( x i +1 / − ξ ) S ( x j +1 / − ξ ) dξ − (cid:21) , (68)the last equality due to the sum rule, Eq. (9), and (cid:80) µ q µ = 1 . The discrete analog of Eq. (23) is the conditionthat the sum over each row (or column) of the covariance matrix is zero; we have (cid:88) j ∆ C i +1 / ,j +1 / = 1 N p (cid:88) j ∆ ˆ S ( x i +1 / − ξ ) S ( x j +1 / − ξ ) dξ − = 1 N p ∆ (cid:88) j ˆ S ( x i +1 / − x j +1 / ) − = 0 , (69)where ˆ S is the convolution of S with itself. Since S satisfies the sum rule and assumptions (3)–(7), so does ˆ S ,therefore, the quantity in Eq. (69) sums to zero and we obtain the analogous discrete result as in Eq. (23). Asin the continuous case, this identity says that the vector u i = 1 is associated with zero eigenvalue, implying thatthe covariance matrix is singular. In the next section we will show exact calculations of these covariance matrixelements for specific particle shapes.The negative correlations of Eq. (21) also appear in Eq. (68). In order to understand these correlations, letus compare a case in which we pick the number N i +1 / of particles in the cells independently and identicallydistributed (iid) from a Poisson distribution with parameter λ = N ppc = N p ∆ , with mean λ and variance λ . Here, N ppc is the expected number of particles per cell. Again, recall that in this discussion, we are assuming ρ ( x ) = 1 ;thus the mean of ρ e,i +1 / = N i +1 / /N ppc is unity and its variance is /N ppc . Indeed, from the iid assumption,the off-diagonal terms are zero and using the property of the variance, Var ( ρ e,i +1 / ) = Var ( N / ) /N ppc , weobtain C i +1 / ,j +1 / = λN ppc δ ij = 1 N ppc δ ij . (70)Referring to Eqs. (8), (68), for a kernel of width h ∼ ∆ , the diagonal (the variance) is ∼ / ∆ N p = 1 /N ppc ,comparable to the dominant part of the diagonal in Eqs. (68), (69). However, the negative correlations − /N p in Eqs. (21) and (68) are not contained in the Poisson model. These negative correlations are traced to the fact17hat the total number of particles N p in Eqs. (21) and (68) is fixed. In a particle code, the total number ofparticles can be assumed fixed at each time step. (The particle number may also be constant throughout thesimulation for certain type of boundary conditions such as periodic, for example.) The fixed number of particlesis in contrast to the Poisson case in which the expected number of particles per cell is N ppc (expected total numberof particles = N p ). Intuitively, when the total number of particles is fixed, if one cell has more than the expectednumber of particles N ppc , other cells must necessarily have fewer particles, leading to negative correlations. (Seealso Ref. [35] where negative correlations between numbers of particles in different cells were obtained workingfrom the multinomial distribution.) We will discuss in a later section the effect of these negative correlations onthe calculation of the electric field. C i +1 / ,j +1 / examples In this section we present covariance matrix calculations with particle shapes from Table 3. The simplest particleshape is the boxcar (top-hat) function. For i (cid:54) = j the overlap integral in Eq. (68) is zero and we find C i +1 / ,j +1 / = − N p ( j (cid:54) = i ) giving C i +1 / ,j +1 / = 1 N ppc δ ij − N p . (71)The first term in Eq. (71) is equal to the value in Eq. (70), and the second term is recognized as the negativeterm of Eq. (68). We conclude indeed that the − /N p term is due to the fact that the total number of particlesis fixed. The condition (69) is obviously satisfied in this example.For the linear particle shape we find C i +1 / ,j +1 / = N p (cid:2) ∆ ´ S ( ξ ) dξ − (cid:3) = N p (cid:2) − (cid:3) = N ppc − N p ( j = i ) , N p (cid:2) ∆ ´ S ( ξ − S ( ξ ) dξ − (cid:3) = N ppc − N p ( j = i ± , − N p otherwise. (72)Note that the condition (69) holds for this example as well. This condition highlights the importance of the − /N p correlations and shows how the variances ( i = j terms) decrease as the particle width increases.For the quadratic particle shape we have C i +1 / ,j +1 / = N ppc − N p ( j = i ) N ppc − N p ( j = i ± N ppc − N p ( j = i ± − N p otherwise. (73)The property (69) is again clearly seen for all these cases. Note again that the variances, i.e. the i = j terms,decrease further as the particle widths increase, with their values distributed to more neighboring bands, withthe negative term − /N p on diagonal as well as off-diagonal terms. In this section we analyze the statistical properties of the electric field on a grid, again assuming uniform density, ρ ( x ) = 1 . The continuous form of these relations was presented in Sec. 3.2. In the previous section, the densitywas assumed to reside at cell centers ρ e,i +1 / , i.e., at x i +1 / . If the density and the electric field are taken tobe staggered, with E i residing at vertices x i for i = 0 , · · · , N g , to obtain a second order accurate differencingscheme, we write E i +1 − E i = ∆ ρ q,i +1 / , (74) This discretization can also be derived using finite elements. .0 0.2 0.4 0.6 0.8 1.0 x E l e c t r i c f i e l d RWBBOU
RW: Random walkBB: Brownian bridgeOU: Ornstein-Uhlenbeck
Figure 5: Illustration of three types of random behavior of the electric field depending on the boundary conditions,as discussed in Appendix A. Both RW and BB have a starting point E (0) = 0 (see inset); the BB case hasadditionally E (1) = 0 , and the OU case is shifted down, satisfying the zero potential difference condition.where ρ q,i +1 / = 1 − ρ e,i +1 / is the charge density and, as before, assuming uniform and immobile ions of unit(neutralizing) density. This relation leads to E i = E + E ,i , where E ,i = ∆ i − (cid:88) j =0 ρ q,j +1 / . (75)The condition E = E N g ( x = 0 , x N g = 1 ) follows from N g − (cid:88) j =0 ρ q,i +1 / = 0 (76)and Eq. (67). (It is instructive to revisit the derivation of the discrete property (69): this result can be obtaineddirectly by writing (cid:80) j (cid:104) ˜ ρ e,i +1 / ˜ ρ e,j +1 / (cid:105) = (cid:104) ˜ ρ q,i +1 / (cid:80) j ˜ ρ q,j +1 / (cid:105) , which is seen to vanish by Eq. (76).) Thecondition in Eq. (25) of having zero applied potential across a period takes the discrete form N g − (cid:88) i =0 ∆ E i = ∆ N g E + N g − (cid:88) i =1 E ,i = 0 . (77)This is a condition on E ; we find E = − ∆ N g − (cid:88) i =0 i − (cid:88) j =0 ρ q,j +1 / = − ∆ N g − (cid:88) j =0 ( N g − j − ρ q,j +1 / . As in Sec. 3.2, the terms in Eq. (75) lead to four distinct terms in the covariance matrix for the noise in theelectric field, (cid:104) ˜ E i ˜ E j (cid:105) = C E + C E ,i + C E ,j + C E ,ij , where C E = (cid:104) ˜ E (cid:105) , (78)19 E ,i = (cid:104) ˜ E ,i ˜ E (cid:105) , (79)similarly for C E ,j , and C E ,ij = (cid:104) ˜ E ,i ˜ E ,j (cid:105) . (80)For the first of these we find C E = ∆ N g − (cid:88) j =0 ( N g − j − N g − (cid:88) k =0 ( N g − k − C j +1 / ,k +1 / , (81)where C j +1 / ,k +1 / = (cid:104) ˜ ρ q,j +1 / ˜ ρ q,k +1 / (cid:105) . For the next term (and similarly for C E ,j ) we have C E ,i = − ∆ i − (cid:88) j =0 N g − (cid:88) k =0 ( N g − k − C j +1 / ,k +1 / . (82)Finally, the last term is C E ,ij = ∆ i − (cid:88) k =0 j − (cid:88) l =0 C k +1 / ,l +1 / . (83)For this covariance matrix we find C E ,ij = 1 N p (cid:0) ∆ min( i, j ) − ∆ ij (cid:1) = 1 N p (min( x i , x j ) − x i x j ) . The first term is the value for the Poisson case and the second is the Brownian bridge contribution from theconstant correlations − /N p in Eq. (71). The remaining terms arise from the Ornstein-Uhlenbeck bridge; wehave C E ,i = 1 N p (cid:18) − ∆ ( N g − i + ∆ i ( i − N g ( N g − i (cid:19) → N p (cid:18) x i − x i (cid:19) , the latter limit as ∆ → with ∆ N g =1. The quantity C E ,j is computed similarly and C E = 1 N p (cid:32) ∆ ( N g − N g (2 N g − − ∆ N g ( N g − (cid:33) → N p and in the same limit, ∆ → . We see that the limiting case is in agreement with Eq. (33).The three different types of random behavior of the electric field, depending on the boundary conditions,is illustrated in Fig. 5. The random walk curve is a Brownian motion with E (0) = 0 without imposing anyboundary condition at x = 1 , the Brownian bridge curve reflects the extra boundary condition E (0) = E (1) = 0 ,and the Ornstein-Uhlenbeck curve satisfies both E (0) = E (1) = 0 and the zero potential difference condition.The latter means that if E ( x ) makes an excursion into the positive half plane, it must do so in negative halfplane as well so that E ( x ) integrates to zero – see the discussion in terms of position, velocity, and accelerationin Appendix A.In the next section we verify numerically our theoretical conclusions. Considerations dictating the choice of a charge deposition rule in a particle algorithm were discussed in Sec. 5.1,where two important factors were identified – smoothness and width. In sections 3 and 4 we have found theparticle width to be more relevant to our discussion: it reduces the variance for uniform density and minimizesthe error in a non-uniform density via the BVO process, independent of grid resolution. For discretization on agrid, we have also emphasized the importance of obeying the sum rule.20ne strategy for applying our theory in practice is to use particles that have width closest to the optimal,i.e., h = i ∆ + δ ≈ h opt with h opt given by (52), i an integer number, and ∆ chosen according to a desired gridresolution. One may further decide on a particle that is an integer number of cells wide or add the additionalcorrection δ . Using particles with exact width equal to h opt is obviously not possible in general simulations wherethe exact density and its gradients are unknown and may vary in time. However, the estimated quantities canbe used as a guideline or a lower resolution simulation may be done to probe for these and other properties ofa physical system. As discussed previously, spline functions of order i > are probably not the most efficientchoice and custom particles may be better suited. If one uses high grid resolution (small ∆ ), one may be ableto approximate h opt sufficiently well with an integer number of cells, h = i ∆ ≈ h opt ; However, one advantage ofusing a fractional width shape with a correction δ is that the same charge deposition can be used and adjusted“in real time,” depending on the values of ρ ( x ) and ρ (cid:48)(cid:48) ( x ) . This provides “fine tuning” ability, which may bepreferable to changing the type/width of particle in the course of a simulation and may help to avoid introducingundesired effects or difficulties.In the following sections, our focus will be on theory comparison and verification, which is why we will notbe concerned with the requirement of grid resolution. Instead, we will use the grid spacing to adjust the widthof the same particle shape, thus applying the scaling transform K ( x ) = K f ( x/h ) /h . When using the scalingmethod, we will not be dealing with fractional width particles; therefore, to change the width of a given particle,we will vary the grid spacing by integer numbers: for example, when N g = 15 , the width of the three-cell-widequadratic spline equals h = 3∆ = 3 /
15 = 0 . , for N g = 30 its width equals h = 3 /
30 = 0 . , etc. Recall that wework in the domain [0 , ; a range of N g (cid:39) . . . will prove to provide a sufficient range of particle widths. Ina separate set of simulations, we present results on a fixed grid but varying the fractional width of the particledefined in (65). We present as a first example numerical computations of the density covariance matrix with uniform true density ρ ( x ) = 1 on a uniform periodic grid on [0 , and we use the linear particle shape from Table 3 for our chargedeposition. Recall that the covariance matrix for this case is given by Eq. (72); however, more convenientquantities to test are the products (cid:101) C i + ,i + ≡ C i + ,i + × N ppc = 23 − ∆ , (84) (cid:101) C i + ,i + ± ≡ C i + ,i + ± × N ppc = 16 − ∆ (85)since they are independent of N ppc and N p .The theory was developed in the limit of infinite number of samples by integrating over a continuous densitydistribution but clearly, numerically we can only use a finite number of samples in averages. The present resultsaim to verify the developed theory as well as to inform us of the number of samples needed in simulations in thenext section. Because in this section we deal with uniform density, the correlations are expected to be the samefor every grid point. Using this fact, we can obtain better statistics by averaging correlations over the wholegrid, i.e., ¯ C i + ,i + = (1 /N g ) (cid:80) N g C i + ,i + for every i = 1 . . . N g . Therefore, for a number of M local samples,such averaging amounts to performing M = N g × M local local cell samples; to the end of this section we omit theover-bar. Each particle sample is drawn from a uniform distribution by drawing N p random numbers R ∈ [0 , and setting particle positions ξ µ = R , µ = 1 , , . . . N p .Simulation results on a fixed grid with N g = 25 ( ∆ = 0 . ) are listed in Table 5. Theory predicts (cid:101) C i + ,i + =0 . . . . − .
04 = 0 . . . . and (cid:101) C i + ,i + ± = 0 . . . . − .
04 = 0 . . . . . We have also intentionally taken N p M (cid:101) C i + ,i + (cid:101) C i + ,i + ± theoretical numerical theoretical numerical250 . × . × . × Table 5: Simulation results for uniform density distribution. The theoretical values have the fixed quantity ∆ subtracted. In the numerical values, this quantity is not subtracted explicitly.the product N p × M = const. in order to examine the role of number of samples versus number of particles. The21umerical values in the table agree with the theoretical values of the correlations and are most accurate (to significant figures) when using the largest number of samples M = 2 . × – the first row in Table 5. Using alarger number of particles and smaller number of samples shows good agreement as well, albeit with a somewhatlarger error in the third significant figure. We conclude that a number of samples of order should be sufficientfor comparison with theory, with three significant figures and a small error bar in the third significant figure. This set of simulations aims to compare numerical results and theory for the local error in a non-uniform density, ρ ( x ) (cid:54) = const. We compare the minimum error and optimal particle width for four different shapes, all havingthe same support h = 3∆ : the boxcar shape scaled to , the quadratic spline and the trapezoidal shape fromTable 3, and the Epanechnikov kernel scaled to , K E ( x ) = 1∆ (cid:40) (cid:16) − (cid:0) x ∆ (cid:1) (cid:17) , (cid:12)(cid:12) x ∆ (cid:12)(cid:12) ≤ otherwise . (86)These four shapes also have different smoothness, which is another point of comparison. We note that the widerboxcar shape satisfies the sum rule just as the ∆ -wide NGP shape does but deposits constant amount of charge onthree grid points instead of one. Unlike the first three shapes, the Epanechnikov kernel (86) does not satisfy thesum rule and therefore does not conserve charge on the grid, i.e., is not a particle shape in the sense of Sec. 5.1.For this reason, the Epanechnikov kernel is not recommended for practical applications; however, for our purposeof statistical calculations, charge conservation is not required. The trapezoidal shape (which does satisfy thesum rule) is also unusual and since we have not tested it in full particle simulations, we do not recommend it forpractical use at this point. The purpose of present comparisons is to (i) verify the theoretical prediction that theminimal error depends primarily on the particle width and not significantly on the specific particle shape (andsmoothness); and (ii) to verify the general dependence of Q min and h opt on the particle shape. For all simulations in this section the true density distribution is given by the periodic function ρ ( x ) = 1 + a cos(2 πmx ) , (87)where a < is a constant amplitude and m is an integer mode number. Since (87) satisfies ´ ρ ( x ) dx = 1 and ρ ( x ) ≥ for x ∈ [0 , , it is a probability density distribution. For all simulations in this section we have chosen a = 1 / , and m = 2 . We show results for two locations: x = 1 / , where ρ (1 /
2) = 3 / and ρ (cid:48)(cid:48) (1 /
2) = − π ,and x = 1 / , where ρ (1 /
3) = 3 / and ρ (cid:48)(cid:48) (1 /
3) = − π . The number of samples is M = 10 and the numberof particles is N p = 10 , unless otherwise stated. Drawing a sample from the density (87) was done by by thetransformation method [36], by which the cumulative distribution function of ρ ( x ) is inverted; this process isrepeated µ = 1 . . . N p times for the particles in the sample. The whole process is repeated for a total of M samples. We note that by connecting particle widths to integer number of cells, we do not reach the absolutetheoretical minimum h opt since N g in a domain of fixed size can only be an integer number and ∆ = 1 /N g canonly take discrete values.Before discussing the numerical simulations, let us look qualitatively at the error minimum for the density(87), based on the order of magnitude estimates (50). We assume the kernel width is a multiple of the gridspacing and for simplicity take h = ∆ , i.e., we estimate the density with a boxcar shape. Using the scalingtransformation (8), we vary the grid in the range N g ∈ [16 , , which corresponds to kernel widths in the range h ∈ [0 . , . . We need to compare h with, the density gradient scale length, which for our case ofEq. (87) at x = 1 / gives l = (cid:112) . / π (cid:39) . . The qualitative bias and variance error curves are illustrated inFig. 6, showing that we should expect a minimum in this range of simulation parameters. Notice that the lowestvalue of h corresponds to N g = 48 and is an acceptable grid resolution for our cosine density profile, having grid points per period; it does not, however, correspond to h opt , the error being about above the actualminimum. If we were to set up a grid corresponding to h opt (using the same particle, with h = ∆ ), it wouldcorrespond to N g (cid:39) and we would be somewhat under-resolving the cosine period with only about gridpoints per period. We shall not further discuss this issue but emphasize that grid resolution is a factor that mustbe chosen based on grid truncation errors and independently from the particle shape, whose width is chosen tominimize Q = V + B , where the variance V is a measure of the statistical error. Finally, the quantity N h fromEq. (51) is equivalent to the number of particles per cell, N ppc , in the discretized system. For a total of N p = 10 particles, we find N ppc ∈ [312 , for N g ∈ [16 , , and N ppc (cid:39) at the intersection of the two curves, which We remind the reader that Q is the mean-square error and the actual error is √ Q .
22s close to the minimum of Q = V + B . Note that in the intersection region the bias curve depends morestrongly on h than the variance curve, suggesting that a more efficient way to “locate” (and follow) the minimumof the error curve is by adjusting h rather than by adjusting the number of particles in a simulation. One shouldkeep in mind that the exact numerical parameters discussed above may differ from the exact figures by a factorof a few, coming from the shape coefficients C and C , which were not included in the estimates (50) and (51). h E rr o r ( v a r i a n c e ) V x ) N p h E rr o r ( b i a s ) B ( hl ) Figure 6: Qualitative comparison of the terms representing the bias and variance errors in Eq. (50).Numerical results for the above setup are shown in Fig. 7 at the location x = 1 / . In this simulation weuse the scaling transform, Eq. (8), which is accomplished by changing the grid spacing ∆ (or N g ) and thus theparticle width h = 3∆ = 3 /N g . In the figure each point corresponds to a different N g but we remind the readerthat the grid resolution is irrelevant to BVO calculations and thus we should not expect interference with theresults for h opt and Q min . As predicted by theory, a minimum of the error is achieved at some value h opt . Further,the value of the minimal error changes very little between the four different shapes. The value of h opt is seento vary, with the exception of the trapezoidal shape and the Epanechnikov kernel, which show similar values forboth h opt and Q min . To understand the similarities and differences better, we use the values of the constants C and C from Table 2 to calculate the theoretical Q min and h opt from Eqs. (52) and (53). The results in Table 6confirm our observations and also show good agreement between theoretical and numerical calculations.One important point concerning comparisons between numerical and theoretical results should be made. Thetheoretical calculations involved certain approximations such as termination of the Taylor expansion in Eq. (41)and neglecting terms ∼ /N p in Eq. (44). In contrast, the numerical results are in this sense “exact” since noapproximations are made; the only inaccuracy in the numerical results stems from the finite number of samplesused in the averages. It is plausible that certain density profiles may require keeping some of the other neglectedterms to better describe the optimal width and error. This situation may occur, for example, when at isolatedlocations, x , the second derivative of the density vanishes, ρ (cid:48)(cid:48) ( x ) = 0 . Such occurrences are exceptions ratherthan the general situation and of course, the averaged theory, Eqs. (55), (56), remains unaffected. For thedensity profile (87) we will see below that at the specified above spatial locations, including one extra term inthe theoretical calculations further improves the agreement by a few percent.In Figure 8 we show comparison between local theory, Eqs. (52), (53), averaged theory, Eqs. (55), (56), andsimulations for the two locations, x = 1 / (top panel) and x = 1 / (bottom panel), for the same quadraticspline particle shape (cf. Table 3). The lines labeled “Theory (local)” are plots of formula (49); the lines labeled“Theory (average)” are plots of the integrated (from to ) Eq. (49). Both theory and averaged theory use thevalues of C and C from Table 2. We see that the averaged theory does not change between the top and bottompanels since h opt and Q min do not depend on x . (The apparent difference is due to the slightly different plottingrange.) In order to obtain even better agreement between local theory and simulations, we have included thesmall correction − ρ ( x ) /N p in Eq. (44), which has been neglected so far; that yields an improvement of about . Including that term also explains the slight difference in the theoretical values of h opt and Q min that thediscerning eye would observe in the legend of Fig. 8 (bottom panel) versus the table values in Table 6 (Quadr.spline).Lastly, we present simulations of bias-variance optimization with the fractional width particle shape, Eq. (65).Recall that the fractional particle shape is not a simple scaling transform as the cases considered thus far but is23 shape transform , which additionally changes the measure of its support. The theory developed in Sec. 4 wasfor a scaling transform only, keeping the particle shape unchanged. Therefore that theory cannot be used forquantitative comparison with the following simulations but can nevertheless serve as a guideline to understandingthe numerically observed behavior of Q min and h opt .These simulations were performed at x = 1 / , for two fixed grid sizes, N g = 8 and , with N p = 1000 .The smaller number of particles was chosen to increase (for clarity) the relative numerical error due to the finitenumber of samples: for particles, the value of Q min is expected to be about / ≈ . times larger whilethe value of h opt to be about / ≈ . times larger. As a guideline for the values of the minimal errorand optimal width, we take the averages of Q min and h opt between the two limiting shapes ( ∆ -wide boxcar and -wide linear tent); these average values are Q min ∼ . and h opt ∼ . . The simulation results in Fig. 9show minimum error Q min ≈ . and optimal width h opt ≈ . , consistent with the above predictions. Theobserved discontinuity in the Q ( h ) curve is also expected and easy to explain. For the range . < h < . we use ∆ = 0 . (or N g = 16 ), and for . < h < . , we use ∆ = 0 . (or N g = 8 ). The fractional widthparticle becomes a ∆ -wide boxcar shape for δ = 0 and a -wide linear shape for δ = ∆ . Therefore, on the leftside of the discontinuity the density is estimated by a linear particle shape of width . while on the right itis estimated by a boxcar shape of the same width. Estimating the density by two different shapes (of the samewidth) is expected to produce different errors Q because of the different values of the shape coefficients C and C (see Table 2).From the results in Fig. 9 we conclude that the fractional width particle shape can indeed be used to attainthe BVO minimum error without having to change the type of particle (or charge deposition rule). Particle width, h Sq u a r e d e rr o r , Q = V + B BoxcarQuadratic splineTrapezoidalEpanechnikov
Figure 7: Local bias-variance optimization comparison of four different shapes at x = 1 / , from numericalcomputations. Although the value of h opt changes noticeably between the different shapes, the value of Q min depends little on the shape of the particle, in agreement with the results summarized in Table 2. Shape Q min h opt theoretical numerical theoretical numericalBoxcar . . . . Quadr. spline . . .
139 0 . Trapezoidal . . .
107 0 . Epanechnikov . . .
103 0 . Table 6: Comparison of theoretical values, Eqs. (52) and (53), and numerical values of the optimal particle widthand the error minimum. The more accurate value Q min = 0 . of the minimal error for Epanechnikov kernelis indeed slightly lower than the value of the trapezoidal shape; the table entry has been rounded off to threesignificant figures. 24 .06 0.08 0.10 0.12 0.14 0.16 0.18 0.20 Particle width, h Sq u a r e d e rr o r , Q = V + B Theory (local), h opt = 0.160, Q min = 0.000895Theory (average), h opt = 0.147, Q min = 0.00130Quadratic spline, h opt = 0.167, Q min = 0.000910 Particle width, h Sq u a r e d e rr o r , Q = V + B Theory (local), h opt = 0.138, Q min = 0.00208Theory (average), h opt = 0.147, Q min = 0.00130Quadratic spline, h opt = 0.136, Q min = 0.00198 Figure 8: Bias-variance optimization comparison between exact local theory, averaged theory, and simulationsfor the quadratic spline particle shape. Top panel: x = 1 / . Bottom panel: x = 1 / . + Q m i n = B + V N g = 8 N g = 16 Figure 9: Local bias-variance optimization at x = 1 / with the fractional width particle (65). The particle width h = ∆ + δ varies by changing δ , but at the step discontinuity you have to change both δ and ∆ . (See the text.)25 Summary and conclusions
We have presented analyses of the noise in particle methods used to study electrostatic models in one dimension.We have described kernel density estimation for continuous x , expressing the kernel in terms of a fundamentalkernel K f of width , namely K ( x ) = (1 /h ) K f ( x/h ) . In this form its shape ( K f ) and its width h are representedseparately. Restricting our attention to uniform true electron density ρ ( x ) for these initial studies (and immobileions of uniform, fixed density ρ i ( x ) throughout the paper), we have computed the covariance matrix C ( x, y ) of thenoise in the estimated electron density ρ e ( x ) . There are positive off-diagonal elements of the covariance matrixrelated to the width h of the kernel. But, more importantly, there are constant negative elements, on and off thediagonal (the latter related to negative correlations). These negative matrix elements arise from the fact that thetotal number of particles is fixed at each time step. That is, for example, if the particles are concentrated in onearea (higher density estimate), they will be necessarily more sparse (lower density estimate) in other areas. Thesenegative correlations lead to the property ´ C ( x, y ) dy = 0 , i.e. C ( x, y ) has a zero eigenvalue, ´ C ( x, y ) u ( y ) dy = 0 for eigenfunction u ( y ) = 1 . We compute the estimated electric field from the estimated density by Gauss’s law, ∂E/∂x = ρ i − ρ e , using E (0) = E (1) by charge neutrality, E (0) = 0 by periodicity. We also assume that theapplied potential across the system ´ E ( x ) dx is zero. These boundary conditions and the negative correlationsin C ( x, y ) lead to properties of the noise in the electric field related to a process called the Ornstein-Uhlenbeckbridge , described in Appendix A. The covariance matrix of the electric field C E ( x, y ) is significantly reducedrelative to that of the commonly known Brownian process or the related Brownian bridge, improving the fidelityof simulation results. Because of the assumed periodic boundary conditions on [0 , , the covariance matrices C and C E have translational invariance properties, i.e. C ( x, y ) = C ( x − y ) and C E ( x, y ) = C E ( x − y ) . The latteralso has an eigenfunction, namely u ( x ) = 1 , with zero eigenvalue.We have also investigated cases with non-constant density ρ ( x ) , but still with continuous x . We haveconsidered the total error in the estimated density and analyzed it in terms of bias-variance optimization (BVO.)Small kernel widths have too few particles within their support leading to too much variance; for kernels withwidths that are large compared to a characteristic density gradient scale length, the actual density is smoothedexcessively. The optimum between these two limits is found by BVO. The analysis also shows that this optimumis weakly dependent on the kernel shape for kernels of equal widths. We find that the scaling of the minimalerror Q min with the total number of particles N p is modified to Q min ∼ N − / p compared to the well knownvariance scaling V ∼ N − p .We have analyzed these properties for a grid of discretized x values. In this case the charge depositionrule is expressed in terms of a particle shape . We have discussed an important property to be preserved in thediscrete system: the exact preservation of the net electron charge ´ ρ ( x ) dx ; the discrete version of this propertyis necessary to ensure that the discretized electric field obeys the periodic boundary conditions. If we assumethat the particle shape obeys a sum rule (cid:80) i ∆ S ( x i − ξ ) = ´ S ( x − ξ ) dx = 1 , saying that the discretized integralover the shape equals the exact integral, then exact preservation of charge holds. The particle obeys such a sumrule, if it is the convolution of two kernels, S ( x ) = ´ K ( y ) ˆ K ( x − y ) dy , where one of the kernels satisfies thesum rule. An example of this occurs when one of the kernels is a particle shape, ˆ K = S ( x − y ) , or when finiteelements are used so that ˆ K ( x − y ) = Ψ( x − y ) with Ψ( x − y ) having width equal to an integral multiple of thegrid spacing, ∆ , and satisfying the sum rule. The convolution formula appears naturally in variational particlemethods [4, 12]. It is important to note that the sum rule property for the particle shape does not require thatthe kernel width is an exact multiple of the grid spacing ∆ .We have relaxed the approximations of the analytic calculations of BVO optimization, doing numericalcomputations of the total error as a function of the particle width. The results show good agreement with theanalytic results over a range of particle shapes.As practical applications of the results in this work, we have provided evidence that noise correlations canbe reduced by using sufficiently wide particle shapes, decreasing finite number of particle numerical effects. Innon-uniform density, guidelines for the design, construction, and implementation of computationally efficientparticle shapes that take advantage of the BVO is proposed. In particular, for large values of h opt and small gridspacing ∆ we recommend custom designed particles of low (polynomial) order but sufficiently wide extent as acomputationally efficient alternative to the traditional spline shape functions. When h opt ∼ ∆ , we recommendfractional width particle shapes, which can also be used to follow the optimal particle width in the course of asimulation while keeping the charge deposition rule unchanged and providing computational efficiency.The bias-variance trade-off was discussed and shown to be important in the context of PIC simulations inplasmas in Ref. [37]. In the present paper we stress the analytical development of this idea as applied to particle-based numerical methods [38–40] and present detailed analysis of particle shapes related to the width of thebias-variance minimum. We also discuss exact charge conservation and the negative correlations due to a fixed26umber of particles and their influence on the statistical properties of the electric field. Acknowledgments
Sandia National Laboratories is a multimission laboratory managed and operated by National Technology &Engineering Solutions of Sandia, LLC, a wholly owned subsidiary of Honeywell International Inc., for the U.S.Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paperreviews objective technical results and analysis. Any subjective views or opinions that might be expressed donot necessarily represent the views of the U.S. Department of Energy or the U.S. Government.The work of EGE was supported in part by NASA WV EPSCoR Grant
A Appendix: Brownian bridge and Ornstein-Uhlenbeck bridge
Consider a random walk on < t < T (continuous time), with dv w dt = r w ( t ) , (88)where v w is the velocity of a Brownian particle and r w is its random acceleration, with (cid:104) r w ( t ) r w ( s ) (cid:105) = V τ c δ ( t − s ) .Here, V is the variance of r w ( t ) and τ c is the correlation time. We start by taking v w (0) = 0 , and find v w ( t ) = ´ t r w ( τ ) dτ , leading to (cid:104) v w ( t ) v w ( s ) (cid:105) = ˆ t dτ ˆ s dσ (cid:104) r w ( τ ) r w ( σ ) (cid:105) = V τ c min ( t, s ) , (89)the standard result [30]. In the analogy with the results of Sec. 3.2, time takes the place of the distance x , theacceleration takes the place of the density, the velocity takes the place of the electric field, and the displacementtakes the place of the electrostatic potential.Now for each realization of the noise r w ( t ) , consider the modified process v b ( t ) = v w ( t ) − tv w ( T ) T . (90)At this stage we still have v b (0) = 0 . The usual random walk has (cid:104) v w ( T ) (cid:105) = 0 but this modified process has v b ( T ) = 0 for each realization of the noise. This is the analog of the condition E (1) = E (0) of Sec. 3.2. We find dv b dt = r ( t ) = r w ( t ) − v w ( T ) T (91)and v b (0) = 0 , with the result C ( t, s ) ≡ (cid:104) r ( t ) r ( s ) (cid:105) = V τ c (cid:20) δ ( t − s ) − T (cid:21) . (92)This is proportional to the covariance in Eq. (22). Notice the stationarity condition C ( t, s ) = C ( t − s ) , theanalog of the translation invariance condition in the spatial context, and ´ T C ( t − s ) ds = 0 , as in Eq. 35. Wealso find the covariance matrix for the the Brownian bridge velocity v b ( t ) , (cid:104) v b ( t ) v b ( s ) (cid:105) = ˆ t dτ ˆ s C ( τ, σ ) dσ = V τ c [ min ( t, s ) − ts ] . (93)It is clear that the process of subtracting tv w ( T ) /T in Eq. (90) is equivalent to integration of dv b /dt = r ( t ) withthe covariance matrix r ( t ) given in Eq. ((92)).The final step is to consider the Ornstein-Uhlenbeck bridge [27] by starting with the system dv/dt = r ( t ) asin Eq. (91), for now just relaxing the requirement v (0) = 0 . We find v ( t ) = v + v ( t ) = v + ´ t r ( τ ) dτ , leadingto (cid:104) v ( t ) v ( s ) (cid:105) = C + C ( t ) + C ( s ) + C ( t, s ) , (94)27here C = (cid:104) v (cid:105) , C ( x ) = (cid:104) v v ( t ) (cid:105) , C ( y ) = (cid:104) v ( s ) v (cid:105) , and C ( t, s ) = ˆ t dτ ˆ s dσ (cid:104) v ( τ ) v ( σ ) (cid:105) , the exact analog of Eqs. (29)-(32). The net displacement x ( t ) of the particle is found by integrating dx/dt = v ( t ) ,with x (0) = 0 , so that x ( t ) = ´ t v ( τ ) dτ . The Ornstein-Uhlenbeck bridge modification is this: for each randomwalk, we choose v so that the net displacement at t = T is zero, ´ T v ( τ ) dτ = 0 . This zero net displacementcondition is the analog of the zero potential difference requirement of Eq. (25)). The covariance matrix (cid:104) v ( t ) v ( s ) (cid:105) is obtained by the methods outlined in Sec. 3.2 and in Appendix B, and illustrated in Fig. 2. Finally, the Ornstein-Uhlenbeck bridge for smooth correlations, (cid:104) r w ( t ) r w ( s ) (cid:105) = V τ c δ ( t − s ) → V τ c K ( t − s ) is treated in Appendix Band in Fig. 2.It is interesting to note that in analogy with Eq. (90), we can relate the Ornstein-Uhlenbeck displacementvariable x ( t ) to the displacement for the Brownian bridge variable x b ( t ) having v = 0 by x ( t ) = x b ( t ) − tx b ( T ) T .
The Brownian Bridge defined above has the requirement that the particle velocity return to zero at t = T .The Ornstein-Uhlenbeck bridge has the further requirement that the particle displacement return to its originalposition, i.e. that the average velocity ´ T v ( t ) dt/T be zero. B Appendix: Electric field covariance matrix for general kernels
In this appendix we derive the covariance matrix for the electric field from Eq. (28), relaxing the special case ofEq. (22) to a general kernel, C ( x, y ) = K ( x − y ) − , (95)where the /N p factor has been suppressed and K ( x ) = (1 /h ) K f ( x/h ) . We take K ( x ) → K ( x ) to be extendedto be periodic of period , so that it is even about x = 1 / . As in Eq. (28), we conclude C E ( x, y ) = C + C ( x ) + C ( y ) + C ( x, y ) . (96)For simplicity we pick the fundamental kernel K f ( x ) to be the boxcar. We find XC E = ( ∂ x + ∂ y ) C E = 0 ,which implies translational invariance C E ( x, y ) = C E ( x − y ) , and symmetry C E ( x, y ) = C E ( y, x ) implies C E ( x, y ) = C E ( | x − y | ) . Finally, these relations imply that C E ( x, y ) is periodic in x − y with period .An alternate approach begins with C ( x, y ) = (cid:104) ˜ ρ ( x )˜ ρ ( y ) (cid:105) = (cid:104) ˜ E (cid:48) ( x ) ˜ E (cid:48) ( y ) (cid:105) . This leads to C ( x, y ) = ∂ ∂x∂y C E ( x, y ) = ∂ ∂x∂y C E ( x − y ) , (97)with the last step following from translational invariance. Eq. (97) leads to: C ( x − y ) = − ∂ x C E ( x − y ) . (98)For C ( x ) = δ ( x ) − (ignoring /N ppc factor), we find ∂ x C E ( x ) = − C ( x ) = − δ ( x ) + 1 implies C E ( x ) = D − | x | + x . (99)This satisfies C E (1) = C E (0) = D or C E ( − /
2) = C E (1 / , related to the boundary condition on the electricfield from overall charge neutrality. The condition ´ C E ( x ) dx = 0 , from the zero applied potential condition ´ E ( x ) dx = 0 , leads to D = 1 / . These results are in agreement with Eqs. (33) and (34) when the factor /N p is reinstated.Finally, for the linear kernel, C ( x ) = K fL ( x/h ) /h − , we solve ∂ x C E ( x ) = − K fL ( x/h ) /h + 1 . The neutralityrequirement C E ( − /
2) = C E (1 / is easily seen to be satisfied. These results as well as those from h = 0 (Eqs. (99), (33), (34)) are plotted in Fig. 2. Note that the cusp at x = y for h = 0 is smoothed and relation ´ C E ( x, y ) dy = 0 is found to hold. With the /N p factor, these results agree with the results derived by themethod above (see Eq. (96)). 28 Appendix: Scaling of the kernel
The information specific to a given kernel is contained in the shape coefficients C and C [cf. Eq. (48)].Sometimes it may be more convenient to work with the scaled kernel (8) instead of the fundamental kernel.Additionally, published literature may define the fundamental kernel with a width different from unity. Tomake a connection between scaled or differently defined kernels and the fundamental kernel as defined in ourpresentation, we examine how the coefficients C and C scale under a scaling transformation h → αh forarbitrary h and scaling factor α > . For example, to obtain the linear particle shape in Table 3 from the linearfundamental kernel in Table 1, we use Eq. (8) with h = α ∆ with α = 2 ; similarly, for the quadratic splineand trapezoidal particles we use α = 3 , etc. (We remind the reader that not all fundamental kernels allow for ascaling transform leading to a particle shape satisfying the sum rule, e.g., the Epanechnikov fundamental kernel.)The kernel scales as K ( x ) = 1 h K f (cid:16) xh (cid:17) → αh K f (cid:16) xαh (cid:17) = 1 α K f (cid:18) ζα (cid:19) , where ζ = x/h . It is easy to verify that the scaled kernel is also normalized to unity [cf. Eq. (3)]. Now wecalculate the scaled coefficients, C s, and C s, . We have C s, = ˆ dζ α K f (cid:18) ζα (cid:19) = 1 α ˆ dγ K f ( γ ) = 1 α C ,C s, = ˆ dζ ζ α K f (cid:18) ζα (cid:19) = α ˆ dγ γ K f ( γ ) = α C , with γ = ζ/α . Again, formulas (48) are based on the fundamental kernel and therefore yield the values of C and C on the right hand sides above, as seen in Table 2. Calculating the scaled values of Q min , h opt , and W Q we get Q s , min ∼ (cid:16) C s, (cid:112) C s, (cid:17) / = (cid:18) C α (cid:112) α C (cid:19) / = (cid:16) C (cid:112) C (cid:17) / ∼ Q min ,h s , opt ∼ W s , Q ∼ (cid:32) C s, C s, (cid:33) / = (cid:18) α C α C (cid:19) / = 1 α (cid:18) C C (cid:19) / ∼ α h opt ∼ α W Q . We see that Q min remains unchanged , while h opt and W Q scale inversely with the scaling factor α . References [1] F. H. Harlow. The particle-in-cell computing method for fluid dynamics.
Methods in Computational Physics ,3:319–343, 1964. 1[2] R. W. Hockney. Computer experiment of anomalous diffusion.
Phys. Fluids , 9(9):1826–1835, 1966. 1[3] A. B. Langdon and C. K. Birdsall. Theory of plasma simulation using finite-size particles.
Phys. Fluids ,13(8):2115–2122, 1970. 1[4] H. R. Lewis. Energy-conserving numerical approximations for Vlasov plasmas.
J. Comput. Phys. , 6(1):136–141, 1970. 1, 5.1, 7[5] J. U. Brackbill and H. M. Ruppel. Flip: A method for adaptively zoned, particle-in-cell calculations of fluidflows in two dimensions.
Journal of Computational Physics , 65(2):314–343, 1986. 1[6] J. M. Dawson. Particle simulation of plasmas.
Rev. Mod. Phys. , 55(2):403–447, April 1983. Publisher:American Physical Society. 1[7] R. W. Hockney and J. W. Eastwood.
Computer Simulation Using Particles . Taylor & Francis Group, NewYork, 1988. 1, 2[8] C. K. Birdsall and A. B. Langdon.
Plasma Physics via Computer Simulation . CRC Press, New York u.a.Taylor & Francis, 1 edition edition, October 2004. 1, 2, 2, 5, 5.1, 5.1[9] G. Chen, L. Chacón, and D. C. Barnes. An energy- and charge-conserving, implicit, electrostatic particle-in-cell algorithm.
J. of Comput. Phys. , 230(18):7018 – 7036, 2011. 12910] S. Markidis and G. Lapenta. The energy conserving particle-in-cell method.
J. Comp. Phys. , 230(18):7037– 7052, 2011. 1[11] J. Squire, H. Qin, and W. M. Tang. Geometric integration of the Vlasov-Maxwell system with a variationalparticle-in-cell scheme.
Phys. Plasmas , 19(8):084501, August 2012. 1[12] E. G. Evstatiev and B. A. Shadwick. Variational formulation of particle algorithms for kinetic plasmasimulations.
J. Comput. Phys. , 245:376–398, July 2013. 1, 2, 5.1, 5.1, 7[13] B. A. Shadwick, A. B. Stamm, and E. G. Evstatiev. Variational formulation of macro-particle plasmasimulation algorithms.
Phys. Plasmas , 21(5):055708, May 2014. 1[14] A. B. Stamm, B. A. Shadwick, and E. G. Evstatiev. Variational Formulation of Macroparticle Models forElectromagnetic Plasma Simulations.
IEEE Trans. Plasma Sci. , 42(6):1747–1758, June 2014. 1[15] P. J. Morrison. The Maxwell-Vlasov equations as a continuous hamiltonian system.
Physics Letters A ,80A(5):383–386, 1980. 1[16] A. Weinstein and P. J. Morrison. Comments on: The Maxwell–Vlasov equations as a continuous Hamiltoniansystem.
Phys. Lett. , 80A(4):235–236, 1981. 1[17] M. Kraus, K. Kormann, P. J. Morrison, and E. Sonnendrücker. GEMPIC: geometric electromagneticparticle-in-cell methods.
J. Plasma Phys. , 83(4), August 2017. Publisher: Cambridge University Press.1[18] C. K. Birdsall and D. Fuss. Clouds-in-clouds, clouds-in-cells physics for many-body plasma simulation.
J.Comp. Phys. , 3(4):494–511, April 1969. 1, 5.1[19] A. B. Langdon. Kinetic theory for fluctuations and noise in computer simulation of plasma.
Phys. Fluids ,22(1):163–171, January 1979. 1[20] M. Kotschenreuther. Bull. Amer. Phys. Soc., 1988. 1[21] R. E. Denton and M. M. Kotschenreuther. Delta-f algorithm.
Journal of Computational Physics , 119:283–294, 1995. 1[22] S. E. Parker and W. W. Lee. A fully nonlinear characteristic method for gyrokinetic simulation.
Physics ofFluids B: Plasma Physics , 5(1):77–86, January 1993. Publisher: A. Institute Phys. 1[23] G. Hu and J. A. Krommes. Generalized weighting scheme for delta-f particle simulation method.
Phys.Plasmas , 1(4):863–874, April 1994. Publisher: A. Institute Phys. 1[24] S. Brunner, E. Valeo, and J. A. Krommes. Collisional delta-f scheme with evolving background for transporttime scale simulations.
Phys. Plasmas , 6(12):4504–4521, November 1999. Publisher: A. Institute Phys. 1[25] C. M. Bishop.
Pattern Recognition and Machine Learning . Information Science and Statistics. Springer-Verlag, New York, 2006. 1, 2, 4.1[26] C. Sammut and G. I. Webb, editors.
Encyclopedia of Machine Learning and Data Mining . Springer US, 2edition, 2017. 1, 4.1[27] Y. Chen and T. Georgiou. Stochastic bridges of linear systems.
IEEE Transactions on Automatic Control ,2016. 1, 3.2, A[28] S. Corlay. Properties of the ornstein-uhlenbeck bridge. arXiv preprint arXiv:1310.5617 , 2013. 1[29] A. Mazzoloa. Constraint ornstein-uhlenbeck bridges.
J. Mathematical Physics , 2017. 1[30] D. Revuz and M. Yor.
Continuous Martingales and Brownian Motion . Springer, Berlin ; New York, 3rdedition, December 1999. 1, A[31] V. A. Epanechnikov. Non-Parametric Estimation of a Multivariate Probability Density.
Theory Probab.Appl. , 14(1):153–158, January 1969. Publisher: Society for Industrial and Applied Mathematics. 2, 4.1[32] R. Mansuy and M. Yor.
Aspects of Brownian Motion . Springer, Berlin, Heidelberg, 2008. 3.23033] Q. Li and J. S. Racine.
Nonparametric Econometrics . Princeton University Press, December 2006. 4.1[34] E. B. Becker, G. F. Carey, and J. T. Oden.
Finite Elements: An Introduction , volume 1. Prentice-Hall,Inc., Englewood Cliffs, N.J., 1981. 5.1[35] S. Ross.
A First Course in Probability, p. 364 . Prentice Hall, Upper Saddle River, NJ, 7 edition edition,2006. 5.2[36] W. Press, S. Teukolsky, W. Vetterling, and B. Flannery.
Numerical Recipes in Fortran . Cambridge UniversityPress, 2 edition, 1993. 6.2[37] W. Wu and H. Qin. Reducing noise for PIC simulations using kernel density estimation algorithm.