Fundamental Photoemission Brightness Limit from Disorder Induced Heating

We determine the limit of the lowest achievable photoemitted electron temperature, and therefore the maximum achievable electron brightness, due to heating just after emission into vacuum, applicable to dense relativistic or nonrelativistic photoelectron beams. This heating is due to poorly screened Coulomb interactions equivalent to disorder induced heating seen in ultracold neutral plasmas. We first show that traditional analytic methods of Coulomb collisions fail for the calculation of this strongly coupled heating. Instead, we employ an N-body tree algorithm to compute the universal scaling of the disorder induced heating in fully contained bunches, and show it to agree well with a simple model utilizing the tabulated correlated energy of one component plasmas. We also present simulations for beams undergoing Coulomb explosion at the photocathode, and demonstrate that both the temperature growth and subsequent cooling must be characterized by correlated effects, as well as correlation-frozen dynamics. In either case, the induced temperature is found to be of several meV for typical photoinjector beam densities, a significant fraction of the intrinsic beam temperature of the coldest semiconductor photocathodes. Thus, we expect disorder induced heating to become a major limiting factor in the next generation of photoemission sources delivering dense bunches and employing ultra-cold photoemitters.


Introduction
Beam brightness is a principle figure of merit for relativistic photoelectron sources for use in high brilliance linear accelerators or for ultrafast electron diffraction (UED) experiments. It is qualitatively defined as the average particle flux per phase space volume. For high brilliance linear accelerators used for x-ray production, the brightness of the x-ray beam is directly determined by that of the electron beam. For UED setups, the brightness of the electron beam is main beam parameter that determines the visibility of diffraction pattern per electron pulse.
Such electron sources are comprised of a photoemitting material placed in an accelerating gradient, where both direct current (DC) and radio frequency (RF) accelerating fields are used for various applications. For either x-ray or electron diffraction experiments, it is often the 4 dimensional transverse normalized brightness that is most pertinent. For a given beam current I, we can define the "micro-brightness" as the phase space density itself, ρ = dI dV 4 , where dV 4 = dxdydp x dp y is the phase space volume element. The normalized total beam brightness can then be defined as a statistical average of the micro-brightness: For a bunched beam with some bunch repetition rate f , the average current can be written as I av = qf , where q is the charge of the bunch. In a previous work [1], it was shown that the maximum achievable beam brightness (either microbrightness or total) can be written: where E acc is the accelerating electric field directly at the photocathode, which sets the maximum supportable charge density at the photocathode. The minimum divergence is set by kT , the temperature of photoemitted electrons, an intrinsic property of the choice of photocathode and laser wavelength. A nonzero temperature arises fundamentally from the electron momentum spread inside the photoemitting material, and can then be significantly increased due to excess laser energy above the photoemission threshold, as well as electron scattering off of imperfections in the emitter.
As it is one of the two independent parameters of the maximum achievable beam brightness per bunch, the photoemission temperature has been the focus of much work in photoelectron sources. Great progress has been made by those working with negative electron affinity semiconductor photocathodes and those pursuing electron emission from laser cooled atoms. Several semiconductor photocathode groups have measured photoemission of equivalent temperatures well below kT = 25 meV, or the thermal energy at room temperature, either via the cryogenic cooling of the photocathode [2], or via the maintenance of a pristine photoemissive conditions in ultra high vacuum, allowing the low effective mass of conduction electrons to produce an effect equivalent to Snell's law in which the velocity spread of electrons is drastically reduced at the vacuum interface [3]. Furthermore, using laser cooled atoms which are photoionized, there has been production of electrons at temperatures equivalent to ∼ 1 meV [4].
The number of electrons per bunch varies widely across various applications, with a range approximately between 10 6 → 10 9 electrons, which for moderate to high flux applications often corresponds to densities in the range of n 0 = 10 17 → 10 20 m −3 . For a near-zero temperature bunch with such high density, we expect some contribution of individual stochastic Coulomb interactions from close encounters just after emission into vacuum to add to the total effective photoemission temperature. In this work, we determine this amount of stochastic heating as a function of beam density and initial temperature, as well as the nature of its evolution in time.
We expect this effect to be most prevalent when the electrostatic potential energy of neighboring particles is comparable to their thermal energy, that is kT ∼ e 2 /4π 0 a, where a is the Wigner-Seitz radius, a = (3/4πn 0 ) 1/3 . Thus for a given density, we expect the heating to be of the order e 2 /4π 0 a, and should thus scale with the cubic root of the density, and should be independent of the number of particles in the bunch. For a rough estimate of the importance of the effect, the plasma coupling parameter Γ can be used, defined as the ratio of kT to the average pair interaction potential. It ranges from Γ = e 2 /4π 0 akT = 0.2 → 2 given an electron temperature of 5 meV. Thus, for applications requiring a large charge density, we expect this stochastic heating to serve as a hard limit to the lowest attainable electron temperature, and limiting the maximum attainable beam brightness.

Failure of Analytic models for Coulomb Collisions
We will now describe how traditional collisional methods in beam physics, familiar to many accelerator practitioners, fail for the case of a cold dense beam. Readers familiar with the inability of such methods to accurately describe our strongly correlated case can bypass this section. A simple model for a bunch that has just been photoemitted into vacuum is a static, uniform, randomly distributed electron sphere with very small initial temperature kT ∼ 0, in a constant accelerating field. We first draw a sharp distinction between the stochastic heating in question and the effects of space charge, which is the collective, mean field effect of Coulomb repulsion. The space charge approximation, applied in most beam physics calculations, self-consistently calculates the interparticle interaction based on the local single particle beam density, ∇ 2 Φ(r) = en(r)/ 0 , where Φ is the total electrostatic potential of a particle at r. This approximation requires that the individual electron interaction is heavily screened, meaning that the Debye screening length λ = 0 kT /n 0 e 2 is much larger than the interparticle separation. However, for the lowest temperatures attained in semiconductor photoemission kT ∼ 5 meV [3], and the lowest of the above densities n 0 = 10 17 m −3 , the Debye length is already on the order of the interparticle separation. In this case, the collective field will describe the overall density evolution, but cannot capture the growth of the stochastic component of the velocity.
The effect of Coulomb interactions in particle beams has been treated analytically in various schemes. Perhaps the most famous is the diffusive Fokker-Planck method. The Fokker-Planck method assumes that the effects of Coulomb collisions can be treated via the calculation of effective velocity diffusion and dynamical friction terms. This approach requires that the shifts in velocity due to Coulomb collisions are small compared to the overall velocity spread of the beam. This assumption is maximally violated for the case of a cold dense beam as described above, in which transverse velocities begin near zero.
Non diffusive, two-particle methods have also been developed that do not require the assumption of Debye screening. A clear presentation of these is given in [5]. The method applicable over the largest parameter space is the Extended Two-Particle Approximation (ETPA), developed by Jansen. This method relies on the calculation of the displacement (either velocity displacement, or position displacement) of a test particle in the presence of a single field particle over some collision time δt c . If we calculate the velocity displacement ∆v of a particle pair initially at rest, the ETPA allows the formation of the probability distribution ρ(∆v), by averaging ∆v over all possible separations and initial velocities in the beam. Each of these encounters is assumed to be statistically and dynamically independent. The second moment of the velocity distribution ∆v 2 is then a measure of the temperature. We will proceed with a sketch of this calculation to highlight the failure of some of its assumptions in the cold dense beam case.
For two particles initially at rest (i.e. kT ∼ 0) with initial separation r i , using the dimensionless variablesr = r/r i andt = t er 2π 0 m the equation of motion for their separation is given by: This equation is integrable for the functiont(r), which is not analytically invertible, but is trivial to invert numerically to obtainr(t). With a global choice of δt c , and with a test particle chosen at the origin, we can obtain ∆v as a function of r i by replacing the scaling factors. Then, we average over the entire distribution of r i : wheret is also evaluated at δt c and r i . The velocity kick ∆v falls off sufficiently fast with large separation, and the integral measure r 2 dr ensures that the averaging does not diverge at small r, and thus we can integrate over all space. The equation 4 is the statistical average of all two particle interactions in a beam for during some time δt c .

Unbound Trajectories
The expression given in 4 and similar uncorrelated two-particle collision methods fail in the cold dense beam case for two reasons. First, it assumes that each pair-wise interaction is statistically independent from each other. There can be no dynamic correlation between separate pair interactions. However in the cold dense beam case a large contribution to the final temperature can be given by simultaneous 3-body (or higher) interactions. The assumption of statistical and dynamical independence allows for the unbound expansion particles with very small initial separation, which are the most pertinent interactions in the cold dense beam case. Even from the definition of the scaled coordinates, it is clear that any choice of δt c corresponds to a some electron separation r c below which all collisions taking place in that time will be sufficiently complete collisions, in which all potential energy is converted to kinetic energy. Alternatively put, given δt c there will always be some electron seperation small enough to make τ arbitrarily large. If we assume the initial distribution of electron separations to be uniform over all length scales, there is no unambiguous choice of cutoff in for r c to avoid such unphysical free expansion, and thus no inherent timescale for two particle collisions across the entire bunch. For perspective, in other Coulomb collision calculations, this ambiguity is seen the calculation of the socalled "Coulomb Logarithm", defined as ln (b max /b min ), or the logarithm of the ratio of the maximum to minimum impact parameters over the whole bunch. However, as it appears under the logarithm in such problems, the minimum distance is often not considered to be a sensitive parameter [6].
In a similar method to the ETPA, in [7] Massey et al. calculate the uncorrelated root mean square fluctuation in the interaction force, and from it they obtain a stochastic energy spread. The authors directly impose a minimum interaction distance which corresponds to those collisions for which half of the potential energy is converted to kinertic energy over a certain time. Collisions with more of a fraction of potential energy release (i.e. smaller sepration) are ignored. They readily acknowledge the ambiguity of this choice, and argue that the effect is small for the beams in their study. However, for a cold beam just after photoemission, it is this fast-timescale release of potential energy as close neighbors rearrange that we are interested in calculating. Thus, an average over independent two particle interactions is not sufficient here.

Scaling with Density
Furthermore, even if one does make a choice of δt c based on some other reasoning, the fact that this scheme involves taking a moment of the single particle distribution means one will always find v 2 ∼ n 0 , whereas the effect we desire to calculate should have v 2 ∼ n 1/3 0 , argued above. This scaling of the velocity spread with density occurs in the Fokker-Planck method (as shown in [6] in Eq. 5.243), in the ETPA (as shown in [5] in Eq. 5.5.11), and in the work by Massey (reference [7] in Eq. 19).
A method that produces the correct scaling with density, but is also flawed, is again given by Jansen in [5]. In what he calls the "thermodynamic limit", Jansen takes the difference of of an the total electrostatic potential enery of an initially uniform distribution of charge with no screening, and a final state of a Debye screened distribution with some kT , and sets this difference equal to the the total thermal energy, 3 2 N kT . However, instead of using a single particle density, Jansen implicitly uses the two-particle correlated density for a Debye screening: Here U i and U f stand for the initial and final potential energy in the system, respectively given by the first and second terms of the integral. φ(r) is the interaction potential of two electrons seperated by r, and is assumed to have the form: φ(r) = e 2 exp(−r/λ)/4π 0 r. The factor n exp(−φ(r)/kT ) is the final two particle correlated density. This equation has a solution of the form: Where α is a dimensionless number determined by numerical solution of the above, α = 0.08702 ‡ This number corresponds to an coupling factor at thermodynamic equilibrium of Γ eq = 3.53. However, if one evaluates the number of particles in the Debye sphere, one finds N D = 4 3 πλ 3 n 0 ≈ 0.03, whereas the Debye approximation requires that the number of particles in the Debye sphere must be large. Thus, we may not apply the Debye/Yukawa form for φ(r) nor for the two particle density function used in 5. It is the use of a two particle correlated density function that provides the correct scaling with density, however, it is difficult to analytically compute the correlation in general. It must be found by some other numerical means [8].

Disorder Induced Heating
The heating associated with the relaxation of a random, near-zero temperature distribution of charges is well known to the ultracold neutral plasma (UNP) community. In such systems, a cold gas is laser ionized, and after a time on the order of the τ = 2πω −1 p = 2π (n 0 e 2 /m e 0 ) −1/2 . In the traditional plasma physics terminology, this effect is referred to as disorder induced heating (DIH), and the effect is seen for Γ of order unity or larger [9], and we will argue that an exact analog of this effect is present in practical photoemisison. Disorder induced heating has been experimentally observed for the ions in a neutral plasma [10]. In cold neutral plasmas, the electrons equilibrate much faster than the ionic ω −1 p , and then serve to screen ion-ion interaction. An expression for the amount of DIH in ions was first given for Yukawa systems (for which electron screening, but not electron-ion recombination, is considered) in [11]. The initial ionic distribution is uncorrelated, and can be defined as the zero energy state. As the ions relax, the ion distribution begins to develop order, and the resulting correlations have an associated ‡ An incorrect value of alpha is quoted for α in [5], which explains the use of the Debye relations, though the final value of λ shows that they are not applicable. binding energy. This correlation binding energy can be calculated from the two particle density function g(r), and is only a function of Γ and the electron screening parameter κ = λ/a: U c N kT ≡Ū = 1 2 This binding energy is balanced by an increasing ion temperature. Calculating the increase in temperature thus only requires knowledge ofŪ for a given electron screening. Analytic calculation of g(r) for strongly coupled systems is difficult, and alternatively, U has been tabulated via molecular dynamics (MD) simulations [11]. It has been demonstrated that a trapped, charged plasma, such as an electron bunch, is equivalent to the one-component plasma (OCP) model, in which the charges exist in a uniform neutralizing background [12]. The external containing potential allows the initial uncorrelated state to be viewed as zero total energy. As electrons relax, correlations develop, and the presence of a "Coulomb hole" in g(r) for r < a creates an effective correlation binding energy equal to that of a one component plasma. In the case of an OCP, there is no second species to provide additional screening, and thus U is only a function of Γ. Owing to this simplification, MD data has been fit to a power series relation forŪ = aΓ + bΓ 1/3 + c, where the coefficients are given in [12], a = −0.90, b = 0.590673, c = −0.26569. We may now write down the expression for the the final coupling Γ f , and thus the final equilibrium temperature, for a fully confined, initially uncorrelated (Ū i = 0), uniform distribution of charges with initial coupling Γ i , analogous to that presented in [11]: It is important to note the power series relation is quoted only for Γ > 1, however we find that for Γ < 1, the above expression is well approximated by Γ i = Γ f , or the limit of no DIH, as we expect. The final coupling given in 8 is plotted in figure 1. For an ideal bunch with zero initial temperature, Γ i → ∞, 8 can be solved to give Γ f ≈ 2.23.
Practical relativistic photoelectron sources contain a number of complicating factors. Inside the emitting material (either a cold gas or crystal), the electron coupling, and thus DIH, will be dramatically reduced by the presence of the ionic/nuclear potential, and thus here we only need to consider DIH developing in vacuum. The beam density during emission will vary in time and with position across the bunch. The temperature associated with DIH should be reached on the order of ω −1 p after electron emission into vacuum, which for typical beam densities is on the picosecond scale. On this scale, we can neglect relativistic effects, such as the 1/γ 2 damping of the interaction force [6]. The coupling given by Γ f = 2.23 ≡ Γ eq , for a given local beam density, can then be viewed as the fundamental photoemission temperature limit. We consider the effects of a changing local beam density, as well as the nonrelativistic effects of acceleration during the evolution of DIH below, and we find it to correspond well with the prediction of 8. Rewriting it in terms of practical units, the heating induced in a zero temperature bunch is given by kT [eV] = 1.04 × 10 −9 (n 0 [m −3 ]) 1/3 . Along with 2, this forms the fundamental phoemission brightness limit.
To test the prediction of 8, we have chosen to use a tree-algorithm electrostatic simulation package, with nonrelativistic, spherical, uniform particle bunches. Traditional particle-in-cell integrators have an intrinsic length scale-the spatial grid on which local fields are calculated. However, tree algorithms lack a spatial grid, and thus are better suited to our initial conditions, where both close interactions and long range correlations are significant. Several plasma treecodes have been developed, see for instance [13], and one has been included in a relativistic particle accelerator code [14]. However, considering the simplicity of the problem, we elected to directly modify a code originally written for gravitational interaction [15].
We consider a randomly distributed spherical bunch of 10 5 particles with uniform density, and vary only the density (10 17 → 10 20 m −3 ) and initial Gaussian velocity distribution (kT ∈ {0.25, 20} meV or kT = 0) , where we use the definition kT = where m e is the electron mass, and i is a Cartesian coordinate. There is no breaking of spherical symmetry in our simulations, thus equipartition will always hold. For a fully contained uniform bunch, x i v i ≈ 0, and the temperature is just a measure of average electron kinetic energy, with no spatial dependence. Containment is done in simulation via the application of a radial external force equal and opposite to the SC field. The parameters of the interparticle force calculation are the Barnes-Hut opening angle θ, and the leap frog integrator time step dt. The opening angle is the force calculation accuracy parameter that determines whether an electron is treated individually or lumped together with other similarly distant particles [16]. No effect was seen from a force-softening distance p, such that the force between two particles is f 12 ∝ (r 12 + p) −2 , from p = p = a/10 6 up to p = a/100. For all simulations presented, we have demonstrated convergence using dt = τ /120, θ = 23 • , N = 10 5 .
For an initial distribution with zero temperature, the final temperature for several applicable beam densities is plotted in figure 2. We note that for extremely dense beams, n 0 > 10 19 m −3 , the DIH is comparable to the temperature of photoemission from the coldest semiconductor cathodes kT ∼ 5 meV. Furthermore, we note the oscillation of the electron temperature at 2ω p , which has been measured experimentally in UNP systems [10]. These oscillations last far longer than those seen in UNPs, as the bunch remains uniform in simulation (rather than Gaussian), and there is thus no spatial temperature smearing.
For bunches with nonzero initial temperature, the amount of additional heating can be calculated in the final coupling, plotted along with the prediction of 8 in figure 1. All simulations here had a constant density n 0 = 10 20 m −3 , with varied initial Maxwellian velocity distribution. The duration of the temperature oscillations visible in figure 2 is reduced with increasing initial temperature, though the exact dependence was not extracted in this study. In practical systems, initial acceleration in a voltage gap V will have cooled any photocathode emission temperature along the acceleration direction by a factor kT i /eV , making kT || ∼ 0 almost instantaneously in high gradient electron sources, whereas the transverse velocity spread will remain unaffected [6]. Simulations were performed for contained bunches of multiple densities with kT x = kT y = 0 and kT z = 0. Equivalent dynamics to figure 2, as well as agreement with 8 were found if the average of the temperatures in all three directions is used to calculate Γ i and Γ f , as shown in Fig 1. The timescale of equipartition here also depends on kT i and n 0 , but for the densities considered above and initial temperatures comparable to DIH heating, it was seen that equipartition occurs between τ and 3τ . A bunch instantaneously emitted into a uniform accelerating field will not have the DIH evolution altered, a fact that was also verified in simulation.
Real electron bunches just after emission are very infrequently fully transversely contained in the sense described above, and are not emitted instantaneously. The overall space charge force can double the beam radius on the timescale of ω −1 p , and acceleration can significantly lengthen the pulse afterwards. Thus, we must consider the process of DIH in bunches with time dependent density reduction. A system of fully coupled charges will have an equilibrium temperature that decreases as kT ∝ 1/Γ eq a, as given by above. However, it is well known that system of noninteracting charges expanding under linear forces will have an equilibrium temperature that decreases as kT ∝ 1/a 2 [6], as a consequence of RMS emittance conservation, which therefore leaves the phase space density or brightness unchanged. During the explosion, the beam is in a highly nonequilibrium state, but our definition of temperature above still applies, where here rv r = 0.
Simulations were performed to determine the temperature as a function of time for a bunch undergoing full Coulomb explosion. Multiple densities were considered, however the previous scaling with n 1/3 0 was seen. Thus, we present only results for n 0 = 10 20 m −3 in figure 3. The total overall expansion and velocity growth is plotted in the inset, with the analytic prediction, showing excellent agreement. The temperature shows an initial increase due to DIH at a time of t max ≈ ω −1 p , to a value of kT (t = t max ) = e 2 /4π 0 a(t max ), or a value a(t max )/a 0 ≈ 1.6 times smaller than for a fully contained beam. The bunch continues to expand, attempting to equilibrate to Γ eq = 2.23, and the temperature continues to decrease. However, at longer times, the temperature falls as a −2 , as in an uncoupled system. Thus, there must be a time at which correlation ceases to grow with decreasing temperature. Such behavior has been noted in the adiabatic expansion of UNPs due to kinetic pressure [17]. Thus, to model the temperature as a function of time, we can write the temperature as a product of the correlated growth, and the decoupled expansion after correlation ceases. Since the growth of the radius is small during the time of DIH, we may presume T = T c (t)  describes the expansion of independent particles. The rate of DIH should be proportional to the how far the system is from the equilibrium temperature at the current density, and so: The only free parameter of this model is decoupling time τ 0 , or the time when the correlation ceases to increase. Integrating this numerically, and fitting to the temperature data, we find τ 0 = 0.325τ , which is also plotted in figure 3. This time is non-negligibly greater than t max , and suggests that the bunch is cooled not only via decoupled expansion, but coupled expansion as well just after t max . We can verify this prediction by looking at the pair correlation as a function of time. This is computed in figure 4, via averaging over 5×10 3 particles randomly selected from the distribution. The development of the "Coulomb hole" and an accompanying shock profile due to violent repulsion is clearly visible, as has also been seen in UNP simulations [18]. Indeed the correlation ceases to change after a time of approximately τ 0 = 0.3τ . Partial confinement effected by linear focusing in the source would yield an altered a(t), and thus a longer τ 0 , whereas acceleration that lengthens the bunch significantly near the DIH heating timescale would yield a shorter τ 0 , but in either case the DIH heating timescale should remain unchanged.
The exploding beam case justifies post-factum the validity of the temperature limit of 8 for practical photoemission, though beams in many applications do not have constant volume. In the exploding case, near τ 0 the temperature scaling shifts from a(t) −1 as given in Eq. 2, to the familiar scaling of σ 2 v ∝ a(t) −2 , as given by emittance/brightness conservation. A more detailed calculation of DIH for practical cases in which the bunch shape is both time and space dependent should involve averaging over a temperature distribution given by 8 and 9 using the local density as determined by space charge tracking.

Conclusion
In summary, we have characterized the fundamental temperature limit of photoemission from the disorder induced heating of electrons due to poorly screened Coulomb interactions at all length scales, where the analytic two particle models fail. We have shown that the tabulated thermodynamic quantities of one component plasmas are sufficient to explain both fully contained and Coulomb exploding instances of DIH, and have verified two simple relations that describe the temperature evolution. Furthermore, many interesting effects of DIH in UNP, such as temperature oscillation and correlation decoupling, should also be present in such cold beams, yielding the possibility of rich interdisciplinary study.
Practically, for next-generation ultracold dense electron sources we have shown that the temperature of photoemission, and thus the maximum beam brightness, cannot be arbitrarily improved. Furthermore, given the rapid progress of photocathode temperature reduction, we anticipate such a limit to be approached in the next generation of high brightness electron sources producing intense beams. . RMS velocity spread of a spherical bunch undergoing coulomb expansion, treecode data (circles), and prediction of 9. The inset shows the overall expansion of the bunch in normalized units, both analytic prediction (solid line), given by Gauss's Law, and treecode data (dots). Uniform bunches remain uniform, and the normalized bunch size R(τ )/R 0 is independent of density.