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.

## I. INTRODUCTION

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 $ \xi = c t \u2212 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 ramp^{6–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 \pi \sigma x \sigma y$, must be as large as $ \u223c 10 34 \u2009 cm \u2212 2 s \u2212 1$, where *N* is the number of particles in each colliding bunch, *f* is the repetition rate of collisions, and $ \sigma 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 $ \u223c 100 \u2009 nm = \u03f5 n , x \u03f5 n , y$, where $ \u03f5 n , x , \u03f5 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.

## II. ION MOTION PARAMETERS AND THE GAUSSIAN PHENOMENOLOGICAL MODEL

### A. Ion motion parameter

^{14,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 \u226b n 0$ (

*n*

_{0}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

*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 = \omega p 0 / c$ is the plasma wavenumber, where $ \omega p 0 = n 0 e 2 \u03f5 0 m$ is the plasma frequency,

*n*

_{0}is the plasma density,

*e*is the electron charge, and

*m*is the electron mass.

*ξ*= 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,

*L*is the length of the witness beam. Qualitatively, ion motion is important when $ \Phi b \u2273 1$. Typically, $ k p 0 L b \u223c 1$, so this occurs when $ Z n b 0 / n 0 2 M / m \u2273 1$.

_{b}It is worth noting that for a beam with a longitudinal Gaussian density profile, the length of the beam *L _{b}* is not well defined. However, we can use the effective bunch length $ L b = 2 \pi \sigma 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 $ \pi / 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 ( $ \pi 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.

### B. Gaussian phenomenological model

*r*(Ref. 16) (other functions can be used as well),

*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

*A*

_{0}and

*σ*. The corresponding focusing force is

_{ion}## III. METHODS

### A. Single particle simulation

*F*and

_{x}*F*are predetermined, and their form can be specified, based on the problems we are interested in. We generally ignore acceleration (changes in

_{y}*γ*). 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

*F*= 0). We also do not push the particles in

_{y}*ξ*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 fields

^{23}was also recently developed and is being widely utilized.

### B. Self-consistent simulation

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 \u2261 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 \u2261 s$ is fixed, and a “3D” part of the code where the beam is advanced using *s*. There is thus a 2D “time step,” $ \Delta \xi $, and a 3D time step, $ \Delta 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 $\varphi $. The fields, charge, and current defined on an *r*-*z* grid are expanded into azimuthal modes in the gridless $\varphi $ 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 $\varphi $ position. This so-called quasi-3D^{26–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 WAKE^{29} 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.

## IV. EMITTANCE EVOLUTION IN A UNIFORM PLASMA

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.

### A. Immobile ions

*n*

_{0}, 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*,

*σ*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.

_{m}*β*and

*α*to parameterize the beam, which are defined as follows:

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 *p _{x}* are normalized to $ \u03f5 n i \gamma i k \beta $ and $ \gamma i k \beta \u03f5 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}

### B. Mobile ions

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 \u2212 3$. The drive beam's density profile is tri-Gaussian, $ n b = n b 0 \u2009 exp ( \u2212 x 2 2 \sigma x 2 ) exp ( \u2212 y 2 2 \sigma y 2 ) exp ( \u2212 z 2 2 \sigma z 2 )$, with a symmetric transverse spot size $ \sigma x = \sigma y = 0.8 k p 0 \u2212 1 = 13.5 \u2009 \mu m$, longitudinal bunch length $ \sigma z = 2 k p 0 \u2212 1 = 23.8 \u2009 \mu 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 $ \Phi b = 0.12$ (calculated using an effective bunch length $ L b = 2 \pi \sigma 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 *I _{head}* = 27.2 kA at the head and

*I*= 17.6 kA at the tail, so the accelerating field

_{tail}*E*is nearly flattened across the witness beam ( $ E z \u223c \u2212 0.56 m c \omega p 0 / e$).

_{z}^{36,37}The length of the witness beam is $ L = 2 k p 0 \u2212 1 = 33.7 \u2009 \mu m$ (so the total charge of the witness beam is 2.5 nC), and its head is located $ 5 k p 0 \u2212 1$ behind the center of the drive beam. The witness beam has an initial energy of 10 GeV, and an initial normalized emittance of $ \u03f5 n i = 1 \u2009 \mu m = 0.06 k p 0 \u2212 1$ for $ n 0 = 10 17 \u2009 cm \u2212 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 , \xi )$ simulation box ( $ \xi = c t \u2212 z$) with 7340 cells in *r* and 2640 cells in *ξ*. The resolution is $ \Delta r = 1.36 \xd7 10 \u2212 3 k p 0 \u2212 1$ and $ \Delta \xi = 5 \xd7 10 \u2212 3 k p 0 \u2212 1$. The witness beam consists of ×10^{6} numerical particles, which are pushed every 3D time step $ \Delta s = 10 c / \omega p 0$. Both plasma electrons and ions are initialized at four locations uniformly spaced per $ r , \xi $ 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)] ( $ \sigma m = 0.025 k p 0 \u2212 1 = 0.41 \u2009 \mu m$, $ n b , head / n 0 \u2248 5.4 \xd7 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 ( $ \Phi b \u223c 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.

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 $ \u03f5 RMS \u2261 \u03f5 n , x \u03f5 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 \u2009 cm \u2212 3$, the tri-Gaussian drive beam has $ 3 \xd7 10 10$ electrons (4.8 nC), normalized emittance 1 $ mm$, and bunch length $ \sigma z = 30 \u2009 \u2009 \mu m \u2009 ( 1.8 k p 0 \u2212 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: $ \sigma r = 10 \u2009 \mu m \u2009 ( 0.6 k p \u2212 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 $ \Phi b = 0.18$ (calculated with an effective bunch length $ L b = 2 \pi \sigma z$). The witness beam has 10^{10} electrons (1.6 nC), an initial energy of 25 GeV, an initial energy spread of 2%, and an initial projected emittance of $ \u03f5 n i = 0.1 \u2009 \mu m \u2009 ( 0.006 k p 0 \u2212 1 )$. It has a tri-Gaussian density profile with bunch length $ \sigma z = 10 \u2009 \mu m \u2009 ( 0.6 k p 0 \u2212 1 )$. The ion parameter for the witness bunch is $ \Phi b \u2248 6$. The distance between the centers of the two beams is $ 115 \u2009 \mu m \u2009 ( 6.8 k p 0 \u2212 1 )$.

The $ ( r , \xi )$ simulation box has 32 768 × 1024 cells with a resolution of $ 3.6 \xd7 10 \u2212 4 k p 0 \u2212 1$ by $ 1.8 \xd7 10 \u2212 2 k p 0 \u2212 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 , \xi $) cell with 16 particles distributed in the azimuthal direction. The witness beam consists of 10^{7} numerical particles.

The beam particles are pushed every 3D time step $ \Delta s = 10 c / \omega 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 ( $ \sigma r = 0.1 \u2009 \mu m = 0.006 k p 0 \u2212 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.

### C. Phenomenological description of transverse phase space evolution from ion collapse using a single particle code

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 $ \sigma ion = 0.005 k p 0 \u2212 1$. The beam slice has 10 GeV energy with $ \u03f5 n = 0.0594 k p 0 \u2212 1 ( 1 \u2009 \mu m )$, and no initial energy spread. There is also no longitudinal acceleration included in the simulation. We use 10^{5} 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 *p _{y}* = 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.

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 $ \sigma ion = 4 \xd7 10 \u2212 4 k p 0 \u2212 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 $ \sigma ion = 10 \u2212 3 k p 0 \u2212 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)].

Note that it is shown in Ref. 16 (although in 1D) that a larger *A*_{0} 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.

## V. MATCHING USING ADIABATIC DENSITY RAMPS

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.

### A. Without ion motion

^{9,12,38}defined as

*z*is the propagation distance. Equation (14) is analogous to the relationship between

*α*and

*β*. An adiabaticity parameter can be defined as $ A = | \alpha 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

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}

*α*[defined in Eq. (14)] is positive. Therefore, we design the ramp with

_{m}*α*starting from $ \alpha 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),

_{m}*β*evaluated at the density plateau, and $ L = 180 \u2009 00 k p 0 \u2212 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.

_{m}In the QPAD simulation, we use a $ ( r , \xi )$ simulation box with dimensions $ 10 k p 0 \u2212 1 \xd7 12 k p 0 \u2212 1$, and cell size $ 2.5 \xd7 10 \u2212 3 k p 0 \u2212 1 \xd7 10 \u2212 2 k p 0 \u2212 1$. The witness beam consists of 10^{6} numerical particles. The beam particles are pushed every 3D time step $ \Delta s = 10 c / \omega p 0$. Both plasma electrons and ions are initialized at four locations uniformly spaced per $ r , \xi $ 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.

### B. With ion motion

#### 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 \u2264 \xi \u2264 2 k p 0 \u2212 1$, where *ξ* = 0 is at the beam head and $ \xi = 2 k p 0 \u2212 1$ is at the beam tail. The thickness of the slices is chosen to be $ \Delta \xi = 0.1 k p 0 \u2212 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 $ \u223c 1 %$ net growth at the exit of the plasma downramp. For the beam slice at the head of the beam (blue, $ \xi = 0.05 k p 0 \u2212 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, $ \xi = 1 k p 0 \u2212 1$) and tail (green, $ \xi = 1.95 k p 0 \u2212 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.

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 [ $ \Phi b \u223c 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 ( $ \xi = 0.05$), middle ( $ \xi = 1.0 k p 0 \u2212 1$), and tail ( $ \xi = 1.95 k p 0 \u2212 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

*p*values are normalized to their matched values at a given

_{x}*z*in the absence of ion motion.

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 $ \Phi b \u223c 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 ( $ \sigma \gamma / \gamma \xaf$) decreases to 1% at the exit, and the projected emittance only grows by 1%.

#### 2. Single particle simulation using Gaussian phenomenological model

*r*and

*z*:

*A*

_{0}at the plasma density plateau.

^{5}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.

At the entrance of the plasma upramp, the beam's distributions in *x* and *p _{x}* 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

*p*[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}*x*and

*p*, 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.

_{x}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.

## VI. LC PARAMETERS

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

### A. Quasi-adiabatic ramps

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 \u2212 3$ length for the plateau region (uniform acceleration section) which is kept at $ 4 \xd7 10 4 \u2009 k p 0 \u2212 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 \xd7 10 10$ electrons (4.8 nC), a normalized emittance $ \u03f5 n = 1 \u2009 \u2009 mm$, and a tri-Gaussian density profile with $ \sigma r = 10 \u2009 \mu m \u2009 ( 0.6 k p 0 \u2212 1 )$, $ \sigma z = 30 \u2009 \mu m \u2009 ( 1.8 k p 0 \u2212 1 )$. It is thus matched to the uniform acceleration section. The drive beam's peak density $ n b 0 = 5.9 \u2009 n 0$ is high enough to produce a fully blown-out wake but not trigger any ion motion ( $ \Phi 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, $ \u03f5 n = 100 \u2009 \u2009 nm$, and a trapezoidal longitudinal current profile ranging from $ I b , head = 25.26 \u2009 kA$ at the head to $ I b , tail = 6.42 \u2009 \u2009 kA$ at the tail. Its length is $ 1.8 \u2009 k p 0 \u2212 1$ ( $ 30.3 \u2009 \u2009 \mu m$), so the total charge is $ 1.0 \xd7 10 10$ electrons (1.6 nC).

The head of the witness beam is located at a distance of $ 6.25 \u2009 k p 0 \u2212 1$ ( $ 105.2 \u2009 \u2009 \mu 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 = \u2212 1.2 \u2009 m c \omega 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.

To smoothly match the witness beam into the uniform acceleration section (plateau), we use a $ L = 3 \xd7 10 4 \u2009 k p 0 \u2212 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 \u2212 4 n 0$, such that the corresponding $ \alpha m i \u2248 1$. We note that this does not strictly satisfy the adiabatic condition, which requires $ \alpha m i \u226a 1$. Thus, we refer to this as a quasi-adiabatic matching section. Strict compliance with $ \alpha m i \u226a 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 ( $ \Phi b$), and the length of the ramp will be discussed later in Sec. VII. We match the witness beam to the plasma entrance with $ \sigma = \sigma m i$ and $ \alpha = \alpha m i$, where *σ _{mi}* is Eq. (9) evaluated using the density at the entrance

*n*.

_{i}In the QPAD simulation, we use a $ ( r , \xi )$ simulation box with dimensions $ 12 k p 0 \u2212 1 \xd7 15 k p 0 \u2212 1$ (see Fig. 11), and cell sizes $ 5 \xd7 10 \u2212 4 k p 0 \u2212 1 \xd7 10 \u2212 2 k p 0 \u2212 1$ ( $ 8.42 nm \xd7 1.54 \xd7 10 \u2212 2 \mu m$). The witness beam consists of 10^{7} numerical particles. The beam particles are pushed every 3D time step $ \Delta s = 10 c / \omega p 0$. Both plasma electrons and ions are initialized with four locations for particles per $ r , \xi $ 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 , \xi $ 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 $ \Delta \xi = 0.1 k p 0 \u2212 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 $ \u223c 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 ( $ \sigma \gamma / \gamma \xaf \u223c 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 $ \beta m 0$ becomes $ \gamma f / \gamma i = 2$ times, where *γ _{i}* and

*γ*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.

_{f}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.

### B. Realistic ramps

This profile has the following properties: $ n ( 0 ) = n \u2032 ( 0 ) = 0$, $ n ( L ) = n 0 ,$ $ n \u2032 ( L ) = 0 , \u2009 n ( L / 2 ) = n 0 / 2 , \u2009$ and $n\u2032(z)$ is maximum at $ z = L / 2$. We choose the length of the ramp to be the same as before, $ L = 3 \xd7 10 4 k p 0 \u2212 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.

^{44}

*β*and

*α*, we numerically integrate Eq. (19) backwards along the upramp

^{22}(see dashed curves in Fig. 14) to obtain estimates for $ \beta i , \alpha 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 $ \beta i , \alpha 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 $ \beta , \alpha $ 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.

## VII. 100 GeV WITNESS BEAM SIMULATION

*n*to be low. Second, the plasma density ramp should ideally be adiabatic (a small

_{i}*α*). However, for a given length of the ramp, improving one aspect will deteriorate the other. This trade-off can be seen by substituting

_{mi}*z*= 0 in Eq. (16) as follows:

*L*,

*n*and

_{i}*α*change in different directions. If we want to have less initial ion motion by decreasing

_{mi}*n*, 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), $ \beta m 0$ is larger. The constraint from Eq. (20) tells us we have to either increase

_{i}*α*(so the ramp becomes less adiabatic), or increase

_{mi}*n*(so initially more ion motion is triggered). This makes the emittance preservation for progressively higher energy beams more challenging.

_{i}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 $ \gamma \u2212 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 \u2212 1$ it is straightforward to show that $ n b / n 0 \u223c \gamma / \u03f5 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 $ \alpha m i \u223c 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 \xd7 10 \u2212 5$, ensuring the same initial matched spot size,

*σ*, as the 25 GeV case. However, the use of a lower density leads to a larger $ \alpha m i \u223c 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.

_{mi}## VIII. CONCLUSIONS

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 \u2009 cm \u2212 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 $ \u03f5 n \u223c 1 \mu m$ (10 GeV/1 *μ*m). For this case, the ion motion parameter was $ \Phi b = 2.2$. This led to $ \u223c 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 $ \u03f5 n \u223c 100$ nm (25 GeV/100 nm). This led to $ \u223c 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 \beta $ 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 $ \u2248 0.25 %$ for the 10 GeV/1 *μ*m case and $ \u2248 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 $ \u2248 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 $ \Phi b \u223c \gamma 1 / 4 / \u03f5 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}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: ADIABATIC RAMP DESIGN

*α*linearly decreases from

_{m}*α*to zero. Assume the plasma density upramp is in the range of $ 0 \u2264 z \u2264 L$, where

_{mi}*L*is the length of the upramp, then we have

*a posteriori*). From $ \beta m ( z ) = 2 \gamma \xaf c \omega p ( z )$, we get

*n*, is related to

_{i}*α*and

_{mi}*L*by the following constraint:

*n*,

_{i}*α*, and

_{mi}*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

*n*), and also an adiabatic plasma ramp (small

_{i}*α*), then we have to use a longer upramp (large

_{mi}*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.

### APPENDIX B: ION MOTION PARAMETER FOR A LINEARLY MATCHED BEAM

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

*L*, then the peak density is given by

_{b}*σ*into Eq. (9) and obtain

_{r}*n*

_{0}is the background unperturbed electron density and $ Z n i 0 = n 0$. If the beam is Gaussian along

*z*with a bunch length

*σ*, it follows that

_{z}### APPENDIX C: NONLINEAR MATCHING THE SPOT SIZE IN A UNIFORM PLASMA

*et al.*We start from the definition of the beam spot size ( $ \u27e8 \u27e9$ means the ensemble average),

*z*to get

*F*is linear in

_{x}*x*, $ F x = \u2212 k \beta x$, then we obtain Eq. (7). However, even if

*F*is nonlinear, useful expressions can still be obtained. For example, consider the 1D phenomenological model for

_{x}*F*including ion motion,

_{x}*σ*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:

For the given ion collapse model parameters (*A*_{0} and *σ _{ion}*) and the given beam parameters (

*ϵ*and

_{x}*γ*), we can find a numerical solution to

*σ*for the above equation, which is a reasonable estimate of beam's steady state spot size.

_{x}## REFERENCES

*r*–

*z*and a gridless description in

*ϕ*into OSIRIS

*Accelerator Physics*