Strategies for particle resampling in PIC simulations
A. Muraviev, A. Bashinov, E. Efimenko, V. Volokitin, I. Meyerov, A. Gonoskov
SStrategies for particle resampling in PIC simulations
A. Muraviev a , A. Bashinov a , E. Efimenko a , V. Volokitin b , I. Meyerov b ,A. Gonoskov c a Institute of Applied Physics, Russian Academy of Sciences, Nizhny Novgorod 603950,Russia b Lobachevsky State University of Nizhni Novgorod, Nizhny Novgorod 603950, Russia c Department of Physics, University of Gothenburg, SE-41296 Gothenburg, Sweden
Abstract
In particle-in-cell simulations, excessive or even unfeasible computational de-mands can be caused by the growth of particle number in the course of prolificionization or cascaded pair production due to the effects of quantum electro-dynamics. Here we discuss how one can arrange a dynamic shrinking of themacro-particle ensemble to maintain an acceptable sampling of arbitrary par-ticle distribution. The approaches of merging and thinning as well as theirvariants are discussed and the aspects of their use are considered.
Keywords: particle-in-cell, resampling, merging, thinning, QED cascades
1. Introduction
In particle-in-cell (PIC) simulations and some other statistical computationsthe use of so-called macro-particles, that sample the distribution of some quan-tity, can be supplemented by the process of adding new macro-particles. Forexample, in PIC plasma simulations this can be done to account for the con-tinuous ionization of matter [1, 2] or for the electron-positron pair productiondue to the effects of quantum electrodynamics (QED)[3, 4, 5, 6, 7, 8]. Whensuch particle sources become prolific the number of macro-particles can growsignificantly. This can slow down the simulation or/and exhaust the memoryavailable for the allocation of macro-particles data. The natural solution is toresample the modeled distribution by a smaller number of macro-particles with
Preprint submitted to Elsevier June 16, 2020 a r X i v : . [ phy s i c s . c o m p - ph ] J un ncreased weight of their contribution. This procedure, referred to as down-sampling, can be repeatedly applied to combat the rise of the macro-particlesor, alternatively, to reduce the computational expenses for highly populated re-gions of phase space, where the representation became excessive with time. Inthe latter case, the released resources can be used to reduce the computationalnoise in underpopulated regions by introducing more macro-particles for bettersampling, i.e. performing so-called up-sampling.The implementation of down-sampling has been considered by several au-thors and a number of methods have been proposed. One can distinguish threemain approaches. According to the first approach, referred to as merging or co-alescing , one (or two) macro-particle is introduced to replace a subset of macro-particles that are close in the phase space. According to the second approach,referred to as thinning , we do not introduce new macro-particles, but removeone or several macro-particles and redistribute their weight among the others,either globally or locally within the given subset. Finally, one can totally re-place the given subset of macro-particles with a new subset of appropriatelyintroduced new macro-particles. This is the third approach, which is referredto as complete resampling in this article. In all cases, for the selected subset ofclosely located macro-particles the procedure can potentially lead to the changeof intrinsically conserved quantities, such as the total mass (weight), charge,energy or momentum. In addition, the procedure can potentially lead to thesudden change of grid values for the charge and current density, leading to ar-tificial noise, heating or other systematic effects. That is why the possibilityof preserving such quantities has been considered with special care by manyauthors.Within the approach of merging Lapenta and Brackbill proposed a methodfor coalescing two particles into one so that the charge assignment to the gridnodes is preserved [9]. Merging several macro-particles of dense clusters selectedwith the use of the Voronoi diagram has been suggested by Luu et al. [10]. Toconserve both energy and momentum Vranic et al. proposed to perform merginginto a pair of particles with the appropriately chosen momenta [11]. This can2e arranged using the selection of highly populated volumes in the phase space.A similar method with some modification has been used in Ref. [8].The approach of thinning provides various options, including rather straightforward ones. For example, Timokhin developed a procedure that repeatedly se-lects a random particle, deletes it and uniformly redistributes its weight amongthe others of the same kind [12]. Nerush et al. used a similar global thinningbut with the redistribution of mass, charge and energy [3]. Although this pro-cedure preserves the mentioned quantities globally, it implies their stochasticlocal sudden variations. One way of preventing such variations is to restrictthe redistribution to a dense cluster or a highly populated volume of the phasespace. A way to perform thinning with conservation of several arbitrary particleand grid quantities is proposed in Ref. [13].Within the approach of complete resampling, Lapenta and Brackbill pro-posed a method of replacing the macro-particles in a given cell with a newsubset, preserving the contributions to the grid quantities and also maximizingthe uniformity of the distribution of these quantities over the new subset [14]. Away of doing such resampling with the conservation of grid values for charge andcurrent density, as well as of the total energy of the resampled macro-particles,has been proposed by Assous et al. [15] and further developed in Ref. [16]. Pfeif-fer at al. proposed a statistical method that conserves momentum and energy[17]. Faghihi et al. reported on the development of an algorithm for both down-and up-sampling that preserves any number of particle and grid quantities [18].Recently, the methods of down-sampling became highly demanded for thenumerical studies of QED cascades that will be inherent for the upcoming ex-periments at the next generation high-intensity laser facilities [19, 20, 21, 22].The numerical studies indicate that, apart from the drastic increase of particlenumber by several orders of magnitude, the physics in strong laser fields in-cludes a variety of new phenomena [23, 24, 25, 26, 27, 28, 29]. The absence ofprior knowledge about the minimal scales of new phenomena raises a new diffi-culty for the implementation of down-sampling: merging macro-particles withindense clusters or volumes of predetermined scale may erase or affect smaller3cale peculiarities that are essential for the modelled phenomena. Although thecoordinates have the cell size as a natural limiting scale, the momentum doesnot have any natural resolution limit according to the PIC method. At the sametime, shrinking the momentum gate for merging in a given cell, increases theneeded density for the selection of several particles to be possible.In this article we consider how one can use the thinning approach or modifythe merging procedure to combat the outlined difficulty. In addition, we considerthe aspect of reducing the difference between the weights of particles as a wayto increase the efficiency of sampling. For our study we develop, compare andanalyse several methods that we release as open-source tools available withinthe hi- χ famework [30]. The article is arranged as follows. In section 2 weintroduce the principle of agnostic down-sampling for developing the thinningprocedures that are applicable without prior knowledge about minimal scalesattributed to the modelled process. In section 3 we describe several thinningmethods that comply this this principle and also discuss how to modify themerging methods to improve their applicability in this context. In section 4we perform a basic comparison of the methods. In section 5, we show realisticexamples. We conclude in section 6.
2. The principle of agnostic down-sampling
When several macro-particles are merged we assume that all these particlessample a uniform part of the particles distribution without any complex trendsor peculiarities represented by some of them. For example, the macro-particlesshown in Fig. 1 are assumed to sample the distribution shown with the solidcurve, not the ones shown with the dashed or dotted curves. Although one canprobably mitigate the risks of such misinterpretation by restricting the choiceof particles in a cluster to sufficiently small phase space volumes, let us consideran alternative down-sampling methodology that is applicable without any priorknowledge and thus called agnostic here.Firstly, note that simple merging into the average position of the selected4 hase space coordinate w e i g h t macro-particles Figure 1: An illustrative clarification of the difficulty in arranging non-destructive mergingwithout prior knowledge about the peculiarities in the particle distribution and their minimalscales. The merging of four macro-particles (grey rectangular peaks) can be non-destructiveif they sample the distribution show with black solid curve. However it can affect the mod-elled processes if the sampled distribution has a more complex shape. The examples of suchpotential shapes are shown with dashed blue and dotted green curves. macro-particles will likely reduce the variance of the peaks that are more nar-row than the selection scale. To mitigate this one can introduce a probabilisticmismatch with the variance determined over the merged macro-particles. How-ever, such a procedure will likely cause filling the gap in the distribution shownwith the dashed curve in Fig. 1. Looking at the distribution shown with thedotted curve we can see that the only way to not introduce any particles inany potentially empty phase space region is to use only the existing particles.In other words, we should try to restrict our action to the removal of one orseveral macro-particles in combination with the change of the weights of theother macro-particles.Secondly, we should not affect any distribution functions. In order to dothis, we can arrange a probabilistic procedure so that the chance of removingany particle in the given subset is compensated by the chance of increasing itsweight. If the expectation value for the weights of each particle exactly equal toits initial weight, all the possible distribution functions remains unchanged onaverage (see rigorous consideration below). This means that if some peculiarityis erased in one subset, it has a chance to be increased in a different subset.5n this way, the procedure only increases the statistical variation but does notremove peculiarities at any scale.We can now formulate the principle of down-sampling that is based on thisidea. A down-sampling is called agnostic if it is restricted to the probabilisticchange of weights so that (1) at least one macro-particle receives zero weightand can be removed and (2) the expectation value for the weight of each macro-particle is exactly equal to its initial weight.If we have initially n macro-particles and the weight of i -th is w i , the principleof agnostic down-sampling implies determining a number of outcomes so thatthe probability p k of choosing the k -th outcome provides (cid:104) ˆ w i (cid:105) = (cid:88) k ∈ outcomes w ki p k = w i for i = 1 , ...n, (1)where (cid:104) ˆ w i (cid:105) is the resulting weight averaged over the outcomes, and w ki is theweight assigned to the i -th particle in the k -th outcome.Let us demonstrate that any agnostic resampling has the property of pre-serving all the distribution functions on average, i.e. any distribution averagedover the outcomes of probabilistic resampling procedure coincides exactly withthat we had before the resampling. In order to show this, we consider an ar-bitrary distribution F = ∂N/∂f , which can be numerically represented as asequence of values F j : F j ( A ) = 1 V ( d j ) (cid:88) f ( x i , p i ,σ i ) ∈ d j w i , (2)where f ( x i , p i , σ i ) is the quantity (can be multidimensional) over which thedistribution is computed, σ i is a generalized vector of the particle’s parametersother than coordinate and momentum (e.g. spin, polarization, etc.), A is thestate of ensemble, d j denotes the sub-regions used to discretize the space of f values and V ( d j ) is the volume of j -th sub-region. In other words, F j ( A ) isthe number of real particles, for which the quantity f falls into d j , divided bythe volume of the element d j . As one can see, in such a way we can definecoordinate-, momentum-, energy-, angular- and other distributions on uniform6s well as non-uniform grids defined by d j . We can now formally compute thevalue F j averaged over outcomes of the resampling procedure that turns A into R ( A ): (cid:104) F j ( R ( A )) (cid:105) R = 1 V ( d j ) (cid:88) f ( x i , p i ,σ i ) ∈ d j (cid:104) ˆ w i (cid:105) R = F j ( A ) . (3)Here the first equality is provided by the fact that an agnostic resampling neitherchanges internal state ( x i , p i , σ i ) of nor adds new macro-particles, whereas thesecond equality follows from the fact that it preserves the weight on average,i.e. ∀ i : (cid:104) ˆ w i (cid:105) R = w i . As one can see, we proved the statement without requiringany knowledge about (1) the quantity f , (2) the numerical intervals d j and (3)the distribution of macro-particles (either over the position in phase space orover the values of their weights). In this context, the term agnostic indicatesthat the procedure preserves all the distribution functions independently of theoutlined entities.
3. Strategies of agnostic down-sampling
In this section we propose a number of methods that comply with the prin-ciple of agnostic down-sampling. We start with the simplest methods and thendescribe more advanced ones and their potential benefits. We quantify the rateof resampling by a parameter k being the target ratio of number of macroparti-cles in the ensemble before resampling to their number after resampling. Eachmethod was given a shortened name in parenthesis for designation on graphics.1. Simple thinning (simple).According to this method each macro-particle is either removed, with equalprobability p = 1 − k − , or has its weight increased by a factor of k . This methodis agnostic because (cid:104) ˆ w i (cid:105) = 0 · (cid:0) − k − (cid:1) + kw i · k − = w i . The method does notstrictly conserve any quantities such as total weight, energy or momentum, butconserves all quantities on average as any agnostic method. The method can beapplied to subsets of any size, which makes it possible to apply it to very smallvolumes in phase space. This method is the easiest to implement and analyze7heoretically. If the initial total number of macroparticles is n , the number ofmacroparticles remaining after resampling is approximately n/k .2. Leveling thinning (leveling).In this method, we first calculate the average weight ¯ w of particles in a givencell . Then, for all particles with weight w i < k ¯ w the weight k ¯ w is assigned withprobability w i / ( k ¯ w ) and the particle is removed otherwise. This method isagnostic because (cid:104) ˆ w i (cid:105) = k ¯ w · w i / ( k ¯ w ) + 0 · (1 − w i / ( k ¯ w )) = w i . It is clearthat this procedure gets rid of macro-particles with weight below kW . Thismay help to balance and optimize the distribution of computational resources.The method does not strictly conserve any quantities. The actual number ofremaining macroparticles can be less or greater than n/k depending on theinitial weight distribution among macroparticles. If all macroparticles initiallyhave the same weight, the number of macroparticles remaining after resamplingis approximately n/k .3. Global leveling thinning (globalLev) is a modification of the previouslydescribed leveling thinning method. The method works similarly, except it usesall particles in the entire calculation area to compute the average weight W .4. Number-conservative thinning (numberT).In this method we select random macro-particles with probability propor-tional to their weight, i.e. w i /W , where W = (cid:80) w i . We repeat this selection m times and count the number of times c i we have chosen the i -th macro-particle.After that the particles that have not been selected even once ( c i = 0) are re-moved and the others are assigned with a new weight equal to ˆ w i = c i W/m .It is clear that the i -th macro-particle will be selected (cid:104) c i (cid:105) = mw i /W times onaverage and thus the mathematical expectation of the change in the macropar-ticle’s weight is zero: (cid:104) ˆ w i (cid:105) = ( (cid:104) c i (cid:105) W/m ) = w i , i.e. the method is agnostic.This procedure strictly conserves the total weight of macro-particles in a cell:ˆ W = (cid:80) ˆ w i = (cid:80) c i W/m = W . In addition, it favors the removal of macro-particles with small weight. This may also contribute to the efficiency of sam-pling. The number of macro-particles ˆ n after this procedure is probabilistic butobviously cannot exceed m . The average number of remaining macro-particles8s given by the expression ˆ n = n (cid:88) i =1 (cid:16) − (cid:16) − w i W (cid:17) m (cid:17) . (4)If we assume that in our distribution macro-particles have similar weights wecan estimate ˆ n ≈ n (cid:0) − (cid:0) − n − (cid:1) m (cid:1) . Assuming also that n is large, we canestimate that ˆ n = n/k is achieved for m ≈ − n ln (cid:0) − k − (cid:1) . (5)This means, for example, that for large n if we need to remove roughly half ofthe macro-particles we need m ≈ ln (2) n . This method is useful when the totalcharge/number of particles needs to be strictly conserved.5. Energy-conservative thinning (energyT) is a modification of the previouslydescribed number-conservative thinning . According to this method we also selecta random macro-particle m times, but do this with probability proportional toenergy, i.e. the i -th macroparticle is selected with probability e i w i /E , where e i is the energy of the particle represented by the i -th macroparticle and E = (cid:80) e i w i . The macroparticles that have been selected c i (cid:54) = 0 times are assignedwith the weight c i E/ ( e i m ) and the others are removed. One can check thatthis procedure complies with the principle of agnostic down-sampling and alsostrictly conserves the total energy E in each cell. However, this method doesnot strictly conserve the total weight W . This method is useful when the totalkinetic energy of particles needs to be strictly conserved.6. Conservative thinning (conserv). This method, described in [13], canbe configured to preserve several invariants simultaneously. Each invariant ( A )can be represented by a linear equation: A = (cid:80) a i w i , where a i amd w i are thecontribution and the weight of the i -th particle, respectively. The conservationof several invariants sets a system of linear equations, where the number ofvariables (weights w i ) can be controlled by the number of particles involved inthe thinning procedure. If the number of particles is greater than the numberof invariants, the system is undetermined. It turns out that it is possible tofind two solutions with one of the weights being equal to zero and others being9ositive, so that the probabilistic choice of one of these solutions results in anagnostic resampling that reduces the number of macroparticles by one. Theprocedure can be repeated several times for a given set of macroparticles toreduce the number of particles down to n/k (assuming that it is still larger thanthe number of invariants). This method is useful when a number of physicalproperties of particles need to be strictly conserved.For comparison we also consider methods that do not comply with the prin-ciple of agnostic down-sampling. Most of these methods revolve around mergingof dense clusters.7. Merging to averaged location (mergeAv). According to this method, foreach cell we determine n cell /k (where n cell is the number of particles in a cell),but at least 3, clusters using the k-means method with respect to location ofparticles in the momentum space. Next we replace all particles in each clus-ter with a new macro-particle that has the mean coordinate and momentum ofparticles in the cluster and weight equal to the total weight W . The numberof selected clusters determines the number of remaining particles. Since thecomplexity of the k-means method is ( n cell /k ), the algorithm may require sig-nificant computational resources. This method is useful when the phase spacecan be adequately represented by a number of dense clusters. If n cell /k >> n/k . Due tocomputational time restrictions, the recommended value of k for merge-basedmethods is such that no more than 30 clusters are formed.8. Merging to random particle (merge). As we mentioned earlier, the merg-ing procedure can naturally result in the systematic relocation of particles to-wards denser regions. An indicative example is the case of a bulk of particles (ora particle beam) with a narrow distribution in coordinate space (as comparedto the coordinate scale of clusters used for merging). In this case merging natu-rally favors bringing macro-particles to the peak of that distribution, removingthe macro-particles in its tails. We can avoid this by introducing the followingmodification to the previous method. The weight of particles in each determinedcluster is brought to a random particle in the cluster. Note, however, that this10oes not prevent the reduction of particles’ spread in the momentum space.As a general note we would like to highlight the following. Any down-sampling results in the loss of information since the amount of unique macro-particles decreases. The inevitable consequence of this is the increase of noise inthe distributions of particles. The more macro-particles are used the less noisewe can expect and vise versa. In practice this leads to the trade-off betweenthe accuracy of results and the computational demands. The goal of arrangingappropriate resampling is to avoid systematic deviations and minimize compu-tational demands for reaching the accuracy necessary for the problem of interest.To not spend computational resources for resampling when it is not needed,we use the trigger for starting of the resampling procedure: the number ofmacroparticles in a shared-memory computational domain must reach a certainthreshold value. All methods except global leveling are applied to each cellindependently.
4. Comparison of methods on model problems
In this section we present the comparison of the resampling methods de-scribed above using model problems where QED effects are negligible: a steady-state plasma and the development of a Weibel instability in two counter-streamingplasma flows.
In our numerical experiments we observe an additional decrease of plasmatemperature caused by the ensemble resampling. To clarify the reason for thiseffect, we start from a brief phenomenological analysis.For simplicity we consider the resampling of electrons, whereas the ions(or positrons) are modelled by a uniform positively charged background. Wecan outline two basic reasons for the change of temperature to happen due toresampling. Firstly, if the applied resampling does not preserve the total kineticenergy of particles and the method is not agnostic, then we can potentially11ave an asymmetric net acquisition of energy mismatches, either positive ornegative. Secondly, even if the used resampling is agnostic, the change of weightseffectively means that we draw the plasma out of the equilibrium state. Therecould be two cases.If the Poissons equation is solved at each iteration (as it happens in somespectral codes), this abrupt local relocation of charges builds up an additionallocal variation of electric field, the energy of which can eventually add up tothe temperature increase. Note that if the resampling does not preserve chargewithin each cell, the local variation of charge density can build up a non-zeroglobal electric field (especially in the 1D case), which can have a significantenergy. This explains why preserving charge may be beneficial.In case the Poissons equation is not solved and the field is only affectedby the charge currents, the local removal of macro-particles would effectivelymean adding compensating charges that are fixed in space (the added positivecharge is compensated by the increase of weight of the remaining macroparti-cles). In this case the plasma will tend to a new equilibrium state relative tothe positively charged background with corresponding local variations of chargedensity. Again, if the charge is not preserved locally, a global effective potentialcan be formed and the placement of compensating charges can cause a signifi-cant energy change. Now, let us imagine a situation when the plasma is leavingsome part of computational region. The added effective positive charges in thispart will show up as unchanged noise in electric field. Since this field wouldcontain a strictly positive energy, we can conclude that this energy is effec-tively deducted from the thermal motion of charges due to the application ofresampling. Hence, resampling can cause an effective cooling of plasma. Thisis observed in our numerical experiments.In order to compare how strongly different resampling strategies can affectsimulations we choose to compare the change of temperature of a steady stateplasma after the application of resampling using different methods.Particularly, we consider a 3D region represented by 32 × ×
32 cells withperiodic boundary conditions filled with homogeneous quasineutral electron-12 igure 2: Change of equilibrium temperature for a plasma with initially N = 100 ppc aftera single instance of resampling depending on the target resampling coefficient k , agnostic methods. Dashed lines: methods’ results, Solid: linear fit. positron plasma with initial temperature T = 0 . mc (here m is the electroncharge, c is the speed of light), cell size equal to 2 Debye radius, and physical density derived from these values. The considered values of the initial numberof macro-particles are N = 100 and N = 1000 particles per cell (ppc). Thetime step is set to 1 /
128 of the period of cold plasma oscillations. After t = 1oscillation period we perform resampling, at t = 10 oscillation periods we calcu-late the temperature of the plasma (as the average kinetic energy of particles)in comparison to the initial temperature T . This procedure was performed forvalues of N mentioned above for every method of resampling and for a set ofresampling coefficients k, which indicate the target decrease ratio in the amountof macroparticles, equal to (1.1; 3; 10; 30; 100; 300; 1000) where possible dueto method limitations. To identify the temperature decrease induced by resam-pling the temperature decrease in the case without resampling is subtractedfrom the value obtained using various methods.For all agnostic methods the results are similar and lie within a narrow rangeof each other. These results are shown also in Fig. 2 and 3 for the cases with N = 100 and N = 1000, respectively.13 igure 3: Change of equilibrium temperature for a plasma with initially N = 1000 ppc aftera single instance of resampling depending on the target resampling coefficient k , agnostic methods. Dashed lines: methods’ results, Solid: linear fit. Interestingly, the results show a linear trend in ∆ T ( k ). In addition, closetrends with respect to N /k values (see Fig. 2 and 3) indicate that the tem-perature decrease depends solely on the number of macroparticles per cell re-maining after resampling rather than on the initial ppc number and the coef-ficient k separately. The results of agnostic methods can be roughly fitted by∆ T = − (0 . /ppc f ) ∗ T , where ppc f is the final number of macroparticles percell after resampling. Actual values of ∆ T for each method vary by about 10%from this rough estimate depending on the particular method in question andthe value of k.Merge-based methods, on the other hand, show considerably poorer perfor-mance according to our chosen metric (see Fig. 4 and 5). Even in the best-casescenario where merge methods perform the closest to agnostic methods, ∆ T shown by merge methods is approximately 15-20 times greater than ∆ T shownby agnostic methods. Particularly, while agnostic methods show a linear trendtowards a 10% temperature drop at k = N (which means the target numberof particles after resampling is N /k = 1 ppc), merge methods yield a whole65% temperature decrease already at ppc=3. The curve for merge methods is14 igure 4: Change of equilibrium temperature for a plasma with initially N = 100 ppc aftera single instance of resampling depending on the target resampling coefficient k , all methods.Dashed lines: methods’ results, Solid: linear fit.Figure 5: Change of equilibrium temperature for a plasma with initially N = 1000 ppc aftera single instance of resampling depending on the target resampling coefficient k , all methods.Dashed lines: methods’ results, Solid: linear fit. slightly concave up, so for lower k the result is even worse in relative comparisonto the agnostic methods.We conclude that from the point of view of numerical cooling, merge-basedmethods are at least an order of magnitude worse (i.e. these methods yield15umerical cooling an order of magnitude stronger) than agnostic methods. The second test problem is the development of Weibel instability [31] incounter-streaming plasma flows [32]. This instability results in an exponentialgrowth of perturbations in plasma density, current and magnetic field along thedirection transverse to the plasma stream. To make our experiment robust weintroduce a periodic modulation of density in the transverse direction to act asa systematic seed for the instability. For each method we carry out an individ-ual simulation and perform a single resampling procedure near the beginningof the growth. In such a way we intend (1) to see whether the introducedrandom perturbations in density can affect the process, and (2) compare differ-ent resampling methods in terms of the introduced perturbations. To quantifythe strength of these perturbations we measure the variance of plasma densitycomputed for individual cells of the computational grid.Let us first note that the reduction of the number of macro-particles shouldnaturally result in the increase of variance for the number of real particles ineach cell. The extent of this increase, however, depends on the method. Forexample, number conservative thinning does not change the number of realparticles (although the variance can grow after the migration of particles be-tween the cells). In this cases the impact of resampling on macro-particlesis coordinated within each cell. To estimate the worst case scenario, let usconsider the case of simple thinning , for which the impact is totally uncoor-dinated. With probability p = 1 /k the particles weight w is increased by afactor of k to w = kw , otherwise (with probability p = 1 − k − ) the macro-particle is deleted. The expected value of the number of physical particles N phys among different realizations of this random process must remain un-changed: E [ N phys ] = k w = w (hereafter by E [ · ] we denote the value averagedover the realization of resampling). For the contribution of individual macro-particles, we can compute the variance D [ N phys ] = E [ N phys ] − E [ N phys ] ,and E [ N phys ] = (cid:80) p i w i , where p i is the probability of the i -th outcome and16 physi = w i is the number of physical particles in that outcome. For thesimple thinning we can write E [ N phys ] = k − w + (1 − k − ) · kw . Fi-nally, the variance is D [ N phys ] = kw − w = ( k − w . Since the varianceis additive, considering a cell with N macroparticles with factor w , we obtain: D [ N phys ] = N ( k − w . We are interested in macroparameters, such as particledensity n = N phys ∆ V − = N w ∆ V − , where ∆ V is the cell volume. We cancalculate D [ n ] = D [ N phys ](∆ V ) = N ( k − w (∆ V ) = n ( k − w ∆ V = n ( k − N (6)Although this expression is the variance of n over different realizations of randomevents, the independence of such events in different cells allows us to use it tocalculate the dispersion of n over coordinate space.For our study, we performed a series of 2D simulations of the Weibel instabil-ity development in counter-streaming flows of quasineutral electron-ion plasma.We considered the following parameters: initial density n = 10 cm − , plasmaflow velocity V ± = ± . c , which corresponds to a Lorentz-factor of γ =100, where c is the light velocity and the ”+” and ”-” signs denote the streamsdirected along and opposite to the x axis in our simulations. The initial densityof ions in both streams was set to be uniform. The initial density of electrons andtheir local momentum were modulated harmonically across the transverse di-rection ( y axis): n ± = ± δn cos( k y y ), p y, ± = ± δm e ω p V ( k Γ) − c − sin( k y y ),where δ = 0 .
02 is the modulation amplitude, k y = πL , L = 2 . · − cm is thespatial scale of the modulation, ω p = (cid:112) πe n /m e is the plasma frequency (oftotal density of both streams), m e and e are the electron mass and charge. In theconsidered case of small-scale spatial modulation, the growth rate of the Weibelinstability is Γ ≈ ω p γ − / V /c . The modulation of electron density leads to theelectromagnetic field variation of the following form: (cid:126)B = 8 πδ en k V + c sin( k y y ) (cid:126)z , (cid:126)E = − πδ Γ ck en k V + c cos( k y y ) (cid:126)y .The size of the simulation region was 2 µ m × µ m and the region was rep-resented by 96 ×
384 cells. The boundary conditions were periodic. The timestep was set to 1 / (128Γ). The computation time was set to match the duration17 igure 6: ∆ D ( t ) for different methods of resampling and k = 50. Group 1: simple (darkblue), leveling (green), global leveling (blue); Group 2: number thinning (red), energy thinning (cyan); Group 3: merge (purple), mergeAv (yellow). of the linear regime during which the plasma density perturbation is negligiblecompared to the plasma density itself, which was the case until t ≈ . / Γ.The computation was performed for each method of resampling except conserv and each value of the resampling coefficient k from the set (1.1;2;5;10;20;50), aswell as for the case without resampling. In every computation the resamplingprocedure was applied once at t = 1 . / Γ. In order to identify the dispersioninduced by the resampling procedure, for each method and each value of k the difference ∆ D ( t ) = D [ n ]( t ) − D [ n ]( t ) was calculated, where D [ n ]( t ) is thetime dependence of the variance of particle density n in that calculation and D [ n ]( t ) ∼ e t is the density variance in the simulation, which was entirelyperformed without resampling.Let us first compare the results of different methods using equal values of k = 50, see Fig. 6. According to these results, the methods can be dividedinto 3 groups. Within each group the results are very similar. Specifically, for The methods merge and mergeAv could not complete for coefficients k = 1 . k = 2due to computational time restrictions. simple, leveling, Globallev ) ∆ D ( t ) takes the form close toa step function. In other words, resampling causes a single leap in ∆ D ( t ) atthe time of resampling, after that the difference in dispersion compared to thecase without resampling stays constant up until the end of the linear regime,despite the fact that the variance itself grows exponentially. This result confirmsthat in this case the influence of resampling on the variance can be consideredindependently of other (physical) processes affecting it.For methods in Group 2 ( numberT, energyT ) the results are of the sameorder of magnitude and follow a similar trend, but these results show a notableperiodic oscillation in ∆ D ( t ), which is not negligible, but still of low amplitudecompared to the value itself.Group 3 methods ( merge, mergeAv ) show a considerably lower ∆ D ( t ) (es-pecially the merge method), but oscillations are of the same order of magnitudeas the value itself.Now let us compare the results of each method depending on the value of k .In order to assess that, we present values of ∆ D ( t ) immediately after resampling.Since before resampling all computations in our series are identical, this valueis exactly the leap in dispersion caused by resampling. We present this valuedepending on the value of the resampling coefficient k for all resampling methodsexcept conserv , as well as the estimate (6) for the simple thinning method, seeFig. 7.As evident, Group 1 results follow the estimate (6): a linear dependenceproportional to k −
1. The results in group 2 have a slightly lower slope: com-pared to Group 1 the variance increase is overall lower, with except of somerange of low values of k . The merge method performed close to linear in therange of coefficients where we were able to complete it, and yields notably lowerincrease than all other methods: the slope is ∼ / mergeAverage method stands out as highly nonlinear, yieldinghigher increase of variance at low coefficients, but for k = 50 the increase islower than that of Group 1.The lowest increase of variance is systematically provided by the method19 igure 7: ∆ D ( k ) immediately after resampling for different methods. Group 1: simple,leveling, global leveling ; Group 2: number thinning, energy thinning ; Group 3: merge (purple), mergeAv (yellow). merge . This can be related to the fact that this method provides the mostadvanced coordinated replacement of macroparticles giving, in some way, mostefficient representation of the ensemble. Nevertheless, we should highlight thatthis is possible because we know the coordinate and momentum scales of thesimulated processed and ensure that they are larger than that of the mergeprocedure. Apart from this, we must note that even for the highest coefficient k = 50 the computational time of merge/mergeAv resamping exceeds that ofother methods by at least a factor of 7. For k = 10 a single application ofresampling doubles the total execution time of the whole simulation, and for k = 5 resampling takes up 80% of the computational time. Such characteristicsare highly dependant on the particular problem at hand, of course, but it canbe firmly said that merging methods may be highly demanding as compared tothe other methods under similar accuracy requirements.
5. Comparison of methods on pertinent physical problems
A large class of tasks that require resampling methods is the study of cas-caded production of electron-positron pairs and high-energy photons in laser20elds of high intensities. The description of these processes is based on therates of transitions between quantum states computed within quantum electro-dynamics (QED) and thus the cascades are commonly referred to as QED cas-cades. One of the most widely used numerical approaches for the simulation ofQED cascades is based on extended PIC codes, also known as QED-PIC codes.With the development of an electrodynamic cascade, the number of particlescan increase exponentially in time and significantly increase the computationaldemands. This means that some problems are impossible to compute withoutthe use of resampling, so there is no benchmark to compare results to. In thiscase we have to rely entirely on results acquired using one or more methods ofresampling.A cascade can have two stages: linear and nonlinear. At the linear stage,the density of the generated electron-positron plasma is not sufficient to signifi-cantly affect the structure of the electromagnetic field and its intensity. At thenonlinear stage, on the contrary, the generated plasma has a higher density andsignificantly affects the electromagnetic field. Therefore, at the linear stage re-sampling can affect only the particle distribution function, affecting, potentially,the local rate of cascade development, or the particle distribution function. Atthe nonlinear stage resampling may also affect the structure and magnitude ofthe fields, changing the plasma-field dynamics. Below, using several examples,we consider various stages of the QED cascade and show how different typesof resampling can affect the simulated processes. For simulations, we used thePICADOR code with the
Adaptive Event Generator module described in [7],which, when necessary, subdivides the time step in order to resolve the QEDcascade.
To investigate the operation of resampling methods at the linear stage of thecascade, we chose a well-studied problem of QED-cascade development in thefield of a standing linearly polarized plane wave (see, for example, [33, 34, 35]).In order to have a benchmark for comparison - results of a simulation without21he use of resampling (’w/o’) - we chose a relatively small wave amplitude E =1000 mω ce , where ω = 2 . × s − is the laser frequency for the wavelength λ = 0 . µ m. The wave has only E z and B y non-zero field components, whichcauses particles to move along the x and z axes. At this wave amplitude,electrons and positrons tend to the vicinity of the nodes of the electric field,but due to the stochasticity of photon emission and a photon decay into pairs,electrons and positrons can reach the vicinity of the electric field antinode [36].Initially, electrons and positrons with the density of approximately 10 cm − are uniformly distributed in a λ × λ × λ simulation box represented by 128 × × T ,where T = 2 .
66 fs is the wave period. The time step was dt = 1 . × − s.The same parameters were used for all resampling methods: every seconditeration, starting from zero, if the amount of macroparticles of any type exceedsthe resampling threshold of 5 × particles, the particles of that type undergoresampling with k = 2, and the amount of these macroparticles decreases byapproximately half.To analyze the effect of resampling on the accuracy of each simulation, weconsider the temporal evolution of the total number of electrons and positrons N e (Fig. 8(a)). We used the following parameters for our analysis. The first oneis the cascade growth rate Γ = ln Ne ( t =7 T ) Ne ( t =3 T ) T , the second one is the relative mean-square deviation of N e ( t ) acquired with resampling from N e ( t ) w/o - acquiredwithout resampling - η = (cid:115) (cid:82) T Ne,res − Ne,w/o )2 N e,w/o dt T , where res denotes a certaintype of resampling. In order to calculate Γ, the moment t = 3 T was used asthe initial one, since at t = 0 there are no photons and it takes about 2 . T forsteady exponential growth to establish.Note, that the principle of agnostic resampling ensures that the growth rateis not affected on average, but for this we need to know that the process is22 igure 8: Comparison of different resampling types via simulations of the QED cascade devel-opment in a standing linearly-polarized plane wave. (a) Time dependence of the total numberof electrons and positrons N e in the simulation box. (b) The cascade growth rate Gamma .(c) The relative mean square deviation η of N e ( t ) acquired with use of resampling from N e ( t )acquired without resampling. The colors in figures (b) and (c) correspond to the color of thelines in figure (a). exponential and compute it using the geometric mean (instead of just mean )value of density over several identical simulations (with different random seed).That is why we here mimic the resampling-related error in the computation ofoutlined parameters without any prior knowledge about the process. This alsoquantifies possible resampling-related distortions in more complex processes aswe explicitly observe in other examples (see section 5.3.).It should be noted that at a given wave amplitude, electrons and positronsemit plenty of photons, but only a small fraction of all photons decays intoelectron-positron pairs. Therefore, during the whole simulation resampling wasinitiated about 100 times for photons, but only once for electrons and positrons(Table 5.1). Thus, hefty macroparticles are added to the ensemble of electronsand positrons as the result of photon decay. The simulation without resamplingshows that Γ T = 0 .
445 and the number of electrons and positrons increased byabout a factor of 12 during the whole simulation (Fig 8(a)).Based on the comparison of Fig. FCsGr (a), (b) and (c), all consideredresampling methods can be divided into 4 groups. The first group includesthe most accurate thinning method leveling with η = 0 .
006 and an error inthe value of Γ of about 0.3%. The second group also includes fairly accurate23hinning methods globalLev , conserv and energyT , for which η ≈ . merge , mergeAv and thinning method numberT with η ≈ .
08; 0 .
08; 0 . η = 0 . merge and mergeAv methods is significantly higher - up to about 11000 seconds. Thisis primarily due to the use of the ’k-means’ method for merging of particles,and also due to a more frequent triggering of resampling. The merge methodsyield a speedup of about 2 times in comparison with the ’w/o’ simulation, andthey have the same or even worse accuracy than thinning methods (the method’simple’ being an exception). Thus, in this setup, thinning methods, providingat least equal accuracy (or much greater accuracy in some cases), noticeablybenefit in performance in comparison with the merge methods. In addition,thinning methods trigger the resampling procedure less frequently. It is espe-cially worth noting that the most accurate method for the considered problemis the leveling method. It provides a fractions of a percent accuracy and a morethan twentyfold simulation speedup. 24 able 1: The influence of different types of resampling on run time and frequency of resam-pling. Type ofresampling
Run Time, s
Resampling ofPhotons Resampling of e − , e + conserv globalLev energyT leveling mergeAv merge numberT simple In this problem we perform a full 3D3P simulation of irradiation of a seedplasma target by two colliding half-infinite linearly-polarized ultraintense laserpulses. As shown previously [37], in such a configuration the electromagneticcascade [3, 4, 5, 6, 7, 8] can result in a rapid generation of plasma on thecascade’s linear stage followed by the formation of ultra-thin current sheets onthe nonlinear stage. We perform the simulation for this problem using differentmethods of resampling.It is important to note that since in this setup plasma density can increasesignificantly over a single half-period of the electromagnetic pulse, the setupis highly sensitive to the phase of the standing electromagnetic wave at themoment in time when the plasma density becomes sufficient to affect the fieldstructure of the standing wave. Consequently, the probabilistic nature of elec-tromagnetic cascading may lead to a discrepancy in results (in physical andnumerical experiments alike) over different realizations of probabilistic physical events. We stress that a slight discrepancy in observed parameters in calcula-25 igure 9: Total Field Energy as a function of time for different methods of resampling tions with the same initial conditions may represent different realizations of aprobabilistic physical process and it alone does not indicate that the resamplingprocess leads to the distortion of the result.The chosen initial parameters for this problem are: laser pulse wavelength λ = 800 nm, field amplitude a = 3500 in relativistic units, width D = 5 λ ,initial seed plasma density at 10 − of the critical plasma density N cr .Below we present the total electromagnetic field energy in the 7 λ × λ × λ simulation box for different methods of resampling as a function of time sincethe start of the nonlinear stage, during which the current sheets begin to form.Agnostic methods were compared using the same default value k = 2. Thisvalue of k is well outside the recommended range for merge-based methodsdictated by computational time restrictions, so here merge-based methods werenot considered.As shown in Fig. 9, agnostic methods excluding simple yield qualitativelysimilar results varying within 5%, thus cross-confirming each other’s results.For the method simple we performed five attempts with different seeds of therandom function. Nevertheless, all these attempts have been terminated at anearly stage with a clearly unphysical surge in field energy.26eeper investigation of the performance of the method simple shows thatalthough average density is conserved well (up to a point), maximal density isnot (also see Sec. 5.3). Since each particle’s weight is probabilistically increasedor zeroed independently, simple resampling with coefficient k may result in thetotal weight in a certain cell increasing by a factor of k . This phenomena is a lotmore likely to occur if much of that cell’s weight is carried by a single particle orseveral particles. If k is large enough, or if the procedure is applied many times(our case), this can result in a significant increase of the total particle weightin some cells (effectively at the cost of other cells). Since on average propertiesare conserved, intuition may suggest this should not lead to physically incorrectresults, apart from some increased numerical noise. However, it can be shownthat this is not the case. The reason behind this lies in the discrete natureof the particle-in-cell code. To operate correctly, all physically relevant timescales have to be resolved by the PIC time step. The artificial increase oftotal weight (and thus physical density) in some cells leads to the increase ofthe local plasma oscillation frequency. If this frequency exceeds the frequencyresolved by the PIC code’s time step, an unphysical instability develops. Thecurrent generated by a very hefty particle or cluster of particles creates anelectric field that inverts the momentum of these particles in the next iteration.The process is additionally fed by new particles created in an electromagneticcascade, which results in an unphysical exponential growth in particle density,energy and field values, resulting in the termination of the computation, oftenin as few as several iterations. Other methods either strictly conserve certainvalues or exclude particles with a large weight from the procedure, helping toavoid this problem.To summarize, we see that the simple thinning may not be applicable if alocal stochastic variation of density can cause numerical instability. All otheragnostic methods show adequate results similar to each other in this setup.27 .3. Pinching of electron-positron plasma in a multi-10PW dipole wave To study the effect of various resampling techniques on the dynamics ofthe particle ensemble in another pertinent problem, we have investigated theinteraction of a multi-10PW-level laser radiation with plasma targets. Such aproblem is of great interest due to the fact that this kind of laser systems willsoon be ready to be used in experiments [38, 39, 40] and even more powerfulsystems are being developed [41, 42]. These experimental setups will allow thestudies of QED cascade development in converging fields of petawatt powerlevel, which is discussed in the previous section. Moreover, these laser systemswill be capable of driving various self-consistent nonlinear regimes during theinteraction of a QED-produced electron-positron plasma with ultraintense fields.As shown earlier [28] and [29], in this configuration there are two mainnonlinear regimes of interaction. In the first regime, which is realized when thelaser power is less than 20 PW, thin electron-positron plasma sheets are formedas a result of the development of a current instability [28]. If the thresholdpower exceeds 20 PW, pinching of electron-positron plasma is possible as aresult of current contraction [29]. In this paper in order to study the influenceof resampling method we consider the interaction of an ideal dipole wave withtotal power P = 27 P W with a low-density plasma target, which acts as a seedfor the development of a QED cascade. The choice of this value of power canbe justified by two factors. First, the development of a QED cascade in sucha configuration was discussed in detail in [29], and, second, such a formulationallows studying the dynamics of the system at both (linear and nonlinear) stagesof evolution.In the numerical simulation, performed with the PICADOR Particle-in-Cellcode, the interaction of an ideal 27 PW dipole wave with a plasma target witha 3 µ m diameter and a density of 10 cm − was simulated. The size of thesimulation area was 4 x 4 x 4 µ m with a grid size of 512 x 512 x 512, whichfor a radiation wavelength of λ = 0.9 µ m corresponds to a resolution of 115points per wavelength. The time step was 0.015 fs, which corresponds to 200steps per period of laser radiation. This resolution is sufficient to resolve the28ynamics of the electron-positron plasma. The processes of photon productionand their decay into an electron-positron plasma were modeled with steps splitaccordingly to resolve the corresponding time scales. In this section we employedthe previously discussed methods to resample the particle ensemble using eitherthinning or merging techiniques. All methods except mergeAv were comparedamong each other. The resampling coefficient k = 2 was used in all schemes.In papers [28] and [29], the globalLev method was used. For completeness ofpresentation, this method is also compared with other methods discussed in thispaper.During the interaction the system evolves through several stages thoroughlydiscussed in related papers. First, the target is compressed towards the centerof symmetry, which can be seen as the peak in maximum density at 7–8 T inFig. 10 (a). The next stage, if initial target density is low enough, is a linearQED cascade in a given field. During this stage the maximum pair density andtotal number of particles increase exponentially. This stage is marked in bluecolor in Fig. 10. When pair plasma density becomes comparable with the criticaldensity, the transition to the nonlinear regime occurs, which manifests itself inthe form of saturation of dependence of the total particle number and plasmadensity on time in Fig. 10. The origin of such behaviour is high absorption inoverdense electron-positron pair plasma which leads to a significant drop of thefield amplitude. This region is marked in red in Fig. 10.As was noted previously, such a formulation of the problem allows us to studyboth linear and nonlinear stages of the QED cascade. It should be expected thatthe choice of the resampling algorithm should have the most significant impacton the dynamics of the system in the linear regime, since upon transition to thenonlinear stage the growth rate significantly decreases and resampling of theensemble of particles occurs less frequently.According to the results of this section, the resampling methods can bedivided into several groups, results are summarized in Table 2.The first group constitutes of the simple method. As discussed previously,disadvantages of this method include the formation of macro-particles with large29 igure 10: Temporal evolution of QED cascade in the field of a converging 27 PW dipolewave. (a) The maximum electron-positron plasma density during the development of theQED cascade. (b) The total number of electrons (positrons) in the cylinder with a diameterand height equal to λ . Methods with similar results are grouped accordingly.Figure 11: Spatial distribution of electron-positron plasma density at t = 14 T for differentresampling methods: (a, e) simple ; (b, f) leveling, globalLeveling, conservative ; (c, g) merge,energyThinning ; (d, h) numberThinning . The density of electron-positron plasma is plottedto a logarithmic scale. weights. Figure 10(a) shows the maximum density of the electron-positronplasma during the development of the QED cascade. It is clear that for thismethod the maximum value of pair plasma density exceeds the value obtainedfor all other methods and is exceptionally noisy. At the same time, the totalnumber of particles is almost identical to other methods (see Fig. 10(b)) whichindicates that the problem is the formation of super-particles. It can also beseen that the presence of huge macroparticles leads to an earlier transition to the30 able 2: Comparison of resampling methods for the problem of pinching of electron-positronplasma in a 27 PW dipole wave. Run time, amount of resampling instances, and the averagenumber of particles for the time interval of 100 iterations at linear and nonlinear stages indomain are given for the central process. Simulation of the nonlinear regime for method simple failed due to a numerical instability. Linear regime Nonlinear regimeType Γ T Time,sec × Time,sec × simple globalLev leveling conserve merge energyT numberT simple method is used at the nonlinear stage. Also, the cascade growthrate is underestimated by about 10%. The difference is not as high as for theslow linear cascade in Sec 5.1, but it may still lead to incorrect results. Spatialdistributions of pair plasma at the intermediate stage in comparison with othermethods are shown in Fig. 11 (a, e), the formation of a large number of densesuper-particles is apparent.The second group of methods includes the globalLev , conserve , and leveling .These methods behave almost identically to each other, both at the linear andnonlinear stages of development of the QED cascade. The dependencies of thetotal number of particles are shown in Fig. 10(b). The growth rate yielded bythese methods differs by about 10 − . The plasma distribution in this case is31ymmetric with respect to the azimuthal angle, as can be seen from the Fig. 11(b, f). In this problem due to a high growth rate, unlike the problem in Sec5.1, it is impossible to complete the simulation without resampling. However,taking into account the conclusions made in that section regarding the directcomparison of the growth rate with and without resampling, we can assumethat this group of methods gives the best possible accuracy.The third group includes the methods merge and energyT . These methodsunderestimate the growth rate of the cascade at the linear stage by about 3%,and the number of particles over time withing this group is also identical. Atthe same time, surprisingly, these methods lead to the appearance of spatialmodulation of plasma density, which is absent in the methods considered earlier(see Fig. 11). At the nonlinear stage, these modulations are amplified andstart to significantly affect the dynamics of the system. The nature of thesemodulations requires a separate investigation. It should also be noted that thecharacteristic simulation time of the merge method in this task is almost anorder of magnitude higher than the simulation time of all other methods, seeTable 2. This, combined with other factors, makes it less convenient for practicaluse.The last group of methods consists of the numberT method, which can beconsidered somewhat intermediate. The plasma modulation is present, but itis not as pronounced, and the growth rate is underestimated less than for thethird group (see Fig. 11 (c, g)).
6. Conclusion
The principle of agnostic down-sampling that is applicable without any priorknowledge about the problem was formulated and several resampling methodscomplying with the agnostic principle were presented. Results acquired withuse of these methods were compared among each other, and also to theoreticalresults, results acquired without resampling (where possible) and to the re-sults acquired by some non-agnostic methods. The comparison was performed32rst using simple model problems, and then using pertinent problems involv-ing generation of plasma via QED cascade and thus often requiring extensiveresampling.It was shown that the relative accuracy of various methods highly dependson the problem at hand and the criteria for determining this accuracy. There-fore we conclude that there is no universal method of resampling which wouldshow the best performance in all cases. However, it can be noted that severalmethods provide stable performance on all problems that we have considered.These methods are leveling , globalLev , conserv , and to a lesser extent energyT and numberT . The methods simple and merge/mergeAv have at least one ex-ample where the method in question significantly alters the physical outcomeeven though in certain other conditions these methods might be the most ad-vantageous. Most of the considered methods are released as open-source toolswithin the hi- χ framework [30].
7. Acknowledgements
The work was funded by Russian Foundation for Basic Research and thegovernment of the Nizhny Novgorod region of the Russian Federation, grantNo. 18-47-520001.
References [1] D. L. Bruhwiler, D. A. Dimitrov, J. R. Cary, E. Esarey, W. Leemans,R. E. Giacone, Particle-in-cell simulations of tunneling ionization effectsin plasma-based accelerators, Physics of Plasmas 10 (5) (2003) 2022–2030. doi:10.1063/1.1566027 .URL https://doi.org/10.1063/1.1566027 [2] M. Chen, E. Cormier-Michel, C. Geddes, D. Bruhwiler, L. Yu, E. Esarey,C. Schroeder, W. Leemans, Numerical modeling of laser tunneling ioniza-tion in explicit particle-in-cell codes, Journal of Computational Physics 236332013) 220–228. doi:10.1016/j.jcp.2012.11.029 .URL https://doi.org/10.1016/j.jcp.2012.11.029 [3] E. N. Nerush, I. Y. Kostyukov, A. M. Fedotov, N. B. Narozhny, N. V. Elk-ina, H. Ruhl, Laser field absorption in self-generated electron-positron pairplasma, Phys. Rev. Lett. 106 (2011) 035001. doi:10.1103/PhysRevLett.106.035001 .URL http://link.aps.org/doi/10.1103/PhysRevLett.106.035001 [4] N. V. Elkina, A. M. Fedotov, I. Y. Kostyukov, M. V. Legkov, N. B.Narozhny, E. N. Nerush, H. Ruhl, QED cascades induced by circularlypolarized laser fields, Physical Review Special Topics - Accelerators andBeams 14 (5) (May 2011). doi:10.1103/physrevstab.14.054401 .URL https://doi.org/10.1103/physrevstab.14.054401 [5] I. V. Sokolov, N. M. Naumova, J. A. Nees, Numerical modeling of radiation-dominated and quantum-electrodynamically strong regimes of laser-plasmainteraction, Physics of Plasmas 18 (9) (2011) 093109. doi:10.1063/1.3638138 .URL http://aip.scitation.org/doi/10.1063/1.3638138 [6] C. Ridgers, J. Kirk, R. Duclous, T. Blackburn, C. Brady, K. Ben-nett, T. Arber, A. Bell, Modelling gamma-ray photon emis-sion and pair production in high-intensity lasermatter interac-tions, Journal of Computational Physics 260 (2014) 273 – 285. doi:https://doi.org/10.1016/j.jcp.2013.12.007 .URL [7] A. Gonoskov, S. Bastrakov, E. Efimenko, A. Ilderton, M. Marklund,I. Meyerov, A. Muraviev, A. Sergeev, I. Surmin, E. Wallin, Extendedparticle-in-cell schemes for physics in ultrastrong laser fields: Review anddevelopments, Phys. Rev. E 92 (2) (2015) 023305. doi:10.1103/PhysRevE. .URL http://link.aps.org/doi/10.1103/PhysRevE.92.023305 [8] Y. W. X. Y. Q. B. CHANG Hengxin, XU Zheng, Study on extreme plasmadynamics by quantum electrodynamic particle-in-cell simulations, CHI-NESE JOURNAL OF COMPUTATIONAL PHYSICS 34 (5) (2017) 526.URL http://cjcp.org.cn/EN/abstract/abstract3417.shtml [9] G. Lapenta, J. U. Brackbill, Dynamic and selective control of the number ofparticles in kinetic plasma simulations, Journal of Computational Physics115 (1) (1994) 213–227. doi:10.1006/jcph.1994.1188 .URL https://doi.org/10.1006/jcph.1994.1188 [10] P. T. Luu, T. Tckmantel, A. Pukhov, Voronoi Particle Merging Algorithmfor PIC Codes, Comput. Phys. Commun. 202 (arXiv:1504.00636) (2015)165–174, comments: 11 figures. doi:10.1016/j.cpc.2016.01.009 .URL https://cds.cern.ch/record/2006714 [11] M. Vranic, T. Grismayer, J. L. Martins, R. A. Fonseca, L. O. Silva, Particlemerging algorithm for PIC codes, Computer Physics Communications 191(2015) 65–73. arXiv:1411.2248 , doi:10.1016/j.cpc.2015.01.020 .[12] A. N. Timokhin, Time-dependent pair cascades in magnetospheresof neutron stars i. dynamics of the polar cap cascade with no par-ticle supply from the neutron star surface, Monthly Notices of theRoyal Astronomical Society 408 (4) (2010) 2092–2114. arXiv:http://mnras.oxfordjournals.org/content/408/4/2092.full.pdf+html , doi:10.1111/j.1365-2966.2010.17286.x .URL http://mnras.oxfordjournals.org/content/408/4/2092.abstract [13] A. Gonoskov, Agnostic conservative down-sampling for optimizing statisti-cal representations and pic simulations (2016). arXiv:arXiv:1607.03755 .3514] G. Lapenta, J. Brackbill, Control of the number of particles in fluid andMHD particle in cell methods, Computer Physics Communications 87 (1-2)(1995) 139–154. doi:10.1016/0010-4655(94)00180-a .URL https://doi.org/10.1016/0010-4655(94)00180-a [15] F. Assous, T. P. Dulimbert, J. Segr´e, A new method for coalescing particlesin PIC codes, Journal of Computational Physics 187 (2) (2003) 550–571. doi:10.1016/s0021-9991(03)00124-4 .URL https://doi.org/10.1016/s0021-9991(03)00124-4 [16] D. Welch, T. Genoni, R. Clark, D. Rose, Adaptive particle managementin a particle-in-cell code, Journal of Computational Physics 227 (1) (2007)143–155. doi:10.1016/j.jcp.2007.07.015 .URL https://doi.org/10.1016/j.jcp.2007.07.015 [17] M. Pfeiffer, A. Mirza, C. D. Munz, S. Fasoulas, Two statistical particle splitand merge methods for Particle-in-Cell codes, Computer Physics Commu-nications 191 (2015) 9–24. doi:10.1016/j.cpc.2015.01.010 .[18] D. Faghihi, V. Carey, C. Michoski, R. Hager, S. Janhunen, C.-S. Chang,R. Moser, Moment preserving constrained resampling with applications toparticle-in-cell methods (2017). arXiv:arXiv:1702.05198 .[19] ELI-NP, , accessed: 2020-01-09.[20] ELI beamlines, , accessed: 2020-01-09.[21] CoReLS, https://corels.ibs.re.kr/html/corels_en/ , accessed: 2020-01-09.[22] SULF, http://english.siom.cas.cn/Newsroom/hotnews/201907/t20190710_212831.html , accessed: 2020-01-09.[23] A. Gonoskov, A. Bashinov, I. Gonoskov, C. Harvey, A. Ilderton, A. Kim,M. Marklund, G. Mourou, A. Sergeev, Anomalous Radiative Trapping inLaser Fields of Extreme Intensity, Physical Review Letters 113 (1) (2014)3614801. doi:10.1103/PhysRevLett.113.014801 .URL http://link.aps.org/doi/10.1103/PhysRevLett.113.014801 [24] L. L. Ji, A. Pukhov, I. Y. Kostyukov, B. F. Shen, K. Akli, Radiation-Reaction Trapping of Electrons in Extreme Laser Fields, Physical ReviewLetters 112 (14) (2014) 145003. doi:10.1103/PhysRevLett.112.145003 .URL https://link.aps.org/doi/10.1103/PhysRevLett.112.145003 [25] A. Gonoskov, A. Bashinov, S. Bastrakov, E. Efimenko, A. Ilderton, A. Kim,M. Marklund, I. Meyerov, A. Muraviev, A. Sergeev, Ultrabright GeV pho-ton source via controlled electromagnetic cascades in laser-dipole waves,Physical Review X 7 (4) (oct 2017). arXiv:1610.06404 , doi:10.1103/PhysRevX.7.041003 .[26] M. Tamburini, A. Di Piazza, C. H. Keitel, Laser-pulse-shape controlof seeded QED cascades, Sci. Rep. 7 (1) (2017) 5694. doi:10.1038/s41598-017-05891-z .URL [27] M. Vranic, O. Klimo, G. Korn, S. Weber, Multi-GeV electron-positronbeam generation from laser-electron scattering, Scientific Reports 8 (1)(2018) 4702. doi:10.1038/s41598-018-23126-7 .URL [28] E. S. Efimenko, A. V. Bashinov, S. I. Bastrakov, A. A. Gonoskov, A. A.Muraviev, I. B. Meyerov, A. V. Kim, A. M. Sergeev, Extreme plasma statesin laser-governed vacuum breakdown, Scientific Reports 8 (1) (2018) 2329. doi:10.1038/s41598-018-20745-y .URL [29] E. S. Efimenko, A. V. Bashinov, A. A. Gonoskov, S. I. Bastrakov, A. A.Muraviev, I. B. Meyerov, A. V. Kim, A. M. Sergeev, Laser-driven plasmapinching in e − e + cascade, Phys. Rev. E 99 (2019) 031201. doi:10.1103/ hysRevE.99.031201 .URL https://link.aps.org/doi/10.1103/PhysRevE.99.031201 [30] hi-chi project at Gihub, https://github.com/hi-chi .[31] E. S. Weibel, Spontaneously growing transverse waves in a plasma dueto an anisotropic velocity distribution, Phys. Rev. Lett. 2 (1959) 83–84. doi:10.1103/PhysRevLett.2.83 .URL https://link.aps.org/doi/10.1103/PhysRevLett.2.83 [32] B. D. Fried, Mechanism for instability of transverse plasma waves, ThePhysics of Fluids 2 (3) (1959) 337–337. arXiv:https://aip.scitation.org/doi/pdf/10.1063/1.1705933 , doi:10.1063/1.1705933 .URL https://aip.scitation.org/doi/abs/10.1063/1.1705933 [33] J. G. Kirk, A. R. Bell, I. Arka, Pair production in counter-propagatinglaser beams, Plasma Physics and Controlled Fusion 51 (8) (2009) 085008. doi:10.1088/0741-3335/51/8/085008 .URL https://doi.org/10.1088%2F0741-3335%2F51%2F8%2F085008 [34] V. F. Bashmakov, E. N. Nerush, I. Y. Kostyukov, A. M. Fedotov, N. B.Narozhny, Effect of laser polarization on quantum electrodynamical cas-cading, Physics of Plasmas 21 (1) (2014) 013105. arXiv:https://doi.org/10.1063/1.4861863 , doi:10.1063/1.4861863 .URL https://doi.org/10.1063/1.4861863 [35] T. Grismayer, M. Vranic, J. L. Martins, R. A. Fonseca, L. O. Silva, Seededqed cascades in counterpropagating laser pulses, Phys. Rev. E 95 (2017)023210. doi:10.1103/PhysRevE.95.023210 .URL https://link.aps.org/doi/10.1103/PhysRevE.95.023210 [36] A. Bashinov, P. Kumar, A. Kim, Quantum electrodynamic cascade struc-ture in a standing linearly polarised wave, Quantum Electron. 48 (2018)833–842. doi:https://doi.org/10.1070/QEL16665 .URL http://mr.crossref.org/iPage?doi=10.1070/QEL16665 doi:10.1134/S0021364015150060 .URL https://link.springer.com/article/10.1134/S0021364015150060https://link.springer.com/article/10.1134/S0021364015150060