Fundamental photoemission brightness limit from disorder induced heating

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 [1]. For UED setups, the brightness of the electron beam is the main beam parameter that determines the resolving power of the instrument [2].

Such electron sources are comprised of a photoemitting material placed in an accelerating gradient, where both direct current and radio frequency accelerating fields are used for various applications. We will restrict our focus in this work to emitters without nano- or micro-scale structure, whose photoemission probability is uniform across the cross sectional area where light illuminates its surface. For either x-ray or electron diffraction experiments, it is often the four dimensional transverse normalized brightness that is most pertinent. For a given beam current I, we can define the 'microbrightness' as the phase space density itself, $\rho =\frac {{\mathrm { d}} I}{{\mathrm { d}} V_4}$, where dV₄ = dx dy dp_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:

$$\begin{equation} \mathcal{B}_{\mathrm{n,4D }}=\frac{I^{-1}}{\left(mc\right)^2}\int \rho(x, y, p_x, p_y)^2 \,{\rm d}V_4. \end{equation} \tag{ 1 }$$

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 [3], it was shown that the maximum achievable beam brightness (either microbrightness or total) can be written as

$$\begin{equation} \mathcal{B}_{\mathrm{n,4D }} \bigg|_{\mathrm{max}}=\frac{mc^2 f\epsilon_0 E_{\mathrm{acc}}}{2\pi kT}, \end{equation} \tag{ 2 }$$

where E_acc is the accelerating electric field directly at the photocathode, which sets the maximum supportable charge density at the photocathode. For such sources, the minimum divergence is set by kT, the temperature of photoemitted electrons, an intrinsic property of the choice of photoemitter and laser wavelength. In the absence of emitter imperfections such as surface roughness, 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. This is achieved either via the cryogenic cooling of the photocathode [4], 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 [5]. Furthermore, using laser-cooled atoms which are photoionized, there has been production of electrons at temperatures equivalent to ∼ 1 meV [6].

The number of electrons per bunch varies widely across various applications, with a range approximately between 10⁶ and 10⁹ electrons, which for moderate to high flux applications often corresponds to densities in the range of n₀ = 10¹⁷ → 10²⁰ m⁻³. For example, the Cornell energy recover linac injector prototype operates with a charge density of 5 × 10¹⁷ m⁻³ as described in [7], whereas photoguns operating in the blowout regime may have densities higher than 10¹⁹ m⁻³, calculated from the parameters given in [8]. 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²/4π₀a, where a is the Wigner–Seitz radius, a = (3/4πn₀)^1/3. Thus for a given density, we expect the heating to be of the order e²/4π₀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²/4π₀akT = 0.2–2 given an electron temperature of 5 meV, and for the density ranges described above. 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.

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, ∇²Φ(r) = en(r)/₀, 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 $\lambda =\sqrt {\epsilon _0 k T/n_0 {e}^2}$ is much larger than the interparticle separation. However, for the lowest temperatures attained in semiconductor photoemission kT ∼ 5 meV [5], and the lowest of the above densities n₀ = 10¹⁷ m⁻³, 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 [9, 10]. 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.

Nondiffusive, two-particle methods have also been developed that do not require the assumption of Debye screening. A clear presentation of these is given in [11]. 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²〉 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 variables $\tilde {r}= r/r_{\mathrm {i}}$ and $\tilde {t}=t \left ( e r_{\mathrm {i}}^{-3/2}/\sqrt {2\pi \epsilon _0\,m}\right )$ the equation of motion for their separation is given by

$$\begin{equation} 1=\left(\frac{{\rm d}\tilde{r}}{{\rm d}\tilde{t}}\right)^2+\frac{1}{\tilde{r}}. \end{equation} \tag{ 3 }$$

This equation is integrable for the function $\tilde {t}(\tilde {r})$, which is not analytically invertible, but is trivial to invert numerically to obtain $\tilde {r}(\tilde {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:

$$\begin{equation} \langle \Delta v ^2 \rangle = \int n(\mathbf r)\Delta v(r, \delta t_{\rm c})^2\, {\rm d}^3\mathbf r =\frac{2 {e}^2 n_0}{\epsilon_0\,m} \int r^2{\rm d}r \left( \frac{{\rm d}\tilde{r}(\tilde{t})}{{\rm d}\tilde{t}}\frac{1}{\sqrt{r}} \right)^2, \end{equation} \tag{ 4 }$$

where $\tilde {t}$ 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² dr ensures that the averaging does not diverge at small r, and thus we can integrate over all space. Equation (4) is the statistical average of all two particle interactions in a beam during some time δt_c.

2.1. Unbound trajectories

2.2. 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²〉 ∼ n₀, whereas the effect we desire to calculate should have 〈v²〉 ∼ n^1/3₀, argued above. This scaling of the velocity spread with density occurs in the Fokker–Planck method (as shown in [9] in equation (5.243)), in the ETPA (as shown in [11] in equation (5.5.11)), and in the work by Massey (see [12] in equation (19)).

A method that produces the correct scaling with density, but is also flawed, is again given by Jansen in [11]. In what he calls the 'thermodynamic limit', Jansen takes the difference of the total electrostatic potential energy 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 total thermal energy, $\frac {3}{2}NkT$. However, instead of using a single particle density, Jansen implicitly uses the two-particle correlated density for Debye screening:

$$\begin{equation} \frac 3 2 nV k T=U_{\mathrm{i}}-U_{\mathrm{f}}=\frac 1 2 n^2\,V \int^\infty _0 {\rm d}^3 r \frac{e^2}{4 \pi \epsilon_0 r}\left(1- \exp[-\phi(r)/kT] \right). \end{equation} \tag{ 5 }$$

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 separated by r, and is assumed to have the form: $\phi (r) =e^2\exp (-r/\lambda )/4\pi \epsilon _0 r$. The factor $n\exp (-\phi (r)/kT)$ is the final two particle correlated density. This equation has a solution of the form

$$\begin{equation} kT = \frac{e^2}{4 \pi \epsilon_0} \left(4 \pi \alpha^2\right)^{1/3} n^{1/3}, \end{equation} \tag{ 6 }$$

where α is a dimensionless number determined by numerical solution of the above, α = 0.087 02.⁵ This number corresponds to a coupling factor at thermodynamic equilibrium of Γ_eq = 3.53. However, if one evaluates the number of particles in the Debye sphere, one finds $N_{\mathrm { D}}=\frac {4}{3}\pi \lambda ^3 n_0 \approx 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 correlated density in general. It must be found by some other numerical means [13].

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 $\tau =2\pi \omega _{\mathrm {p}}^{-1}=2 \pi \left (n_0 e^2/ m_{\mathrm {e}} \epsilon _0 \right )^{-1/2}$ the temperature rises to an equilibrium nonzero value higher than the initial gas temperature. 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 [14], and we will argue that an exact analogue of this effect is present in practical photoemission.

DIH has been experimentally observed for the ions in a neutral plasma [15]. In cold neutral plasmas, the electrons equilibrate much faster than the ionic ω⁻¹_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 [16]. 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 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:

$$\begin{equation} \frac{U_{\mathrm{c}}}{NkT} \equiv \bar U= \frac{1}{2}\frac{e^2 n_0}{4\pi \epsilon_0 } \int \frac{ g(r, \Gamma, \kappa)\,{\rm d}^3r}{r}. \end{equation} \tag{ 7 }$$

This binding energy is balanced by an increasing ion temperature. Calculating the increase in temperature thus only requires knowledge of $\bar U$ for a given electron screening. Analytic calculation of g(r) for strongly coupled systems is difficult, and alternatively, $\bar U$ has been tabulated via molecular dynamics (MD) simulations [16].

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 [17]. 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 OCP. In the case of an OCP, there is no second species to provide additional screening, and thus $\bar U$ is only a function of Γ. Owing to this simplification, MD data has been fit to a power series relation for $\bar U=a\Gamma +b\Gamma ^{1/3}+c$, where the coefficients are given in [17], a = −0.90, b = 0.590 673, c = −0.265 69. We may now write down the expression for the final coupling Γ_f, and thus the final equilibrium temperature for a fully confined, initially uncorrelated ($\bar U_{\mathrm { i}} =0$), uniform distribution of charges with initial coupling Γ_i, analogous to that presented in [16]:

$$\begin{equation} \bar U_{\rm f} =\frac{3}{2} \left(\frac{\Gamma_{\rm f}}{\Gamma_{\rm i}}-1 \right)= a\Gamma_{\rm f}+b\Gamma_{\rm f}^{1/3}+c. \end{equation} \tag{ 8 }$$

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.

Zoom In Zoom Out Reset image size

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 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 ω⁻¹_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/γ² damping of the interaction force [9]. 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⁻⁹(n₀ (m⁻³))^1/3. Along with (2), this forms the fundamental photoemission 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 tree-algorithms have been developed, see for instance [18], and one has been included in a relativistic particle accelerator code⁶. However, considering the simplicity of the problem, we elected to directly modify a code originally written for gravitational interaction [19], simply by inverting the sign of the interaction force and changing its scale.

We consider a randomly distributed spherical bunch of 10⁵ particles with uniform density, and vary only the density (10¹⁷ → 10²⁰ m⁻³) and initial Gaussian velocity distribution (kT∈{0.25,20} meV or kT = 0), where we use the definition $kT=m_{\mathrm {e}}\left (\langle v_i ^2\rangle -\langle x_i v_i \rangle ^2/\langle x_i^2 \rangle \right )$, where m_e is the electron mass and i is a Cartesian coordinate. This definition of temperature subtracts position-correlated motion, as in the bulk flow or expansion of charge, from the total rms velocity distribution width. For a fully contained uniform bunch there is no correlated flow, thus 〈x_iv_i〉 ≈ 0, and the temperature is just a measure of average electron kinetic energy, with no spatial dependence. Bunches are contained 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 numerical integration time step dt. In any Barnes–Hut simulation, the opening angle θ = s/D is the force calculation accuracy parameter that determines whether an electron in a cubic volume of size s, at a distance D from the test particle is treated individually or lumped together with other electrons in that volume [20]. An opening angle of zero corresponds to a fully point-to-point summation of the interaction force, whereas the large angle θ = 1 corresponds to lumping electrons together in a box as wide as it is distant from the test particle. No effect was seen from a force-softening distance p, such that the force between two particles is f₁₂∝(r₁₂ + p)⁻², from p = a/10⁶ up to p = a/100. While p is not exactly equivalent to the r_c of two particle methods, the insensitivity of simulation results to this parameter does illustrate the significance of multiple strong interactions for closely spaced particles, which contain their initially rapid expansion. For all simulations presented, we have demonstrated convergence using dt = τ/120, θ = 23°, N = 10⁵.

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₀ > 10¹⁹ m⁻³, 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 [15]. These oscillations last far longer than those seen in UNPs, as the bunch density remains uniform in simulation (rather than evolving into a Gaussian), and there is thus no spatial temperature smearing.

Zoom In Zoom Out Reset image size

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₀ = 10²⁰ m⁻³, 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 [9]. 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 figure 1. The timescale of equipartition here also depends on kT_i and n₀, 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 ω⁻¹_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/Γ_eqa, as given 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² [9], 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₀ was seen. Thus, we present only results for n₀ = 10²⁰ m⁻³ in figure 3. The total overall expansion and velocity growth is plotted in the inset, with the prediction of the equation of motion given by Gauss' law (integrated numerically), showing excellent agreement. The temperature shows an initial increase due to DIH at a time of t_max ≈ ω⁻¹_p, to a value of kT(t = t_max) = e²/4π₀a(t_max), or a value a(t_max)/a₀ ≈ 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⁻², 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 [21].

Zoom In Zoom Out Reset image size

Thus, to phenomenologically model the temperature as a function of time, we can factorize the temperature as a product of the correlated growth and the decoupled expansion after correlation ceases. Thus, $T=T_{\mathrm { c}}(t)\frac {a_0^2}{a(t)^2}$, where T_c is the correlated temperature dependence, and the term $\frac {a_0^2}{a(t)^2}$ describes the expansion of independent particles. Since the growth of the radius is small during the time of DIH, at early times T(t)∝T_c, and later times as T(t)∝a(t)⁻² with T_c unchanging or changing slowly. To model T_c, we know that the rate of DIH should dependent on the difference of the current temperature T(t) to the equilibrium temperature at that radius, $T_{\mathrm { eq}}=e^2/\left (4\pi \epsilon _0 a(t) \Gamma _{\mathrm { eq}}\right )$. The simplest ansatz one can make is a linear dependence of the DIH temperature growth rate on the temperature difference between the current temperature to the equilibrium temperature:

$$\begin{equation} \frac{{\rm d}T_{\rm c}}{{\rm d}t}=\frac{T_{{\rm c},{\rm eq}}-T_{\rm c}(t)}{\tau_0}=\frac{1}{\tau_0}\left( \frac{e^2}{4 \pi \epsilon_0 a(t) \Gamma_{\rm eq}} -T_{\rm c}(t) \right). \end{equation} \tag{ 9 }$$

The only free parameter of this model is decoupling time scale τ₀, or the time scale when the correlated heating ceases. Integrating this numerically and fitting to the temperature data, we find τ₀ = 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³ 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 [22]. Indeed the correlation ceases to change after a time of approximately τ₀ = 0.3τ. Partial confinement effected by linear focusing in the source would yield an altered a(t), and thus a longer τ₀, whereas acceleration that lengthens the bunch significantly near the DIH heating timescale would yield a shorter τ₀, but in either case the DIH heating timescale should remain unchanged.

Zoom In Zoom Out Reset image size

The above case demonstrates that even for practical dense bunches that do not have constant volume, DIH will still play a large role. For a changing beam density, the temperature is well described by (8) until the beam density decreases enough to decouple the interaction at τ₀. After τ₀ the temperature scaling shifts away from a(t)⁻¹ as given in (8), to the familiar scaling of σ²_v∝a(t)⁻², 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.

In summary, we have characterized the fundamental temperature limit of photoemission from unstructured emitters caused by DIH 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.

It is important to note that it may be possible circumvent the above limit with engineered photoemitters that produce electrons in a pre-ordered state, such as has been suggested in metal photocathodes [23], and in cold gas electron sources [24]. The full degree to which the effects of DIH can be mitigated remains an open question. As such structured emitters are arrays of smaller sub-emitters, unless each sub-emitter emits only one electron at a time, the DIH limit applies to each one individually. Furthermore, it may be possible to create beams of higher brightness downstream of the source by emitting a low density beam for which DIH effects are negligible, and employing transverse and longitudinal bunch compression after relativistic acceleration, for which the interparticle interaction force will be relativistically reduced. However, the above brightness limit should hold in the nonrelativistic regime near an unstructured source.

In practical terms, we have shown that the brightness of a next-generation cold electron unstructured photoemission source cannot arbitrarily increase its brightness. 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.

The authors thank two anonymous reviewers for helpful suggestions that clarified the presentation of the work. The authors also thank Joshua Barnes, the original author of the treecode package modified for use in this work. This work was supported by the National Science Foundation grant no. DGE-0707428, and by the Department of Energy grant no. DE-SC0003965. This work was also supported by the Director, Office of Science, Office of Basic Energy Sciences of the US Department of Energy, under contract numbers DE-AC02-05CH11231, KC0407-LSJNT-I0013 and DE-SC0005713.