{'fraction_non_alphanumeric': 0.033738148841043855, 'fraction_numerical': 0.009290674445888956, 'mean_word_length': 4.119805194805195, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A new method to dispatch split particles in Particle-In-Cell codes

Particle-In-Cell codes are widely used for plasma physics simulations. It is often the case that particles within a computational cell need to be split to improve the statistics or, in the case of non-uniform meshes, to avoid the development of fictitious self-forces. Existing particle splitting methods are largely empirical and their accuracy in preserving the distribution function has not been evaluated in a quantitative way. Here we present a new method specifically designed for codes using adaptive mesh refinement. Although we point out that an exact, distribution function preserving method does exist, it requires a large number of split particles and its practical use is limited. We derive instead a method that minimizes the cost function representing the distance between the assignment function of the original particle and that of the sum of split particles. Depending on the interpolation degree and the dimension of the problem, we provide tabulated results for the weight and position of the split particles. This strategy represents no overhead in computing time and for a large enough number of split-particles it asymptotically tends to the exact solution.


Introduction
The multi-scale character of laboratory and astrophysical plasmas is ubiquitous. Examples include, but are not limited to, the turbulence in the solar wind or in fusion experiments [1], the reconnection of magnetic field lines [2], the formation and propagation of shocks [3]. In all these systems, spatial and temporal scales can span many orders of magnitude and impose stringent constraints on grid-based numerical codes. In situations where small scales are spatially localized, a solution to circumvent some of the constraints is to adaptively refine the mesh in the regions of interest. These Adaptive Mesh Refinement (AMR) technique are widely used in fluid codes (see e.g. [4], [5]) and have also begun to be successfully implemented in Particle-In-Cell (PIC) codes ( [6], [7], [8], [9], [10], [11], [12]).
In PIC codes, a large collection of physical particles is described by a smaller set of computational particles or "macro-particles" [13]. Such a model is acceptable as long as the statistical properties of the set of the macroparticles is close to that of the physical particles. It means that the moments resulting from integration in velocity space of the particle distribution and of the macro-particle distribution are very close, whatever their order. For the sake of readability, we shall use in the rest of the paper the word "particle" instead of "macro-particle". We remind that the moment of order n of a distribution function is the integral over velocity space of this distribution function multiplied by the velocity at power n. Density, bulk velocity and pressure are the moments of order 0, 1 and 2 respectively, and are the most widely used.
A characteristic of the particles used in PIC codes is their finite size, whose spatial profile is given by a continuous function called the "assignment function" or "shape factor". This function has bounded support and it is thus zero outside of it. Importantly the size of the particle (their support) depends on the local spatial resolution as well as on the order of interpolation (see Ref. [13] for an extensive review). In the context of AMR codes this is crucial: when a particle moves from a position where the grid size is ∆ to enter a refined region where the grid size is ∆/r, the size of the particle should also decrease by a factor r, where r is the refinement factor. When entering this region of smaller grid size, the "parent" particle has to be split into two or more "children" particles for at least two reasons: the first is to avoid spurious self-forces that are associated with the time derivative of the assignment function [14,15], the second is to maintain a sufficiently large number of particles per cell [16,17,11], [18] in order to insure stability and accuracy. In that respect, particle splitting is also useful in codes using a uniform mesh. Therefore an important question that needs to be addressed is: how many children particles should the parent particle be split into? And as a corollary question, once this number is fixed, where should the children particles be dispatched, with which velocities and with which weights? So far, only empirical answers to these questions can be found in the literature.
Instead, this paper aims at addressing these questions rigorously by using an optimization technique that minimizes the difference between the assignment function of the parent particle and the sum of the assignment functions of the children particles.
The paper is structured as follows. In section 2, we review the important features of basis-splines (B-splines), as these are used to represent macroparticles in the vast majority of PIC codes. In section 3, we discuss the dichotomy between the exact solution of the splitting problem and the approximate one we are proposing. In section 4, we present some general considerations on splitting requirements, and the possible splitting patterns depending on the dimension of the problem. In section 5, we discuss the way to evaluate the accuracy of the method, and provide the optimal parameters of the split particles. In section 6, we discuss the consequences of the splitting on the level of density fluctuation on the refined grid. In section 7, we discuss possible implications of our conclusions for the inverse problem of merging process.

B-splines for Particle-In-Cell method
The particles' assignment functions used in PIC codes are generally Bspline functions [19]. These functions can be of any degree: the higher the degree, the larger the particle size, i.e. the length of the support of the assignment function and the smoother the particle shape. Examples of Bsplines of degree one to three are shown in Fig. 1. For clarity, we will treat in this paper the linear (p = 1), quadratic (p = 2) and cubic (p = 3) cases.
The p = 0 case has been proven to be unstable [13], while values of p > 3 are very diffusive, both cases are rarely used in actual simulations. A B-spline of degree p is the union of p + 1 polynomials of degree p on p + 2 knots.
To handle B-splines in a more general way, we shall use the notation S p (x) for a B-spline of degree p, defined on a grid of mesh size ∆ centered at x = 0 and s p (x − δ) the B-spline of degree p defined on a grid of mesh size ∆/r centered on x = δ. These represent respectively the parent and children particles. By construction, B-splines have their integral equal to unity. As an illustration, on the right panel of Fig. 1, the B-spline s 3 (x) for r = 2 is depicted in dark gray, with its associated knots as white bullets. When a particle leaves a position where the grid size is ∆ to reach a new one where and right panel is the B-spline for p = 3. In each panels, light gray is used to represent the B-splines with a refinement factor r = 1 (i.e. the parent particle). In addition, an example of a B-spline with a refinement factor r = 2 (i.e. the child particle), is shown on right panel in dark gray. Black knots are for r = 1 and white knots are for r = 2.
the grid size is ∆/r, the support of the B-splines has to be divided by r. For convenience, the coarse and fine grid have been displayed in the right panel of Fig. 1, It is clear that when moving in a domain of smaller mesh size, the spatial profile of the assignment function of a particle can not be preserved if this particle is substituted by only one with a narrower support (see right panel of Fig. 1). The obvious consequence is that a parent particle needs to be split into several child particles. One may hope that the error introduced when splitting the parent particle somehow becomes unimportant when splitting many parent particles. However, we shall demonstrate that this is not the case and that the number and placement of child particles may have dramatic consequences on the conservation of the associated distribution function. For convenience we restate here the analytic form of B-splines with ∆ = 1. The linear B-spline is given by the quadratic B-spline is defined by and the cubic B-spline is defined by Without loss of generality, we shall focus in this paper on the special case (for the numerical values) where ∆ = 1 for the parent and r = 2 for the children, but we shall keep ∆ and r in the notations. A particle k is thus defined by its weight w k , position x k and velocity v k , so the continuous density function at position x where f (x, v, t) is the distribution function, is approximated by the sum In the 2D and 3D cases, the assignment function is the product of the 1D assignment function defined for each directions, for example:

Exact and approximate solution
In a pioneering study by Lapenta [19], it was emphasized that the rezoning of particles (splitting or merging) can not be exact for uniform meshes. While this assertion is true for a uniform grid, it no longer applies if the children are dispatched to a mesh of different grid size from that of their parent. In this case the size of the support of the assignment function of the children is smaller (by a factor r) than that of their parent. Indeed, an exact solution exists when r is an integer [20]. However when applied to PIC codes, the solution comes at the cost of a large number of split particles: given the dimension d of the grid and a refinement factor r = 2, the number of children In section 5, we provide the parameters for the exact splitting in one dimension, for B-splines of degree p = 1, p = 2 and p = 3. More specifically, we provide the w i and δ i satisfying the relation The 2D and 3D generalization is straightforward using Eq. (6). To our knowledge, in AMR codes using macro-particles (e.g. [7], [8], [10]) the splitting of particles is neither exact nor optimized. In addition, these methods are all based on split particles having the same velocity as their parent. While arbitrary, this choice simplifies the problem and it is the one Each panel is for d = 1.
we adopt here. We now discuss in more details the splitting methods existing in the literature.
• In the study of Fujimoto [7], the symmetry in the 2D case is preserved with four children. The children are located at a distance δ ⋆ = ∆/ √ N P P C from the parent where ∆ is the grid size associated to the parent and N P P C is the number of particle in the cell where the parent is located. With ∆ = 1 and N P P C ∼ 100 (which is a commonly used value), one obtains δ ⋆ ∼ 0.1. of Ref. [10] shows that, in our notation, these children are uniformly dispatched along each direction with a spacing 1 / 2 (∆/r) = 1 / 4 for ∆ = 1 and r = 2. Hence, in the 2D case, δ ⋆ = 0.25 for four children (keeping in mind that a larger, even number of split particles can be used with this same δ ⋆ value).
• In the study of Muller et al. [8] the symmetry is not preserved. The parent is split into two children independently of the dimension of the problem. Furthermore, the position δ ⋆ of the two children is not precisely given, except that the independent shift of each of them, relative to the parent, is "small". We shall consider that δ ⋆ is of the order of 0.1.
The star notation for δ ⋆ is intended to outline the fact that these values are defined on the fine grid (see right panel of Fig. 1) and not on the coarse one.
Making a clear difference with the choices of these studies, we stress that it is crucial to ensure that the assignment function of the parent and of the set of children are as close as possible, at every points of their support.
Hence, the optimization process must rely on a procedure considering these assignment functions as continuous functions, and not solelly considering a finite number of their moments (density, momentum, energy) as in the studies cited above.
To make clear the close relation between the position of the children relative to their parent and the associated changes in the density profile, we randomly distribute 100 parents in a single, 2D cell of unitary mesh size (some more particles are also deposited on the adjacent cells to avoid the drop of density on the border of the cell). We Focused on the 2D case as this case is nowadays more prevalent than the 1D case. Nonetheless, to ease the representation, we decided to display a line out along the X direction, obtained at a given Y position. As an illustration, we take an assignment function of degree p = 1. The density profile n(x), for a given y, calculated The first two panels correspond to the two values of δ ⋆ used in previous studies [7], [8], [10]. The last panel shows the method presented in this paper.
It is clear from Fig. 3 that the choice of δ ⋆ is crucial in order for the density distribution of the children to be as close as possible to that of the parents.
In addition, the density profiles show larger total variation depending on the choice of δ ⋆ , which may have consequences on the density level of fluctuations.
This point will be deepened in section 5.

Constraints on the splitting strategy
The problem treated in this paper is to determine for a single parent of weight unity located at the origin, the number of children N, their weight . . N} so that the associated assignment functions are as close as possible. If the dimension d of the problem is larger than one, then δ ⋆ i are vectors. It is reasonable to think that the larger N, the smaller the associated error. The yet unspecified definition of this error will be discussed in section 5, but it essentially quantifies how far the assignment function of the children are from the assignment function of the parent. In order to control the accuracy we want to reach, we keep N as a free parameter.
We focus on the position of the children because, as already mentioned, their velocity is the one of their parent.
As already said, δ and δ ⋆ are the position of the children, relative to their parent, defined on the coarse (of grid size ∆) and fine (of grid step ∆/r) grid, respectively. Hence, δ ⋆ = rδ. Yet, we also use δ, defined on the coarse grid, in some cases, as in Eq. 7. From a technical point of view, refined particles are dispatched on the refined grid, so their positions have to be defined on this refined grid. The method for r = 2 is also applicable, but needs to recalculate the w i and δ ⋆ i as discussed in the next section. Such calculations can be done using the dedicated optimization code [21].
Among all the sets of w i and δ ⋆ i , we want to pick the ones for which the assignment function of the children is the closest to the one of the parent. This is illustrated in Fig. 4 for the simple 1D case, with p = 1 and N = 2, where N is the number of split particles on the fine grid for a single particle on the coarse grid. On the left panel, the two children are depicted by a black line, while the parent is in gray. One notes that the support of the assignment function of the children being half the one of the parent, its maximum value is twice the one of the parent because B-splines are normalized functions.
In the left panel of Fig. 4, the two children are deposited with a weight 1 / 2 in order to conserve mass. On the middle panel, we display the sum of the assignment functions of the two children (solid black line), as well as the one of the parent (also in gray). On the right panel, we display the difference between the assignment function of the parent and the one of the children. Red is used when this difference is positive and blue is used when it is negative. The best solution is the one for which the total (red plus blue) surface is as small as possible. We can then define the cost function in order to determine, for each d, p and N values, the w i and δ i for each child i. For the 1D case, the cost function is: For the 2D case, in order to preserve spatial symmetry (which is a neces-sary condition to reach the best solution), there are two ways to split a single parent in four children of equal weight : (i) the four particles are located at (x = ±δ, y = ±δ) or (ii) two are located at (x = 0, y = ±δ) and the two others are located at (x = ±δ, y = 0). This first solution is depicted by purple bullets in middle panel of Fig. 5 and is called pattern of type 2, while the second one is depicted by pink bullets and is called pattern of type 1.
One can notice that pink children are located at a distance δ from its parent while the purple ones are located at δ √ 2. This farther location will have consequences on the δ value as well as on the accuracy that we will discuss in next section. If a single child of weight w 0 (associated to pattern of type

Accuracy of the splitting strategy
For a subset of children, the associated pattern (defining the number and space distribution of children) is called τ j where j is the index of this subset.
Naming N j the number of children for the j subset of children with pattern τ j , the assignment function of this subset of particles is noted In order to make the notations as explicit as possible, index j of T is the index of the subset of children (associated to type τ j ), while index i refers to a given child of this subset. The meaning of the type τ j index can be clearly explained for the 3D case. A child i belonging to the pattern of type τ j is dispatched at a position where a, b and c belong to {-1, 0, +1} andx,ŷ,ẑ are unit vectors. It is so because the δ i differ one from the other because of the different {a, b, c} values, but rely on the same δ j value for the pattern of type τ j . The type of the pattern j is defined as τ j = a 2 + b 2 + c 2 so the associated pattern contains all possible values of {a, b, c} satisfying this relation. In Fig. 5, pink, purple and lime pattern are then associated to τ j = 1, 2 and 3, respectively. The naming of the τ j = 0 pattern is then obvious. The full knowledge of a subset j of children is then given by the triplet (τ j , w j , δ j ). This can be extended straightforwardly to 2D and 3D cases.
Eq. (8) can hence be written in a more compact way as In Eq. (11), M is the number of subsets (each one being associated to a pattern type τ j ) so that the total number of childern in Eq. (8) is given by It is clear that a pattern type τ j for j ∈ [1, M] can have zero or multiple occurences, with different w j and δ j values (see right panel of Fig. 2). Because of the piecewise definition of B-splines, the analytical form of this cost function is not easy to obtain. However, optimal parameters can be obtained numerically and tabulated so to be used in codes at run time. Innocenti et al. [10], which is off the chart for d = 3.
We calculate Q p N defined by Eq. (8) for each combinations of p = 1, 2, 3 and d = 1, 2, 3. The values obtained for the single free parameter δ ⋆ 1 (reminding that the star notation is associated to the fact that this value is defined on the refined grid) are depicted in Fig. 6. The left panel is for d = 1 and two children (N 1 = 2), the middle panel is for d = 2 and four children (N 1 = 4)), and the right panel is for d = 3 and six children (N 1 = 6). The curves correspond to the τ 1 = 1 pattern, indicated by the pink positions in Fig. 5, whatever the value of d. In each panels, the red curve corresponds to p = 1, green is for p = 2 and blue is for p = 3. One notes that it exists for each set of p and d values, an absolute minimum and no local minima, i.e. an optimal δ ⋆ 1 value to dispatch the split particles. In Fig. 6, we also displayed, using bullets, the positions of the minima for each curves.
The minimum values Q p N are increasing with d, but decreasing with p. This last point results from the fact that the larger p, the smaller the total variation of the function. Moreover, the optimum δ ⋆ 1 value is increasing with d, whatever the p value. This results from the increasing values of N 1 with d, hence the decrease of the weight of each child: being lighter, they need to be dispatched further from the parent to fulfill the tail of the assignment function of the parent. Furthermore, the support of a B-spline is increasing with p. The support of a child being half the one of its parent, this child has to be dispatched farther from its parent. As a consequence, δ ⋆ 1 increases with p.
As already discussed in section 3, the suggested δ ⋆ 1 value in Ref. [7] and [8] are close to 0.1. This value is off chart, very far from the optimal δ ⋆ 1 value and it is associated to a very large error. In Ref. [10], the suggested value of . It clearly appears that these two values are not the best one in order to satisfy our accuracy criterion. Moreover, as can be observed on Fig. 3, these δ ⋆ 1 values also increase the level of fluctuations of the density (and higher order moments) profile, which generally compromises the stability of the code. The numerical increase of the level of fluctuation will be provided and discussed in the next section.
With the power 2 in Eq. (8), one recognizes the use of the L 2 norm. We tried to use a different power, namely 1 and 3. The derivative of Q p N is smaller for higher power values in the norm, meaning that the associated minimum value is approximately at the same position, but in a shallower potential well. The good stability for the location of the minima then suggests that the results given here do not depend on the choice of the norm.
We now carry out the same optimization with two degrees of freedom,  respectively, and columns are for p = 1, 2, 3 from left to right, respectively. We use a different scaling for the Y -axis, depending on the dimension d, but the same at a given d value whatever the p value. The reason is that the relation w 0 = 1 − N w 1 has to be satisfied as well as the positiveness of the w i . As the number of children N increases with d, the range of possible values for w 1 thus decreases. The scaling of the X-axis is the same for all the nine panels, but the origin of the X-axis depends on the degree p of the B-splines: the larger p, the larger the origin of the X-axis. As for Fig. 6, this is a consequence of the increasing support of B-splines with their degree p.
feature from these panels is that with two degrees of freedom, the property of having a single absolute minimum is preserved, whatever the p and d  values. We can conjecture it is so for larger M value, even if this point is not mandatory for the numerical resolution of this problem.
We pursue these calculations to find the w j 's and δ ⋆ j 's values that minimize Q p N for larger values of N using Eq. (11). We report these w j values in Tab. 1, 3 and 5 for d = 1, d = 2 and d = 3, respectively. The associated δ ⋆ j are in Tab. 2, 4 and 6 for d = 1, d = 2 and d = 3, respectively. Note that in Tab. 1 to 6, the first row contains the total number N of children for this configuration, while the second contains the type τ j of the pattern used for this subset j (one per line) of children that have to be considered in the summation of Eq. (11). For the interested reader, the associated python code for this optimization [21] (using the NLopt optimization library [22]) can be easily used and extended for larger p, different r and/or larger N. It is also very important (for the implementation on an AMR code) to remind that the δ ⋆ j values are given on the refined mesh in units of ∆/r = 1 / 2 .  The most important for PIC simulations is to decrease as much as possible the value of Q p N . In Fig. 8 we show the optimal Q p N as a function of the number of children N.     The parameters for exact splitting in two and three dimensions are not tabulated (except for d = 2 and p = 1) because they involve some more complex patterns, that would need a different notation. But these values can be deduced very straightforwardly from the ones at one dimension from Eq. (6).

Density fluctuation level when splitting many particles
While the values of Q p N are important to evaluate the accuracy of the method, it is not that easy to interpret the associated consequences in a PIC code. The level of density fluctuations is the root mean square of the density over a large enough number of grid points (for an acceptable statistics) in an homogeneous system. It has a strong impact on the stability of non-AMR PIC codes [13,16] but in the AMR PIC context, a too large increase of the density fluctuation level can dramatically compromised the stability of such simulation [17]. We hence evaluate how the density fluctuation level is modified by splitting a set of parent particles using the optimal splitting described above.  Tab. 1 to 6, we split each parent particle, and calculate the associated density value at the same location, i.e. at the center of each cell. As for the parent, we then calculate over this sample of 4000 points the standard deviation σ ⋆ 100 .
Finally, Fig. 9 depicts the relative increase (in percentage) of the density fluctuation level.
As a first remark, the level of fluctuation of the childern is always larger than the one of the parents. Without clear evidence, one could suspect that the larger total variation of the assignment function of the children compared to the one of the parent could play a significant role in this feature. Furthermore, Fig. 9 clearly shows that when increasing the number of children  In Tab

Conclusions
We presented a new method to find the optimum way to dispatch particles in PIC codes when a splitting strategy is needed in a AMR context. Splitting techniques have already been investigated and discussed, but not necessarily in the framework of AMR. In this specific case, the change of the grid mesh for the split particles has strong consequences. Contrary to most of the previous studies, we emphasized the fact that the way to dispatch split particles must depend on the degree p of the B-spline used as assignment function for the particles. We also demonstrated that such splitting is optimum for a symmetrical distribution of the split particles, meaning that the number of split particles is constrained by this requirement. This symmetry is preserved by adding a particle in the set of split particles at the parent's particle locus. We then provide in this paper the loci of the split particles in Fig. 5, and the associated weights and positions of these split particles in Tab. 1 to 6, for first, second and third order of interpolation, and for both 1D, 2D and 3D cases. The Python code [21] we use for this optimization can then be run for refinement factor larger than two, not tabulated in this paper.
For the 1D cases, the exact solution is affordable as this type of simulations is generally light enough. For the 2D case, one already reaches a very good solutions for N = 9. It means that the number of split particles is We emphasized the fact that the assignment function of the set of split particles has to be as close as possible to the assignment function of the par-ticle to be split. This constraints is more important than the one generally considered [23], [24], only associated to its integrated values: mass/charge, impulsion, energy. This point has already been pointed out [25] while conserving the full distribution function needs a heavier numerical effort. From a computational point of view, the method we present has the advantage of having absolutely no overhead associated to the calculations of the coefficients (weight and locus) if these values are calculated and tabulated. In previous studies, the splitting process can be aborted if at least one of the split particle is out from the cell (on the fine grid) of the parent locus. Such a constraint has the drawback to violate the symmetrical constraint on the number and loci of the children, discussed in this paper. This requirement seems to us unsuitable, and should be re-evaluated in the light of the method presented here.
We have also emphasized the importance for a set of split particle to have an associated assignment function as close as possible as the one of the particle which is split. In Maxwell's equations, these particles are playing a role through their zero order (charge density) and first order (current density) moment of their distribution function. The definition of these quantities being linear with the assignment function, the errors on this quantities will evolve in the same way as the one presented for the assignment function. This is clearly the case for the charge density, but also for the current density as all the split particles have the same velocity as the one of the particle that is split.
Up to now, all the calculations have been conducted with the assumption that the grid is isotropic, meaning that the mesh size is the same in all directions. Actually, it is not a constraint, and the above w i and δ i values are also the optimum ones for anisotropic grid, whatever the degree of anisotropy.
For the two imensional cases, a simple homothety centered on the parent position with a ratio ∆ Y /∆ X only in the y direction is the link between the isotropic case with mesh size ∆ X and the anisotropic one with mesh size (∆ X , ∆ Y ). As this can be checked with the optimization code [21] (for both two and 3D cases), this transformation does not modify the optimization results, meaning that the optimum parameters are not depending on the degree of anisotropy of the grid.
These conclusions should also be considered in the opposite process of merging. Up to now, most of the existing methods identify a set of particles as close as possible in phase space, and merge two particles in a single one [26], or a larger set of particles in a pair [27]. In the merging strategy, the degree of the B-spline of the assignment function could also play an important role as in the splitting strategy, so the way to merge particles might also depend on this degree. To be more specific, the set of particles to be merged should be as close as possible in velocity space, but at an optimal finite distance in position space.
An important conclusion of this work is that when splitting a particle, the split particles should be dispatched at a given finite distance from the original particle, this distance depending only on the degree p of the B-spline and on the dimension d of the problem. We can draw an interesting parallel with a recent study on the merging problem by Luu et al. [28]. This study uses a Voronoi diagram to identify particles close enough in phase space to be merged. This algorithm needs a tolerance parameter, which is somewhat the threshold value below which one can consider that the particles to be merged are close enough. In this study by Luu et al., they show that this value has to be very small in velocity space, but can be much larger in position space.
While not a proof, it suggests the importance of the finite distance between particles involved in merging or splitting processes.
We also outlined the existence of an exact splitting at the cost of a large number of split particles. Considering the existence of efficient rezoning algorithm [25], one can wonder how the approximate splitting we present in this paper can compare with the exact one, followed by a rezoning procedure. Both options could be used in order to have the same number of split particles. The first option will be the cheapest in term of CPU, but their efficiency to preserve the distribution function and keep the fluctuations at the same level should be evaluated in a future work.