Diamagnetic to paramagnetic oscillations in an exploding Coulomb gas

We examined the dynamics of electron SCPs exploding from rest within a uniform magnetic field strength of $B_{ext}$ over a wide variety of initial conditions (not shown) using $N$-particle simulations, and compared the $N$-particle results to envelope equations. We term the direction along the magnetic field axis the longitudinal direction and the two directions perpendicular to this direction the transverse direction(s). $N$-particle simulations were carried out using the Velocity-Verlet algorithmVerlet (1967) with the electric field calculated using the Fast Multipole MethodGimbutas and Greengard (2015) at every step.

Before discussing specific results, we briefly introduce our spheroidal envelope model and notation. For simplicity, we assumed cylindrical symmetry (perfect Pearson correlation between $x$ and $y$), with $R^{2}=\frac{1}{2}(s_{x}^{2}+s_{y}^{2})$ representing the transverse “radial” variance, with $s_{d}^{2}$ representing the variance along $d$ and $Z=s_{z}$ the longitudinal standard deviation; the parameters $R$ and $Z$ are often termed the RMS envelope of the SCP as we will do here. Under the assumption of perfect Pearson correlation between $x$ and $y$ and no Pearson correlation with $z$, we define the radial ($\eta_{R}=\frac{s_{x,v_{x}}+s_{y,v_{y}}}{2R^{2}}$), angular ($\omega_{\rm{r}}=\frac{s_{x,v_{y}}-s_{y,v_{x}}}{2R^{2}}$), and longitudinal ($\eta_{Z}=\frac{s_{z,v_{z}}}{s_{z}^{2}}$) chirps by the line of best fit for the velocity, i.e. ${\vec{{v}}}_{bf}(r,z)={{\eta}}_{R}r\hat{R}+{{\omega}}_{\rm{r}}r\hat{\phi}+{{\eta}}_{Z}z\hat{Z}$, where $s_{x,v_{x}}$ represents the covariance of $x$ and $v_{x}$, etc. We use tilde above a parameter to indicate the dimensionless version of the parameter and $0$ in the subscript to indicate its initial value. The analysis can be made dimensionless by scaling the spatial length by the appropriate initial length and the time by the electron’s cyclotron frequency associated with the external magnetic field $\omega_{\rm{c}}=\frac{eB_{ext}}{m}$ where $e$ is the elementary charge and $m$ is the mass of the electron. Three examples of this are $\tilde{R}=\frac{R}{R_{0}}$, $\tau=\omega_{\rm{c}}t$, and ${\tilde{\eta}}_{R}=\frac{\eta_{R}}{\omega_{\rm{c}}}$. The resulting dimensionless system of equations describing the evolution of the phase space coordinates starting from rest follows:

	$\displaystyle\frac{d\tilde{R}}{d\tau}$	$\displaystyle={\tilde{\eta}}_{R}\tilde{R}$		(1a)
	$\displaystyle\frac{d\tilde{Z}}{d\tau}$	$\displaystyle={\tilde{\eta}}_{Z}\tilde{Z}$		(1b)
	$\displaystyle\frac{d{\tilde{\eta}}_{R}}{d\tau}$	$\displaystyle=\frac{1}{2}{\tilde{\omega}}_{\text{\rm{p}3,0}}^{2}\Gamma_{1}\left(\alpha\right)\frac{1}{\tilde{R}^{2}\tilde{Z}}-\frac{1}{4}+\frac{1}{4\tilde{R}^{4}}-{\tilde{\eta}}_{R}^{2}$		(1c)
	$\displaystyle{\tilde{\omega}}_{\text{r}}$	$\displaystyle=\frac{1}{2}-\frac{1}{2{\tilde{R}}^{2}}$		(1d)
	$\displaystyle\frac{d{\tilde{\eta}}_{Z}}{d\tau}$	$\displaystyle=\tilde{{\omega}}_{\text{\rm{p}3,0}}^{2}\Gamma_{2}\left(\alpha\right)\frac{1}{\tilde{R}^{2}\tilde{Z}}-{\tilde{\eta}}_{Z}^{2}.$		(1e)

where $\omega_{\text{\rm{p}3,0}}^{2}=\frac{e^{2}}{m\epsilon_{0}}\frac{N}{\frac{4}{3}\pi R_{0}^{2}Z_{0}}$ is the initial plasma frequency of the SCP (in three dimensions, hence the p$3$ in the subscript), $\alpha=\frac{Z}{R}$ is the aspect ratio, and $\Gamma_{1}(\alpha)=\frac{\alpha}{\alpha^{2}-1}\left(\alpha-\frac{\cosh^{-1}\alpha}{\sqrt{\alpha^{2}-1}}\right)$ and $\Gamma_{2}(\alpha)=\frac{\alpha}{\alpha^{2}-1}\left(\frac{\cosh^{-1}\alpha}{\sqrt{\alpha^{2}-1}}-\frac{1}{\alpha}\right)$ are functions that determine the effect geometry has on the force. We provide a fluid derivation of this system of equations in Section XXX of the supplemental. Similar envelope models are well established in the literatureLund et al. (2014); in fact, the philosophy for presenting this model here is to have minimal complexity for understanding the physics while maintaining the agreement between the model and the dynamics of the situation.

We note that Eq. (1d) is not typically reported as part of the envelope equations; however, equivalent forms of this expression can be found in the literature describing the equilibrium distribution of SCPs (for instance, Eq. (5.36) ofReiser (1994)). This equation is always derived from a consideration of the conservation of angular momentum, as we have done in Section XXX of the supplement. However, our use of this equation is somewhat different from the standard use. In most cases where this equation may be found, the researcher is looking at a plasma at equilibrium. Therefore, in addition to Eq. (1d), equilibrium analyses also assume $R$ is constant in time, i.e. $\eta_{R}=0$ and $\eta_{Z}=0$ or more generally $(\eta_{X},\eta_{Y},\eta_{Z})=(0,0,0)$. Substantial effort has been made to describe such distributions of interacting SCPs using these constraints [refs]; namely, we know “equilibrium beams” satisfies Eq. (1d) and $R=$ a constant by adopting a non-uniform profile when the generalized angular momentum is non-zeroBogema Jr and Davidson (1970). However, this well-established observation need not necessarily be true for far-from-equilibrium SCPs and beams. Namely, if the radius is not constant, the distribution may be undergoing oscillations while remaining nearly uniform. We examine this possibly controversial statement further in the next section.

For our first representative simulation, we initialized $10,000$ electrons with spatial coordinates sampled from a spherical-uniform distribution with radius $2~{}\mu$m and evolved this distribution from rest under the influence of the particles’ pair-wise repulsion and an external magnetic field strength of $4$ T. Fig. 1 shows the evolution of the longitudinal magnetic moment from both the $N$-particle simulation and the envelope equations starting from rest. As can be seen, the $N$-particle results and envelope prediction are in near perfect agreement. As the longitudinal magnetic moment reverses direction twice during every oscillation, the magnetic moment is aligned with the external magnetic field, i.e. is diamagnetic, and opposite to the external magnetic field, i.e. is paramagnetic, for periodic intervals during this oscillation. We note this oscillation is consistent with the Bohr-Van Leeuwen theoremVan Leeuwen (1921) as our system is far-from-equilibrium; likewise, it is not inconsistent with the standard attribution of magnetically confined plasmas as being diamagneticO’Neil (1995) as such a description again assumes equilibrium-like conditions.

As the longitudinal magnetic moment, $m_{z}$, in our notation is given by

$\displaystyle m_{z}$

$\displaystyle=-\frac{e}{2}\omega_{\text{r}}R^{2}\hat{z},$

(2)

we see that the sign of the magnetic moment of the SCP is determined by the direction of rotation captured by the angular chirp. The oscillation of the angular chirp for our first characteristic simulation is displayed in Fig. 2a; again, these results are in excellent agreement. As the radial motion is captured in the model through the evolution of the parameter $R$, which is the sole time-dependent parameter in Eq. (1d), we see that $\omega_{\text{r}}$ couples to the radial dynamics so that $\omega_{\text{r}}$ and $R$ oscillate at the same frequency. The envelope predicted and $N$-particle simulated oscillation of $R$ and the unbound expansion of $Z$ are seen in Figs. 2b and 2c, respectively. Notice ${\tilde{\omega}}_{\text{r}}$ is negative for any value of $R$ that is less than $R_{0}$ for simulations starting from rest; that is, the transition between the diamagnetic and paramagnetic states occurs when the transverse radius crosses its original size as can be seen by comparing Figs. 1 or 2a with Fig. 2b. We note that the evolution of $Z$ seen in Fig. 2c is typical of unconfined Coulomb explosion; namely, the rapid acceleration of the growth in the statistic followed by a period of almost linear growth leading to a consistent decrease in the density and hence a weakening of space charge effects.

As this effect is in $3D$ and the expanding longitudinal dimension results in a highly dynamic density, it is not apparent that this effect is the same for every magnetic field. As the only frequencies relative to this situation are the cyclotron and plasma frequencies, we consider this effect as a function of $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$. In Fig. 3a, we show how the oscillatory dynamics respond to changes in $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ when all other parameters are held constant. The amplitude of the oscillations decreases as the ratio between the frequencies is increased; interestingly, though, the size of the oscillation’s amplitude within the paramagnetic regime increases to a point at which further increasing of the magnetic field likewise begins to decrease the size of the oscillation in the paramagnetic region. In other words, at very low values of $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$, the self-magnetic field oscillation does not enter very far into the paramagnetic region and such oscillation is effectively turned on by increasing the magnetic field.

To obtain better intuition into this effect, we quantified the relative size of the paramagnetic and diamagnetic oscillations by looking at the ratio of the minimum to the maximum self-magnetic moments, i.e. $\frac{m_{z,\rm{para,max}}}{|m_{z,\rm{dia}}|_{\rm{max}}}$. As can be seen in Fig. 3a, the fist diamagnetic peak is somewhat larger than subsequent peaks, so this ratio in our envelope treatment effectively uses the size of this first peak in the diamagnetic region to gauge the size of the subsequent peaks in the paramagnetic region. Moreover, the subsequent peaks within the paramagnetic region slowly increase toward an asymptote as can be seen in Fig. 3b — this is due to the SCP becoming less and less dense as it expands longitudinally resulting in the space-charge effect becoming less and less important in determining where the oscillating SCP will reverse its inward motion. On the other hand, as we increase the magnetic field, the number of oscillations occuring before the space-charge effect becomes negligible increases. As we are comparing simulation results at the same multiple of the the cyclotron period, this means that it takes the space-charge effect will emerge later in simulation with relative larger cyclotron frequencies. This is the origin of the peak seen in Fig. 3a.

As can be seen, increasing the ratio $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ does appear to turn this oscillation on — specifically, the turning on occurs when $\omega_{\rm{c}}^{2}/\omega_{\rm{p}3,0}^{2}$ is in the range $[0.1,3]$. Note that the “Brillouin point” occurs when the initial angular motion is rotating at $\frac{1}{2}\omega_{\rm{c}}$ and $\omega_{\rm{c}}^{2}/\omega_{\rm{p}3,0}^{2}=2$, but here there is no initial angular motion and the oscillation is driven entirely by the initial expansion of the electron SCP; therefore, the occurence of $2$ in this range is unrelated to the Brillouin point.

Comparing the results of the RMS envelope equations to the $N$-particle simulations, the general trends are captured by the model in both Figs. 3a and 3b and are almost exact at earlier times. Notice that the deviations in the ratio monotonically increase during the simulation and are larger for stronger external magnetic fields, and a close inspection of the magnetic moment oscillation for higher magenetic fields can be seen in Fig. 3a. Coincidentally, these are the cases that have larger emittance growth, so we infer these deviations are due to accumulated emittance growth that are not included in Eqs. (1). We will model the effect of RMS emittance in this situation in a later section, but such treatment will still consider the emittance as being conserved — just non-zero — but still not dynamic as it is in the simulations.

As mentioned in the previous section, substantial work has been done to describe the steady-state properties of confined plasmas[refs]. Initial work by Brillouin described a uniform steady-state where the angular chirp is given by the Larmor frequency, $\omega_{\rm{r}}=\omega_{\rm{L}}=\frac{1}{2}\omega_{\rm{c}}$Brillouin (1945), and subsequent work has extended this understanding to non-uniform steady-state distributions rotating with angular chirps in the range $[0,\omega_{c}]$ known as rigid-rotor equilibrium.

To our knowledge, there is no clear quantitative definition of rigid-rotor equilibrium; so we develop one here. If we were to calculate the angular kinetic energy of a rigid-rotor, we would expect that the amount of angular kinetic energy predicted by the linear chirp, i.e. $\frac{N}{10}m\omega_{\rm{r}}^{2}R^{2}$, should be much larger than the amount of angular kinetic energy associated with random motions. Consider the angular velocity of the $i^{th}$ particle, $v_{\phi,i}$. We know that the residual for the $i^{th}$ particle between the chirp-predicted and measured angular velocity can be written as

$\displaystyle\delta_{v_{\phi},i}$

$\displaystyle=v_{\phi,i}-\omega_{\rm{r}}r_{i}.$

(3)

where $r_{i}$ and $v_{\phi,i}$ are the radial and angular velocity coordinates for particle $i$. Therefore the variance of the residuals can be written as

$\displaystyle s_{\delta_{v_{\phi}}}^{2}$

$\displaystyle=s_{v_{\phi}^{2}}+\omega_{\rm{r}}^{2}s_{r}^{2}-2\omega_{\rm{r}}s_{v_{\phi}r}.$

(4)

Note that $s_{r}\neq R$ as $r>0$ for all $r$ whereas $\bar{x}$ and $\bar{y}$ are approximately $0$. As $\frac{N}{2}ms_{\delta_{v_{\phi}}}^{2}$ represents the portion of the kinetic energy that is associated with motions that are not linearly correlated and is often attributed to being random, it follows that the ratio

$\displaystyle\nu$

$\displaystyle=\frac{1}{5}\frac{\omega_{\rm{r}}^{2}R^{2}}{s_{v_{\phi}^{2}}+\omega_{\rm{r}}^{2}s_{r}^{2}-2\omega_{\rm{r}}s_{v_{\phi}r}}$

(5)

can quantify the rigid-rotor approximation; we call $\nu$ the rigid-rotor statistic. Indeed if the rigid-rotor statistic is larger than $1$, the majority of the angular kinetic energy is associated with the motion described by the angular chirp.

For a second characteristic simulation, we show the rigid-rotor statistic as a function of time in Fig. 4a. This second simulation again had $10,000$ electrons but was initially sampled from the spherical uniform distribution with initial radius $R_{0}=5~{}\mu$m from rest with a magnetic field strength of $1$ T. As can be seen in this figure, the rigid-rotor statistic is much larger than $1$ except when it crosses $R_{0}$ at which time the SCP is in the process of reversing its direction of rotation; at this time, the angular chirp is zero and therefore the kinetic energy associated with the coherent portion of the rotation captured by this statistic is zero as well. Such dynamics should not be described as rigid-rotor-like. However, as the radius continues to decrease, the angular chirp becomes negative and the rigid-rotor statistic once again become greater than $1$. This means the distribution begins to rigidly rotate in the opposite direction after briefly passing through a state dominated by random angular motion. Due partially to the fact that the expansion beyond $R_{0}$ is relatively much larger than the contraction within $R_{0}$, the rigid-rotor statistic for this latter region does not get as large. Moreover, the rate at which the rigid rotor statistic evolves within the larger than $R_{0}$ and smaller than $R_{0}$ regions differ as can be seen in Figs. 4a and 4b.

It is unclear from $N$-particle simulations alone to what extent each of non-zero emittance, space-charge, and emittance evolution effects contribute to evolution of the rigid-rotor statitics seen in Fig. 4a especially once the SCP has expanded significantly. To obtain insight into such dynamics, we re-simulated the last two transverse oscillations without inter-particle forces starting with the same configuration produced from the interacting simulation of Coulomb explosion. We show the comparison between the results of these non-interacting and interacting simulations in Fig. 4c. We note that the evolution of the rigid-rotor statitic is very similar in both cases suggesting that non-interacting-like behavior dominates towards the end of the transverse oscillations; such behavior would also suggest that emittance growth is responsible for the small deviations. On the other hand, we would expect that emittance effects should dampen the rigid-rotor statistic, but what is seen is the opposite — that is, the rigid-rotor statistic for the interacting simulation is larger than the non-interacting simulation.

As the rigid-rotor statistic is correlated with the transverse radius of the SCP, one explanation for this effect could be if the interacting simulation were to obtain a larger radius than the non-interacting simulation thus affecting the numerator of the rigid-rotor statistic. In Fig. 5a, the evolution of the non-interacting and interacting radii are shown. As can be seen, the simulated SCP with no interaction predicts a less extreme transverse oscillation. Thus this explanation is consistent with the observed difference between the maxima of the rigid-rotor statistics in Fig. 5c; however, a small deviation is also observed at the minimum radius, and the rigid-rotor statistic does not display the complementary deviation in the rigid-rotor statistic during this time period. So something more is going on near the transverse minima that compensates for the radial effect seen at the maximum.

We explored the densities predicted near the final simulated transverse minimum to understand if a space-charge shift in the bunch profile could be compensating for the radial-induced expected deviation near the minima. As the second-order moments are sensitive to profile changes, such a profile change should result in a change in the rigid-rotor statistic that we would then attribute to the use of second order moments. We show the transverse radial density at the $14^{th}$ transverse radial minimum for both simulations in Fig. 5b. As can be seen, the simulated radial densities do in fact deviate substantially from the theoretically uniform distribution, but they do not deviate from one another. Specifically, a clear edge effect that is well discussed in the equilibrium literature [refs] arises when the bunch reaches its minimum transverse radius, but again, this occurs in both the interacting and non-interacting simulations. Therefore, this density effect cannot be causing the difference in the rigid-rotor statistic. Further, as the only portion of the non-interacting dynamics that is driving distribution change is the emittance, we conclude that the edge effect is solely an emittance effect in both cases.

A second possible mode of compensation would be if our statistical estimate of the kinetic energy in the angular mode from the angular chirp could differ between the two simulations if non-linearities were to arise. We show the angular phase-space at the $14^{th}$ transverse radial minimum for both simulations in Fig. 5c. As in the radial density, a clear edge effect, which could be described as a shear, can be seen; however, again the same effect is seen in both the interacting and non-interacting simulations. Therefore, we conclude that the compensation in the rigid-rotor statistic is not driven by this shear effect and further that the shear is likewise a result of emittance.

A third possible compensatory mechanism is if the denominator in Eq. (5) were to somehow differ between the two simulations. To examine this denominator, we point out that normalized RMS emittance, standardly defined [refs] as

$\displaystyle\varepsilon_{x}$

$\displaystyle=\sqrt{s_{x}^{2}s_{v_{x}}^{2}-s_{x,v_{x}}^{2}}$

(6)

and defined transversely in the presence of cylindrical-symmetry about a magnetic field byLund et al. (2014)

$\displaystyle\varepsilon_{\perp}$

$\displaystyle=\sqrt{\frac{1}{4}(s_{v_{x}}^{2}+s_{v_{y}})^{2}-R^{2}(\eta_{R}^{2}+\omega_{r}^{2})},$

(7)

is effectively the velocity spread times the spatial spread when the bunch is projected to a single dimension [refs]. Therefore, we multiplied the velocity spread, or the square root of the residual variance from Eq. (4), by $R$ to define an “angular emittance”,

$\displaystyle\varepsilon_{\phi}$

$\displaystyle=Rs_{\delta_{v_{\phi}}}.$

(8)

As can be seen in Fig. 6a, this measure is very roughly conserved during the non-interacting simulation. Interestingly, it is also roughly the same value as the cylindrically-symmetric emittance seen in Fig. 6b. The small emittance growth seen in both the angular and transverse emittances appears to contribute to the fact that the rigid-rotor statistic is smaller at the minimum of the transverse width than should be expected from the evolution of the radius alone.

We note that this roughly conserved angular emittance has not been described in the literature but can be argued to be the complement of the round-beam thermal emittance defined by Lund et alLund et al. (2014) presented in Eq. (7). Also interestingly, obtaining the Cartestian covariances and re-sampling the distribution to make an equivalent initial condition [not shown] neither captures the correct angular emittance nor the evolution of the rigid-rotor statistic. This suggests that the standard statistics are not capturing all of the spatial correlations present in the SCP.

We now consider cases where the SCP starts with some temperature, which we quantify by the dimensionless emittances defined by

	$\displaystyle\varepsilon_{\perp}^{2}$	$\displaystyle=\frac{R^{2}}{R_{0}^{4}\omega_{\rm{c}}^{2}}\left(\frac{s_{v_{x}}^{2}+s_{v_{y}}^{2}}{2}-\left(\eta_{R}^{2}+\omega_{\rm{r}}^{2}\right)R^{2}\right)$		(9a)
	$\displaystyle\varepsilon_{z}^{2}$	$\displaystyle=\frac{1}{R_{0}^{4}\omega_{\rm{c}}^{2}}\left(s_{z}^{2}s_{v_{z}}^{2}-s_{z,v_{z}}^{2}\right)$		(9b)

where $\varepsilon_{\perp}$ is a dimensionless version of Eq. (7). Following the formulation of the envelope equations by Lund et alLund et al. (2014), we add standard definitions of emittance to our equations with the additional assumption that we are starting with non-zero angular motion, $\omega_{\rm{r},0}$, resulting in changes to Eqs. (1c), (1d), (1e):

	$\displaystyle\frac{d{\tilde{\eta}}_{R}}{d\tau}$	$\displaystyle=\frac{1}{2}{\tilde{\omega}}_{\text{\rm{p}3,0}}^{2}\Gamma_{1}\left(\alpha\right)\frac{1}{\tilde{R}^{2}\tilde{Z}}-\frac{1}{4}$
		$\displaystyle\quad+\left(\left({\tilde{\omega}}_{\rm{r},0}-\frac{1}{2}\right)^{2}+\varepsilon_{\perp}^{2}{\tilde{\omega}}_{\text{\rm{p}3,0}}^{2}\right)\frac{1}{\tilde{R}^{4}}-{\tilde{\eta}}_{R}^{2}$		(10a)
	$\displaystyle{\tilde{\omega}}_{\text{r}}$	$\displaystyle=\frac{1}{2}-\left({\tilde{\omega}}_{\rm{r},0}-\frac{1}{2}\right)\frac{1}{{\tilde{R}}^{2}}$		(10b)
	$\displaystyle\frac{d{\tilde{\eta}}_{Z}}{d\tau}$	$\displaystyle=\tilde{{\omega}}_{\text{\rm{p}3,0}}^{2}\Gamma_{2}\left(\alpha\right)\frac{1}{\tilde{R}^{2}\tilde{Z}}+\frac{\varepsilon_{z}^{2}}{\tilde{Z}^{4}}-{\tilde{\eta}}_{Z}^{2}$		(10c)

Our addition of the transverse and longitudinal emittance terms can be argued analogous to Lund et al. or derived completely from a statistic standpoint (in preparation); angular emittance is not added to Eq. (10b) as the equation is a statement of conservation of angular momentum. The addition of the dimensionless normalized RMS emittance terms captures the effect of motion that deviates from the linear chirps, sometimes described as thermal effects [refs], that result in an outward-like force on the RMS envelope parameters. This is analogous to temperature pressure in fluid models [refs], but is formulated from a discrete statistical vantage point where no theoretical approximations are needed.

To examine the effect of the emittance on the diamagnetic/paramagnetic oscillation, we conducted $N$-particle simulations of initially spherically symmetric Coulomb explosion, initially with $\varepsilon_{z,0}=\varepsilon_{\perp,0}$, and integrated the envelope equations. For the initial conditions we held $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ fixed and varied $\varepsilon_{\perp,0}$. The evolution of some such characteristic simulations can be seen in Fig. 7a. What is apparent is that increasing the initial emittance leads to larger oscillations, shifts the frequency of oscillations higher and decreases the size and duration of the magnetic moment within the paramagnetic region — in a way, this is nearly the opposite of the effect of increasing the confining magnetic field. A plot of the ratio of the maximum self-magnetic moments in each region with such initial conditions viewed as a function of $\varepsilon_{\perp,0}$ can be seen in Fig. 7b. The envelope data is shifted rightward relative to the $N$-particle simulated data as the $N$-particle simulations undergo disorder induced heating that increases the emittance. As can clearly be seen, increasing the dimensionless emittance beyond about 0.1 decreases the paramagnetic/diamagnetic statistic, and the effect is again turned off once the dimensionless emittance is increased past about 2.

Our envelope model predicts the behavior of the $N$-particle simulation well especially before the emittance has grown significantly. We have also produced a figure detailing the interaction between $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ and the emittance using just the envelope predictions within the first $6.5$ plasma periods. The competitive nature of these two effects can be seen in Fig. 8; namely, an increase in the emittance requires larger values of $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ before a substantial paramagentic state is observable. Further, increasing the emittance delays at what ratio of $\frac{1}{{\tilde{\omega}}_{\rm{p}3,0}^{2}}$ this effect arises as can be seen by the fact that the curves with significant emittance in Fig. 8 appear to be scaled versions of the successive evolution of this statistic seen in Fig. 3b. As this delay enables additional effects, like cyclotron radiation and other electromagentic effects, outside of the scope of this paper to occur and further drive the SCP toward equilibrium, sufficiently large emittance will effectively turn this effect off.

So far, our analysis has been principally of numerical results. The spheroidal envelope model does not have an exact analytical solution; however, simpler models do allow more analysis. If we were to match such a model to an experiment where the space charge effect remained important even once the bunch had become significantly oblate, it would require us to fit the initial parameters of the model to the parameters observed analogous to our non-interacting treatment. We analyze this continuous-beam model.

The continuous model is very similar to the spheroidal model. The principal differences are: (1) there is no longitudinal dimension to model and (2) the geometric factor is absorbed into the plasma frequency giving a cylindrically symmetric plasma frequency. We denote this new plasma frequency by $\omega_{\text{p2,0}}=\sqrt{\frac{e^{2}}{m\epsilon_{0}}\frac{n_{\text{z}}}{\pi R_{0}^{2}}}$ where $n_{z}$ is the number density along the length of the beam. The subsequent radial oscillation is determined by

$\displaystyle\frac{d{\tilde{\eta}}_{R}}{d\tau}$

$\displaystyle=\frac{1}{2}{\tilde{\omega}}_{\text{p2,0}}^{2}\frac{1}{\tilde{R}^{2}}-\frac{1}{4}+\left(\left({\tilde{\omega}}_{\text{r,0}}-\frac{1}{2}\right)^{2}+\varepsilon_{\perp,0}^{2}{\tilde{\omega}}_{\text{p2,0}}^{4}\right)\frac{1}{\tilde{R}^{4}}.$

(11)

Note the evolution of the angular frequency is only changed by the third dimension only through its dependence on $R$ as we assume angular momentum is completely along the axis of symmetry. We note that many equivalent models to this can be found in the literatureDubin (1993); Reiser (1994). While this model depends on three parameters, ${\tilde{\omega}}_{\text{p2,0}}^{2}$, ${\tilde{\omega}}_{\text{r,0}}$, and the initial radial velocity, the freedom to assign the initial phase of the oscillation allows us to set one of the initial velocities, by convention the initial radial velocity, to zero. Thus the properties of the oscillations predicted by this model can be understood within the $\left({\tilde{\omega}}_{\text{r,0}},{\tilde{\omega}}_{\text{p2,0}}^{2}\right)$ parameter space when $\varepsilon_{\perp,0}$ is held fixed, and we will conduct our analysis within this parameter space for various values of $\varepsilon_{\perp,0}$.

We note that this parameter space is standard in descriptions of rigid-rotor equilibria; namely, this parameter space is typically used to visualize the parabola of zeros of the effective force that describe the well known rigid-rotor equilibriaGould (1995). In the accelerator literature where continuous beam behavior near rigid-rotor states is discussed, small radial oscillations about a rigid-rotor equilibrium are described as RMS-mismatched induced transverse ripples along the longitudinal coordinateReiser (1994); however, such descriptions focus solely on the effect of $R$ and not its coupling to $\omega_{\text{r}}$, which we emphasize as being important here.

A diagram denoting all of the regions of parameter space where the initial conditions of the continuous-beam model will experience the phase transition for various values of $\varepsilon_{\perp,0}$ is displayed in Fig. 9. One interesting aspect of this diagram is the fact that the two shaded regions, one inside the rigid-rotor parabola and a second describing the half-plane with $\omega_{\text{r,0}}<0$, describe the same oscillations. That is, our convention of no initial radial velocity is ambiguous as there exists two times per period with zero radial velocity: (1) the time where the radius is maximum (corresponding to the entire region of parameter space inside the rigid-rotor parabola) and (2) the time at which the radius is minimum (that region’s complement). Deducing the fact that all initial conditions within the half plane where $\omega_{\text{r,0}}<0$ produce the paramagnetic oscillations is simple; we know that the point in the oscillation that corresponds to the maximum radius is inside the parabola and has $\omega_{\text{r}}>0$, so all such initial states in this region must change their direction of rotation at some time. The initial states within the rigid-rotor equilibrium parabola that experience the phase transition requires more analysis, and we present the mathematical description of this region’s boundary in Section VI of the supplement.

We have described the dynamics of Coulomb explosion of an electron SCP within a magnetic field where we have found that the rotation of the SCP can spontaneously reverse direction due to a combination of the space-charge effect and the conservation of angular momentum. This reversal of the rotational direction of the SCP results in the SCP periodically switching between diamagnetic and paramagnetic states. This observation is surprising given the established diamagnetic nature of SCPs in ubiquitous Brillouin flow models. We showed that this behavior is turned on by increasing the magnetic field, and turned off by increasing the emittance or thermal kinetic energy in the SCP. Moreover, we showed that this behavior arises from coupling between the radial and angular motion introduced by the magnetic field, and presented equations that describe this coupling. We modeled the radial and hence angular motion of such SCPs with mean-field particle-particle interactions through envelope equations and showed that the envelope equations captures the angular dynamics to the same extent that the radial dynamics are captured.

It is clear that SCPs particularly in traps rapidly evolve toward an equilibrium distribution. During this time, the plasma typically undergoes disorder induced heatingMurillo (2006) as well as cooling through cyclotron radiationO’Neil (1980). However, a full kinetic picture of how such a plasma evolves from a possibly far-from-equilibrium initial state with substantial inter-particle correlation to the eventual equilibrium distribution where particle inter-particle correlations are treated randomly is not apparent within the literature. Our Coulomb explosion model provides insight into some possible early dynamics that could further evolve the inter-particle correlations if relativistic effects were to be properly considered in future work.

We note the parameters we use have statistical definitions. We have been calling the description of the evolution of such parameters statistical kinematics although they are sometimes known as moment equations in the literatureSacherer (1971). These statistical definitions have the benefit of being applicable to arbitrary distributions; in this work, this is true for the radial measure, $R$, angular measure, $\omega_{\rm{r}}$, and magnetic moment, $m_{z}$. However, to capture the time dependence of one of these parameters, specifically the radial measure, we needed to apply mean field theory using the uniform distribution to avoid space-charge non-linearities for which general techniques have yet to been developed. So despite the fact that Eq. (1d) relates the dependence of $\omega_{\rm{r}}$ on $R$ in any generic distribution, we can only capture the time evolution through the use of the uniform distribution for which an analytic expression for the space-charge force is available.

The insight into the radial/angular coupling adjusts the typical depiction of the SCP breathing mode as a periodic inward and outward motionDubin and Schiffer (1996). Instead, we’ve shown that the angular motion oscillates coherently with the radial motion. That is, instead of moving in and out, the SCP also rotates with rate of the rotation specified by the size of the radius. The microscopic origin of this is of course the cyclotron motion of the particles, and we have provided a simplified schematic of this motion for the non-interacting distribution in the supplemental documentation.

The fundamental theory can be tested experimentally by careful accounting for the mechanical angular momentum of a beam. We note that generalized Courant-Snyder theoryQin et al. (2013) provides tools to track the rotation angle of the envelope matrix in the presence of transverse-coupled dynamics, and this rotation can alternatively be conceptualized as the angular momentum of the beam as we’ve done here. Therefore, that theory and the tools we present here may be used to do such accounting. More specifically, a long solenoid can be used to introduce mechanical angular motion to a prolate SCP in an accelerator setting; the amount of angular momentum can be determined using Eq. (1d) and (11) and knowledge of the angular momentum kick across the solenoid’s fringe. As the SCP evolves in a region without a magnetic field, these equations may again be employed with $\omega_{\text{L}}=0$ and the $\omega_{r,0}$ introduced from the solenoid. The SCP could then be passed through a second long solenoid. By tuning the length of the first solenoid, the magnetic field strengths, and the separation distance between the solenoids for an SCP of a known length and number of electrons, it should be possible to induce a near equilibrium SCP once it enters the second solenoid.

Work by R.M.V. and O.Z. was in part supported by a Packard Fellowship in Science and Engineering from the David and Lucile Packard Foundation. Work by P.D. and B.Z. was supported by NSF Grant numbers RC1803719 and RC108666.

O.Z. first noticed the phase transition. O.Z. and B.Z. developed theoretical model, wrote code, and analyzed the simulation and model data. O.Z., R.M.V., P.M.D., and B.Z. interpreted the data and wrote the manuscript.

Data available on request from the authors. Simulation and envelope code available on request from authors.

The authors declare no competing interests.