Phase-space moment-equation model of highly relativistic electron-beams in plasma-wakefield accelerators

Annals of Physics 356 (2015) 306–319

Contents lists available at ScienceDirect

Annals of Physics journal homepage: www.elsevier.com/locate/aop

Phase-space moment-equation model of highly relativistic electron-beams in plasma-wakefield accelerators R.E. Robson a,∗ , T. Mehrling b , J. Osterhoff b a

College of Science, Technology and Engineering, James Cook University, Townsville 4811, Australia

b

Deutsches Elektronen-Synchrotron DESY, 22607 Hamburg, Germany

article

info

Article history: Received 27 November 2014 Accepted 6 March 2015 Available online 14 March 2015 MSC: 82D10 35Q83 Keywords: Plasma physics Plasma-based acceleration Vlasov equation Relativistic electron-beams

abstract We formulate a new procedure for modelling the transverse dynamics of relativistic electron beams with significant energy spread when injected into plasma-based accelerators operated in the blow-out regime. Quantities of physical interest, such as the emittance, are furnished directly from solution of phase space moment equations formed from the relativistic Vlasov equation. The moment equations are closed by an Ansatz, and solved analytically for prescribed wakefields. The accuracy of the analytic formulas is established by benchmarking against the results of a semi-analytic/ numerical procedure which is described within the scope of this work, and results from a simulation with the 3D quasi-static PIC code HiPACE. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Plasma-based accelerator technology [1] provides gradients for the acceleration of charged particles in excess of 10 GV/m, which surpasses the fields in conventional accelerators by three orders of magnitude. For this reason, plasma accelerators promise a dramatic miniaturization of electron accelerator technology and accordingly raise hopes for a significant reduction of their cost. The field of plasma acceleration has made remarkable progress in the past decade with the demonstration of mono-energetic spectral features [2–4], the breaking of the GeV-energy barrier [5], controlled

∗

Corresponding author. E-mail address: [email protected] (R.E. Robson).

http://dx.doi.org/10.1016/j.aop.2015.03.004 0003-4916/© 2015 Elsevier Inc. All rights reserved.

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

307

bunch-injection techniques for improved beam quality [6–9], improved stability [10] and high efficiency [11]. Also, first applications for photon generation from plasma-generated electron beams have been demonstrated [12,13]. Despite this tremendous development, relativistic electron beams from plasma accelerators are not yet achieving the phase-space quality required for more demanding applications such as free-electron lasers and high-energy physics colliders. Further progress in this area critically depends on obtaining a better understanding of, and thus, control over the key processes, such as those influencing transverse beam properties. It is here that theoretical modelling can play a significant role. PIC (particle in cell) codes [14] have hitherto provided the most detailed insight, yet, require significant computational resources. Analytical modelling offers in many cases a complementary approach without the need for extreme computation. The procedures outlined in this paper have been developed to that end, and provide an analytical way of describing the evolution of the transverse phase space of an electron beam propagating in a plasma wave along a longitudinally homogeneous plasma-profile. In addition, a numerical method is proposed, which models the transverse dynamics of an electron beam in a plasma wakefield in inhomogeneous plasma in a computationally efficient manner. We start with the observation that the quantities of experimental interest may be expressed as phase space averages or ‘‘moments’’ of the electron phase space distribution function. The latter may be found as the solution of the relativistic Vlasov kinetic equation, but we prefer to work instead solely with equations for the moments themselves, which are found by forming appropriate averages of the Vlasov equation. In other words, we seek to find the averages without the need to first find the distribution function. The price of taking this ‘‘short cut’’ (see Fig. 1) is that it produces more unknowns than equations, the familiar problem of closure, and the challenge is then to find an Ansatz (or postulate) which closes the moment equations while still producing acceptably accurate solutions. As with fluid modelling in plasma physics [15], an Ansatz built on sound physical reasoning, rather than ad hoc supposition, is an essential ingredient for success. In this paper we consider a relativistic electron bunch with energy spread which is injected into wakefields, driven by a laser or particle beam in the blowout regime [16]. The field configuration in this regime is featured by focusing fields which are proportional to the radius from the propagation axis of the driving beam [17,18]. Moment equations are built from first principles, and solved rigorously for the case where the driver beam does not evolve and the focusing fields remain constant along the direction of driver propagation. This is a valid assumption if the injected beam is ultra-relativistic, such that the betatron wavelength is long compared to the transition from the vacuum to the plasma, and if the driver beam does not evolve significantly during the time over which the expected quality degradation of the witness beam occurs. Analytic expressions for all physical properties of interest, including the emittance, are thereby obtained. In the case where the field varies along the axis of propagation, a semi-analytic/numerical approach is used to furnish the required expressions. The aim is to obtain a computationally economic method for modelling the transverse beam dynamics. The analytic results are benchmarked against the semi-analytic/numerical procedure and results from the quasi-static three-dimensional (3D) PIC code HiPACE [19]. Other procedures aimed at providing an analytic perspective [20,21] start at the single particle level, and proceed to an ensemble average through a series of steps involving an assumed form of the phase space distribution function. In contrast, we consider average properties from the outset, as furnished by phase space moments of the Vlasov equation [22] without stipulating the shape of the phase space distribution explicitly. It is emphasized that the procedure is completely self-contained and we do not need to solve the Vlasov equation for the distribution function itself. In passing we note a similar treatment by Kumar for low energy, non-relativistic electron ‘‘swarms’’ in gases starting with the classical Boltzmann kinetic equation [23]. The structure of this paper is as follows: In Section 2, we discuss the formalism for the phase space moment method and establish the general form of the moment equation. This is followed by discussions of the solution of the equations, both general aspects and the specific details, including the semi-analytic/numerical procedure and the Ansatz, in Sections 3 and 4 respectively. In Section 5, we compare the analytic results with semi-analytic/numerical values and results obtained from a PIC simulation, thereby establishing the credentials of the moment method and the accuracy of the Ansatz. Section 6 summarizes the results and points the way to further applications.

308

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

Fig. 1. Two possible ways of obtaining quantities of physical interest as phase space averages: (a) Vlasov’s kinetic equation is solved for the electron phase space distribution function f (r, p, t ) and the required moments then are found by integration (upper path); (b) The ‘‘short cut’’ method, discussed in this paper, which furnishes the averages directly from low order moment equations which are formed from the Vlasov equation, plus a closure Ansatz, without any need to know f (lower path).

Although the idea is conceptually straightforward, there is good deal of mathematical working required to arrive at the desired result. In order to avoid any possibility of distraction from the main goal, most of these details have been placed in the Appendix. 2. Phase space averaging procedures 2.1. The relativistic Vlasov kinetic equation The properties of an electron beam (or ‘‘bunch’’) in phase space at time t may be completely described by a density distribution function f = f (r, p, t ), from which the quantities of experimental interest may be found by taking appropriate averages or ‘‘moments’’. If binary interactions between electrons can be neglected, then f is the solution of the Vlasov (or collisionless Boltzmann) kinetic equation,

(∂t + v · ∇ + ∂p · F) f = 0,

(1)

1 + (p/mc )2 , c is the speed of light, and where v = p/(γ m), m is the electron rest mass, γ = F = −e [E + v × B] is the force acting on an electron due to electric and magnetic fields E and B. In this work, it is assumed that self-consistent fields are negligible, and the fields are externally prescribed. A system of coordinates (x, y, z) is chosen such that the z-axis is parallel to the direction of laser propagation, and properties in the transverse direction are specified by the vector x⊥ = (x, y, 0). The Vlasov equation may be transformed by introducing a dimensionless momentum u = p/mc, and a coordinate ξ ≡ z − ct measured with respect to reference frame whose origin is at s = ct:



(∂s + (βz − 1) ∂ξ + β⊥ · ∇⊥ + ∂u · F/mc 2 ) f = 0.

(2)

Here βz = vz /c and β⊥ = v⊥ /c are dimensionless velocities in the longitudinal and transverse directions respectively, ∇⊥ = (∂x , ∂y , 0), and f = f (x⊥ , ξ , u, s) is now to be regarded as a distribution in the transformed phase space coordinates x⊥ , ξ , and u. 2.2. Averaging and moment equations 2.2.1. Moment equations We now averaging over all momenta u and transverse spatial coordinates x⊥ . Thus  consider  dN = dξ dx⊥ du f (x⊥ , ξ , u, t ) gives the number of electrons, while

 N =

dξ



 dx⊥

du f (x⊥ , ξ , u, s),

(3)

is the corresponding particle number. Since the beams under consideration are highly relativistic with γ ≫ 1, then βz = vz /c → 1, and the term in (2) involving the derivative with respect to ξ vanishes. The Vlasov equation reduces to

(∂s + β⊥ · ∇⊥ + ∂u · F/mc 2 ) f = 0,

(4)

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

309

where, since there is effectively one less coordinate, f = f (x⊥ , u, s). Note that integration of Eq. (4) over all x⊥ and u gives

∂s N = 0

(5)

showing that N is a conserved quantity. Multiplying Eq. (4) by Φ (x⊥ , u) and integrating over x⊥ and u furnishes the moment equation

∂s ⟨Φ ⟩ =



x′⊥

∂ · ∇⊥ Φ + F⊥ /mc · Φ + Fz /mc Φ ∂ u⊥ ∂ uz 



2

∂





2

 (6)

where an average is defined by 1





du Φ (x⊥ , u) f (x⊥ , u,s), n and the slope of a particle trajectory in the transverse directions is given by

⟨Φ ⟩ (s) =

dx⊥

x′⊥ ≡ u⊥ /uz = β⊥ /βz ≃ β⊥ .

(7)

(8)

2.2.2. Cylindrically symmetric case If the electron bunch is symmetric about an axis which coincides with the direction of laser propagation, i.e. the z-axis, the x and y-directions are equivalent, and only one coordinate, say x, is needed to describe variation in the transverse direction. Similarly only the x-components of the transverse vectors x′⊥ and F⊥ , denoted by x′ = px /pz and Fx respectively, enter into the equations. The general moment equation (6) then simplifies to

  ∂s ⟨Φ ⟩ = x′ ∂x Φ +



Fx mc

   Fz ∂ Φ + ∂ Φ . u u x z 2 2

(9)

mc

2.2.3. Specification of the force terms The force terms are assumed to be as follows: Fx

Fz

= f∥ (10) mc 2 where both kx and f∥ are independent of x, but may, in general, be functions of s. After making these substitutions we obtain, mc 2

= f⊥ = −kx x ;

      ∂s ⟨Φ ⟩ = x′ ∂x Φ − kx x ∂ux Φ + f∥ ∂uz Φ

(11)

a moment equation which forms the basis for all subsequent calculations. 3. The fundamental moment equations and emittance 3.1. Moment equations for small relative variations of uz The quantities of physical interest are the averages x2 , xx′ and x′2 (or equivalently, Courant– Snyder parameters [24]) and the corresponding moment equations are found by setting Φ = x2 , xx′ , x′2 in succession into Eq. (11):

  

    ∂s x2 = 2 xx′







(12)



    x ∂s xx′ = x′2 − kx

2

uz



 − f∥

xx

′

uz

 ′  ′2    xx x ∂s x′2 = −2kx − 2f∥ . uz

uz

(13)

(14)

When analysing the magnitude of the terms on the right hand side  of Eqs. (13) and (14), the following can be noted. In a blow-out regime plasma wave |f⊥ | ∼ f∥ , and moreover for the slope of the

310

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

transverse trajectories of beam particles x′ = βx ≪ 1. Hence,

   ′     2       xx   ≪  f⊥ x  = kx x  f∥  u    u  uz  z z   ′2    ′    ′   x    x  xx  f∥       u  ≪ f⊥ u  = kx u  . z

z

(15)

(16)

z

Thus, the terms containing f∥ in Eqs. (13) and (14) may be neglected, and hence:

 2     x ∂s xx′ = x′2 − kx

(17)

 ′   xx . ∂s x′2 = −2kx

(18)

uz

uz

In addition, we assume ⟨uz ⟩ to experience only an insignificant relative change during the considered propagation length L,

∂s ⟨uz ⟩ L ≈ 0. ⟨uz ⟩

(19)

This assumption corresponds to a configuration with a pre-accelerated electron beam, which has a significant energy at injection compared to the energy gain during the acceleration. Such a scenario is realistic for staged plasma-based acceleration schemes. The assumption will be implicit in what follows when discussing the solution of the three basic moment equations (12), (17) and (18). 3.2. Normalized transverse trace-space emittance and its evolution The phase space quality of an electron bunch may be characterized by the normalized transverse trace-space emittance (see e.g. [20]), ϵn = ⟨uz ⟩ ϵ , where

ϵ=

  

x2 x′2 − ⟨xx′ ⟩2



(20)

      is the transverse trace-space emittance, which may be evaluated after solving for x2 , xx′ and x′2 . Note that a general expression for the evolution  of  the emittance   can  be obtained by differentiation of (20) with respect to s and substitution for ∂s x2 , ∂s xx′ and ∂s x′2 from Eqs. (12), (17) and (18):        xx′ k x  ′  x2 ∂s ϵn ≈ ⟨uz ⟩ ∂s ϵ = ⟨uz ⟩ xx − x2 . (21) ϵ uz uz 3.3. Closure through an Ansatz The three moment (12),  equations   (17) and (18), are not closed, and therefore cannot be solved, until the terms x2 /uz and xx′ /uz are specified by an Ansatz. Thus for example, we might assume that an average of a quotient is the quotient of averages, i.e.,



x2 uz

 2

 ≈

x

⟨uz ⟩

 ;

xx′ uz



 ≈

xx′



⟨uz ⟩

(22)

to obtain a closed set of equations, which could be solved for any kx (s). However, this Ansatz effectively   assigns the same value of uz to each electron and thereby ignores any ‘‘spread’’ (uz − ⟨uz ⟩)2 in uz about the mean ⟨uz ⟩. This leads to unphysical results, e.g., substitution of Eq. (A.1)   in Eq.  (21) gives  ∂s ϵn ≡ 0, i.e., a non-evolving emittance, regardless of the explicit expressions for x2 , xx′ and x′2 . It is therefore important to correctly account for the spread, which may be measured by the parameter,

 δ=

(uz − ⟨uz ⟩)2 ⟨uz ⟩

 (23)

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

311

and to formulate an Ansatz with δ > 0. The naive Ansatz (A.1) assumes δ = 0 and cannot by itself reproduce the correct physics. Nevertheless, it can be used to provide the basis for a useful numerical procedure, as outlined below. 3.4. A semi-analytic/numerical procedure Here we outline a semi-analytic/numerical procedure to compute the averages of beams with significant energy spread for all times and for arbitrary profiles of kx (s). In Eqs. (12), (17) and (18) the average ⟨⟩ is understood to be over all x, x′ and uz . Suppose now that we consider a subset of particles which all have the same uz , or  equivalently, that the averaging does not yet include uz . The corresponding averages, written as x2 u etc., are formed at a particular value z of uz and satisfy moment equations which formally follow from (A.1), with ⟨uz ⟩ = uz :   d  2 x u = 2 xx′ u z z ds

 

x2 u   d  ′ z xx u = x′2 u − kx z z ds uz



(24)



xx′ u d  ′2  z x u = −2kx . z ds uz These must in general be solved numerically for each value of uz , for a given kx (s). Now suppose that the uz are distributed about a specified mean value ⟨uz ⟩ according to an assumed normalized distribution function f (uz ). Multiplying Eq. (24) by f (uz ) and integrating over all the allowed values of uz gives Eqs. (12), (17) and (18) again, with

 2



x

=

x2



duz f (uz ) x2

 

(25)

uz

and



uz

 =

 2 duz f (uz )

x

uz

uz

and similarly for the other quantities. In general, the integrals (25) must be evaluated numerically. In summary, the procedure allows x2 , etc., (average over all variables) to be found by first solving

 

the differential equations (24) for x2 u (zero spread, average over only some variables) and then z applying the rule (25). In a sense, the latter is equivalent to an Ansatz. This approach is compared later with the analytic results obtained below. 4. Solution of the moment equations 4.1. Analytic solution for constant kx In plasma-based accelerators operated in the blow-out regime, kx is only a function of the plasma density. Hence, for a plasma target with longitudinally homogeneous profile, kx is constant. In this case, the solution of the moment equations for beams with energy spread can be found by employing a Laplace transformation (see Appendix for full details):

   = A + exp −(δ kβ s)2 /2 B cos 2kβ s + C sin 2kβ s  ′    xx = kβ exp −(δ kβ s)2 /2 C cos 2kβ s − B sin 2kβ s  ′2      x = k2β A − exp −(δ kβ s)2 /2 B cos 2kβ s + C sin 2kβ s  2 x

(26)

where

 kβ =

kx

⟨uz ⟩

(27)

312

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

 

corresponds to the betatron wavenumber. The constants are given in terms of initial conditions x2 0 ,





 ′

xx

0

and x

 ′2

by

 20

A=

x



0

x′ 2

+

2 x

−

2

 C =

 0

xx′

0

2k2β

 2 B=



x′ 2

 0

(28)

2k2β



kβ

0

.

Substitution of Eq. (26) in Eq. (20) gives for the transverse trace-space emittance:

 ϵ = kβ A2 − exp[− (δ kβ s)2 ](B2 + C 2 ).

(29)

The accuracy of these expressions, and hence the Ansatz itself, is determined by benchmarking against numerical results which are obtained by use of the method described in 3.4. Some general points to note: (a) For the special case of ‘‘matching’’ initial conditions,

  2 x

0

=

x′ 2



2

kβ

0



;

xx′

 0

=0

(30)

it follows from Eq. (28) that B = C = 0, and hence all quantities remain constant. (b) In the limit of zero spread, δ → 0, exp[−(δ kβ s)2 /2] → 1 and Eq. (26) then describe pure undamped harmonic behaviour. Eq. (29) also shows that emittance is constant. The same results could be obtained more directly by applying the Ansatz (A.1) to the moment equations (12), (17) and (18). As noted earlier, a spread in momentum spread is essential for producing an emittance which evolves in time. (c) In the limit where both kx and kβ → 0, Eqs. (26) and (29) then yield the expression for the moments,

      = x2 0 + 2 xx′ 0 s + x′2 0 s2  ′  ′   xx = xx 0 + x′2 0 s  ′2   ′2  x = x 0  2 x

(31)

and emittance,

ϵ=

   x2

0

x′ 2

 0

− ⟨xx′ ⟩20

(32)

respectively. The same expressions could also be obtained directly from Eqs. (12), (17) and (18) with kx = 0. Note that in the asymptotic region, quantities should vary according to (31) and (32), although the ‘‘initial conditions’’ correspond now to the values attained emerging from the transition region, and are therefore not prescribed as such. 5. Numerical calculations Fig. 2 shows a comparison between moments calculated from the analytic expressions (26), corresponding values found from the numerical procedure described in Section 3.4, and results obtained from a 3D Particle-In-Cell (PIC) simulation with the quasi-static PIC code HiPACE [19] for typical parameters kx and δ. The simulation parameters for the 3D PIC simulation with a reference density of n0 = 1017 cm−3 are chosen as follows. The dimensions of the co-propagating simulation box in the longitudinal and transverse directions are 201.7 × 235.3 × 235.3 µm3 with 512 × 256 × 256 cells with 8 particles per

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

313

Fig. 2. Evolution of beam phase-space moments after injection into a blow-out plasma wakefield. Solid lines depict the evolution according to the analytic description, dashed lines show the evolution according to results obtained by use of the semi-analytical approach described in 3.4, and dash-dotted lines show the result  of  a PIC simulation with the code HiPACE.

 

Initial beam parameters were chosen as follows: x2

0

  = (5 µm)2 , xx′ 0 = 0, x′ 2

The betatron wavenumber in this configuration was chosen to be kβ = 0.94 mm n0 = 1017 cm−3 (see e.g. [25]).

0

−1

= (0.5 mrad)2 , ⟨uz ⟩ = 2000, δ = 0.1.

, in accordance to a plasma density of

cell for the drive beam, a total number of 5 × 104 particles for the witness beam and 8 particles per cell for the plasma in the vicinity of the drive beam, down to 1 particle per cell for the plasma near to the box boundary. The current density deposition is linear in the longitudinal direction and quadratic in the transverse directions. The drive beam has a Gaussian profile with a peak electron density of 4.0 × n0 and longitudinal and transverse rms values of 11.8 µm, 8.4 µm and 8.4 µm. The mean energy corresponds to a Lorentz factor of 4 × 103 , the relative energy spread is 0.001, and uncorrelated, and the initial normalized trace-space emittance is 13.3 µm in both transverse directions. The witness beam has a Gaussian profile with a peak electron density of 0.5 × n0 and longitudinal and transverse rms values of 1.68 µm, 5 µm and 5 µm. The initial mean energy of the witness beam corresponds to ⟨uz ⟩0 = 2000, the initial relative energy spread is δ0 = 0.1 and uncorrelated, and the initial normalized trace-space emittance is 5 µm in both transverse directions. The spatial offset between driver and witness beam is 58.8 µm. Such parameters are close to what can be obtained by modern electron accelerators. The plasma density is uniform at 1017 cm−3 transversely and along the propagation axis. The agreement between the analytically and semi-analytically/numerically obtained curves in Fig. 2 is excellent over the entire range of s considered. Similarly, the normalized emittance calculated from the analytic expression Eq. (29) agrees very well with the semi-analytic/numerical procedure, apart from some oscillatory structure. This fine oscillatory structure evident in the numerical calculations of emittance is absent in the analytic solution because of a small approximation made when inverting the Laplace transform (see A.4). Similar excellent agreement has been found for a wide range of values of kx and δ .

314

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

While no change in the energy spread and no change of the mean energy were assumed in the basis for the analytic approach and (although possible) for the computation of the semi-analytic/numerical approach, these effects are included in the PIC simulation. The longitudinal momentum of the witness beam after 5 cm of acceleration in the PIC simulation was at ⟨uz ⟩ = 2795, and δ = 0.0777. The effect 2 of the energy gain can be seen in the curves for ⟨x′ ⟩ and ⟨xx′ ⟩ obtained in the PIC simulation, which are adiabatically changed since x′ = ux /uz . The witness beam in the PIC simulation does not load [26] the plasma wave significantly, and has a finite length, hence different parts of the beam are accelerated at a different rate and the energy spread increases during the acceleration, while the relative energy spread δ in this example decreases. This results in the deviations in the decoherence rate seen in the curves for the moments and for the normalized emittance from the PIC simulation and from the analytic and semi-analytic/numerical approach. Despite the neglect of these effects in the analytic approach, it should be noted that the curves for the normalized emittance ϵn from the PIC simulation, from the analytic approach and from the semi-analytic/numerical approach are close with respect to each other. 6. Discussion and concluding remarks In this paper the main aim has been to introduce the moment equation method for calculation of the properties of a relativistic beam of electrons with energy spread. The Vlasov kinetic equation for the phase space distribution function f has been integrated to furnish a set of moment equations, which must, however, be closed through an Ansatz. The aim is to obtain quantities of physical interest directly as solutions of the moment equations, without the need to find f . To simplify the analysis as much as possible we have assumed that: (a) Mutual interactions between electrons are negligible; (b) Space-charge effects are negligible, and any fields are externally prescribed; and (c) The beam is axially symmetric. In order to obtain an analytic solution, it has also been assumed that: (d) The radial focusing field varies with radial position, but not with axial position s. The accuracy of the analytic expressions for moments and normalized emittance thus obtained has been established by comparison with the results of an independent semi-analytic/numerical procedure described in Section 3.4 and simulation results obtained with a three dimensional PIC code. The above limitations can be relaxed: for example, space-charge fields can be included by solving the moment equations in conjunction with a correspondingly averaged form of Maxwell’s equations. However, the priority as we see it is to relax the restriction (d), and allow for a focusing field which varies with s, and there is particular interest in the case where kx (s) decreases with s through a transition region, tending a small value asymptotically. In this asymptotic region, where kx (s) → 0, we know that quantities must vary with s according to Eqs. (31) and (32), but the ‘‘initial conditions’’ are prescribed by the emergence of the beam from the transition region, and are therefore not known at the outset. The method of Laplace transformation outlined in this paper for the case of constant kx produces accurate analytic expressions. When benchmarked against results from semi-analytical studies and results from PIC simulations, this analytic procedure shows an excellent agreement of the normalized emittance evolution for a variety of parameters. It therefore allows for the quantitative prediction of the emittance growth of electron beams, externally injected into a plasma wakefield operated in the blow-out regime. This is of vital interest for the study of the quality-preservation of beams with finite energy spread in such scenarios and hence for the improvement of the beam-quality in plasma-based accelerators. Acknowledgements We gratefully acknowledge the support of the Alexander von Humboldt Foundation, Helmholtz Association VI-503 and EuCARD2 WP7. One of us (RER) wishes to warmly thank members of the FLA group at DESY for their hospitality during the period 2011–2014.

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

315

Appendix. Analytic procedure for solving the moment equations for constant kx Here we outline the procedure solving the three basic Eqs. (12), (17) and (18) for ⟨x2 ⟩, ⟨xx′ ⟩, ⟨x′2 ⟩, under the assumption that both kx and ⟨uz ⟩ are constant. The equations are closed by an Ansatz, which allows for the all-important spread in longitudinal momenta. Since the full analysis is quite long and algebraically involved, only the most important details are given. A.1. A reformulation We start with a reformulated set of moment equations. Let U be an arbitrary constant, and substitute Φ = x2 eU /uz , xx′ eU /uz and x′2 eU /uz successively in Eq. (11), thereby generating moment equations for X (s, U ) ≡ ⟨x2 eU /uz ⟩;

Y (s, U ) ≡ ⟨xx′ eU /uz ⟩;

Z (s, U ) ≡ ⟨x′2 eU /uz ⟩.

(A.1)

Since

 U / uz  e ∂X = x2 ; ∂U uz

  ∂Y eU /uz = xx′ ∂U uz

(A.2)

the new set of moment equations may be written as

∂s X = 2Y ∂s Y = Z − kx

∂X ∂U

(A.3)

∂Y ∂s Z = −2kx ∂U

where terms involving the longitudinal force have been neglected as before. These equations are to be solved for X (s, U ), Y (s, U ), Z (s, U ) for arbitrary U, with assumed initial conditions (herein lies the Ansatz), and finally one sets U = 0 to obtain the quantities of physical interest, i.e., X (s, 0) ≡ ⟨x2 ⟩;

Y (s, 0) ≡ ⟨xx′ ⟩;

Z (s, 0) ≡ ⟨x′2 ⟩.

(A.4)

It is not necessary to return to the original moment equations (12), (17) and (18), although these could be regained simply by setting U = 0 in (A.3). Likewise, although it is not now necessary, we could obtain expressions for the unknown terms in Eqs. (17) and (18), by setting U = 0 in Eq. (A.2):



x2 uz



 =

∂X ∂U



 ; U =0

xx′ uz



 =

∂Y ∂U



.

(A.5)

U =0

To this point, the discussion is valid for any kx (s). However, in the following we will assume the case of a constant kx . In this case an analytic solution is possible. A.2. Laplace transformed moment equations for constant kx We start by taking the Laplace transform of Eqs. (A.3): pX = 2Y + X0 pY = Z − kx pZ = −2kx

∂X + Y0 ∂U

∂Y + Z0 ∂U

where an overhead bar denotes transformed quantities, X = X ( U , p) ≡

∞



e−ps X (U , s) ds 0

(A.6)

316

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

and a subscript 0 denotes the initial value, X0 = X0 (U ) ≡ X (U , 0)

(A.7)

and similarly for Y and Z . Note that these initial conditions are in fact known only for U = 0, through the relationship with the known initial values of the actual physical moments of interest, e.g., X0 (0) = ⟨x2 ⟩0

(A.8)

and we have to make a postulate for X0 (U ) for U ̸= 0. The solution of (A.6) may be found exactly, e.g., X ( U , p) =

X0 ( U ) 2p

p2 U

 + exp −



U



′

dU exp −

K ⟨uz ⟩

p2 U ′ K ⟨uz ⟩



f0 (U ′ )

(A.9)

where (here and elsewhere) the indefinite integral is implied in the right hand side, f0 ( U ) =

1 K ⟨uz ⟩



p 2

X0 (U ) + 2Y0 (U ) +

2 p



Z0 (U )

(A.10)

and K ≡ 2k2β = 4kx /⟨uz ⟩. The Laplace transform of the moment of physical interest then follows by setting U = 0 in (A.9). Using (A.8) and making the substitution U ′ = iu in the integral, we find

⟨x2 ⟩ = X (0, p) =

⟨x2 ⟩0 2p



0

+i

 du exp −

ip2 u



K ⟨uz ⟩

f0 (iu).

(A.11)

The corresponding expressions for Y (U , p) and Z (U , p), and hence ⟨xx′ ⟩ and ⟨x′2 ⟩ respectively, may be found similarly. Eq. (A.11) is exact and, in principle, one has only to invert the Laplace transform to get an exact expression for ⟨x2 ⟩. However f0 (iu) is not specified, since the initial values X0 (iu), Y0 (iu) and Z0 (iu) are not known for arbitrary u. One has to somehow connect these with the known initial values ⟨x2 ⟩0 , ⟨xx′ ⟩0 , ⟨x′2 ⟩0 of the actual physical moments, and that requires a knowledge of the initial distribution of uz .The nature of the Ansatz has thus been shifted from a requiring postulate for the unknown terms in Eqs. (17) and (18) to one relating to initial conditions. A.3. The Ansatz The unknown quantity f0 (iu) involves averages like X0 (iu) = ⟨x2 eiu/uz ⟩0

(A.12)

which are effectively integrals of eiu/uz with the initial probability  distribution in uz . The latter is assumed to be a Gaussian centred on ⟨uz ⟩, sharply peaked so that ⟨(uz − ⟨uz ⟩)2 ⟩0 is small compared with ⟨uz ⟩ and hence we may make the approximation eiu/uz ≈ eiu/⟨uz ⟩−iu(uz −⟨uz ⟩)/⟨uz ⟩ . 2

The average (A.12), which is an integral over uz , is thus effectively the Fourier transform of a Gaussian, which itself is a Gaussian in u: X0 (iu) = ⟨x2 eiu/uz ⟩0 = exp[iu/⟨uz ⟩ − u2 /2β 2 ]⟨x2 ⟩0

(A.13)

where

⟨uz ⟩2 β=  . ⟨(uz − ⟨uz ⟩)2 ⟩0

(A.14)

While the rational for the Ansatz (A.13) is clear, its accuracy can only be tested subsequently by appropriate benchmarking (see Section 5).

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

317

Similar reasoning applies for the other moments, and hence Y0 (iu) = ⟨xx′ eiu/uz ⟩0 = exp[iu/⟨uz ⟩ − u2 /2β 2 ]⟨xx′ ⟩0 Z0 (iu) = ⟨x′2 eiu/uz ⟩0 = exp[iu/⟨uz ⟩ − u2 /2β 2 ]⟨x′2 ⟩0 . Thus, after substitution into Eq. (A.10), we find f0 (iu) =

exp[iu/⟨uz ⟩ − u2 /2β 2 ]



K ⟨uz ⟩

p 2

2

′

2



′2

⟨x ⟩0 + 2⟨xx ⟩0 + ⟨x ⟩0

(A.15)

p

and hence Eq. (A.11) gives i

⟨x2 ⟩0

⟨x2 ⟩ =



+

2p

p 2

⟨x2 ⟩0 + 2⟨xx′ ⟩0 + 2p ⟨x′2 ⟩0

 (A.16)

K ⟨uz ⟩

0



du exp iu 1 + p2 /K /⟨uz ⟩ − u2 /2β 2

 

×





.

(A.17)

Using the identity (Appendix B of Ref. [27]) the integral may be conveniently written as



iβ

0

dk exp iu 1 + p /K /⟨uz ⟩ − u /2β



2



2



2



= √ Z 2



1 + p2 /K



√

2δ

(A.18)

where

δ≡

⟨uz ⟩ = β



⟨(uz − ⟨uz ⟩)2 ⟩ ⟨uz ⟩

(A.19)

is a parameter measuring the spread in longitudinal momenta, and

π

2

∞

e−x

−∞

x−ζ



1 Z (ζ ) ≡ √

dx

(A.20)

is the plasma dispersion function [28]. Eqs. (A.17) and (A.18) together give

⟨x2 ⟩ =

F0 (p) − √ Z 2p 2δ

⟨x2 ⟩0



1 + p2 /K

 (A.21)

√

2δ

where F0 (p) =

1



K

p 2

 2 ⟨x2 ⟩0 + 2⟨xx′ ⟩0 + ⟨x′2 ⟩0 .

(A.22)

p

Eq. (A.21) is mathematically exact, to the extent that the Ansatz (A.13) holds, but no further approximations have been made. Expressions for the other moments may be similarly expressed in terms of 1+p2 /K ). Z( √ 2δ

A.4. Inverting the Laplace transforms From (A.20) we have

 Z

1 + p2 / K

√



2δ

where K ′ ≡ K (1 −

√  2 −K 2δ ∞ e−x = √ dx 2 ′ π −∞ p + K

√

2δ x).

The last expression holds to a very good approximation, since the spread δ is small, and the integrand is significant only if x . 1. After factorization of the integrand, 1 p2 + K ′

1

=− √

2i K ′



1

√

p + i K′

−

1

√

p − i K′



318

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

the Laplace transform may be readily inverted,

 ∞  √ ′ √  Kδ 2 ′ = √ √ e−x e−i K s − ei K s dx 2δ 2i π −∞ √  ∞  √ √ √ √ √ √  Kδ 2 e−x e−i K s+ixδ K s/ 2 − ei K s−ixδ K s/ 2 dx. ≈ √ √ 2i π −∞ √ √ 1 In the last step we have made the approximation (1 − 2δ x) 2 ≈ (1 − δ x/ 2), which can be justified on the grounds that δ is small, and that the integrand is significant only for x < 1.  

L−1 Z

1 + p2 /K



√

Since



∞

2 e−x ± ixδ

√

√

K s/ 2

dx =

√ π exp(−K δ 2 s2 /8)

−∞

then L

−1

  Z

1 + p2 /K



√

2δ

√ √ 2 2 ≈ − 2K δ e−K δ s /8 sin K s.

Similarly



L−1 pZ L−1



1 p





√ √ 2 2 ≈ − 2K δ e−K δ s /8 cos K s



√ √ 2 2 ≈ − 2δ [1 − e−K δ s /8 cos K s].

√

2δ

 Z

1 + p2 /K 1 + p2 /K

√

2δ

These expressions are sufficient to enable us to find the inverse of the Laplace transformed moment (A.21). Laplace transforms of the other moments may be similarly inverted, and the complete solution Eq. (26) follows. References [1] E. Esarey, P. Sprangle, J. Krall, A. Ting, IEEE Trans. Plasma Sci. 24 (2) (1996) 252–288. http://dx.doi.org/10.1109/27.509991. [2] C.G.R. Geddes, C. Toth, J. van Tilborg, E. Esarey, C.B. Schroeder, D. Bruhwiler, C. Nieter, J. Cary, W.P. Leemans, Nature 431 (7008) (2004) 538–541. http://dx.doi.org/10.1038/nature02900. [3] S.P.D. Mangles, C.D. Murphy, Z. Najmudin, A.G.R. Thomas, J.L. Collier, A.E. Dangor, E.J. Divall, P.S. Foster, J.G. Gallacher, C.J. Hooker, D.A. Jaroszynski, A.J. Langley, W.B. Mori, P.A. Norreys, F.S. Tsung, R. Viskup, B.R. Walton, K. Krushelnick, Nature 431 (7008) (2004) 535–538. http://dx.doi.org/10.1038/nature02939. [4] J. Faure, Y. Glinec, A. Pukhov, S. Kiselev, S. Gordienko, E. Lefebvre, J.-P. Rousseau, F. Burgy, V. Malka, Nature 431 (7008) (2004) 541–544. http://dx.doi.org/10.1038/nature02963. [5] W.P. Leemans, B. Nagler, A.J. Gonsalves, C. Toth, K. Nakamura, C.G.R. Geddes, E. Esarey, C.B. Schroeder, S.M. Hooker, Nature Phys. 2 (10) (2006) 696–699. http://dx.doi.org/10.1038/nphys418. [6] J. Faure, C. Rechatin, A. Norlin, A. Lifschitz, Y. Glinec, V. Malka, Nature 444 (7120) (2006) 737–739. http://dx.doi.org/10. 1038/nature05393. [7] A.J. Gonsalves, K. Nakamura, C. Lin, D. Panasenko, S. Shiraishi, T. Sokollik, C. Benedetti, C.B. Schroeder, C.G.R. Geddes, J. van Tilborg, J. Osterhoff, E. Esarey, C. Toth, W.P. Leemans, Nature Phys. 7 (11) (2011) 862–866. http://dx.doi.org/10. 1038/nphys2071. [8] A. Pak, K.A. Marsh, S.F. Martins, W. Lu, W.B. Mori, C. Joshi, Phys. Rev. Lett. 104 (2010) 025003. http://dx.doi.org/10.1103/ PhysRevLett.104.025003. [9] C.E. Clayton, J.E. Ralph, F. Albert, R.A. Fonseca, S.H. Glenzer, C. Joshi, W. Lu, K.A. Marsh, S.F. Martins, W.B. Mori, A. Pak, F.S. Tsung, B.B. Pollock, J.S. Ross, L.O. Silva, D.H. Froula, Phys. Rev. Lett. 105 (2010) 105003. http://dx.doi.org/10.1103/ PhysRevLett.105.105003. [10] J. Osterhoff, A. Popp, Z. Major, B. Marx, T.P. Rowlands-Rees, M. Fuchs, M. Geissler, R. Hörlein, B. Hidding, S. Becker, E.A. Peralta, U. Schramm, F. Grüner, D. Habs, F. Krausz, S.M. Hooker, S. Karsch, Phys. Rev. Lett. 101 (2008) 085002. http://dx.doi.org/10.1103/PhysRevLett.101.085002. [11] M. Litos, E. Adli, W. An, C.I. Clarke, C.E. Clayton, S. Corde, J.P. Delahaye, R.J. England, A.S. Fisher, J. Frederico, S. Gessner, S.Z. Green, M.J. Hogan, C. Joshi, W. Lu, K.A. Marsh, W.B. Mori, P. Muggli, N. Vafaei-Najafabadi, D. Walz, G. White, Z. Wu, V. Yakimenko, G. Yocky, Nature 515 (7525) (2014) 92–95. http://dx.doi.org/10.1038/nature13882. [12] M. Fuchs, R. Weingartner, A. Popp, Z. Major, S. Becker, J. Osterhoff, I. Cortrie, B. Zeitler, R. Horlein, G.D. Tsakiris, U. Schramm, T.P. Rowlands-Rees, S.M. Hooker, D. Habs, F. Krausz, S. Karsch, F. Gruner, Nature Phys. 5 (11) (2009) 826–829. http://dx.doi.org/10.1038/nphys1404.

R.E. Robson et al. / Annals of Physics 356 (2015) 306–319

319

[13] S. Kneip, C. McGuffey, J.L. Martins, S.F. Martins, C. Bellei, V. Chvykov, F. Dollar, R. Fonseca, C. Huntington, G. Kalintchenko, A. Maksimchuk, S.P.D. Mangles, T. Matsuoka, S.R. Nagel, C.A.J. Palmer, J. Schreiber, K.T. Phuoc, A.G.R. Thomas, V. Yanovsky, L.O. Silva, K. Krushelnick, Z. Najmudin, Nature Phys. 6 (12) (2010) 980–983. http://dx.doi.org/10.1038/nphys1789. [14] C. Birdsall, A. Langdon, The Adam Hilger Series on Plasma Physics, McGraw-Hill, 1985. [15] R.E. Robson, R.D. White, Z.L. Petrović, Rev. Modern Phys. 77 (2005) 1303–1320. http://dx.doi.org/10.1103/RevModPhys.77.1303. [16] A. Pukhov, J. Meyer-ter Vehn, Appl. Phys. B 74 (4–5) (2002) 355–361. http://dx.doi.org/10.1007/s003400200795. [17] I. Kostyukov, A. Pukhov, S. Kiselev, Phys. Plasmas (1994-present) 11 (11) (2004) 5256–5264. http://dx.doi.org/10.1063/1. 1799371. [18] W. Lu, C. Huang, M. Zhou, W.B. Mori, T. Katsouleas, Phys. Rev. Lett. 96 (2006) 165002. http://dx.doi.org/10.1103/ PhysRevLett.96.165002. [19] T. Mehrling, C. Benedetti, C.B. Schroeder, J. Osterhoff, Plasma Phys. Control. Fusion 56 (8) (2014) 084012. http://dx.doi.org/ 10.1088/0741-3335/56/8/084012. [20] K. Floettmann, Phys. Rev. ST Accel. Beams 6 (2003) 034202. http://dx.doi.org/10.1103/PhysRevSTAB.6.034202. [21] P. Michel, C.B. Schroeder, B.A. Shadwick, E. Esarey, W.P. Leemans, Phys. Rev. E 74 (2006) 026501. http://dx.doi.org/10. 1103/PhysRevE.74.026501. [22] A. Vlasov, J. Phys. USSR 9 (1945) 25–40. [23] K. Kumar, J. Phys. D: Appl. Phys. 14 (12) (1981) 2199. URL http://stacks.iop.org/0022-3727/14/i=12/a=008. [24] E. Courant, H. Snyder, Ann. Physics 3 (1) (1958) 1–48. http://dx.doi.org/10.1016/0003-4916(58)90012-5. [25] E. Esarey, P. Catravas, W.P. Leemans, AIP Conf. Proc. 569 (1) (2001) 473–486. http://dx.doi.org/10.1063/1.1384377. [26] M. Tzoufras, W. Lu, F.S. Tsung, C. Huang, W.B. Mori, T. Katsouleas, J. Vieira, R.A. Fonseca, L.O. Silva, Physics of Plasmas (1994-present) 16 (5) (2009) http://dx.doi.org/10.1063/1.3118628. [27] R. Robson, Introductory Transport Theory for Charged Particles in Gases, World Scientific, 2006. [28] B.D. Fried, S.D. Conte, The Plasma Dispersion Function, Academic Press, 1961, pp. 1–419. http://dx.doi.org/10.1016/B9781-4832-2929-4.50001-0. URL http://www.sciencedirect.com/science/article/pii/B9781483229294500010.