Periodic Poisson model for beam dynamics simulation

A method is described to solve the Poisson problem for a three dimensional source distribution that is periodic into one direction. Perpendicular to the direction of periodicity a free space (or open) boundary condition is realized. In beam physics, this approach allows us to calculate the space charge field of a continualized charged particle distribution with periodic pattern. The method is based on a particle-mesh approach with equidistant grid and fast convolution with a Green’s function. The periodic approach uses only one period of the source distribution, but a periodic extension of the Green’s function. The approach is numerically efficient and allows the investigation of periodicand pseudoperiodic structures with period lengths that are small compared to the source dimensions, for instance of laser modulated beams or of the evolution of micro bunch structures. Applications for laser modulated beams are given.


I. INTRODUCTION
The fundamental problem of the analysis of beam dynamics is to solve the equation of motion for a multiparticle system (N ∼ 10 3 Á Á Á 10 9 ) in the presence of external forces (caused by magnets, cavities etc.) and self-forces.A particular difficulty is to determine the self-fields of the dynamic many-particle system.The coarsest approximation to solve this problem is the space charge or Poisson approach.It assumes that all particles are in uniform motion with identical velocity vectors.Therefore the problem is equivalent to an electrostatic problem, and in principle one has to add for each observer particle the one over distance contribution from all other sources.This point-to-point interaction would allow the slightly better approach of individual uniform motion per particle, but its effort scales as N 2 .This is avoidable by a continualization so that the point particles are represented by a smooth density function.Additionally the continualization solves a second problem that is usually not in the focus of interest: the near scattering of particles that approach much closer than the averaged particle distance.A proper tracking with near scattering needs not only the point-to-point interaction but also a very fine time step related to the minimal distance between particles.
Continuous density functions and fields are numerically represented by superpositions of basis functions.Often this is done by interpolating field values on a mesh.In this report we use a three dimensional equidistant cartesian mesh for binning and compute the static potential by a fast convolution with a Green's function as described in [1,2].This Green's function fulfills free space (or open) boundary conditions.Also the electric field is calculated on the mesh, by numerical differentiation.The field is interpolated to the particle positions and transformed into a collective force as it appears in the equation of motion.The numerical integration of the equation of motion is not described here.
The Poisson-particle-mesh approach is used widely and successfully in beam dynamics, for instance, in straight sections with constant energy, but also in particle sources or accelerating sections and even in sections with weak curvature.Sometimes the particle distribution has a fine substructure which has to be taken into account, for instance if the particles perform plasma oscillations [3][4][5], or a micro-modulation builds up (microbunch-instability) or they interact in an undulator with a laser.Such high resolution computations need very fine meshes and high particle (or macroparticle) numbers so that they become numerically very expensive.
This numerical effort can be dramatically reduced if the approach assumes and utilizes periodicity: only the particles of one period have to be tracked and meshed.We describe such an approach where the periodicity of the source distribution is interchanged with the periodicity of the Green's function.Except for the calculation of the Green's function, the established method (Poisson-particlemesh with equidistant mesh and fast convolution) can be used completely unchanged.We consider periodicity into one direction.The direction is arbitrary and does not need to coincide with a mesh axis.Perpendicular to this direction free space (or open) boundary conditions are realized.
The memory requirements are reduced due to the lower number of particles and the smaller volume of the problem domain.The reduction factor 1=R is approximately 1=R N ≈ 1=R V , with R N the ratio of the number of particles in the bunch N b to the number of particles per period, and R V the ratio of the volume for the full bunch V b to volume of one period V p .The effort for numerical operations (per step of particle tracking) can be split into the effort for the computation of the Green's function and the effort for all other operations (binning to mesh, convolution, field interpolation, and particle motion).The effort for the first part is increased by the factor E Gp =E G , as the effort E Gp to compute the periodic Green's function is larger than the effort E G for the nonperiodic function.This enhancement is typically 10, as ascertained in Sec.III.Therefore the effort for the first part scales as E Gp =ðE G R V Þ while the effort for the second part decreases with 1=R.The periodic approach is highly effective if the ratio R is large and one is interested in the "pure" periodic behaviour, or in one or few "periodsamples" along the bunch.For a bunch with charge q ¼ 1 nC, peak current Î ¼ 100 A and period length λ p ¼ 1 μm the ratio is approximately R ¼ cq=ðλ p ÎÞ ≈ 3000.
The report is organized as follows.Section II is about general aspects of periodicity.In Sec.III we describe the approach, first the standard method (Poisson, particlemesh, convolution with Green's function), then its modification with periodic Green's function.For a MATLAB implementation of our method, called QField, the numerical effort to calculate the periodic Green's function in a given volume is compared with the effort to determine the nonperiodic function in the same volume.In Sec.IV we demonstrate the method with examples: the FLASH seeding section [6][7][8] with different particle densities and modulations, and parasitic heating after the LCLS laser heater [9,10].For comparison of field calculation methods, we did tracking simulations with Astra [11] and QField.

II. PERIODICITY
Some aspects of spatial periodicity and pseudoperiodicity are illustrated in Fig. 1.The particles in panel (a) are randomly distributed while panel (b) shows a periodically repeated set S p of random particles.Distribution (a) is caused by applying a periodic transport map to random particles that are uniformly distributed into the direction of periodicity.Strictly periodic distributions as in (b) are used to approximate pseudoperiodic distributions as in (a).
For the strictly periodic approach we consider only a set S p of nonrepeated particles.As it can be seen in panels (b-e) there are many sets S p that describe the same periodic distribution.It is easy to manipulate such sets by shifting particles by a multiple of the repetition vector r p .The sets in (c) and (d) are chosen so that all particles are in a volume slice, where the front plane corresponds to the back plane shifted by r p .These slices might have any orientation (or shift), but there are orientations so that the volume of a box (or sphere) that includes all particles is minimal.The distribution in (b) are particles from one volume slice, that have been mapped by a periodic transport map.(Particles with nearly the same spatial coordinates, but with different momenta may be mapped to very different spatial coordinates.)For numerical reasons it might be advantageous to use sets without sudden change of density, as it can be seen in panel (e).This set has been composed similar as (c), but some particles close to the left boundary are shifted by r p (beyond the right boundary) and some particles close to the right boundary are shifted into the opposite direction.The probability of shifting depends on the distance to the boundary.
To keep periodicity, the set S p has to represent a distribution in the six dimensional phase space X ¼ ½x; y; z; p x ; p y ; p z which is periodic with respect to P, the six dimensional periodicity vector.The required property for the transport TðXÞ → Y is TðX þ nPÞ → TðXÞ þ nTðPÞ.This is for instance fulfilled if only one spatial component of P is nonzero and T is a linear function of that component.
Real distributions are always pseudoperiodic since they are finite.Sometimes only a small part of a (long) bunch behaves periodically, as for instance in Figs. 8, 9.The phase space volume covered by the beam is usually small and the transport can be linearized, or at least locally (for one or few periods) linearized.Therefore the Poisson approach can be used for distributions with finite extension in momenta space, although the uniform motion approach assumes identical momentum of all particles.In the same sense the periodical Poisson approach is applicable to distributions with nonzero momentum components of P, although this would mean particles with arbitrarily large momentum.An example is a bunch with energy chirp and micro modulation.
In many cases the effect of field periodicity appears even for charge distributions with very localized periodic pattern, as in Fig. 2. To estimate this, we consider a distribution ρðx; y; zÞ of finite length with localized pseudoperiodic behavior into z-direction.We split this distribution into the averaged, slowly varying part ρ av ðx; y; zÞ −L p =2 ρðx; y; z þ zÞdz and the remaining modulation ρðx; y; zÞ ¼ ρðx; y; zÞ − ρ av ðx; y; yÞ.L p is the period length.We want (a) to show that the contribution of the modulation is localized and (b) to estimate the effect of the slowly varying part in comparison to the periodic modulation.
(a) The total charge per period of ρ is nearly zero, so that each single period of ρ causes essentially a dipole field.
[Due to the pseudoperiodic nature of ρ, the monopole part is not exactly zero, but it is small compared to the field of ρ av ðx; y; zÞ that is estimated in (b).]The longitudinal electric dipole field of one period of ρ declines as Longitudinal field E z on axis of a bunch with round Gaussian transverse profile.The longitudinal profile is rectangular with length L, superimposed by two short periodic bursts.Both bursts contain 5 dirac pulses with equidistant spacing.The bunching factor in the periodic intervals is 0.1, or 10 percent of the charge per period is in the dirac pulses.The period length in the bursts is 0.5σ r =γ and 2σ r =γ.The field is normalized to E n ¼ Λ=ð2πε 0 σ r γÞ with Λ ¼ q=L.(a) Full bunch, (b) detail with burst in the middle, (c) detail with burst close to the edge.
for transverse offsets x, y of the order of the typical transverse size L t .The influence of all periods, that are further away from the point of interest than NL p , scales as The term in square brackets describes the relative effect.It is obvious that distant periods are negligible for L t =ðγNL p Þ ≪ 1.The first example in Sec.IVA demonstrates that about 4 periods in the middle of the distribution are enough to realize quasiperiodic conditions, compare Figs.11 and 12. On the other hand, the pattern about 4 μm from the center is not regular enough.Alternatively the complete laser pulse could be calculated and simulated with a super-period of about 20 laser periods.Such superperiods are also helpful for other broadband processes as the excitation and amplification of microbunching due to shot noise.(b) The field of the slowly varying part is the macroscopic interaction (over lengths that are large compared to the period length) in a bunch without microstructure.For bunches that are long compared to the transverse dimension L t =γ, it can be estimated by the one-dimensional longitudinal space charge field that can be calculated from the longitudinal charge profile ΛðzÞ and the field E δ ðzÞ of a charged disc with the transverse charge profile ηðx; yÞ.For instance the longitudinal field on axis of a round Gaussian distribution with radius σ r is and the transversely averaged field hηðx; yÞE z ðx; y; zÞi is [12,13] with The amplitude of the periodic field can be estimated by ĨZðc 0 2π=L p Þ with Ĩ the amplitude of current modulation and ZðωÞ the longitudinal impedance per length.The impedance per length corresponding to Eq. ( 5) is [12,13]: A macroscopic field with quasiperiodic pattern of microstructures can be seen in Fig. 2. The difference between the curves for finite and infinite bunch length length is the field of the finite distribution without microstructure.The bunch has to be long enough so that macroscopic field is negligible.This also means that in the reverse approach, to calculate periodic behavior by a nonperiodic method, the macroscopic distribution must be sufficiently long, especially if the modulation is weak.

III. METHOD A. Poisson approach and Lorentz force
The electromagnetic field caused by a set S of discrete point particles with charge q i , position r i ðtÞ and velocity v i ðtÞ is estimated based on the assumption of collective uniform motion (of all particles) with the velocity v c .The particle positions are represented by orthonormal coordinates wherein the z-axis points in the direction of motion e c ¼ v c =‖v c ‖ and the directions e ⊥1 and e ⊥2 are perpendicular.After Lorentz transformation to the rest frame, the coordinates are The electrostatic field may be computed from point-to-point interaction as with excluded self-interaction.Therefore S i is the subset of all particles without particle i.The continualization replaces the particle distribution fq i ; r i;c g by a continuous charge density function ρ c ðrÞ and the electrostatic field is the negative gradient of the electrostatic potential Finally we transform back to the moving frame and calculate the Lorentz force

B. Particle-mesh-method
A simple and robust continualization method is to approximate the discrete charges by a piecewise continuous distribution that is constant in the cells of an equidistant mesh.The mesh points r m ðj; k; lÞ ¼ jh x e ⊥1 þ kh y e ⊥2 þ lh z e c are in the middle of the cells and the charges per cell are q c ðj; k; lÞ.For simplicity we have chosen the direction of motion e c to coincide with one mesh axis.This has numerical advantages, as the resolution requirement in longitudinal direction is usually very different from that into transverse directions, but it is not necessary.The choice of the perpendicular directions e ⊥1 , e ⊥2 is free.
The electrostatic potential caused by the charge density ρ ¼ q c ð0; 0; 0Þ=ðh x h y h z Þ in the cell which includes the origin, is with the Green's function This integral can be solved as [2]: . The asymptotic behavior is obvious from Eq. ( 15).To compute the potential at all mesh points ϕ m ðj; k; lÞ, caused by all mesh charges, one has to perform the summation ϕ m ðj; k; lÞ with gðj; k; lÞ the Green's function at the mesh points r m ðj; k; lÞ.This can be done efficiently by a threedimensional fast convolution.The components of the electric field are determined as difference quotients (f.i.E c ¼ ½ϕ m ðj; k; lÞ − ϕ m ðj; k; l þ 1Þ=h z ).These components, that are allocated on shifted meshes (as r m ðj; k; l þ 1=2Þ), are interpolated to the Lorentz transformed particle positions r i;c to find the fields E i;c in the rest frame.
For an efficient computation of the Green's function on the mesh, one utilizes the asymptotic behavior for points far from the origin and calculates only near points by Eq. ( 16).Therefore one calculates first the antiderivative Hð½j − 1=2h x ; ½k − 1=2h y ; ½l − 1=2h z Þ in one octant (j, k, l ≥ 0) of the near-volume and then the Green's function itself.The field in other octants follows from symmetry.Technically, the change to the asymptotic approximation is realized with a smooth switch function as "switched" Green's function The volume of pure asymptotic behavior follows from the condition In Fig. 3 the Green's function is compared with the asymptotic approximation, and the ratios G a =G and G s =G are plotted for cubic mesh cells (h x ¼ h y ¼ h z ) and for mesh cells with large aspect ratio (h x ¼ 10h y ¼ 10h z ).The relative deviation of the asymptotic function is below 0.001 for x= maxðh x ; h y ; h z Þ > 10.Note that the case of infinite aspect ratio (maxðh y ; h z Þ=h x → 0) is very similar to the case h x ¼ 10h y ¼ 10h z in the volume r ≥ 2h x .The switch function is calculated for the parameters C 2 ¼ 2C 1 ¼ 10.The relative deviation of the "switched" function is below 0.001 for the whole volume.

C. Periodic source and periodic Green's function
Each particle of the set S p is infinitely repeated with the spatially periodic shift r p , see Fig. 1(b-c).The continualization of the particles of the set gives the smooth density function ρ c;p ðrÞ and the periodic extension is The particles in S p and the continuous function ρ c;p ðrÞ generate the periodic behavior, but they are not restricted to a certain volume, as the volume of one period.
As for two dimensional charge distributions with infinite size into the third dimension, the potential of periodic distributions is usually not limited to a finite range.Therefore the Poisson integral Eq. ( 10) diverges for all observation points of interest.It is possible to avoid this and to achieve convergence, by normalizing the potential to be zero at a certain point r 0 .The gradient of φc is identical to that of ϕ c .In the following we skip the tilde and write for the potential of the periodic distribution The summation term in the second integral can be interpreted as Green's function of a periodically repeated point charge.The Green's function of one mesh cell follows from the integration This can be expressed by the nonperiodic Green's function as To compute the potential at mesh points ϕ m ðj; k; lÞ, caused by all mesh charges, one has to perform the summation with q c;p ðj; k; lÞ the integrated charge of density ρ c;p in mesh cell ðj; k; lÞ and g p ðj; k; lÞ the periodic Green's function at the mesh points.Therefore the field calculation for periodic source distributions can be done by the same method as for a nonperiodic source, only the Green's function and the source distribution have to be replaced by the periodic extension G p and the nonrepeated representation q c;p .

D. Numerical calculation of periodic Green's function
For the numerical evaluation we set r 0 ¼ 0 and split G p into the terms FIG. 3. The Green's function G, the asymptotic approximation G a and the "switched" approximation G s for switch parameters C 2 ¼ 2C 1 ¼ 10: (a) and (c) comparison of G and its asymptotic approximation G p ; (b) and (d) ratios G a =G and G s =G; (a) and (b) for cubic mesh cells; (c) and (d) for mesh cells with large aspect ratio (h x ¼ 10h y ¼ 10h z ).The Green's function is only calculated for discrete mesh points ( The finite sum is computed directly, using Eq. ( 21).The number M has to be chosen to be large enough to fulfill the condition ‖r AE Mr p ‖ > C 2 maxfh x ; h y ; h z g for the applicability of the asymptotic approximation, compare Eq. ( 22).Especially for meshes with small stepwidth, an early transition to the asymptotic approach is possible.Usually M is not much larger than 2.
In order to calculate the infinite sums, we use the asymptotic Eq. ( 18) and the auxiliary function aðn > 0; p; qÞ In the following we describe two approximations for the auxiliary function.The first, with a Fourier series, converges well for large values of jqj.The second, for small values of jqj, sums the lowest terms of the infinite sum and uses a truncated Taylor expansion for the rest.Approximation 1: The auxiliary function fðp; q; 0Þ is periodic in p and can be written as Fourier series Approximation 2: For small values of q the auxiliary function is calculated with help of a Taylor expansion and the Taylor polynomial with S kþ1;M t ¼ −ψðk; M t Þ=ðk!Þ and ψ the polygamma function [13].For small values of M it is necessary to split the function into a finite sum and a Taylor approximation for a sufficiently high number M t fðp; q; MÞ ≈ To compute fðp; q; MÞ to a desired level of accuracy for all possible arguments, it is necessary to choose a transition value q t [to decide if Eq. ( 34) or (37) is more efficient] and to determine M t and M f , the number of summands in Eq. ( 33) to be taken into account.The values for q t , M t , M f are selected by comparison with "high precession" results f hp ðp; q; 0Þ which have been calculated to numerical accuracy by using a Fourier series to sufficiently high harmonics or by a sufficiently high M t for the Taylor expansion.
For q t ¼ 0.5, M t ¼ 16 and M f ¼ 8 the numerically calculated function fðp; q; 0Þ and the deviation Δfðp; q; 0Þ ¼ jfðp; q; 0Þ − f hp ðp; q; 0Þj are plotted in Fig. 4. The highest values of Δf appear at q ≃ 0.5 where the type of expansion is switched.The accuracy in the range −0.5 ≤ p ∈ ≤ and 0 ≤ q ≤ 1.2 is better then 2 × 10 −12 .
The numerical effort to compute the periodic Green's function is higher than for the nonperiodic function.To estimate this enhancement of effort, we take the CPU times T G , T Gp to compute G, G p on the same equidistant grid as measure for the efforts E G , E Gp .We used a MATLAB implementation with the program version "MATLAB-R2012b" on a processor "Intel(R) Core(TM) i5-200 CPU @ 3.30 GHz" to determine CPU times.Figure 5 shows the enhancement ratio E Gp =E G for meshes with equidistant stepwidth Δx ¼ Δy ¼ γΔz and different mesh volume N x ¼ N y , N z as function of the aspect ratio N z =N x .For aspect ratios greater than 0.3 the enhancement factor is approximately 10 or lower.It should be noted, that the mesh volume for a full-bunch calculation is usually much larger than for one period.The numerical effort for full bunch simulations with high number of particles, for instance the macroparticles coincide with real particles, is usually not driven by the computation of the Green's function, but by the effort for operations as binning to the mesh, convolution, field interpolation and particle motion.This effort scales linear with the number of particles or the mesh volume and can be significantly reduced by a periodic approach.

E. Mesh-periodicity and symmetry
It is obvious from Fig. 1 that the choice of S p is not unique.Nevertheless the convolution of the smoothed density with the periodic Green's function might result in the same values ϕ m ðj; k; lÞ of the potential on the mesh.This happens if the mesh supports the periodicity, which means the vector r p coincides with a mesh vector jh x e ⊥1 þ kh y e ⊥2 þ lh z e c as in Fig. 6(b).
For practical simulations it is desirable to use a mesh that supports periodicity, and to shift the particles of the original set by a multiple of r p so that they fill a one-period-slice as in Fig. 1(b) and (c).It stands to reason, to align the slice with respect to the mesh as in Fig. 1

(b). A small volume of filled mesh cells is preferable.
Meshes that do not support periodicity can be favorable to avoid extreme aspect ratios of mesh cells, as it would happen for supported periodicity if the ratio of nonzero components of r p is extreme.Such a mesh, as in Fig. 6(d 0Þ and absolute error Δfðp; q; 0Þ of the numerical calculation.As f is singular at the origin (p ¼ q ¼ 0), it is truncated for the plot to values below 10.The transition from the Taylor approximation (for small jqj) to the Fourier series is at jqj ¼ 0.5.The Taylor approximation is used for the upper part of the sum from M t ¼ 16 to infinity.The harmonic expansion is truncated after 8 terms.gives reasonable results supposed the generating set S p is sufficiently smooth.(Sudden changes of density have to be resolved by the discretization.)The potential has to be computed in a volume that overlaps partially with the shifted volume.In overlapping regions the potential is not quite identical, due to nonperiodic discretization errors.Therefore an appropriate interpolation between potentials in overlapping regions (while entering one region and leaving the other) helps to avoid discontinuous potentials and infinite fields.

IV. APPLICATIONS
A. FLASH seeding section, with 800 nm energy modulation of a 1.5 kA bunch In the seeding section of the FLASH free electron laser [6][7][8] an electron bunch is energy modulated by an external laser and sent through three chicanes and several undulators as it can be seen in Fig. 7. To investigate microbunching effects the radiator undulator (after chicane 1) and the sFLASH undulators (after chicane 2) have been switched off.The simulation starts directly after the modulator (longitudinal position 162 m) and ends 10 m behind after the transverse deflecting structure TDS (at longitudinal position 204.2 m).
The initial beam (after the modulator) has the following properties: energy E 0 ¼ 585 MeV, charge 0.3 nC, normalized emittance 1.5 μm, peak current 1.5 kA, uncorrelated energy spread σ E ¼ 150 keV and Gaussian longitudinal and transverse profiles.Therefore, the rms length is 80 fs.The beam has been energy modulated in the modulator undulator by a laser with wavelength λ ¼ 800 nm.The duration of the laser pulse of minimal 30 fs (FWHM) is short compared to bunch length, but long enough to permit pseudoperiodic behavior.As the optical beam in the undulator is much wider than the particle beam, the energy modulation amplitude Ê (of 250 keV) is nearly transverse offset independent.
Effects in drift-quadrupole-lattices are the modification of particle energies, the modification of transverse optics and plasma oscillations.For the given parameters, the energy effect is strongest, transverse defocusing is weak and the length for a full plasma oscillation, in longitudinal phase space, is about L p ¼ 120 m [5].Therefore the longitudinal phase space distribution is nearly not altered on the short distance to the first chicane.The longitudinal dispersion of the first chicane R 56 ¼ 220 μm couples relative energy deviations δ E ¼ ðE − E 0 Þ=E 0 to longitudinal displacements δz ¼ R 56 δ E , and creates density modulation due to the initial energy modulation.The effect of the second chicane is weak due to its small dispersion R 56 ¼ 3 μm.More important is the interchange between density and energy modulation, as it is typical for plasma oscillations, on the distance between chicane 1 and 3 of about L p =6.The bunch current with its modulation is shown in Fig 8(a).Up to this position, the energy modulation has grown by more than a factor of ten to an amplitude of about 4 MeV.This amplification is beyond the linear regime and this will be exacerbated in the last chicane with R 56 ¼ 170 μm.This longitudinal dispersion is sufficient to shear particles, that were originally (during modulation) in one period, over several periods as it can be seen within the longitudinal phase space in Fig. 9(a).At rollover "points" in longitudinal phase space, where the band with particles is folded back, current spikes appear as can be seen in Fig. 8(b).These spikes cause strong space charge fields and lead to a complicated energy pattern after a further drift of about 15 m, see Fig. 9(b).
Figure 9 has been calculated by the established space charge code Astra [11] with a three dimensional solver based on the particle-mesh method.For comparison we did nonperiodic and periodic simulations with our implementation, called QField.The spatial resolution of all computations is 20 meshlines per period in longitudinal direction and 5 meshlines per rms-width in transverse direction.The results of full-bunch simulations with one million macro particles with Astra and QField-none-periodic can be seen in Figs. 9 and 10, showing the longitudinal phase space after chicane 3 and 15 m downstream.As both methods use the same approach for field calculation and the same spatial resolution, it is not surprising that the results are in good agreement.
To verify the periodic approach, we did a high resolution full-bunch simulation with 20 million macro particles with QField-none-periodic, and a periodic simulation with same macroparticle density in only one period.The full-bunch simulations was done with longer laser modulation (60 fs, FWHM) to increase the length of the pseudo periodic range.Figs.11 and 5) for an averaged beta-function of 10 m and the section length of about 40 m is 0.35 MeV.In the range of the laser modulation (AE8 μm) the macroscopic wake is even below 0.17 MeV.

B. FLASH seeding section, with 266 nm energy modulation of a 45 A bunch
This example compares the three-dimensional periodic model with a one-dimensional impedance model, for a periodic distribution that is modulated and compressed to achieve current spikes of maximal amplitude.Significant self-effects happen on the drift-quadrupole-lattice after the four-magnet compressor chicane.The one-dimensional impedance model averages the longitudinal field of a round Gaussian beam versus the transverse offset, as derived in [12], see also [5].A periodic version of the one dimensional model was used for the comparison.(The periodic repetition can be realized either by a summation of the spatially shifted field, or by a summation of the harmonics of the impedance in frequency domain.) The setup is the same as before, but with different parameters and only one active chicane: beam energy E 0 ¼ 700 MeV, bunch charge 0.3 nC, normalized emittance 1.5 μm, peak current 45 A, uncorrelated energy spread σ E ¼ 3 keV and Gaussian profiles.The particle beam is energy modulated by a laser of wavelength λ ¼ 267 nm with the modulation amplitude Ê ¼ 0.2 MeV.The dispersion R 56 ≈ 145 μm of the first chicane is adjusted for maximum bunching and maximum current as it can be seen Fig. 13.The second and third chicane are switched off.
Directly before the first chicane, the phase space distribution is almost not changed by self-effects, due to the small beam current, the high particle energy and the short traveled distance to this position.In contrast to the previous example, the particles are tracked through the magnetic fields of the chicane, but the effect of one-respectively three-dimensional self fields are so small, that Fig. 13 looks identical for both cases, and it is not to distinguish from the corresponding figure for a calculation without self-fields.The sharpness of the current spikes is not only determined by uncorrelated energy spread and longitudinal dispersion, but also by transverse emittance and second order coupling to the longitudinal direction.The current spikes are very short and the bunching factor b n ¼ j P expðin2πz ν =λÞj=N is above ten percent even for spectral lines beyond the 30th harmonic.The required resolution of the longitudinal mesh has to be much better than λ=30.It is set to 0.6 nm, which is 450 meshlines per period of the fundamental modulation and thus below the shortest possible rms length of R 56 σ E =E 0 that would appear if the particles are ideally (linear chirp, no nonlinear effects) compressed to minimal length.There are about 220000 electrons per period which are directly simulated.A macroparticle number equal (or higher) than the number of electrons can be easily handled with the periodical approach.
The particle distribution with sharp density spikes is tracked about 23 m through quadrupoles and drifts to the exit of the (inactive) third chicane.Figures 14, 15 show the longitudinal phase space, the bunch current and the bunching for the one-and three-dimensional field model.The current spikes are diverged and the high harmonic content is significantly reduced, in good agreement with both models.It can be seen that the slope of the initial sawtooth-like modulation has been reduced, similar to the reduction of the energy modulation during a plasma oscillation, after a length that is short compared to the period length.In the one-dimensional model all particles in one slice feel the same longitudinal field.Therefore the only mechanism to increase the slice energy spread is individual longitudinal particle motion.The longitudinal field in the spikes is strongest and causes an energy modulation comparable to the initial laser modulation and a large energy spread.But the space charge distribution significantly violates the requirements for the one dimensional approach of a long bunch: γσ z ≫ σ ⊥ .The length of the current spike of few nanometers times γ is not large 13.Application B: longitudinal phase space, current and bunching after 1st chicane.FIG.14. Application B: longitudinal phase space, current and bunching after 3rd chicane, calculated with one dimensional impedance.Chicanes 2 and 3 are off.compared to the typical transverse dimension of about 100 μm (that follows from the emittance ϵ ¼ 1.5 μm=γ and the beta function in Fig. 7).The vectorial field of the three dimensional model depends on the transverse offset and is integrated along betatron trajectories.(For the other extreme of flat bunches γσ z ≪ σ ⊥ the longitudinal field is proportional to the transverse density, as ∂E z =∂z ≈ divE ¼ ρ=ε 0 .)This increases the slice energy spread as it can be seen in Fig. 15.

C. Parasitic heating after LCLS laser heater
At LCLS a laser heater is used to increase the uncorrelated energy spread of the electron beam to counteract a microbunch instability [14].Therefore the electron beam is energy modulated in an undulator that is in the middle of a short chicane, see Fig. 16.The periodic modulation directly after the undulator is (in principle) visible in longitudinal phase space.This is different after the dogleg comprising of the last two chicane magnets: the angular divergence before the dogleg, due to emittance, causes spatial longitudinal divergence due to the dispersion parameter r 52 .This zdivergence is large compared to modulation wavelength so that the periodic behavior is obscured, see ðΔE; zÞ-diagram in Fig. 19.Similarly, in spatial space appears no pattern of the harmonic excitation, see ðx; zÞ-diagram in Fig. 19.Therefore it was unexpected to measure energy spectra with rms widths that are not proportional to amplitudes of the exciting laser field.The position of the spectrometer can be seen in Fig. 16.This effect is a clear indication for self-effects and has been discussed, explained and estimated in [10].
We give a simplified explanation for this effect.The transport equation from the exit of the laser heater undulator to an arbitrary position after the dogleg is with x 0 , x 0 0 , z 0 , δ 0 ¼ δ cosðkz 0 Þ the phase space coordinates after the undulator and x, z the spatial coordinates (at the arbitrary position) after the dogleg and the relative modulation amplitude δ.The energy spread before the undulator is neglected and for the beginning we neglect also the contribution of x 0 .Figure 17 illustrates the transformation by Eq. ( 38): an equidistant orthogonal grid in the ðx 0 0 ; z 0 Þplane is transformed to the ðx; zÞ plane.The chosen numbers are typical for the given setup.It is obvious that the original area has been sheared and modulated and spatial density oscillations have been created.Also particles with the same initial coordinate z 0 but different slope x 0 0 are distributed horizontally, they are still located on lines.To achieve a homogeneous spatial mixing, a horizontal spread of x 0 is necessary that is coupled by r 11 into the horizontal direction again.Supposed the fluctuation of r 11 x 0 is sufficiently large (compared to the horizontal distance between lines with Δz 0 ¼ λ), then the density is uniform, as for instance in Fig. 19.This is different in the vicinity of zero crossings of r 11 as in Fig. 18 with r 11 and optical functions β x , β y .At these positions appear spatial density fluctuations and space charge fields cause an additional energy modulation, which is correlated with the original modulation.
Figures 19, 20, and 21 are results of a numerical simulation with the periodic space charge solver of the setup described in [10] (E ¼ 135 MeV, σ E ¼ 2 keV, peak current 37 A, charge 250 pC, normalized emittance 0.4 μm, modulation wavelength λ ¼ 758 nm).The beam is heated to σ E ¼ 8 keV.The blue particles in the ðΔE; zÞ-and ðx; zÞdiagrams are particles that have been in the same longitudinal period during modulation, the gray particles are a periodic repetition.Figure 19 refers to a position short after the laser heater chicane.The particles behave as expected by a theory neglecting self-effects: the oscillating behavior is not visible, the particles are longitudinally smeared over a length larger than the period λ.The spectral density agrees with the predicted double bump structure [15].heating (from initially 2 keV) to 8 keV as in the previous figures.The parasitic heating effect can be clearly seen.The characteristic of this heating is described in Fig. 23 showing the rms energy spread at the spectrometer entrance versus the initial energy spread.It is remarkable that parasitic heating changes not only the rms width but also the shape of the spectrum, see Figs. 19, 20, and 21 with the energy spectrum in the upper right plot.This calculation has been performed with one million macro particles per period although the number of electrons per period is approximately 550 000.The longitudinal resolution is Δz ¼ λ=50, the transverse resolution is Δx ¼ Δ y ¼ γΔz, so that the Lorentz transformed mesh has the same resolution in all directions.The phase space pictures have been computed separately, with a resolution of 20 meshlines per period.
Figure 24 shows the published measurement [10] together with a periodic and nonperiodic simulation with QField.Although the shapes agree qualitatively, the results are not in good agreement.The first maximum of the measured curve describes a heating from 7.5 keV ("with laser only") to 28 keV.The maximum in the Qfield diagram is a heating from 7 keV ("with laser only") to 12 keV.For details about the curve "theory (laser þ LSC)" in diagram (a) please refer to the original publication [10].

V. CONCLUSION
In this paper we developed a method to calculate the potential and electromagnetic field of a periodic charge distribution that moves in free space with constant velocity.The charge distribution is three dimensional and periodic in one direction of space.The numerical technique is a modification of the well-known particle-mesh-method with fast convolution of the discretized source with the Green's function of one cell.An efficient method has been described to compute the periodic Green's function so that only particles, representing one period, have to be considered for the convolution.
The method has been demonstrated with examples of pure spatial periodicity into the direction of motion.For example (A) a comparison with nonperiodic simulations is shown.Although the quasiperiodic part of the bunch is quite short, a good agreement is achieved.A nonperiodic simulation for example (B) was beyond the available computational capacities.For example (C) some points of the curve in Fig. 24(b) have been verified with large effort by nonperiodic simulations with lower spatial resolution.The CPU time requirement for each point was about 1000 times larger despite the lower resolution.
Not all capabilities of the method have been demonstrated by these examples.For instance a density modulated beam with energy chirp represents a periodicity in geometric-and momentum-space.Due to longitudinal and transverse dispersion the periodicity vector r p changes its length and direction and deviates from the direction of motion.This happens in bunch compressor chicanes.
The periodic approach is numerically much more efficient than nonperiodic simulations of bunches with quasiperiodicity, as the number of particles and the volume of the problem domain is reduced by orders of magnitude.
The generalization of the Poisson approach to a Maxwell approach with coherent synchrotron radiation (CSR) is a topic that has not finally been solved, see f.i.[16].Particularly the conditions for a periodic approach are more restricted, as CSR effects include long range interactions, while the shape of the source distribution might change on that distance.On the other hand, simplified CSR models, that describe the interaction by a time-dependent, offsetindependent longitudinal wake, are widely used and can be considered in a periodic approach.In the linear regime of weak microbunching this is done in with a local steady state impedance [17,18] and a transient impedance [19].

FIG. 1 .
FIG. 1. Periodic particle distributions and periodic shift vector r p : (a) pseudoperiodic random distribution with periodic behavior; (b) generating set S p and periodically repeated particles in blue resp.gray; (c) equivalent generating set in slice volume between two planes shifted by r p ; (d) other equivalent set in different volume slice; (e) equivalent generating set without sharp truncation.
cosð2πnpÞ:ð33ÞThe harmonic sum converges well for large values of jqj, due to the expð−xÞ= ffiffi ffi x p behavior of the modified Bessel function for large arguments.The sum is usually truncated after M f terms, see below.Function values for positive M can be calculated by subtracting the first terms of the summation in Eq. (33) fðp; q; MÞ ¼ f f ðp; qÞ − X

FIG. 4 .
FIG.4.Function fðp; q; 0Þ and absolute error Δfðp; q; 0Þ of the numerical calculation.As f is singular at the origin (p ¼ q ¼ 0), it is truncated for the plot to values below 10.The transition from the Taylor approximation (for small jqj) to the Fourier series is at jqj ¼ 0.5.The Taylor approximation is used for the upper part of the sum from M t ¼ 16 to infinity.The harmonic expansion is truncated after 8 terms.

FIG. 5 .
FIG. 5. CPU time ratio EGp =E G to compute the Green's function with and without periodicity.The functions are calculated on an equidistant mesh with N x × N y × N z points and stepwidth Δx ¼ Δy ¼ γΔz.The periodic shift is r p ¼ N z Δze z .

FIG. 6 .
FIG.6.Equivalent generating sets S p for the same periodic distribution: (a) volume is not limited to a one-period-slice; (b) particles in one-period-slice and mesh that supports periodicity; (c) particles in equivalent one-period-slice; (d) particles without sharp truncation and mesh that does not support periodicity.
FIG. Seeding section of FLASH: optical functions and layout.

FIG
FIG. 19.Application C: longitudinal phase space ðΔE; zÞ, top view ðx; zÞ and energy spectrum for the position 11 m, short after the chicane.