Benchmarking of 3D space charge codes using direct phase space measurements from photoemission high voltage DC gun

We present a comparison between space charge calculations and direct measurements of the transverse phase space for space charge dominated electron bunches after a high voltage photoemission DC gun followed by an emittance compensation solenoid magnet. The measurements were performed using a double-slit setup for a set of parameters such as charge per bunch and the solenoid current. The data is compared with detailed simulations using 3D space charge codes GPT and Parmela3D with initial particle distributions created from the measured transverse and temporal laser profiles. Beam brightness as a function of beam fraction is calculated for the measured phase space maps and found to approach the theoretical maximum set by the thermal energy and accelerating field at the photocathode.


I. INTRODUCTION
The generation of high-brightness electron beams remains the principle challenge for a number of linear accelerator based projects including the Energy Recovery Linac (ERL). Design of the electron sources relies heavily on the use of space charge simulations [1,2,3,4,5,6]. Still, there remain a number of open questions with regards to understanding and modeling of space charge dominated bunched beams in photoemission guns. Different space charge codes use varying degree of approximations to capture the most significant physics relevant to beam dynamics in photoinjectors. Certain codes allow self-consistent inclusion of complex conducting boundaries at the expense of considerable increase in the required computation time [7,8], while many of the mainstream codes widely used in the design of photoinjectors assume open boundary conditions everywhere except in the vicinity of the photocathode. Due to an inability to use direct self-force calculations in a bunch with ∼ 10 9 particles, artificial smoothing of the space charge potential is employed when tracking macroparticles, which represent the actual bunch, either through meshing of the electron cloud or via introduction of an effective size to the macroparticles in a point-to-point calculation. As a result, and depending on the simulation parameters, the space charge force in simulations can either be overly smooth or grainy as compared to the actual self-force in space charge dominated beams. Additional assumptions are employed by different codes with respect to details of modeling the emission process from the photocathode, whether or not the velocities of individual electrons in the rest frame of the bunch are treated as negligible, 3D * Electronic address: ib38@cornell.edu or 2D or lower dimensionality nature (e.g. uniform cylinders of HOMDYN) of the space charge, etc. The validity of these assumptions must be evaluated for each individual case. Yet, data comparing directly measured beam conditions, detailed phase space distributions in particular, with simulations is sparse for space charge dominated bunched beams such as found in either DC or RF photoinjectors [9,10,11]. Even a relatively simple configuration involves a number of "free" parameters that need to be varied within the uncertainty of the measurement such as the RF phase in an RF gun in order to obtain good agreement between simulations and measurements [11].
We present direct measurements of the transverse phase space distribution using a simple beamline consisting of a DC gun followed by an emittance compensation solenoid. The number of potential variables affecting beam performance is reduced to the bare minimum in such a setup. Careful characterization of the initial conditions such as the laser transverse and temporal profiles and thermal emittance of the photocathode allows us to carry out a direct cross-checking between the measurements and 3D space charge simulations using the codes Parmela3D [1] and GPT [3]. Phase space distributions are then used to calculate beam brightness vs. the included beam fraction, which is compared to the theoretical limit set by the thermal transverse energy and accelerating field at the photocathode.
The paper is organized as following: Section II details the experimental setup and beam diagnostics used in this work as well as our experimental procedures. Section III presents simulation details and data processing procedures used to extract information from the measured data such as second moments and rms emittances. Comparison between data and simulations follows in Section IV. Finally, we conclude with a discussion and outlook for future work.

II. EXPERIMENTAL SETUP
A. Beamline Fig. 1 shows the experimental setup, which consists of a high-voltage DC gun [12] followed by a solenoid used for emittance compensation located at 0.335 m from the photocathode to its center. Fig. 2b shows the distribution of the magnetic field in the solenoid as calculated by POISSON [13] and the actual measured values using a Hall-probe. Both the emittance measurement system (EMS) and an insertable viewscreen are positioned 1.244 m from the photocathode. Additionally, the beamline is equipped with a deflecting RF cavity [14] used to characterize the initial temporal profile of photoemitted electrons, two beam scanners and a Faraday cup all integrated into a data acquisition system for direct phase space measurement. Details on the EMS are in the next subsection. Two different types of materials have been employed for viewscreens: high sensitivity BeO used in temporal measurements with very low bunch charges and the deflecting cavity, and CVD diamond used with average beam currents of up to 100 µA. Each viewport is equipped with 12-bit CCD camera interfaced to the control system. The HV DC gun, initially designed for 750 kV, was operated at 250 kV. The field distribution is shown in Fig. 2a. The gun has reached 420 kV voltage during high voltage processing in the year 2007. Since then, however, we had to limit the gun voltage to a conservative value below 300 kV due to field emission problems. Upon disassembly of the gun we found a considerable amount of dust coming from the ceramic resistive coating, which is believed to be the primary reason behind the strong field emission. Work is underway to eliminate this source of dust from the gun.
The laser system has been detailed elsewhere [15,16]. The laser spot-size was monitored using a 12-bit CCD camera at the location of a virtual cathode.A Pockels cell was used to reduce the 50 MHz train of pulses, with an average power of about one Watt and 520 nm wavelength, to a lower duty factor for beam measurements. A typical average current during phase space measurements was between 10 to 100 µA.

B. Emittance Measurement System
Considerable care is required when designing a system for direct measurement of the phase space of space charge dominated beam [17,18]. The EMS employed in this study is a double-slit system, with non-moving parts; see Fig. 3. The beam motion is achieved with a pair of corrector coils, designed with a vanishing sextupole component in order to provide uniform kick across most of the vacuum beampipe cross-section (see Fig. 4). Each of the two scanners, one prior to the 1st slit, and the other between the two slits consists of a pair of identical coils with opposite direction of excitation current. The coils in each beam scanner have been measured to cancel each other to better than 1%. Thus, each scanner changes only the position of the beam without affecting its divergence. Additionally, a pair of horizontal and vertical  steering coils before each of the slits allows correction for yaw/pitch alignment errors, leading to overall relaxed tolerances for the device. Finally, a small solenoid (about ±1 • rotation angle) is positioned between the two slits to allow for roll compensation, although its use proved unnecessary in practice.
The 1st EMS slit consists of an armor slit with a 200 µm opening followed by a 20 µm vertically selecting precision slit brazed to the water cooled armor slit; see Fig. 3b. With this design most of the beam power is intercepted by the armor slit. ANSYS analysis shows that the EMS system is able to perform without significant deformation of the precision slit (< 10%) with 1 kW of incident beam power.
Phase space measurements can be carried out using either single slit and viewscreen or double-slit and Fara-day cup configurations. Monte-Carlo analysis on scattered radiation using GEANT4 [19] has been carried out for the full EMS system showing excellent signal to noise performance of the system with either of the two configurations [20].
All of the measurements reported in this work have been done using a double-slit method. The maximum scan rate for the beam scanners was 200 Hz. The signal from both slits has been detected using the Faraday cup connected to a low noise current amplifier. A solenoid positioned just before the Faraday cup was used to focus the beamlet to the center of the Faraday cup for increased charge collection efficiency. A typical transverse phase space scan of 100×100 points would take on the order of one minute.
The slit opening size and distance between the two slits have been determined by solving coupled beam envelope equations for a beamlet selected by the 1st slit with inclusion of the space charge force: Here I is the beam peak current after passing through the 1st slit, I 0 = 17 kA is the Alfven current, and (βγ) is the normalized momentum. After a vertically selecting slit with opening d small compared to the beam size, one has σ y = d/ √ 12 and normalized rms emittance ǫ n,y = ǫ n,y0 (d/ (12))/σ y0 , where ǫ n,y0 and σ y0 are the emittance and vertical size of the full beam prior to the slit. By solving Eq. 1 for the beamlet size σ y (L) at the location L of the 2nd slit (or the viewscreen), and comparing σ y (L)/L to the uncorrelated divergence σ ′ y0 at the 1st slit position, one can gauge the effectiveness of the slit system in the presence of space charge. Here we note that the simple requirement that the beamlet after the 1st slit be emittance dominated [17] is necessary but not sufficient. Because the ratio of the space charge and emittance terms in Eq. (1) scales as fast as ∝ σ 3 y , meeting such condition at the location of the 1st slit does not ensure that the beamlet stays emittance dominated all the way through to the location of the 2nd slit. Numerically solving the coupled equations Eq. 1, on the other hand, allows proper characterization of the slit performance. Fig. 5 shows the results of emittance overestimation by the double-slit method in the case of 80 pC bunches with the parameters as indicated in the figure. The separation between the two slits was chosen to be 38 cm leading to the overestimation in emittance measurement to be less than 10% of 0.3 µm normalized rms emittance at 0.5 MeV kinetic energy.

C. Experimental Procedures
Measurements have been taken at 3 different bunch charges: 80 pC, 20 pC and 0.5 pC. The measured laser intensity stability was 2% rms. The laser spot was initially magnified, then passed through a 2.6 mm diameter aperture, which was 1:1 imaged onto the photocathode. The laser pointing stability was 60 µm rms in each transverse direction. Each data set involved taking multiple images of the laser spot on the virtual photocathode, and an image representing the average centroid position was chosen for simulations as detailed in the next section. A typical transverse laser spot profile is shown in Fig. 6.
The laser pulses were temporally stacked using three birefringent crystals [16]. Direct measurement of the initial temporal distribution of electrons was performed with the deflecting cavity with negligible charge per bunch. See Fig. 7. The resolution of the temporal measurement in this case is 1.5 ps rms as limited by the RF to laser synchronization and finite electron beam spot size. To obtain the actual temporal profile of the electron distribution, the data was fitted with 8 Gaussians, then each Gaussian was assigned 1.0 ps sigma corresponding to the value found in an autocorrelation measurement for an unshaped laser pulse [15]. Both the fit to the data and reconstructed temporal profile used in simulations are shown in Fig. 7.
Early on in the measurements, we were able to observe asymmetric transverse phase space distribution, e.g. see Fig. 18. To eliminate possible causes for such an occurrence, careful beam based alignment was carried out before each data collection. A small laser aperture (0.25 mm) placed concentrically with the larger one was used to create a small beam with negligible charge per bunch. The beam centroid vs. the solenoid current data was fit to obtain both the angle and the offset of the magnet's magnetic axis with respect to the beam. The solenoid has been physically adjusted so that its axis coincides with the beam axis to within a few 10's µm and µrad in offset and angle respectively. Similarly, the center of the aperture was imaged to coincide with the electrostatic center of the gun to about 10 µm, ensuring that the central orbit is well aligned throughout the system.

A. Simulation Parameters And Conditions
As the quantum efficiency (QE) of the photocathode in the region of interest was found to have about 10% peakto-peak fluctuations (in absolute terms the QE for GaAs was about 6% in this work), the measured laser transverse profiles and the temporal shape shown in Fig. 7 were used to created 3D distributions for simulations. The thermal emittance of GaAs photocathodes has been investigated previously [21] and for an rms laser spot size σ ⊥ was found to be ǫ n,th = σ ⊥ kT ⊥ /mc 2 with kT ⊥ = 120 ± 8 meV for 520 nm light, mc 2 is electron rest energy.
We have used two 3D space charge codes: Parmela3D which employs a fast Fourier Transform method for solving the Poisson equation on a 3D grid [22] and General Particle Tracer (GPT) with a non-equidistant mesh solver for the space charge force calculation [23]. The same field maps for both the gun and the solenoid magnet were used in both codes. Convergence of the calculation results has been checked for 20k, 100k, and 500k particle distributions and different mesh sizes. Results for the beam envelope are presented in Fig. 8. Fig. 9 shows the results for emittance vs. longitudinal position for the case of 80 pC bunches. An additional difference between Parmela and GPT is that the former reports the relevant beam parameters as a function of time, whereas 3D coordinates of the bunch were projected to a given longitudinal position for the latter. We observe that sufficient convergence is demonstrated with 100k macroparticles for both Parmela3D and GPT. The mesh size was set to 64 × 64 × 64 for Parmela3D and 50 × 50 × 50 for the nonequidistant mesh Poisson equation solver in GPT with the bounding box size set to 5σ in each dimension. Additionally, the Poisson equation solver in GPT provides a choice of 3 different boundary conditions at the bounding box: Dirichlet with zero potential, an open boundary, and an approximate boundary in which the potential at the bounding box is assigned analytically computed values from a uniform elliptical cylinder with rms dimensions set equal to those of the actual bunch. (Zero potential is assigned at the cathode to include image charge  Fig. 10 for emittance vs. longitudinal position. Beam envelope calculations for these boundary conditions produce essentially identical results for these parameters, and therefore are not shown. Open boundary conditions were used for subsequent comparisons with the measurements.

B. Data Processing Procedures
In order to extract second moments from the measured beam profiles and phase space distributions, the data requires appropriate noise subtraction. The general approach follows the general notion of self-consistent unbiased rms emittance analysis (SCUBEEx) [24]: i) a certain contour delineates the signal plus noise region from the noise only region; ii) the average intensity of the outside region represents a noise bias; iii) the noise bias is subtracted from the data while the outside region is assigned 0 intensity; iv) the contour is grown, and the parameter of interest should not change significantly as a function of the included area once all signal is accounted for if uniform random noise is present. We have used two types of contours: circular type for the viewscreen data, and a special boundary detection technique for the measured phase-space distributions. The boundary detection technique is based on the following observation: a binary im- age obtained by applying a threshold to the phase space 2D distribution is likely to form a continuous region for the signal, and many individual islands for the noise. See Fig. 11. The boundary detection algorithm proceeds as follows: i) the data is convolved with a n × n square (image blurring); ii) the smallest threshold is found that generates a single continuous island; iii) n is incremented and step i) is repeated. The process stops when the island starts to include chunks of noise region, which becomes clearly visible. Once the boundary has been found, noise subtraction is verified by growing/shrinking the contour. The contour growing is stopped when less than half of the whole image area becomes available for noise estimation. Fig. 12 illustrates the procedure further. The change in the parameter of interested (e.g. rms emittance) vs. the included area delineating the signal from noise regions represents the uncertainty in the measurement due to the noise subtraction. Emittance or rms values so calculated correspond to 100% of the beam.

IV. COMPARISON OF MEASUREMENTS WITH SIMULATIONS
A. Transverse Phase Space Distributions Fig. 13 shows the comparison of measured transverse phase space at the location of the 1st slit (z = 1.244 m) for 0.5 pC charger per bunch with Parmela3D and GPT simulations. The resolution of the measured transverse phase space is 90 × 90 steps. In the case of simulations, each image is produced using a 300 × 300 2D histogram with additional convolution (blurring) with a 3×3 square. Each image is normalized to the same maximum intensity value. An identical color map to that of Fig. 6 is used throughout. As expected, the calculated rms normalized emittance is in good agreement with the thermal emittance value for this case. Fig. 14 and 15 show the comparison for 20 pC and 80 pC charge per bunch respectively vs. solenoid lens strength. The streak features seen in the measured phase space are due to the motion of the laser spot at the photocathode. Good qualitative agreement can be seen for 20 pC/bunch data, while some discrepancy in the shape of transverse phase space distributions can be seen at larger solenoid current values for 80 pC/bunch. Fig. 16 shows a comparison of vertical rms beam size (a) at the location of the viewscreen (z = 1.244 m) and rms normalized vertical emittance for 100% of the beam (b) as a function of solenoid current for 20 pC bunches. Excellent agreement is seen for the spot size comparison and good overall agreement for the emittance values, although the measured rms emittance appears to be systematically smaller for the 20 pC/bunch case. Fig. 17 shows similar results for 80 pC charge bunches. Different sets of curves for EMS and viewscreen simulations correspond to different laser spots as registered for the two data sets. Good agreement between simulations and measurements of the beam size is seen for 80 pC/bunch before the formation of a beam waist at the location of measurement, while the agreement at larger solenoid currents appears to be less conclusive.

C. Asymmetric Phase Space Distributions
Asymmetric phase space distributions have been measured on multiple occasions for space charge dominated beam conditions in our setup. For example, see Fig. 18. No such asymmetry was observed for low bunch charge running under otherwise identical operation conditions. It was important to understand the origin of this behavior for its subsequent mitigation. We were able to reproduce similar phase space distributions in 3D simulations for laser spots with noticeable asymmetry. To elucidate the mechanism for this tail formation, the particles comprising the tail have been tagged and their portion of the distribution is shown in Fig. 19 for the transverse profile at the location of the photocathode, z = 0 (b), and inside the solenoid, z = 0.35 m (c). The asymmetry in the laser spot causes the space charge forces to push these particles away from the central axis (three is about 2 mm difference between the top and the bottom edges of the transverse distribution in Fig. 19c) so that the particles experience a stronger focusing kick from the solenoid lens. These particles then undergo a cross-over and form the observed phase space tail. Improving the transverse laser shape would reduce the asymmetry in the phase space.

D. Theoretical Limit To Beam Brightness
It is instructive to consider rms emittance as a function of the contributing beam fraction [18]. For example, Fig. 20 shows rms normalized emittance vs. beam fraction as measured and simulated using GPT for 20 pC/bunch for a solenoid current of 3.8 A. Fig. 21 shows corresponding results for 80 pC/bunch for a solenoid current of 3.6 A. Dot-dashed lines show the expected curve for a Gaussian distribution in the phase space with the same rms emittance as 100% of the actual beam. It can be seen that in the case of 80 pC/bunch, the beam strongly deviates from Gaussian distribution, having a substantially brighter core. Core emittance, defined as ǫ n,y,core ≡ dǫ n,y (ξ = 0)/dξ with ǫ n,y (ξ) being the normalized rms emittance as a function of beam fraction 0 ≤ ξ ≤ 1, is given in both figures along with the core fraction ξ core defined as the fraction of the beam with the emittance equal to the core emittance value: ǫ n,y (ξ core ) = ǫ n,y,core .
The beam brightness available from photoinjectors forms through an interplay of several phenomena such as space charge dominated beam dynamics in the presence of time transient and position dependent external fields. The upper limit, however, is set by the thermal emittance of the photocathode and the available accelerating gradient. Consider a short laser pulse illuminating a photocathode placed in the accelerating field E cath . The electron bunch after the emission will assume a pancake shape provided that the laser pulse duration is sufficiently short: σ t ≪ σ ⊥ m/eE cath , m and e are electron mass and charge respectively. This condition is satisfied in most operating photoinjectors with bunched beams. The maximum charge density that can be supported by the electric field is then given by Note the inclusion of the image charge. The average (normalized) beam brightness can be defined as a ratio of average current I avg over its 4D-volume A 4 defined for (x, p x /mc, y, p y /mc) coordinates: This quantity would be related to x-ray brightness for a properly matched undulator, for example. Beam brightness normalized per single bunch is given by B n,avg /f = q/A 4 , with f being the repetition rate, and q the charge contained in the 4D-volume A 4 . E.g. the 4D volume A 4 = dxdp x dydp y /(mc) 2 for a 4D-hypercuboid element with sides dx, dp x , dy, and dp y . The charge contained in A 4 can be written as with σ p being the rms value of the transverse momentum (assumed to be isotropic for both transverse directions), which is σ p = √ mkT ⊥ for Maxwell-Boltzmann distribution of velocities. The dimensionless coefficient κ depends on details of momentum distribution, e.g. κ = 4π corresponds to a uniform circular distribution in p x and p y with a diameter 4σ p , while κ = 2π corresponds to the peak of a 2D Gaussian distribution. Combining Eq. 4 with the charge density as given by Eq. 2, we find the maximum beam brightness normalized per single bunch This result shows the maximum beam brightness available from a photoinjector to be independent from the bunch charge, and is determined by the accelerating field E cath and transverse thermal energy kT ⊥ of the electrons leaving the photocathode. To compare Eq. 5 with the measured data, we compute the beam brightness per single bunch as where is 100% charge, and ǫ n,y (ξ) is emittance vs. beam fraction curve (cf. Fig. 21. Eq. refeq:brmeas additionally assumes an axially symmetric beam with uniform phase space distribution inside an equivalent ellipse with 4πǫ n,y (ξ) area. Fig. 22 shows brightness normalized per single bunch vs. the beam fraction. The theoretical (average) brightness maximum as given by Eq. 5 is shown as well (κ = 4π). Additionally, (B) n,avg /f is computed for a beam that has the same 100% rms emittance as the actual beam but adopts a 2D Gaussian distribution in the phase space. It is seen that an equivalent Gaussian beam does a poor job of describing the 80 pC/bunch beam as the core is substantially brighter for the measured beam approaching the limit given by Eq. 5.

V. DISCUSSION AND OUTLOOK
Benchmarking of 3D space charge codes has been performed with the direct measurements of the transverse phase space for a bunched beam in the spacecharge dominated regime from a DC gun. Overall, good agreement has been found between the measurements and simulations. We observe that an equivalent Gaussian beam assigned the measured 100% rms emittance poorly describes the peak brightness available in the beam at 80 pC bunch charge due to the presence of a substantially brighter core. In particular, for the case of the minimum measured 100% rms normalized emittance ǫ n,y = 1.8 ± 0.2 µm, the core emittance is found to be ǫ n,y,core = 0.31 ± 0.04 µm with 60% beam fraction contained in the core. Additionally, the comparison of the measured beam brightness vs. beam fraction shows that it approaches the maximum theoretical brightness as set by the available accelerating gradient and transverse thermal energy of the photocathode. While the peak brightness of the beam cannot be improved without changing kT ⊥ and E cath , it should be possible to bring a larger portion of the electron beam to approach the brightness limit given by Eq. 5 through proper control of the space charge forces. For example, Fig. 23 shows simulated possible rms normalized emittance (100% of the beam) from the same beamline as used in this experiment for 80 pC bunches as a function of gun voltage. A uniform cylindrical laser distribution has been used in these calculations, with the (same) thermal transverse energy corresponding to GaAs illuminated by 520 nm. The laser pulse duration in these simulations was 12 ps rms, about 50% longer that what was used in the measurements reported here. In addition to employing a longer laser pulse and continuing the work to reach the gun design voltage of 750 kV, noticeable improvements are sought for the transverse laser shape as well as the pointing stability. E.g. based on simulations, the worst 10% emittance beam fraction in our measurements can be mostly attributed to the less than ideal transverse laser profile (see Fig. 6). Finally, we note that despite a significantly more complicated setup and beam dynamics in the full ERL injector [25], where the bunch undergoes acceleration to over 10 MeV, bunch compression and matching into the linac, a much simpler beamline such as the one considered in this study allows exploration of the best beam brightness achievable from the whole photoinjector. In particular, simulations for the full injector [26], where substantially shorter bunches are produced though subsequent drift bunching (∼ 3 ps rms) indicate rms normalized emittances at ∼ 11 MeV which are about 50% lower but otherwise very similar to that shown in Fig. 23. Thus, to continue the work on improving the HV DC gun design and pushing for lower emittances, it is sufficient in many ways to have a simpler setup with beam diagnostics dedicated to such research similar to the one described in this work.

VI. ACKNOWLEDGMENTS
We acknowledge control system support by John Dobbins, John Barley, and Mike Forster, laser support by Dimitre Ouzounov and Heng Li, magnetic measurements of the solenoid by Alexander Temhykh. The following individuals are acknowledged for many useful discussions and their interest in this work: Tsukasa Miyajima, Georg Hoffstaetter, and Dave Rice. This work is supported by the NSF grant PHY-0131508 and NSF/NIH-NIGMS award DMR-0225180.