Dynamical bunching and density peaks in expanding Coulomb clouds

Expansion dynamics of single-species, non-neutral clouds, such as electron bunches used in ultrafast electron microscopy, show novel behavior due to high acceleration of particles in the cloud interior. This often leads to electron bunching and dynamical formation of a density shock in the outer regions of the bunch. We develop analytic fluid models to capture these effects, and the analytic predictions are validated by PIC and N-particle simulations. In the space-charge dominated regime, two and three dimensional systems with Gaussian initial densities show bunching and a strong shock response, while one dimensional systems do not; moreover these effects can be tuned using the initial particle density profile and velocity chirp.

One application of these theories that is of particular current interest is to high-density electron clouds used in next-generation ultrafast electron microscopy (UEM) development [40][41][42]. The researchers in the UEM and the ultrafast electron diffraction (UED) communities have conducted substantial theoretical treatment of initially extremely short bunches of thousands to ultimately hundreds of millions of electrons that operate in a regime dominated by a virtual cathode (VC) limit [36,40,[43][44][45] which is akin to the Child-Langmuir current limit for beams generated under the steady-state conditions [18]. These short bunches are often generated by photoemission, and such bunches inheret an initial profile similar to that of the driving laser pulse profile. Typically, the laser pulse has an in-plane, "transverse" extent that is of order one hundred microns and a duration on the order of fifty femtoseconds, and these parameters translate into an initial electron bunch with similar trans-verse extents and sub-micron widths [40]. After photoemission, the electrons are extracted longitudinally using either a DC or AC field typically in the 1-10 MV/m [46][47][48][49] through tens of MV/m [50][51][52] ranges, respectively. However, the theoretical treatments of such "pancakelike" electron bunch evolution have largely focused on the longitudinal dimension [32][33][34][35]44], and the few studies looking at transverse dynamics have either assumed a uniform-transverse distribution [35] or have looked at the effect of a smooth Gaussian-to-uniform evolution of the transverse profile on the evolution of the pulse in the longitudinal direction [34,37]. Of specific note, only one analytic study found any indication, a weak longitudinal signal, of a shock [34].
On the other hand, an attractive theoretical observation is that an ellipsoidal cloud of cool, uniformly distributed charged particles has a linear electric field within the ellipsoid which results in maintenance of the uniform charge density as the cloud spreads [23]. In the accelerator community, such a uniform distribution is a prerequisite in employing techniques such as emittance compensation [53] as well as forming the basis of other theoretical analyses. It has long been proposed that such a uniform ellipsoid may be generated through proper control of the transverse profile of a short charged-particle bunch emitted from a source into vacuum [44], and experimental results have shown that an electron cloud emitted from a photocathode and rapidly accelerated into the highly-relativistic regime can develop into a final ellipsoidal profile characteristic of a uniform charge distribution [54]. Contrary to expectations from the free expansion work but consistent with the longitudinal analyses, this shadow lacks any indication of a peripheral region of high-density shocks. However, recent work has indicated that a substantial high density region may indeed form in the transverse direction [55], and N-particle simulation results, as demonstrated in Fig. (1), demonstrate a rapidly-developed substantial ring-like shock circumscribing the median of the bunch when the bunch starts from sufficient density. Moreover, this shock corresponds to a region of exceedingly low brightness, or conversely, high, local temperature, and that experiments show that ) are x-z and x-y projections, respectively, of just the portion of the distribution within 20% of the standard deviation of the median value of y and z, respectively. Notice the ring-like substructure that is evident in the "slices", (c.) and (d.) but absent from the full distribution projections, (a.) and (b.). (e.) N-particle radial-distributions obtained near the longitudinal median plane of the bunch at various times. Density is calculated by binning 1000 macroparticles and assigning the resulting density to the average radial position of those particles. The initial distribution is sampled from a Gaussian, and the square-like nature of the plot results from the discreteness of the bins. The sub-graph in the upper left corner shows the position of maximum density as a function of time, which is non-zero at initial time due to binning resulting in a non-zero minimum radial position. The sub-graph in the upper right shows the ratio of the maximum density to the density at the minimum r value. Notice, the "phase transition" in the 45-80 ps range where the location of the non-zero, non-stochastic peak first appears well away from the origin.
removal of this region results in a dramatic increase in the bunch brightness [55], which we term "Coulomb cooling" as it is similar to evaporative cooling in the fact that the "hottest" charged particles are removed from the distribution's edge thus leaving behind a higher-quality, cooler bunch.
To understand Coulomb cooling, we first investigate this transverse shock. Here we demonstrate the formation of a ring-like shock within N-particle simulations [56,57] of electron bunches with initial transverse Gaussian profile and offer an explanation of why this phenomena has not been noted previously within the UED literature. We then utilize a Poisson fluid approach to derive analytic predictions for the expansion dynamics in planar (1D), cylindrical, and spherical geometries, and we derive conditions for the emergence of density peaks distinct from any initial density maximum. We show that peak formation has a strong dependence on dimension, with one dimensional systems less likely to form shocks, while in cylindrical and spherical geometries bunching is more typical. Particle-in-cell (PIC) methods, utilizing WARP [58], and N-particle simulation are then used to validate the analytical predictions for peak emergence.

II. OBSERVATION OF TRANSVERSE SHOCK
One reason that a transverse shock has not been seen previously in N-particle simulations is apparent in Fig. (1). We consider pancake electron bunches typical of 100keV ultrafast electron microscopy, and we consider the thin direction of the bunch to be the zaxis. Previous studies of the expansion dynamics of these bunches, including our own work, have looked at the projection of the particle density distribution to the x-z plane [38,44,54,59,60]. Fig. (1) shows that by projecting the distribution in this manner, and with the ability to statistically discern density fluctuations at about the 10% level, results in what appears to be a uniform distribution; however, by restricting the projection to only electrons near the median of the bunch, a restriction that can only be done computationally presently, results in evidence of a transverse ring-like density substructure near the median longitudinal (z) position. To better understand when this shock emerges, simulations with the same distribution parameters (σ r ≈ 100 µm and ∆z ≈ 0.4 µm) but various numbers of electrons were run. The average radial density was calculated for 30 instances of bunches with 1 thousand, 10 thousand, and 100 thousand electrons. As can be seen in Fig. (2), the shock emergence is present for bunches with 100 thousand electrons but not for those with 1 thousand electrons. The case of bunches with 10 thousand electrons suggests the emergence of the shock, but the shock becomes less defined at later times. Fig (3) shows nine density profiles at time 10 τ p , where τ p represents the plasma period, with τ p = 2π m 0 ne 2 and n = N πσ 2 r σz . The dots in each figure indicate average densities in cylindrical rings, originating from three randomly chosen initial conditions (figures in each row) and for three values of the total number of electrons N : 1 thousand, 10 thousand, and 100 thousand electrons; for the top, middle, and bottom rows, respectively. These representative density plots support the conclusion that the shock is only present in the case of bunches with 100 thousand electrons; where the spline fit to the data indicates a significant peak removed from the center of the bunch in all instances examined. As expected, the density of bunches with one thousand electrons is noisy due to low statistics both from the small number of electrons in the simulation and the large proportion of electrons that spread beyond the analysis region due to the initial velocity spread; and the density profile of bunches with 10 thousand electron has a consistent general shape that fits well to the spline fit but lacks significant emergent peaks; which are the indicators of shock formation.
We define the emergence time as the time at which peaks indicative of a shock emerge in the dynamics of Coulomb clouds. Fig. (4) shows the dependence of the emergence time, and its variability, on the number of electrons in a bunch. It also shows very clearly that the emergence time is proportional to the plasma period, τ p . As can be seen in this figure, the spread in the emergence time is large for a bunch with 10 thousand electrons, but this spread decreases as the density of the bunch increases. For bunches with N ≥ 100, 000 the spread in the emergence time is small, moreover the emergence time appears to converge toward approximately 5 τ p at large N (for Gaussian initial distributions). We note here that for Gaussian pulses with similar spatial and temporal extents, simulations at and above 10 million electrons, a goal of the community [61], result in relativistic velocities as a result of the stronger space-charge effects. As the discussion here focuses on non-relativistic physics, we present data for up to 1 million electrons, where the velocity obtained from the self-field remains non-relativistic.
The results presented in Figs. 2-4 are a second reason that shock formation has not been seen previously in studies of electron bunches. Specifically, most work has been conducted using ≤ 10, 000 electrons with a transverse standard deviation of 100 µm, and in the regime where there is no consistent emergence of a shock. Moreover, the fact that the non-relativistic evolution of the bunch profile has a time scale proportional to the plasma period, a fact that we derive under special geometries later in this manuscript, means that higher density bunches result in faster, more consistent evolution of the transverse profile. In other words, the emergence time of a shock happens earlier as the density of the bunch is increased. Specifically, a transverse shock emerges at on the order of 50 ps for an initially Gaussian profile (σ r = 100 µm with sub-micron length) with 10 6 electrons, which is the number of electrons which is the current goal for the diffraction community [61]. This implies that for modern bunches, this transverse shock is happening well within the photoemission gun before the onset of the relativistic regime. The goal of 10 8 electrons for the imaging community needs to be further examined as the transverse velocity spread will be relativistic, but we expect to find this effect there as well, however we expect that it occurs at short times, of order a few picoseconds. For the density at 10 τp, 30 cylindrical shells of equal volume and length σz partitioned the distribution out to 0.6 mm, and the numbers of electrons in each of these shells were used to calculate a density at the shell's average radius. Due partially to the different numbers of electrons and partially due to the fact that the longer simulations, namely the simulation with N = 1000, resulted in significantly more electrons migrating out of the analysis region as a result of the initial velocity spread , the density scales are different for the three rows in the figure: for the top row, 0.1 (0.01mm) 3 for the middle row, and for the bottom row. Red dashed lines represent splines of order 3 with 10 knots. Notice the clear presence of a shock for the case N = 100, 000, an ambiguous shock at N = 10, 000, and essentially noise at N = 1, 000.

III. 1D MODEL
As noted in the introduction, formation of a shock in the longitudinal direction of an expanding pancake pulse has not been observed, and the analysis of Reed [34] demonstrates that this is true for cold initial conditions. Here we re-derive this result using an elementary method, which enables extension to include the possibility of an initial chirp; and we find chirp conditions at which shock formation in the longitudinal direction can occur.
Consider the non-relativistic spreading of an electron bunch in a one dimensional model, which is a good early time approximation to the longitudinal spreading of a pancake-shaped electron cloud generated at a photocath-FIG. 4. The emergence time divided by the plasma period as a function of the number of electrons in the initial Gaussian profile with σr ≈ 100 µm and σz ≈ 0.4µm. Emergence time was determined as the first time the density away from the inner-most-value exceeded the inner-most-value by 2%. Notice that the emergence time converges to about 5 τp for high densities, but at low densities the emergence time has high variability with a median shifted to higher multiples of the plasma period.
ode. In one dimensional models, the density, ρ, only depends on one coordinate, which we take to be z. We also take ρ to be normalized so that its integral is one. For the sake of readability, denote the position of a particle from the Lagrangian perspective to be z = z(t) and z 0 = z(0). The acceleration of a Lagrangian particle is where q is the charge of the particle (e.g. electron), m is its mass, Q tot is the total charge in the bunch, and The key observation enabling analytic analysis is that if the flow of electrons is lamellar, so that there is no crossing of particle trajectories, then these integrals and the acceleration calculated from them are time independent and hence may be determined from the initial distribution. Therefore, we denote a(z; t) = a 0 , ρ 0 = ρ(z; 0), and δσ = z0 −z0 ρ 0 (z)dz. Moreover, due to the fact that for any particle trajectory, the acceleration is constant and given by a 0 , the Lagrangian particle dynamics reduces to the elementary constant acceleration kinematic equation where v 0 is the initial velocity of the charged particle that has initial position z 0 . Notice that both v 0 and a 0 are functions of the initial position, z 0 , and we shall see later that the derivatives of these parameters, v = dv0 dz0 and a = da0 dz0 are important in describing the relative dynamics of Langrangian particles starting at different initial positions. Moreover, the special case of v 0 = 0 everywhere, which we will call the cold-case, is commonly assumed in the literature, and we now examine this case in detail.
First we consider the speading charge distribution within the Eulerian perspective, where z is an independent variable instead of it describing the trajectory of a particle. We denote the charge distribution at all times to be Q tot ρ(z; t) with ρ(z; t) a unitless, probability-like density and Q tot the total charge per unit area in the bunch. Since particle number is conserved, we have so that in the non-relativistic case derived above Notice that the derivative of the acceleration with respect to the initial position is directly proportional to the initial distribution, so that for initial distributions that are symmetric about the origin, Plugging Eq. (6) into Eq. (5), we get For the cold-case, Eq. (7) reduces to the density evolution equation derived by Reed [34] using different methods. Also, in the cold-case, Eq. (8) simplifies into a proportionality between the initial slope of the distribution and the slope of the distribution at any later time. Therefore, a charge distribution that is initially at rest and unimodal, i.e only a single initial location has non-zero density with ρ 0 = 0, never develops a dynamically generated second maximum. This explains why we should not expect to see an emergent shock in the longitudinal direction; provided the 1D model is applicable and cold initial conditions are valid.
However, if particles in the initial state have an initial velocity that depends on initial position, i.e. v 0 (z 0 ), then density peaks will emerge at z when t = This occurs at positive time if v 0 ρ 0 > v 0 ρ 0 . In the special case v 0 ρ 0 = v 0 ρ 0 , the distribution may be reframed as a cold-case distribution starting from t = − mc1 0 qQtot for some z 0 -independent constant c 1 with velocity units when mc 2 1 0 qQtot < 1 or a distribution starting from a singularity with velocity distributionṽ 0 = c 1 ρ 0 − 1 2 δσ . As noted earlier, the function a 0 (z 0 ) is monotonically increasing as a function of distance from the center of the pulse, which means that electrons at the edges of the bunch always have larger accelerations away from the center of the pulse than electrons nearer the pulse center. Thus crossover, where an inner electron moves past an outer electron, cannot occur unless the initial velocities of inner electrons overcome this relative acceleration.
A practical case where crossover may be designed is where the initial distribution has an initial velocity chirp, i.e. v 0 = cz 0 where c has units of inverse time. Intuitively we can expect that the velocity chirp needs to be negative in order for crossover to occur. To find the crossover time, we consider the time at which two electrons that were initially apart, are at the same position at the same time. In this case it is straight forward to find the time at which crossover occurs by considering an electron at initial position z 0 , and a second electron at position z 0 +δz 0 . Before either of these electrons experiences a crossover Eq. (3) is valid, and setting z(z 0 , t x ) = z(z 0 + dz 0 , t x ) reduces to, where A = q 2m 0 ρ 0 (z 0 ) and B = dz dz0 . Solving the quadratic equation leads to the crossover time given by Since A is always positive, the square root is real only if B 2 is larger than 4A. Moreover the time is only positive if B is negative. Therefore crossover only occurs if the chirp has a negative slope, as expected on physical grounds. The conditions for tuning the chirp to produce crossover in 1D are then The results above are applicable to the spreading in the longitudinal direction of non-relativistic pancake bunches, because the expression Eq. (3) is linear in acceleration. In that case, the position of a charged particle at any time can be calculated from a superposition of the contribution from the space-charge field and any external constant field such as a constant and uniform extraction field. In that case, the space charge field leads to spreading of the pulse, while the extraction field leads solely to an acceleration of the center of mass of the entire bunch. In that case the center of mass and spreading dynamics are independent and can be decoupled. The extension of the description above to asymmetric charge density functions is also straightforward, as is the inclusion of an image field at the photocathode. Moreover, inclusion of these effects does not change Eq. (7), Eq. (8), nor the conclusions we have drawn from them. These results apply generally to all times before the initial crossover event within the evolution of the bunch, and once crossover occurs, the distribution can be reset with a new Eq. (7) to follow further density evolution.
As we show in the next section the one-dimensional results do not apply, even qualitatively, to higher dimensions, as the constant acceleration situation is not valid and crossover can occur even with cold initial conditions, as demonstrated in the simulations presented in previous section. In the next section we present fluid models in higher dimensions where the origin of these new effects is evident.

IV. CYLINDRICAL AND SPHERICAL MODELS
The methodology for the cylindrical and spherical systems is similar so we develop the analysis concurrently.
Consider a non-relativistic evolving distribution Q tot ρ(r, t), where ρ(r, t) is again taken to be the unitless particle distribution and Q tot is again the total charge in the bunch. In a system with cylindrical symmetry, the mean field equation of motion for a charge at r ≡ r(t) is given by where p r is the momentum of a Lagrangian particle, λ(r, t) is the cummulative distribution function (cdf) λ(r, t) = r 0 2πrρ(r, t)dr (14) and Q tot λ(r; t) is the charge inside radius r. Analogously in a system with spherical symmetry, the equation of motion for an electron at position r is given by where P (r, t) = r 0 4πr 2 ρ(r, t)dr is the cdf and Q tot P (r; t) is the charge within the spherical shell of radius r. Notice, r in Eq. (13) denotes the cylindrical radius while r in Eq. (15) represents the spherical radius. In both cases, before any crossover occurs, the cdf of a Lagrangian particle is constant in time. For simplicity we write λ(r, t) = λ(r 0 , 0) ≡ λ 0 and P (r, t) = P (r 0 , 0) ≡ P 0 for a particle starting at r 0 ≡ r(0). In other words, since Q tot λ 0 and Q tot P 0 can be interpreted as the charge contained in the appropriate Gaussian surface, if we track the particle that starts at r 0 , these contained charges should remain constant before crossover occurs. It is convenient to also define the average particle density to beρ 0 = λ0 πr 2 0 in the cylindrically symmetric case andρ 0 = 3P0 4πr 3 0 in the spherically symmetric case. Notice that these average particle densities are a function solely of r 0 , and we will use these parameters shortly. Eq. (13) and Eq. (15) may now be rewritten as, for the cylindrical and spherical cases respectively dp r dt = qQ tot λ 0 2π 0 r , which apply for the period of time before particle crossover. Note that unlike the one dimensional case, in two and three dimensional systems the acceleration on a Lagrangian particle is not constant, and has a time dependence through the time dependent position r = r(r 0 : t) term in the denominator. Since Eq. (16) and Eq. (17) represent the force on the particle in the cylindrical and spherical contexts, respectively, we can integrate over the particle's trajectory to calculate the change in the particle's energy. Integrating from r 0 , 0 to r, t gives for the cylindrical and spherical cases respectively where the term on the right side of the equality can be interpreted as the change in the potential energy within the self-field of the bunch. These expressions are fully relativistic, and in the nonrelativistic limit, we can derive implicit position-time relations for the particle by setting the energy difference equal to the non-relativisitic kinetic energy mv 2 /2, and integrating. The details of this derivation have been placed in Appendix B, and the resulting expressions in the cold-case for the cylindrical and spherical systems are respectively t =τ where F (·) represents the Dawson function andτ p,0 represents the plasma period determined from the initial conditions:τ p,0 = 2π m 0 qQtotρ0 = 2π ω0 , indicating that the appropriate time scale is the scaled plasma period as seen in Figs. (2) and (4) for the case of pancake bunches used in ultrafast electron diffraction systems. Eq. (21) and its derivation is equivalent to previous time-position relations reported in the literature [26,62] although the previous work did not identify the plasma period as the key time-scale of Coulomb spreading processes and cylindrical symmetry was not discussed (Eq. (20)).
The time-position relations detailed in the equations above depend solely on the amount of charge nearer to the origin than the point in question, i.e. Q totρ0 , and not on the details of the distribution. Notice however, that it is the difference between the time-position relationships of different locations where the details of the distribution become important and may cause neighboring particles to have interesting relative dynamics; leading to the possibility of shock formation in the density.
To translate the Lagrangian particle evolution equations above to an understanding of the dynamics of the charge density distribution, we generalize Eq. (5) to where d is the dimensionality of the problem, i.e. 1 (planar symmetry), 2 (cylindrical symmetry), or 3 (spherical symmetry). The factor in the denominator, dr dr0 , may be determined implicitly from the time-position relations above, and the details are presented in Appendix C. The resulting expressions for the density dynamics, in the cold case, for d = 2 (cylindrical) and d = 3 (spherical) cases are where is a function only of the initial position. The functions f d for cylindrical systems is given by while for systems with spherical symmetry we find Note that these are functions of the ratio r/r 0 . The functions f d can also be written as mixed functions of r and t, specifically f 2 r r0 = r0 r ln r r0 ω 0 t and However, care must be used when using these mixed forms as r is implicitly dependent on t. Here we work with these functions in terms of relative position, r r0 . Substituting Eq. (23) into Eq. (22), we find that the density evolution in systems with cylindrical (d=2) and spherical symmetry(d=3) can be compactly written as This expression is general and can be applied to arbitrary, spherically symmetric or cylindrically symmetric initial conditions. Analogously to the 1D case the condition dr dr0 < 0 results in particle crossover. However, as detailed in Eq. (23), the sign of dr dr0 depends on the sign of 1 + D d (r 0 )f d r r0 . It is very interesting to note that D d (r 0 ) is the deviation from a uniform distribution function, so that the D functions are solely functions of the initial conditions and are positive at locations where the local density is larger than the average density at r 0 , and negative when the local density is smaller than the average density at r 0 . On the other hand, the functions f d are functions of the evolution of the Langrangian particle. One immediate consequence of Eq. (27) is that for a uniform initial density distribution, for either cylindrical or spherical systems, the corresponding D function is zero at every location where the original density is defined. Thus, the uniform density evolution in Eq. (22) reduces to the generally recognized expressions: ρ(r, t)πr 2 = ρ 0 πr 2 0 for the cylindrical case; and ρ(r, t) 4 3 πr 3 = ρ 0 4 3 πr 3 0 in the spherical case. We provide additional details for the uniform distribution in the next section. However, Eq. (23) is general for any distribution before particle crossover, not just the uniform distribution.
For a particle starting at position r 0 and having a deviation from uniform function D d (r 0 ), crossover occurs when the particle is at a position, r, that satisfies f d r r0 = −1/D d (r 0 ). Since every particle moves toward positive r, every particle will have a time for which it will assume every value of the function f r r0 ). The character of the two and three dimensional f 's are similar as can be seen in Fig (5) where the value of the function is plotted against r r0 . Specifically, both functions increase to a maximum and then asymptote towards 1 from above. This means that all density positions eventually experience uniform-like scaling since lim r→∞ f d ( r r0 ) = 1 results in Eq. (22) simplifying to ρπr 2 = ρ 0 r0) in the cylindrical and spherical cases, respectively, for large enough r. Notice, this uniform-like scaling does not mean that the distribution goes to the uniform distribution, which is what happens in 1D but need not happen under cylindrical and spherical geometries.
The main difference between the cylindrical and spherical symmetries are that the cylindrical function's maximum is larger than the spherical function's maximum; and we find max(f 2 ) ≈ 1.28 while max(f 3 ) ≈ 1.07. Moreover the maximum of the cylindrical function occurs at a larger value of r0 r than that of the spherical function; specifically r ≈ 9.54r 0 instead of r ≈ 8.27r 0 , respectively. The first observation means cylindrical symmetry is more sensitive to the distribution than the spherical case, while the second observation indicates that if crossover is go- ρ 0 (r 0 ) − 1 as is done in (c-d.). Specifically, crossover may occur at some time when dr dr 0 < 0. For negative values of D d crossover first occurs at the maximum of f d and the value of D d at which that occurs is given in the schematics of (c,d). For values of D d that are more negative than this value, crossover will occurs at some time in some parts of the distrubiton. For the uniform distribution, D(r0) = 0 for all points inside the distribution. The dashed blue line indicates expansion less quickly than the uniform distribution, while the green line indicates more rapid expansion than for the uniform distribution.
ing to occur for a specific particle, it will occur before the r value for which the corresponding f function is maximum (i.e. r ≈ 9.54r 0 or r ≈ 8.27r 0 ), otherwise the particle will never experience crossover. From this reasoning, we obtain the earliest time for crossover by minimizing the time taken for a trajectory to reach the maximum of the function f d , with the crossover constraint dr dr0 = 0. This may be achieved by using Lagrange multipliers or by running calculations for a series of values of r/r 0 to find the position at which crossover happens first. The mean field theory is valid before the minimum crossover time, and the results presented below are well below this time.

V. UNIFORM AND GAUSSIAN EVOLUTIONS: THEORY AND SIMULATION
In this section, the mean field predictions are compared to the N-particle and PIC simulations. First we present the evolution of the initially-at-rest cylindrically-and spherically-symmetric uniform distribution of 1.875×10 7 and 2 × 10 4 electrons within radii's of 1 mm (see Fig.  (6(a,b)). Note, in this fairly trivial case, crossover should not occur and the analytic results should be valid meanfield-results for all time. Since ρ 0 =ρ 0 in this case, D d (r 0 ) = 0 and Eq. (27) reduces to Notice that r can be solved for a specific time using Eq. (20) or Eq. (21), depending on whether we are examining the cylindrically-or spherically-symmetric case, respectively, and due toρ 0 's independence from r 0 , these equations need only be solved once for a give time to describe all r. Therefore, we may write r = α(t)r 0 ≡ αr 0 , where α is independent of r 0 , and we immediately see that Eq. (28) can be written as ρ(r; t) = α d ρ 0 (r 0 ) suggesting that the density simply scales with time as generally recognized by the community. We solve for α at 6 times, and present a comparison with both PIC and N-particle cylindrically-symmetric and spherically-symmetric simulations in Fig. (6(a,b)). As can be seen, despite the presence of initial density fluctuations arising from sampling, the simulated results follow the analytic results exceedingly well. Specifically, the distributions simply expand while remaining essentially uniform, and the analytic mean field formulation correctly calculates the rate of this expansion. While this comparison is arguably trivial, it is reassuring to see that our general equation reduces to a form that captures these dynamics. Less trivial is the evolution of Gaussian distributions. We simulated 3.75 × 10 7 and 10 5 electrons for the cylindrical and spherical cases, respectively, using σ r = 1 mm. Solving for the minimum crossover time, we get approximately 44 ns for each distribution. Therefore, we simulate for 37.5 ns, which is well before any crossover events.
For the Gaussian distributions we introduce the scaled radius variables s = r √ 2σr and s 0 = r0 √ 2σr , so that from Eq. (25) for the cylindrical and spherical cases we have, where erf is the well known error function. Putting these expressions into Eq. (27) we find for the cylindrical and shows the analytic position of max density as a function of time, and the sub-graph in the upper right of (c.d.) shows the analytic ratio of the max density to the density at the minimum r value both. The corresponding analytic ratio for the uniform distribution is shown in this sub-graph as a dashed horizontal line at 1. Unsurprisingly, the PIC results and the analytical results, both mean-field models, are in almost perfect agreement, and the N-particle results are in surprisingly good agreement as well. Notice that the models predict peak formation on a time-scale dependent on the initial plasma frequency similar to the peak formation seen in the N-particle disc-like density evolution seen in Fig. (1(e)) and detailed in Fig. (4).
To find r(t)/r 0 we solve Eq. (20) or Eq. (21), depending on whether we are examining the cylindrically-or spherically-symmetric cases, respectively; and for r r0 , and for every time step, we calculate the predicted distribution at 5000 positions, r, corresponding to 5000 initial positions, r 0 , evolved to time t. As can be seen in Fig. (6), both the cylindrically-and spherically-symmetric Gaussian distributions develop peaks similar to those seen in the simulations of expanding pancake bunches described in the first section of the paper. As can be seen in Fig.  (6), both the PIC and the N-particle results match the analytical results very well. Notice, the primary differences between the cylindrically-and spherically-symmetric evolutions is in their rate of width expansion and the sharpness of the peak that forms, and both of these facets are captured by the analytic models.

VI. CONCLUSIONS
In this work, we have shown that a shock occurs in the transverse, but not longitudinal, direction during expansion of pancake-like charged particle distributions typical of those use in ultrafast electron microscope (UEM) systems.
Fluid models for arbitrary initial distributions, Eq. (7), a generalization of a model already in the literature, showed that the formation of such a shock should not occur for any cold initial distribution in one dimension. This result is consistent with the finding that typically no shock is visible in the longitudinal direction dynamics of UEM bunches; however by tuning the initial velocity distribution it should be possible to generate a dynamic shock.
We generalized the fluid theory to cylindrical and spherical symmetries deriving implicit evolution equations for the charge density distributions Eq. (27). We analyzed these models for the advent of particle crossover, which occur for some distributions even when the initial distribution is cold due to the behavior of the Coulomb force in higher dimensions; and we found that the time scales associated with the space charge expansion are proportional to the plasma period. One interesting detailed observation is that in the case of cylindrical symmetry, the pre-factor of τp π of Eq. (20) is roughly 0.3 while for the spherical symmetric case, corresponding prefactor in Eq. (21), 3 2 τp 2π , is roughly 0.2 plasma periods. Interestingly beam relaxation has been independently found to occur at roughly 0.25 the plasma period [63], which falls directly in the middle of our cylindrically and spherically symmetric models. The analytic theory predicts that emergence of a shock is distribution dependent, and as expected, a uniform initial distribution does not produce a shock. However we showed that electron bunches that are initially Gaussian distributed produce a shock well before the advent of particle crossover indicating that the emergence of a shock is well described by fluid models presented here. This is consistent with the observation of a shock in N-particle simulations of the transverse expansion of UEM pancake bunches (see Figs. 1-4).
To our knowledge, we have presented the first analytic derivation of the cold, single-species, non-neutral density evolution equations for cylindrical and spherical symmetries. These equations are general enough to handle any distribution under these symmetries, and can be used across specialties from accelerator technology, to electronics, to astrophysics. While simulation methods, like the N-particle and PIC codes used here are general tools, the insights provided by these simple analytical equations should provide fast and easy first-approximations for a number of calculations; while providing physical insights and parameter dependences that are more difficult to extract from purely computational studies.
The analysis presented here has been carried out for the non-relativistic regime; which is only valid for cases of sufficiently low density where the shock occurs prior to the electrons achieving relativistic velocities. For higher densities or other physical situations where the bunch becomes relativistic more quickly than the formation of this shock, a relativistic analysis is needed. On the other hand, for sufficiently high densities, i.e. approaching 10 7 or more electrons in the pancake geometry used in this manuscript and typical in the UE field, relativistic effects in the transverse direction become important and need to be considered. The extension to fully relativistic cases will be addressed in future work.
We point out that Child-Langmuir current should not have these dynamic shocks except at the onset of the current before the steady-state condition sets in. Previous studies note the "hollowing" of a steady-state beam due to fringe field effects [64], but a steady state Child-Langmuir current is largely independent of emission parameters; so that this hollowing effect is not dynamical, but part of the continuous emission process itself, and is therefore a very different mechanism than the dynamic shocks we see here. It would be interesting to study the combined effects of steady state beam hollowing and dynamic shock formation in pancake bunches to determine if the combination of these processes provides new opportunities for optimization of beam properties.
The analytic models presented here treat free expansion whereas most applications have lattice elements to confine the bunches. Substantial work, in particular the particle-core model, has been very successful at predicting transverse particle halos of beams [65,66]. This model assumes a uniform-in-space beam-core density called a Kapchinsky-Vladimirsky (KV) distribution due to its ease of theoretical treatment. Such an assumption is supported by the analysis presented here as we find that the distribution within the shock is nearly uniform. However, the particle-core models do not treat the initial distribution as having a large density on the periphery. It would be interesting to revisit such treatments with this new perspective although we would like to point out that the main effect the particle-core model attempts to capture, halos, occur even after aperturing the beam [65]. Specifically, it should be possible to examine the effect of radialfocussing fields on the evolution of the three-dimensional distributions we have investigated here.
The experimental work that motivated this analysis, [55], not only predicted a shock but also a correlated decrease in brightness near the periphery. We emphasize that the mean-field equations used here explain the density shock only and do not provide a quantitative theory of the emittance and the Coulomb cooling achieved by removing the electrons in the shock. Specifically, the true emittance in the analytic models presented here remains zero for all time as all particles at a radius r have velocity 2 3 r0 ωp,0 1 − r0 r resulting in zero local spread in velocity space. This perfect relationship between velocity and position means that the true emittance is zero even if the relation is non-linear; however, in such a nonlinear chirp case, the rms emittance will not remain zero despite the true emittance being zero. Moreover, the analytic model does capture some of the rms emittance growth as a change in the distribution has a corresponding change on the variance measures used to determine the rms emittance. Specifically, a Gaussian distribution should have especially large emittance growth due to its evolution to a bimodal distribution, a distribution that is specifically problematic for variance measures. Such a large change in the emittance of the transverse Gaussian profile has been seen by Luiten [44] experimentally and us computationally [37,38]. On the other hand, the perfectly uniform distribution does not change its distribution throughout its evolution and therefore should have zero rms emittance growth as the chirp exactly cancels out the expansion of the pulse at all times. Moreover, Luiten et al. found experimentally that the uniform distribution does have an increase in emittance although less than the Gaussian casel [44], an observation that is corroborated by our own work with PIC and N-particle calculations [37,38]. The analytic formulation of meanfield theory presented here provides new avenues to treating emittance growth, by treating fluctuations to these equations in a systematic manner. This analysis will be presented elsewhere. where u = aλ 0 ln r r0 ,w =ũ √ aλ0 , and w = ln r r0 . The remaining integral,