Chromatic Dynamics of an Electron Beam in a Plasma Based Accelerator

We present a theoretical investigation of the chromatic dynamics of the witness beam within a plasma based accelerator. We derive the single particle motion of an electron in an ion column within a nonlinear, blowout wake including adiabatic dampening and adiabatic variations in plasma density. Using this, we calculate the evolution of the beam moments and emittance for an electron beam. Our model can handle near arbitrary longitudinal phase space distributions. We include the effects of energy change in the beam, imperfect wake loading, initial transverse offsets of the beam, and mismatch between the beam and plasma. We use our model to derive analytic saturation lengths for the projected, longitudinal slice, and energy slice emittance under different beam loading conditions. Further, we show that the centroid oscillations and spot sizes vary between the slices and the variation depends strongly on the beam loading. Next, we show how a beam evolves in a full plasma source with density ramps and show that the integral of the plasma density along the ramp determines the impact on the beam. Finally, we derive several simple scaling laws that show how to design a plasma based injector to produce a target beam energy and energy spread.


I. INTRODUCTION
Accelerators, in the form of high energy colliders and light sources, have proven to be important tools for a diverse range of research fields. Unfortunately, the size and cost of these machines are prohibitive, especially at the energy frontier. Plasma based accelerators are promising, compact alternatives that have been shown to produce accelerating gradients two to three orders of magnitude larger than conventional radio frequency accelerators. Great progress has been made in demonstrating low energy spread beams and high efficiency acceleration [1][2][3]. Colliders and light sources, however, require beams with high brightness and thus place strict limits on the beam's transverse emittance. Current plasma based accelerators struggle to meet these strict requirements. The emittance of the accelerated beam, called the witness beam, typically grows considerably as the beam traverses an accelerating stage.
Multiple mechanisms contribute to the beam's emittance growth within the plasma stage. Recent work has primarily focused on emittance growth due to mismatch between the plasma and the beam [4][5][6][7][8][9][10][11] and beam instabilities [12][13][14][15][16][17][18][19][20][21]. Yet the emittance growth of a transversely offset beam-due to chromatic phase spread-has only been briefly considered: Refs. [22,23] derived simple expressions for the saturated emittance and initial growth rate, Ref. [24] considered emittance in the presence of a laser driver, and Ref. [25] worked out the saturated emittance for an injection mismatch in a conventional accelerator. In the case of mismatch, emittance growth due to energy gain in the presence of plasma ramps, or with the inclusion of wake loading, has not been considered. Typically, the longitudinal phase space is assumed to take on a simple form and previous approaches cannot handle arbitrary distributions. Further, the longitudinal slice emittance has only been briefly in-vestigated [24,26,27], while the energy slice emittance has yet to receive serious attention despite its importance in a transverse gradient undulator [28][29][30].
We derive the projected emittance, slice emittance (longitudinal and energy), and moment evolution of an electron beam travelling in a nonlinear, blowout wake including the effects of adiabatic plasma ramps, energy gain, beam loading, initial transverse offsets, and the initial longitudinal phase space distribution. We start by deriving the single particle motion of an electron, including energy change, in an ion column with adiabatically varying density. Next, we derive the evolution of the beam moments and the emittance in a way that allows straightforward evaluation of the slice and projected beam parameters. We separate the effect of the beam's initial longitudinal phase space into a single parameter that can be analytically evaluated in simple cases. For complicated phase space distributions, this parameter can be evaluated numerically while retaining the rest of the analytic formulation.
To show the broad applicability of our model, we present several examples. First, we consider a beam in a uniform plasma with continuous energy gain. We use two simple models for the longitudinal accelerating field to represent a beam that overloads the wake and one that does not sufficiently load the wake. In both cases, we derive analytic saturation lengths for the projected, longitudinal slice, and energy slice emittance. We show that beam loading causes particles to mix between energy slices leading to growth in the energy slice emittance. We then calculate the transverse offset of the different longitudinal slices and the spot size of the different energy slices. We show that the slice parameters depend on the beam loading -a result with experimental consequences. Second, we calculate the witness beam evolution through a full plasma source with density ramps, and show that the integral over the ramp density determines how much of an impact the ramp has on the beam. We show why it is the ramp shape, and not the length, that is important. Third, we calculate the beam evolution in a plasma based injector and derive simple, analytic expressions for designing an injector to produce a beam with a target energy and energy spread. These examples demonstrate the generality of our approach. It combines multiple effects in a straightforward analytic framework.

II. SINGLE PARTICLE MOTION
In a plasma based accelerator operating in the blowout regime, all of the plasma electrons are evacuated from the center of the wake leaving a column of ions behind. If the plasma has a transversely uniform density and ion motion is neglected, the ions produce a linear focusing force on the witness beam that is independent of ξ = ct− s, where s is the distance along the accelerator [31,32]. The equations of motion in the transverse directions for an electron in the witness beam are decoupled and given by [9,33] where a prime denotes a derivative with respect to s. The focusing strength is given by In general, K is a function of s through the local plasma density and the particle energy. The plasma frequency ω p is defined as ω 2 p (s) = n(s)e 2 /(m e ǫ 0 ), γ b (s) is the relativistic factor of the electron, c is the speed of light, n(s) is the plasma density, e is the elementary charge, m e is the mass of the electron, and ǫ 0 is the permittivity of free space. The betatron wavenumber is defined as k 2 β (s) = K(s).
An approximate solution to the equation of motion can be derived for a plasma density that varies adiabatically in s. The adiabatic condition is defined as [10]: where α m is one of the matched Courant-Snyder (CS) parameters defined for a single particle as: β m = 1/k β , α m = −β ′ m /2, and γ m = (1 + α 2 m )/β m . β m sets the natural length scale of the transverse evolution. In a uniform plasma with no energy gain, β m is constant and α m = 0. The matched CS parameters are functions of s due to the variation in plasma density and particle energy along the accelerator.
In the absence of energy gain, the transverse motion of a single particle in an adiabatically varying plasma density is given by [10,11] x =x 0 β m β m0 where the subscript 0 indicates the initial value of a variable at s = s 0 and x ′ = dx/ds = p x /p s , p x and p s are the longitudinal and transverse momentum of the particle, respectively. The betatron phase advance is defined as We are interested in the general case where the particles can gain and lose energy. Energy gain (loss) results in a reduction (increase) in the amplitude of the oscillations due to adiabatic dampening. We use the ansatz that the single particle motion has an additional s dependent amplitude A(s): where x c is the position of a single particle with constant energy given by Eq. (4). x c satisfies the differential equation x ′′ c + k 2 β x c = 0. Inserting the ansatz into Eq. (1) and requiring that A can take on any value gives a differential equation for A: The solution of which is A(s) = γ b0 /γ b . For Ax c to approximately satisfy the equation of motion, Eq. (1), the relative change in energy over one betatron period must be small The transport matrix M defines the motion of the particle based on its initial conditions Combining the adiabatic dampening term A(s) with x c gives the transport matrix for a particle in an adiabatic plasma with energy change: In a uniform plasma, β m = β m0 γ b /γ b0 . Inserting β m into Eq. (9) and assuming the energy gain over a single betatron period is small, α m ≈ 0, we recover the expression for single particle motion in a uniform plasma from [9,33]. Fig. 1 shows a comparison between Eq. (8) and a numerical solution to Eq. (1) for a single electron travelling through a plasma based accelerator. The plasma stage is designed to double the particle's energy and has identical adiabatic density ramps on each end. The entrance ramp reduces the incoming beta function by a factor of 10, and the beam undergoes approximately 12 betatron periods within the uniform plasma section. k βu is the betatron wavenumber at the start of the uniform plasma section, located at sk βu = 200. The parameters used here are similar to those found in current beam driven plasma wakefield accelerator experiments such as those at FACET-II [34,35] and FLASHForward [36,37]. The analytic solution shows excellent agreement with the numerical solution; it accurately captures both the adiabatic dampening and focusing in the ramp. In the next section, we use this solution to derive general expressions for the evolution of the beam's moments and emittance.

III. EVOLUTION OF THE BEAM MOMENTS AND PROJECTED EMITTANCE
The transverse beam quality can be quantified using the normalized emittance [38]: where the beam sizes are σ 2 x − p x 2 and the correlation is σ xpx = xp x − x p x . In the ultra-relativistic limit the normalized emittance simplifies to where σ 2 δ = ( γ 2 b − γ b 2 )/ γ b 2 and ǫ is the geometric emittance defined as with σ 2 In most cases, the first term in Eq. (11) is very small and the geometric emittance dominates, giving the more familiar formula ǫ n ≈ γ b ǫ.
Emittance growth results from the γ b dependence of k β . Different energy slices of the beam oscillate with different frequencies in the ion channel. Over time the different slices dephase and no longer overlap in transverse phase space as shown in Fig. 2. Even if the beam has no initial energy spread, this chromatic dephasing will occur if imperfect wake loading induces energy spread in the beam during acceleration. This same effect dampens out the centroid oscillations of an offset beam. To calculate the growth explicitly requires the evaluation of the moments of the beam distribution.
The beam moments can be evaluated under the assumption that there is no initial correlation between the longitudinal distribution and the transverse distribution. In this case the phase space distribution of the witness beam at s 0 can be written as where δ = (p − p 0 )/p 0 parameterizes the energy spread in the beam and p 0 is the momentum of the reference particle. Here, we are ignoring the vertical transverse dimension because it is decoupled and evolves independently in an analogous fashion. We let f ⊥ describe the witness beam without offset; i.e., the beam centroid is located at the origin in transverse phase space: Further, we assume that the distribution is normalized to 1.
If the witness beam is then given an initial offset of ∆x in x and ∆x ′ in x ′ , the beam distribution at s 0 can be FIG. 2. When the witness beam is mismatched or offset transversely with respect to the ion column, it will undergo emittance growth through chromatic dephasing. Different energy slices of the beam will rotate at different frequencies in transverse phase space, leading to growth in the projected emittance (the area enclosed by the white dashed line grows to the area enclosed by the black dashed line).
Single particle evolution is described by Eq. (8) and is a function of s, , where the ξ 0 and δ 0 dependence enters through the relativistic factor γ b (s, ξ 0 , δ 0 ). The functional form of γ b depends on the acceleration model chosen, we show several in the examples. Making the change of variables to u = x − ∆x and u ′ = x ′ − ∆x ′ , we can write the single particle evolution in terms of u and u ′ : x(s, u 0 + ∆x, u ′ 0 + ∆x ′ , ξ 0 , δ 0 ). The first beam moment is given by: The other beam moments, x 2 , x ′ , x ′2 , and xx ′ take the same form.
To make the integral tractable, we make the assumption that the energy spread is small: δ ≪ 1. This allows us to treat the matched CS parameters and the adiabatic dampening term, A = γ b0 /γ b , in Eq. (9) as constant with respect to ξ 0 and δ 0 . We use bars to denote values for the reference particle:γ b = γ b (s, ξ 0 = 0, δ 0 = 0). Using γ b =γ b in A, β m , α m , and γ m for all particles defines a matched set of CS parameters for the beam; this is in contrast to the approach taken in Ref. [10], where different matched CS parameters are defined for each energy slice. The only dependence on ξ 0 and δ 0 remaining in Eq. (9) is in the energy dependence of the betatron phase φ.
We have to integrate over u and u ′ in Eq. (13) first because x depends on the longitudinal coordinates. The u and u ′ integrals can be written in terms of moments of the witness beam's initial (non-offset) phase space distribution. The second central moments of the phase space distribution, σ x , σ x ′ , and σ xx ′ , are defined as: . These moments can be expressed in terms of the CS parameters: β = σ 2 x /ǫ, γ = σ 2 x ′ /ǫ, and α = −σ xx ′ /ǫ. Evaluating the u, u ′ integrals, and writing the beam moments in terms of the beam's initial CS parameters gives the following expressions for the moments of the offset beam: where we have kept terms of order α 2 m . The a, b, and d terms are all unit-less constants that depend on the beam's initial conditions: The first three constants have physical meaning; d 1 and d 2 are the normalized phase space coordinates for the motion of the beam centroid and, as we will show later, a is the relative emittance growth at saturation. The s dependence of the moments enters throughγ b directly (adiabatic dampening) and through the matched CS parameters (plasma focusing), which are functions of γ b ; however, the primary dependence on s is given by the C and S terms which capture both the chromatic phase spreading and the betatron oscillations. These terms are integrals over the longitudinal phase space: The two dimensional longitudinal phase space distribution can always be reduced to a one dimensional distribution of betatron phase advance f φ (φ). The integrals can then be written in the form: where we have only written the C 1 integral for brevity. The s dependence enters through the variation of the f φ distribution as the beam propagates through the plasma. At this point, the saturated emittance can be determined.
As the different energy components of the beam dephase, the width of the f φ distribution will grow. At saturation, the width is much larger than 2π and the C and S integrals tend towards zero. Setting C and S to zero in the moments and calculating the emittance gives the saturated emittance where terms higher than α m have been dropped. The pre-factor 1 + σ 2 δ can be set equal to one if the final energy spread is sufficiently small.
The detailed emittance evolution can be found by rewriting the C and S integrals. We start by writing the betatron phase as φ =φ + ∆φ(ξ 0 , δ 0 ). Further, we note that the four C and S integrals can be written as the real and imaginary parts of two complex integrals. After removing the e iφ term from the them, the integrals generically evaluate to a pair of complex numbers with real amplitudes H 1 , H 2 , and arguments ψ 1 , ψ 2 : The C and S integrals are then given by If one of three conditions is met or approximately met, the expression for the emittance simplifies significantly. The first condition is if In all the examples we consider later, either the first or third condition is satisfied. The emittance then evolves according to where H 2 terms describe emittance growth due to mismatch and H 1 terms describe emittance growth due to transverse offsets in position and angle. In addition, H 1 , when combined with the adiabatic dampening, describes the dampening of the beams centroid oscillations. The H(s) functions describe the amount of chromatic coherence in the beam, they get smaller as the beam propagates and dephases. ψ 1 is the betatron phase difference between the projected beam and the reference particle. In Fig. 3 we compare the results of numerical particle tracking with Eq. (23) for a uniform plasma without energy gain. The witness beam has an initial energy spread and is either mismatch to the plasma (red), offset transversely from the drive beam (blue), or both mismatched and offset (black). In this simple case, Eq. (5) gives the betatron phase advance as Here, φ is independent of ξ 0 because the accelerating field is constant (in this case zero) in ξ 0 . The betatron phase advance of the reference particle isφ =k β s, is the betatron wavenumber of the reference particle. Expanding to first order in δ 0 simplifies the expression to φ =φ −φδ 0 /2. Since φ depends linearly on δ 0 , the I 1 and I 2 integrals are Fourier transforms of the energy distribution with frequencies ω δ =φ/2 and 2ω δ =φ, respectively. If the distribution of energy is , and ψ 1 = ψ 2 = 0. Here,f δ is the Fourier transform of the energy distribution f δ . Eq. (23) then fully describes the evolution of the projected emittance.
In the case of a uniform energy spread between −∆δ/2 and ∆δ/2,f δ (φ) = sinc∆Φ andf δ (φ/2) = sinc(∆Φ/2), where ∆Φ is the range of φ spanned by the energy spread ∆Φ = ∆δφ. In the absence of transverse offset, this formula reduces to that given in Ref. [8]. The emittance growth of such a beam is shown in Fig. 3(a) for various initial mismatches and transverse offsets. The ana- lytic expression shows excellent agreement with the numerical particle tracking code. The small oscillations of the numerical result about the theoretical emittance from Eq. (23) are due to treating β m and γ b as constants and are of order σ 2 δ (for more details see Ref. [10]). There are two length scales for the emittance growth. The growth due to mismatch saturates when ∆Φ = π giving a saturation length of s 2 = π/(∆δk β ), while the emittance growth due to transverse offset requires twice the length to saturate: s 1 = 2π/(∆δk β ).
In the case of a Gaussian distribution of energy, Without any transverse offset, ∆x = ∆x ′ = 0, the solution reduces to that given in Ref. [9]. The theoretical solution is compared to numerical particle tracking in Fig. 3(b). As can be seen in the figure, the saturated emittance is the sum of the emittance growth due to an offset and mismatch.

IV. EVOLUTION OF THE SLICE MOMENTS AND SLICE EMITTANCE
Although the projected emittance is normally taken as the figure of merit, both the longitudinal slice emittance and the energy slice emittance are of practical interest. We start with the longitudinal slice emittance, which is of importance for both light sources and colliders. For an individual slice at ξ, the I 1 integral takes the form: where N (ξ) = dδ 0 f (ξ, δ 0 ) is a normalization factor and δ D is the Dirac delta function. I 2 has the same form with ∆φ replaced by 2∆φ. The expressions for the moments and emittance of each slice are given by Eqs. (14)- (18) and Eq. (23), but with H and ψ replaced by the ξ dependent expressions from Eq. (24). Further accuracy can be achieved by usingγ b and β m calculated for each slice rather than the beam as a whole. Similarly, the energy slice emittance is found by replacing the I 1 and I 2 integrals with where γ b is a function of ξ 0 and δ 0 , . Again, I 2 has the same form as I 1 with ∆φ replaced by 2∆φ. Growth in the energy slice emittance is due to particle exchange between energy slices as a consequence of imperfect loading of the wake.
Energy slice emittance is important for two reasons: First, projected emittance growth can be reversed down to the energy slice emittance using a suitable apochromatic beam line. An example of this technique is presented in Ref. [39]. Second, the energy slice emittance influences the results of some emittance measurements, such as those taken using the butterfly technique. Third, the performance of a transverse gradient undulator depends on the energy slice emittance [28][29][30].

V. EMITTANCE GROWTH OF A BEAM WITH IMPERFECT WAKE LOADING
In the previous sections, we develop general expressions for the projected and slice emittance growth. In this section, we work out analytic solutions for the projected and slice emittance for two situations where the witness beam does not perfectly load the wake.
Perfect loading requires the witness beam to have a trapezoidal current profile to produce a uniform accelerating field [40]. Existing accelerators and PWFA experiments tend to use beams that have Gaussian or other non-trapezoidal current profiles. Depending on the length and current of the witness beam, the accelerating field can be approximated around the beam centroid as a linear or quadratic function of ξ. A linear function is a good approximation when the witness beam does not shows a beam that is slightly too short, the variation in the accelerating field is well approximated by a quadratic function. In both cases the field from −3σ ξ0 to 3σ ξ0 was used for fitting.
have sufficient current to flatten the wake and a quadratic function is a good approximation if the witness beam has too much current. To demonstrate this, we ran a series of particle in cell (PIC) simulations using the code VSim [41]. A 0.5 nC witness beam and 2 nC drive beam were propagated through a 3.5×10 16 cm −3 plasma. Both beams had Gaussian current profiles. Simulations were run for witness beams of different lengths, and thus, different peak currents. Fig. 4 shows the longitudinal electric field in the region of the witness beam for two beams with the same charge, but bunch lengths of σ z = 4.0 µm and σ z = 8.0 µm. The electric field around the shorter beam is well represented by a quadratic function centered on the beam centroid while the field around the longer beam can be approximated as linear. We proceed to work out the projected emittance, longitudinal slice emittance, and energy slice emittance growth for both cases.
In the linear case, the accelerating field is given by Here, E 0 is the accelerating field at ξ = 0 and E 1 /∆ξ describes the slope of the electric field; ∆ξ is an arbitrary length scale, we typically set it to the bunch length. The energy of a particle in the beam as a function of s depends on ξ 0 and δ 0 according to . It is straightforward to evaluate Eq. (5) to get the betatron phase advance of each particle If we make the reasonable assumption that E 1 ξ 0 /(E 0 ∆ξ) ≪ 1, then φ can be expanded in the aforementioned quantity giving Even for simple distributions, the C and S integrals are not closed form without further expansion of φ. We use the fact that γ b ≈γ b to expand the square root in the above expression in E 1 ξ 0 /(E 0 ∆ξ) and δ 0 : where we have dropped terms higher than first order. The betatron phase advance of the reference particle is Because φ depends linearly on ξ 0 and δ 0 as ∆φ = −δ 0 ω δ −ξ 0 ω ξ , the I integrals are related to the 2D Fourier transform of the distribution: where Consider a beam with a Gaussian longitudinal distribution with length σ ξ0 and energy spread σ δ0 , The distribution is symmetric sof is real and ψ 1 = ψ 2 = 0. H 1 and H 2 are given by and the projected emittance and moment evolution is given by Eq. (23) and Eqs. (14)- (18), respectively. The longitudinal slice emittance of a Gaussian beam is straightforward to calculate. Evaluating Eq. (24) gives δ0 ω 2 δ /2 and ψ 1 = −ω ξ ξ. H 2 and ψ 2 are given by replacing ω δ with 2ω δ and ω ξ with 2ω ξ and the emittance growth is given by Eq. (23). The longitudinal slice emittance has no dependence on ξ because all the slices have the same initial energy spread. The phase term ψ 1 shows that the beam offset will vary sinusoidally along the bunch as shown in the top of Fig. 5 because ψ 1 ∝ ξ 0 . The oscillation frequency is given by ω ξ .
When finding the energy slice emittance, it will be convenient to use the RMS energy spread:  Fig. 4(a), and the bottom shows the shorter beam from Fig. 4(b). In both cases the initial energy spread is σ δ0 = 4%. The magnitude of the oscillations is dampened by chromatic phase spreading due to the initial energy spread within each slice. The oscillations decohere due to the energy differences the slices pick up as a result of imperfect wake loading. In (a) the accelerating field varies significantly across the slices and the oscillations rapidly decohere. In (b) the accelerating field along the beam is more uniform and the beam oscillates nearly coherently.
Using Eq. (25), we calculate H 1 and ψ 1 for each energy slice: H 2 and ψ 2 are given by replacing ω δ with 2ω δ and ω ξ with 2ω ξ . The emittance growth is given by Eq. (23) and the moments are given by Eqs. (14)- (18). As with the longitudinal slice emittance, the energy slice emittance is the same for all slices because H 1 is independent of δ and the betatron phase offset is linearly proportional to the slice energy because ψ 1 is proportional to δ. Chromatic dephasing cannot occur within an energy slice unless particles are able to mix between slices. The variation in accelerating field along the wake causes this mixing to occur.
There are several emittance saturation length scales. The transverse offset and mismatch driven emittance growth are determined by H 1 and H 2 respectively. We define the saturation length as the distance s when the argument of the exponential in H equals -1. The saturation length for the longitudinal slice emittance growth due to mismatch is given by 2σ 2 δ0 ω 2 δ = 1; this length also quantifies the contribution of the initial energy spread to the emittance growth. The length scale for the emittance growth due to imperfect wake loading is 2σ 2 ξ0 ω 2 ξ = 1. Finally, the saturation length for the energy slice emittance growth is given by The solution for s for each of these expressions is summarized in Table I. To calculate the saturation length for the energy slice emittance, we assumed σ 2 δ0 γ 2 b0 ≪ σ 2 ξ0 s 2 E 2 1 e 2 /(∆ξ 2 m 2 e c 4 ) (i.e. the initial energy spread is smaller than the energy spread induced by the wake), to simplify the expression for σ δ . This assumption is typically well satisfied for a plasma based accelerator with a single stage of reasonable length.
The saturation lengths reveal several properties about the emittance growth. First, the energy slice emittance always has a saturation length longer than that of the longitudinal slice emittance. Consequently, the energy slice emittance is smaller than the longitudinal slice emittance at every s. Second, both slice emittances are approximately independent of σ ξ0 and E 1 ; they depend only on the initial energy spread σ δ0 . Third, if the initial energy spread is too low, the slice emittances will never fully saturate because acceleration reduces the relative energy spread faster than dephasing can occur. Mathematically, this appears as a divergence in the saturation lengths of the slice emittances. Saturation occurs because the spread in betatron phase becomes larger than 2π; the phase spread in a longitudinal slice is ∆φ = σ δ0 ω δ , which does not grow without bound. The maximum phase spread is found by taking the large s limit: The maximum emittance is found by inserting this limit into the expression for H 1 and H 2 . This also means the magnitude of the transverse oscillations of the longitudinal slices are not dampened to zero. The above conclusions apply to the energy slice emittance as long as E 1 is sufficiently large (σ 2 δ0 γ 2 b0 ≪ σ 2 ξ0 s 2 E 2 1 e 2 /(∆ξ 2 m 2 e c 4 )). Figure 6 compares the growth of the longitudinal slice, energy slice, and projected emittance. The accelerating field is the same as that shown in Fig. 4(a), but the beam length has been doubled to σ ξ0 = 16.0 µm to exaggerate the difference between the projected and slice emittance saturation lengths. The initial energy spread is σ δ0 = 0.05. The saturation length due to the variation in the accelerating field is shorter than that due to the beam's initial energy spread. As a result, the saturation length due to imperfect wake loading adequately describes the projected emittance growth. The difference between the numerical and analytic solution for the TABLE I. Saturation lengths for the projected and slice emittances. The projected emittance growth is driven by contributions from both the initial energy spread and the energy spread induced by imperfect wake loading. The shorter of the two saturation lengths dominates the projected emittance growth. In contrast, both energy and longitudinal slice emittance growth depend only on the initial energy spread, assuming the variation in accelerating field is sufficiently strong.

Accelerating field linear in ξ Saturation Length -Offset Beam Saturation Length -Mismatched Beam
Imperfect Loading Contribution a s 1ξ = E0e Accelerating field quadratic in ξ projected emittance arises due to the large final energy spread in this example. This growth can be handled analytically using the approach presented in Ref. [10]; however, for typical experimental parameters the energy spread is small enough that the correction is negligible. For the quadratic case, the accelerating field is given by E z = −E 0 − E 2 ξ 2 /∆ξ 2 , where E 2 /∆ξ 2 is the quadratic fitting parameter. As before, ∆ξ is an arbitrary length scale typically set to the bunch length. The energy of a particle in the beam is whereγ b γ b0 +E 0 se/(m e c 2 ), the same as before. Following the same procedure as before, the betatron phase advance of the particle is given by The reference phase advanceφ is given by Eq. (28). Unlike before, φ does not depend linearly on ξ 0 and thus the C and S integrals are no longer Fourier transforms. To keep the math simple, we define the w ξ as while ω δ is given by Eq. (30). We again assume the Gaussian longitudinal phase space distribution of Eq. (31). In this case, the I integrals are straightforward and the evolution is described by H 2 and ψ 2 are given by making the replacement ω δ → 2ω δ and w ξ → 2w ξ . In this case ψ does not strictly satisfy the requirements for Eq. (23) to be valid; however, 2σ 2 ξ0 w ξ must be small when the emittance is still growing (the emittance saturation length for mismatch is given by σ 2 ξ0 w ξ = 2), allowing us to expand the arctan and satisfy the condition that ψ 2 = 2ψ 1 . To calculate the saturation lengths, we let 1 + 16σ 4 ξ0 w 2 ξ −1/4 ≈ e −1 for mismatch (H 2 ) which corresponds to σ 2 ξ0 w ξ = 2. For transverse offset we assume σ 2 ξ0 w ξ = 4 in H 1 in order to calculate the saturation length. The saturation lengths are shown in Table I. As before, all longitudinal slices have the same slice emittance. The evolution of each longitudinal slice is given by H(ω δ ) = e −σ 2 δ0 ω 2 δ /2 and ψ(w ξ ) = −w ξ ξ 2 . The saturation length is the same as in the linear loading case. If the bunch starts offset, the ξ 2 dependence of the  Table I: black is the contribution from the nonuniform accelerating field to the projected emittance, blue is the initial energy spread contribution as well as the longitudinal slice, and magenta is the energy slice. The projected emittance is dominated by contribution from the non-uniform accelerating field.
phase leads to a large region in the center of the bunch that undergoes transverse oscillations in phase. This is evident in the bottom of Fig. 5.
The energy slice emittance is analytically tractable, but the solution is cumbersome and the H and ψ functions are not easily extracted. To compare to the linear case, we use numerical particle tracking to propagate a witness beam in the accelerating field from Fig. 4(b) and then numerically evaluate Eq. (25) to get the theoretical prediction. For this example, the beam is mismatched but not transversely offset. Fig. 7 shows a comparison of the energy slice spot size evolution between the linear and quadratic cases. In both cases the spot size varies across the energy slices at the exit of the plasma. Unlike the linear loading case, the emittance in the quadratic loading case varies across the energy slices. This variation is described by the H 2 function which is visible in the figure as a δ dependent dampening of spot size oscillation. The low energy slices of the beam disappear because particles can only move to slices with larger δ (if the initial energy spread is uncorrelated).
The variation in beam parameters with energy could be used to measure the wake loading. If the beam is intentionally mismatched into the plasma, the C-S parameters and emittance of each energy slice will have a dependence on the wake loading. An imaging spectrometer can then indirectly measure the loading by looking at the variation in σ x with energy.  Fig. 4(a) and the bottom shows the shorter beam from Fig. 4(b). In both cases the initial energy spread is σ δ0 = 4%. Chromatic dephasing cannot occur within a slice; the dampening of the β function oscillations is due to particle exchange between slices driven by wake loading. If Ez is linear in ξ (top), the emittance grows uniformly for all slices and the beam size only varies with the β function of each slice. If Ez is quadratic in ξ (bottom), the emittance varies across the slices leading to more complex dynamics.

VI. EMITTANCE EVOLUTION IN A PLASMA SOURCE WITH DENSITY RAMPS
In the previous examples, we have only considered plasma sources with uniform density. Here, we analytically calculate the quantities necessary to find the emittance growth in a plasma source with adiabatic density ramps at the entrance and exit. We assume the accelerating field varies with plasma density according to the simple model where η = n e /n eu , n eu is the uniform density in the bulk, E 1 η/∆ξ is the slope of the wakefield in the blowout regime [31,32] and E 0 √ η − 2η describes the variation in longitudinal phase and maximum wake amplitude with plasma density [42]. The fields can be approximated from first principles as where κ is a multiplier accounting for the current of the drive beam. κ typically varies between 1 and 2. Beam loading is ignored in this example to keep the math tractable.
The energy of a single particle is then given by where To find the betatron phase advance in the ramp we assume the energy spread and energy gained in the ramp is small compared to the beam's centroid energy throughout the ramp. We can then expand k β as and ω 2 pu = n eu e 2 /m e ǫ 0 is the plasma frequency of the uniform density region. As before, the betatron phase advance is linear in δ 0 and ξ 0 : ∆φ = −δ 0 ω δ −ξ 0 ω ξ . The I integrals are given by Eq. (29) with The betatron phase advance of the reference particle is where (46) The emittance growth and beam moments are then found by getting H 1 , H 2 , ψ 1 , and ψ 2 from the I integrals and then using Eq. (23) to find the emittance and Eqs. (14)- (18) to find the moments. Without carrying out the full calculation, it is apparent that minimizing the G and D integrals will minimize ω ξ , ω δ , and ∆γ b , i.e., minimize the undesirable impacts the ramp has on the beam.
As the beam enters the uniform density region, it is accelerated significantly and we can no longer assume ∆γ b ≪ γ b0 . Instead, we use the solution presented in Sec. V with the addition of an initial phase and an initial energy of γ 0 → γ bl =γ b (s = l), where l is the length of the ramp. The exit ramp is handled the same way as the entrance ramp except the initial energy is γ bL =γ b (s = l + L), where L is the length of the uniform plasma. The  8. Evolution of the spot size (top) and divergence (bottom) of a mismatched beam entering a plasma source with entrance and exit ramps. The beam is initially focused by the entrance ramp into the bulk of the plasma source where it undergoes chromatic dephasing and becomes matched to the plasma. As the beam is accelerated, the divergence continues to decrease due to adiabatic dampening before the beam is defocused by the exit ramp.
resulting piece-wise functions forφ, ω ξ , and ω δ are fairly cumbersome and are given in Appendix B.
Qualitatively, the solution is similar to that of a uniform plasma presented in Sec. V. The primary difference is the conversion of beam size into divergence by the entrance ramp (focusing) and divergence into beam size by the exit ramp (defocusing). This is shown in Fig. 8 where σ x and σ x ′ are plotted for a mismatched beam propagating through a plasma source with ramps. In this example, the ramps are short and only a small amount of chromatic dephasing occurs in them. The majority of the emittance growth occurs in the bulk plasma.
If the ramps are long or poorly shaped, they can have significant impacts on the witness beam. The integrals G 1 , G 2 , D 1 , and D 2 determine the amount of phase spread, betatron oscillation, and, in the case of a beam driven wake, drive beam energy loss in the ramp. All the integrals are reduced if the integrated plasma density-G 2 -is minimized. The adiabatic ramp that minimizes the impact on the beam therefore has the highest density gradient possible while remaining adiabatic. Solving Eq. (3) gives the optimal ramp as η = 1/(1 − 2α m s) 2 . From experience, the ramp needs |α m | < 0.1 for the ramp to be well described by adiabatic theory. In Fig. 9, we compare several different adiabatic ramp shapes that all focus the beam by the same amount. The ramp shapes that minimize the density integral induce smaller impacts on the witness beam. The full ramp shape is important and the impact of the ramp on the beam cannot be described using only the half width and the adiabatic parameter. Half Widths e −s 2 /2σ 2 , β1 = 64.8, β2 =48.8 e −s/η , β1 =49.7, β2 =29.8 1/(1 +ηs) 2 , β1 = 40.2, β2 =20.0 FIG. 9. Three plasma ramps with different profiles that all reduce the beta function by a factor of 5. The chromatic phase spread, betatron oscillations and drive beam energy loss in the ramp all depend on the integral over the plasma density (G2 shaded area). The first ramp (blue, bottom) is the most efficient adiabatic ramp possible while the Gaussian (red, top) is barely adiabatic at the tail. The ramp shape can have a significant impact on the beam, and the ramp half width is not a good indicator of the effectiveness of the ramp.

VII. PARTICLE INJECTION INTO A LINEAR ACCELERATING FIELD
Consider an injector where a low charge beam is generated within the wake. The following discussion is agnostic to the details of the injection scheme; for example, it could be applied equally well to either ionization injection or plasma photo-cathode injection [43,44]. We can use our theory to extend the theoretical treatment presented by Ref. [33] to include arbitrary distributions of injected electrons and offset of the injected beam with respect to the center of the wake.
Assume at time τ 1 a group of electrons is injected into the wake and rapidly accelerated until each electron reaches a phase locked longitudinal position ξ 0 with energy γ b0 . After the particles are phase locked, the longitudinal distribution of the particles is given by f (ξ 0 , τ 1 ), and the transverse distribution of the particles is described by an initial set of CS parameters β 0 , α 0 , and γ 0 and an initial transverse emittance ǫ n0 . Assume particles ionized at a different time τ 2 are described by the same CS parameters and emittance once they are phase locked. The longitudinal phase space of particles injected at τ 2 is given by f (ξ 0 , τ 2 ).
Further, assume a small amount of charge is injected into the wake, and thus the longitudinal electric field E z varies linearly along the length of the injected beam, E z = −E 0 − E 1 ξ/∆ξ. The energy of an injected particle depends on the injection time and the phase within the wake according to whereγ b = γ b0 + E 0 se/(m e c 2 ) and γ b0 is the energy at the moment of phase locking. The wake is moving at c, so a particle injected at τ is phase locked into position ξ 0 at a location of s 0 = τ c . Using the same approach and assumptions as in Sec. V to evaluate Eq. (5) gives the betatron phase advance of each particle whereφ is given by Eq. (28). The C and S integrals are now in terms of τ rather than δ, we only write C 1 for brevity: Because φ depends linearly on τ and ξ, Eq. (29) holds with δ replaced by τ . ω ξ remains the same as in Eq. (30) while ω τ =k β c.
As pointed out in Ref. [33], the projected emittance initially rises rapidly during injection before decreasing to a minimum and finally growing to saturation. We can solve for the propagation distance that minimizes the emittance growth by finding the value of s where H 1 and H 2 are maximized. Further, we can solve for the propagation distance that minimizes the energy spread.
In an optimal injector, the energy spread and emittance will reach their minimal values simultaneously at the exit of the accelerator.
It is straightforward to solve for the energy spread starting with Eq. (47). Dropping moments of the longitudinal distribution higher than order 2 gives: where σ ξ and σ τ are the standard deviations of the longitudinal distribution about the mean and σ ξτ = ξτ − ξ τ . As long as sE 0 e/m e c 2 ≫ γ b0 , the constant term γ b0 inγ b can be dropped, giving the explicit s dependence of the relative energy spread as The first term is the asymptotic energy spread resulting from imperfect wake loading. The second term is the energy spread induced by the finite injection time. This energy spread is fixed; thus, its contribution to the relative energy spread is suppressed by a factor of 1/γ b as the beam accelerates. The final term accounts for any initial correlation between injection time and energy, particles that are injected later (earlier) experience larger (smaller) accelerating fields, thus flattening the longitudinal phase space. The minimum energy spread is given by The beam reaches its minimum energy spread at which requires σ ξτ > 0 in order for s > 0, which is equivalent to an initial positive chirp. The energy spread can be minimized by creating a large correlation between injection time and injection position. Calculating the value of the minimum emittance requires assuming an initial longitudinal distribution. Take as an example an injection process that injects particles into a blowout wake at a uniform rate from τ = −∆τ /2 to τ = ∆τ /2. After becoming phase locked, the particles injected at a given τ are longitudinally distributed in the wake with a Gaussian distribution centered at vτ . Here, v describes the change in the longitudinal position of the particles with injection time. Positive v means particles injected at later times are injected closer to the back of the wake. The distribution f (ξ 0 , τ ) is separable and can be written as where T = s/c + ∆τ /2 for 0 ≤ s < c∆τ /2 and T = ∆τ for c∆τ /2 ≤ s. The width T accounts for the change in the distribution during injection. The instantaneous injection length σ ξ0 is the bunch length of particles injected at a given τ after they become phase locked. The full bunch has a length of σ 2 ξ = σ 2 ξ0 T 2 v 2 /12. If the injection time is sufficiently long, the energy spread of the beam can rise to to 100% (σ δ ≈ 0.5 for our distribution); then, we can no longer assume β m and γ b are approximately equal for all particles. The minimum achievable energy spread, however, is small and our approach can be used to describe the beam around this minimum. Using Eq. (29) we get Since 2ψ 1 = ψ 2 , the emittance evolution is described by Eq. (23).
In the case of large v and small σ ξ0 (strictly 1 ≪ T v/σ ξ0 + T ω τ /(σ ξ0 ω ξ )), we can approximate the exponential in H 1 and H 2 as 1 and solve for s where both H functions are maximized: which is same distance where the minimum energy spread occurs. This injection distribution works well for the accelerating field assumed here. For even moderate injection duration, however, the minimum emittance is only marginally smaller than the saturated emittance.
Using the estimates for E 0 and E 1 from Eq. (41), we can derive useful formulas for designing a plasma injector. In this case κ can also be used to describe the wakefield phase the particles are injected into. We assume the length of the injector is chosen to minimize the final energy spread of the beam, combining Eq. (53) with Eq. (41) gives The final energy spread of the injected beam is given by and the final energy of the beam is given by The final energy at the point of minimum energy spread is independent of the plasma density and only depends on the correlation between τ and ξ. This occurs because the length of the plasma scales inversely to the plasma density and the correlation v. The final energy spread scales with √ n e σ ξ0 . Thus, the injection region needs to be reduced in proportion to the skin depth to maintain a given energy spread. Some injection schemes have a minimum attainable σ ξ0 resulting in a trade-off between emittance and energy spread because space charge effects, and thus the initial emittance, are reduced as plasma density increases. This trade-off can potentially be mitigated by appropriately loading the wake to reduce the energy spread.
As an example, let us design an injector to produce 2 GeV beams with sub 1% energy spread. Assume κ = 1.5 and γ b0 = 25. Immediately we can find the required correlation is v = c/350, this is reasonable for currently proposed injection schemes [26,27]. The plasma density might depend on the target emittance or the drive beam available; we use a typical experimental value of n e = 5 × 10 17 cm −3 . The plasma should be 12.4 mm long to reach minimum energy spread; to achieve sub 1% energy spread, the particles should be injected such that σ ξ0 ≤ 0.35 µm; we want to emphasize that this is not the final bunch length, but the phase locked length of particles injected at a given τ . Fig. 10 shows the evolution of the longitudinal phase space, energy spread, and projected emittance within the injector. The wake loading effectively cancels out the initial energy chirp of the beam. For parameters of interest to practical injector designs, the minimum emittance is approximately equal to the saturated emittance, thus length is not a concern for ǫ n . The initial spike in emittance is primarily due to the very large energy spread present while injection is still occurring. Notice that while the injector length is not particularly important for the final emittance, it has a significant impact on the final energy spread.

VIII. CONCLUSION
We have demonstrated an analytic approach to calculating the evolution of the witness beam in a plasma based accelerator operating in the nonlinear blowout regime. We included the effects of energy change, loading of the wake, and adiabatic variations in plasma density. We developed our approach to describe the chromatic dephasing of the beam. This dephasing will cause emittance growth of a beam if the beam is either transversely offset or mismatched to the plasma. The growth will saturate if the plasma is sufficiently long. The saturated emittance is the sum of the contribution from the offset and the mismatch, with the transverse offset requiring a longer distance to saturate. In addition, we showed how to calculate both the energy slice and longitudinal slice emittance evolution and saturation values.
For simple cases, the projected and slice emittances can be calculated analytically, letting us investigate general properties of the emittances. Because the particles are phase locked in the wake, the longitudinal slice emittance depends only on the initial energy spread within the slice and grows more slowly than the projected emittance. In the presence of imperfect beam loading, the variation in accelerating field along the beam causes particles to mix between energy slices, leading to growth of the energy slice emittance. Depending on the details of the wake loading, the energy slice emittance will vary across the slices.
In addition to the emittance, our approach provides the beam moments and thus the transverse offset and spot size of the projected beam and the longitudinal/energy slices. Chromatic dephasing leads to a dampening of any transverse offset on the same time frame as the emittance growth. For the energy slices, the mixing process results in an energy dependence of the beam spot size at the exit of the plasma. This dependence is sensitive to the details of the beam loading. This will impact the signal the electron beam generates in an imaging spectrometer, which can be used to indirectly measure the beam loading.
We showed two examples of how our general approach can be applied to specific situations. First, we considered a full plasma accelerator with ramps but with insufficient charge in the witness beam to load the wake. In this case, the energy spread produced by the variation in the accelerating field is sufficient for the emittance to reach saturation regardless of the initial energy spread. We also showed that it is the integral of the plasma density ramp profile that determines how much the ramp perturbs the beam. Second, we considered a plasma injection scheme. For this example, we derived some simple scaling laws for the final energy spread and optimal length of the injector to simultaneously minimize the energy spread and emittance.
In many of our examples, we considered low charge beams that do not significantly load the wake. Our approach, however, is capable of handling more complex loading situations if the longitudinal dependence of the accelerating field can be written analytically. Even if it cannot, the integrals can be solved numerically. Depending on the longitudinal phase space, this may be faster than particle tracking.
The particles are initialized in action-angle variable space (J and ψ) using the distribution The particle's initial position in real space is calculated from J and ψ using x 0 = 2Jβ 0 cos ψ + ∆x, (A5) The longitudinal positions are initialized based on the distribution of interest.