Three-Dimensional Analysis of Wakefields Generated by Flat Electron Beams in Planar Dielectric-Loaded Structures

An electron bunch passing through dielectric-lined waveguide generates $\check{C}$erenkov radiation that can result in high-peak axial electric field suitable for acceleration of a subsequent bunch. Axial field beyond Gigavolt-per-meter are attainable in structures with sub-mm sizes depending on the achievement of suitable electron bunch parameters. A promising configuration consists of using planar dielectric structure driven by flat electron bunches. In this paper we present a three-dimensional analysis of wakefields produced by flat beams in planar dielectric structures thereby extending the work of Reference [A. Tremaine, J. Rosenzweig, and P. Schoessow, Phys. Rev. E 56, No. 6, 7204 (1997)] on the topic. We especially provide closed-form expressions for the normal frequencies and field amplitudes of the excited modes and benchmark these analytical results with finite-difference time-domain particle-in-cell numerical simulations. Finally, we implement a semi-analytical algorithm into a popular particle tracking program thereby enabling start-to-end high-fidelity modeling of linear accelerators based on dielectric-lined planar waveguides.


I. INTRODUCTION
Next generation multi-TeV high-energy-physics lepton accelerators are likely to be based on nonconventional acceleration techniques given the limitation of radio-frequency (rf) normal-conducting [1] and superconducting [2] structures. Non-conventional approaches based on the laser plasma wakefield accelerator have recently demonstrated average energy gradients one order of magnitude higher than those possible with state-ofthe art conventional structures [3]. The applicability of laser-driven techniques to high-energy accelerators is currently limited as attaining luminosity values similar to those desired at the International Linear Collider would demand a laser with power approximately four orders of magnitude larger than the most powerful lasers currently available [4]. Another class of non-conventional accelerating techniques includes beam-driven methods which rely on using wakefields produced by high charge drive bunches traversing a high-impedance structure to accelerate subsequent witness bunches [7]. Such an approach has the advantage of circumventing the use of an external power source and can therefore operate at mm and submm wavelengths. Structures capable of supporting wakefield generation include plasmas [5], and dielectric-loaded waveguides [10]. The possible use of plasma-wakefield accelerators (PWFAs) as the backbone of a a multi-TeV electron-positron linear collider is limited by plasma ions motion due to the intense electromagnetic field of the bunch [6]. Dielectric wakefield accelerators (DWFAs) are not prone to similar limitations.
In this paper we concentrate on the colinear DWFA. In such a configuration a highly-charged drive bunch propagates through a dielectric-lined waveguide (DLW) and excites an electromagnetic wake [8][9][10]. A delayed witness bunch moving on the same path as the drive bunch can experience an accelerating field.
Recent experiments [13] confirm that DLW can support accelerating fields in excess of a GV/m thereby making DWFA a plausible candidate for the next generation high-energy-physics linear accelerators [14] or compact short-wavelength free-electron lasers [15].
In cylindrically-symmetric DLW structures, the electric field amplitude of the wakefield is approximately inverse proportional to the aperture radius. Given the linear charge-scaling of the field, peak electric field amplitude of the order of Gigavolt-per-Meter can be obtained by either using high-charge drive beams (> 100 nC) in mm-sized DLW structures or by focusing sub-nC bunches in micron-sized DLW structures [9,11,12].
To date, cylindrically-symmetric DLW's have been extensively studied both theoretically [9,[16][17][18] and experimentally [19][20][21][22]. Since the angular divergence of the beam also increases during focusing it is hard to maintain a low transverse size of the beam over a long propagation distance. Therefore, the design of the DLW structures must compromise between a small transverse size to maximize the intensity of the wakefield, and a longer interaction length to maximize the energy gain of the test charge. An appealing solution consists of using flat electron beams (σ x ≫ σ y ) passing through slab-geometry DLWs [11,[23][24][25]. This possibility has become more attractive since the recent advances toward generating flat beams directly in photoinjectors [26,27].
In the following Sections we present an analysis of the generation of wakefields in rectangular DLW structures excited by drive beams with arbitrary three-dimensional charge distributions. The main goal of this paper is to extend the formalism introduced in Ref. [11] by including all types of modes excited in a planar DLW. This paper is pedagogical in the sense that it is about the method of derivation rather than the results, which have been derived in previous papers, most extensively in Ref. [24]. However, this paper provides close-form formulae for the eigenfrequencies and electromagnetic fields excited in planar DLWs. These analytical results are benchmarked against three-dimensional finite-difference time-domain (FDTD) simulations. Finally, the model is implemented in a popular particle-in-cell (PIC) beam dynamics program. The latter provides a fast and highfidelity model enabling start-to-end simulation of DLWbased linear accelerators.

II. WAKEFIELD GENERATION
The geometry of the problem analyzed in this paper is depicted in Fig. 1. Transient effects resulting from the injection of the electron bunch in the structure are not included (the structure is assumed to be infinitely long) and the drive bunch is taken to be ultra-relativistic with its Lorentz factor γ much greater than unity. The building block of a real drive bunch is assumed to consist of a linear charge distribution oriented along the x-axis which moves in the z-direction with velocity v and has offset y 0 in the vertical direction. To simplify, it is also assumed that the linear charge distribution is symmetric in the x-coordinate and it vanishes at the ends of the dielectric structure. Under these assumptions the charge distribution can be written as a Fourier series in the horizontal direction [11] ρ(x, y, z) = m λ m cos(k x,m x)δ(y − y 0 )δ(z − vt), (1) where λ m is a constant and each term is indexed by the integer m = 0, 1, · · · defined such that k x,m ≡ (2m + 1) π Lx . In the charge-free vacuum region the electromagnetic wakefield satisfies the wave equation: In this paper we consider only the propagating modes with harmonic longitudinal and temporal dependencies of the form E(x, y, z, t) = E(x, y)e i(ωt−kz z) [11,25] where k z = k.ẑ. In addition to the propagating-mode solutions, Eq. 2 when subject to the boundary conditions also admits evanescent modes with imaginary eigenfrequencies ω [25]. Since these evanescent modes are present over short distances behind the drive bunch, they do not contribute to the beam dynamics of a subsequent witness bunch. These short-range fields are consequently ignored in the remaining of this paper.
In the ultra-relativistic regime the synchronism condition is achieved because the wakefield phase velocity v ϕ ≡ ω kz = v ≈ c is very close to the velocity of the witness beam. Since the transverse wave vector k ⊥ ≡ k 2 − k 2 z ≃ k/γ = 0 the wave equation 2 can be written only in terms of the partial derivatives with respect to the transverse coordinates: The symmetry of the drive-beam charge distribution and the boundary conditions at x = ±L x /2, determine (up to constants) the analytical expression of the fields. In the vacuum region the electric field components are simple combinations of trigonometric and hyperbolic functions: The axial field E z associated with this set of solution is symmetric with respect to the horizontal axis [i.e. E z (x, −y, z) = E z (x, y, z)] we henceforth refer this set to as "monopole" modes. Similarly, a set of fields with an antisymmetric axial field [E z (x, −y, z) = −E z (x, y, z)] is obtained by substitution of sinh(k y y) with cosh(k y y) and vice versa. This latter set is termed as "dipole" modes in the remaining of this paper.
Inside the dielectric, the transverse wave number cannot be neglected where ǫ r is the relative electric permittivity of the dielectric medium. The boundary conditions at y = ±b and at x = ±L x /2 determine the trigonometric form of the field expressions. Inside the dielectric region, there is no distinction between monopole and dipole modes and the electric field components are given by The corresponding expressions for the magnetic field in the vacuum and dielectric regions can be easily obtained following a similar prescription.
Inspection of the z-component of the fields indicates that the normal modes cannot be categorized in the usual transverse electric or transverse magnetic sets. The reason is that the separation surface between the dielectric and vacuum regions is in the x − z plane unlike the case of the uniformly-filled waveguides where this surface is in the transverse x − y plane. Therefore, it is natural to categorize the modes depending on the orientation of the fields with respect to the dielectric surface. Following the definition introduced in Ref. [29] we classify the mode as Longitudinal Section Magnetic (LSM) and Longitudinal Section Electric (LSE) modes corresponding respectively to the case when the magnetic and electric field component perpendicular to the dielectric surface vanishes. In the configuration shown in Fig. 1, LSE and LSM modes correspond respectively to E y = 0 and H y = 0 at the vacuum-dielectric interface.

A. Longitudinal Section Magnetic (LSM) modes
To obtain the normal mode frequencies it is convenient to use the hertzian potential vector method [29]. In the source-free region, magnetic and electric fields can be expressed in terms of the vector potential Π: where Π is the solution of the homogeneous Helmholtz equation in the vacuum region. To determine the frequencies and the fields for the modes when H y = 0 (LSM) it is convenient to choose Π ≡ Π e = ψ e (x, y)e i(ωt−kz z)ŷ where the the scalar function ψ e (x, y) must be a solution of (∇ 2 ⊥ + k 2 ⊥ )ψ e = 0. From Eqs. 7 and 8 all fields can be expressed in terms of the unknown function ψ e (x, y): The expression for ψ e that satisfies Eqs. 9, 5, and 7 is given by: where A and B are constants. The boundary conditions completely determine the eigenfrequencies of the LSM modes. Since H x (∝ ǫψ e ) and E z (∝ ∂ψe ∂y ) are continuous at y = a the constants A and B can be eliminated giving the dispersion equation Therefore, for each discrete value of k x,m there is an infinite set of discrete k y,n values where n is an integer. So, the eigenfrequencies are indexed by the integer couple (m, n) and verify It is worthwhile noting that the boundary conditions do not completely determine either the function ψ e or the fields. The reason being that, up to this point, the field source terms (ρ, j) were not taken into account although their symmetries were invoked. Still, it is straightforward to show that all fields (and also ψ e ) associated to a certain mode depend on a common normalization constant, the amplitude E 0;m,n , which remains to be determined. The expressions of the fields are given by

B. Longitudinal Section Electric (LSE) modes
The case of the LSE modes (E y = 0) can be treated in the same way as the LSM modes. The hertzian vector electric potential is replaced by a vector magnetic potential Π h which is related to the fields: As in the previous case, it is convenient to factor out the t and z-dependencies of the Π h and to define a scalar function ψ h (x, y): It is straightforward to derive the dispersion equation for the LSE modes and the expressions for the electromagnetic-field components:

IV. WAKEFIELD AMPLITUDES
Although we have obtained the general expression for the electromagnetic field, the constants E 0;m,n still remained to be determined. An often used method to find the constants E 0;m,n in Eqns. 13 and 17, is to evaluate the Green function of the wave equation with sources included. The latter is usually a cumbersome procedure and we instead chose to determine E 0;m,n 's based on energy balance considerations. The electromagnetic energy stored in the DLW equals the mechanical work performed on the drive bunch. For a linear wakefield, the decelerating field E d acting on a drive point-charge is half of the wakefield amplitude, i.e. E d = E0;m,n 2 as a consequence of the fundamental wakefield theorem [30]. . Suppose the drive charge moves over an infinitely short distance δz. The work performed by the wakefield on the drive charge should equal the energy stored in the field m,n δ(y − y 0 ) E z;m,n (z = vt) 2 dxdyδz where the integration extends over the transverse plane. It is straightforward to evaluate the full expressions of the wakefield amplitudes in terms of the wave numbers k x,m and k y,n from Eqns. 18, 13 and 17: The field amplitudes for the dipole modes can be obtained from the previous equations by substituting sinh(k x,m u) for cosh(k x,m u) where u takes on the val-ues of a and y 0 . As expected, the wakefield amplitude scales linearly with the drive bunch charge and inverse proportionally with the transverse size of the structure. To this point this model is fully three-dimensional and the only limitation stems from the assumed symmetry of the transverse drive charge distribution with respect to the y axis.

V. COMPARISON WITH THREE-DIMENSIONAL FDTD SIMULATIONS
To evaluate the wakefield associated to a drive bunch, an integration over the full three-dimensional continuous charge distribution must be performed. In practice the integration is replaced by numerical summations of discrete charge distributions similar to those described by Eq. 1. This process is similar to the charge discretization procedure used in standard PIC algorithms [32]. An important feature of this model is that the integration over x-direction is already performed through the Fourier expansion of the charge distribution. For the charge distributions considered in this paper, only a few Fourier terms (< 10) are sufficient to obtain an accurate representation of the distribution along the x direction. This is significantly less than the number of grid points in x-direction needed by most PIC codes to evaluate the 3D-collective effects and external fields.
The integration over the vertical direction is straightforward and in most cases, when the drive charge distribution is also symmetric with respect to the horizontal axis, the contribution of the dipole modes cancels out.
Since the phase velocity of the wakefield is the same as the velocity of the drive beam, causality principle requires that the wakefield vanishes ahead of the drive charge. Therefore, the integration over the longitudinal direction extends only from the observation point to the actual drive charge position. The wakefield assumes the form where W stands for any of the electromagnetic field components and W m,n are the corresponding field component associated to the LSM m,n and LSE m,n modes given respectively by Eq. 13 and Eq. 17. The summation is performed over all excited modes.
The results obtained from Eq. 20 are benchmarked against simulations performed with vorpal a conformal FDTD (CFDTD) PIC electromagnetic solver [28]. vorpal is a parallel, object-oriented framework for three dimensional relativistic electrostatic and electromagnetic plasma simulation. The DLW model implemented in vorpal is fully three dimensional; see Fig. 2. The model consists of the rectangular DLW surrounded by perfectly-conducting boundaries (PCBs), The lower and upper z planes are terminated by perfectly matched layers (PMLs) that significantly suppresses artificial reflections of incident radiation [31]. A particle source located on the surface of the lower z plane produced macroparticles uniformly distributed in the transverse plane (x, y) and following a Gaussian longitudinal distribution with total charge Q described by where ζ ≡ z − vt, w x (w y ) is the full width transverse horizontal (vertical) beam size, and σ z the longitudinal root-mean-square (rms) length. The function H(...) is the Heaviside function. The longitudinal Gaussian distribution is truncated at ±3σ z . Finally, a "particle sink" at the upper z plane allows macroparticles to exit the computational domain without being scattered or creating other source of radiation. In order to precisely benchmark our model, a DLW which supports both LSM and LSE modes is chosen. The parameters of the structure and driving bunch are gathered in Tab. I and the frequency and the amplitude associated to the first few modes appear in Fig. 3. The drive-bunch energy is set to E = 1.0 GeV consistent with the ultra-relativistic approximation used in the analytical model. The longitudinal component of the electric field simulated with vorpal is shown in Fig. 4 as a twodimensional projection in the y − z plane.  Fig. 6 (right)] transverse coordinate at a given axial location z = 18.6 mm corresponding to the minimum E z shown in Fig. 5 (i.e. maximum accelerating field). All plots show a decent agreement (relative discrepancy < 25 %) between our model and vorpal simulations. The notable disagreement observed in Fig. 5 for the transverse fields in the vicinity of the driving charge (i.e. z ≃ 27.5 mm) is rooted in the absence of velocity fields in our model (only radiation field contributes to the wakefield) while the vorpal simulations include both velocity and radiation fields. In fact given the bunch charge and duration, the amplitude of the velocity field can be evaluated by convolving the charge distribution with the electric field generated by an ultra- relativistic particle E(r, ζ) = q/(2πǫ 0 )(r ⊥ /r 2 ⊥ )δ(ζ) where e and ǫ 0 are respectively the electronic charge and vacuum permittivity, and r ⊥ ≡ (x, y) and ζ are respectively the transverse and longitudinal coordinates of the observation point referenced to the particle's location. Such a convolution with Eq. 21 yields the transverse electric field components at (x, y, ζ) = (w x , w y , 0) to be ∼ 1 MV/m for the bunch parameter listed in Tab. I. The latter value is in agreement with the observed difference between the vorpal and theoretical models; see E x and E y components in Fig. 6.

VI. THE TWO-DIMENSIONAL LIMIT
An interesting limiting case occurs when L x ≫ L y . In this "two-dimensional limit", the dependence of the fields on the horizontal coordinate x is weak and completely vanishes when L x −→ ∞ (so that k x ≃ 0 and m = 0). In such a case the field amplitudes associated to the LSM and LSE modes are respectively , and E LSE 0;0,n ≃ 0, where Λ is charge per unit length in the horizontal direction. The LSE modes are suppressed and the latter equation is in agreement with the results of Ref. [11]. This limit case can be practically reached by using structures with large aspect ratios (L x ≫ L y ) driven by flat electron beams (tailored such that σ x ≫ σ y ). Flat beams can be produced in photoinjectors by using a round-to-flat beam transformation [33,34]. In such a scheme, a beam with large angular-momentum is produced in a photoinjector [35]. Upon removal of the angular momentum by applying a torque with a set of skew quadrupole, the beam has its transverse emittance repartitioned with a tunable transverse emittances ratio [36]. Flat beams with transverse sizes of ≈ 100 µm and aspect ratio of ∼ 20 have been produced, even at a relatively-low energy of 15 MeV [27]. Based on the previous experience in producing flat beams and preliminary simulation of the Advanced Superconducting Test Accelerator (ASTA) currently in construction at Fermilab [37]. It is reasonable to consider a 3-nC flat beam generated from the ASTA photoinjector to have parameters tabulated in Tab. II when accelerated to 1 GeV. Considering a structure with a = 100 µm, b = 300 µm and ǫ = 4.0 would yield a maximum axial wakefield amplitude of ∼ 300 MV/m; see Fig. 7. The fundamental LSM frequency is f 0,0 = 193 GHz and the nearest LSE mode amplitude is ∼ 136 times lower. In Fig. 7 the wakefield is computed with the asymptotic limit provided in Eq. 22 and is in excellent agreement with the FDTD simulations.

VII. IMPLEMENTATION IN A PARTICLE-IN-CELL BEAM DYNAMICS PROGRAM
Although FDTD simulations provide important insight that can aid the design and optimization of the DLW geometry, their use to optimize a whole linear accelerator would be time and CPU prohibitive. Therefore it is worthwhile to include a semi-analytical version of the model developed in this paper in a well-established beam dynamics program impact-t [38]. Our main motivation toward this choice stems from the availability of a wide range of beamline elements models. In addition, impactt takes into account space-charge forces using a threedimensional electrostatic solver. The algorithm consists in solving Poisson's equation in the bunch's rest frame and Lorentz-boosting the computed electrostatic fields in the laboratory frame.
In a typical particle-tracking PIC program, like impact-t, an electron bunch is described by a set of "macroparticles" arranged to mimic the bunch phase space distribution. Each macroparticle represents a large number of electrons (typically 10 3 in our simulations). To evaluate the electrostatic fields in the bunch's rest frame the macroparticles are deposited on the cells of a threedimensional grid. As a result of this charge deposition algorithm, the initial charge distribution is approximated by a set of point charges, located at the nodes of the three-dimensional grid.
The line charge density corresponding to certain y and z-coordinates is ρ(x) = x0 q(x 0 )δ(x − x 0 ) and provided this distribution is symmetric the Fourier coefficients λ m from Eq. 1 are given by: In practice the linear charge distribution may not be "exactly" symmetric. In this case the charge at an arbitrary grid point x 0 is set at the average of the values given by the charge deposition algorithm at positions x 0 and −x 0 . The general expressions of the wakefields at a given position have the following form where i ≡ x, y, z, n is the mode index, f i 's are geometry dependent functions of the type described in Eq. 5 and q(x 0 , y 0 , z 0 ) cos(k x,m x 0 ) cosh(k x,m y 0 ). (25) In the latter equation the substitution cosh ↔ sinh should be made when considering dipole modes.
The simulated electromagnetic field components are compared with vorpal simulations in Fig. 8 for the same case as presented in Fig. 5. In the Fig. 8 the impactt simulation are performed with and without activating the space-charge algorithm. When accounting for spacecharge forces, impact-t is in very good agreement with vorpal.
The convolution over the z-summation in Eq. 24 can be be performed with a numerical Fast Fourier Transformation (FFT). Therefore the total number of operations needed to evaluate the wakefields is ∝ N modes N 2 x N 2 y N z log N z . The total number of modes is twice the product between the modes allowed for each of the transverse wave numbers: N modes = 2N kx N ky . The factor of 2 comes from the inclusion of the dipole modes along with the vertically symmetric monopole modes. For the data generated in Fig. 8, the impact-t simulations are more than two orders of magnitude faster than the vorpal ones.
The altered version of impact-t was used to explore a possible DLW experiment at the ASTA facility as a 1-GeV electron beam is injected in a DLW. The beam and structure have the parameters displayed in Tab. II . For these simulations, the electron beam is an idealized cold beam without energy spread nor divergence. Figure 9 compares the longitudinal phase spaces simulated with vorpal (top row) and impact-t (bottom row) at two axial locations. The agreement on the phase space structure developing as this non-optimized beam propagates in the DLW is excellent. The mean and rms energies evolution for the witness and drive bunches reported in Fig. 10 confirms the quantitative agreement between the FDTD and semi-analytical models. Both models give mean and rms energies in agreement within 5 %.

VIII. CONCLUSIONS
In this paper we presented a three-dimensional model to evaluate wakefields in slab-symmetric dielectric-lined waveguides. The only limitation of this model stems from the assumed charge distribution symmetry with respect to the vertical axis [ρ(x) = ρ(−x)]. The model was successfully validated against three-dimensional FDTD PIC simulations performed with vorpal and was implemented in the popular beam dynamics tracking program impact-t. The added capability to impact-t enables start-to-end simulation of linear accelerators based on DLW accelerating structures. Furthermore, because impact-t includes a space charge algorithm, the upgraded version provides a valuable tool for investigating the performances of DLW acceleration when the dynamics of either, or both, of the drive and witness bunches is significantly impacted by space charge effects. The observed good agreement between the developed algorithm and simulations performed with the vorpal FDTD PIC program demonstrates that our model strikes an appropriate balance between, efficiency, accuracy, and simplicity.