Grid induced noise and entropy growth in 3d particle-in-cell simulation of high intensity beams

The numerical noise inherent to particle-in-cell simulation of 3D high intensity bunched beams is studied with the TRACEWIN code and compared with an analytical model by Struckmeier. The latter assumes that entropy growth can be related to Markov type stochastic processes due to temperature anisotropy and the artificial collisions caused by using macro-particles and calculating the space charge effect. The resulting noise can lead to growth of the six-dimensional rms emittance and a suitably defined rms entropy. We confirm the dependence of this growth on the bunch temperature anisotropy as predicted by Struckmeier. However,we also find an apparently modified mechanism of noise generation by the non-Liouvillean effect of the Poisson solver grid, which exists in periodic focusing systems even in the complete absence of temperature anisotropy. Our findings are applicable in particular to high current linac simulation, where they can provide guidance to estimate noise effects and help finding an effective balance between the number of simulation particles and the grid resolution.


INTRODUCTION
In modelling of high intensity beams by particle-in-cell (PIC) computer simulation it is of importance to understand whether any observed emittance growth or beam loss is caused by real physical processes, or by non-physical processes like numerical discretisation and grid effects from the Poisson solver, which can be also seen as pseudocollisions with the grid charge distribution. Our emphasis is on the numerical noise generated by the discreteness of the spatial grid used for the Poisson solver and the finite number of particles. We focus on beam parameters typical for high intensity linear accelerators, but our findings can also be extended to circular accelerators with very different ratio of transverse to longitudinal parameters.
It is commonly accepted that real physical processes leading to amplitude and emittance growth can be initial mismatch; mismatch at structure transitions; errors; external or internal resonances; but also the fact that with space charge strictly periodic stationary distribution functions are not explicitly known in 3d (for a detailed discussion see Ref. [2]). In storage rings intra-beam scattering is often observed as another kind of real physical process, which leads to a total emittance growth under certain circumstances. Numerical noise, on the other hand, can have a similar effect on the beam as real intra-beam scattering. It is one of the challenges of high intensity beam simulation to be able to distinguish between physical and numerical growth effects. Towards this end an in-depth understanding and parametrisation of numerical noise is crucial.
Collision or noise effects can, in principle, be associated with entropy growth. The concept of entropy in accelerator beams was first discussed in general terms by Lawson et al. [3]. Later on the Fokker-Planck equation -see Ref. [4] for a detailed discussion of its applicability to intra-beam scattering -was employed to analyse mismatch and turbulent heating of a sheet beam by Bohn [5]. A significant step ahead in this direction has been the noise and entropy growth model by Struckmeier [6,7], which provides a useful theoretical framework based on second order moments of the Vlasov-Fokker-Planck equation. It assumes collisional behaviour and temperature anisotropy to drive 6d rms emittance growth, which is used to explain rms entropy growth.
Evidence of such an emittance growth by grid and particle number related "spurious collisions" in a linac environment -using the PICNIC space charge routine within the PARMILA code -was given in Ref. [8]. The present computational study makes use of the analytical framework by Struckmeier, examines its validity for 3d bunched beams under space charge conditions typical for high intensity proton or ion linear accelerators, and presents different parametric dependencies to serve as basis for optimization. For our study we use the TRACEWIN PIC code [9] developed primarily for linear accelerators. It also employs the PICNIC space charge routine with an rz as well as an xyz option. Detailed comparison with results from other codes or simulation codes more suitable for circular accelerators is left to future work.

RMS APPROACH TO ENTROPY GROWTH
It is worth reminding the fact that in analytical descriptions of high intensity beam phenomena it is generally assumed that the flow of particles in 6d phase space can be described by a Hamiltonian with -in general -time-dependent external forces and a time-dependent self-consistent space charge potential. Liouville's theorem applies, which implies that the volume of a phase space element occupied by particles remains invariant in time. In such a system there is no entropy growth, if infinite phase space resolution is assumed and the Boltzmann entropy as logarithm of probabilities in 6d or µ space is evaluated accordingly.
The same applies to a Hamiltonian flow in the 6Ndimensional phase space (Γ space), which includes sys-tems with Coulomb collisions. Growth of the Gibbs entropy (logarithm of probabilities in Γ space) requires introduction of a "coarse-grained" phase space as opposed to the "fine-grained" infinite resolution phase space, which of course breaks the Hamiltonian nature of the flow. There is a certain analogy between coarse-graining in Γ space and the PIC-code technique of distributing particles on a 3d grid to solve Poisson's equation. The latter also breaks the Hamiltonian nature of the flow in 6d phase space and is a source of entropy growth. Under certain circumstances the interaction with this grid can be also looked at as "collision" with the grid charge distribution.
The idea of a probability based approach to entropy by using the logarithm of the rms emittance was introduced by Lawson et al. [3]. For a time-independent Kapchinskij-Vladimirskij distribution in 4d phase space they thus obtained where k is the Boltzmann constant. An important step towards dynamically evolving distributions has been the approach by Struckmeier [6,7]. He demonstrated that the rms -second order moments -approach originally introduced by Sacherer [10] to derive envelope equations, and later generalized to include the field energy in 6d (see Ref. [11]), can be extended to the full Vlasov-Fokker-Planck equation. The thus obtained equation applies to the product of the three rms emittances, which can be understood as a six-dimensional rms emittance, ǫ 6d ≡ ǫ x ǫ y ǫ z , hence where s = βct measures the distance and k f ≡ β f /βcγ, with β f the dynamical friction coefficient, which is proportional to the Coulomb logarithm. The r nm are rms based "temperature ratios" or ratios of intrinsic spreads of velocities given here in non-relativistic approximation: For upright ellipses we simply have the familiar rms expressions and similar for the remaining ratios. Some caution is necessary when using a concept like "temperature". Strictly speaking, it cannot be defined properly for our beams, which are not in an equilibrium state. In the above model the concept of locally near-isotropic temperatures is formally adopted as basis for the validity of the Einstein relation connecting isotropic diffusion and friction coefficients, which is only valid for not too large deviation from thermal equilibrium.
It is important to note that Eq. 2 suggests a separation of a "friction term" (given by k f ) and a "temperature anisotropy term" I A defined by the bracket on the r.h.s. of Eq. 2. For the present study we shall not be concerned further about a more rigorous justification of k f , which is left to a forthcoming study. Instead, we adopt a phenomenological approach with the goal to demonstrate that this separation is a practical ansatz. The details of the numerical scheme and of the Coulomb logarithm are assumed to enter into the coefficient k f ), which is separated from the "driving" anisotropy term. Note that -as one may expect from extended plasmas -no growth of the rms emittance is predicted from Eq. 2, if all temperatures are identical everywhere and the r.h.s. vanishes. This issue will be examined and discussed further below.
As the r.h.s. of Eq. 2 is either zero or positive, ln ǫ x ǫ y ǫ z cannot decrease. Noting that the probability of decoupled events is the product of individual probabilities and that entropy should be an "extensive" (i.e. additive) quantity, it is justified to relate this logarithm to the entropy in an rms sense analogous to the static KV-case in Ref. [3]: In the remainder of this work growth of ǫ 6d is thus used as a synonym to entropy growth. It is, however, necessary to apply some caution here. An entropy definition based on rms quantities cannot be applied to collective or resonant processes beyond second order, which may also cause growth of rms emittances, but at the same time a decrease of ǫ 6d . This will be discussed in more detail in the following section.

3D BUNCHES IN PERIODIC SOLENOID WITH RZ POISSON SOLVER
We employ the TRACEWIN code [9] with a bunched beam, and in all examples of this section a periodic lattice with thin solenoid lenses and thin rf gaps at the same location. Due to the rotational symmetry the rz Poisson solver option of TRACEWIN/PICNIC is sufficient. It is characterized by the number n cr of radial cells from the origin to R max ; and n cz as radial cells from −Z max to Z max . For these maximum grid extent values we have assumed throughout the values 3σ, where σ is the rms size in each direction. Beyond the grid limits PICNIC determines fields analytically assuming a Gaussian density distribution of the same rms sizes.
Typical matched periodic envelopes over 5 cells -assuming a length of 1 m per cell -for a spherical (breathing) bunch with equal emittances in x, y, z and k 0x,y,z = 60 o are shown in Fig. 1. In all cases -unless mentioned otherwise -we set the option "steps of calculation" to 20/10 in TRACEWIN, which implies 20 lattice calculations and 10 space charge steps per meter of drift space, in addition to one lattice calculation and one space charge step per solenoid. Due to a rounding up algorithm this amounts to 14 space charge steps per lattice cell, which is found as sufficient. As starting distribution function we employ either the TRACEWIN standard option "6d-ellipsoid" (randomly generated in the six-dimensional phase space ellipsoid with parabolic density profiles) or the option "Gaussian" (Parmila type 36).

Distinction of "coherent" and "noise dominated" regimes
As mentioned above the validity of the rms entropy evolution is limited to incoherent noise effects. A distinction between essentially three regimes is recognized in Fig. 2, where the noise effect is enhanced by choosing a particle number as low as 1000. An initial anisotropy with ǫ x,y /ǫ z = 0.35 is assumed for a zero current phase advance k 0x,y,z = 60 o , which results in space charge depressed tunes of k x,y,z = 32/32/38 o . The starting distribution function is the "6d-ellipsoid" option. Results of the emittance evolution are shown in Fig. 2. Three regimes of emittance behaviour are identified: 1. Coherent regime: It is noted that in less than 20 cells a partial exchange occurs between transverse and longitudinal emittances. This is caused by the coherent space charge driven fourth order resonance with a stop-band around k z /k x,y = 1 as shown in the stability chart of Fig. 2 (details see Ref. [12]). The selfconsistently calculated tune foot print (evolving from larger to smaller values of k z /k x,y due to the dynamically changing emittances) indicates the rapid initial changes of the tune ratio over the first few cells. During this resonant phase there is a steep increase of ǫ 6d , followed by pronounced oscillations, and ǫ 6d cannot be interpreted as entropy. This is consistent with the observation that the underlying dynamics takes place beyond second order. Note that the only partial approach to equipartition in this regime is related to the starting tune ratio at the right edge of the stop-band in Fig. 2.

Anisotropy induced noise:
The remaining anisotropy slowly vanishes in the subsequent noise and anisotropy driven regime, which gradually leads to equal emittances and tunes at about 500 m, where k x,y,z ≈ 36 o . Hence the beam becomes fully equipartitioned, and the vanishing r.h.s. of Eq. 2 predicts no further emittance or entropy growth.

Grid induced noise:
The continuing growth of ǫ 6d beyond the point of equipartition is unexpected in the frame of Eq. 2. We assume this growth is caused by a violation of the assumption of a locally isotropic collision type diffusion process -the underlying assumption in the derivation of Eq. 2 -in our periodic focussing system. It must be assumed that this violation is actually a combined effect of several mechanisms: • Coulomb interaction is via charges deposited on the grid rather than direct particle-particle interaction; the distribution of charges on a grid introduces a "coarse-graining" component, which is not subject to Hamiltonian mechanics, hence Liouville doesn't apply and entropy may grow. • The periodic modulation of focusing introduces a "coherent" streaming against the grid -even in a fully isotropic situation. • The modulation of focussing is not slow compared with "collision times".
This suggests that it may be appropriate to introduce a purely grid-related noise term I GN in Eq. 5, which is independent of the temperature anisotropy term: Note that we have introduced a k ⋆ f to take into account that these grid related effects might also have an effect on k f , besides the offset term I GN .
Derivation of such a modified equation from first principles is beyond the scope of this paper. Instead, we take this phenomenological approach and present further support to it in the following sections.

Scaling with intensity
We first study the relative increase in ǫ 6d over 1000 cells for a strictly spherical bunch with Gaussian distribution function, ǫ x,y /ǫ z = 1, k 0x,y,z = 60 o and k x,y,z ≈ 33 o . An example for N = 4000 is shown in Fig. 3 for an rz Poisson solver setting with 16x16 grid cells within 3σ. Note that in TRACEWIN the notion of n cells within mσ implies that n radial or axial grid cells are counted over a radial or axial semi-axis of extent mσ. Beyond this extent -practically in the beam halo region -a different method is applied: fields are determined analytically by assuming a Gaussian charge distribution within mσ. Also note that Fig. 3 indicates a relatively large jitter due to the relatively low number of particles, besides the linear increase of ǫ 6d with distance. The explicit dependence of this grid-induced noise on space charge is shown in Fig. 4, where the linac current is varied. Here, as in all following graphs, we have introduced a relative growth of ǫ 6d relative to its initial value, i.e. ∆ǫ 6d /ǫ 6d . Throughout this study ∆ǫ 6d is normalized to 1000 lattice cells, which makes use of a practically linear evolution of ǫ 6d in s; except for relatively small values of N , where the initial gradient of ǫ 6d was determined and extrapolated to the same value of 1000 cells. We thus obtain from Eq. 6: On the current scale of Fig. 4 the example of Fig. 3 with a space charge depressed tune of k x,y,z ≈ 33 o is equivalent to I=5 mA. The effect of space charge results in a

Scaling with N and grid in rz
While the emittance exchange in the coherent regime of Fig. 2 changes only slightly for a larger number N of macro particles, the rates in the noise dominated regime are large for small N . Noting that the charge per simulation particle is ∝ 1/N , we expect a corresponding decrease for large N .
A more complete picture is shown in Fig. 5, where the relative growth of ǫ 6d , hence the rms entropy, is plotted against 1/N as measure of the charge per macro particle (note the double-logarithmic scale, with the dashed line indicating a strictly linear fit through the origin). For the 16x16 grid a linear dependence of ∆ǫ/ǫ is noted in the range from 1000 up to ≈20.000 particles: For higher N we find that ∆ǫ/ǫ bends off, however. We also show cases with grid resolution of only 12x12 and find that the departure from a linear law occurs even with less particles, and even more for 8x8 grids. The data indicate that for sufficiently large N we enter into a grid resolution limited region, where further increasing of N doesn't help much to reduce growth of ǫ 6d , unless the grid resolution is increased as well. This transition region is characterized by Figure 5: Relative growth of ǫ 6d for N from 1000 to 128.000 macro particles indicating transition to grid resolution limited regions.
a typical average number of particles per cell. With the average number of particles per (toroidal) cell in ∆r∆z given as N/(n 2 c * 3.14/4) we find that the transition occurs typically for 80-100 particles in a toroidal cell (for example for 16.000 particles for n c = 16). This suggests that for significantly less particles statistical fluctuations of charges on the grid become large, hence more particles will help to reduce the noise effect. Using significantly more particles, instead, has no further benefit unless the grid resolution is increased. This connection should be considered if one wishes to optimize the efficiency of a simulation with regard to CPU-time.
In Fig. 6 we show the rms emittance and entropy growth as function of the number of grid cells n c , which is chosen as equal in r and z. As observed above, a resolution above 16 cells has no benefit for only 16000 particles. The 128000 particles case, instead, shows progressive benefit for increasing n c , hence this region is particle number limited. On the other hand, a n c as low as 8 shows only a factor 2 reduction of ǫ 6d by increasing N from 16000 to 128000 due to the lack of grid resolution. In practice, the noise level argument may not be the only one for selecting N . Another important consideration is very small fractional beam loss at the level of 10 −4 to 10 −6 . A reliable statistical level can be reached only with an adequately large N .
In summary the observed grid-induced noise appears to be a combined effect of "pseudo-collisions" of particles with a discretised grid and a statistically fluctuating pop-

TEMPERATURE ANISOTROPY DRIVEN NOISE
Justification of the separation of k ⋆ and a "kinetic" term I A + I GN needs careful examination. It is not a priori clear that I GN is fully independent of, for example, space charge, if we notice that the collision distance relative to the focusing period is not independent of space charge. In the following we study the entropy growth as function of the transverse:longitudinal temperature ratio for given N.

Comparison of numerical anisotropy effects with theoretical predictions
In the first example we assume k 0xyz = 60 o , N=4000, a 16x16 grid with the "6d-ellipsoid" initial distribution. The linac current is chosen again such that k xyz = 33 o for equal emittances. Different emittances in transverse and longitudinal directions are realized in such a way that the product ǫ x ǫ y ǫ z remains invariant. The initial temperature anisotropy on the abscissa of Fig. 7 is determined via Eq. 4 using actual emittances and space charge dependent tunes k from TRACEWIN. In Fig. 7 we show ∆ǫ 6d /ǫ 6d determined from the initial gradient of the 6d emittance from TRACEWIN data (taken over the first 100 or few hundred cells): The simulation results confirm the existence of a minimum emittance growth of ∆ǫ 6d /ǫ 6d ≈ 0.03 at the point, where the temperature ratios r in Eq.3 are all equal to unity. In order to check the applicability of the analytical I A in Eq. 2 we need to fit k f and I GN to the numerical data using Eq.7. The thus resulting "theoretical" curve in Fig. 7 shows good agreement with the simulation data.
We have also tested this comparison for split k 0xy and k 0z to 60/60/47 o -this time using a Gaussian input distribution -and find reasonably good agreement as shown in

Scaling with N
The separation of N −dependence and of anisotropy for the same focusing and initial distribution is further examined in the double-logarithmic graphs of Fig. 9. For large N (above ≈ 10 4 in our examples) we retrieve the typical behaviour of deviation from linear scaling in N −1 already observed in Fig. 5 and interpreted as effect of insufficient lattice resolution in spite of sufficient particle statistics. It is noted that in this linear scaling regime the curves for different anisotropy are basically shifted by constants, which supports that N doesn't enter into I A , but only into k f .

3D BUNCHES WITH XY Z POISSON SOLVER.
The findings obtained for rz−symmetry have also an equivalent in the use of the xyz Poisson solver option. The n cx,y,z are understood here as half number of cells between the maximum grid extent values. We also consider variation of the number of space charge "steps of calculation"keeping the number of lattice calculations fixed at 20/m.

XYZ versus RZ Poisson solver
We first consider the identical periodic solenoidal case used in Figs. 1, 5, 6 and replace the rz Poisson solver by the xyz option. This is shown in Fig. 10 employing (as before) 14 space charge steps per lattice cell. However, we find that for 16.000 particles ∆ǫ ǫ saturates around 0.08mor than 5 times larger than for the rz solver. For 128.000 particles, instead, the xyz-solver results are roughly only a factor of 2 larger than the corresponding values for the rz-solver. For N =16.000 we note that increasing the number of space charge kicks from seven to eleven per cell leads to only half the growth of ǫ 6d for values of n c ≈ 10. Less frequent space charge kicks are "harder", and it is understandable that they enhance entropy growth. This difference, however, is less pronounced for larger N . The observation of "particle statistics limited" and "grid resolution limited" is even more obvious than for the rz Poisson solver: • Grid resolution limited: For small n c there is only weak benefit by increasing N from 16.000 to 128.000 due to insufficient grid resolution • Particle statistics limited: Increasing n c above 6 enhances again the emittance growth for the case N =16.000; in this case the average number of particles per grid cell falls under 10 and statistics appears insufficient. Raising N to 128.000 has visible benefit and no saturation is seen, although the improvement above n c ≈ 8...10 is only minor.
• Optimum grid resolution: There is evidence that for each N an optimum grid resolution exists, where the emittance growth is smallest; the corresponding n c appears to be the higher the larger N .

FODO lattice
We adopt an "equivalent" periodic FODO lattice with rf gaps and cell length again 1 m as shown in Fig. 11. The same phase advance per m is chosen as in the previous examples as well as identical emittances, which results in approximately the same space charge depressed tunes. The alternating focusing causes a strong modulation and local imbalance of "temperatures", which is described in Ref. [6] as source of entropy growth. We expect that this mechanism amplifies respectively adds to the already mentioned "non-Markov" effects found for the periodic solenoid. This is shown in Fig. 12 as function of n c and for different N as well as space charge steps. It is noted that the difference by decreasing the number of space charge steps per cell from 15 (10/m in TRACEWIN) to 11 (3/m in TRACEWIN) is minor. Both cases, low and high particle numbers, show again that an optimum n c exists, which is the higher the larger N . As before, large N is efficient only if the grid resolution is sufficiently large.
Comparing with Fig. 10 for the periodic solenoid we find an enhanced noise level. The results of Fig. 12 also differs from the result reported in Ref. [8], where low and high particle number (transverse) emittance growth is nearly identical around n c = 8. For this grid resolution we find that 16.000 particles lead to practically four times the growth in ǫ 6d than 128.000 particles. Theoretically, following the result of Fig. 5 and Eq. 8, the difference could be even as large as a factor eight. In practice, such high resolution is not feasible -in particular jointly with large N .

DTL lattice
In TRACEWIN -when applied to linac design -it is customary to employ the compressed and faster option of DRIFT TUBE LINAC CELLS. It simplifies space charge calculation to one single step of calculation per DTL cell, or optionally to three steps (marked by 1/cell and 3/cell). We apply this method to the cell shown in Fig. 11 and obtain the results shown in Fig. 13 The results indicate that the frequently used 1/cell option implies a not insignificant additional noise level.

CONCLUSION AND OUTLOOK
We have explored space charge and grid induced noise in 3d PIC simulation with the TRACEWIN code, and compared results with the analytical rms entropy growth equation by Struckmeier based on a collisional approach. Our findings confirm that the 6d rms emittance -defined as product of the three 2d rms emittances -is a suitable measure for the numerical noise and entropy growth generated by space charge and grid effects. It can thus be used as "quality measure" for optimizing simulation of high intensity beams. In particular, we have found that an optimum combination of N and grid resolution exists: further increasing of the grid resolution is useful only, if the simulation particle number (per grid cell) is increased as well. Likewise, a larger number of simulation particles has little benefit unless there is sufficient grid resolution.
We also find that our simulations are in a grid effect dominated regime, which differs from the assumption of a collisional regime assumed in the work by Struckmeier. Thus, we obtain entropy growth even in fully isotropic cases with no temperature differences. Further work is needed to explore in more detail the transition between two distinct regimes: the "grid effect dominated regime" -claimed here -where particles interact in a non-Liouvillean way with the charge distribution on a grid; and a "particle collision regime", where (Markov type) particle-particle collisions are resolved. The latter can be assumed to be increasingly relevant at much higher grid resolution than employed here. In the grid dominated regime case it should also matter that our time step is not sufficiently small compared with a typical transit time across a grid cell. Also, the presence of "coherent flow" in periodic focusing -as opposed to incoherent temperature-like motion typical for steady state plasmas -must in a suitable way enter into a future more detailed analysis.
Further work is required to extend this analysis to 3d bunches with substantially different oscillation frequencies -like slow synchrotron motion -in order to explore PIC simulation noise for circular accelerators in a similar framework.