Density shocks in the relativistic expansion of highly charged one component plasmas

TL;DR

This paper investigates how relativistic effects influence density shock formation in highly charged, one-component plasmas during expansion, revealing new behaviors and the enhancement of shocks under relativistic conditions.

Contribution

It demonstrates that relativistic coupling alters plasma expansion dynamics, leading to novel density shock phenomena not present in non-relativistic models.

Findings

Relativistic coupling causes new behaviors in plasma expansion.

Density shocks are enhanced under relativistic conditions.

A uniform cold bunch can develop shocks in relativistic regimes.

Abstract

In a previous paper we showed that dynamical density shocks occur in the non-relativistic expansion of dense single component plasmas relevant to ultrafast electron microscopy; and we showed that fluid models capture these effects accurately. We show that the non-relativistic decoupling of the relative and center of mass motions ceases to apply and this coupling leads to novel behavior in the relativistic dynamics under planar, cylindrical, and spherical symmetries. In cases where the relative motion of the bunch is relativistic, we show that a dynamical shock emerges even in the case of a uniform bunch with cold initial conditions; and that density shocks are in general enhanced when the relative motion becomes relativistic. Furthermore, we examine the effect of an extraction field on the relativistic dynamics of a planar symmetric bunch.

Equations

$$D_{01}=\frac{{\bar{\rho }}_{01}}{{\rho }_{01}}-1$$

$$D_{0d}=\frac{{\rho }_{0d}}{{\bar{\rho }}_{0d}}-1$$

$$E_0(z_0)=E_{01}=\frac{{\Sigma }_{tot}{P}_{01}}{2{ϵ}_0}=E_{T1}{P}_{01}$$

$$P_{01}=P_{01}(z_0)={\int }_0^{z_0}2{\rho }_{01}(z)dz$$

$$E_0(r_0)=E_{02}=\frac{{\Lambda }_{tot}{P}_{02}}{2\pi {ϵ}_0{r}_0}=E_{T2}{P}_{02}$$

$$P_{02}={\int }_0^{r_0}2\pi r{\rho }_{02}(r)dr$$

$$E_0(r_0)=E_{03}=\frac{Q_{tot}{P}_{03}}{4\pi {ϵ}_0{r}_0^2}=E_{T3}{P}_{03}$$

$$P_{03}={\int }_0^{r_0}4\pi r^2{\rho }_{03}(r)dr$$

$$\frac{dp}{dt}=qE$$

$$p=\gamma mv$$

$$\gamma =\frac{1}{\sqrt{1-(\frac{v}{c}{)}^2}}$$

$$E_1(z)=E_{01}$$

$$E_2(r)=E_{02}\frac{r_0}{r}$$

$$E_3(r)=E_{03}{\left(\frac{r_0}{r}\right)}^2$$

$$K=(\gamma -1)m{c}^2$$

$$\Delta U_1=q{E}_{01}(z_0-z)$$

$$\Delta U_2=q{E}_{02}{r}_0\ln \left(\frac{r_0}{r}\right)$$

$$\Delta U_3=q{E}_{03}{r}_0^2\left(\frac{1}{r}-\frac{1}{r_0}\right)$$

$$(\gamma -1)m{c}^2=-\Delta U$$

$$v=\frac{\sqrt{{\left(1-\frac{\Delta U}{m{c}^2}\right)}^2-1}}{1-\frac{\Delta U}{m{c}^2}}c.$$

$$t-t_0={\int }_{y_0}^y\frac{1}{v}dy$$

$${\rho }_{01}d{z}_0={\rho }_1(z;t)dz$$

$${\rho }_{02}2\pi r_{0}d{r}_0={\rho }_2(r;t)2\pi rdr$$

$${\rho }_{03}4\pi r_0^{2}d{r}_0={\rho }_3(r;t)4\pi r^{2}dr$$

$${\rho }_d(y;t)$$

$${\omega }_p=\sqrt{\frac{q^{2}n}{{ϵ}_{0}m}}$$

$${\omega }_{01}=\sqrt{\frac{q{\Sigma }_{tot}{\overline{\rho }}_{01}}{{ϵ}_{0}m}}=\sqrt{\frac{q{E}_{01}}{m{z}_0}},$$

$${\omega }_{02}=\sqrt{\frac{q{\Lambda }_{tot}{\overline{\rho }}_{02}}{{ϵ}_{0}m}}=\sqrt{\frac{2q{E}_{02}}{m{r}_0}},$$

$${\omega }_{03}=\sqrt{\frac{q{Q}_{tot}{\bar{\rho }}_{03}}{{ϵ}_{0}m}}=\sqrt{\frac{3q{E}_{03}}{m{r}_0}}.$$

$$\frac{{\omega }_{0d}^2}{c^2}=\frac{d}{y}\frac{q{E}_{0d}}{m{c}^2}$$

Full text

I Introduction

The expansion dynamics of highly charged plasmas is a fundamental problem in areas ranging from astrophysics to nanotechnology to beam physics. Previous analytic work has focused on initial conditions where a highly charged plasma is cold and has uniform densityJansen (1988); Reiser (1994); Batygin (2001); Bychenkov and Kovalev (2005); Grech et al. (2011); Kaplan et al. (2003); Kovalev and Bychenkov (2005); Last et al. (1997); Eloy et al. (2001); Krainov and Roshchupkin (2001); Morrison and Grant (2015); Boella et al. (2016); Bychenkov and Kovalev (2011) However, the vast majority of this work, with the exception of Bynchenkov and KovalevBychenkov and Kovalev (2011), have assumed non-relativistic conditions. In ultrafast electron microscopy (UEM) and some beam physics applications, electron sources are used to produce dense bunches of charged particles within an intense extraction field that is used to accelerate the distribution to near-luminal speeds. In addition, for sufficient densities, the bunch self-field can result in relativistic velocities within the frame of the bunch. These concerns indicate that a relativistic theory is required in many practical cases, and here we present the relevant theory.

The typical analytic approach to relativistic expansion dynamics of such systems is a treatment based on envelope equations that are predicated on the conservation of emittance and the use of uniform density distributionsReiser (1994). Uniform ellipsoidal distributions are particularly amenable to analysis as the self-electrostatic field in these distributions is linear and the expansion dynamics results in a simple power law growth of the ellipsoid axes – at least in the non-relativistic regime. Furthermore, it is fairly straightforward to show that uniform ellipsoids conserve emittance as long as all the particles can be treated as having identical Lorentz factors. Moreover, analysis of beam dynamics such as emittance oscillationAnderson (1987), emittance compensationRosenzweig et al. (2006), and the beam halo Gluckstern (1994) generally assume similar uniform-like conditions. However, for electron injectors utilizing photoemission the initial conditions of the bunch is often Gaussian, or at the very least non-uniform, and it has long been known that charge redistribution from the non-uniform to the uniform bunch is one of the major sources of emittance growthWangler (1991) suggesting that uniform distributions are at best an idealization that miss much of the physics present in the typical situation. However, we have recently shown that dynamics similar to Wangler’s charge redistribution for freely expanding bunches leads to an opportunity of “Coulomb cooling” — a mechanism we believe employs the intense Coulomb fields to carry off heat from non-neutral plasmasWilliams et al. (2017); Zerbe et al. (2018). Our analysis is directed at a better understanding of the relativistic expansion dynamics of uniform and non-uniform systems to study mechanisms of emittance growth near the particle source but also with the goal of understanding and optimizing Coulomb cooling. Here we concentrate on characterizing the relativistic density dynamics for both uniform and non-uniform initial conditions.

Numerous works within the UEM literature have already looked at various aspects of the evolution of non-uniform distributionsDegtyareva et al. (1998); Luiten et al. (2004); Musumeci et al. (2008); Morrison et al. (2013); Li and Lewellen (2008); Siwick et al. (2002); Qian and Elsayed-Ali (2002); Reed (2006); Collin et al. (2005); Gahlmann et al. (2008); Tao et al. (2012); Portman et al. (2013, 2014); Michalik and Sipe (2006); Zerbe et al. (2018). Reed presented a fluid model that described the dynamics of non-uniform bunch expansion under the assumptions that the bunch could be treated as having planar symmetry in the non-relativistic regimeReed (2006). In our previous paperZerbe et al. (2018), we showed that the expansion dynamics in the planar case differs greatly from symmetries in higher dimensions and verified analytic description of the dynamics utilizing both N-particle and particle-in-cell (PIC) simulations. The analytic descriptions depend on terms that can be written as functions multiplied by the quantity

[TABLE]

for planar geometries or by the quantity

[TABLE]

for cylindrical ($d=2$) and spherical ($d=3$) geometries, where $\rho_{0d}$ is the initial density at location $r_{0}$ and ${\bar{\rho}}_{0d}$ is the average density within that location. For the uniform distribution, the local density and the average density are the same so that $D_{0d}=0$ everywhere and the density evolution’s dependence on the aforementioned functions vanishes reducing the dynamics to the uniform dynamics utilized extensively in the literature. Otherwise, these functions play a large role in the density evolution leading to differences in the dynamics of distributions with different initial density profiles. Specifically, we showed that a density shock, seen in Coulomb explosion studiesGrech et al. (2011); Kaplan et al. (2003); Kovalev and Bychenkov (2005); Last et al. (1997); Murphy et al. (2014); Reed (2006); Degtyareva et al. (1998), was present in the analytic density evolution of initially Gaussian distributions under both cylindrical and spherical symmetries but absent under planar symmetry unless an appropriate initial chirp in phase space is present. Moreover, such a shock is absent in non-relativistic dynamics of uniform systems with cold initial conditions suggesting that uniform distribution evolution is unique in this regardZerbe et al. (2018).

We noted in our previous work that the density evolution seen in PIC simulations using an electromagnetic (EM) solver and relativistic particle pusher, which should capture all relativistic effects if the initial field is accurate, do not significantly differ from our analytic description for the densities analyzed thereZerbe et al. (2018). Moreover, we have seen that the PIC simulations using an EM solver do not significantly differ from those using an electrostatic (ES) solver with relativistic particle pusher for much higher densities than those examined in our previous work, or in the work described here. Assuming the accuracy of the initial fields, we conclude that the relativistic effects are then adequately captured within the relativistic particle pusher, which is equivalent to simply including the relativistic momentum in the analysis. Precisely such an analysis of the relativistic free-expansion of a spherically-symmetric, cold uniform charge distribution was completed by Bychenkov and KovalevBychenkov and Kovalev (2011), and part of what we do in this manuscript is extend this analysis to non-uniform cases as well as additional symmetries.

Here we treat charge distributions with general initial spatial distributions starting from rest under planar, cylindrical and spherical geometries introducing a novel length scale that is associated with each symmetry. First in Section II.1, we present general results applicable to all cases. In Section II.2 we derive expressions for planar symmetry for any arbitrary initial spatial distribution and examine these expressions in the non-relativistic and highly relativistic limits. We then introduce $M$-shell simulations, which are simulations of $M$ equally charged planes in 1D, and show that these simulations reproduce the density evolution derived analytically. In Sections II.3 and II.4, we derive relativistic density evolution expressions under cylindrical and spherical symmetries, respectively, for arbitrary initial distributions, and we show that these expressions are consistent with PIC calculations using an EM solver and relativistic particle pusher utilizing the well-known package warpFriedman et al. (2014). Further, we show that $M$-shell simulations, which track $M$ equally charged cylindrical- and spherical- shells in 2- and 3- dimensions, respectively, also capture the same density evolution. We validate these expressions against their non-relativistic and uniform relativistic counterparts, and we examine the expressions in the highly relativistic limit.

In Section III we introduce an external extraction field in the case of a planar electron bunch, and we point out that the acceleration from self-fields and external fields do not decouple in the relativistic case. Though the analysis is captured by a straightforward extension of the analysis used for the planar case with no extraction field, the physical effects are quite interesting and relevant to ultrafast electron microscopy. In Section IV, we demonstrate how to apply the time dependent distributions to calculate the statistical width of the distribution. Section V contains discussions and conclusions; noting in particular the emergence of density shocks due to relativistic effects, even in planar uniform systems where a shock does not emerge in the non-relativistic limit.

II Relativistic density evolution

In this section we consider planar, cylindrical and spherical geometries in the case of cold initial conditions and where there is no external electromagnetic field. We treat the dynamics using the relativistic treatment of momentum and energy, however we treat the forces between electrons using electrostatics. We then check the latter approximation under cylindrical and spherical geometries using PIC calculations using a full EM solver and find excellent agreement. We start with a general analysis and then specialize in the later three subsections.

II.1 General considerations

II.1.1 General formulation

We consider cold symmetric initial charge distributions so that the electric field at position $z_{0}$ in a planar geometry is given by,

[TABLE]

where $\Sigma_{tot}=Nq/A$, with $A$ the pulse area, $N$ the number of particles in the bunch, $q$ the charge of each particle, $E_{T1}=\frac{\Sigma_{tot}}{2\epsilon_{0}}$ is the total planar field produced by the particles, and $P_{01}$ is a position specific scalar representing the fraction of the electrons in the bunch contributing to the net electric field on the particle at position $z_{0}$. The physical significance of the total field $E_{T1}$ will be discussed in detail in Section III. The direction along which the charge density varies is $z$, and we assume the initial charge distribution is symmetric about the origin for the sake of simplicity. We define $\rho_{1}(z;t)$ to be the number density per unit length at position $z$ and $\rho_{01}=\rho_{01}(z_{0})=\rho_{1}(z_{0};t=0)$ to be the initial number density at position $z_{0}$. As the initial density $\rho_{01}$ is symmetric about the origin, the cumulative probability is given by,

[TABLE]

Notice that the quantity $\Sigma_{tot}P_{01}$ represents the charge per unit area after integrating the number density per unit length $\rho_{01}$ over the range $[-z_{0},z_{0}]$.

In systems with cylindrical symmetry, we have,

[TABLE]

where $\Lambda_{tot}$ is the total charge per unit length along the cylindrical charge distribution, $E_{T2}=\frac{\Lambda_{tot}}{2\pi\epsilon_{0}r_{0}}$ is the electric field a particle at $r_{0}$ would feel if the entire distribution were distributed cylindrically within $r_{0}$, and $P_{02}$ is the cumulative probability given by,

[TABLE]

where $\rho_{2}(r;t)$ is the number density in two dimensions (number per unit area) and where we define $\rho_{02}=\rho_{02}(r_{0})=\rho_{2}(r_{0};t=0)$. Notice that the quantity $\Lambda_{0}P_{02}$ represents the charge per unit length inside radius $r_{0}$.

In systems with spherical symmetry we have

[TABLE]

where $Q_{tot}$ is the total charge in the system, $E_{T3}=\frac{Q_{tot}}{4\pi\epsilon_{0}r_{0}^{2}}$ and is the electric field a particle at $r_{0}$ would feel if the entire distribution were distributed spherically within $r_{0}$, and $P_{03}$ is the cumulative probability given by,

[TABLE]

where $\rho_{3}(r;t)$ is the number density in three dimensions (number per unit volume) and where we define $\rho_{03}=\rho_{03}(r_{0})=\rho_{3}(r_{0};t=0)$. Again notice that $P_{03}$ represents the fraction of the particles that lie inside radius $r_{0}$ and $Q_{tot}P_{03}$ gives the charge inside radius $r_{0}$.

To make analytic progress we make the laminar fluid approximation, which states that there is no mixing of the charged particle trajectories. As a result, the symmetries of the charge distributions are conserved. If we consider a particle of charge $q$ and rest mass $m$ starting from rest (cold initial conditions); at position $z_{0}$ (planar case) or $r_{0}$ (cylindrical or spherical cases), we want to find the position and velocity of the particle at later times. In a planar system we thus want to find $z(z_{0},t)=z$ and $v(z_{0},t)=v$; while in the cylindrical and spherical cases we want to find the radial position and radial velocity $r(r_{0},t)=r$; $v(r_{0},t)=v$. One approach is to use the relativistic form of Newton’s second law,

[TABLE]

where we are considering the $z$ component of these vectors in planar geometries and the radial component in cylindrical and spherical geometries. We use the relativistic expression for the momentum,

[TABLE]

where $m$ is the rest mass and

[TABLE]

Due to the laminar fluid property, the electric field experienced by a particle at position $z(z_{0},t)$ in planar geometries is the same as the electric field this particle experienced at it’s initial position, so that

[TABLE]

In cylindrical geometries using Gauss’ law we find,

[TABLE]

and analogously in spherical geometries we have,

[TABLE]

These expressions may be used with Newton’s second law to solve for the particle dynamics. Alternatively in an energy formulation, conservation of energy requires that the change in kinetic energy equals the change in potential energy. We use the relativistic form of the kinetic energy,

[TABLE]

and the change in kinetic energy is equal to $K$ as the bunch starts from rest. The change in potential energy is found by integrating the force $qE$, which for the planar case gives,

[TABLE]

while for the cylindrical case we find,

[TABLE]

and for the spherical case,

[TABLE]

Notice that $\Delta U\leq 0$ for all time due to the physics as the potential energy decreases as the electron bunch expands. Setting the sum of the changes in potential and kinetic energy to zero, we have,

[TABLE]

where the appropriate expression for $\Delta U$ must be utilized. From this expression a general relation between the velocity and position is found to be,

[TABLE]

Moreover, since $\Delta U$ only depends on the position, this equation may be integrated to find an expression relating time and position,

[TABLE]

where $y=z$ for planar cases and $y=r$ for cylindrical and spherical cases.

Finally to obtain an expression for the time evolution of the density, we use the conservation of the charge density under laminar conditions. This conservation can be stated for planar systems as

[TABLE]

for cylindrical geometries as

[TABLE]

and for spherical systems as

[TABLE]

In general, this results in the relationship between the density and the initial density of

[TABLE]

where again $y=z$ for planar cases and $y=r$ for cylindrical and spherical cases, $d$ is $1$, $2$, and $3$ for these symmetries, respectively, and ${}^{\prime}\equiv\frac{d}{dy_{0}}$ with the $d$’s in this last expression representing differentiation – not dimensionality of the problem.

II.1.2 Fundamental parameters

Some fundamental parameters need to be considered in the discussion of high density single component plasmas. The first is the plasma frequency,

[TABLE]

which describes the frequency of coherent plasma oscillations, where $n$ is the number density (number of particles per unit volume) and $q$ the particle charge. We note that relativistic effects affect the plasma frequency, but as our distribution is starting from rest, it is sufficient to consider the non-relativistic plasma frequency; however, the plasma frequency for different symmetries is not apparent from Eq. (26). We define average initial densities $\overline{\rho}_{01}=\frac{P_{01}}{2z_{0}}$, $\overline{\rho}_{02}=\frac{P_{02}}{\pi r_{0}^{2}}$, $\overline{\rho}_{03}=\frac{P_{03}}{\frac{4}{3}\pi r_{0}^{3}}$ which are the average densities inside distance $z_{0}$ (planar case), or inside radius $r_{0}$ for the cylindrical and spherical cases. These definitions are used to define initial plasma frequencies as follows for planar systems,

[TABLE]

cylindrical systems,

[TABLE]

and spherical systems

[TABLE]

This can be summarized by

[TABLE]

for $d\in\{1,2,3\}$ and where $y$ is $z_{0}$ when $d=1$ and $r_{0}$ for $d=2$ and $d=3$.

As will be seen below, the time $\tau_{0d}$ defined as,

[TABLE]

sets the timescale for the relativistic expansion of high density charge clouds; as was found in the non-relativistic casesZerbe et al. (2018).

In addition to the plasma frequency, we find it advantageous to define the related 1D-number density as

[TABLE]

where $\rho_{r0d}$ has units of inverse length. We call $\rho_{r0d}$ the relativistic crossover density for planar, cylindrical, and spherical symmetries for $d=1,2,3$, respectively. The physical interpretation of the relativistic crossover density is that it provides a scale for the potential energy as $\frac{\Delta U_{1}}{mc^{2}}=\rho_{r01}(z_{0}-z)$, $\frac{\Delta U_{2}}{mc^{2}}=\rho_{r02}r_{0}\ln\left(\frac{r_{0}}{r}\right)$, and $\frac{\Delta U_{3}}{mc^{2}}=\rho_{r03}r_{0}\left(\frac{r_{0}}{r}-1\right)$. The relativistic crossover density is related to the plasma frequency through

[TABLE]

The relativistic length scale, $l_{r0d}$ is related to the relativistic density through,

[TABLE]

where $l_{r0d}$ is seen to be independent of the initial distribution. $l_{r0d}$ can be thought of as the distance a particle experiencing the force obtained by the full distribution at the given coordinate needs to travel before having kinetic energy of $mc^{2}$. Notice, that $l_{r01}$ is a constant and is specifically independent of $z_{0}$; however, $l_{r02}\propto r_{0}$ and $l_{r01}\propto r_{0}^{2}$.

II.2 Planar symmetry

In this case, Eq. (20) becomes

[TABLE]

where $\rho_{r01}$ is from Eq. (33). The integral in Eq. (21) may be carried out to find,

[TABLE]

which can be inverted to find $z(z_{0},t)$ as

[TABLE]

where

[TABLE]

Taking the time derivative of Eq. (37), the velocity as a function of time becomes,

[TABLE]

From Eq. (25), we find the density dynamics,

[TABLE]

where

[TABLE]

The non-relativistic limit occurs when $\rho_{r01}ct<<1$ or equivalently when $t<<t_{x}$ where

[TABLE]

and in this limit the expressions above reduce to the known results, i.e. $z=z_{0}+qE_{01}t^{2}/2m$ and

[TABLE]

where $\rho_{1NR}(z;t)$ is the density in the non-relativistic limit and $\omega_{01}$ is the plasma frequency defined in Eq. (27). Reed (2006); Zerbe et al. (2018).

The highly relativistic limit is when $\rho_{r01}ct>>1$ or equivalently $t>>t_{x}$. Note that this second inequality implies that any point in the distribution except the center point at $z_{0}=0$ becomes highly relativistic for sufficient time; this is part of the nature of the planar symmetries, and we find similar nature for the cylindrical symmetries below. In this limit, we find,

[TABLE]

where the sign of the luminal velocity is determined by on which side of $z_{0}=0$ the particle originated and

[TABLE]

or equivalently

[TABLE]

where $\rho_{1HR}(z;t)$ is the density distribution in the highly relativistic limit. The interpretation of this result is interesting. First, the majority of the distribution essentially becomes two pulses traveling at near luminal speeds away from one another. Second, as the particles within the distribution reach luminal speeds, the density no longer significantly changes as the particles propagate to the left or right; that is, the density evolves toward an “asymptotic density” determined by Eq. (45). If $\rho_{r01}<<\rho_{01}$, then $\rho_{HR}\to\frac{(P_{01})^{2}}{l_{r01}}$; however, if $\rho_{01}>>\rho_{r01}$, then on the edges $\rho_{HR}\to\rho_{01}$ whereas as you go further in the distribution transitions to $\frac{(P_{01})^{2}}{l_{r01}}$. This behavior for the uniform and Gaussian distributions for various ratios of $\frac{l_{r01}}{L_{0}}$, where $L_{0}$ indicates the original width, may be seen in Fig. 1.

Analytically, for the case of an initial uniform distribution, $\rho_{01}=\frac{1}{L_{0}}$ and $P_{0}=\frac{2z_{0}}{L_{0}}=2z_{0}\rho_{01}$ where $L_{0}$ is the initial width of the distribution. In this case,

[TABLE]

Thus, the shape of this asymptotic distribution is determined entirely by the length scale, $l_{r01}$, and the initial width, $L_{0}$. For any point $z_{0}<<l_{r01}$, including the entire distribution if $L_{0}<<l_{r01}$, this asymptotic density is essentially parabolic with zero density at the center and $\frac{1}{l_{r01}}$ at the edge. This case can be seen in Fig. 3a. For extremely dense distributions where $L_{0}>>l_{r01}$, the asymptotic density at the edges approaches the original density, $\rho_{01}$. There is also a period of transition between the parabolic and original density when the length scale is much smaller than the original width. Both asymptotic behaviors can be seen in Fig. 1 for both the uniform and Gaussian cases.

The mechanism for the relativistic peak emergence may be further seen in Fig. 2 where the density associated with different Lagrangian particles within the distribution are tracked and shown to asymptote to various density values predicted by Eq. (45). One way to describe this mechanism is to notice that all particles, excepting the center particle, in a planar model will asymptote to the speed of light. As the density is physically smooth, the particles’ velocities in the neighborhood of the Lagrangian particle asymptote similarly to the speed of light. In other words, the relative velocity of the particles in the Lagrangian particle’s neighborhood asymptotes to zero, and the particles cease to spread in the $z$ dimension. As the $z$ dimension is the only dimension in which the density is spreading in the planar model, this is the same as freezing the density to a constant value – an asymptote. Moreover, as particles toward the edge of the distribution have larger accelerations, these particles asymptote earlier than particles farther in. These differences in “freezing” time result in the middle of the distribution expanding, and becoming less dense, before the onset of the relativistic regime. Coupled with the initial distribution, this results in the formation of density peaks toward the edge of the distribution, as is seen in both the uniform and Gaussian distributions in Fig. 3. In the non-relativistic limit, there is no Coulomb shock in planar bunches with cold initial conditions; while in the relativistic limit a strong shock emerges and an initial bunch described by either uniform or a Gaussian density distribution evolves to a two peak structure with one bunch moving to the right and the other to the left (see Fig. 3).

Also apparent in Fig. 3 is the fact that stochastic effects are initially strong in simulated density profiles. However, at long times the theoretical density and simulated density agree well. This is a real effect. Specifically, consider the inter-particle distance between the $i^{th}$ and $(i+1)^{th}$ shells denoted as $d_{i}$. For a uniform distribution, order statistics tells us that $d_{i}(0)=\frac{L}{M+1}+\epsilon$ where $L$ is the total width of the distribution, $M$ is the total number of shells, and $\epsilon$ is a stochastic factor roughly of the size $\frac{L}{M+1}$. Thus, due to stochastic fluctuations, we’d expect some sheets to be bunched together giving a higher local density than the average and likewise other sheets to be further apart giving a lower local density than the average. This is precisely what is seen with the initial distribution in Fig. 3. However, as these sheets evolve, the relative non-relativistic acceleration is $\frac{2qE_{T1}}{mM}$, so $d_{i}(t)=d_{i}(0)+\frac{qE_{T1}}{mM}t^{2}$. Given sufficient time, $t>>\sqrt{\frac{mL}{qE_{T1}}}$, $d_{i}(t)\approx\frac{qE_{T1}}{mM}t^{2}$. That is, the inter-particle distance (and hence the distribution) is dominated by the space-charge effect and converges to the space-charge predicted distribution everywhere. Of course, if the bunch enters the relativistic regime prior to this smoothing, the stochastic effects will be preserved. We will see such behavior once we add an extraction field, but such behavior requires extremely dense bunches that may not be physically possible in free expansion experiments.

II.3 Cylindrical symmetry

Now we consider the expansion of an initially cold charged particle cloud with cylindrical symmetry. In this case, Eq. (20) becomes

[TABLE]

where $y^{2}=\ln\left(\frac{r}{r_{0}}\right)$, $\zeta^{2}=\frac{r_{0}\rho_{r02}}{2}=\frac{r_{0}P_{02}}{2l_{r02}}=\frac{r_{0}^{2}\omega_{02}^{2}}{4c^{2}}$, with $\rho_{r02}$ coming from Eq. (33). As $l_{r02}\propto r_{0}$, it should be apparent that $\zeta$’s dependence on $r_{0}$ is completed determined by $P_{02}(r_{0})$.

From Eq. (48) and (21), we find the implicit relation between time and radial position through the integral,

[TABLE]

To make the connection with previous work, we introduce a generalized Dawson function $\mathcal{F}$ through the definition,

[TABLE]

where $\zeta$ can be written as a function of $x$. Thus the time-spatial relation may be expressed as

[TABLE]

where

[TABLE]

When $g(\zeta,y)=1$, we reproduce the Dawson function, $F(x)=\mathcal{F}(1,x)$. Specifically when we are in the non-relativistic regime, we have $2\zeta y<<1$ and $g(\zeta,y)\approx 1$, so Eq. (51) reduces to

[TABLE]

which is the result we derived previously in the non-relativistic case. We can write down the derivative of the generalized Dawson function by applying the Leibniz rule

[TABLE]

Note that in the non-relativistic limit, $g(\zeta,y)=1$, and Eq. (54) reduces to the normal Dawson function derivative $\frac{dF}{dx}=-2xF(x)+1$.

Following the same reasoning as our previous workZerbe et al. (2018), we can obtain an analytic form for the time dependent density, i.e. the density evolution expression (see Eq. (25)). Evaluating $r^{\prime}=\frac{dr}{dr_{0}}$ by taking a derivative of Eq. (51) with respect to $r_{0}$, we find,

[TABLE]

where $\mathcal{F}$ is shorthand for $\mathcal{F}(g(\zeta,y),y)$ $\mathcal{F}_{\partial}$ is shorthand for $\mathcal{F}\left(\frac{\partial g}{\partial\zeta},y\right)$, and $D_{02}$ is from Eq. (2). Note $D_{02}$ measures the deviation from a uniform cylindrically-symmetric distribution, and for the uniform cylindrically-symmetric distribution case it is zero for all values of $r_{0}$ where $\rho_{0}$ is not [math].

From the above analysis, the density evolution is found to be,

[TABLE]

In Fig. 4, we compare the predictions of Eq. (56) to simulations for both uniform and Gaussian initial distributions. We choose the initial radius and radial standard deviation, respectively, to be 1 cm for $N=1\times 10^{13}$ electrons/cm. We again simulate with Warp using the EM solver as well as the 2D version of $M$-shell simulations. For the $M$-shell simulations, the initial radius of the $M$ cylindrical shells are sampled and then evolved according to Eq. (48) and Eq. (53) but with $\omega_{02}$ replaced by $\sqrt{\frac{3q\Lambda_{s}}{\pi r_{s,0}^{2}m\epsilon_{0}}}$ where $\Lambda_{s}$ is the charge per unit length contained in the cylindrical shell and $r_{s,0}$ is the initial radius of the shell. As can be seen in Fig. 4, the theory and both simulations agree on the evolution of both the uniform and non-uniform initial distributions. Similar to the planar case, the initial variance about the predicted value can be seen to decrease as the simulations evolve. Again, this indicates that the inter-shell distances are dominated by the space-charge effects resulting in the later simulations having less statistical variation from the expected distribution.

In the non-relativistic regime, $2\zeta y<<1$; $\mathcal{F}\to F$ and $\zeta\mathcal{F}_{\partial}\to 0$. Thus Eq. (56) reduces to

[TABLE]

which is the expression we found in our earlier, non-relativistic workZerbe et al. (2018).

Similar to the planar symmetric case, we are interested in the density the distribution should evolve under specific limits. We were unable to analytically obtain a limit analogous to the limit we found under planar symmetry in Eq. (45) as doing so requires evaluating the value of the modified Dawson function as $\ln\left(\frac{r}{r_{0}}\right)$ goes to infinity. We were able to see the freezing of the dimension in the extremely dense limit where $\zeta>>1$ as we have shown in Appendix C . Specifically, the evolution at the edge of the distribution can be approximated by $\rho_{2}(r;t)=\frac{r_{0}}{r}\rho_{0}$, which is the evolution of the uniform distribution under non-relativistic conditions in one-dimension lower, i.e. 1D. This situation is analogous to the high density 1D case that causes the edges to essentially immediately become relativistic likewise resulting in evolution of the uniform distribution under non-relativistic conditions in one-dimension lower, i.e. 0D or constant. However, this condition, $\zeta>>1$, is analogous to the 1D case when the entire distribution is essentially in the highly relativistic limit. We will shortly show that even in this case, the spherically symmetric evolution can be shown to freeze out a dimension; however, we believe that this freezing happens for cylindrically symmetric distributions regardless of the size of $\zeta$.

II.4 Spherical symmetry

Now we consider the expansion of an initially cold charged particle cloud with spherical symmetry. In this case, Eq. (20) becomes

[TABLE]

where $x^{2}=1-\frac{r_{0}}{r}$, $g_{1}(x)=1+\zeta^{2}x^{2}$, $g_{2}(x)=1+2\zeta^{2}x^{2}$, $\zeta^{2}=\frac{r_{0}\rho_{r03}}{2}=\frac{r_{0}P_{03}}{2l_{r03}}=\frac{r_{0}^{2}\omega_{03}^{2}}{6c^{2}}$, and $\rho_{r03}$ is from Eq. (33). As $l_{r03}\propto r_{0}^{2}$, it should be apparent that $\zeta\propto r_{0}P_{03}(r_{0})$.

From Eq. (48) and (21), we find the implicit relation between time and radial position through the integral,

[TABLE]

where $T(x)=\tanh^{-1}\left(\sqrt{\frac{g_{1}(1)}{g_{1}(x)}}x\right)$. Note that the $1$ inside the $g$ functions corresponds to $x$ at infinitely long times, i.e. $\lim_{\frac{r}{r_{0}}\to\infty}x=1$, so $g_{1}(1)=1+\zeta^{2}$ and $g_{2}(1)=1+2\zeta^{2}$. This expression is essentially the same expression as derived by Bychenkov and Kovalev, who first derived it for the case of uniform initial density distributions Bychenkov and Kovalev (2011). Our expression differs only in the interpretation of $\omega_{03}$ as ours can be dependent on $r_{0}$ whereas their $\omega_{03}$ is a constant, which is the correct interpretation for the uniform distribution. This difference in interpretation allows us to treat general initial distributions but requires additional consideration when determining the derivative of Eq. (59) with respect to $r_{0}$ as $\omega_{03}^{\prime}=\frac{3\omega_{03}}{2r_{0}}D_{03}$, where ${}^{\prime}\equiv\frac{d}{dr_{0}}$, with the $d$’s in this last expression representing differentiation – not dimensionality of the problem, and $D_{03}$ is from Eq. (2).

We follow the same reasoning as our previous work Zerbe et al. (2018) in order to obtain the density evolution expression. After taking the derivative of Eq. (59) with respect to $r_{0}$, we can solve for $r^{\prime}$ giving

[TABLE]

where

[TABLE]

and

[TABLE]

Plugging Eq. (60) into Eq. (25) we obtain the evolution of the density distribution

[TABLE]

In Fig. 5, we compare the prediction of Eq. (63) to simulations for both a uniform and Gaussian initial distributions of $N=1\times 10^{13}$ electrons. The simulation have $R=1$ cm (uniform) or $\sigma_{r}=1$ cm (Gaussian). We again simulate with Warp using the EM solver and as well as the 3D version of $M$-shell simulations. For the $M$-shell simulations, the initial radius of the $M$ spherical shells are sampled and then evolved according to Eq. (58) and Eq. (59) but with $\omega_{03}$ replaced by $\sqrt{\frac{3qQ_{s}}{4\pi r_{s,0}^{3}m\epsilon_{0}}}$ where $Q_{s}$ is the charge contained in the shell and $r_{s,0}$ is the initial sampled radius of the shell. As can be seen in Fig. 5, the theory captures the evolution of both the uniform and non-uniform initial distributions. Similar to both the planar and cylindrical cases, the initial variance around the theoretical value primarily seen in the uniform distribution decreases as the distribution expands.

For further validation, we compare Eq. (63) to the expression derived by Bychenkov and Kovalev. Their expression detailed the relativistic density evolution for the uniform distributionBychenkov and Kovalev (2011), $\rho_{unif}(r;t)$, which should be equivalent to our expression when $D_{03}=0$. In this case, Eq. (63) reduces to

[TABLE]

This expression for the density evolution for uniform initial conditions is identical to the expression published in the English translation of Bychenkov and Kovalev except for an obvious typo in that workBychenkov and Kovalev (2011).

Next, we compare this expression to our previous, non-relativistic expression. In the non-relativistic regime $2\zeta^{2}<<1$. Unlike the planar and cylindrical cases, the spherical model need never enter the relativistic regime and therefore this model may be relevant for all time. In this non-relativistic regime, Eq. (63) reduces to

[TABLE]

which is identical to the non-relativistic expression we previously derived but with $D_{03}=D$ in our previous notationZerbe et al. (2018).

Again we would like to analyze specific limits of the density evolution; fortunately, under spherical symmetry we can analyze the long time limit. In Appendix B we show that

[TABLE]

where $\rho_{x3}(r_{0})=\frac{1+3\zeta^{2}+2\zeta^{4}}{1+\frac{3}{2}D_{03}}\rho_{03}$ is entirely determined by the initial conditions. Notice, the pre-factor in $\frac{1+3\zeta^{2}+2\zeta^{4}}{1+\frac{3}{2}D_{03}}$ is essentially $1$ in the center where $\zeta\approx 0$ and $D_{03}\approx 0$, but that this value increases as $r_{0}$ increases. The time evolution of $\rho_{x3}$ and the predicted asymptote for this quantity can be seen in Fig. 6. For the uniform distribution, the increase in $\rho_{x3}$ as a function of $r_{0}$ is quartic as $D_{03}=0$ for all values of $r_{0}$. In real distributions, though, there should be a value for $r_{0}$ where $D_{03}=-\frac{2}{3}$, and we see that $\rho_{x3}$ has a zero in the denominator. This violates the assumptions made in the derivation of $\rho_{x3}$, and inspection of Appendix B shows that $r^{\prime}$ becomes [math] in the locality of $D_{03}=-\frac{2}{3}$ suggesting a violation of the laminar fluid assumption. For the Gaussian distribution, roughly 80% of the distribution is contained within the radius where $D_{03}(r_{0})=-\frac{2}{3}$ suggesting that at least the majority of the distribution is captured by this theory. Furthermore, the the precise shape for $\rho_{x3}(r_{0})$ for a uniform and Gaussian distribution may be seen in Fig. 7; however, in the 1D case, $\rho_{01}$ truly asymptotes whereas here $\rho_{03}$ continues to decrease eventually with the uniform-like behavior of $\frac{r_{0}^{3}}{r^{3}}$. This difference is largely due to the fact that all particles asymptote to the same velocity, $c$, in the planar case but different velocities in the spherical case.

This analysis leads to the second limit of $\lim_{\zeta\to\infty}$, an unphysical limit, analogous to the 1D and cylindrically-symmetric cases where we saw the density at the edge lose dimensionality. Likewise, in the spherical-symmetric case we show in Appendix C that

[TABLE]

which is again the uniform density evolution of a symmetric distribution in one dimension less than the one being considered. While appearing unphysical, this does have some physical significance, though. This suggests that such distribution evolve toward $\rho_{x}(r_{0})$ in such a way that the factor in the denominator exactly cancels out the factor of $\frac{r_{0}}{r}$, i.e. $r^{\prime}\approx 1$ early on. However, as time progresses, the evolution shifts toward the decay of the uniform distribution in the appropriate dimension.

III Extraction field in the planar model

It is straightforward to introduce an extraction field in the planar model, and this is relevant to dynamics of electron density distributions in the pancake bunches used in ultrafast electron microscopy. The equations for this case are identical to the equations derived for the planar model in the absence of an extraction field with the single replacement

[TABLE]

or equivalently

[TABLE]

where $\rho_{r01}$ is from Eq. (33) and $\rho_{a}=\frac{eE_{a}}{mc^{2}}$, which also can be interpreted as a new length scale, $l_{a}=\frac{mc^{2}}{eE_{a}}$, associated with the extraction field. We choose the applied field to be in the positive $z$ direction. For applied fields, $E_{a}$, with $E_{a}>E_{T1}$ ($l_{a}<l_{r01}$) the applied field is sufficiently strong to overcome the space charge field throughout the bunch, hence accelerating all of the particles in the same direction. For smaller applied fields, $E_{a}<E_{T1}$ ($l_{a}>l_{r01}$), particles in the negative $z$ regions of the initial charge distribution may experience a stronger intrinsic space charge field than can be overcome by the applied field. In this case the initial distribution breaks up into two bunches moving in opposite directions. This is the virtual cathode limit defined by Valfells et al.Valfells et al. (2002). We also point out that $l_{r01}$ corresponds to the length scale of this limit. The fraction of charge in the bunch that moves in the positive and negative $z$ direction is simply $\frac{1}{2}\left(1\pm\frac{E_{a}}{E_{T1}}\right)$, respectively . The form of the two bunches is given by Eq.(40), with the substitution given in Eq. (68). Taking the relativistic limit of this expression, we obtain the asymptotic form of the two bunches

[TABLE]

and the asymptotic form for the uniform and Gaussian distributions for various applied fields are demonstrated in Fig. 8. Eq. (70) can be written in terms of the plasma period, $\omega_{01}$, and other terms, but we find that both the applied field scale, $E_{T1}$, and the associated length scale, $l_{r01}$, are more apparent in this formulation.

For the case of a uniform initial distribution $\rho_{01}=\frac{1}{2L}$ on the domain $[-L,L]$, the field of a particle at position $z_{0}$ is given by

[TABLE]

Setting the total field to zero gives the point at which the pulse breaks up into two pulses,

[TABLE]

Notice that as $E_{a}$ goes above $E_{T1}$, $z_{x}$ goes below $-L$ indicating that there is no split in the pulse consistent with prior analysis. As long as $E_{a}<E_{T1}$, the peak density of each pulse after it has gone relativistic can be calculated from Eq. (70) and is

[TABLE]

again with the positive, rightward pulse corresponding to the $+$ and the negative, leftward pulse corresponding to the $-$.

The effect of the extraction field on the time-dependent density evolution can be seen in Figs. (9) and (10), which show the evolution of initially uniform and Gaussian distributions, respectively, in the presence of various extraction fields. First, notice that the inclusion of a non-zero $E_{a}$ breaks the symmetry of the left and right pulses and that we can see that the double pulses are replaced by a single pulse as the applied field crosses the virtual cathode limit, $E_{a}=E_{T1}$. As $E_{a}$ is increased beyond $E_{T1}$, all Lagrangian particles eventually become relativistic, and the density “lifts” away from the axis. Eventually (not shown), the extraction field should be strong enough that no appreciable expansion occurs and the initial distribution is simply displaced at the speed of light; this can be shown to occur when $E_{a}>>E_{T1}$.

Also as can be seen in Figs. (9) and (10), the initial stochastic variation in the density is lost for simulations of sufficiently low extraction field but is retained for $E_{a}\geq 10E_{T1}$. This is due to the same effect discussed in the planar model without an electric field; however, the relativistic time scale needs to be adjusted, namely

[TABLE]

When $\tau_{rel}>>\tau_{exp}$, we again have the case where the expansion dynamics dominate and the inter-particle spacings essentially are equivalent to the inter-particle spacings determined by theory. However, once $\tau_{rel}<<\tau_{exp}$, the inter-particle spaces do not expand sufficiently to overcome the initial stochastics and the variance is preserved. The new wrinkle is that $\tau_{rel}$ can be reduced by simply increasing the extraction field. Therefore, we do in fact see a distribution evolve that retains the initial variance, i.e. $E_{a}=10E_{T1}$ in Fig. 9, as for that simulation $\tau_{rel}<<\tau_{exp}$.

Moreover, the influence of the extraction field is important in the highly relativistic regime not only for influencing the time scale but also influencing the asymptotic distribution. Specifically, the effect of the extraction field in the 1D model is apparently not to accelerate the front of the distribution, i.e. all simulation had the front of the bunch traveling near the speed of light, but instead to shape the eventual distribution as can be seen in Figs. (9) and (10). The asymptotic densities for the initially uniform and Gaussian distributions and for various extraction fields can be seen in Fig. 8 where we have used Eq. (45). Specifically, every point besides the point corresponding to $P_{01}+\frac{E_{a}}{E_{T1}}=0$ will eventually have $|P_{01}+\frac{E_{a}}{E_{T1}}|\frac{ct}{l_{r01}}>>1$ and thus the density corresponding to such points will eventually become a constant. However, while $|P_{01}+\frac{E_{a}}{E_{T1}}|\frac{ct}{l_{r01}}$ is not much larger than $1$, the density of the point will decrease toward the eventual constant value.

IV Application of the 1D distribution

In this section, we demonstrate one use of the spatial distributions; specifically, we calculate the width evolution of UEM-relevant planar-symmetric distributions as a function of time. Specifically, we define the rms width of the distribution as

[TABLE]

Theoretically $<a>=\int_{-\infty}^{\infty}a\rho_{1}dz=\int_{-\infty}^{\infty}a(z_{0})\rho_{01}dz_{0}$, and in simulation $<a>=\frac{1}{N}\sum_{i=1}^{N}a_{i}$ where $N$ is the number of particles in the simulation and $a_{i}$ is the value of $a$ for the $i^{th}$ particle.

We consider a Gaussian bunch with transverse radius of $100~{}\mu$m and longitudinal width of $\sigma_{r}=0.1~{}\mu$m, and we consider both $N=10^{6}$, relevant for diffraction studies, and $N=10^{8}$, relevant for imaging studies. We treat the longitudinal expansion with the planar model both using the non-relativistic distribution, Eq. (43), as well as the relativistic version, Eqs. (40) and (69). We calculate the theoretical expectation numerically for initially Gaussian distributed planar-symmetric distributions for various values of $E_{a}$ as well as the non-relativistic width prediction and compare the results to the standard deviation of $M=10^{4}$-shell simulations with the same parameters. Results of this treatment may be seen in Fig. 11.

As can be seen in Fig. 11, the theory and simulations result in the same width evolution. It is worth noting that $E_{T1}\approx 0.3MV/m$ for $N=10^{6}$ and $E_{T1}\approx 30MV/m$ for $N=10^{8}$. As can be seen in the figures, the width growth does not vary much from the unaccelerated case until an extraction field is increased beyond $\approx 10~{}E_{T1}$, that is the expansion dynamics will dominate the width determination until we are far beyond the “total” field within this 400 ps time frame. Also apparent is that, within this time frame, the dynamics of the un-accelerated bunch does not differ from the non-relativistic model; on the other hand, the higher density dynamics do differ suggesting that relativistic expansion occurs in the $N=10^{8}$ planar model. Of course, as this bunch expands its longitudinal length will quickly become larger than the transverse width suggesting higher-dimensional dynamics will become important. Specifically after the transition to higher-dimensional dynamics, the planar model overestimates the longitudinal width and underestimates the transverse width. However, if a sufficient extraction field is obtained, i.e. around $100E_{T1}$, the asymptotic longitudinal width is of sufficiently small size to result in planar dynamics for the bunch meaning that the bunch width can be modeled with the planar model in that regime.

V Discussion and Conclusions

In this work, we extended our previous density evolution analysis into the relativistic regime. Specifically, we showed that the uniform distribution in any dimension develops density shocks as the outer portion of the distribution becomes relativistic; we also found expressions for such peaks in other distributions which occur in competition with the non-uniform Coulomb mechanism that leads to peaks in such distributionsZerbe et al. (2018). We showed that the analytic results accurately predicted 1D-like $M$-shell simulation results under all symmetries and PIC results under cylindrical and spherical symmetries. The PIC simulations conducted here were completed using an EM solver with an initial ES solve used to initialize the fields. As these simulations agree with the theory that is essentially based on electrostatics, it is apparent that EM effects beyond electrostatics are not significantly affecting the density evolution for the problems examined.

We emphasize that the mechanism for the relativistic shock development is distinct from non-relativistic shock development seen in the Gaussian distributionZerbe et al. (2018). Previously, we demonstrated that shocks arise in non-planar non-uniform distribution evolutions due to the initial distribution leading to non-linear Lagrangian particle velocities that lead to inner Lagrangian particles catching up to outer Lagrangian particles. On the other hand, shock in a relativistic bunch are caused by the “shrinking” of one-dimension of the density as the Lagrangian particles approach the luminal speed limit. This can be seen by considering the energy of continuum particles within the distribution. Specifically, as the particles expand, their kinetic energies increase according to Eq. (19). In the planar and cylindrically symmetric models, this increase is linear and logarithmic in their position (and eventually time), respectively, as can be seen by Eqs. (16) and (17). This leads to all particles in a neighborhood approaching the speed of light resulting in the “freezing” of the expansion along the expanding dimension; that is, all planar symmetric distributions eventually asymptote to a time independent density while all cylindrically symmetric distributions eventually expand “uniform-like” but with one dimension less than that being considered, i.e. $\rho_{2}\to\frac{r_{0}}{r}A\rho_{02}$ where $A$ is some parameter determined from the initial conditions. On the other hand, in the spherically symmetric case, the kinetic energy is bounded by $\frac{qQ_{tot}P_{03}}{4\pi\epsilon_{0}r_{0}}$, which is finite. As this kinetic energy is dependent on $r_{0}$ through $\frac{P_{03}(r_{0})}{r_{0}}$, the asymptotic velocity of neighboring continuum particles differs. This is why the density at long times for highly relativistic portions of the distribution drops in a uniform-like manner, i.e. $\rho_{3}\to\frac{r_{0}^{3}}{r^{3}}A\rho_{03}$ with $A=\frac{1+3\zeta^{2}+2\zeta^{4}}{1+\frac{3}{2}D_{03}}$ (see Eq. (66)). While both the planar and cylindrically symmetric cases lost a power to the uniform-like evolution, i.e. planar cases evolve like $\frac{1}{r^{0}}$ and cylindrically symmetric cases evolve like $\frac{1}{r^{1}}$, the spherically symmetric case’s uniform-like evolution retains $\frac{1}{r^{3}}$. This difference arises as the particles in the spherically symmetric case have finite potential energy and thus asymptote towards a velocity that is a little less than the speed of light. Now neighboring Lagrangian particles can have very small differences in their asymptotic velocity leading to the expansion in the radial direction being slower than what is seen non-relativistically; specifically, this is what leads to $A>1$. However, there does remain a non-vanishing small velocity difference meaning that the distribution continues to expand in all dimensions. Nonetheless, these effects lead to density peaks forming towards the edges of the distribution in all cases; regardless, we find it interesting that the behavior in each dimension is qualitatively unique.

Furthermore, we demonstrated that if the distribution is given enough time to expand, the stochastic effects in the initial distribution are overwhelmed by the space charge effects. This means that under such conditions, repeated instances of similar bunches should look more or less the same. However, if the bunch is quickly accelerated into the relativistic regime, the initial stochastic fluctuations are preserved.

While these results are somewhat surprising, we do need to emphasize that these conclusion are based on analyzing the symmetric models at long times and that some of the physical assumptions inherent in the models should be violated at some point. The three assumptions for these model are (1.) the temperature is small compared to the kinetic energy delivered to the particle due to Coulomb interaction, (2.) the distributions remain laminar, and (3.) the symmetry under consideration represents the physical situation. While these assumptions all break down to some extent at some point, we emphasize that in the simulations we have conducted that the model almost exactly matches with the PIC simulations.

Note that even as the distribution becomes more and more diffuse, if a region of the distribution had much less heat than the kinetic energy delivered to it by space-charge effects, we’d expect the particles’ trajectories to not be drastically altered from the trajectory determined by the space-charge effects alone — as long as the potential energy is quickly converted to kinetic energy. Nonetheless, there is always a portion of the center of the distribution that does not meet this assumption. In the planar and cylindrically symmetric cases where the kinetic energy is unbounded, this portion of the distribution is always shrinking; on the other hand, in the spherically symmetric case, there is a portion of the distribution that will never become space-charge dominated as the kinetic energy transferred to the particles in this region will never overcome the energy associated with the initial temperature. In other words, in real world situations, the center of the distribution is emittance dominated regardless of the fact that farther out in the distribution the particles may be space-charge dominated.

The second assumption of laminar behavior is surprisingly robust. Obviously, having a higher temperature should lead to issues with this assumption, but for the cases we’ve examined within the temperature range where space-charge dominated fluid is present, this does not seem to be much of an issue — at least early on. The biggest success of the laminar fluid assumption is the planar symmetry case as the acceleration of successive sheets is monotonically increasing making laminar fluid assumption violating events impossible unless the initial velocity distribution is correctly tuned. The real issue with this assumption, though, is with the non-uniform bunches under cylindrical and spherical symmetries. As we have discussed in our previous paper, crossover events that violate the laminar fluid assumption occur when $r^{\prime}=0$. The density shock that ends up forming in the evolution of a non-uniform bunch can be thought of as occurring in region(s) of substantial initial density that have $r^{\prime}\to 0$ relatively quickly. In other words, successive cylindrical or spherical shells begin to bunch up as they expand resulting in a relatively higher density in those regions. In the non-relativistic case, these shells eventually cross-over resulting in a violation of the laminar fluid assumption (although the model still predicts the density evolution fairly well even past such events). However, considering relativity, if the initial distribution is of sufficient density, the expansion may be able to “freeze” before this point. — at least in the cylindrical case. On the other hand, relativity does not help the spherically symmetric case as complete freezing never occurs. Instead, as can be seen by analyzing Eq. (78), for the region of the distribution where $D_{03}\approx-\frac{2}{3}$ all of the terms may be relevant. As the term in front of the $\tanh^{-1}$ function is negative in this region whereas the rest of the expression is positive, there should be some value of $r_{0}$ that leads to $r^{\prime}=0$. This tells us that at some time we should expect the laminar fluid assumption to be violated by a spherically symmetric distribution. Of course for truly uniform distributions, $D_{03}=0$ everywhere and this crossover does not happen, but for any realistic distribution, all values of $D_{03}<0$ are present and crossover should occur. Specifically, roughly 20% of the Gaussian distribution has $D_{03}(r_{0})<-\frac{2}{3}$ and this crossover in this region apparently does not drastically change the evolution of the distribution for the times we’ve considered in this paper although further examination of the behavior of the model in this region is warranted. Regardless, before the time of crossovers, we are confident that the dynamics of the distribution are captured by the expressions we have presented here.

In the UEM community, the planar symmetric model is applied to a bunch that is thin along one axis with much larger widths along the other dimensions; we denote this as $L_{0}<<R_{0}$ where $L_{0}$ represents the initial width of the thin dimension and $R_{0}$ the initial widths of the other two (equivalent) dimensions. If planar symmetric dynamics are present, at some time $L\approx R$, and the planar symmetric model should no longer apply instead requiring a higher dimensional description. The time scale for the expansion of the bunch is $\tau_{exp}\approx\frac{1}{\omega_{01}}$; on the other hand, the time scale described by Eq. (73) indicates the time at which we would expect the edges of the distribution to have energy equivalent to the rest energy of the particle. As we assume $L_{0}<<R_{0}$, we’d expect that if these two timescales are of the same order or the relativistic timescale is shorter than we’d expect relativistic effects described by the models presented here to occur. This occurs when $l_{r01}\leq L_{0}$. Likewise, for the cylindrical case in fields like accelerator physics, it is generally assumed $R_{0}<<L_{0}$; which again breaks down when $L\approx R$. Nonetheless, we again expect the cylindrically symmetric dynamics described here to be apparent if $l_{r02}\leq R_{0}$.

It is straightforward to add an extraction field to the laminar theory. For $10^{8}$ electrons in a uniform bunch of radius $100~{}\mu$m, $E_{T1}\approx 30MV/m$; thus an acceleration field of $100MV$, which is the upper limit of the UEM community used at the Stanford Linear AcceleratorMusumeci et al. (2010); Weathersby et al. (2015); Murooka et al. (2011), is only about $3.5\times$ this quantity. As we saw in Figs. (9) and (10), in this range we would still expect space-charge effects that enact substantial expansion and distortion of the initial distribution. On the other hand, table top UEM devices typically have extraction fields up to $5MV/m$Srinivasan et al. (2003); Ruan et al. (2009); vanOudheusden et al. (2010); Sciaini and Miller (2011) , which is only slightly more than the total internal field of $10^{7}$ electrons in a pancake with a radius of $100~{}\mu$m, $E_{T1}\approx 3MV/m$ and is far below $E_{T1}$ for $10^{8}$ electrons. Thus $10^{8}$ electron bunches are beyond the capability of such table top devices, and $10^{7}$ electron bunches should expand immensely within the extraction field making them very difficult to work with.

Notice that in previous treatments of the evolving densityReed (2006); Zerbe et al. (2018), the extraction field was left out of the analysis. This is accurate as the density evolution in the non-relativistic limit, Eq. (43) is independent of the effective field. However, this is not true in the general case as relativistic effects make the electric field couple to the dynamics. This leads to an interesting opportunity to control the density through this coupling effect. Specifically, in 1D, the density freezing leads to the concept of asymptotic density, which is a density that no longer evolves in time. We showed that this asymptotic density can be manipulated through the inclusion of an extraction field, $E_{a}$. Specifically, the initial density is essentially the asymptotic density when $E_{a}>>\frac{\Sigma_{tot}}{\epsilon_{0}}$; however, when the extraction field is not sufficiently large, the asymptotic density can be significant different from the initial density. This suggests that if we are accelerating a bunch well into the relativistic regime, we may need to consider this asymptotic density when determining optimal criteria. Namely, in the relativistic regime, an initially uniform distribution should no longer be the distribution with the smallest emittance as relativistic considerations introduce non-linearities in the phase space that may be absent from correctly chosen initial distributions. We will develop this idea further in future work.

Appendix A Cylindrical symmetric density evolution in the highly relativistic regime

Assuming $2{\tilde{\beta}}^{2}>>1$ and analyzing all but $\mathcal{F}$ and $\mathcal{F}_{\partial}$, Eq. (55) becomes

[TABLE]

However,

[TABLE]

and

[TABLE]

Placing these approximations back into Eq. (79) results in $r^{\prime}\approx 1$.

Appendix B Long-time limit

Consider spherical symmetry. Notice that $\frac{r_{0}}{r}r^{\prime}$ is a function of $\frac{r_{0}}{r}$, so in the limit $\lim_{\frac{r}{r_{0}}\to\infty}$, these terms go to zero. As a result, $x\to 1$ in the expression for $r^{\prime}$, so

[TABLE]

where ${\tilde{\beta}}=\frac{r_{0}{\bar{\omega}}_{03}}{\sqrt{6}c}$ and where the second term is lost since inverse hyperbolic tangent goes to infinity logarithmically which is slower than $\frac{r}{r_{0}}$.

Appendix C High density limit

At high densities, the edges of the planar symmetric distribution do not significantly evolve, and therefore the distribution is essentially preserved in this region. This occurs when $\frac{2\omega_{01}L_{0}}{c}>>1$. We now extend this to the other symmetries.

C.1 Cylindrical symmetry

Assuming $\tilde{\beta}>>1$, where ${\tilde{\beta}}=\frac{r_{0}{\bar{\omega}}_{02}}{2c}$, and analyzing all but $\mathcal{F}$ and $\mathcal{F}_{\partial}$, Eq. (r’) becomes

[TABLE]

However,

[TABLE]

and

[TABLE]

Placing these approximations back into Eq. (79) results in $r^{\prime}\approx 1$. That is, $r^{\prime}$ cancels out the factor $\frac{r}{r_{0}}$ term.

C.2 Spherical symmetry

Assuming $\tilde{\beta}>>1$, where ${\tilde{\beta}}=\frac{r_{0}{\bar{\omega}}_{03}}{\sqrt{6}c}$, we see that Eq. (78) is approximately [math]; therefore we need to return to the full expression and expand in terms of $\frac{r}{r_{0}}$. We find only keeping the highest order of ${\tilde{\beta}}$

[TABLE]

So like the planar and cylindrical symmetric cases, super-highly relativistic densities result in essentially the loss of one dimension during expansion.