Nonlinear equilibria and emittance growth in plasma wakefield accelerators with ion motion
C. Hansel W. An W. Mori J. B. Rosenzweig
physics.acc-ph





claire.hansel@ucla.edu




 Department of Physics and Astronomy, University of California, Los Angeles, California 90095 






The plasma wakefield accelerator may accelerate particles to high energy in a future linear collider with unprecedented acceleration gradients, exceeding the GeV/m range. Beams for this application would have extremely high brightness and, subject to the intense plasma ion-derived focusing, they would achieve densities high enough to induce the plasma ions to collapse into the beam volume. This non-uniform ion density gives rise to strong nonlinear focusing which may lead to deleterious beam emittance growth. The effects of ion collapse and their mitigation has been investigated recently through particle-in-cell simulations, which show that by dynamically matching the beam to the focusing of the collapsed ion distribution, one may avoid serious emittance growth. We extend this work by exploring the near-equilibrium state of the beam-ion system reached after the ions have collapsed, a condition yielding the emittance growth mitigation observed. We show through PIC simulations  <cit.> and analytical theory that in this case a dual electron beam-ion Bennett-type equilibrium distribution is approached. Here, the beam and ion distributions share nearly the same shape, which generates nonlinear transverse electromagnetic fields. We exploit a Bennett-type model to study beam phase space dynamics and emittance growth over time scales much longer than permitted by PIC simulations through use of a 2D symplectic tracking code with Monte Carlo scattering based on Moliere's theory of small angle multiple scattering. We find that while phase space diffusion due to parametric excitations of the beam size due to plasma non-uniformity is negligible, scattering from collapsed ions gives rise to manageable emittance growth in the case of a linear collider. The implications of these results on experiments planned at FACET-II are examined.



Nonlinear equilibria and emittance growth in plasma wakefield accelerators with ion motion
    J. B. Rosenzweig
    December 30, 2023
==========================================================================================





§ INTRODUCTION


Due to radio-frequency cavity breakdown <cit.> and related phenomena, current conventional particle linear accelerators are limited to a maximum practical acceleration gradient not notably in excess of 100 MeV/m. To diminish the size of high energy accelerators in applications exemplified by the linear collider (LC), a number of advanced acceleration techniques have been proposed and are in various stages of development.  One such technique is plasma wakefield acceleration, which employs waves in a plasma - a medium chosen to naturally evade breakdown issues - excited by an intense drive beam, to accelerate a trailing beam. Plasma wakefield accelerators (PWFAs) have already demonstrated acceleration gradients in excess of 50 GeV/m <cit.>, orders of magnitude above what is achievable with conventional accelerators. Further scenarios have been identified that extend this gradient reach to beyond a TeV/m <cit.>. To explore the suitability of the PWFA for LC applications, proposals <cit.> have been put forward that analyze the use of an “afterburner" at the end of a conventional LC injector, which have the goal of doubling the beam energy available for high energy physics exploration. There is considerable worldwide research interest in development of the PWFA for this high impact application, with major new facilities now being commissioned to explore the physics issues supporting PWFA development such as FACET-II <cit.> and FLASHForward <cit.>.

Current experimental and theoretical research on PWFAs has focused on the nonlinear “blowout" regime due to its favorable properties for acceleration and focusing. In the blowout regime the electron beam density is much greater than the ambient plasma electron density, and the collective fields of the beam eject (or “blow out") the plasma electrons from the region near the beam path, and form a bubble of negligible electron density <cit.>. An  electromagnetic wave is trapped inside of this bubble that provides acceleration in uniform phase fronts, as in standard relativistic electron accelerators. Further, in this scenario the plasma ions left behind, if undisturbed, provide a uniform charged column that yields strong, linear (emittance preserving) focusing. In this way, one may achieve high quality, low energy spread acceleration without emittance growth due to geometric aberrations. However, in the proposed PWFA afterburner case, the low emittance, high current LC beams acting under the influence of strong ion focusing will move significantly during beam passage. The plasma ions, previously treated as stationary with constant density n_0, can then rearrange dramatically, collapsing into a tight column with much enhanced, non-uniform density <cit.>. This in turn causes the beam to focus further under the influence of large, nonlinear forces, potentially inducing notable emittance growth. 

The challenge of controlling such emittance growth has been known for over a decade, but recently a robust solution has been found through an analysis based on particle-in-cell (PIC) simulations <cit.>. This work demonstrates that that by matching the initial beam phase space to the eventual focusing arising from the collapsed ion distribution, a quasi-equilibrium scenario develops, in which emittance growth is mitigated to a nearly ignorable level. This quasi-equilibrium has been investigated further by Benedetti, et al., <cit.> in a transient regime, where its nature is not clear, as the result of the analysis is a complicated analytical description of the beam-plasma state. Here, we extend these previous explorations by examining the fundamental characteristics of this matched, nonlinear equilibrium, in the limit where it is established (collisionlessly, through nonlinear phase mixing) in a near steady state over the majority of the beam. The analytical form of this equilibrium is available through a Maxwell-Vlasov description, and is found to be a Bennett-type distribution where the ion and electron distributions follow nearly the same radial profile. The establishment of this equilibrium is described theoretically and verified through PIC simulations. 

Due to the computationally intensive nature of these simulations, the evolution of the beam's phase space can only be modeled on a short time scale of a few betatron oscillations. As a consequence, PIC simulations cannot model effects that take place over longer timescales such as diffusion induced by the nonlinear focusing or scattering due to the high plasma ion density. We have developed a tracking code that exploits the analytically known focusing forces in equilibrium to investigate the electron beam behavior over timescales much longer than those that can be feasibly simulated with PIC codes. In these tracking simulations, we are concerned with two effects: diffusion in phase space due to parametric excitations and Coulomb scattering due to the very dense (over 10 times atmospheric density) plasma ions. We use this code to estimate the observability of phase space dilution and concomitant emittance growth in a relevant LC case, and for the parameters of the E314 Ion Motion experiment currently under development at the SLAC wakefield research facility FACET-II. We discuss aspects of this experiment in light of the physics understanding developed from simulation and theory. 

With these issues in mind, we organize this paper as follows. In Sec. (<ref>) we derive the Maxwell-Vlasov equilibrium equations that yield the Bennett-type distributions, and give expressions for the the Bennett-type profile parameters expected in terms of beam and plasma initial conditions. In Sec. (<ref>) we display and discuss particle-in-cell simulations that illustrate the collapse to equilibrium and validate the description of the equilibrium as having Bennett-type form. In Sec. (<ref>) we describe the tracking code that permits evaluation of long-time phase space dynamics in the beam. In Sec. (<ref>) we discuss the possibility, based on the use of this tracking code, of observing diffusion due to chaos and parametric resonances driven by plasma density fluctuations. Finally in Sec. (<ref>) we discuss conclusions, implications for present experiments, and future concerns related to a PWFA-based LC. Additionally, Appendix <ref> contains a derivation of the multiple scattering angular distribution as well as a description of how the tracking code was validated.



§ BENNETT-TYPE EQUILIBRIA
 

In this section, we develop the analysis of the Maxwell-Vlasov equilibrium that arises due to this two species – electron beam and collapsed ion population – interaction. It is in many ways an extension and generalization of the magnetically self-focused beam scenario originally analyzed by Bennett in 1934.



 §.§ Equilibrium Conditions


We begin with a brief description of the Maxwell-Vlasov analysis leading to the emergence of a Bennett-type equilibrium. Let n_e and n_i be the number densities of beam electrons and plasma ions respectively which are assumed to be cylindrically symmetric. The beam is taken to travel with approximately constant velocity β c ẑ, i.e. predomininantly in the z-direction. The plasma ions are assumed to be in a thermal equilibrium without a net directed flow in z. The total charge density is thus given by ρ = -e n_e + Z e n_i and the total current density by J⃗ = -e β c n_e ẑ. The electromagnetic fields are then obtained in excellent approximation by using Gauss's law and Ampère's law, as is customary in ultra-relativistic beam analyses. From these fields, the forces on a beam electron and on a plasma ion are obtained from the Lorentz force relation, and the resulting expressions are simplified using the assumptions ṙ≪ c and ż≃β c for beam electrons and ṙ≪ c and ż≪ c for plasma ions. The total force on a beam electron is thus


    F⃗_e(r⃗) = -e^2/ϵ_0 r∫_0^r (Z n_i(r) - n_e(r)/γ^2) r dr r̂.


By assuming Z n_i ≫ n_e / γ^2, which holds for GeV beams, this expression can be simplified to


    F⃗_e(r⃗) = -Ze^2/ϵ_0 r∫_0^r n_i(r) r dr r̂.


The force on a plasma ion is, on the other hand, 


    F⃗_i(r⃗) = -Ze^2/ϵ_0 r∫_0^r ( n_e(r) - Z n_i(r) ) r dr r̂.


The Hamiltonian yielding such forces on a beam electron or plasma ion can now be written:


    ℋ_α = p_x^2/2 γ_α m_α + p_y^2/2 γ_α m_α - ∫ F_α, r(r) dr,


where the subscript α indicates either e for beam electron or i for plasma ion. Note that γ_e = γ and γ_i ≃ 1 since the ions are assumed to move non-relativistically. To obtain the equilibrium distribution function f_α (q⃗, p⃗), rather than solving the Maxwell-Vlasov equation directly as done in Ref. <cit.>, we exploit the fact that this thermal (globally uniform temperature) form can be generally written as f_α(q⃗, p⃗) ∝ e^-ℋ_α(q⃗, p⃗)/τ_α where τ_α is the transverse temperature of the species. Thus


    f_α(x, y, p_x, p_y) = C_α e^-p_x^2 + p_y^2/2 γ_α m_ατ_α + 1/τ_α∫ F_α, r(r) dr.


 We now explicitly place f_α in in separable form by performing a coordinate transformation to yield


    f_α(r, θ, p_x, p_y)    = λ_α f_r, α(r) f_θ, α(θ) f_p_x, α(p_x) f_p_y, α (p_y) 
    
    f_r, α(r)    = C_α e^1/τ_α∫ F_α, r(r) dr
    
    f_θ, α(θ)    = (2π)^-1
    
    f_p_x, α(p_x)    = 1/√(2π)σ_p_x, αe^-p_x^2/2σ_p_x, α^2
    
    f_p_y, α(p_y)    = 1/√(2π)σ_p_y, αe^-p_y^2/2σ_p_y, α^2


where λ_α is the number of particles of species α per unit longitudinal length, σ_p_x, α = σ_p_y, α = √(γ_α m_ατ_α), and f_α is normalized such that


    ∫_0^∞f_r, α(r)r dr = ∫_0^2πf_θ, α(θ) dθ = 
        = ∫_-∞^∞ f_p_x, α(p_x) d p_x = ∫_-∞^∞ f_p_y, α(p_y) d p_y = 1.


The equilibrium condition for each species is obtained from the fact that the radial component off the equilibrium distribution must be proportional to the number density of that species. Mathematically,


    n_α(r) = C_α e^1/τ_α∫ F_α, r(r) dr,


where C_α has been redefined to absorb the constant of proportionality. 



 §.§ Bennett Profile
 

Now we examine the Bennett-type equilibrium obtained from this analysis. Combining Eq. (<ref>) with Eq. (<ref>) and taking α = e gives an integral equation for the beam electron density n_e. Similarly combining Eq. (<ref>) with Eq. (<ref>) and taking α = i gives an integral equation for the plasma ion density n_i. It is straightforward to verify that the densities


    n_e = ρ_e/(1+(r/a)^2)^2


and


    n_i = ρ_i/(1+(r/a)^2)^2


solve this system of integral equations, provided the temperatures are given by


    τ_e = Z e^2 a^2 ρ_i/8 ϵ_0


and


    τ_i = Z e^2 a^2 (ρ_e - Zρ_i)/8 ϵ_0.


Substituting Eq. (<ref>) into Eq. (<ref>) yields the force on a beam electron,


    F⃗_e(r) = -Z e^2 ρ_i/2 ϵ_0r⃗/1+(r/a)^2.


Similarly, substituting Eq. (<ref>) into Eq. (<ref>) yields the force on a plasma ion


    F⃗_i(r) = -Ze^2(ρ_e - Z ρ_i)/2 ϵ_0r⃗/1+(r/a)^2.


Substituting these forces into into Eq. (<ref>) yields the Hamiltonian, which for beam electrons is


    ℋ_e    = p_x^2/2 γ m_e + p_y^2/2 γ m_e + 
       + Z e^2 a^2 ρ_i/4 ϵ_0ln( 1 + (r/a)^2 )


and for plasma ions is


    ℋ_e    = p_x^2/2m_i + p_y^2/2m_i + 
       + Ze^2 a^2(ρ_e - Z ρ_i)/4 ϵ_0ln( 1 + (r/a)^2 ).


The beam electron and ion distribution functions are given in this analysis by


    f_α(r, θ, p_x, p_y) = λ_α×(2a^2/(a^2 + r^2 )^2) ×( 1/2 π) ×
       ×( 1/√(2π)σ_p_x, αe^-p_x^2/2σ_p_x, α^2) ×( 1/√(2π)σ_p_y, αe^-p_y^2/2σ_p_y, α^2)


where


    λ_α = π a^2 ρ_α


and


    σ_p_x, α = σ_p_y, α = √(γ_α m_ατ_α).




 §.§ Modified Bennett Profile
 

As discussed in the following section, PIC simulations show that while the beam density is accurately described by a Bennett profile, the ion density is not so satisfactorily described. Indeed, while the center of the ion column is relatively well approximated by a Bennett profile, the density profile shown in Fig. <ref> exhibits “wings" at a distance of a few Bennett radii off axis as well as a constant background ion density “pedestal" for r ≫ a. While accounting for the wing shape is outside the scope of the following analysis, it is important to explore an analytical description of the resulting changes to the beam profile due to the constant background ion density that extends outside the collapsed region.

We begin our analysis by defining the modified ion Bennett profile as the original ion Bennett profile Eq. (<ref>) plus a uniform background ion density n_0:


    n_i = ρ_i/(1+(r/a)^2)^2 + n_0.


To derive the modified electron Bennett profile we start with the equilibrium Eq. (<ref>) with α = e. Substituting in the force Eq. (<ref>), we obtain an equation for n_e. We make the assumption that the temperature τ_e in this equation is independent of the background ion density n_0 – the ions are approximately in a global thermal equilibrium after the collapse process. 

With this assumption, and since the modified electron Bennett profile should reduce to the unmodified electron Bennett profile Eq. (<ref>) for n_0 → 0, τ_e can be determined. The modified electron Bennett profile is thus given by


    n_e = ρ_e e^-r^2/2 σ_^2/(1+(r/a)^2)^2


where


    σ_ = a/2√(ρ_i/n_0)


and


    τ_e = Z e^2 a^2 ρ_i/8 ϵ_0.


The modified Bennett profiles, in contrast to the original Bennett form, are only an approximate equilibrium because they do not self-consistently solve for the ion equilibrium. However, as the errors in the modified description are quantitatively scaled to the undisturbed ion density, which is small compared to the final peak ion density. Thus the error arising from this approach should be negligible.  To further provide justification of this assertion, we point to the fact the the form of the modified ion Bennett profile Eq. (<ref>) agrees well with the results of the PIC simulations (except for the distribution wings), as shown in Fig. <ref>.

The existence of the ion density pedestal in equilibrium results in a major conceptual benefit: while Eq. (<ref>) possesses a divergent second moment ⟨ x^2 ⟩,  Eq. (<ref>) yields a bounded second moment, and the rms emittance of the electron beam in the modified Bennett case is no longer infinite. This is also a helpful property for the tracking code discussed in Sec. (<ref>). If the rms emittance diverges, the emittance computed by the tracking code yields an unphysical dependence on the number of particles tracked. Despite this important conceptual difference, the radial dependence of modified Bennett distribution is nearly indistinguishable in practice from the Bennett profile, as illustrated in Fig. <ref>. This is owed to the fact that only the  tails of the distribution are affected by the ion pedestal. In this regard, the feature that the force is linear in r far off axis rather than decaying as r^-1 makes phase space diffusion less of a concern. Such diffusion processes are discussed at length in Sec. (<ref>).



 §.§ Determination of Bennett Radius


The goal of this section is to derive expressions for the three parameters a, ρ_e, and ρ_i of the beam electron and plasma ion Bennett distributions Eqs. (<ref>,<ref>) derived from the initial conditions of the system. Three equations are required to determine these three parameters.

The first expression comes from the conservation of the number of beam electrons per unit longitudinal length λ_e during the process of collapse. This, along with Eq. (<ref>) yields 


    λ_e,  = π a^2 ρ_e.


Similarly, the second equation comes from the fact that the number of plasma ions per unit longitudinal length should be constant during the process of collapse. However since the plasma ions are initially uniformly distributed, the initial number of ions per unit longitudinal length is, before considering which ones may be involved in the interaction, infinite. To resolve this issue, we include as relevant to the analysis only plasma ions with an initial distance from the axis less than some value r_i. The initial number of ions per unit length is thus π r_i^2 n_0 and from this along with Eq. (<ref>), we obtain 


    n_0 r_i^2 = ρ_i a^2.


The value of r_i is estimated through the following considerations. The equation of motion for a plasma ion can be obtained utilizing the force given by Eq. (<ref>):


    r⃗” + Z e^2 (ρ_e - Z ρ_i)/2 ϵ_0 m_i β^2 c^2r⃗/1 + (r/a) ^2 = 0,


where derivatives are respect to z = β c t. For r ≪ a the motion is simple harmonic with angular wave-number


    k_i^2 = Z e^2 (ρ_e - Z ρ_i)/2 ϵ_0 m_i β^2 c^2.


We use this frequency to quantify the time needed for nonlinear phase mixing that yields the equilibrium. For this purpose, we take as distance for ion equilibration to be approximately two linear ion oscillation periods. Further, We take r_i to be as the maximum initial radius of an ion that can fall to a final radius of one Bennett radius in this distance. Assuming negligible angular momentum and initial radial velocity, an expression for r'(z) of a particle starting at r = r_i can be derived from Eq. (<ref>) employing conservation of energy. Defining r̃≡ r / a, r̃_i ≡ r_i / a, and z̃≡ k_i z, this expression becomes


    dr̃/dz̃ = √(ln( 1 + r̃_i^2/1 + r̃^2)).


Since r̃ goes from r̃ = r̃_i to r̃ = 1 in a distance z̃ = 4π,


    4 π = ∫_1^r̃_i1/dr̃/dz̃ dr̃ = ∫_1^r̃_i( ln( 1 + r̃_i^2/1 + r̃^2) ) ^ -1/2 dr̃.


Solving this expression numerically we obtain r̃_i = 10.24.

The third and final equation is obtained by assuming there is no dilution of the distribution density in the center of the beam's transverse phase space during the process of collapse to equilibrium. This condition will be enforced in the case of the stronger condition that there is negligible transverse emittance growth during the establishment of the equilibrium. While the lack of distribution density dilution is clearly not always a valid assumption, it has been shown in Ref. <cit.> that it is possible to match the initial beam size so that emittance growth due to collapse is nearly ignorable. With negligible beam-core phase space dilution, the value of the distribution function f(r⃗, p⃗) at r⃗ = p⃗ = 0 is approximately conserved. If we take the initial beam spatial distribution to be described by a cylindrically symmetric bi-gaussian function, the initial value of f(0⃗, 0⃗) is


    f_e, (0⃗, 0⃗) = λ_e, /4 π^2 σ_x^2 σ_p_x^2 = m_e^2 c^2 λ_e, /4 π^2 ϵ_n^2


where we have used the fact that the normalized rms emittance of a cylindrically symmetric bi-gaussian distribution is ϵ_n = σ_x σ_p_x / m_e c. From Eq. (<ref>), the value of f(0⃗, 0⃗) after collapse is


    f_e, final(0⃗, 0⃗) = λ_e, final/2 π^2 a^2 σ_p_x, eσ_p_y, e = 4 ϵ_0 ρ_e/πγ m_e Z e^2 a^2 ρ_i


where we have used Eqs. (<ref>), (<ref>), and (<ref>). To obtain the third and final equation we set Eqs. (<ref>) and (<ref>) equal to eachother to get


    ρ_e/a^2 ρ_i = r_e m_e^4 c^4 Z γλ_e, /4 ϵ_n^2.


Combining Eqs. (<ref>), (<ref>), and (<ref>) we obtain


    ρ_e = λ_e r̃_i/2 ϵ_n√(r_e γ Z n_0/π)



    ρ_i = n_0 r̃_i^2



    a = ( 4 ϵ_n^2/π r_e γ Z n_0 r̃_i^2)^1/4


where


    r̃_i ≈ 10.24


This is a physically intuitive result, as it is similar to the expression for the square of rms beam size σ^2 in the case of linear focusing; this spot area is proportional to the the rms emittance ϵ_n/γ multiplied by the inverse of the linear-focusing betatron wave-number k_β. In the case of ion focusing, k_β=√(2π r_e Z ρ_i/γ), and the scaling a∝σ∝√(ϵ_n/γ k_β) is manifested in Eq. (<ref>).



§ ION COLLAPSE SIMULATION: DEMONSTRATION OF NEAR-EQUILIBRIUM DISTRIBUTIONS
 








PIC simulations were performed in order to verify that an approximate equilibrium is reached in which the beam electron and plasma ion transverse densities are described by a Bennett-type profile after ion collapse. Simulations were performed using QuickPIC, a collisionless 3D parallel quasi-static particle-in-cell code <cit.>. Because the collapse to near-equilibrium is very fast compared to the full time scale of the PWFA interaction, we use the PIC simulations only to examine this short (few 10's of femtosecond) transient period, leaving longer time-scale issues such as emittance growth due to scattering after near equilibrium has been reached for the tracking code studies described below. Similar to the approach of <cit.>, the ion channel is created prior to the arrival of the beam involved in the collapse dynamics. 

The beam and plasma parameters of these ion collapse PIC simulation are given in Table <ref> and are based on what is likely achievable in FACET-II experiments with beam optics solutions using permanent magnetic quadrupoles <cit.>. The plasma density and bunch length are chosen to give ample interaction time for near-equilibrium to develop, with  an estimate of k_iσ_z≃ 4.2. After the collapse, the value of k_iσ_z is increased further, to above 25, as seen from the results indicated below. 

While we have developed a theoretical framework to understand the dual Bennett-type equilibrium in this system, the context of this present work is also experimentanl. In this regard, the final matching of the beam to the extremely strong ion focusing after collapse is experimentally very challenging in the FACET-II context due to final focus limitations. As such, in the simulations the  initial beam spot size was matched to the linear focusing scenario, a condition that is experimentally in reach.  This of course implies that there will be emittance growth present in this simulation which can be avoided with more careful matching. Because the beam is not nonlinearly matched, the values of ρ_e, ρ_i, and a in the simulation cannot be compared to the theoretical expressions in Eqs. (<ref>), (<ref>), and (<ref>) – one may point to their consistency given the degree of observed emittance growth. Further,  matching the beam to mitigate emittance growth during ion collapse using the method developed in Ref. <cit.> requires a very time-intensive optimization, for the simulation. As such, the comparison to nonlinear matched simulations awaits a future investigation undertaken with enhanced computational resources. 

The nearly equilibrated beam electron and plasma ion densities inside the PWFA bubble at the end of the simulation relevant to FACET-II are shown in Fig. <ref>. The beam and ion density line-out profiles at ξ = 150 μ m, The ion density line-out clearly displays the oscillation period expected, near 60 μm. This figure also shows the fits of the radial distribution functions of the ions and electrons to both the unmodified Bennett profile discussed in Sec. (<ref>) and the modified Bennett profile discussed in Sec. (<ref>).  The fits were performed using nonlinear least squares fitting. The parameters of the fit to modified Bennett profiles as functions of ξ are shown in Fig. <ref>. 

It can be seen from these figures that the beam electron distribution resembles both the unmodified and modified Bennett distribution, which are nearly identical. The plasma ion distribution, unlike the beam electron distribution, displays some notable deviations from the (modified) Bennett equilibrium. This is because, even though we have not perfectly matched the beam to the final size, the beam electrons begin the interaction as a spatially gaussian, thermal distribution, which is quite similar to a Bennett distribution up to a relatively modest difference in spatial extent due to use of the linear matching conditions.  The initial state of the plasma ions is a cold, uniform distribution, however, which is very far from the near-equilibrium formed by fast collection (collapse) and nonlinear phase-mixing. It is observed in Fig. <ref> that the center of the ion distribution follows a modified Bennett profile fairly well while the residual effects of the non-thermal initial conditions are displayed as artifacts in the wings found at r> a. These artifacts can also be seen as a wake effect in Fig. <ref>.

The maximum ion density after collapse in the simulation is found to be ≈ 50 times the initial background ion density, a number consistent with theoretical predictions, which state that up a factor of 100 enhancement of density is expected in a nonlinear-matched scenario. The slightly smaller observed increase in density indicates that the final beam size is slightly increased due to phase space filamentation in the unmatched case that is to be experimentally encountered. Finally, we note that the large ion density present after collapse may give rise to significant emittance growth due to scattering, an issue we examine quantitatively below.



§ PARTICLE TRACKING IN BENNETT PROFILE-DERIVED FIELDS
 



 §.§ Code Overview


The modified Bennett equilibrium formed by the ions generates electrostatic fields with strongly nonlinear radial dependence. In addition, the ions are very dense, well in excess of atmospheric gas density. In order to establish the combined effects of these two phenomena encountered by beam electrons, simulations are utilized. The goal of these simulations is to quantify the emittance growth of the electron beam due to both multiple Coulomb scattering of beam electrons off of plasma ions and diffusion driven by plasma density fluctuations

Because emittance growth from these two sources take place on a timescale orders of magnitude longer than the ion collapse, a process already computationally intensive to simulate with PIC codes, it is not feasible to use a PIC approach to simulate these effects. Instead we have developed a custom 2D tracking code. When a simulation is run with this code, the initial positions and momenta of the beam electrons are randomly sampled from the modified Bennett distribution function. These electrons are then symplectically tracked through the analytically known electromagnetic fields which are derived assuming the system is in equilibrium. At the end of each step, for each particle, an angle is sampled from the scattering angle distribution predicted by Molière's theory, and the particle's momentum is deflected by this angle. Statistics describing the electron distribution are computed and saved at multiple z positions for later analysis. Acceleration is accounted for  by adiabatically adjusting the beam energy and parameters which scale with the energy at each step. A sinusoidal plasma density modulation or gaussian noise in the ion density can be included to drive parametric resonances. The tracking simulation code is written in C++ and is parallelized using the Boost MPI library. Analysis and plotting is performed with Python 3. Simulations are run on UCLA's Hoffman2 computing cluster. The code is available on Github [Equilibrium Plasma Accelerator Multiple Scattering Simulator (EPAMSS): <https://github.com/clairehansel3/EPAMSS>.].



 §.§ Tracking: Macroscopic Forces
 

As noted above, the PIC approach to understanding the behavior of the beam-plasma interaction is very useful for shorter time scales, but to study the longer term evolution of the beam, another approach based on tracking in using model Hamiltonian may be employed. This is analysis is intended to reveal slower, yet critically important phenomena such as phase space diffusion. We now derive the equations of motion used in the tracking code. To obtain the force on a beam electron, we begin with the modified plasma ion Bennett profile Eq. (<ref>). In order to explore parametric resonances we wish to vary n_i as a function of z in our dynamics model. This is done by multiplying n_i by a modulation factor m(z). The equations of motion are then easily obtained from Eq. (<ref>). Introducing normalized variables: t̃≡ω_p t = z̃≡ k_p z, r̃⃗̃≡ k_p r⃗, r̃≡ k_p r, ρ̃_̃ĩ≡ k_p^-3ρ_i, and ñ_̃0̃≡ k_p^-3 n_0, and assuming as before ultra-relativistic (β≃ 1) and paraxial (p_≪ p_z) electron motion, the equations of motion are


    d^2 r̃⃗̃/d z̃^2 + 1/γd γ/d z̃dr̃⃗̃/d z̃ +Z m(z̃)/2 γ(1 + ρ_i/n_0/1+ (r̃/ã)^2) r̃⃗̃ = 0.


Our intention is to use the symplectic velocity-Verlet method to solve the equations of motion numerically. However, due to the damping term in the above expression, the velocity-Verlet method cannot be employed. To remedy this we perform the coordinate transformation r̃⃗̃ = γ^-1/2ũ⃗̃ which gives


    d^2 ũ⃗̃/d z̃^2   + [ 1/4 γ^2 ( d γ/d z̃ )^2 - 1/2 γd^2 γ/d z̃^2 +  . 
       +  . Z m(z̃)/2 γ(1 + ρ̃_̃ĩ/ñ_̃0̃/1+ 1/γ(ũ/ã)^2) ] ũ⃗̃ = 0


which is now in the correct form to be solved numerically using the velocity-Verlet method. Constant acceleration is accounted for by adiabadically updating γ and a with z. Since a scales as γ^-1/4, for constant acceleration we can write


    γ(z̃) = γ_0 + dγ/dz̃z̃


and 


    ã(z̃) = ã_̃0̃(1 + 1/γ_0dγ/d z̃z̃) ^ -1/4.


where dγ/dz̃ is constant. Validation of the tracking algorithm is discussed in Appendix <ref>. 



 §.§ Tracking: Coulomb Scattering
 

As the subject of relativistic particle Coulomb scattering in an ion environment is not commonly discussed, we discuss here the details of the computational approach to modeling the relevant scattering dynamics. The Coulomb scattering algorithm used in the tracking simulation is based on Molière's theory of small angle multiple scattering <cit.>. At the end of each computational step after a particle's phase space state has been tracked, a scattering angle is sampled from the scattering angle distribution which is given by Molière's theory. The velocity of the particle is then deflected by this angle. The step size and scattering angle are assumed to be small enough that the contribution of the scatter to the particle's position can be ignored. Because Molière's theory requires constant transverse density, it is assumed that the step size Δ z is small enough so that the plasma ion density along the particle's trajectory is constant. Note that Molière's theory assumes no ion recoil takes place during scattering.

Molière's original theory is based off a differential cross section which includes atomic screening, and so it must be modified to account for scattering off ions, which lack atomic screening. For this purpose, we use the Mott-Born differential cross section <cit.>.  Because the total cross section for the Mott-Born differential cross section diverges, a minimum angle θ_min is introduced and the cross section is assumed to vanish for θ < θ_min:


    dσ/dΩ(θ) = 
    Z^2 r_e^2 (1 - β_0^2sin^2(θ/2))/4 γ_0^2 β_0^4 sin^4(θ/2)   θ≥θ_min,
    
        0    θ < θ_min.


The code computes θ_min from the total cross section σ which is defined as σ≡π r_σ^2 where the cross sectional radius r_σ is an adjustable parameter of the simulation. Modifying the Molière theory by using the Mott-Born differential cross section (Eq. <ref>) gives the scattering angle distribution, with the following relations applied:


    f(ϑ) = 2ϑ e^-ϑ^2 + ∑_n=1^∞f^(n)(ϑ)/B^n;



    f^(n)(ϑ) = ϑ/n!∫_0^∞ e^-η^2/4( η^2/4ln( η^2/4) )^n J_0(ηϑ) η dη;



    B = -W_-1( -e^-2(1 - γ_E)/Ω_0);



    θ = (θ_min√(Ω_0 B)) ϑ;



    θ_min = 2 Z r_e/γ_0 β_0^2 r_σ;



    Ω_0 = n_i σΔ z.


Here W_-1(x) is the lower branch of the Lambert W-function <cit.>, γ_E is the Euler-Mascheroni constant, and n_i is the density of plasma ions at the position of the particle, Δ z is the step size of the simulation, and Z is the ion atomic number. The function f(ϑ) is normalized such that ∫_0^∞ f(ϑ) dϑ = 1. The results stated above are derived in Appendix <ref>. Before a simulation is initiated, a table of values for f^(n)(ϑ) is created by evaluating this function at evenly spaced points between zero and a cutoff value ϑ_max. The evaluation is performed by using an adaptive Gauss-Kronrod quadrature to numerically integrate (<ref>). This permits  f^(n)(ϑ) to be evaluated during the simulation by using cubic B-spline interpolation on the stored tables. The simulation truncates the sum in Eq. <ref> at a maximum order which is specified as a simulation parameter. In order to evaluate scattering at an angle θ, the code samples the scaled angle ϑ from Eq. <ref> using rejection sampling on an initial uniform distribution from 0 to ϑ_max, and then the angle θ is computed using Eq. (<ref>). If the particle is far enough off axis, the number of scatters is small enough that Molière's theory does not apply. As discussed in Appendix <ref>, Molière's theory requires Ω_0 ≳ 25. If a particle is far enough off-axis such that this condition is not met, the effects of scattering for that particle are assumed to be negligible and ignored. To determine the effects of scattering, the simulation is run twice with the same starting particles, once with scattering turned on and once with scattering off. The validation of the scattering algorithm is discussed in Appendix <ref>.



 §.§ Tracking: Results
 









The tracking code was first run with parameters based on what is achievable at FACET-II, although the plasma was simulated to be 100 m to have more distance to permit observation of possible emittance growth. The tracking run parameters are shown in Table <ref>. The observed emittance growth is plotted in Fig. <ref>. From this plot it is clear that the emittance growth due to scattering in this particular case is negligible.

The tracking code was then run with parameters associated with a hypothetical LC. These parameters are shown in Table <ref>. The emittance growth is plotted in Fig. <ref>. From this plot we can see that the emittance growth in this case is small but non-negligible.



§ CHAOTIC BEHAVIOR AND PARAMETRIC RESONANCES



While scattering is of primary importance in giving phase space diffusion and emittance growth, the nonlinearity of the restoring force can give rise to other relevant processes. In particular, one may excite parametric resonances due to periodic changes in the focusing ion distribution population from plasma density fluctuations or breathing of the bulk beam-ion distribution. In this regard, we note that there is a range of effective locally resonant frequencies due to the non-uniformity of the ion distribution, ranging from √(2 π r_e Z n_0 / γ) at large amplitude, to a near-axis frequency of √(2 π r_e Z (ρ_i + n_0) / γ). In order to evaluate possible phase space diluting effects, we have again employed the tracking code.



To demonstrate resonances and island formation, simulations were performed in which a sinusoidal modulation was added to the ion density. This was done by employing the modulation function m(z) as defined in Sec. (<ref>) to


    m(z) = 1 + A_sin(k_ z).


For these simulations, a single particle was tracked through phase space. This particle was initialized at x = 5a, y =0 and p_x = p_y = 0 at the beginning of the plasma. Simulations were performed with the scattering process turned off and with the LC parameters shown in Table <ref> except l_ = 1 m was used. The relative amplitude of the modulation A_ was chosen to be 0.1. Poincare plots were created by sampling the phase space coordinates at values of z which are multiples of the near axis betatron wavelength. Poincaré plots obtained from these tracking simulations show island formation in cases where the modulation wave-numbers k_ were equal to or near 1, 1/2, or 1/4 times the near axis betatron wavenumber. The Poincaré plots for these three cases are shown in Fig. <ref>.



In order to provide a physical seed for observing diffusive emittance growth due to excitation of parametric resonances, white noise was added to the ion density. This was done by at each step z_i sampling from a gaussian with σ = 0.01. This sampled value, added to unity, was used as the value of m(z_i). The LC simulation from Sec. (<ref>) was repeated without scattering or sinusoidal modulation. The emittance growth is this scenario is shown in Fig. <ref>; it is significant but not unmanageably large. It is of similar importance as the scattering-induced diffusive emittance growth. In the future, a Fokker-Planck approach might be used to better understand diffusion.



§ CONCLUSIONS AND OUTLOOK
 

In this paper we added to the analytical understanding of nonlinear equilibria in PWFAs with ion motion. We utilized and extended previous research, using both theoretical and computational approaches, to confirm that these equilibria are described by Bennett-type profiles, and given an estimate of the associated dimensions such as the Bennett radius a. We have, through tracking simulations, quantified the expected emittance growth due to scattering in a hypothetical PWFA based linear collider that includes strong ion motion and found it to be acceptably small. Finally we discussed diffusion due to parametric resonances in the nonlinear focusing fields of the Bennett-type profile, with similar conclusions reached considering possible problems encountered in LC scenarios.  These results thus significantly extend previous results obtained from  PIC simulations that demonstrated emittance mitigation for PWFAs with ion motion <cit.>. The results of the current work give further reason to be optimistic in this regard.

The issue of ion motion has, due to the demanding conditions needed to explore the relevant physics, to date not been addressed experimentally. Our results give context which identifies the important physics to be investigated in firt experiments. Indeed, we have discussed above a scenario that is to be explored at FACET-II by the E-314 experimental collaboration. In this experiment, by utilizing cutting edge methods in photoinjector electron sources and beam preparation, in tandem with very high gradient final focusing permanent magnet quadrupoles, we expect to be able to create the conditions <cit.> for achieving a collapsed-ion equilibrium at the FACET-II interaction point. These experiments should be able to demonstrate propagation of a joint Bennett-type equilibrium.

In these experiments, the plasma needed is within the current state-of-the-art <cit.>. It has a nominal plasma density of 10^18 cm^-3, in order to provide high phase advance for the ions within the beam. This is also accomplished using a slightly relaxed bunch length of σ_z=40 μm. Thus, as seen in the simulations above, we anticipate a Bennett-type equilibrium will develop over nearly the entire beam. This model situation is ideal for studying the physics of ion collapse.

The beam at final focus will have a spot size that is well sub-optical. We will need to utilize both appearance intensity <cit.> and ionization dynamics (as in the E-317 proposal on high field atomic physics at FACET-II, see Ref. <cit.>), In order to resolve these spot sizes (1̃00 nm) one may use advanced imaging methods <cit.>, including coherent imaging approaches <cit.>, as well as  other tools developed at SLAC at previous test facilities, e.g. laser-wire measurements <cit.>. Electron beam-derived forces produce ions having kinetic energy up to the 100 keV range that may be collected to give a signature of ion motion, as well as a  measure of the electron beam size. It is planned to build a short ion transport line to extract the ions away from the electron beam path, and to measure their energies with a compact magnetic spectrometer made from a permanent magnet dipole <cit.>. 

The betatron radiation spectra emitted by the electron beam for the parameters given above show unique signatures of gamma-ray production. For the parameters used in simulation, photons are produced well into the MeV range, and the spectrum is broadened by the particular nonlinear focusing of the Bennett-type profile’s fields. The Lenard-Wiechert field-based simulations are now underway at UCLA <cit.> to provide detailed theoretical predictions of the spectral shape.  Measurement of the MeV-scale double-differential spectrum will be accomplished  by use of an Compton-pair spectrometer under development at UCLA <cit.>.

A central goal of the experiment should also verify the predicted emittance growth (near a factor of two) under ion collapse conditions, which should be due predominantly to the phase space filamentation during Bennett-type near-equilibrium formation. This entails measurement of the electron beam downstream of the plasma source. Sensitive emittance measurement techniques are now being developed for the E-310 Trojan Horse Injection experiment <cit.>, for deducing emittance growth through beam profile measurements downstream of the plasma exit.

In summary, the physics of critically important phenomena in PWFA have been explored in a theoretical and simulational approach placed in the context of both FACET-II and LC experimental concerns. The issues raised here provide for a rich physics study of ion collapse and beam-ion equilibrium formation. Successful exploration of the physical effects involved should give confidence to the developing view of a PWFA-based future linear collider. 



§ ACKNOWLEDGMENTS


This work was performed with support of the US Dept. of Energy, Division of High Energy Physics, under contract no. DE-SC0009914, and the National Science Foundation, under grant no. PHY-1549132, an NSF Science and Technology, the Center for Bright Beams.

This work used computational and storage services associated with the Hoffman2 Shared Cluster provided by UCLA Institute for Digital Research and Education’s Research Technology Group.
 




§ DERIVATION OF MULTIPLE SCATTERING ANGLE DISTRIBUTION



The derivation of the multiple scattering angular distribution is similar to that found in Ref. <cit.> and Ref. <cit.>. The major difference is that the Mott-Born differential cross section with a cutoff at angles below a minimum angle θ_min is used instead of a cross section with atomic screening.

Let θ be the polar angle and β be the azimuthal angle. ϕ_x and ϕ_y are the projected angles in x and y respectively. Assuming scattering angles are small, these angles are related by


    [     ϕ_x       = θcos(β);     ϕ_y       = θsin(β) ].


Let F(θ, β) be a scattering angle distribution written as a function of the polar and azimuthal angles. The scattering angle distribution as a function of the projected angles is given by F_p(ϕ_x, ϕ_y) = F(θcos(β), θsin(β)). These distributions are normalized so that


    ∫_0^2 π∫_0^∞ F(θ, β) θ dθ dβ = 1


and


    ∫_-∞^∞∫_-∞^∞ F_p(ϕ_x, ϕ_y) dϕ_x dϕ_y = 1.


Note that when the small angle approximation is made, the domain of θ is extended from [0, π] to [0, ∞) and the domains of ϕ_x and ϕ_y are extended from [-π, π] to (-∞, ∞). Now let F_p,M(ϕ_x, ϕ_y) be the distribution of scattering angles after M scatters. The distribution of scattering angles after M+1 scatters is given by the convolution




    F_p,M+1(ϕ_x, ϕ_y) = ∫_-∞^∞∫_-∞^∞ F_p,M(ϕ'_x, ϕ'_y) F_p,1(ϕ_x - ϕ'_x, ϕ_y - ϕ'_y) dϕ_x dϕ_y.
 

In Fourier space, this equation is


    F_p,M+1(ξ_x, ξ_y) = F_p,M(ξ_x, ξ_y) F_p,1(ξ_x, ξ_y)


where


    F_p,M(ξ_x, ξ_y) = ∫_-∞^∞∫_-∞^∞ e^i(ξ_x ϕ_x + ξ_y ϕ_y) F_p,M(ϕ_x, ϕ_y) dϕ_x dϕ_y


and


    F_p,M(ϕ_x, ϕ_y) = 1/(2π)^2∫_-∞^∞∫_-∞^∞ e^-i(ξ_x ϕ_x + ξ_y ϕ_y)F̃_p,M(ξ_x, ξ_y) dξ_x dξ_y.


From Eq. (<ref>) it can be seen that


    F_p,M(ξ_x, ξ_y) = F_p,1(ξ_x, ξ_y)^M.


We now define F̃_M(ξ,α) as 


    F̃_M(ξ,α) ≡F̃_M(ξcos(α), ξsin(α)).


Substituting in Eq. (<ref>),


    F̃_M(ξ,α)    = ∫_0^2 π∫_0^∞ e^iξθ(cos(α) cos(β) + sin(α) sin(β) F_p,M(θcos(β), θsin(β) θ dθ dβ
       = ∫_0^2 π∫_0^∞ e^iξθcos(β - α) F_M(θ, β) θ dθ dβ
       = ∫_0^∞ F_M(θ) [∫_0^2 π e^i ξθcos(β - α) dβ] θ dθ
       = 2 π∫_0^∞ F_M(θ) J_0(ξθ) θ dθ




where we have exploited the fact that the angular distribution is cylindrically symmetric and thus F_M(θ, β) = F_M(θ). Note that from Eqs. (<ref>), F̃_M(ξ,α) = F̃_M(ξ) is independent of α. Eq. (<ref>) is a Hankel transform of order zero, and the inverse transform is


    F_M(θ) = 1/2π∫_0^∞F̃_M(ξ) J_0(θξ) ξ dξ.


The projected multiple scattering angle distribution can be written in terms of F_p,M(ϕ_x, ϕ_y), the projected scattering angle distribution for M scatters, and p_M, the probability of a particle scattering M times, as


    F_p(ϕ_x, ϕ_y) = ∑_M=0^∞ F_p,M(ϕ_x, ϕ_y) p_M.


To find p_M, we imagine a particle which travels an infinitesimal distance dz in a material with number density n of scattering targets each with total cross section σ. The probability that the particle scatters is σ n dz. The probability that it does not scatter is 1 - σ n dz. If we now let the particle travel a finite distance z in the, we can split z into N infinitesimal distances of length dz = z / N. p_M is is given by


    p_M    = lim_N →∞NM( σ n z/N)^M ( 1 - σ n z/N)^N - M
       = (σ n z)^M/M! e^-σ n z = Ω_0^M/M! e^-Ω_0


where we have introduced the dimensionless parameter Ω_0 ≡σ n z. Now we Fourier transform both sides of Eq. (<ref>) to give


    F̃_p(ξ_x, ξ_y) = ∑_M=0^∞F̃_p,M(ξ_x, ξ_y) p_M.


Substituting in Eqs. (<ref>) and (<ref>), we obtain


    F̃_p(ξ_x, ξ_y)    = ∑_M=0^∞F̃_p,1(ξ_x, ξ_y)^M Ω_0^M/M! e^-Ω_0
       = e^(F̃_p,1(ξ_x, ξ_y) - 1) Ω_0.


Now using Eq. (<ref>) we have


    F̃(ξ)    = F̃(ξ, α) 
       = F̃_̃p̃(ξcos(α), ξsin(α) 
       = e^(F̃_p,1(ξcos(α), ξsin(α)) - 1) Ω_0
       = e^(F̃_1(ξ, α) - 1) Ω_0
       = e^(F̃_1(ξ) - 1) Ω_0.


The angular distribution after a single scatter F_1(θ) is proportional to the differential cross section. The constant of proportionality is determined by the normalization condition Eq. (<ref>). F_1(θ) is given by


    F_1(θ) = 1/σdσ/dΩ(θ)


where it is evident from the definition of total cross section that F_1(θ) obeys the normalization condition Eq. (<ref>). In summary, the multiple scattering angular distribution is given by


    F(θ) = 1/2π∫_0^∞ e^(F̃_1(ξ) - 1) Ω_0 J_0(θξ) ξ dξ,


where


    F̃_1(ξ) = 2π∫_0^∞1/σdσ/dΩ(θ) J_0(ξθ) θ dθ.


We now evaluate these equations using the Mott-Born differential cross section with a cutoff at angles below a minimum angle θ_min. The differential cross section is


    dσ/dΩ(θ) = 
    Z^2 r_e^2 (1 - β_0^2sin^2(θ/2))/4 γ_0^2 β_0^4 sin^4(θ/2)   θ > θ_min
    
        0    θ < θ_min,


where Z is the atomic number of the ions and r_e is the classical electron radius. The total cross section is


    σ   = ∫_0^2π∫_0^πdσ/dΩ(θ) sin(θ) dθ dβ .
       = π r_e^2 Z^2/β_0^4 γ_0^2[( θ_min/2)^2 + ln( sin ( θ_min/2) ) ], 
       ≃4π r_e^2 Z^2/β_0^4 γ_0^2 θ_min^2,


where we have assumed θ_min≪ 1. We write the total cross section in terms of the cross-sectional radius r_σ≡√(σ / π). In terms of r_σ, the minimum angle is


    θ_min = 2 Z r_e/γ_0 β_0^2 r_σ.


The series expansion of dσ/dΩ(θ) about θ = 0 is


    dσ/dΩ(θ) ≃ 
    r_σ^2 θ_min^2
    1/θ^4 + 2 - 3β_0^2/12 θ^2 + O(θ^0)   θ > θ_min
    
        0    θ < θ_min.


From Eqs. (<ref>) and (<ref>),


    F̃_1(ξ)    = 2π/σ∫_0^∞dσ/dΩ(θ) J_0(ξθ) θ dθ
       = 2θ_min^2 ∫_θ_min^∞(1/θ^3 + 2 - 3β_0^2/12 θ) J_0(ξθ) dθ
       = 2θ_min^2 ξ^2 ∫_ξθ_min^∞J_0(u)/u^3 du 
       + (2-3β_0^2)θ_min^2 ξ/6∫_ξθ_min^∞J_0(u)/u du 
       = 2x^2 I_1(x) + (2-3β_0^2)θ_min/6 x I_2(x)


where x ≡θ_minξ and the integrals I_1(x) and I_2(x) are defined as


    I_1(x) ≡∫_x^∞J_0(u)/u^3 du,


and 


    I_2(x) ≡∫_x^∞J_0(u)/u du.


Integrating by parts twice and using the relations (J_0(x))' = -J_1(x) and (x J_1(x))'= x J_0(x), I_1(x) can be written


    I_1(x) = J_0(x)/2x^2 - J_1(x)/4x - 1/4I_2(x).


Eq. (<ref>) now becomes


    F̃_1(ξ)    = J_0(x) - x/2 J_1(x) 
       + ((2-3β_0^2)θ_min x/6 - x^2/2) I_2(x).


Eq. (<ref>) can be simplified to


    I_2(x)    = 2/π∫_x^∞∫_0^π/2cos(u cos(ϕ))/u dϕ du 
       = 2/π∫_0^π/2∫_x cos(ϕ)^∞cos(v)/v dv dϕ
       = -2/π∫_0^π/2Ci(x cos(ϕ)) dϕ
       = -[ln(x e^γ_E/2) + ∑_n=1^∞(-1)^n x^2n/2 n 4^n (n!)^2]


where γ_E = 0.57721... is the Euler-Mascheroni constant and where we have used


    J_0(x) = 2/π∫_0^π/2cos(x cos(ϕ)) dϕ,



    Ci(x) = -∫_x^∞cos(u)/u du,


and


    Ci(x) = γ + ln(x) + ∑_n=1^∞(-x^2)^n/2n (2n)!.


From Ref. <cit.> and Ref. <cit.> the important contribution to the integral Eq. (<ref>) occurs when x ≲ 1. Substituting the series expansions of J_0(x), J_1(x), and I_2(x) in Eq. (<ref>) and ignoring the O(x^2) and larger terms, and using the fact that θ_min≪ 1, we obtain


    F̃_1(ξ) = 1 + x^2/2ln(x e^γ_E/2) + O(x^2).


Let η≡√(Ω_0 B) x where B is a constant defined later.


    F̃_1(ξ) - 1    = η^2/4 Ω_0 Bln( η^2  e^2 γ_E/4 Ω_0 B e^2) 
       = η^2/4 Ω_0 B( ln( e^2 γ_E/Ω_0 B e^2)  + ln( η^2/ 4) )


If we choose B such that


    1/Bln( e^2γ_E/Ω_0 B e^2) = -1


or equivalently


    B = -W_-1( -e^-2(1 - γ_E)/Ω_0)


where W_-1 is the lower branch of the Lambert W function, Then by combining Eqs. (<ref>), (<ref>), and (<ref>), we have


    F̃(ξ)    = e^- η^2/4  + η^2/4 Bln( η^2/ 4)
       ≃ e^-η^2/4∑_n=1^∞1/n!(η^2/4 Bln( η^2/4))^n


Plugging this into Eq. (<ref>) and defining  ϑ≡θ / (θ_min√(Ω_0 B)), the probability density function is given by 


    f(ϑ) = 2ϑ e^-ϑ^2 + ∑_n=1^∞f^(n)(ϑ)/B^n;


where 2π f(ϑ) ϑ dϑ = 2π F(θ) θ dθ and


    f^(n)(ϑ) = ϑ/n!∫_0^∞ e^-η^2/4( η^2/4ln( η^2/4) )^n J_0(ηϑ) η dη.




§ EQUILIBRIUM SIMULATION CODE TRACKING VALIDATION
 



By definition, if a beam is sampled from the equilibrium distribution function and tracked according to the equations of motion, it will remain distributed according to the equilibrium distribution function. Thus the tracking code can be validated by comparing the distribution of beam after it has been tracked to the theoretical distribution function. Including acceleration in the tracking allows the gamma scaling of the equilibrium distribution parameters to be checked.

In Fig. <ref>, we compare the equilibrium distribution function to the simulated particle distribution at the end of the plasma in the linear collider simulation described in Sec. (<ref>). Because the simulated distribution matches the theoretical distribution in the figure, we can conclude that both the tracking algorithm correctly tracks particles through accelerating fields.



§ EQUILIBRIUM SIMULATION CODE SCATTERING VALIDATION
 








The first step taken to validate the scattering algorithm was to check that the code correctly samples the scattering angles from the distribution given by Eq. (<ref>). This is necessary due to the complexity of the algorithm which is described in Sec. (<ref>). First, to ensure the values of Eq. (<ref>) computed by the code were correct, they were compared to values computed using Wolfram Mathematica. These were found to be in agreement for various values of ϑ. The next step was to check that the sampled angles are distributed according to the probability density function Eq. (<ref>). This was done by arbitrarily choosing B = 10 and sampling 1,000,000 values of ϑ and plotting a histogram of these values against the probability density function. This is shown in  Fig. <ref> and it is clear that the sampled values are distributed according to the correct probability density function.

For the second test of the scattering algorithm, the code was modified to track electrons initialized on axis with zero transverse momentum through a uniform hydrogen gas of density n_0 with no forces due to electromagnetic fields. To zeroth order in this situation, Molière's theory predicts that the projected angles ϕ_x and ϕ_y are normally distributed with standard deviation given by <cit.>:


    σ_ϕ = 26.6/γ√(m_p n_0 z/X_0)( 1 + 0.038 ln(m_p n_0 z/X_0) )


where X_0 is the radiation length which is 630.4 m^2/kg for hydrogen. Besides those such as the integration tolerance used in the calculation of f(ϑ), there are two unphysical parameters in the scattering algorithm: the cross sectional radius r_σ and the number of steps. First the modified code was run with 1000 particles at 10 GeV tracked in 10000 steps through a 100 m plasma of density 10^20cm^-3 with the cross sectional radius scanned from 10^-10m to 10^-7m. The simulation results are plotted against theory in Fig. <ref>. From these results we can see that the scattering algorithm is more or less in agreement with the theoretical estimate in Eq. (<ref>) for r_σ≈ 10 nm. Next the same simulation was run except this time the cross sectional radius was fixed at 10 nm and the number of steps was scanned from 10 to 1,000,000. The simulation results are plotted against theory in Fig. <ref>. From this plot we can see that the growth in σ_θ was only slightly affected by changing the number of steps by four orders of magnitude. The number of steps in the simulations performed for the FACET-II and LC cases are approximately 3,000,000. From the results shown in Fig. <ref>, we can say the scattering algorithm is a mild underestimate of the scattering predicted by Eq. (<ref>). However not too much faith should be vested in Eq. (<ref>) as it is a zeroth order approximation.



