Propagation of numerical noise in particle-in-cell tracking

Particle-in-cell (PIC) is the most used algorithm to perform self-consistent tracking of intense charged particle beams. It is based on depositing macro-particles on a grid, and subsequently solving on it the Poisson equation. It is well known that PIC algorithms occupy intrinsic limitations as they introduce numerical noise. Although not significant for short-term tracking, this becomes important in simulations for circular machines over millions of turns as it may induce artificial diffusion of the beam. In this work, we present a modeling of numerical noise induced by PIC algorithms, and discuss its influence on particle dynamics. The combined effect of particle tracking and noise created by PIC algorithms leads to correlated or decorrelated numerical noise. For decorrelated numerical noise we derive a scaling law for the simulation parameters, allowing an estimate of artificial emittance growth. Lastly, the effect of correlated numerical noise is discussed, and a mitigation strategy is proposed.


I. INTRODUCTION
The necessity of studying space charge effects via selfconsistent simulation codes has increased with the advent of new projects, which require the generation of high intensity beams. The difficult scenario of storing a high intensity bunched beam for a time of the order of seconds raised significant computational challenges. This situation is found for the future SIS100 synchrotron of the FAIR project, and also for the LIU project for the CERN accelerator complex [1,2]. The suspect that artificial effects due to the PIC method would enhance emittance growth has constrained the first simulation studies to frozen models, which rely on noise-free tracking schemes [3,4]. There, the space charge is computed assuming that the beam remains frozen, which allows an analytic description of the space charge force. This approach surely avoids artificial emittance growth. However, this advantage is obtained at the expenses of the self-consistency, which is excluded in simulations with frozen space charge models. Long term simulations of high intensity beams with PIC codes allow a self-consistent treatment, but require an understanding of the origins and the propagation of numerical noise. Hence, simulation parameters have be set efficient-wise in order to mitigate noise induced artifacts and to keep the simulation feasible. The problem of evaluating the artificial effects created by a PIC algorithm while tracking is not new, and studies are made from several points of view. For example in Ref. [5] a Fokker-Planck approach to the beam evolution is taken, and in recent works [6,7] the original concept is extended. However, these approaches gives us no information on the single particles, and non-Markovian regimes (correlations) are neglected. Differently from previous approaches, we discuss here the effect of PIC induced noise from the point of view of the typical beam dynamics code. We start with an analysis of the numerical noise due to the PIC algorithm and its conjoint effect with the integration scheme of a particle tracker. Following this approach, the depen-dence of the rms-emittance growth on the simulation parameters is derived. Further, it is found that the particle's rotation in phase space creates correlations in the numerical noise. Therefore, the generated noise may not purely be white such that anomalous emittance growth appears. This finding is new, and not discussed in previous studies. The paper is organized as follows. In Section II we model the effect of random PIC noise on the dynamics of a single particle, while the effect on the whole ensemble of particles is treated in Sec. III. In Section IV we introduce the concept of correlated and decorrelated numerical noise. We conclude our studies with a summary in Sec. V. The effect of the particle's rotation in phase-space on the numerical noise is discussed in Appendix A, whereas the effect of the periodic random walk on the diffusion of macro-particles is treated in Appendix B.

II. NOISE AFFECTING THE SINGLE PARTICLE DYNAMICS
Before studying the effect of numerical noise on a large ensemble of particles, we focus on the dynamics of a single particle that is exposed to noise. This section is meant to illustrate the modeling of the effect of a fluctuating force field on the dynamics of a single particle by use of the random walk theory. A particle i at the longitudinal position s is described by its transverse phase space coordinates (x i , x ′ i , y i , y ′ i ). Its dynamics in a linear lattice is given by the Hill equations with k x , k y the restoring force terms. In the following only the x-plane is discussed, since the y-plane can be treated equivalently. A particle i with a charge q and mass m 0 of a high intensity beam is exposed to a force due the electric field E x (x i , y i , s) of the beam. The effect of the force can be modeled via space charge kicks [8], where the j-th kick is applied at a longitudinal position s j . The longitudinal distance in between two consecutive space charge kicks ∆s = s j+1 − s j is (for simplicity) kept constant during the tracking. A kick refers to a sudden change in momentum, as with c the speed of light and relativistic β and γ. If the positions of all charged particles are known, the resulting space charge field can be calculated as a superposition of Coulomb fields. We refer to this exact solution as E x,0 (x i , y i , s j ). However, if approximations are used to reduce the computational cost of a space charge simulation, imprecisions are induced, see e.g. Ref. [9][10][11].
The result of the electric field calculation is thus not only fixed by physical parameters, but varies due to the numerical method. The strengths of this imprecision is quantified by the standard deviation of the electric field The effect of electric field fluctuations is studied in an easy mathematical model, in which the standard deviation is randomly subtracted or added the to the exact solution, i.e.
This yields an average of E x,0 (x i , y i , s j ) and a standard deviation of δE x (x i , y i , s j ), which correctly models the fluctuations for a large number of space charge kicks. The fluctuations δE x (x i , y i , s j ) are assumed to be random, i.e. for every longitudinal position s j , where a space charge kick is applied, the fluctuation δE x (x i , y i , s j ) does not depend on the fluctuation at any other longitudinal position s j+k = s j + k · ∆s for any k ∈ Z. In the following, we refer to this kind of fluctuations as decorrelated numerical noise. The set of Equations 1, 2, and 3 describes the single particle tracking with space charge in the presence of decorrelated numerical noise. The evolution of the single particle trajectories is randomly distorted at each position s j , where a space charge kick is applied. After the application of the space charge kick, the single particles are tracked for a distance ∆s to the next phasespace position, where a space charge kick is applied. The evolution of the single particles can be described by a linear transformation of the phase space coordinates, given by a map M , as with the random kick ∆x ′ defined as Here, the positions (x i (s j+1 ), x ′ i (s j+1 )) are the single particle coordinates at s j+1 for a noise-free tracking, i.e. δE x (x i , y i , s j ) = 0, with space charge. Due to the map M , i.e. the transfer matrix, both phase space coordinates x, x ′ are affected by the random kick ∆x ′ . To check the effect of the decorrelated numerical noise on the phase space coordinates, we do a numerical experiment. A test-particle is initialized at (x, x ′ ) = (0, 0) and N M = 1000 macro-particles are distributed according to a matched Gaussian probability density function. The beam is tracked for 1000 turns in a constant focusing channel of length L = 1m. One space charge kick is applied per turn. The tune is set to Q x ≈ Q y ≈ 0.064, while the space charge induced tune shift is set to ∆Q x ≈ ∆Q y ≈ −0.005. The selfconsistent two-dimensional Poisson solver uses a mesh of N G × N G = 64 × 64 grid points. The tracking is repeated n = 10000 times for the same physical conditions, but with different initialization of the particle distribution. Due to the decorrelated numerical noise the particle is found at different final positions. In Figure 1, we plot the distribution of final positions of the single particle after the tracking. The distribution is normalized to the number of particles in the center. The results presented in Fig. 1 show that decorrelated numerical noise equally effects both planes x and β x x ′ , if a large number of space charge kicks is applied. This insight can be retrieved directly from a mathematical model developed from Eq. 4, which is discussed in detail in Appendix A.
To model the effect of decorrelated numerical noise on N M macro-particles at p space charge kick positions, we define two sets of independent variables Z, Z ′ ∈ {−1, +1} NM×p which are Z ij = −1 or Z ij = 1 (resp. Z ′ ij = −1 or Z ′ ij = 1) with equal probability. Using Equation 4, we model the effect of random kicks on the phase space positions (with linear tracking) as The additional random elements Z ij account for the independence of both planes x and x ′ , while the normalization factor 1/ √ 2 can be derived directly from the mathematical model, see Appendix A. The model given by Eq. 6 is our random walk model [13]: At each longitudinal position s j , where a space charge kick is applied on the single particles, we substitute a random kick as described by Eq. 6 to its phase space coordinates.
A. Application: Single particle emittance For an analysis of the effect of the random walk model on the single particle emittance, we have to specify the spacial dependence of the fluctuations δE x (x i , y i , s j ). As a simplified approach, we consider with δE x (0, 0, s j ) = δE x (0, 0) the electric field fluctuation in the center of the beam, which are (for simplicity, only in this first approach) assumed to be equal for any s j . As we will discuss in Sec. III, the spacial dependence, as given by Eq. 7, approximates the fluctuations for round Gaussian beams with standard deviation σ r for the radial coordinate r 2 = x 2 + y 2 . A particle i is initialized at the phase space position (x i,0 , x ′ i,0 ) and tracked following Eq. 1 and 2. At every position where the particle gets a space charge kick, the random walk model is applied. Following the standard theory of random walks, the most probable position of the particle (x i,p , x ′ i,p ) after a series of p random kicks is given by the position for a noiseless tracking simulation, i.e. with δE x (x, y, s) = 0. This can also be seen in Fig. 1, as the maximum of the distribution of final positions is in the center. The variances of the particle's phase space positions are given after p random kicks as: To each particle can be associated a single particle emittance with Twiss parameters α x , β x , γ x . The effect of random kicks on the initial single particle emittance ε SP,xi,0 can be studied by calculating the expectation value of the single particle emittance after a series of p random kicks ε SP,xi,p = ε SP,xi,0 + ∆ε SP,xi,p . Using Equations 8 and 9, we find the average growth rate for noise induced single particle emittance growth to be ∆ε SP,xi,p ∝ pδE x (0, 0).
As intuitively expected, decorrelated random fluctuations in the space charge force field unavoidably lead to a single particle emittance growth. For a detailed analysis of the effect on the whole ensemble of particles, we have to discuss the origins and properties of the numerical noise. These studies are presented in the next sections.

III. NOISE AFFECTING A DISTRIBUTION OF PARTICLES
In this section the numerical noise for the electric self-field of a Gaussian particle distribution is quantified with systematic numerical studies using the random start approach. By applying the random walk model for decorrelated numerical noise, as described in Eq. 6, on each particle of a distribution, the evolution of the rmsemittance for a coasting beam is derived. As an application, the artificial emittance growth for a self-consistent space charge simulation for the SIS100 synchrotron at FAIR [1] is presented.

A. Scaling law for field fluctuations
We consider a Gaussian beam, i.e. the particle's phase space positions are randomly defined according to a Gaussian probability density function (p.d.f.), given as where σ x , σ x ′ , σ y , σ y ′ are the standard deviations. In the following, we use round Gaussian beams with σ 2 r := 2σ 2 x = 2σ 2 y the standard deviation for the radial coordinate r 2 = x 2 + y 2 . The associated radial space charge field at a longitudinal position s j can be derived as where E 0 (s j ) is proportional to the field gradient in the vicinity of x = y = 0. The transverse electric field of a coasting beam, that determines the self-consistent space charge forces, is computed in long-term simulation most efficiently by using a two-dimensional particle-in-cell (PIC) scheme. There, to achieve a high computational efficiency, the Poisson equation in the transverse plane is solved on a finite set of N G × N G grid points and the particle distribution is approximated by N M macroparticles. Since the number of macro-particles N M is much smaller than the number of physical particles N p , an artificial granularity of the distribution is introduced. This is one of the sources of numerical noise.
To study the impact of the simulation parameters N G and N M on the numerical noise, we adopt the random start technique (opposite to the quiet start). In the random start approach the positions of N M macro-particles are randomly initialized according to a p.d.f, as given by Eq. 11 for round Gaussian beams. The resulting electric field is therefore not only determined by the p.d.f. itself, but varies according to the random initialization of macro-particles with an amplitude δE x (x, y, s). The fluctuations δE x (x, y, s) are quantified by the standard deviation of the electric field for multiple random initializations of macro-particles. The dependence of δE x (x, y, s) on the simulation parameters can be retrieved by systematically changing N G , N M and σ r . By fitting the simulation data, we quantify the electric field fluctuations of a round Gaussian beam as Here, is a factor to quantify the standard deviation of the electric field, that is independent on the parameters N G,0 , N M,0 and σ 2 r,0 . It is determined by calculating the standard deviation in the center of the beam, given by δE x (0, 0, s), for multiple random start initializations for σ 2 r,0 the rms beam size, N M,0 the number of macro-particles, and N G,0 × N G,0 the number of grid points. It has thus to be calculated only once for given machine and beam parameters. The factor ξ(x, y, s) describes the effect due to the bi-linear interpolation of the electric field in between grid points. This interpolation causes a grid texture in the standard deviation. As an example, the standard deviation for N G = 16 and N M = 2000 and a Gaussian beam is given in Fig. 2. The statistics is done with 500 random start initializations. The scaling law, Eq. 13, is valid as long as the beam is reasonably resolved, which is the case for N ≥ 16 grid points within [−2σ x , 2σ x ] (and resp. in [−2σ y , 2σ y ]) of the Gaussian particle distribution. The effects of a low mesh resolution are recently studied in Ref. [7]. It is shown that due to grid heating, the emittance growth is artificially enhanced. In this section, we derived the electric field fluctuations δE x (x, y, s) for a two-dimensional beam distribution by using the method of the random start. These results can thus be used for an analysis of the numerical noise in coasting beam simulations. However, the same method can also be used to find the dependence of the fluctuations δE x (x, y, s) on s for a bunched beam in the 2.5D scheme, see Ref. [14].

B. Effect of numerical noise on the rms-emittance
The quantification of the fluctuations of the electric self-field, given by Eq. 13, includes the dependence on x and y. The strength of fluctuations is thus known for any particle i at position x i , y i as δE x (x i , y i , s). Therefore, the effect of decorrelated numerical noise on the full ensemble of particles can be studied. In the following, we derive the evolution of the rms-emittance in the x-plane, defined as As before, only the x-plane is discussed, since the yplane can be treated equivalently. The effect of a random space charge kick on the rms-emittance is described by using random walk model, Eq. 6, developed for the single particle dynamics. The emittance after the application of a space charge kick (with decorrelated numerical noise) at s j and linear tracking over a distance ∆s is given by ε x (s j+1 ). By applying the random walk model, see Eq. 6, on the definition of the rms-emittance, Eq. 15, we derive All fluctuating terms (proportional to ∆x ′ ) are averaged out to zero after many kicks. For the random variables we use (Z ij ) 2 = (Z ′ ij ) 2 = 1 ∀ i, j. Further, it can be easily shown that Z ij Z ′ ij (∆x ′ ) 2 = 0, since Z ij and Z ′ ij are statistically independent. Since (∆x) 2 (∆x ′ ) 2 is much smaller than any other term in Eq. 16, it can be neglected. For a matched beam, with x 2 = β 2 x x ′2 , we therefore find In the following, the electric field fluctuations δE x (x, y, s), as given by Eq. 13, are used in a simplified way, where the fluctuations due to the bi-linear interpolation are set to zero, i.e. ξ(x, y, s) = 1. This is valid for a fine grained mesh, since the effect is averaged out to zero. We derive with For a large number of macro-particles N M the expectation values · can be evaluated by an integration over the p.d.f, which yields Equation 20 gives the averaged artificial emittance growth for decorrelated numerical noise. It depends on the number of macro-particles N M , the number of gridpoints N G , the integration length ∆s and the variance σ 2 r of the distribution.
C. Tune shift due to numerical noise Space charge forces lead to a defocussing of the beam and thus to a change of the betatronic tune Q x . For a Gaussian particle distribution (un-bunched) the associated tune shift is given by where r 0 = e 2 /(4πǫ 0 m 0 c 2 ) is the classical particle radius. Since during a simulation the initial emittance ε x (t = 0) changes while the number of physical particles N p is kept constant, the associated tune shift changes as For the evolution of the emittance ε 2 with the growth rate ∆ε 2 x given by Eq. 20.

D. Application to the SIS100
The scaling law, see Eq. 20, can be applied to predict the emittance growth in a tracking simulation. As an application, the SIS100 heavy ion synchrotron (without non-linear elements) is chosen. The tracking is done with the MICROMAP [15] library. For the PIC solver [11] we used a mesh of N G × N G = 128 × 128 grid points. The current is set to I = 30 mA while the emittances are set to ε x = ε y = 30 mm mrad. The resulting space charge tune shift was ∆Q x = −0.15644 and ∆Q y = −0.15704. The integration step ∆s is chosen to have 12 space charge kicks per betatronic-wavelength. The simulation parameters are chosen in a way to emphasize the numerical noise in order to assure it is the dominant effect causing the emittance growth. The results are presented in Fig. 3, where the red dots show the emittance growth from the simulation. The blue line is the theoretical prediction from Eq. 20, i.e. the random walk model applied to the full beam. No fitting is done to retrieve the blue line.   3: (Color) Emittance growth due to numerical noise for a simulation with NM = 1000 macro-particles (top) and NM = 5000 macro-particles (bottom) for a beam tracked in the SIS100 lattice. The beam parameters are given in the text. The theory (blue line) predicts the artificial emittance growth (red dots). The strong fluctuations in the total emittance ε = ε 2 x + ε 2 y are due to the randomness of the process.
As it can be seen in Fig. 3, the theory developed in this section recovers the artificial emittance growth gener-ated from the simulation of a complex machine. For decorrelated numerical noise the random walk model can thus be used in order to a) understand the effect of noise on single particles and b) make predictions on the evolution of the whole ensemble of particles. Since the scaling law Eq. 20 is derived, this theory can be used to choose reasonable simulation parameters. For example, if in a simulation with N M1 macro-particles an artificial emittance growth of ∆ε x,NM 1 is observed, we predict that for a simulation with N M2 macro-particles The scaling law Eq. 20 can also be used to find a sufficient number of macro-particles for a different number of grid points N G . Until now, the integration length ∆s has been taken fixed. The analysis gets more complex if the integration length ∆s is changed, since this may cause a correlation in the electric field fluctuations. A study dedicated to this effect is presented in the following section.

IV. CORRELATED NUMERICAL NOISE
For decorrelated numerical noise the random walk model works well. However, this description is no longer valid in the presence of noise correlations, since randomness (white noise) is assumed. If a coasting beam is tracked through some circular lattice, depending on the integration length ∆s and the machine tunes Q x and Q y , the particle ensemble will return (close) to its initial positions after a certain number m * of space charge kicks. Correlations in the numerical noise are then created, since each macro-particle i will see the same electric field as m * kicks before, including the same electric field fluctuations δE x (x i , y i , s). Thus, the fluctuation of the electric field will reoccur periodically and will therefore no longer be random. Macroparticles will thus be periodically pushed in positive or negative x-direction (depending on the random initialization of the ensemble of particles). These correlations may lead to an enhanced emittance growth that exceeds the estimates of the random walk model. Further, the correlations may cause a deviation from the Gaussianity of our distribution, by e.g. creating tails. We refer to this deviation from the random walk model as a stochastic resonance [16]. Correlations in the numerical noise are potentially created, if where ∆Φ j gives the intra-kick phase advance in between s j and s j+1 . We define m * ∈ N as the order of the stochastic resonance. For a constant focusing channel, where ∆Φ j = ∆Φ is constant for any j, Eq. 25 can be simplified as The effect of a stochastic resonance of order m * is modeled by an extension of our random walk model. When a stochastic resonance is excited, the numerical noise on the first m * space charge kicks is still decorrelated and thus completely random. The fluctuation of the electric field is now subjected to the correlation property The set Ω R contains all natural numbers up to R, while m * R is the number of kicks, for which the numerical noise itself is not weakening the correlations. The effect of these m * kicks can thus be modeled as before. We define two sets of random variables Z, To model the effect of correlated numerical noise, the random elements Z ij and Z ′ ij are applied repeatedly R-times on all particles. Afterwards, when the distribution has (partly) changed (e.g. due to the artificial emittance growth), the field fluctuations are (partly) alternated and the elements Z ij and Z ′ ij are (partly) randomized. The number of repetitions R until the correlations are washed out depends on the beam-, machine-, and simulation parameters. The product m * R can be taken as a quantifier of the persistence of correlations in the numerical noise. The order of the stochastic resonance m * fixes the number of random elements Z ij and Z ′ ij . Thus, for N k space charge kicks in a simulation, if m * ≥ N k we then retrieve the random walk. If instead m * is a small natural number, correlations in the numerical noise are created. Depending on the number of repetitions R, a small m * may enhance the artificial emittance growth. This issue is discussed in the context of a periodic random walk in Appendix B. The persistence of the correlations depends on m * and on the granularity of the distribution that is controlled by N G , N M and σ r . Therefore, a scaling law of artificial emittance growth in the presence of a stochastic resonance, equivalent to Eq. 20 for decorrelated numerical noise, can only be derived by a dedicated study on the damping of the correlations which is beyond the purpose of this paper. The rate of artificial emittance growth additionally depends on the parity of m * . If m * is an even number, then after m * /2 space charge kicks the signs of all particle's phase space positions in the resonant plane are inverted. Therefore, the whole distribution is flipped, which leads to the same field fluctuation (with opposite sign), i.e.
The effect of a stochastic resonance is thus emphasized if m * is an even number, because additional correlations are created. Since the number of random elements in Z and Z ′ is half as much as without this correlation, a stochastic resonance of an even order m * corresponds effectively to a resonance of order m * /2.

A. Effect of correlations on the emittance growth
To analyze correlations in the electric field fluctuations, it is necessary to check the spectrum of the electric field. A coasting beam is tracked in a constant focusing channel with machine tunes Q x ≈ 0.19 and Q y ≈ 0.24. For simplicity we chose ∆s = L. The electric field E x (x, y, s) and E y (x, y, s) at x = y = 0 is calculated for 1024 turns at a longitudinal position s. In the spectrum, retrieved by a Fourier transform of the data set and shown in Fig. 4, the highest peak occurs at the machine tune positions. Other frequencies are excited due to a coupling of both planes. To study the effect of stochastic resonances in one plane, we have to ensure that we do not have any interference of the other plane. This can be guaranteed by setting the tune of the plane that is not studied close to an integer value. Here we chose Q x ≈ 1.001, which sufficiently decouples both planes, see spectrum in Fig. 5. We systematically scan the vertical tune Q y in the range Q y = 0.1...0.4, while we keep Q x ≈ 1.001 fixed (see discussion above). The initial emittance in the y-plane is chosen for any machine tune such that we have a round beam, i.e. σ x = σ y . For each machine tune we The emittance growth rates, as presented in Fig. 6, are confirming the prediction given by the theory developed in this paper. In the absence of stochastic resonances (decorrelated numerical noise), we find a growth rate that can be predicted by the random walk model. The emittance growth rate has a (local) peak whenever a stochastic resonance is hit. The lower the order of the resonance m * , the stronger we see its effect on the emittance growth rate. If the parity of m * is even, the artificial emittance growth is emphasized, as described above. The emittance growth for m * = 6 is thus very similar to the case of m * = 3.

B. Change of Integration Length
For long term tracking studies it is desirable to avoid the occurrence of correlated noise to a) minimize the artificial emittance growth and b) avoid uncontrollable artifacts, like tails in the distribution. In most simulation studies the machine tune Q y is a multiple of the intra-kick phase advance, i.e.
because the integration length ∆s is a fraction of the machine length L. For a simulation close to a machine resonance, i.e. of type nQ y = N ∈ N for n ∈ N, stochastic resonances due to PIC noise will then be excited simultaneously. This can be avoided by choosing an appropriate integration length ∆s. As a demonstration of the influence of the integration length on the occurrence of stochastic resonances, the study presented in Sec. IV A is repeated for an integration length ∆s = L · 0.95. The results are presented in Fig. 7. The blue line in Fig. 7 is the same as in Fig. 6, where ∆s = L, while the red line corresponds to ∆s = L ·0.95.
The peaks due to the stochastic resonances are shifted for the modified integration length. The effect of changing the integration length can most easily be understood for the simplified case of a constant focusing channel. In fact, if the integration length is set to a fraction f of the initial one, i.e. ∆s = L · f , the intra-kick phase advance is reduced by a factor f . Therefore, the resonance condition, given by Eq. 26, is fulfilled if We can thus remove correlated noise from our working point Q y by changing the integration length, if necessary. We can then guarantee that single particles are solely affected by white noise that can be described by a random walk model, whose effect can be mitigated using the scaling law presented in this paper.

V. CONCLUSION
In this paper we construct a model for predicting the effect of white noise on a particle beam by using an equivalent effective field fluctuation acting on all phase space coordinates. We use simulation results to infer the mathematical modeling of the electric field fluctuations on a single particle in the beam, and use this model to develop a treatment of the full beam ensemble. The inferred scaling law as function of the simulation parameters allows to estimate the artificial emittance growth. We find that for special machine tunes, which recall those of a machine resonance, the PIC noise becomes correlated. The effect of these correlations on emittance growth is more dramatic, and we explain this unusual dynamics in terms of stochastic resonances. Our study shows that the intra-kick phase advance plays a crucial importance on numerical noise as it defines the tunes at which the stochastic resonances appear. This understanding is useful for developing a strategy to avoid or minimize correlations in numerical noise, hence avoiding stochastic resonances. The theory developed in this paper can thus be applied to control unwanted noise effects in tracking simulations of high intensity beams. Lastly, we remark that in our studies we decoupled the PIC induced numerical noise in the transverse planes by setting the horizontal tune Q x close to an integer value. However, in more realistic applications, the effect of coupling can not be neglected, and certainly will introduce corrections to the model developed in this study. The investigations on the coupling between the planes will be part of future studies.
Appendix A: Effect of the intra-kick phase advance on numerical noise In the following, the combined effect of linear tracking and numerical noise is studied. For simplicity, we discuss the constant focusing channel. Let us consider a single particle initialized at position (x i , x ′ i ) that is tracked with space charge through a constant focusing channel. Due to the numerical noise, an imprecision ∆x ′ (x i , x ′ i , s j ) is induced at any longitudinal position s j , where a space charge kick is applied. In total K space charge kicks are applied, i.e. the particles are tracked for a distance d = ∆s · K. The numerical noise at position s j affects the particle's phase space coordinates thus as with a ij resp. b ij the random kick on x resp. x ′ due to numerical noise at s = s j . Since we specialize on a constant focusing channel, we can decompose the transfer and R(∆Φ) = cos(2π∆Φ) sin(2π∆Φ) − sin(2π∆Φ) cos(2π∆Φ) . (A3) Using l = K − j, we simplify Eq. A1 as The random elements a ij , b ij can be derived as (A5) A test-particle initialized at (x i (s 0 ), x ′ i (s 0 )) and tracked through the constant focusing channel accumulates the random elements a ij , b ij . For K space charge kicks it can be described as where (x i (s K ), x ′ i (s K )) are the final coordinates of the test particle for a noise-free tracking. The random elements a ij , b ij , as derived in Eq. A5, depend on the intra-kick phase advance ∆Φ. To study the effect of ∆Φ on the cumulation of numerical noise, the situation of a fourth order stochastic resonance is considered, where the random elements Z ′ ij are correlated as Z ij = Z i(j+4) ∀j ∈ K. (A7) The cumulated effect for K = 1000 is given by a i and b i , as described in Eq. A6. The averaged cumulated numerical noise a i and b i for 100 random configurations of Z ′ ij with a correlation as in Eq. A7 is shown in Fig. 8 and Fig. 9 for intra-kick phase advances ∆Φ in the range ∆Φ ∈ [0, 1].
For most intra-kick phase advances ∆Φ the cumulated numerical noise remains small, since the contributions a ij resp. b ij are averaged out to zero, i.e. a i ≈ 0 and resp. b i ≈ 0. This changes if the tune is close to where here m * = 4, since a fourth order stochastic resonance is considered. In this case, the cumulated numerical noise sums up, i.e. | a i | >> 0 and | b i | >> 0.
For the propagation of numerical noise, these results can be interpreted in the following way. If the condition for a stochastic resonance, see Eq. A8, is not fulfilled, then the elements a ij , b ij are random and thus the random walk model can be applied. Since both parameters are decoupled and are in average equivalent, the model given by Eq. 6 is justified. If, however, a stochastic resonance is excited the model is not valid anymore and the noise has to be treated differently, as discussed in Sec. IV.  To explain the effect of a periodic random walk, we study the following simplified case. A single particle with the one dimensional coordinate x is subjected to N k random kicks of strength ∆x = 1 in arbitrary units. We further enforce a periodicity of m * , i.e. we define a set of random numbers X i ∈ {−1, 1} with i ∈ {1, 2, . . . m * }, while the set of random numbers is applied R times on the particle's coordinate. The number of repetitions R is chosen to fulfill N K = m * R, with N K constant. The dependence of the averaged standard deviation σ 2 x on m * is retrieved in a numerical experiment. The averaged standard deviation is evaluated for 100 random initializations of X i for each m * ∈ {1, 2, . . . 2000}. The results are presented in Fig. 10. We find that for m * < N k the standard deviation has a proportionality as while for m * ≥ N k the standard deviation gets independent on m * since we apply only N k kicks in this example. These results show that for a periodic random walk the diffusion gets enhanced. The smaller the periodicity m * , the larger is the diffusion. In these studies it is assumed that the set of random kicks X i remains unchanged during the N k random kicks. When applying these results to the situation of a stochastic resonance created by the PIC tracking, we have to consider that the correlations in the numerical noise are washed out due to the numerical noise itself. The scaling, as given in Eq. B1, may suggest that the periodic random walk could be an upper bound on the diffusion, because the damping of the correlations reduces R. If the correlations are completely washed out, the diffusion will be equivalent to a random walk. This defines the lower bound on the diffusion, as indicated by the blue line in Fig 10.