---
abstract: '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.'
author:
- 'B. S. Zerbe'
- 'X. Xiang'
- 'C.-Y. Ruan'
- 'S. M. Lund'
- 'P. M. Duxbury'
bibliography:
- 'CoulombDynamics.bib'
title: Dynamical bunching and density peaks in expanding Coulomb clouds
---

=1

​

Introduction
============

Non-neutral plasma systems arise in a variety of physical contexts ranging from astrophysics[@Arbanil:2014_charged_star_review; @Maurya:2015_charged_sphere; @Yousefi:2014_dust_aggregates]; accelerator technologies [@Bacci:2014_plasma_acceleration; @Boine:2015_intense_beams; @Whelan:2016_MRI; @Bernal:2016_recirculator]; ion and neutron production [@Bulanov:2002_charged_beam_generation; @Fukuda:2009_species_generation; @Esirkepov:2004_highly_efficient_ion_generation; @Kaplan:2015_preprint; @Parks:2001_neutron_production; @Bychenkov_2015_review]; sources for electron and ion microscopy[@Murphy:2014_cold_ions; @Gahlmann:2008_ultrashort]; to high power vacuum electronics[@Booske:2011_vacuum_review; @Liu:2015_maximal_charge; @Zhang:2016_review]. Understanding of the dynamics of spreading of such systems is critical to the design of next generation technologies, and simple analytic models are particularly helpful for instrument design. As a result, substantial theoretical efforts have already been made in this vein[@Jansen:1988_book; @Reiser:1994_book; @Batygin:2001_self; @Bychenkov:2005_coulomb_explosion; @Grech:2011_coulomb_explosion; @Kaplan:2003_shock; @Kovalev:2005_kinetic_spherically_coulomb_explosion; @Last:1997_analytic_coulomb_explosion; @Eloy:2001_coulomb_explosion; @Krainov:2001_ce_dynamics; @Morrison:2015_slow_down_dynamics; @Boella:2016_multiple_species]. Specifically, free expansion of clouds of charged single-specie particles starting from rest have been well studied both analytically and computationally[@Last:1997_analytic_coulomb_explosion; @Eloy:2001_coulomb_explosion; @Grech:2011_coulomb_explosion; @Batygin:2001_self; @Degtyareva:1998_gaussian_pileup; @Siwick:2002_mean_field; @Qian:2002_fluid_flow; @Reed:2006_short_pulse_theory; @Collin:2005_broadening; @Gahlmann:2008_ultrashort; @Tao:2012_space_charge; @Portman:2013_computational_characterization; @Portman:2014_image_charge; @Michalik:2006_analytic_gaussian], and a number of studies have found evidence of the formation of a region of high-density, often termed a “shock”, on the periphery of the clouds under certain conditions[@Grech:2011_coulomb_explosion; @Kaplan:2003_shock; @Kovalev:2005_kinetic_spherically_coulomb_explosion; @Last:1997_analytic_coulomb_explosion; @Murphy:2014_cold_ions; @Reed:2006_short_pulse_theory; @Degtyareva:1998_gaussian_pileup].

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[@King:2005_review; @Hall:2014_report; @Williams:2017_longitudinal_emittance]. 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[@Valfells:2002_vc_limit; @Luiten:2004_uniform_ellipsoidal; @King:2005_review; @Miller:2014_science_review; @Tao:2012_space_charge] which is akin to the Child-Langmuir current limit for beams generated under the steady-state conditions[@Zhang:2016_review]. 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 transverse extents and sub-micron widths[@King:2005_review]. After photoemission, the electrons are extracted longitudinally using either a DC or AC field typically in the 1-10 MV/m[@Srinivasan:2003_UED; @Ruan:2009_nanocrystallography; @van_Oudheusden:2010_rf_compression_experiment; @Sciaini:2011_review] through tens of MV/m[@Musumeci:2010_single_shot; @Weathersby:2015_slac; @Murooka:2011_TED] ranges, respectively. However, the theoretical treatments of such “pancake-like” electron bunch evolution have largely focused on the longitudinal dimension[@Luiten:2004_uniform_ellipsoidal; @Siwick:2002_mean_field; @Qian:2002_fluid_flow; @Reed:2006_short_pulse_theory; @Collin:2005_broadening], and the few studies looking at transverse dynamics have either assumed a uniform-transverse distribution[@Collin:2005_broadening] 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[@Reed:2006_short_pulse_theory; @Portman:2013_computational_characterization]. Of specific note, only one analytic study found any indication, a weak longitudinal signal, of a shock[@Reed:2006_short_pulse_theory].

  ---------- -- ---------
                
   \[-3mm\]      \[8cm\]
  ---------- -- ---------

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 [@Grech:2011_coulomb_explosion]. In the accelerator community, such a uniform distribution is a prerequisite in employing techniques such as emittance compensation[@Rosenzweig:2006_emittance_compensation] 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[@Luiten:2004_uniform_ellipsoidal], 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[@Musucemi:2008_generate_uniform_ellipsoid]. 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[@Williams:2017_transverse_emittance], and N-particle simulation results, as demonstrated in Fig. (\[fig:distribution substructure\]), 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 removal of this region results in a dramatic increase in the bunch brightness[@Williams:2017_transverse_emittance], 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 [@Berz:1987_cosy; @Zhang:2015_fmm_cosy] 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[@Friedman:2014_warp], and N-particle simulation are then used to validate the analytical predictions for peak emergence.

Observation of Transverse Shock
===============================

One reason that a transverse shock has not been seen previously in N-particle simulations is apparent in Fig. (\[fig:distribution substructure\]). We consider pancake electron bunches typical of 100keV ultrafast electron microscopy, and we consider the thin direction of the bunch to be the z-axis. 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[@Luiten:2004_uniform_ellipsoidal; @Musucemi:2008_generate_uniform_ellipsoid; @Portman:2014_image_charge; @Morrison:2013_measurement; @Li:2008_quasiellipsoidal]. Fig. (\[fig:distribution substructure\]) 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.

  -- -- --
        
  -- -- --

![\[fig:splined\_grid\] Density near the z-median of simulated pancake bunches with transverse Gaussian profiles ($\sigma_r = 100 ~ \mu$m) in an extraction field of 10 MV/m. Each figure is the transverse radial density of a section of width $\sigma_z \approx 0.4~\mu$m for different initial conditions and different numbers of electrons, at time $10 ~ \tau_p$ where $\tau_p = 2 \pi \sqrt{\frac{m \epsilon_0}{n e^2}}$ is the plasma frequency; where $n = \frac{N}{\pi \sigma_r^2 \sigma_z}$. The number of electrons in each horizontal panel is different and equal to $N = 1,000$ (top), $N=10,000$ (middle), and $N = 100,000$ (bottom). For the density at $10 ~ \tau_p$, 30 cylindrical shells of equal volume and length $\sigma_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: $\frac{1}{(0.1 mm)^3}$ for the top row, $\frac{0.1}{(0.01 mm)^3}$ for the middle row, and $\frac{1}{(0.01 mm)^3}$ 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$.](grid_of_splined_2.png){width="45.00000%"}

  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --
   ![\[fig:emergence time\] The emergence time divided by the plasma period as a function of the number of electrons in the initial Gaussian profile with $\sigma_r \approx 100 ~\mu$m and $\sigma_z \approx 0.4 \mu$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~\tau_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. ](emergence_time.png "fig:"){width="45.00000%"}  
  -------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --

To better understand when this shock emerges, simulations with the same distribution parameters ($\sigma_r \approx 100 ~\mu$m and $\Delta z \approx 0.4~\mu$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. (\[fig:average distribution\]), 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 (\[fig:splined\_grid\]) shows nine density profiles at time $10 ~ \tau_p$, where $\tau_p$ represents the plasma period, with $\tau_p = 2 \pi \sqrt{\frac{m \epsilon_0}{n e^2}}$ and $n = \frac{N}{\pi \sigma_r^2 \sigma_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. (\[fig:emergence time\]) 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, $\tau_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\ge 100,000$ the spread in the emergence time is small, moreover the emergence time appears to converge toward approximately $5~\tau_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[@BES_report:2016_electron_sources], 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 $\le 10,000$ electrons with a transverse standard deviation of 100 $\mu$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 ($\sigma_r = 100 ~\mu$m with sub-micron length) with $10^6$ electrons, which is the number of electrons which is the current goal for the diffraction community[@BES_report:2016_electron_sources]. 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.

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 [@Reed:2006_short_pulse_theory] 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 photocathode. In one dimensional models, the density, $\rho$, only depends on one coordinate, which we take to be $z$. We also take $\rho$ 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 $$\begin{aligned}
  a(z;t) = \frac{qQ_\text{tot}}{2 m \epsilon_0} \delta \sigma\end{aligned}$$ where $q$ is the charge of the particle (e.g. electron), $m$ is its mass, $Q_{\text{tot}}$ is the total charge in the bunch, and $$\begin{aligned}
  \delta \sigma 
                       = \int_{-z}^{z} \rho(\tilde{z};t) d\tilde{z}\end{aligned}$$ 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$, $\rho_0 = \rho(z;0)$, and $\delta \sigma =  \int_{-z_0}^{z_0} \rho_0(\tilde{z})d\tilde{z}$. 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 $$\label{eq:non-relativistic position}
  z(t) = z_0 + v_0 t + {1\over 2} a_0 t^2$$ 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' = \frac{dv_0}{dz_0}$ and $a' = \frac{da_0}{dz_0}$ 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_\text{tot} \rho(z; t)$ with $\rho(z; t)$ a unitless, probability-like density and $Q_\text{tot}$ the total charge per unit area in the bunch. Since particle number is conserved, we have $$\begin{aligned}
  \rho(z; t)dz = \rho_0 dz_0 \end{aligned}$$ so that in the non-relativistic case derived above $$\begin{aligned}
  \rho(z,t) &=  \rho_0 \left({dz\over dz_0}\right)^{-1} =\frac{ \rho_0}{1 + v_0' t +\frac{1}{2} a_0' t^2} \label{eq:1D density evolution}\end{aligned}$$ 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, $$\begin{aligned}
\label{eq:first order acceleration}
  a_0' &= \frac{q Q_\text{tot}}{2 m \epsilon_0}  \frac{d{\delta \sigma}}{dz_0}\nonumber\\
          &= \frac{q Q_\text{tot}}{m \epsilon_0}\rho_0\end{aligned}$$ Plugging Eq. (\[eq:first order acceleration\]) into Eq. (\[eq:1D density evolution\]), we get $$\begin{aligned}
  \rho(z,t) &=  { \rho_0 \over  1 + v_0' t +  \frac{q Q_\text{tot}}{2 m \epsilon_0} \rho_0t^2}\label{eq:1D density evolution subbed}\\
  \frac{d}{dz} \rho(z,t) &= \frac{ \rho_0'\left( 1  + v_0' t\right) - \rho_0v_0''t }{\left(1 + v_0' t +  \frac{q Q_\text{tot}}{2 m \epsilon_0} \rho_0t^2\right)^3}\label{eq:1D rho derivative} \end{aligned}$$ where $\rho_0' = \frac{d\rho_0}{dz_0}$ and $v_0'' = \frac{d^2v_0}{dz_0^2}$. A detailed derivation of the second expression is in Appendix \[ap:1D appendix\].

For the cold-case, Eq. (\[eq:1D density evolution subbed\]) reduces to the density evolution equation derived by Reed[@Reed:2006_short_pulse_theory] using different methods. Also, in the cold-case, Eq. (\[eq:1D rho derivative\]) 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 $\rho_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 = \frac{1}{v_0''\rho_0 - v_0' \rho_0'}$. This occurs at positive time if $v_0'' \rho_0 > v_0' \rho_0'$. In the special case $v_0'' \rho_0 = v_0' \rho_0'$, the distribution may be reframed as a cold-case distribution starting from $t = -\frac{m c_1 \epsilon_0}{q Q_\text{tot}}$ for some $z_0$-independent constant $c_1$ with velocity units when $\frac{m c_1^2 \epsilon_0}{q Q_\text{tot}} < 1$ or a distribution starting from a singularity with velocity distribution ${\tilde{v}}_0 = c_1 \left( \rho_0 - \frac{1}{2} \delta\sigma\right)$. 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 = c z_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 + \delta z_0$. Before either of these electrons experiences a crossover Eq. (\[eq:non-relativistic position\]) is valid, and setting $z(z_0, t_x) = z(z_0 + dz_0, t_x)$ reduces to, $$\begin{aligned}
  At^2_x + Bt_x + 1 = 0\end{aligned}$$ where $A = \frac{q}{2m \epsilon_0} \rho_0(z_0)$ and $B = \frac{dz}{dz_0}$. Solving the quadratic equation leads to the crossover time given by $$\begin{aligned}
  t_x &= \frac{-B \pm \sqrt{B^2 - 4A}}{2A}\end{aligned}$$ 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 $$\begin{aligned}
  \frac{dv_0}{dz_0} &< 0\label{eq:negative chirp}\\
  \left| \frac{dv_0}{dz_0}\right| &\ge \sqrt{\frac{2q}{m\epsilon_0}\rho_0(z_0)}\label{eq:sufficient chirp})\end{aligned}$$

The results above are applicable to the spreading in the longitudinal direction of non-relativistic pancake bunches, because the expression Eq. (\[eq:non-relativistic position\]) 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. (\[eq:1D density evolution subbed\]), Eq. (\[eq:1D rho derivative\]), 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. (\[eq:1D density evolution subbed\]) 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.

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_\text{tot} \rho(r,t)$, where $\rho(r,t)$ is again taken to be the unitless particle distribution and $Q_\text{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 \equiv r(t)$ is given by $$\label{eq:2D second derivative}
\frac{d \vec{p}_\text{r}}{dt} = { q Q_\text{tot} \lambda (r,t)\over 2 \pi \epsilon_0 r}\hat{r},$$ where $\vec{p}_\text{r}$ is the momentum of a Lagrangian particle, $\lambda(r,t)$ is the cummulative distribution function (cdf) $$\lambda(r,t) = \int_0^{r} 2\pi \tilde{r} \rho(\tilde{r},t) d\tilde{r}$$ and $Q_\text{tot} \lambda(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 $$\label{eq:3D second derivative}
\frac{d \vec{p}_\text{r}}{dt} = { q Q_\text{tot} P(r,t)\over4\pi \epsilon_0 r^2}\hat{r},$$ where $P(r,t) = \int_0^{r} 4\pi \tilde{r}^2 \rho(\tilde{r},t) d\tilde{r}$ is the cdf and $Q_\text{tot} P(r;t)$ is the charge within the spherical shell of radius $r$. Notice, $r$ in Eq. (\[eq:2D second derivative\]) denotes the cylindrical radius while $r$ in Eq. (\[eq:3D second derivative\]) 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 $\lambda(r,t) = \lambda(r_0,0) \equiv \lambda_0$ and $P(r,t) = P(r_0,0) \equiv P_0$ for a particle starting at $r_0 \equiv r(0)$. In other words, since $Q_\text{tot} \lambda_0$ and $Q_\text{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 $\bar{\rho}_0 = \frac{\lambda_0}{\pi r_0^2 }$ in the cylindrically symmetric case and $\bar{\rho}_0 = \frac{3 P_0}{4\pi r_0^3 }$ 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. (\[eq:2D second derivative\]) and Eq. (\[eq:3D second derivative\]) may now be rewritten as, for the cylindrical and spherical cases respectively $$\begin{aligned}
  \frac{d p_r}{dt} &= { q Q_\text{tot} \lambda_0\over 2 \pi \epsilon_0 r}\label{eq:2D equation of motion}, \\ 
  \frac{d p_r}{dt} &= { q Q_\text{tot} P_0\over 4 \pi \epsilon_0 r^2}\label{eq:3D equation of motion}, \end{aligned}$$ 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. (\[eq:2D equation of motion\]) and Eq. (\[eq:3D equation of motion\]) 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 $$\begin{aligned}
  E(r,t) - E(r_0,0) &= \frac{q Q_\text{tot} \lambda_0}{2 \pi \epsilon_0} ln\left(\frac{r}{r_0}\right)\label{eq:2D energy}\\
  E(r,t) - E(r_0,0) &= \frac{q Q_\text{tot} P_0}{4 \pi \epsilon_0} \left(\frac{1}{r_0} - \frac{1}{r}\right)\label{eq:3D energy}\end{aligned}$$ 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 non-relativistic 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 \[ap:time-location appendix\], and the resulting expressions in the cold-case for the cylindrical and spherical systems are respectively $$\begin{aligned}
  t &= \frac{\bar{\tau}_{p,0}}{\pi} \frac{r}{r_0} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)}\right)\label{eq:2D time}\\
 t &= \sqrt{\frac{3}{2}}\frac{\bar{\tau}_{p,0}}{2\pi} \left( \tanh^{-1} \left( \sqrt{1 -  \frac{r_0}{r}} \right) + \frac{r}{r_0}\sqrt{1 -  \frac{r_0}{r}}\right)\label{eq:3D time}\end{aligned}$$ where $F(\cdot)$ represents the Dawson function and $\bar{\tau}_{p,0}$ represents the plasma period determined from the initial conditions: $\bar{\tau}_{p,0} = 2\pi \sqrt{\frac{m\epsilon_0}{q Q_\text{tot} \bar{\rho}_0}} = \frac{2\pi}{{\bar{\omega}}_0}$, indicating that the appropriate time scale is the scaled plasma period as seen in Figs. (\[fig:average distribution\]) and (\[fig:emergence time\]) for the case of pancake bunches used in ultrafast electron diffraction systems. Eq. (\[eq:3D time\]) and its derivation is equivalent to previous time-position relations reported in the literature[@Boyer:1989_kinetic_energy; @Last:1997_analytic_coulomb_explosion] 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. (\[eq:2D time\])).

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_\text{tot} \bar{\rho}_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. (\[eq:1D density evolution\]) to $$\begin{aligned}
  \rho(r,t) = \rho_0 \left( \left(\frac{r}{r_0}\right)^{d-1} {\frac{dr}{dr_0}}\right)^{-1}\label{eq:density evolution general}\end{aligned}$$ 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, ${\frac{dr}{dr_0}}$, may be determined implicitly from the time-position relations above, and the details are presented in Appendix \[ap:spatial derivatives\]. The resulting expressions for the density dynamics, in the cold case, for $d=2$ (cylindrical) and $d=3$ (spherical) cases are $$\begin{aligned}
  \frac{d r}{dr_0} &= \frac{r}{r_0} \left(1 + D_d(r_0) f_\text{d}\left(\frac{r}{r_0}\right)\right)\label{eq:dr over dr_0 no v_0}\end{aligned}$$ where $$D_d = D_d(r_0) = \frac{d}{2} \left(\frac{\rho_0}{\bar{\rho}_0} - 1\right),$$ is a function only of the initial position. The functions $f_d$ for cylindrical systems is given by $$f_\text{2}\left(\frac{r}{r_0}\right) = 2  \sqrt{\ln\left(\frac{r}{r_0}\right)}F\left(\sqrt{\ln\left(\frac{r}{r_0}\right)}\right),$$ while for systems with spherical symmetry we find $$f_\text{3}\left(\frac{r}{r_0}\right) = \frac{r_0}{r} \sqrt{1 - \frac{r_0}{r}} \tanh^{-1}(\sqrt{1 - \frac{r_0}{r}}) + 1 - \frac{r_0}{r}.$$ 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\left(\frac{r}{r_0}\right) = \frac{r_0}{r}\sqrt{\ln\left(\frac{r}{r_0}\right)}{\bar{\omega}}_0t$ and $f_3\left(\frac{r}{r_0}\right) = \frac{r_0}{r}\sqrt{1 - \frac{r_0}{r}}{\bar{\omega}}_0t$. 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, $\frac{r}{r_0}$. Substituting Eq. (\[eq:dr over dr\_0 no v\_0\]) into Eq. (\[eq:density evolution general\]), we find that the density evolution in systems with cylindrical (d=2) and spherical symmetry(d=3) can be compactly written as $$\begin{aligned}
  \rho(r;t) &= \left(\frac{r_0}{r}\right)^d \frac{\rho_0}{1 + D_d(r_0)  f_\text{d}\left(\frac{r}{r_0}\right)}\label{eq:general density evolution}\end{aligned}$$ This expression is general and can be applied to arbitrary, spherically symmetric or cylindrically symmetric initial conditions.

Analogously to the 1D case the condition $\frac{dr}{dr_0} < 0$ results in particle crossover. However, as detailed in Eq. (\[eq:dr over dr\_0 no v\_0\]), the sign of $\frac{dr}{dr_0}$ depends on the sign of $1 + D_d(r_0) f_\text{d}\left(\frac{r}{r_0}\right)$. It is very interesting to note that $ D_\text{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.(\[eq:general density evolution\]) 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. (\[eq:density evolution general\]) reduces to the generally recognized expressions: $\rho(r,t) \pi r^2=  \rho_0 \pi r_0^2$ for the cylindrical case; and $\rho(r,t) \frac{4}{3} \pi r^3=  \rho_0 \frac{4}{3} \pi r_0^3$ in the spherical case. We provide additional details for the uniform distribution in the next section. However, Eq. (\[eq:dr over dr\_0 no v\_0\]) 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_\text{d}(r_0)$, crossover occurs when the particle is at a position, $r$, that satisfies $f_\text{d}\left(\frac{r}{r_0}\right) = -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\left(\frac{r}{r_0}\right))$. The character of the two and three dimensional $f$’s are similar as can be seen in Fig (\[fig:classifying D\]) where the value of the function is plotted against $\frac{r}{r_0}$. 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 \to \infty} f_{d}(\frac{r}{r_0}) = 1$ results in Eq. (\[eq:density evolution general\]) simplifying to $\rho \pi r^2 =  \rho_0 \frac{\pi r_0^2}{1 + D_\text{d}(r_0)}$ and $\rho \frac{4}{3} \pi r^3=  \rho_0 \frac{\frac{4}{3} \pi r_0^3}{1 + D_\text{d}(r_0)}$ 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_\text{2}) \approx 1.28$ while $max(f_\text{3}) \approx 1.07$. Moreover the maximum of the cylindrical function occurs at a larger value of $\frac{r_0}{r}$ than that of the spherical function; specifically $r \approx 9.54 r_0$ instead of $r \approx 8.27 r_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 going to occur for a specific particle, it will occur before the $r$ value for which the corresponding $f$ function is maximum (i.e. $r \approx 9.54 r_0$ or $r \approx 8.27 r_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 $ {\frac{dr}{dr_0}} = 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.

  -- --
     
     
  -- --

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 \times 10^7$ and $2 \times 10^4$ electrons within radii’s of 1 mm (see Fig. (\[fig:density evolution\](a,b)). Note, in this fairly trivial case, crossover should not occur and the analytic results should be valid mean-field-results for all time. Since $\rho_0 = {\bar{\rho}}_0$ in this case, $D_d(r_0) = 0$ and Eq. (\[eq:general density evolution\]) reduces to $$\begin{aligned}
  \rho(r;t) = \left(\frac{r_0}{r}\right)^d \rho_0(r_0)\label{eq:uniform evolution}\end{aligned}$$ Notice that $r$ can be solved for a specific time using Eq. (\[eq:2D time\]) or Eq. (\[eq:3D time\]), depending on whether we are examining the cylindrically- or spherically- symmetric case, respectively, and due to ${\bar{\rho}}_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 = \alpha(t) r_0 \equiv \alpha r_0$, where $\alpha$ is independent of $r_0$, and we immediately see that Eq. (\[eq:uniform evolution\]) can be written as $\rho(r;t) = \alpha^d \rho_0(r_0)$ suggesting that the density simply scales with time as generally recognized by the community. We solve for $\alpha$ at 6 times, and present a comparison with both PIC and N-particle cylindrically- symmetric and spherically-symmetric simulations in Fig. (\[fig:density evolution\](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 \times 10^7$ and $10^5$ electrons for the cylindrical and spherical cases, respectively, using $\sigma_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 = \frac{r}{\sqrt{2} \sigma_r}$ and $s_0 = \frac{r_0}{\sqrt{2}\sigma_r}$, so that from Eq. (25) for the cylindrical and spherical cases we have, $$\begin{aligned}
   D_2(s_0) &= \frac{(1 + s_0^2) e^{-s_0^2} - 1}{1 - e^{-s_0^2}}\\
  D_3(s_0) &= \frac{(2 s_0^3 + 3 s_0)e^{-s_0^2} - \frac{3 \sqrt{\pi}}{2} \text{erf}(s_0)}{\sqrt{\pi}\text{erf}(s_0) - 2 s_0 e^{-s_0^2}}\end{aligned}$$

where erf is the well known error function. Putting these expressions into Eq. (\[eq:general density evolution\]) we find for the cylindrical and

spherical cases respectively $$\begin{aligned}
  \rho(s;t) &= \frac{\frac{s_0^2}{\pi s^2} e^{-s_0^2}}{1 + 2  \frac{(1 + s_0^2) e^{-s_0^2} - 1}{1 - e^{-s_0^2}} \sqrt{\ln\left(\frac{s}{s_0}\right)} F\left(\sqrt{\ln\left(\frac{s}{s_0}\right)}\right)}\\
  \rho(s;t) &= \frac{\frac{s_0^3}{\pi^{\frac{3}{2}} s^3} e^{-s_0^2}}{1 + \frac{(2 s_0^3 + 3 s_0)e^{-s_0^2} - \frac{3 \sqrt{\pi}}{2} \text{erf}(s_0)}{\sqrt{\pi}\text{erf}(s_0) - 2 s_0 e^{-s_0^2}} \left(\frac{s_0}{s}\sqrt{1 - \frac{s_0}{s}} \tanh^{-1}\left(\sqrt{1 - \frac{s_0}{s}}\right) + 1 - \frac{s_0}{s}\right)}\end{aligned}$$

To find $r(t)/r_0$ we solve Eq. (\[eq:2D time\]) or Eq. (\[eq:3D time\]), depending on whether we are examining the cylindrically- or spherically- symmetric cases, respectively; and for $\frac{r}{r_0}$, 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. (\[fig:density evolution\]), 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. (\[fig:density evolution\]), 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.

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. (\[eq:1D density evolution subbed\]), 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. (\[eq:general density evolution\]). 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 $\frac{\tau_p}{\pi}$ of Eq. (\[eq:2D time\]) is roughly 0.3 while for the spherical symmetric case, corresponding prefactor in Eq. (\[eq:3D time\]), $\sqrt{\frac{3}{2}}\frac{\tau_p}{2\pi}$, is roughly 0.2 plasma periods. Interestingly beam relaxation has been independently found to occur at roughly 0.25 the plasma period[@Wangler:1985_emittance_relaxation], 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[@Luginsland:1996_child_langmuir_2d], 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[@Gluckstern:1994_analytic_halo; @Wangler:1998_particle_core_review]. 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[@Gluckstern:1994_analytic_halo]. Specifically, it should be possible to examine the effect of radial-focussing fields on the evolution of the three-dimensional distributions we have investigated here.

The experimental work that motivated this analysis, [@Williams:2017_transverse_emittance], 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 $\sqrt{\frac{2}{3}}\frac{r_0}{{\bar{\omega}}_{p,0}} \sqrt{1 - \frac{r_0}{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 non-linear 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[@Luiten:2004_uniform_ellipsoidal] experimentally and us computationally[@Portman:2013_computational_characterization; @Portman:2014_image_charge]. 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[@Luiten:2004_uniform_ellipsoidal], an observation that is corroborated by our own work with PIC and N-particle calculations[@Portman:2013_computational_characterization; @Portman:2014_image_charge]. The analytic formulation of mean-field 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.

***Acknowledgment*** This work was supported by NSF Grant 1625181, the College of Natural Science, the College of Communication Arts and Sciences, and the Provost’s office at Michigan State University. Computational resources were provided by the High Performance Computer Center at MSU. We thank Martin Berz and Kyoko Makino for their help with employing COSY as an N-particle code.

1D Density Derivative {#ap:1D appendix}
=====================

In the main text, we argued $$\begin{aligned}
  \rho(z)=  \rho_0(z_0) \left({dz\over dz_0}\right)^{-1}\end{aligned}$$ To determine the slope of the density, we take the derivative with respect to the z coordinate $$\begin{aligned}
  \frac{d}{dz} \rho(z) &= \frac{d}{dz_0} \left(\rho_0(z_0) \left(\frac{dz}{dz_0}\right)^{-1}\right)\left(\frac{dz}{dz_0}\right)^{-1}\nonumber\\
                                &= \frac{d}{dz_0} \left(\rho_0(z_0)\right) \left(\frac{dz}{dz_0}\right)^{-2} -  \rho_0(z_0) \left(\frac{dz}{dz_0}\right)^{-3} \frac{d^2z}{dz_0^2}\nonumber\\
                                &= \frac{\frac{d}{dz_0} \left(\rho_0(z_0)\right) \frac{dz}{dz_0} - \rho_0(z_0)\frac{d^2z}{dz_0^2}}{\left(\frac{dz}{dz_0}\right)^{3}}\label{eq:generic z derivative}\end{aligned}$$ For the sake of conciseness, denote $\rho_0 = \rho_0(z_0)$, $\rho' = \frac{d}{dz} \rho(z)$, $\rho_0' = \frac{d}{dz_0} \rho_0$, $v_0' = \frac{d v_0}{d z_0}$, and $v_0'' = \frac{d^2 v_0}{d z_0^2}$. From the main text, we have $$\begin{aligned}
  \frac{dz}{dz_0} &= 1 + v_0' t +  \frac{q}{2 m \epsilon_0} \rho_0(z_0)t^2\end{aligned}$$ and from this it is straightforward to show $$\begin{aligned}
  \frac{d^2z}{dz_0^2} &= v_0'' t +  \frac{q}{2 m \epsilon_0} \rho_0't^2\end{aligned}$$ Subbing this back into Eq. (\[eq:generic z derivative\]), we get $$\begin{aligned}
  \rho' &= \frac{\rho_0'\left(1 + v_0' t +  \frac{q}{2 m \epsilon_0} \rho_0t^2\right) - \rho_0\left(v_0'' t +  \frac{q}{2 m \epsilon_0} \rho_0't^2\right)}{\left(1 + v_0' t +  \frac{q}{2 m \epsilon_0} \rho_0t^2\right)^3}\nonumber\\
          &= \frac{ \rho_0'\left( 1  + v_0' t\right) - \rho_0v_0''t }{\left(1 + v_0' t +  \frac{q}{2 m \epsilon_0} \rho_0t^2\right)^3}\label{eq:1D rho derivative appendix} \end{aligned}$$

Derivation of Time-location Relations {#ap:time-location appendix}
=====================================

Integral form
-------------

Starting with the relativistic expression for change in particle energy derived in the main text $$\begin{aligned}
  \text{cyl:  }E(t) - E(0) &= \frac{q Q_{tot} \lambda_0}{2 \pi \epsilon_0} ln\left(\frac{r}{r_0}\right)\label{eq:2D energy}\\
  \text{sph: }E(t) - E(0) &= \frac{q Q_{tot} P_0}{4 \pi \epsilon_0} \left(\frac{1}{r_0} - \frac{1}{r}\right)\label{eq:3D energy}\end{aligned}$$ we approximate the energy change with a change in non-relativistic kinetic energy starting from rest $$\begin{aligned}
  \text{cyl: }\frac{1}{2} m v^2 &= \frac{q Q_{tot} \lambda_0}{2 \pi \epsilon_0} ln\left(\frac{r}{r_0}\right)\label{eq:2D kinetic energy}\\
   \text{sph: }\frac{1}{2} m v^2 &= \frac{q Q_{tot} P_0}{4 \pi \epsilon_0} \left(\frac{1}{r_0} - \frac{1}{r}\right)\label{eq:3D kinetic energy}\end{aligned}$$ where $v = \frac{d r}{dt}$ are the velocity of the particle at time $t$ in the two or one of the three dimensional models, respectively, with the appropriate definition of $r$. Solving these equations for the velocity at time $t$, we get $$\begin{aligned}
  \text{cyl: } \frac{d r}{d t} &= \sqrt{\frac{q Q_{tot} \lambda_0}{\pi m \epsilon_0} ln\left(\frac{r}{r_0}\right)}\label{eq:2D velocity}\\
  \text{sph: } \frac{d r}{d t} &= \sqrt{\frac{q Q_{tot} P_0}{2 \pi m \epsilon_0} \left(\frac{1}{r_0} - \frac{1}{r }\right)}\label{eq:3D velocity}\end{aligned}$$ Separating the variables and integrating, we obtain $$\begin{aligned}
  \text{cyl: } t &= \int_{r_0}^{r} \frac{d\tilde{r}}{ \sqrt{\frac{q Q_{tot} \lambda_0}{\pi m \epsilon_0} ln\left(\frac{\tilde{r}}{r_0}\right)}}\label{eq:2D time integral}\\
  \text{sph: } t &= \int_{r_0}^{r} \frac{d\tilde{r}}{\sqrt{\frac{q Q_{tot} P_0}{2 \pi m \epsilon_0} \left(\frac{1}{r_0} - \frac{1}{\tilde{r}}\right)}}\label{eq:3D time integral}\end{aligned}$$ Defining $a = \frac{q Q_\text{tot}}{\pi m \epsilon_0}$, we rewrite Eq. (\[eq:2D time integral\]) and Eq. (\[eq:3D time integral\]) as $$\begin{aligned}
  \text{cyl: }t &= \int_{r_0}^{r} \frac{d\tilde{r}}{ \sqrt{a \lambda_0 ln\left(\frac{\tilde{r}}{r_0}\right)}}\label{eq:2D time integral subbed}\\
  \text{sph: }t &= \int_{r_0}^{r} \frac{d\tilde{r}}{\sqrt{\frac{aP_0}{2r_0}}\sqrt{1 - \frac{r_0}{\tilde{r} }}}\label{eq:3D time integral subbed}\end{aligned}$$

Cylindrically-symmetric integral solution
-----------------------------------------

We solve the sylindrically-symmetric integral first. Define $\tilde{u} = \sqrt{a \lambda_0 ln\left(\frac{\tilde{r}}{r_0}\right)}$. Solving this equation for $\tilde{r}$ in terms of $\tilde{u}$, we see that $\tilde{r} = r_0 e^{\frac{\tilde{u}^2}{a \lambda_0}}$. It is also straightforward to see that $$\begin{aligned}
  d \tilde{u} &= \frac{1}{2} \frac{1}{ \sqrt{a \lambda_0 ln\left(\frac{\tilde{r}}{r_0}\right)}} \frac{a \lambda_0}{\tilde{r}} d \tilde{r}\nonumber\\
                  &= \frac{1}{ \sqrt{a \lambda_0 ln\left(\frac{\tilde{r}}{r_0}\right)}} \frac{a \lambda_0}{2 r_0} e^{\frac{-\tilde{u}^2}{a \lambda_0}} d \tilde{r} \nonumber\end{aligned}$$ Applying this change of coordinates to Eq. (\[eq:2D time integral subbed\]), we get $$\begin{aligned}
  \text{cyl: } t &= \int_{0}^{u} \frac{2 r_0}{a \lambda_0}e^{\frac{\tilde{u}^2}{a \lambda_0}} d \tilde{u}\nonumber\\
                    &= \frac{2 r_0}{a \lambda_0}e \int_{0}^{u} e^{\frac{\tilde{u}^2}{a \lambda_0}} d \tilde{u}\nonumber\\
                    &= \frac{2 r_0}{\sqrt{a \lambda_0}} \int_{0}^{w} e^{\tilde{w}^2} d \tilde{w}\label{eq:2D time integral in w}\end{aligned}$$ where $u = \sqrt{a \lambda_0 ln\left(\frac{r}{r_0}\right)}$, $\tilde{w} = \frac{\tilde{u}}{\sqrt{a \lambda_0}}$, and $w = \sqrt{ln\left(\frac{r}{r_0}\right)}$. The remaining integral, $\int_{0}^{w} e^{\tilde{w}^2} d \tilde{w}$ can be written in terms of the well-studied Dawson function, $F(\cdot)$: $$\begin{aligned}
\int_{0}^{w} e^{\tilde{w}^2} d \tilde{w} &= e^{w^2}F(w) \nonumber\\
                                                           &= \frac{r}{r_0} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)}\right)\label{eq:integral as Dawson}\end{aligned}$$ Subbing Eq. (\[eq:integral as Dawson\]) back into Eq. (\[eq:2D time integral in w\]) gives us our time-position relation $$\begin{aligned}
\label{eq:2D time vs r}
  \text{cyl: } t &= \frac{2 r}{\sqrt{a \lambda_0}} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\nonumber\\
                     &= 2 \frac{r}{r_0} \sqrt{\frac{\pi r_0^2 m \epsilon_0}{q Q_\text{tot}\lambda_0}} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\nonumber\\
                     &= 2  \frac{r}{r_0} \sqrt{\frac{m \epsilon_0}{q Q_\text{tot}\frac{\lambda_0}{\pi r_0^2}}}F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\nonumber\\
                     &= 2  \frac{r}{r_0} \sqrt{\frac{m \epsilon_0}{q Q_\text{tot}{\bar{\rho}}_0}}F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\nonumber\\
                     &= \frac{2}{{\bar{\omega}}_{p,0}}  \frac{r}{r_0} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\nonumber\\
                     &= \frac{{\bar{\tau}}_{p,0}}{\pi} \frac{r}{r_0} F\left(\sqrt{ln\left(\frac{r}{r_0}\right)} \right)\end{aligned}$$ where ${\bar{\rho}}_0 = \frac{\lambda_0}{\pi r_0^2}$ and ${\bar{\omega}}_{p,0} = \sqrt{\frac{q Q_\text{tot}{\bar{\rho}}_0}{m \epsilon_0}} = \frac{2 \pi}{{\bar{\tau}}_{p,0}}$.

Spherically-symmetric Integral Solution
---------------------------------------

We solve the spherically-symmetric integral with an analogous approach. Define $\tilde{u} = \sqrt{1 - \frac{r_0}{\tilde{r} }}$ and solving for $\tilde{r}$ gives $\tilde{r} = \frac{r_0}{1-\tilde{u}^2}$. Thus $$\begin{aligned}
  d \tilde{u} &= \frac{1}{2} \frac{1}{\sqrt{1 - \frac{r_0}{\tilde{r} }}}\frac{r_0}{\tilde{r}^2 } d \tilde{r}\nonumber\\
                  &= \frac{1}{\sqrt{1 - \frac{r_0}{\tilde{r} }}} \frac{(1 - \tilde{u}^2)^2}{2 r_0} d\tilde{r}\nonumber\end{aligned}$$ Applying this change of coordinates to Eq. (\[eq:3D time integral subbed\]) with $u = \sqrt{1 - \frac{r_0}{r}}$, we get $$\begin{aligned}
  \text{sph: }t &= \sqrt{\frac{2 r_0}{aP_0}}\int_{0}^{u} \frac{2 r_0}{(1 - \tilde{u}^2)^2} d \tilde{u}\nonumber\\
                    &=  2 \sqrt{\frac{2 r_0^3}{aP_0}} \int_{0}^{u} \frac{1}{(1 - \tilde{u}^2)^2} d \tilde{u} \nonumber\\   
                    &=  2 \sqrt{\frac{2 r_0^3}{aP_0}} \left( \frac{1}{2} \tanh^{-1} \left(\tilde{u}\right) + \frac{1}{2 } \frac{\tilde{u}}{1 - \tilde{u}^2}\right) \bigg \rvert_{\tilde{u} = 0}^{\tilde{u} = \sqrt{1 - \frac{r_0}{r}}}\nonumber\\
                    &= \sqrt{\frac{2 r_0^3}{a P_0}} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right) + \frac{\sqrt{1 - \frac{r_0}{r}}}{1 - 1 + \frac{r_0}{r}}\right)\nonumber\\
                    &= \sqrt{\frac{2 \pi r_0^3 m \epsilon_0}{q Q_\text{tot} P_0}} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right) + \frac{r}{r_0}\sqrt{1 - \frac{r_0}{r}}\right)\nonumber\\
                    &= \sqrt{\frac{2}{3}} \sqrt{\frac{m \epsilon_0}{q Q_\text{tot} \frac{P_0}{\frac{4}{3}\pi r_0^3}}} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right)\right.\nonumber\\
                    &\quad\quad \left.+ \frac{r}{r_0}\sqrt{1 - \frac{r_0}{r}}\right)\nonumber\\
                    &= \sqrt{\frac{2}{3}} \sqrt{\frac{m \epsilon_0}{q Q_\text{tot} {\bar{\rho}}_0}} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right)\right.\nonumber\\
                    &\quad\quad \left. + \frac{r}{r_0}\sqrt{1 - \frac{r_0}{r}}\right)\nonumber\\
                    &= \sqrt{\frac{2}{3}} \frac{1}{{\bar{\omega}}_{p,0}} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right) + \frac{r}{r_0}\sqrt{1 - \frac{r_0}{r}}\right)\nonumber\\
                    &= \sqrt{\frac{2}{3}} \frac{{\bar{\tau}}_{p,0}}{2 \pi} \left( \tanh^{-1} \left(\sqrt{1 - \frac{r_0}{r}}\right) + \frac{r}{r_0}\sqrt{1 - \frac{r_0}{r}}\right)\label{eq:3D time vs r}\end{aligned}$$ where the solution to the integral was obtained with Mathematica’s online tool[@Mathematica] and where ${\bar{\rho}}_0 = \frac{P_0}{\frac{4}{3}\pi r_0^3}$ and ${\bar{\omega}}_{p,0} = \frac{q Q_\text{tot}{\bar{\rho}}_0}{m \epsilon_0} = \frac{2 \pi}{{\bar{\tau}}_{p,0}}$.

Derivation of Derivatives with Respect to Initial Position {#ap:spatial derivatives}
==========================================================

As noted in the main text, much of the physics of distribution evolution in our models is captured in the term $\frac{dr}{dr_0}$. The procedure to derive the expressions for this derivative is to take the derivative of Eq. (\[eq:2D time vs r\]) and Eq. (\[eq:3D time vs r\]). We do this mathematics here.

The Cylindrically-symmetric Derivative
--------------------------------------

We begin by re-writing $t$ from Eq. (\[eq:2D time vs r\]) as $$\begin{aligned}
  t &= 2  r \sqrt{\frac{1}{a \lambda_0}} F\left(y \right)\end{aligned}$$ where $y = \sqrt{ln\left(\frac{r}{r_0}\right)} = \sqrt{\ln \left(r \right) - \ln\left({r_0}\right)}$. So $$\begin{aligned}
  \frac{dy}{dr_0} &= \frac{1}{2y} \left(\frac{1}{r} \frac{d r}{d r_0} - \frac{1}{r_0} \right)\label{eq:derivative 2D sqrt term} 
 \end{aligned}$$ The Dawson function has the property $\frac{d}{dy} F(y) = 1 - 2 y F(y) = \left(\frac{1}{F(y)} - 2 y\right) F(y)$, and with the chain rule this becomes $\frac{d}{dr_0} F(y) = \left(\frac{1}{F(y)} - 2 y\right)  \frac{dy}{dr_0} F(y)$. Using Eq. (\[eq:derivative 2D sqrt term\]), this becomes $$\begin{aligned}
  \frac{d}{dr_0} F\left(y\right) &= \left(\frac{1}{F(y)} - 2 y\right)\frac{F(y)}{2y} \left(\frac{1}{r} \frac{d r}{d r_0} - \frac{1}{r_0} \right)\nonumber\\
                                              &=F(y) \left(\frac{1}{2 y F(y)} - 1\right)\left(\frac{1}{r} \frac{d r}{d r_0} - \frac{1}{r_0} \right)\label{eq:derivative our Dawson} \end{aligned}$$ Also, note $$\begin{aligned}
  \frac{d}{dr_0}\frac{1}{\sqrt{\lambda_0}} &= -\frac{1}{2}\frac{1}{(\lambda_0)^{3/2}} \frac{d\lambda_0}{d r_0}\nonumber\\
                                                                   &= -\frac{1}{2 \sqrt{\lambda_0}} \frac{d \ln(\lambda_0)}{d r_0}\label{eq:derivative sqrt a lambda_0}
 \end{aligned}$$ So $$\begin{aligned}
  0 &= \frac{dt}{dr_0}\nonumber\\
     &= \frac{t}{r} \frac{dr}{dr_0} - \frac{1}{2} t  \frac{d \ln(\lambda_0)}{d r_0} + t \left(\frac{1}{2 y F(y)} - 1\right)\left(\frac{1}{r} \frac{d r}{d r_0} - \frac{1}{r_0} \right)\nonumber\\
     &= \frac{t}{r} \frac{1}{2 y F(y)}\frac{dr}{dr_0} - t \left(\frac{1}{2} \frac{d \ln(\lambda_0)}{d r_0} + \left(\frac{1}{2 y F(y)} - 1\right) \frac{1}{r_0}\right)\end{aligned}$$ which gives $$\begin{aligned}
  \frac{dr}{dr_0} &= 2 y F(y) r \left(\frac{1}{2} \frac{d \ln(\lambda_0)}{d r_0} + \left(\frac{1}{2 y F(y)} - 1\right) \frac{1}{r_0}\right)\nonumber\\
                         &= \frac{r}{r_0}\left( 1 + \left(\frac{r_0}{2} \frac{d \ln(\lambda_0)}{d r_0} - 1 \right) 2 y F(y) \right)\nonumber\\
                         &= \frac{r}{r_0}\left( 1 + \left(\frac{r_0}{2} \frac{2 \pi r_0 \rho_0}{\lambda_0}  - 1 \right) 2 y F(y) \right)\nonumber\\
                         &= \frac{r}{r_0}\left( 1 + \left(\frac{\rho_0}{\frac{\lambda_0}{\pi r_0^2}}  - 1 \right) 2 y F(y) \right)\nonumber\\ 
                         &= \frac{r}{r_0} \left(1 + \left(\frac{\rho_0}{\bar{\rho}_0} - 1\right) 2  \sqrt{\ln\left(\frac{r}{r_0}\right)}F\left(\sqrt{\ln\left(\frac{r}{r_0}\right)}\right)\right) \end{aligned}$$

The Spherically-symmetric Derivatives
-------------------------------------

We begin by re-writing $t$ from Eq. (\[eq:3D time vs r\]) as $$\begin{aligned}
  t &= \sqrt{\frac{2 r_0^3}{a P_0}} \left(\tanh^{-1} y + \frac{r}{r_0} y \right)\end{aligned}$$ where $y = \sqrt{1 - \frac{r_0}{r}}$. So $$\begin{aligned}
  \frac{dy}{dr_0} &= \frac{1}{2y}\left(\frac{r_0}{r^2} \frac{dr}{dr_0} - \frac{1}{r}\right)\nonumber\\
                          &= - \frac{1}{2yr}\left(1- \frac{r_0}{r} \frac{dr}{dr_0}\right) \end{aligned}$$ Hence $$\begin{aligned}
  \frac{d\tanh^{-1} y}{dr_0} &= \frac{1}{1 - y^2} \frac{dy}{dr_0}\nonumber\\
                                        &= \frac{r}{r_0} \left(- \frac{1}{2yr}\left(1- \frac{r_0}{r} \frac{dr}{dr_0}\right)\right)\nonumber\\
                                        &= - \frac{1}{2yr_0}\left(1- \frac{r_0}{r} \frac{dr}{dr_0}\right)\end{aligned}$$ and $$\begin{aligned}
  \frac{d\left(\frac{r}{r_0} y \right)}{dr_0} &= \frac{y}{r_0} \frac{dr}{dr_0} - \frac{r}{r_0^2}y + \frac{r}{r_0} \frac{dy}{dr_0}\nonumber\\
                                                              &= \frac{1}{r_0}\left( \left(\frac{dr}{dr_0} - \frac{r}{r_0}\right)y -\frac{1}{2y}\left(1- \frac{r_0}{r} \frac{dr}{dr_0}\right)   \right)\nonumber\\
                                                              &= \frac{1}{2 y r_0}\left( 2 \left(\frac{dr}{dr_0} - \frac{r}{r_0}\right)\left(1 - \frac{r_0}{r}\right) - 1+ \frac{r_0}{r} \frac{dr}{dr_0}\right)\nonumber\\ \end{aligned}$$ Therefore $$\begin{aligned}
  \frac{d\left(\tanh^{-1} y + \frac{r}{r_0} y \right)}{dr_0} &= \frac{1}{y r_0}\left( \left(\frac{dr}{dr_0} - \frac{r}{r_0}\right)\left(1 - \frac{r_0}{r}\right)\right.\nonumber\\
                            &\quad\quad \left. - 1+ \frac{r_0}{r} \frac{dr}{dr_0}\right)\nonumber\\ 
                                                                                    &= \frac{1}{y r_0}\left( \frac{dr}{dr_0} - \frac{r}{r_0}\right)\nonumber\\  \end{aligned}$$ Also, similar to Eq. (\[eq:derivative sqrt a lambda\_0\]), $\frac{d}{dr_0}\frac{1}{\sqrt{P_0}} = -\frac{1}{2 \sqrt{P_0}} \frac{d \ln(P_0)}{d r_0}$. Putting this together we have $$\begin{aligned}
  0 &= \frac{dt}{dr_0}\nonumber\\
     &= \frac{3}{2} \frac{t}{r_0} - \frac{t}{2} \frac{d \ln(P_0)}{d r_0} +  \sqrt{\frac{2 r_0^3}{a P_0}} \frac{d\left(\tanh^{-1} y + \frac{r}{r_0} y \right)}{dr_0}\nonumber\\
     &= \frac{3}{2} \frac{t}{r_0} - \frac{t}{2} \frac{d \ln(P_0)}{d r_0} +  \sqrt{\frac{2 r_0^3}{a P_0}} \frac{1}{y r_0}\left( \frac{dr}{dr_0} - \frac{r}{r_0} \right)\end{aligned}$$ Solving for $\frac{dr}{dr_0}$ we get $$\begin{aligned}
  \frac{dr}{dr_0} &= \frac{r}{r_0} - \frac{3y}{2} \left(\tanh^{-1} y + \frac{r}{r_0} y \right)  \nonumber\\
                        &\quad\quad + \frac{y r_0}{2}\left(\tanh^{-1} y + \frac{r}{r_0} y \right)\frac{d \ln(P_0)}{d r_0}\nonumber\\
                         &= \frac{r}{r_0}  + \frac{3y}{2}\left( \frac{r_0}{3}\frac{d \ln(P_0)}{d r_0} - 1\right) \left(\tanh^{-1} y + \frac{r}{r_0} y \right)\nonumber\\
                         &=  \frac{r}{r_0} \left(1 + \frac{3}{2}\left(\frac{r_0}{3}\frac{4 \pi r_0^2 \rho_0}{P_0} - 1\right) \left(\frac{r_0}{r} y \tanh^{-1} y + y^2 \right)\right)\nonumber\\
                         &=  \frac{r}{r_0} \left(1 + \frac{3}{2}\left(\frac{\rho_0}{\frac{P_0}{\frac{4}{3}\pi r_0^3}} - 1\right)\left(\frac{r_0}{r} \sqrt{1 - \frac{r_0}{r}}  \tanh^{-1}\left(\sqrt{1 - \frac{r_0}{r}}\right)\right.\right.\nonumber\\
                         &\quad\quad \left.\left.  + 1 - \frac{r_0}{r} \right)\right)\nonumber\\
                         &= \frac{r}{r_0} \left(1 + \frac{3}{2}\left(\frac{\rho_0}{{\bar{\rho}}_0} - 1\right)\left(\frac{r_0}{r} \sqrt{1 - \frac{r_0}{r}}  \tanh^{-1}\left(\sqrt{1 - \frac{r_0}{r}}\right) \right.\right. \nonumber\\
                         &\quad\quad \left.\left. + 1 - \frac{r_0}{r} \right)\right)\end{aligned}$$
