Plasma-based acceleration (PBA) is being considered for a next generation linear collider (LC). In some PBA-LC designs for the electron arm, the extreme beam parameters are expected to trigger background ion motion within the witness beam, which can lead to longitudinally varying nonlinear focusing forces and result in an unacceptable emittance growth of the beam. To mitigate this, we propose to use quasi-adiabatic plasma density ramps as matching sections at the entrance and exit of each stage. We match the witness electron beam to the low density plasma entrance, where the beam initially has a large matched spot size so the ion motion effects are relatively small. As the beam propagates in the plasma density upramp, it is quasi-adiabatically focused, and its distribution maintains a non-Gaussian equilibrium distribution in each longitudinal slice throughout the process, even when severe ion collapse has occurred. This only causes small amounts of slice emittance growth. The phase mixing between slices with different betatron frequencies leads to additional projected emittance growth within the acceleration stage. A density downramp at the exit of an acceleration section can eliminate much of the slice and projected emittance growth as the beam and ion motion adiabatically defocuses and decreases, respectively. Simulation results from QuickPIC with Azimuthal Decomposition show that within a single acceleration stage with a 25 GeV energy gain, this concept can limit the projected emittance growth to only ∼2% for a 25 GeV, 100 nm emittance witness beam and ∼20% for a 100 GeV, 100 nm normalized emittance witness beam. The trade-off between the adiabaticity of the plasma density ramp and the initial ion motion at the entrance for a given length of the plasma density ramp is also discussed.

Research on plasma-based acceleration (PBA) has made great progress during the past few decades.1 Much of this effort has been motivated to build a considerably less expensive and more compact next generation linear collider (LC) or an x-ray free electron laser (XFEL). In PBA, an intense particle or laser beam is used to form a plasma wake that accelerates a second electron/positron beam that is properly loaded inside the wake. In the blowout regime of PBA,2–5 the drive beam expels all the plasma electrons away, creating a uniform density ion channel with a varying radius behind it. When there is azimuthal symmetry, the wakefield provides an ideal focusing force on an electron beam that is linear in transverse position r and independent of the longitudinal position ξ = c t z, while the accelerating field is independent of r. The transversely uniform accelerating field preserves the beam's slice energy spread, and the longitudinally uniform linear focusing field ensures that all the beam particles (within the same longitudinal slice or across different slices) will oscillate at the same betatron frequency if they have the same energy. These properties remain even when the wakefield is loaded.

For both a collider and XFEL, it is critical that the emittance of the accelerated beam be preserved (kept below a target value). Even in a linear focusing field, preserving the emittance of a beam with a finite energy spread requires matching the beam to the plasma focusing forces. However, the beam's matched spot size in the acceleration stage is typically much smaller than the conventional focusing optics can provide—unless complicated and long beam lines are used, which would lead to greatly reduced averaged gradients. Thus, it is important to design relatively short matching sections at the entrance and exit of a PBA stage. Studies have shown, if the plasma ions are immobile, that by using a tailored short plasma ramp or an adiabatic plasma ramp6–12 it is possible to transition the beam from having a large spot size to having a small matched spot size at the entrance of the acceleration stage, while preserving the beam's emittance.

In order to build a TeV class linear collider, the luminosity, L = f N 2 / 4 π σ x σ y, must be as large as 10 34 cm 2 s 1, where N is the number of particles in each colliding bunch, f is the repetition rate of collisions, and σ x , y is the spot size of the bunch at the interaction point. To achieve such a luminosity, PBA-based linear collider designs have proposed using witness beams with ∼1 nC charge and 100 nm = ϵ n , x ϵ n , y, where ϵ n , x , ϵ n , y are the normalized emittances in each plane.13 For such a beam, the transverse Coulomb field of a matched spot size can pull the ions inward significantly during the transit time of the beam.14,15 In some cases the ions will be attracted by the witness beam and collapse on the axis, forming an ion density spike. The perturbed ions cause a ξ (slice) dependent, transversely nonlinear focusing force, that might result in intolerable emittance growth. It also follows that the accelerating field will vary in r that can might result in energy spread.

Matching the beam's spot size directly to the accelerating stage, which has a uniform density, no longer works well in the presence of ion motion due to the longitudinally dependent and evolving focusing force. The nonlinear nature of the focusing force also causes a Gaussian beam to evolve into a non-Gaussian beam. For typical LC parameters, simulations have shown that the projected emittance grows by around 80% within a few betatron oscillations, and then saturates very quickly.16 

Several ideas have been proposed and studied to mitigate the emittance growth even further in the presence of ion motion using adiabatic processes. Such mitigation is important for multi-stage PBA-LC designs,17 where 80% emittance growth per stage for tens of stages is disastrous. Reference 18 proposed to use a plasma matching section with an adiabatically decreasing ion mass. Reference 19 proposed matching the transverse beam phase space distribution to the nonlinear ion motion-perturbed plasma wakefields slice by slice. However, these ideas may be difficult to realize experimentally. Reference 20 proposed an adiabatic matching procedure where the beam is injected with a low enough energy and thus a large spot size such that ion motion effects are initially small. As the beam is accelerated, the beam adiabatically becomes matched. However, this idea will only work at the initial stage. None of this previous work discussed the use and importance of matching sections at the exit of each stage.

Here, we present a scheme for achievable emittance preservation in the presence of ion motion that is applicable to high energy, high-density electron bunches required in a multi-stage PBA-LC scenario. We use quasi-adiabatic (the ramps do not have to be strictly adiabatic) plasma density upramps and downramps to match the witness beam into and out of a uniform plasma density (plateau) region. We show this method can preserve the witness beam emittance from start to end, even though there is a significant amount of slice dependent ion motion triggered in the plateau. As the beam adiabatically evolves, it naturally finds the slice by slice phase space matching in the plateau, as investigated in Ref. 19. Although there can be small slice and projected emittance growth in the upramp and plateau, we find that much of this can be reversed in a downramp at the exit. We also show that by properly choosing the beam's Courant–Snyder (C-S) parameters at the plasma entrance, the emittance growth can be mitigated even for more general density ramps that are not adiabatic at lower densities.

This paper is organized as follows: In Sec. II, we introduce and review the parameter that quantifies the degree of ion motion. We also introduce a Gaussian phenomenological model for the ion density profile, which will be used to study the evolution of a beam slice. In Sec. III, we introduce the simulation tools used in this paper. In Sec. IV, we study the emittance evolution in a uniform plasma, for cases without and with ion motion. In Sec. V, we study the emittance evolution in a plasma with quasi-adiabatic density ramps using the fully nonlinear quasi-static PIC code, QPAD (QuickPIC with Azimuthal Decomposition),21 for cases without and with ion motion. In this section, the witness beam parameters are close to what we will refer to as FACET II-22 like parameters, so the adiabatic condition of the density ramp is strictly satisfied. In Sec. VI, we use witness beam parameters that are close to LC design parameters. For these parameters, we consider the emittance evolution in both adiabatic density ramps and more realistic, i.e., non-adiabatic density ramps. In Sec. VII, we discuss the trade-off between the adiabaticity of the plasma density ramp and the initial ion motion at the entrance for a given length of the plasma density ramp. We also show simulation results for a 100 GeV witness beam. Finally, conclusions are presented in Sec. VIII.

We note that throughout this paper, we assume the ion motion triggered by the drive beam is negligible, so only the witness beam triggers ion motion. If the drive beam is also adiabatically focused, then the whole process is still adiabatic. We also assume azimuthal symmetry throughout unless otherwise noted.

For a plasma-based (PB) linear collider (LC), the witness beam typically has a very low emittance and spot size, corresponding to a high density with respect to the background density. The resulting large space charge field will induce the motion of background ions significantly during the time it takes the witness beam to pass by an ion. In order to quantify how severe the ion motion is, we follow14,15 and introduce an ion motion parameter, which is defined to be the phase advance of a single ion. To the lowest order, we can estimate the phase advance of an ion using a simple harmonic oscillation model.14,15 We assume that near the axis, the witness beam is simply a uniform cylinder of density n b 0. Since n b 0 n 0 (n0 is the background plasma density) in the ion motion regime, we neglect the self-consistent repelling fields between the ions and only consider the electric field of the witness beam. From Gauss's law the equation of motion for a single plasma ion (ring) initially located inside the beam is
r ( z ) + k i 2 r ( z ) = 0 ,
(1)
where
k i = Z n b 0 e 2 2 ϵ 0 M c 2 = Z n b 0 / n 0 2 M / m k p 0 .
z = ct is the coordinate along the propagation direction, Z is the ion charge state, n b 0 is the witness beam's peak density, ϵ0 is the vacuum permittivity, M is the mass of a single ion, k p 0 = ω p 0 / c is the plasma wavenumber, where ω p 0 = n 0 e 2 ϵ 0 m is the plasma frequency, n0 is the plasma density, e is the electron charge, and m is the electron mass.
In this work, we always use a wide drive bunch with a low peak density to avoid perturbing the background ions. Therefore, here we only consider ion motion triggered by the witness beam. For the slice of the witness beam at the very front, there is no ion motion, and thus the focusing force is linear. However, the space charge forces at the front will begin to pull the ions in, which will affect the later part of the beam. Thus, as we move from the head to tail in the witness beam, the ion motion collapse is stronger. Following Refs. 14 and 15, we can quantify the degree of ion motion through the phase advance parameter Φ for a given slice of the witness beam as follows:
Φ ( ξ ) = k i ξ ,
(2)
where we choose ξ = 0 to be at the head of the beam, so ξ is the longitudinal distance between the beam slice and the beam head. We can also define this ion motion parameter for the entire witness beam to be the ion motion parameter evaluated at the tail of the witness beam,
Φ b = Φ ( L b ) = k i L b = Z n b 0 / n 0 2 M / m k p 0 L b ,
(3)
where Lb is the length of the witness beam. Qualitatively, ion motion is important when Φ b 1. Typically, k p 0 L b 1, so this occurs when Z n b 0 / n 0 2 M / m 1.

It is worth noting that for a beam with a longitudinal Gaussian density profile, the length of the beam Lb is not well defined. However, we can use the effective bunch length L b = 2 π σ z in the equation above to estimate the degree of ion motion.

It turns out Eq. (1) is fairly accurate when the phase advance of the ion is less than π / 2 (before the first ion collapses) even when the response is nonlinear, and it can, therefore, be used to accurately predict the location of the first ion density peak on the axis ( π 2 k i). In Refs. 19 and 20, an Eulerian fluid description was presented for examining ion motion that directly provides a differential equation for the focusing force driven by the space charge of the beam. The ring (Lagrangian) and Eulerian descriptions are essentially identical for linear responses; however, the ring description is still accurate even when strong ion motion has occurred (it fails when rings cross). Both descriptions assume azimuthal symmetry.

For linear collider parameters, Φ b 1, and the first ion density peak resides within the bunch ( π 2 k i < L b). Therefore, the exact transverse profile of the ion density is complicated and longitudinally slice dependent; and it cannot be described by the approaches in Refs. 14, 15, and 19. However, an understanding of the evolution of a beam in the presence of a nonlinear ion response can be obtained through modeling the ion collapse at a given slice as a Gaussian function in r (Ref. 16) (other functions can be used as well),
n ion ( r ) / n 0 = 1 + A 0 exp ( r 2 / 2 σ ion 2 ) .
(4)
The term “1” is included to match the boundary condition as r goes to infinity. We note that this profile does not conserve charge; however, it is reasonable for radii within the beam. This model is parameterized by A0 and σion. The corresponding focusing force is
F r ( r ) = [ r 2 A 0 σ ion 2 1 exp ( r 2 2 σ ion 2 ) r ] m e ω p 0 2 .
(5)
In the later parts of the paper (Secs. IV C and V B 2), we will use this phenomenological model to gain insight for the evolution of a beam slice in a given nonlinear focusing force. A comparison of a Gaussian to the simulation generated ion density profile is given in Sec. IV C.
In order to isolate physics, we study the evolution of a beam slice under a given focusing force using the phenomenological model for ion collapse described in Sec. II B. The equations of motion in the transverse direction for a collection of non-interacting particles in a prescribed force are numerically integrated as follows:
{ x ( z ) = p x ( z ) / γ , p x ( z ) = F x ( x , y , z ) , y ( z ) = p x ( z ) / γ , p y ( z ) = F y ( x , y , z ) .
(6)
The focusing forces Fx and Fy are predetermined, and their form can be specified, based on the problems we are interested in. We generally ignore acceleration (changes in γ). In the absence of ion motion and for azimuthally symmetric wakes the transverse forces are decoupled. For nonlinear focusing forces, the motion in the two transverse planes (x and y) are coupled. Although it is straightforward to include the coupling of the motion in x and y directions, in some cases, to simplify the interpretation and analysis we set y = 0 (where Fy = 0). We also do not push the particles in ξ so each particle remains in a given slice. The numerical integration is done using solve_ivp function from scipy library in Python. The code is also parallelized using mpi4py library. In the later part of the paper, we refer to this method as “single particle simulation” since we simply integrate the equation of motion for each single particle in prescribed fields. We note that a code to study the motion of non-interacting particles in the prescribed fields23 was also recently developed and is being widely utilized.

All the self-consistent simulations are carried out using QPAD, a particle-in-cell (PIC) code that uses the quasi-static approximation (QSA) based on the QuickPIC framework.24,25 In the QSA, it is assumed that the driver evolves on distances much larger than the wavelength of the wake. Thus, the driver is assumed to be static in order to calculate the wakefields. The wakefields are then used to push the beam (driver) forward in z s. To determine the wakefields, there is a “2D” part of the code where the plasma is advanced in ξ on a (x, y) plane assuming z s is fixed, and a “3D” part of the code where the beam is advanced using s. There is thus a 2D “time step,” Δ ξ, and a 3D time step, Δ s. It is fully parallelized through MPI and OpenMP. QPAD uses cylindrical geometry (as opposed to Cartesian geometry as in QuickPIC) with a PIC description in r-z and a gridless description in ϕ. The fields, charge, and current defined on an r-z grid are expanded into azimuthal modes in the gridless ϕ direction, and are truncated at a desired mode number. Particle data are deposited onto the r-z grid for each azimuthal mode based on its ϕ position. This so-called quasi-3D26–28 description is especially efficient for problems with approximate azimuthal symmetry because we only need to keep a small number of modes. When we run with only the m = 0 mode it reduces to the 2D quasi-static model of WAKE29 and the quasi-static version of INF&RNO.30–32 

Modeling beam loading scenarios with matched beams is computationally challenging due to the ion collapse. For such scenarios there is extreme variation in scales between the spot size of the wake structure and the spot size of the witness bunch (and the region of ion collapse). To efficiently simulate this problem it is useful to use mesh refinement,33,34 the quasi-3D approach (QPAD) with high resolution, or mesh refinement on an r-ξ grid.20 For this work we rely on QPAD. Currently, QPAD does not have mesh refinement, so high resolution is used throughout the simulation domain.

Before we investigate the use of matching sections, we investigate how the emittance evolves for an unmatched beam focused on a density plateau with and without ion motion.

To set the stage for matching sections including ion motion, we review the concept of matching in a uniform plasma in the absence of ion motion (linear focusing force). From the beam's envelope equation (longitudinal acceleration is neglected but the formula for the matched spot size remains unchanged if it is included),
σ x = ϵ 2 σ x 3 ( 1 k ¯ β 2 σ x 4 ϵ 2 ) ,
(7)
where ϵ = x 2 x 2 x x 2 is the beam's geometric emittance ( means the ensemble average).
By setting σ x = 0, we can get the matched beam spot size,
σ m = ϵ k ¯ β ,
(8)
where k ¯ β = k p 2 γ ¯, and γ ¯ is the average relativistic factor among all the beam particles. Henceforth we will refer to Eq. (8) as the linearly matched spot size, and matching the beam's spot size to this value is called “linearly matched.”
To explicitly see the dependence on density and energy, Eq. (8) can also be written as follows:
σ m = ( 2 γ ¯ ϵ n 2 k p 0 2 n 0 n ) 1 4 ,
(9)
where ϵ n = 1 m e c x 2 p x 2 x p x 2 is the beam's normalized emittance. For a beam with small energy spread, we have ϵ n γ ¯ ϵ. k p 0 1 is the plasma skin depth corresponding to the reference density n0, and n is the local plasma density. We define a local plasma density, n, to set the stage for considering density ramps. We can see that for small n, σm is large. This means a matched beam's peak density will be low, so the ion motion effects are small. This motivates us to use a plasma density ramp with a low density at the entrance so that the matched beam would trigger negligible ion motion, as we will see in Sec. V B 1.
Throughout the remainder of the paper, we use the Courant–Snyder (C-S) parameters β and α to parameterize the beam, which are defined as follows:
β = x 2 ϵ , α = x x ϵ ,
and α = 1 2 β . The matching condition in a uniform plasma in terms of the C-S parameters is
β = β m , α = 0 ,
(10)
where
β m = σ m 2 ϵ = 1 k ¯ β ,
(11)
is the matched beta function.

For an unmatched beam with energy spread, its emittance will grow due to a chromatic effect;35 particles with different energies rotate at different angular frequencies in the phase space, thus mapping out a larger effective phase space area [Fig. 1(a)]. If the beam is matched, the normalized phase space looks like a circle [Fig. 1(b)], where x and px are normalized to ϵ n i γ i k β and γ i k β ϵ n i, respectively (the subscript i represents the initial value). The phase space area remains the same even if there is phase mixing due to the energy spread, so the emittance is well preserved.8 

FIG. 1.

Phase space plot from a QPAD simulation to illustrate the chromatic effect. (a) Unmatched ( β = 2 β m , α = 0) 10 GeV beam with 2% energy spread, emittance grows due to phase mixing from the energy spread: particles with high energy (blue) and low energy (orange) rotate at different angular frequencies in the phase space. (b) Matched ( β = β m , α = 0) 10 GeV beam with 2% energy spread. Normalized phase space is a circle.

FIG. 1.

Phase space plot from a QPAD simulation to illustrate the chromatic effect. (a) Unmatched ( β = 2 β m , α = 0) 10 GeV beam with 2% energy spread, emittance grows due to phase mixing from the energy spread: particles with high energy (blue) and low energy (orange) rotate at different angular frequencies in the phase space. (b) Matched ( β = β m , α = 0) 10 GeV beam with 2% energy spread. Normalized phase space is a circle.

Close modal

To study the witness beam's emittance evolution in the presence of ion motion, we run self-consistent QPAD simulations. We start with parameters that are not as severe as they would be for a LC design to illustrate some physics. These parameters also overlap somewhat with what might be available at FACET II.22 The background plasma is pre-ionized hydrogen with a uniform density n 0 = 10 17 cm 3. The drive beam's density profile is tri-Gaussian, n b = n b 0 exp ( x 2 2 σ x 2 ) exp ( y 2 2 σ y 2 ) exp ( z 2 2 σ z 2 ), with a symmetric transverse spot size σ x = σ y = 0.8 k p 0 1 = 13.5 μ m, longitudinal bunch length σ z = 2 k p 0 1 = 23.8 μ m, and a peak density of n b 0 / n 0 = 4. This corresponds to 4.4 nC of charge. We assume the drive beam is non-evolving. The drive beam's ion motion parameter is Φ b = 0.12 (calculated using an effective bunch length L b = 2 π σ z) so that the ion motion triggered by the drive beam can be neglected. The witness beam's density profile is transversely bi-Gaussian and longitudinally trapezoidal, with Ihead = 27.2 kA at the head and Itail = 17.6 kA at the tail, so the accelerating field Ez is nearly flattened across the witness beam ( E z 0.56 m c ω p 0 / e).36,37 The length of the witness beam is L = 2 k p 0 1 = 33.7 μ m (so the total charge of the witness beam is 2.5 nC), and its head is located 5 k p 0 1 behind the center of the drive beam. The witness beam has an initial energy of 10 GeV, and an initial normalized emittance of ϵ n i = 1 μ m = 0.06 k p 0 1 for n 0 = 10 17 cm 3. These parameters above are chosen to be the same as in Ref. 20 for convenience. The simulation setup is similar to Fig. 11.

The simulation is performed by using QPAD with only the m = 0 mode (thereby assuming azimuthal symmetry), where we used a ( r , ξ ) simulation box ( ξ = c t z) with 7340 cells in r and 2640 cells in ξ. The resolution is Δ r = 1.36 × 10 3 k p 0 1 and Δ ξ = 5 × 10 3 k p 0 1. The witness beam consists of ×106 numerical particles, which are pushed every 3D time step Δ s = 10 c / ω p 0. Both plasma electrons and ions are initialized at four locations uniformly spaced per r , ξ cell with 16 particles distributed in the azimuthal direction (effectively 64 particles per cell).

First, we consider cases where the witness beam is linearly matched [Eq. (8)] ( σ m = 0.025 k p 0 1 = 0.41 μ m, n b , head / n 0 5.4 × 10 3) with 0 and 2% initial energy spread, respectively (we stress that conventional optics cannot be used to provide such a small spot size). This means the witness beam is unmatched to the nonlinear focusing force triggered by the ion motion ( Φ b 2.2 if we use n b , avg = ( n b , head + n b , tail ) / 2 to estimate the peak density of the beam), leading to emittance growth. This is shown in Fig. 2. We can see that the projected emittance rapidly grows and oscillates during the first few betatron oscillations, after which it asymptotes to a function with a small linear slope. The spot size also eventually asymptotes to a function with a small negative linear slope. The slow growth in emittance has been determined to be due to a combination of acceleration and ion motion. The acceleration leads to an adiabatic decrease in spot size which increases the amount of ion motion. As a result, the shape of the beam's equilibrium profile is modified leading to emittance growth. Benedetti et al.20 also recognized that the slow emittance growth was due to acceleration when studying a different adiabatic matching concept.

FIG. 2.

The emittance (a) and spot size (b) evolution for beams with different initial parameters for QPAD simulations.

FIG. 2.

The emittance (a) and spot size (b) evolution for beams with different initial parameters for QPAD simulations.

Close modal

If the beam is linearly matched, then having an initial energy spread does not lead to extra emittance growth, as we can see from Fig. 2(a).

A possible way to mitigate emittance growth is to focus the beam to a smaller spot size than the linearly matched spot size,16 so that it is closer to some nonlinear equilibrium value after ion motion is triggered. We ran another QPAD simulation with the witness beam's initial spot size to be 85% of the linearly matched spot size and no initial energy spread. In Fig. 2, we see a smaller emittance growth, but now the emittance and spot size oscillate with a larger amplitude. This is because the slices closer to the head of the witness beam experience less ion motion [see Eq. (2)], where the focusing force is just slightly perturbed from the linear focusing force from a fixed ion background. However, the projected spot size is now smaller than for an initially linearly matched spot size. Thus, the spot size oscillates and if there is no energy spread the oscillations will continue for a long time before the weak nonlinearity causes the beam to reach an equilibrium. This leads to the oscillations seen in the projected emittance [green line in Fig. 2(a)].

We also consider a case where the witness beam has the LC relevant parameters used in the simulations in Ref. 16, except here we simulate the entire transverse region of the wake. We note that ϵ RMS ϵ n , x ϵ n , y varies from a few 100 nm to 1 μm in various LC designs.13 Here, we use 100 nm, which is at the low end (which we henceforth refer to as LC parameters), while the FACET II cases use 1 μm, which is at the high end of LC designs.

Due to the smaller emittance, the matched spot size is now smaller and the beam density is higher, which requires finer resolution. Carrying out numerous simulations with sufficient resolution of this type is now possible with the advent of QPAD. The physical and simulation parameters are as follows: The pre-ionized hydrogen plasma density is 10 17 cm 3, the tri-Gaussian drive beam has 3 × 10 10 electrons (4.8 nC), normalized emittance 1  mm, and bunch length σ z = 30 μ m ( 1.8 k p 0 1 ). The drive beam is still set to be moving at the speed of light and non-evolving (it is represented by a specified current profile), and its spot size is set to be linearly matched: σ r = 10 μ m ( 0.6 k p 1 ), so it has a peak density of n b 0 / n 0 = 5.9, which triggers negligible ion motion. The drive beam's ion motion parameter Φ b = 0.18 (calculated with an effective bunch length L b = 2 π σ z). The witness beam has 1010 electrons (1.6 nC), an initial energy of 25 GeV, an initial energy spread of 2%, and an initial projected emittance of ϵ n i = 0.1 μ m ( 0.006 k p 0 1 ). It has a tri-Gaussian density profile with bunch length σ z = 10 μ m ( 0.6 k p 0 1 ). The ion parameter for the witness bunch is Φ b 6. The distance between the centers of the two beams is 115 μ m ( 6.8 k p 0 1 ).

The ( r , ξ ) simulation box has 32 768 × 1024 cells with a resolution of 3.6 × 10 4 k p 0 1 by 1.8 × 10 2 k p 0 1, where ξ = 0 is chosen to be at the longitudinal Gaussian centroid of the witness beam. Both plasma electrons and ions are initialized at four locations uniformly spaced per ( r , ξ) cell with 16 particles distributed in the azimuthal direction. The witness beam consists of 107 numerical particles.

The beam particles are pushed every 3D time step Δ s = 10 c / ω p 0. We only keep the lowest m = 0 azimuthal mode. A few simulations with m = 0, 1, 2 azimuthal modes were carried out, which verified that no hosing seeded by noise occurs with ion motion.

In Fig. 3(a), we observe an ∼80% projected emittance growth from linearly matching the beam to the uniform plasma ( σ r = 0.1 μ m = 0.006 k p 0 1) for these LC relevant parameters. Overfocusing the initial witness beam's spot size to 1/2 of the linearly matched spot size reduces the emittance growth but causes an oscillation in the projected emittance. These results are essentially identical with the results in Ref. 16, giving confidence in QPAD as well as the previous conclusions from the QuickPIC simulations where only a small region near the axis was simulated. However, the QPAD simulations include the self-consistent acceleration of the beam as well as the possibility of the beam hosing.

FIG. 3.

The emittance (a) and spot size (b) evolution for beams with different initial parameters for QPAD simulations. The witness beam has LC parameters.

FIG. 3.

The emittance (a) and spot size (b) evolution for beams with different initial parameters for QPAD simulations. The witness beam has LC parameters.

Close modal

In general, the ion collapse near the axis and the corresponding focusing force is slice dependent and complicated. However, there is essentially no mixing between slices. In this section, we will use single particle simulations to study the evolution of a beam slice. We study the evolution of a beam slice by modeling the density profile of ion collapse using Eq. (4) and begin by choosing A 0 = 100 and σ ion = 0.005 k p 0 1. The beam slice has 10 GeV energy with ϵ n = 0.0594 k p 0 1 ( 1 μ m ), and no initial energy spread. There is also no longitudinal acceleration included in the simulation. We use 105 numerical particles in the beam slice. For simplicity we assume the particles' transverse motion is 1D, where they move under a 1D version of the focusing force in Eq. (5). We simply replace r with x in Eq. (5), which is not strictly valid for a nonlinear focusing force except for y = 0 and py = 0.

If we match the beam with a spot size equal to Eq. (8) (the matched spot size for A 0 = 0), we can see that both the emittance and spot size saturate very quickly within a few betatron oscillations [solid curves in Fig. 4(a)], which is qualitatively similar to what we observed in Fig. 2. The rapid emittance growth at the beginning is due to the nonlinear phase mixing, and Fig. 4(b) shows the spiral phase space caused by phase mixing at an early time of the simulation.

FIG. 4.

(a) Emittance and spot size evolution using a Gaussian phenomenological model. (b) Spiral phase space due to nonlinear phase mixing.

FIG. 4.

(a) Emittance and spot size evolution using a Gaussian phenomenological model. (b) Spiral phase space due to nonlinear phase mixing.

Close modal

We also tried to nonlinearly match the beam's spot size to the steady state spot size of the Gaussian ion model.  Appendix C shows how to estimate the steady state spot size. We can see in Fig. 4(a) that if the beam's initial spot size is closer to the steady state spot size in the nonlinear focusing force, emittance growth can be mitigated significantly. In Ref. 19, matching the beam, including ion motion slice by slice, was self-consistently analyzed for modest amounts of ion collapse.

To make connections to the self-consistent QPAD simulations, we take output from the simulation results presented in Sec. IV B, where the witness beam has LC relevant parameters, for which there is significant ion motion (blue curves in Fig. 3). For comparison, we plot the steady-state transverse ion density lineout at the longitudinal centroid of the witness beam [blue curve in Fig. 5(a)] and the emittance growth of the beam slice at the longitudinal centroid [blue curve in Fig. 5(b)]. We then ran two 2D single particle simulations (the motion in x and y are coupled) with different ion collapse parameters. Both of them have the same on-axis ion density ( 1 + A 0) as the QPAD simulation, but their σion are different. The orange curve in Fig. 5(a) has a smaller σ ion = 4 × 10 4 k p 0 1. While it captures the central part of the ion density lineout from the self-consistent QPAD simulation well, it fails to capture the long tails on both sides, thereby underestimating the slice emittance growth [see the orange curve in Fig. 5(b)]. If we choose a wider σ ion = 10 3 k p 0 1 [green curve in Fig. 5(a)], the core is wider than in the self-consistent result, so the slice emittance growth is overestimated slightly [see the green curve in Fig. 5(b)].

FIG. 5.

Comparison of QPAD simulation and 2D single particle simulations with a Gaussian ion collapse model. (a) Blue: The transverse ion density lineout at the longitudinal centroid of the witness beam (ξ = 0) from QPAD simulation. The dots correspond to the grid points in the simulation. Orange and green: Gaussian ion collapse models with the same on-axis density and different σion. (b) Blue: Emittance evolution of witness beam's centroid slice from QPAD simulation. Orange and green: Emittance of a beam evolving in a fixed focusing force corresponding to the Gaussian ion collapse models shown in (a).

FIG. 5.

Comparison of QPAD simulation and 2D single particle simulations with a Gaussian ion collapse model. (a) Blue: The transverse ion density lineout at the longitudinal centroid of the witness beam (ξ = 0) from QPAD simulation. The dots correspond to the grid points in the simulation. Orange and green: Gaussian ion collapse models with the same on-axis density and different σion. (b) Blue: Emittance evolution of witness beam's centroid slice from QPAD simulation. Orange and green: Emittance of a beam evolving in a fixed focusing force corresponding to the Gaussian ion collapse models shown in (a).

Close modal

Note that it is shown in Ref. 16 (although in 1D) that a larger A0 or a larger σion in the Gaussian ion collapse model will lead to a larger slice emittance growth. Due to the complicated nature of ion collapse, it is not accurately represented by a Gaussian function, especially the “long tail” parts that extend further away from the axis. However, we can still acquire qualitative insight from this model.

In this section, we investigate if emittance growth can be mitigated through the use of an adiabatic density upramp and downramp, placed on each side of a density plateau. The importance of the density downramp at the exit when ion motion is triggered has not been previously emphasized.

In a non-uniform plasma, the matching condition is modified from Eq. (10) and becomes
β = β m , α = α m ,
(12)
where β m , α m are the locally matched Courant–Snyder (C-S) parameters,9,12,38 defined as
β m ( z ) = 2 γ c / ω p ( z )
(13)
and
α m ( z ) = 1 2 β m ( z ) ,
(14)
where z is the propagation distance. Equation (14) is analogous to the relationship between α and β. An adiabaticity parameter can be defined as A = | α m |,9,38 which describes how adiabatic the plasma density is changing from the beam's point of view. The adiabatic condition is given as
A 1.
(15)
If the adiabatic condition is satisfied, it means the betatron frequency k β changes very little within one betatron wavelength. Loosely speaking, the beam feels the plasma density is changing very slowly.

For an adiabatic plasma ramp, it has been shown that if we match the beam to the plasma entrance, the beam will stay matched as it propagates in the plasma, and the emittance can be preserved perfectly even if the beam has a finite energy spread.6,12,39

Here, we show results from a self-consistent QPAD simulation that demonstrates the utility of including upramp and downramp matching sections in the absence of ion motion. The beam parameters are identical to the FACET-II case in Sec. IV B. The plasma density profile is shown in Fig. 6(a). We assume the background plasma is pre-ionized hydrogen plasma with density n 0 = 10 17 cm 3 at the density plateau. The plasma ions are kept fixed. The length of the plateau is 4 × 10 4 k p 0 1 (67.3 cm for n 0 = 10 17 cm 3). Before the density plateau, we use an adiabatic density upramp designed as follows: At the entrance of the plasma upramp, the density is n entrance = 0.01 n 0. In the density upramp, αm [defined in Eq. (14)] is positive. Therefore, we design the ramp with αm starting from α m i = 0.1 at the entrance and linearly decreasing to 0 at the end of the upramp, so that the transition from the upramp to the plateau is smooth. Thus, the ramp is adiabatic throughout. Based on this assumed linear dependence of αm, we can integrate to solve for the density dependence of the upramp (see  Appendix A),
n ( z ) n 0 = 1 ( 1 + α m i ( z L ) 2 β m 0 L ) 2 ,
(16)
where α m i = 0.1 , β m 0 = 2 γ k p 0 1 = 200 k p 0 1 for a 10 GeV beam is βm evaluated at the density plateau, and L = 180 00 k p 0 1 is the length of the upramp. We design the density profile this way so that the adiabatic condition [Eq. (15)] is satisfied everywhere in the upramp, and A is continuous at the end of the upramp. We use a plasma density profile for the downramp, which is completely symmetric to the upramp. The entire plasma density profile is shown in Fig. 6(a). It is worth noting that after being accelerated in the density plateau, the beam's energy is higher in the downramp than in the upramp, so it may be beneficial to carefully design a different downramp that accounts for this energy gain. This will be explored in Sec. VI A. Here, we just choose a symmetric downramp for simplicity. Also, when designing the ramp, we assume there is no longitudinal acceleration and ion motion. In the presence of ion motion (discussed later), the adiabatic condition is not well defined and may not strictly hold due to the enhanced betatron frequency. But in principle, we can always elongate the ramp to make it more adiabatic.
FIG. 6.

(a) Plasma density ramp (blue) with αm (red). (b) Emittance (blue) and spot size (red) evolution for a matched beam in the absence of ion motion.

FIG. 6.

(a) Plasma density ramp (blue) with αm (red). (b) Emittance (blue) and spot size (red) evolution for a matched beam in the absence of ion motion.

Close modal

In the QPAD simulation, we use a ( r , ξ ) simulation box with dimensions 10 k p 0 1 × 12 k p 0 1, and cell size 2.5 × 10 3 k p 0 1 × 10 2 k p 0 1. The witness beam consists of 106 numerical particles. The beam particles are pushed every 3D time step Δ s = 10 c / ω p 0. Both plasma electrons and ions are initialized at four locations uniformly spaced per r , ξ cell with 16 particles distributed in the azimuthal direction. We keep the lowest three azimuthal modes: m = 0, 1, 2. We did this to ensure that for a beam that begins with azimuthal symmetry—there is a lack of hosing growth from noise.40,41 We found that all the quantities for m = 1, 2 indeed essentially vanish.

We match a witness beam with a 2% initial energy spread to the plasma entrance [Eq. (12)]. In Fig. 6(b), we can see that as the witness beam propagates, its spot size evolution strictly follows the locally matched spot size [Eq. (9)], so the beam always stays matched, and its emittance is preserved perfectly.

1. QPAD simulation

In Sec. V A, we showed that in the absence of ion motion, the witness beam emittance can be perfectly preserved using adiabatic plasma density ramps. In this section, we will show that by using the same method, even in the presence of ion motion, the beam emittance can still be well preserved. We consider a simulation with identical parameters to the ones in Sec. V A, except that we now allow the plasma ions to move.

In Fig. 7(a), we show the evolution of the projected emittance (red dashed) and the slice emittance at the head (blue), middle (orange), and tail (green) of the witness beam. The witness beam (with a trapezoidal current profile) is located within the range 0 ξ 2 k p 0 1, where ξ = 0 is at the beam head and ξ = 2 k p 0 1 is at the beam tail. The thickness of the slices is chosen to be Δ ξ = 0.1 k p 0 1. The projected emittance steadily grows in the density upramp, then grows very slowly in the density plateau,20 but eventually decreases in the density downramp, with only 1 % net growth at the exit of the plasma downramp. For the beam slice at the head of the beam (blue, ξ = 0.05 k p 0 1), the slice emittance remains constant, and the slice spot size is just the linearly matched spot size. They are the same as in Fig. 6(b), because the beam slice at the head experiences no ion motion. For slices in the middle (orange, ξ = 1 k p 0 1) and tail (green, ξ = 1.95 k p 0 1) of the beam, the emittance grows in the upramp due to the nonlinear focusing force from the ion motion. However, in the density downramp, the slice emittances decrease. In Fig. 7(b), we can see that the slice spot sizes at the rear of the beam (larger ξ) are also smaller than the linearly matched spot size at the head of the beam (smaller ξ) due to ion motion.

FIG. 7.

(a) Evolution of projected emittance (red dashed) and slice emittance. (b) Evolution of projected spot size (red dashed) and slice spot sizes.

FIG. 7.

(a) Evolution of projected emittance (red dashed) and slice emittance. (b) Evolution of projected spot size (red dashed) and slice spot sizes.

Close modal

To understand the evolution of the projected emittance, we need to investigate how the phase space evolves within individual slices. At the entrance of the plasma upramp, the density is low, and according to Eq. (9), σm is large. The beam's peak density (and self-electric field) is therefore low. As a result, ion motion is initially small [ Φ b 0.69, estimated using Eq. (B3)], and the beam is essentially matched to the unperturbed linear focusing force. In Fig. 8, we present the focusing force at the head ( ξ = 0.05), middle ( ξ = 1.0 k p 0 1), and tail ( ξ = 1.95 k p 0 1) of the witness beam for three propagation distances that correspond to the entrance of upramp, middle of plateau, and exit of downramp. We also show the phase space for slices centered at the same ξ locations (represented by different colored dots). In the upramp region before there is any ion motion, Fig. 8(a) shows that the focusing force at each slice is nearly the same as from the ion column. In Fig. 8(b), the phase space distributions for each slice all overlap as well. Note that x and px values are normalized to their matched values at a given z in the absence of ion motion.

FIG. 8.

(a), (c), and (e) Focusing force in x ̂ at different longitudinal ξ positions for different z propagation distances. (b), (d), and (f) Phase space ellipses corresponding to different longitudinal slices. The colors of each particle correspond to colors of the ξ for the focusing force.

FIG. 8.

(a), (c), and (e) Focusing force in x ̂ at different longitudinal ξ positions for different z propagation distances. (b), (d), and (f) Phase space ellipses corresponding to different longitudinal slices. The colors of each particle correspond to colors of the ξ for the focusing force.

Close modal

As the beam propagates into the ramp, the plasma density slowly increases, causing the beam to be adiabatically compressed. This leads to a higher peak density (thus a stronger self-electric field), so ion motion is adiabatically triggered. In this process, the beam's distribution will slowly and continuously evolve to match the local nonlinear focusing force. Each slice of the beam experiences a different transverse focusing force and, therefore, evolves differently. At the density plateau, the ion motion parameter has increased to Φ b 2.2, so ion motion is significant. We can see the focusing forces at different ξ are different [Fig. 8(c)], leading to different evolution for the phase space ellipse [Fig. 8(d)] of each slice. Due to the nonlinear focusing force there is some emittance growth within each slice. However, in contrast to the growth due to the phase mixing that occurs for initially unmatched beams, this growth is reversible in the downramp.42 Due to the nonlinear focusing force, the distributions of the beam particles are no longer Gaussian. These non-Gaussian profiles lead to larger values for the slice emittance in the density plateau, as shown in Fig. 7(a). We will analyze the slice emittance evolution in greater detail later in Sec. V B 2 using the phenomenological model described in Sec. II B. In addition, slices at the head, middle, and tail of the beam feel different focusing forces, so their corresponding phase space ellipses no longer overlap with each other, resulting in a projected emittance growth that is much larger than for each slice.

In the downramp, the beam's evolution is the reverse of that in the upramp; the beam's spot size slowly expands so the peak density becomes lower, the ion motion retreats, and the focusing force becomes roughly linear and longitudinally independent again, see Fig. 8(e). Therefore, the non-Gaussian distributions that led to a slight increase in the emittance of each slice is reversed. In addition, the different phase space ellipse rotations for different longitudinal slices gradually also reverse such that the slices overlap with each other again at the exit [Fig. 8(f)]. After the entire acceleration stage, the witness beam doubled its energy from 10 to 20 GeV corresponding to a loaded transformer ratio of | E z , acc / E z , dec | = 0.62, while the normalized energy spread ( σ γ / γ ¯) decreases to 1% at the exit, and the projected emittance only grows by 1%.

2. Single particle simulation using Gaussian phenomenological model

In Sec. V B 1, we explained that the projected emittance first increases and then decreases back to near its original value due to the independent evolution of each beam slice. However, the evolution of the slice emittance is more subtle. In this section, we study the evolution of a beam slice in a predetermined focusing force. To take into account the longitudinal plasma density variation, we modify our model for the ion collapse in Eq. (4) and assume the following for the ion density as a function of r and z:
n ion ( r , z ) / n 0 = n ( z ) / n 0 + A ( z ) exp ( r 2 / 2 σ ion 2 ) ,
where A ( z ) = A 0 n ( z ) n ( 0 ) n 0 n ( 0 ), so that the nonlinear part of the ion density is exactly 0 at the entrance, and has a peak density of A0 at the plasma density plateau.
The corresponding focusing force is
F r ( r , z ) = [ r 2 n ( z ) n 0 A ( z ) σ ion 2 1 exp ( r 2 2 σ ion 2 ) r ] m e ω p 0 2 .
(17)
We choose A 0 = 100 , σ ion = 0.005 k p 0 1, and n ( z ) / n 0 to be same as in Fig. 6(a). This predetermined focusing force is used in Eq. (6) to solve for the phase space trajectories of a group of particles. The beam slice still has 10 GeV energy with ϵ n = 0.0594 k p 0 1 ( 1 μ m ), no initial energy spread, and is initially linearly matched at the plasma entrance. The beam has no longitudinal acceleration in the simulation. We use 105 particles in the beam.

In Fig. 9(a), we show the evolution of the emittance and spot size of the beam slice, which are similar to the results from the self-consistent QPAD simulation. The emittance grows in the upramp, remains relatively constant in the plateau, and decreases to near the original value at the plasma exit. To understand this emittance evolution, we need to analyze the beam's dynamics in more detail.

FIG. 9.

(a) Emittance and spot size evolution using a Gaussian phenomenological model. (b) Normalized phase space at the density plateau.

FIG. 9.

(a) Emittance and spot size evolution using a Gaussian phenomenological model. (b) Normalized phase space at the density plateau.

Close modal

At the entrance of the plasma upramp, the beam's distributions in x and px are both Gaussian, which remain matched distributions in a linear focusing force. As the beam propagates in the density upramp, the nonlinear part of the focusing force adiabatically increases. The beam will keep adjusting its profile to match to the local nonlinear focusing force, resulting in non-Gaussian profiles in both x and px [Figs. 10(c) and 10(d)], and a corresponding non-elliptical phase space [Fig. 10(b)], which looks like a rhombus. Particles at the “tips” make a larger contribution to the “spread” in x and px, leading to a larger value of the RMS emittance as compared to a 90% or 95% value. In other words, this RMS measurement of emittance is no longer an accurate reflection of the phase space area when the distribution is not Gaussian. This “emittance increase” is because the beam profile deviates from a Gaussian profile, or the phase space deviates from an ellipse. It is purely a geometrical effect and is reversible. This is different from the emittance growth caused by phase mixing within a slice [Fig. 4(b)], which is difficult to reverse.

FIG. 10.

(a)–(f) The distribution of beam particles' x and px for different z propagation distances.

FIG. 10.

(a)–(f) The distribution of beam particles' x and px for different z propagation distances.

Close modal

As the beam propagates in the density downramp, the nonlinear part of the focusing force adiabatically decreases, and becomes zero at the exit. The evolution of the beam's profile in the downramp is essentially reversed from that in the upramp. Eventually at the exit, the focusing force is back to linear, and the beam's distribution is back to Gaussian [Figs. 10(e) and 10(f)], so the RMS emittance remains unchanged from its original value.

In this section, we consider the utility of quasi-adiabatic matching sections for LC design parameters in the presence of ion motion.43 

For even the most stringent linear collider (LC) relevant witness beam parameters, quasi-adiabatic density ramps are still beneficial for preserving emittance. To demonstrate this, we undertake a QPAD simulation using a preformed hydrogen plasma with the same density n 0 = 10 17 cm 3 length for the plateau region (uniform acceleration section) which is kept at 4 × 10 4 k p 0 1 (67.3 cm) as before. The drive beam has a higher energy, 25 GeV, so the plateau could have been made larger; however, we keep it the same as for the FACET II case for ease of comparison. The drive beam has 3.0 × 10 10 electrons (4.8 nC), a normalized emittance ϵ n = 1 mm, and a tri-Gaussian density profile with σ r = 10 μ m ( 0.6 k p 0 1 ), σ z = 30 μ m ( 1.8 k p 0 1 ). It is thus matched to the uniform acceleration section. The drive beam's peak density n b 0 = 5.9 n 0 is high enough to produce a fully blown-out wake but not trigger any ion motion ( Φ b = 0.18). The drive beam is set to be non-evolving (it is represented by a specified current profile) in the simulation to isolate the physics. The witness beam also has an initial energy of 25 GeV, but with a normalized emittance, ϵ n = 100 nm, and a trapezoidal longitudinal current profile ranging from I b , head = 25.26 kA at the head to I b , tail = 6.42 kA at the tail. Its length is 1.8 k p 0 1 ( 30.3 μ m), so the total charge is 1.0 × 10 10 electrons (1.6 nC).

The head of the witness beam is located at a distance of 6.25 k p 0 1 ( 105.2 μ m) behind the center of the drive beam such that the blowout wake is optimally loaded,36,37 as shown in Fig. 11. The nearly constant accelerating field along the witness beam is around E z , acc = 1.2 m c ω p 0 / e for all beam slices, giving a transformer ratio of | E z , acc / E z , dec | = 1.6. The witness beam is initialized with no energy spread to isolate physics.

FIG. 11.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301, (2023);43 licensed under a Creative Commons Attribution (CC BY) license. A snapshot of the wake and beams in the density plateau. The color bar corresponds to the charge density for the plasma electrons, drive beam, and witness beam. The witness beam has a trapezoidal longitudinal current profile (black) in order to flatten the accelerating field Ez (on-axis lineout in red). The flattened accelerating field for the witness beam is E z , acc = 1.2 m e c ω p 0 / e, and the peak decelerating field for the drive beam is E z , dec = 0.75 m e c ω p 0 / e, giving a transformer ratio | E z , acc / E z , dec | = 1.6.

FIG. 11.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301, (2023);43 licensed under a Creative Commons Attribution (CC BY) license. A snapshot of the wake and beams in the density plateau. The color bar corresponds to the charge density for the plasma electrons, drive beam, and witness beam. The witness beam has a trapezoidal longitudinal current profile (black) in order to flatten the accelerating field Ez (on-axis lineout in red). The flattened accelerating field for the witness beam is E z , acc = 1.2 m e c ω p 0 / e, and the peak decelerating field for the drive beam is E z , dec = 0.75 m e c ω p 0 / e, giving a transformer ratio | E z , acc / E z , dec | = 1.6.

Close modal

To smoothly match the witness beam into the uniform acceleration section (plateau), we use a L = 3 × 10 4 k p 0 1 (50.5 cm) long quasi-adiabatic plasma density upramp with an analytical expression given by Eq. (16) before the density plateau and a symmetric downramp after the density plateau. The entire longitudinal plasma density profile is shown by the black solid curve in Fig. 12(a). We choose a small density at the entrance of the upramp n i = n ( 0 ) = 10 4 n 0, such that the corresponding α m i 1. We note that this does not strictly satisfy the adiabatic condition, which requires α m i 1. Thus, we refer to this as a quasi-adiabatic matching section. Strict compliance with α m i 1 will lead to a much longer ramp, which leads to a lower effective acceleration gradient. The trade-off between the adiabaticity of the ramp, the degree of initial ion motion ( Φ b), and the length of the ramp will be discussed later in Sec. VII. We match the witness beam to the plasma entrance with σ = σ m i and α = α m i, where σmi is Eq. (9) evaluated using the density at the entrance ni.

FIG. 12.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a) The resulting adiabatic plasma density profile (solid black) based on a linear ramp for αm. The profile assumes the witness beam energy remains constant, so the downramp is symmetric to the upramp. The black dashed downramp takes into account the energy doubling of the witness beam, thus it is 2 times as long as the solid black downramp. (b) Evolution of the projected emittance (solid black) and selected slice emittances. The witness beam's head and tail are located at ξ = 0 , 1.8 k p 0 1, respectively. The red dashed curve shows the evolution of the projected emittance when three azimuthal modes, m = 0, 1, 2, are kept in a QPAD simulation. The black dashed curve shows the evolution of the projected emittance for the case where the black dashed downramp in (a) is used.

FIG. 12.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a) The resulting adiabatic plasma density profile (solid black) based on a linear ramp for αm. The profile assumes the witness beam energy remains constant, so the downramp is symmetric to the upramp. The black dashed downramp takes into account the energy doubling of the witness beam, thus it is 2 times as long as the solid black downramp. (b) Evolution of the projected emittance (solid black) and selected slice emittances. The witness beam's head and tail are located at ξ = 0 , 1.8 k p 0 1, respectively. The red dashed curve shows the evolution of the projected emittance when three azimuthal modes, m = 0, 1, 2, are kept in a QPAD simulation. The black dashed curve shows the evolution of the projected emittance for the case where the black dashed downramp in (a) is used.

Close modal

In the QPAD simulation, we use a ( r , ξ ) simulation box with dimensions 12 k p 0 1 × 15 k p 0 1 (see Fig. 11), and cell sizes 5 × 10 4 k p 0 1 × 10 2 k p 0 1 ( 8.42 nm × 1.54 × 10 2 μ m). The witness beam consists of 107 numerical particles. The beam particles are pushed every 3D time step Δ s = 10 c / ω p 0. Both plasma electrons and ions are initialized with four locations for particles per r , ξ cell with 16 particles distributed in the azimuthal direction at each cell location when initializing plasma particles. We only kept the m = 0 mode, so we essentially have 64 plasma particles per r , ξ cell.

In Fig. 12(b), we show the evolution of the projected emittance (black solid) and the slice emittance at the head (blue), middle (orange), and tail (green) of the witness beam. The thickness of the slices is chosen to be Δ ξ = 0.1 k p 0 1. The projected emittance steadily grows in the density upramp, then grows very slowly in the density plateau,20 but eventually decreases in the density downramp, with only a 2 % net emittance growth at the exit of the plasma downramp. After the entire acceleration stage, the witness beam gained 25 GeV while maintaining a very small energy spread ( σ γ / γ ¯ 0.1 % at the exit). In Fig. 13, we show the focusing forces and the phase spaces at the head, middle, and the tail of the witness beam for three propagation distances that correspond to the entrance of upramp, middle of plateau, and exit of downramp.

We also carried out another simulation that included the m = 0, 1, 2 modes. The projected emittance for this simulation is shown in the red dashed curve in Fig. 12(b), which is similar to the m = 0 mode only result (solid black curve). This indicates that there is no hosing growth from noise during the entire stage.

The downramp is symmetric (mirror image) to the upramp for simplicity. However, the witness beam's energy doubled after being accelerated in the density plateau, so β m 0 becomes γ f / γ i = 2 times, where γi and γf denote the witness beam's energy at the entrance and the exit of the density plateau, respectively. From Eq. (16), we can see that in order to keep the same amount of adiabaticity, the length of the downramp should be 2 times that of the upramp. The modified profile is shown by the black dashed curve in Fig. 12(a). Using this downramp instead, the evolution of the projected emittance in the downramp is shown in Fig. 12(b). We can see that the projected emittance growth at the exit of the downramp is even smaller. In general, a longer ramp is more adiabatic and, therefore, better at preserving the beam emittance.

It is worth noting that the total length of the upramp, plateau, and downramp is 2.5 times the plateau along. Thus, the effective gradient of the full stage is 15 GeV/m and not 37 GeV/m. We also note that in the above simulations, the witness beam doubled its energy by gaining 25 GeV. Since the transformer ratio is 1.6 (Fig. 11), which is greater than one, we could have made the plasma density plateau longer, so that the witness beam can gain even more energy before the drive beam depletes its energy.

A ramp which is effectively quasi-adiabatic ( α m < 1 throughout) is extremely challenging to construct, therefore, we consider more general ramp profiles. Specifically, we use the following fifth order polynomial:
n ( z ) = n 0 [ 6 ( z L ) 5 15 ( z L ) 4 + 10 ( z L ) 3 ] ( 0 z L ) .
(18)
FIG. 13.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a), (c), and (e) Focusing force in x ̂ at different longitudinal ξ positions for different z propagation distances. (b), (d), and (f) Phase space ellipses corresponding to different longitudinal slices. The colors of each particle correspond to the colors of ξ for the focusing force. The distribution for px and x are shown for the ξ = 1.8 slice (green curve), and the Gaussian fits that match the standard deviations of the distributions are shown as red dashed lines.

FIG. 13.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a), (c), and (e) Focusing force in x ̂ at different longitudinal ξ positions for different z propagation distances. (b), (d), and (f) Phase space ellipses corresponding to different longitudinal slices. The colors of each particle correspond to the colors of ξ for the focusing force. The distribution for px and x are shown for the ξ = 1.8 slice (green curve), and the Gaussian fits that match the standard deviations of the distributions are shown as red dashed lines.

Close modal

This profile has the following properties: n ( 0 ) = n ( 0 ) = 0, n ( L ) = n 0 , n ( L ) = 0 , n ( L / 2 ) = n 0 / 2 , and n ( z ) is maximum at z = L / 2. We choose the length of the ramp to be the same as before, L = 3 × 10 4 k p 0 1. The density upramp is shown as the blue curve in Fig. 14. The density downramp is still chosen to be symmetric to that of the upramp. Unlike the fully adiabatic ramp given by Eq. (16), the density of the ramp given by Eq. (18) vanishes at the entrance. As a result, both the matched spot size [defined in Eq. (9)] and the adiabaticity parameter αm diverge toward infinity as the density approaches zero.

FIG. 14.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. The density profile for the fifth order polynomial upramp (blue) and the evolution of the C-S parameters. The solid green and orange curves are the β and α of the witness beam from the QPAD simulation, and the dashed curves are from numerical backward-propagation. The red curve shows the projected emittance evolution in the entire plasma density (upramp, plateau, and downramp) profile shown as the blue curve in Fig. 15(a).

FIG. 14.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301 (2023);43 licensed under a Creative Commons Attribution (CC BY) license. The density profile for the fifth order polynomial upramp (blue) and the evolution of the C-S parameters. The solid green and orange curves are the β and α of the witness beam from the QPAD simulation, and the dashed curves are from numerical backward-propagation. The red curve shows the projected emittance evolution in the entire plasma density (upramp, plateau, and downramp) profile shown as the blue curve in Fig. 15(a).

Close modal
Despite the divergence of the adiabaticity at the beginning of the ramp, we find it is still possible to roughly match the beam. In the absence of ion motion and any longitudinal acceleration, the C-S parameters evolve according to44 
1 2 β β 1 4 β 2 + β 2 k β 2 = 1 , α = 1 2 β .
(19)
Starting from the matched C-S parameters at the density plateau: β m = 2 γ k p 0 1, α m = 0 as the initial conditions for β and α, we numerically integrate Eq. (19) backwards along the upramp22 (see dashed curves in Fig. 14) to obtain estimates for β i , α i at the plasma entrance. It is worth noting that for an adiabatic ramp (where the density does not vanish at the entrance), the C-S parameters calculated using this numerical backward-propagation agree with the matched C-S parameters defined in Eq. (13). We then initialize the witness beam's C-S parameters using β i , α i in QPAD. All other drive and witness beam parameters are the same as that used in Sec. VI A, so the initial emittance of the witness beam is still 100 nm.

The solid green and orange curves in Fig. 14 show the actual evolution of the C-S parameters. Initially they follow the backward-propagated β , α because the density is low near the entrance such that the matched beam spot size is large. Therefore, ion motion effects are small and Eq. (19) is valid. Ion motion starts to play a role when the solid curves and the dashed curves deviate from each other where the beam has already entered the adiabatic region. The solid red curve in Fig. 14 shows the emittance evolution from a QPAD simulation for the entire plasma density profile shown (upramp, plateau, and downramp) in the blue curve in Fig. 15(a). As before, the emittance evolution reverses in the downramp, and the overall emittance growth is around 13%. It is likely that this could be improved by experimenting with the C-S parameters at the entrance.

FIG. 15.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301, (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a) Various plasma density profiles used for a 100 GeV witness beam. The downramps are symmetric to the corresponding upramps. The upramp in blue is the same as the upramp in Fig. 14. The density profile in black is the same as that in Fig. 12(a). The αmi is now 2 rather than 1 because the initial energy of the witness beam becomes four times higher (100 GeV compared to 25 GeV). The orange profile has an even lower density at the entrance to such that the witness beam has the same initial σm as for the case in Fig. 12, but at the expense of a worse adiabatic condition ( α m i = 4). (b) Projected emittance evolution for a 100 GeV witness beam when matching to the corresponding ramps in (a).

FIG. 15.

Reproduced with permission from Zhao et al., Phys. Rev. Accel. Beams 26, 121301, (2023);43 licensed under a Creative Commons Attribution (CC BY) license. (a) Various plasma density profiles used for a 100 GeV witness beam. The downramps are symmetric to the corresponding upramps. The upramp in blue is the same as the upramp in Fig. 14. The density profile in black is the same as that in Fig. 12(a). The αmi is now 2 rather than 1 because the initial energy of the witness beam becomes four times higher (100 GeV compared to 25 GeV). The orange profile has an even lower density at the entrance to such that the witness beam has the same initial σm as for the case in Fig. 12, but at the expense of a worse adiabatic condition ( α m i = 4). (b) Projected emittance evolution for a 100 GeV witness beam when matching to the corresponding ramps in (a).

Close modal
In order to preserve the emittance well, there are two important factors. First, at the plasma entrance the matched witness beam should ideally trigger negligible ion motion, which requires the density at the entrance ni to be low. Second, the plasma density ramp should ideally be adiabatic (a small αmi). However, for a given length of the ramp, improving one aspect will deteriorate the other. This trade-off can be seen by substituting z = 0 in Eq. (16) as follows:
n i n 0 = 1 ( 1 + α m i L β m 0 ) 2 .
(20)
From Eq. (20) we can see that for a given ramp length L, ni and αmi change in different directions. If we want to have less initial ion motion by decreasing ni, we need to sacrifice the adiabaticity of the plasma ramp, and vice versa. Another observation is that for a higher energy beam (e.g., in later acceleration stages), β m 0 is larger. The constraint from Eq. (20) tells us we have to either increase αmi (so the ramp becomes less adiabatic), or increase ni (so initially more ion motion is triggered). This makes the emittance preservation for progressively higher energy beams more challenging.

A TeV class PBA-LC will require limiting emittance growth over 20+ stages during which the witness beam energy will be higher for the later stages. As can be seen from Eq. (3) the matched spot size scales as γ 1 / 4. The ion parameter scales as n b / n 0. For a bunch with a given charge and length that scales as k p 0 1 it is straightforward to show that n b / n 0 γ / ϵ n. Therefore, ion motion and emittance growth will be more severe in the later stages. To examine this scaling, we carried out simulations with a 100 GeV witness beam. The longitudinal current profile of the witness beam is the same as that used in Sec. VI A such that the beam loading is still optimized at the density plateau. We still used preformed hydrogen plasma. We used the three different plasma density ramps shown in Fig. 15(a). The blue upramp is for the same profile used in Fig. 14. We match the witness beam through the numerical backward-propagation of Eq. (19). The black profile is the same as the solid black curve in Fig. 12, with the expression of the upramp given by Eq. (16). However, in this case α m i 2 since the witness beam's energy is 100 GeV. Furthermore, since the beam's energy is higher, σmi is smaller, leading to more initial ion motion. The orange upramp is also described by Eq. (16), but with n i / n 0 = 2.5 × 10 5, ensuring the same initial matched spot size, σmi, as the 25 GeV case. However, the use of a lower density leads to a larger α m i 4, such that the adiabatic condition is severely violated. For the black and orange profiles, we match the witness beam directly to the entrance rather than back propagate Eq. (19). Figure 15(b) shows the projected emittance evolution corresponding to the profiles in Fig. 15(a). We can see that the emittance growth of a 100 GeV/100 nm witness beam can be limited within 20% if we appropriately choose the plasma density ramp and match the beam.

In this paper, we examined emittance growth in a PBA-LC that is induced by the ion collapse caused by the intense space charge forces of matched witness beams. We used a combination of fully self-consistent quasi-static PIC simulations (QPAD) together with simulations of non-interacting particles (single particles) in prescribed forces.

Self-consistent QPAD simulations showed that when a witness beam is directly matched into a uniform plasma of density 10 17 cm 3, the resulting ion collapse leads to emittance growth. We considered two cases. In the first, which we refer to as FACET II-like parameters, the witness beam had an energy of 10 GeV and a normalized emittance of ϵ n 1 μ m (10 GeV/1 μm). For this case, the ion motion parameter was Φ b = 2.2. This led to 20 % emittance growth. In the second case, which we refer to as LC like, the beam had an energy of 25 GeV and a normalized emittance of ϵ n 100 nm (25 GeV/100 nm). This led to 80 % emittance growth. We also showed in QPAD simulations that the emittance growth within a slice can be mitigated by focusing the beam on a smaller spot size than the linearly matched spot size, such that it is closer to a nonlinear equilibrium state after ion motion is triggered. This agrees with the results found from single particle simulations, where we showed that nonlinearly matching the beam's spot size to the estimated steady state spot size of the Gaussian ion model can mitigate the emittance growth significantly.

As it is quite challenging to focus a beam on the matched spot size using conventional optics within meter scale distances and to create a beam with slice dependent spot size, we proposed to use quasi-adiabatic plasma density ramps as matching sections at both the entrance and exit, where ion collapse is tolerable at the entrance of the upramp and exit of the downramp but significant within the plateau. The matching sections significantly reduce the emittance growth. We showed that for both FACET II- or LC-like witness beam parameters, the quasi-adiabatic matching sections nearly preserve the witness beam emittance from start to end, even though there is a significant amount of ion motion triggered in the uniform acceleration stage. It is found that the slice emittance grows rapidly with the upramp and then grows slowly within the plateau. This growth is due to the phase space distribution of a slice becoming non-Gaussian. In addition, each slice has a different nonlinear k β and, thus, the projected emittance grows due to phase mixing between slices. Both of these effects are mostly reversible in the downramp leading to a final projected emittance growth of 0.25 % for the 10 GeV/1 μm case and 2 % for the 25 GeV/100 nm case.

It was also shown that the projected emittance growth for a (25 GeV/100 nm) beam is limited to 13% (it was 22% after the upramp) for general upramps and downramps that are not adiabatic at the low density regions. It was found that the best results are obtained by choosing Courant–Snyder parameters for the beam at the entrance based on integrating Eq. (19) from the plateau to the entrance.

Last, we considered LC parameters for a higher witness beam energy as would be the case in the later stages of a PBA-LC. We found that for a 100 GeV/100 nm beam, the projected emittance growth was limited to 25 % for a general ramp. Although it is currently not feasible to simulate the 500 GeV stages, based on the scaling of the ion motion parameter Φ b γ 1 / 4 / ϵ n we can infer that the relative effect of the emittance growth from ion motion for a 400 GeV/400 nm witness beam would scale as a 25 GeV/100 nm beam; and the emittance growth for a 400 GeV/200 nm beam would scale as a 100 GeV/100 nm beam. Thus, it seems feasible that for the electron arm of a PBA-LC (based on an electron beam or laser driver), it is possible to preserve the emittance to a few 100 nm through the use of quasi-adiabatic density upramp and downramp matching sections within each stage. While this work provides a possible path for emittance preservation for the electron arm of an LC, there is still much work to address, including misalignment between the drive beam and the witness beam,45–47 asymmetric beams,16,48 a more detailed assessment of the cumulative emittance growth through multi-stages,49 and maintaining spin polarization.1 

This work was supported by the U.S. Department of Energy (DOE) High Energy Physics (HEP) under Grant No. DE-SC0010064 and a SciDAC FNAL subcontract 644405, the National Science Foundation under Grant No. 2108970, and the National Natural Science Foundation of China (NSFC) Grant No. 12075030. The simulations were carried out on Cori at NERSC.

The authors have no conflicts to disclose.

Yujian Zhao: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Visualization (lead); Writing—original draft (lead). Lance Hildebrand: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing—review & editing (supporting). Weiming An: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing—review & editing (equal). Xinlu Xu: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Writing—review & editing (equal). Fei Li: Investigation (equal); Software (lead); Writing—review & editing (supporting). Thamine N. Dalichaouch: Investigation (equal); Software (equal); Writing—review & editing (supporting). Qianqian Su: Investigation (equal); Writing—review & editing (supporting). Chan Joshi: Conceptualization (equal); Funding acquisition (lead); Resources (equal); Writing—review & editing (equal). Warren B. Mori: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (equal); Project administration (lead); Resources (lead); Supervision (lead); Writing—review & editing (lead).

Raw data were generated at NERSC. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.

In this Appendix, we show how to design the plasma density upramp such that the αm linearly decreases from αmi to zero. Assume the plasma density upramp is in the range of 0 z L, where L is the length of the upramp, then we have
α m ( z ) = α m i z L L .
From α m ( z ) = 1 2 d β m ( z ) d z, we can get
β m ( z ) = β m 0 + α m i ( z L ) 2 L .
We also assume that the beam's energy is a constant when propagating in the density ramp (this assumption can be verified a posteriori). From β m ( z ) = 2 γ ¯ c ω p ( z ), we get
n ( z ) n 0 = ω p ( z ) 2 ω p 0 2 = β m 0 2 β m ( z ) 2 = 1 ( 1 + α m i ( z L ) 2 β m 0 L ) 2 .
The density at the entrance, ni, is related to αmi and L by the following constraint:
n i n 0 = n ( 0 ) n 0 = 1 ( 1 + α m i β m 0 L ) 2 .
So, among ni, αmi, and L, we can choose two parameters, and the third one is determined by the constraint above. Note that if we want the degree of ion motion to be small at the entrance (small ni), and also an adiabatic plasma ramp (small αmi), then we have to use a longer upramp (large L). So, there is a trade-off among the (1) degree of ion motion at the entrance, (2) adiabaticity of the ramp, and (3) length of the ramp. This trade-off is something we need to consider when designing the ramp.

In this Appendix, we calculate the ion motion parameter Φ b when the beam is linearly matched to the plasma. This follows arguments in Refs. 14 and 15.

Assume the density profile of the witness beam is transversely bi-Gaussian and longitudinally a flat-top, with a bunch length Lb, then the peak density is given by
n b 0 = N 2 π σ r 2 L b .
(B1)
If the beam is linearly matched to the plasma, then we substitute σr into Eq. (9) and obtain
n b 0 = γ ¯ 2 N ω p 2 π L b ϵ n c .
(B2)
Substituting the above expression into Eq. (3) leads to
Φ b = k i L b = ( Z n b 0 / n 0 2 M / m k p 0 2 L b 2 ) 1 2 = ( Z 1 2 π 1 2 M / m γ ¯ 2 N k p 0 3 n 0 L b ϵ n n n 0 ) 1 2 ,
(B3)
where n0 is the background unperturbed electron density and Z n i 0 = n 0. If the beam is Gaussian along z with a bunch length σz, it follows that
n b 0 = N ( 2 π ) 3 2 σ r 2 σ z ,
(B4)
Φ b = k i ( 2 π σ z ) = ( Z n b 0 / n 0 2 M / m k p 0 2 ( 2 π σ z ) 2 ) 1 2 = ( Z 1 2 π 1 2 M / m γ ¯ 2 N k p 0 3 n 0 σ z ϵ n n n 0 ) 1 2 .
(B5)
In this appendix, we numerically find the approximated steady state spot size of the beam in a given nonlinear focusing force described by the 1D Gaussian phenomenological model. Details of the derivation will appear in a future publication by Hildebrand et al. We start from the definition of the beam spot size ( means the ensemble average),
σ x = x 2 ,
and then take the first and second derivatives with respect to z to get
σ x = x x σ x
and
σ x = x 2 x 2 x x 2 σ x 3 + x x σ x = ϵ x 2 σ x 3 + x F x γ σ x ,
where we have assumed the beam has no energy spread. At steady state, σ x = 0. If Fx is linear in x, F x = k β x, then we obtain Eq. (7). However, even if Fx is nonlinear, useful expressions can still be obtained. For example, consider the 1D phenomenological model for Fx including ion motion,
n ion ( x ) / n 0 = 1 + A 0 exp ( x 2 / 2 σ ion 2 ) , F x ( x ) = x 2 A 0 σ ion 2 1 exp ( x 2 2 σ ion 2 ) x .
We cannot determine x F x as a function of σ because we do not know the beam distribution function in x. We, therefore, assume the beam's steady state distribution can be approximated by a Gaussian distribution as follows:
f ( x ) = 1 2 π σ x exp ( x 2 2 σ x 2 ) .
Under this assumption,
x F x = + x F x f ( x ) d x = σ x 2 2 A 0 σ ion 2 ( 1 σ ion σ x 2 + σ ion 2 ) ,
(C1)
which leads to the following steady state or nonlinear matching condition:
0 = ϵ x 2 σ x 3 1 γ ( σ x 2 + A 0 σ ion 2 σ x ( 1 σ ion σ x 2 + σ ion 2 ) ) .

For the given ion collapse model parameters (A0 and σion) and the given beam parameters (ϵx and γ), we can find a numerical solution to σx for the above equation, which is a reasonable estimate of beam's steady state spot size.

1.
C.
Joshi
,
S.
Corde
, and
W. B.
Mori
, “
Perspectives on the generation of electron beams from plasma-based accelerators and their near and long term applications
,”
Phys. Plasmas
27
,
070602
(
2020
).
2.
J. B.
Rosenzweig
,
B.
Breizman
,
T.
Katsouleas
, and
J. J.
Su
, “
Acceleration and focusing of electrons in two-dimensional nonlinear plasma wake fields
,”
Phys. Rev. A
44
,
R6189
(
1991
).
3.
A.
Pukhov
and
J. M.
ter Vehn
, “
Laser wake field acceleration: The highly non-linear broken-wave regime
,”
Appl. Phys. B
74
,
355
(
2002
).
4.
W.
Lu
,
C.
Huang
,
M.
Zhou
,
W. B.
Mori
, and
T.
Katsouleas
, “
Nonlinear theory for relativistic plasma wakefields in the blowout regime
,”
Phys. Rev. Lett.
96
,
165002
(
2006
).
5.
W.
Lu
,
M.
Tzoufras
,
C.
Joshi
,
F. S.
Tsung
,
W. B.
Mori
,
J.
Vieira
,
R. A.
Fonseca
, and
L. O.
Silva
, “
Generating multi-GeV electron bunches using single stage laser wakefield acceleration in a 3D nonlinear regime
,”
Phys. Rev. Spec. Top.--Accel. Beams
10
,
061301
(
2007
).
6.
K.
Floettmann
, “
Adiabatic matching section for plasma accelerated beams
,”
Phys. Rev. Spec. Top.--Accel. Beams
17
,
054402
(
2014
).
7.
I.
Dornmair
,
K.
Floettmann
, and
A. R.
Maier
, “
Emittance conservation by tailored focusing profiles in a plasma accelerator
,”
Phys. Rev. Spec. Top.--Accel. Beams
18
,
041302
(
2015
).
8.
X. L.
Xu
,
J. F.
Hua
,
Y. P.
Wu
,
C. J.
Zhang
,
F.
Li
,
Y.
Wan
,
C.-H.
Pai
,
W.
Lu
,
W.
An
,
P.
Yu
,
M. J.
Hogan
,
C.
Joshi
, and
W. B.
Mori
, “
Physics of phase space matching for staging plasma and traditional accelerator components using longitudinally tailored plasma profiles
,”
Phys. Rev. Lett.
116
,
124801
(
2016
).
9.
R.
Ariniello
,
C.
Doss
,
K.
Hunt-Stone
,
J.
Cary
, and
M.
Litos
, “
Transverse beam dynamics in a plasma density ramp
,”
Phys. Rev. Accel. Beams
22
,
041304
(
2019
).
10.
M. D.
Litos
,
R.
Ariniello
,
C. E.
Doss
,
K.
Hunt-Stone
, and
J. R.
Cary
, “
Beam emittance preservation using gaussian density ramps in a beam-driven plasma wakefield accelerator
,”
Philos. Trans. R. Soc. A
377
,
20180181
(
2019
).
11.
X.
Li
,
A.
Chancé
, and
P. A. P.
Nghiem
, “
Preserving emittance by matching out and matching in plasma wakefield acceleration stage
,”
Phys. Rev. Accel. Beams
22
,
021304
(
2019
).
12.
Y.
Zhao
,
W.
An
,
X.
Xu
,
F.
Li
,
L.
Hildebrand
,
M.
Hogan
,
V.
Yakimenko
,
C.
Joshi
, and
W.
Mori
, “
Emittance preservation through density ramp matching sections in a plasma wakefield accelerator
,”
Phys. Rev. Accel. Beams
23
,
011302
(
2020
).
13.
A.
Seryi
,
M.
Hogan
,
S.
Pei
,
T.
Raubenheimer
,
P.
Tenenbaum
,
T.
Katsouleas
,
C.
Huang
,
C.
Joshi
,
W.
Mori
, and
P.
Muggli
, “
A concept of plasma wake field acceleration linear collider (PWFA-LC)
,” in
Particle Accelerator Conference (PAC 09)
(
SLAC National Accelerator Lab
,
Menlo Park, CA
,
2010
), p.
WE6PFP081
.
14.
S.
Lee
and
T.
Katsouleas
, “
Wakefield accelerators in the blowout regime with mobile ions
,”
AIP Conf. Proc.
472
,
524
(
1999
).
15.
J. B.
Rosenzweig
,
A. M.
Cook
,
A.
Scott
,
M. C.
Thompson
, and
R. B.
Yoder
, “
Effects of ion motion in intense beam-driven plasma wakefield accelerators
,”
Phys. Rev. Lett.
95
,
195002
(
2005
).
16.
W.
An
,
W.
Lu
,
C.
Huang
,
X.
Xu
,
M. J.
Hogan
,
C.
Joshi
, and
W. B.
Mori
, “
Ion motion induced emittance growth of matched electron beams in plasma wakefields
,”
Phys. Rev. Lett.
118
,
244801
(
2017
).
17.
C. A.
Lindstrøm
, “
Staging of plasma-wakefield accelerators
,”
Phys. Rev. Accel. Beams
24
,
014801
(
2021
).
18.
R.
Gholizadeh
,
T. C.
Katsouleas
,
P.
Muggli
,
C. K.
Huang
, and
W. B.
Mori
, “
Preservation of beam emittance in the presence of ion motion in future high-energy plasma-wakefield-based colliders
,”
Phys. Review Lett.
104
(
15
),
155001
(
2010
).
19.
C.
Benedetti
,
C. B.
Schroeder
,
E.
Esarey
, and
W. P.
Leemans
, “
Emittance preservation in plasma-based accelerators with ion motion
,”
Phys. Rev. Accel. Beams
20
,
111301
(
2017
).
20.
C.
Benedetti
,
T. J.
Mehrling
,
C. B.
Schroeder
,
C. G. R.
Geddes
, and
E.
Esarey
, “
Adiabatic matching of particle bunches in a plasma-based accelerator in the presence of ion motion
,”
Phys. Plasmas
28
,
053102
(
2021
).
21.
F.
Li
,
W.
An
,
V. K.
Decyk
,
X.
Xu
,
M. J.
Hogan
, and
W. B.
Mori
, “
A quasi-static particle-in-cell algorithm based on an azimuthal Fourier decomposition for highly efficient simulations of plasma-based acceleration: QPAD
,”
Comput. Phys. Commun.
261
,
107784
(
2021
).
22.
C.
Joshi
,
E.
Adli
,
W.
An
,
C. E.
Clayton
,
S.
Corde
,
S.
Gessner
,
M. J.
Hogan
,
M.
Litos
,
W.
Lu
,
K. A.
Marsh
,
W. B.
Mori
,
N.
Vafaei-Najafabadi
,
B.
O'shea
,
X.
Xu
,
G.
White
, and
V.
Yakimenko
, “
Plasma wakefield acceleration experiments at FACET II
,”
Plasma Phys. Controlled Fusion
60
,
034001
(
2018
).
23.
Á.
Ferran Pousa
,
R.
Assmann
, and
A.
Martinez de la Ossa
, “
Wake-T: A fast particle tracking code for plasma-based accelerators
,”
J. Phys.
1350
,
012056
(
2019
).
24.
C.
Huang
,
V.
Decyk
,
C.
Ren
,
M.
Zhou
,
W.
Lu
,
W.
Mori
,
J.
Cooley
,
T.
Antonsen
, and
T.
Katsouleas
, “
QUICKPIC: A highly efficient particle-in-cell code for modeling wakefield acceleration in plasmas
,”
J. Comput. Phys.
217
,
658
(
2006
).
25.
W.
An
,
V. K.
Decyk
,
W. B.
Mori
, and
T. M.
Antonsen
, “
An improved iteration loop for the three dimensional quasi-static particle-in-cell algorithm: QuickPIC
,”
J. Comput. Phys.
250
,
165
(
2013
).
26.
A.
Lifschitz
,
X.
Davoine
,
E.
Lefebvre
,
J.
Faure
,
C.
Rechatin
, and
V.
Malka
, “
Particle-in-cell modelling of laser–plasma interaction using Fourier decomposition
,”
J. Comput. Phys.
228
,
1803
(
2009
).
27.
A.
Davidson
,
A.
Tableman
,
W.
An
,
F.
Tsung
,
W.
Lu
,
J.
Vieira
,
R.
Fonseca
,
L.
Silva
, and
W.
Mori
, “
Implementation of a hybrid particle code with a PIC description in rz and a gridless description in ϕ into OSIRIS
,”
J. Comput. Phys.
281
,
1063
(
2015
).
28.
R.
Lehe
,
M.
Kirchen
,
I. A.
Andriyash
,
B. B.
Godfrey
, and
J.-L.
Vay
, “
A spectral, quasi-cylindrical and dispersion-free particle-in-cell algorithm
,”
Comput. Phys. Commun.
203
,
66
(
2016
).
29.
P.
Mora
and
T. M.
Antonsen
, Jr.
, “
Kinetic modeling of intense, short laser pulses propagating in tenuous plasmas
,”
Phys. Plasmas
4
,
217
(
1997
).
30.
C.
Benedetti
,
C. B.
Schroeder
,
E.
Esarey
,
C. G. R.
Geddes
, and
W. P.
Leemans
, “
Efficient modeling of laser–plasma accelerators with INF&RNO
,”
AIP Conf. Proc.
1299
,
250
(
2010
).
31.
C.
Benedetti
,
C. B.
Schroeder
,
C. G. R.
Geddes
,
E.
Esarey
, and
W. P.
Leemans
, “
Efficient modeling of laser-plasma accelerator staging experiments using INF&RNO
,”
AIP Conf. Proc.
1812
,
050005
(
2017
).
32.
C.
Benedetti
,
C. B.
Schroeder
,
C. G. R.
Geddes
,
E.
Esarey
, and
W. P.
Leemans
, “
An accurate and efficient laser-envelope solver for the modeling of laser-plasma accelerators
,”
Plasma Phys. Controlled Fusion
60
,
014002
(
2017
).
33.
T.
Mehrling
,
C.
Benedetti
,
C.
Schroeder
, and
E.
Esarey
, “
A subgrid algorithm for the efficient modeling of plasma-based accelerators with ion motion using quasi-static particle-in-cell codes
,” in
IEEE Advanced Accelerator Concepts Workshop (AAC)
(
IEEE
,
2018
), pp.
1
5
.
34.
S.
Diederichs
,
C.
Benedetti
,
A.
Huebl
,
R.
Lehe
,
A.
Myers
,
A.
Sinn
,
J.-L.
Vay
,
W.
Zhang
, and
M.
Thévenet
, “
HiPACE++: A portable, 3D quasi-static particle-in-cell code
,”
Comput. Phys. Commun.
278
,
108421
(
2022
).
35.
T.
Mehrling
,
J.
Grebenyuk
,
F. S.
Tsung
,
K.
Floettmann
, and
J.
Osterhoff
, “
Transverse emittance growth in staged laser-wakefield acceleration
,”
Phys. Rev. Spec. Top--Accel. Beams
15
,
111303
(
2012
).
36.
M.
Tzoufras
,
W.
Lu
,
F. S.
Tsung
,
C.
Huang
,
W. B.
Mori
,
T.
Katsouleas
,
J.
Vieira
,
R. A.
Fonseca
, and
L. O.
Silva
, “
Beam loading in the nonlinear regime of plasma-based acceleration
,”
Phys. Rev. Lett.
101
,
145002
(
2008
).
37.
T. N.
Dalichaouch
,
X. L.
Xu
,
A.
Tableman
,
F.
Li
,
F. S.
Tsung
, and
W. B.
Mori
, “
A multi-sheath model for highly nonlinear plasma wakefields
,”
Phys. Plasmas
28
,
063103
(
2021
).
38.
R. L.
Williams
and
T.
Katsouleas
, “
Numerical studies on the ramped density plasma lens
,”
AIP Conf. Proc.
279
,
565
(
1992
).
39.
R.
Ariniello
,
C. E.
Doss
,
V.
Lee
,
C.
Hansel
,
J. R.
Cary
, and
M. D.
Litos
, “
Chromatic transverse dynamics in a nonlinear plasma accelerator
,”
Phys. Rev. Res.
4
,
043120
(
2022
).
40.
L.
Hildebrand
,
W.
An
,
X.
Xu
,
F.
Li
,
Y.
Zhao
,
M. J.
Hogan
,
V.
Yakimenko
,
S. S.
Nagaitsev
,
E.
Adli
,
C.
Joshi
, and
W.
Mori
, “
Mitigation techniques for witness beam hosing in plasma-based acceleration
,” in
IEEE Advanced Accelerator Concepts Workshop (AAC)
(
IEEE
,
2018
), p.
1
.
41.
T. J.
Mehrling
,
R. A.
Fonseca
,
A.
Martinez de la Ossa
, and
J.
Vieira
, “
Mitigation of the hose instability in plasma-wakefield accelerators
,”
Phys. Rev. Lett.
118
,
174801
(
2017
).
42.
P. G.
O'Shea
, “
Reversible and irreversible emittance growth
,”
Phys. Rev. E
57
,
1081
(
1998
).
43.
Y.
Zhao
,
L.
Hildebrand
,
W.
An
,
X.
Xu
,
F.
Li
,
T. N.
Dalichaouch
,
Q.
Su
,
C.
Joshi
, and
W. B.
Mori
, “
Emittance preservation in the presence of ion motion for low-to-high energy stages of a plasma based accelerator
,”
Phys. Rev. Accel. Beams
26
,
121301
(
2023
).
44.
S. Y.
Lee
,
Accelerator Physics
, 2nd ed. (
World Scientific Publishing Co. Pte. Ltd
,
2004
).
45.
M.
Thévenet
,
R.
Lehe
,
C. B.
Schroeder
,
C.
Benedetti
,
J.-L.
Vay
,
E.
Esarey
, and
W. P.
Leemans
, “
Emittance growth due to misalignment in multistage laser-plasma accelerators
,”
Phys. Rev. Accel. Beams
22
,
051302
(
2019
).
46.
C.
Lindstrøm
,
E.
Adli
,
J.
Pfingstner
,
E.
Marin
, and
D.
Schulte
, “
Transverse tolerances of a multi-stage plasma wakefield accelerator
,”
Report No. CERN-ACC-2016-212
,
2016
.
47.
L.
Hildebrand
, “
Beam realignment with emittance preservation in a plasma wakefield accelerator stage
,” Phys. Rev. Lett. (submitted) (
2024
).
48.
Y.
Zhao
, “
Emittance preservation in a plasma wakefield accelerator
,” Ph.D. thesis,
UCLA
,
Los Angeles (Main)
,
2023
.
49.
A. G. R.
Thomas
and
D.
Seipt
, “
Modeling chromatic emittance growth in staged plasma wakefield acceleration to 1 tev using nonlinear transfer matrices
,”
Phys. Rev. Accel. Beams
24
,
104602
(
2021
).