Compact electron acceleration and bunch compression in THz waveguides

We numerically investigate the acceleration and bunch compression capabilities of 20 mJ, 0.6 THz-centered coherent terahertz pulses in optimized metallic dielectric-loaded cylindrical waveguides. In particular, we theoretically demonstrate the acceleration of 1.6 pC and 16 pC electron bunches from 1 MeV to 10 MeV over an interaction distance of 20mm, the compression of a 1.6 pC 1 MeV bunch from 100 fs to 2 fs (50 times compression) over an interaction distance of about 18mm, and the compression of a 1.6 pC 10 MeV bunch from 100 fs to 1.61 fs (62 times) over an interaction distance of 42 cm. The obtained results show the promise of coherent THz pulses in realizing compact electron acceleration and bunch compression schemes.


Introduction
The desire to realize electron acceleration schemes that can surpass the approximately 50 MeV/m maximum acceleration gradient of conventional radio-frequency (RF) technology has spurred much research into the use of alternative regions of the electromagnetic spectrum. Methods that have been investigated include laser-induced plasma acceleration [1], vacuum acceleration [2][3][4][5] using optical pulses and dielectric-based acceleration. Dielectric-based acceleration is achieved either by an external optical laser source [6][7][8] or by the wakefields of another electron bunch (i.e. dielectric wakefield accelerator) [9][10][11]. With the advent of efficient coherent THz pulse generation techniques [12][13][14][15], forays have also been made into the acceleration of electrons in vacuum [16] and in waveguides [17] by coherent THz pulses. This paper demonstrates the capabilities of waveguides optimized for acceleration and/or compression of relativistic electron bunches by coherent THz pulses. The relativistic fewfemtosecond pico-Coulomb electron bunch achieved in the bunch compression scheme has applications in single-shot few-femtosecond electron diffraction [18]. We choose to study dielectric-loaded cylindrical metallic waveguides for their ease of manufacturing and theoretical evaluation. The THz frequency range is chosen as the operation range because it appears to strike a compromise between the large wavelength and low acceleration gradient (due to breakdown limitations) of RF radiation and the small wavelength but high acceleration gradient of optical radiation. Note that a higher acceleration gradient is more favorable for bunch compression and acceleration, but space-charge effects make it difficult to confine a bunch of substantial charge well within a half-cycle if the wavelength is too small. The absence of plasma in a vacuum-core waveguide scheme precludes problems associated with the inherent instability of laser-plasma interactions. Although using a guiding structure leads to intensity limitations, it also increases acceleration efficiency due to a smaller driving energy required and a larger interaction distance.
The high thermal conductivity and breakdown properties of chemical-vapor-deposited diamond at THz frequencies are well-recognized, and has led to its use in waveguides for wakefield acceleration [11] and other applications involving intense terahertz radiation [19]. For this reason, we use diamond for the dielectric throughout this study and assume a relative dielectric constant of  r = 5.5 [20]. We employ the fundamental transverse-magnetic waveguide mode (TM 01 mode) because every field component in this mode vanishes on axis except for the z-directed electric field, so an electron bunch close to the axis will be accelerated mainly in the forward direction. Figure 1 illustrates an example of concurrent compression and acceleration of an electron bunch in our scheme. We present this example before any technical discussion to give some preliminary intuition of the electrodynamics that ensues when a 1 MeV electron bunch (obtained, for instance, from an RF gun) is injected into a coherent THz pulse propagating in a dielectric-loaded cylindrical metal waveguide. Note that the work pursued here differs from the study presented in [17], which discusses the design of a uniformly-accelerating 100 MeV/m coherent THz pulse-driven waveguide accelerator. Here, we study the acceleration as well as bunch compression capabilities of a  coherent THz pulse of finite duration. Moreover, the presented simulation results for coherent THz pulse-driven acceleration and compression cannot be taken for granted or inferred by scaling the results from studies at optical or RF frequencies, because of the non-negligible impact of space-charge. In Sec. 2 and 3, we furnish a technical discussion of the equations upon which our model rests. In Sec. 4, we demonstrate the acceleration of a 1.6 pC electron bunch from a kinetic energy of 1 MeV to about 10 MeV over an interaction distance of about 20mm, using a 20mJ pulse centered at 0.6 THz in a dielectric-loaded metallic waveguide. The implications of using an arbitrarily distant injection point, as well as the prospects of dielectric breakdown and thermal damage for our optimized design are also analyzed. In Sec. 5, we investigate the acceleration of 16 pC and 160 pC 1 MeV electron bunches. In Sec. 6, we optimize the dielectric-loaded metal waveguide design for simultaneous acceleration and bunch compression, achieving a 50 times (100 fs 1.6 pC electron bunch compressed to 2 fs over an interaction distance of about 18 mm) and 62 times (100 fs to 1.61 fs over an interaction distance of 42 cm) compression for 1 MeV and 10 MeV electron bunches, respectively.

Relativistic electrodynamics in a waveguide and simulation algorithms
This section introduces the equations governing the behavior of an electron bunch in the vacuum-filled core of a waveguide, and discusses our approach in modeling this behavior. The electron bunch is made up of N interacting electrons that may be modeled classically as N point charges propagating according to Newton's second law: the electron experiences as a result of its own radiation. In this study, we neglect wf i F  because the relatively short propagation distances and bunch lengths make the effect of wakefields negligible. For acceleration studies involving long propagation distances, or multiple bunches of substantial charge, wakefields should be taken into consideration by implementing formulas derived in previous studies [21]. We also neglect the radiation reaction force since the employed scheme accelerates the electrons primarily via the z-directed component of the electric field, with minimal transverse wiggling. Consequently, radiation losses are negligible. Electrodynamic studies in which the radiation reaction force plays a significant role have commonly employed the Landau-Lifshitz formula [22] for the force.
where q is the electron's charge and i r  its position.
where 0 ε is the permittivity of free space, i v   the acceleration of particle i, is the retarded time along particle i's trajectory corresponding to time t and observation point r  . Given t and r  , the retarded time i t  solves the implicit equation  is the only non-zero term on the right-hand side of Eq. (1), the equation is simply an ordinary differential equation. With inter-particle interaction described by Eq. (3) and Eq. (4), the right-hand side of Eq. (1) becomes a function of i t as well as t, and the equation is no longer an ordinary differential equation. Note that Eq. (4) considers both the velocity field (near-field) and the radiation field (far-field), which are given by the first and second term respectively. If the effect of the radiation field is insignificant and we assume that each particle always travels at its current velocity during each time step, Eq. (4) can be simplified to a function of only t, making Eq. (1) an ordinary differential equation and reducing the computation of inter-particle forces considerably. The formulas that should replace Eq. (4) are then the space-charge formulas obtained by Lorentz-boosting the Coulomb field of each electron from the electron's rest frame to the lab frame. These formulas are used in particle tracer programs like the General Particle Tracer (GPT) [24]. We chose not to use externally-provided software packages in part to ascertain, by implementing Eq. (4), the significance of non-uniform motion and electron radiation in interparticle interaction. It turns out that for the regime investigated in this paper, the use of the exact formulas in Eq. (4) affects overall acceleration and bunch compression results negligibly, and for computational efficiency one may simply revert to the Lorentz-boosted Coulomb fields in modeling inter-particle interaction.
We solve Eq. (1) using a fifth-order Runge-Kutta algorithm with adaptive step-size [25]. If the exact inter-particle fields in Eq. (4) are used, we adapt the Runge-Kutta algorithm to the problem by maintaining a history of i r  and i p  , i=1,…,N, in a ring buffer. At each time t, cubic spline interpolation is applied to compute the retarded time Eq. (5) needed in Eq. (4). Gaussian-distributions of electrons in 6-dimensional phase space are generated by applying the Box-Muller transformation to the normalized output of the rand() function in C, and computations of variance and covariance (required for emittance calculations) are performed using the corrected two-pass algorithm [26]. Multi-core processing capabilities are implemented using OpenMP.
In this study, we are interested in simulating bunches on the order of pCs and tens of pCs, implying that we deal with 10 7 -10 8 electrons. To speed up the computational process, each particle i = 1,…, N is treated as a macro-particle -with the charge and mass of a large number of electrons -instead of a single electron. We can verify that this approach is a good approximation if the solution converges as the number of macro-particles increases while the total number of electrons is kept constant. We have verified this for all results presented in this paper.

The pulsed TM 01 mode in a dielectric-loaded metallic waveguide
We are interested in obtaining an analytical expression that models a coherent THz pulse in the waveguide. This involves integration over the continuous-wave (CW) solutions of the waveguide. The method we use to obtain these solutions is very similar to that detailed in [27] for the optical Bragg fiber, so we only give an overview of the method here. For a general multilayer cylindrical waveguide, the continuous-wave solutions are obtained by solving the Helmholtz equation in cylindrical coordinates [23]: where k  c = 2,  being angular frequency, c the speed of light in vacuum and  the vacuum wavelength. E z CW  e (r) . exp(i(t -z ± l)) and H z CW  h (r) . exp(i(t -z ± l)) are the complex CW z-directed electric and magnetic fields respectively,  is the propagation constant, r the radial coordinate, the azimuthal coordinate, z the direction of propagation along the waveguide, and l a non-negative integer that determines the order of azimuthal variation. According to Eq. (6), a general solution for  e in layer i of an n-layer cylindrical waveguide (the core counts as layer 1) is e; e; e; 1 where r 0  0, r n ∞, and r i for 0 < i < n is the radial position of the boundary between layers i and i+1. J l and Y l are Bessel functions of the first and second kind respectively, A e;i and B e;i are constant complex coefficients within each layer and h i  r;i () r;i ()k    ) 1/2 ,  r;i and  r;i being the dispersive relative permittivity and permeability respectively of the dielectric in layer i. The general solution for  h;i is identical in form to Eq. (7) except that "e" should be replaced by "h" in all subscripts. In the core, it is usually expedient to express Eq. (7) me inverse eal-valued electric of ) with an l = 0 and nents exist. on bunch concentrated at the waveguide axis will experience forces mainly along the direction of propagation. This facilitates longitudinal compression and acceleration of the bunch without significant transversal wiggling, which is undesirable since it tends to increase radiative losses.
To excite the TM 01 mode of the cylindrical waveguide, it would be necessary to apply a radially-polarized (preferably TM 01 ) beam to the waveguide. Studies on coupling linearlypolarized THz pulses into cylindrical metal waveguides show that the dominant modes excited are the TE 11 , TE 12 and TM 11 modes [28], so a linearly-polarized incoming beam is unlikely to serve our purpose. Although THz pulses generated by optical rectification are typically linearly-polarized, the direct generation of radially-polarized THz pulses has been demonstrated [29,30]. Alternatively, a scheme to convert linearly-polarized THz pulses into radially-polarized pulses may be adopted [31].
Equation (8) provides a rigorous way to compute the electromagnetic field at any point in space and time required for an electrodynamic simulation. However, performing a summation over a large number of frequency components at every time step for every macro-particle is computationally expensive. To obtain an analytical approximation for more efficient numerical simulation, notice that in the vacuum-filled core, the CW TM 01 mode is of the form: where q i    r;i () r;i ()k  ) 1/2 is the radial wavevector and   is the vacuum impedance. I 0 and I 1 are the modified Bessel functions of the first kind, of order 0 and 1 respectively. We need to make three more assumptions in the remainder of the formulation: Firstly, variations in propagation constant κ across the spectrum are small enough that their effects on magnitude can be ignored. Secondly, variations in κ are negligible above the second order. Thirdly, the imaginary part of κ(ω) is negligible beyond its 0 th order term, and the quadrature term produced by this imaginary part in E r does not contribute significantly to the field. Hence, Taylor-expanding κ() about central angular frequency  0 we have κ( 0 + ∆) κ 0 -i + κ 1 ∆κ 2 (∆  , where κ i denotes the real part of the i th derivative of κ() with respect to at  0 .  to be physically valid and represents field attenuation per unit distance.
To obtain the approximate analytical field solution, the rightmost expressions of Eq. (9) should be inserted into Eq. (8). Assuming a transform-limited Gaussian pulse at z = 0, we have E z (z = 0, t) ~ exp(-(t/T 0 ) 2 /2)exp(i 0 t), where  0  k 0 c is the central frequency and T 0 is the half-width at 1/e intensity, related to the full-width-at-half-maximum (FWHM) intensity T FWHM as T FWHM = 2(ln2) 1/2 T 0    This is related to the spectral FWHM intensity width ∆ FWHM as ∆ FWHM = 4ln2/T FWHM . Finally, we have where A 0 is an arbitrary complex constant and its role is replaced in the second line of Eq. (10) by |E z0 |, which represents the amplitude of the z-directed field at t = 0 and z = z i = z s , with z i being the initial position of the pulse peak. The precise relationship between |E z0 | and the total pulse energy is complicated and must be obtained by integrating over the Poynting vector in both core and cladding regions. q 1,0 is q 1 evaluated at  0 . z s is the position of the start of the waveguide, where pulse attenuation begins, and before which Eq. (10) does not apply. Note that setting z s ≠ 0 implies that some special pulse, not transform-limited, is being coupled into the waveguide. We set z s = 0 for all simulations in this paper. The carrier phase   is given by where   is a real phase constant. The corresponding E z , E r and H  fields are approximated as   .
where k 0   0 /c. Essentially, Eq. (10) and Eq. (12) furnish an approximate analytical description of a TM 01 pulse moving with an approximate phase velocity and group velocity of v ph =  0 /κ 0 and v g = 1/κ 1 respectively in the vacuum core of a cylindrical waveguide. If z s = 0, the pulse at the start of the waveguide (z = z s = 0) is a transform-limited pulse with a peak longitudinal electrical amplitude of |E z0 |. The primary reason for introducing z i in our formulas is to control when the pulse arrives at the start of the waveguide without having to compromise the intuitive convention of having t = 0 as the initial time (when the simulation begins and the initial electron bunch starts evolving according to Eq. (1)).

Optimization procedure and acceleration results
In this section, we optimize the dielectric-loaded metal waveguide for electron bunch acceleration and perform a rudimentary thermal damage and dielectric breakdown analysis to verify the realism of the scheme. We numerically demonstrate the acceleration of a 1.6 pC electron bunch from a kinetic energy of 1 MeV to one of 10 MeV, using a 20 mJ 10-cycle pulse centered at 0.6 THz. Note that for a 10-cycle pulse, ∆ FWHM /  = 4ln2/(  T FWHM ) = 4ln2/(As will be seen in the results, some longitudinal compression is also inadvertently achieved in the process. Optimizing the dielectric-loaded metallic waveguide for bunch acceleration involves adjusting a large number of parameters, including operating frequency, choice of waveguide mode, waveguide dimensions, THz pulse energy and pulse duration, the type of dielectric, the type of external conductor and initial electron bunch properties. To make this optimization tractable, we fix all parameters in advance based on the available technology except for three degrees of freedom: i) the carrier-envelope phase   , ii) the initial position of the pulse z i (with initial position of electron fixed at z = 0), and iii) the radius of the vacuum core r 1 . In particular, we fix the phase velocity at v ph c and the center frequency at f 0 = 0.6 THz, which limits the dielectric thickness d to specific values depending on r 1 . However, because acceleration results can be very sensitive to small variations in the value of v ph , we take the liberty of treating v ph as an optimization parameter (but ensuring that v ph c) after using v ph c to determine properties of the TM 01 waveguide mode. Therefore, four degrees of freedom are ultimately considered. In practice, after the waveguide has been fabricated according to the optimal specifications, the operating frequency should be perturbed to vary the phase velocity until maximum electron acceleration is achieved. As long as the perturbation is small, the waveguide properties should be very close to those determined for v ph c and f 0 = 0.6 THz. The electron acceleration process is much more sensitive to small variations in v ph than to small variations in any other parameter caused by perturbing the operating frequency alone.  Figure 3(a) shows a color map of the operation frequency as a function of r 1 and d. As noted before, we define the operation frequency as the frequency of the TM 01 mode in the waveguide corresponding to v ph c. Fig. 3(b) shows a color map of the final electron kinetic energy of a single electron of initial kinetic energy 1 MeV, optimized over    z i and v ph (ensuring that v ph c), as a function of r 1 and d. We see that greater electron acceleration is generally achieved at higher operation frequencies. However, choosing a very small wavelength makes it challenging to accelerate a large number of electrons due to smaller waveguide dimensions. As pointed out previously, the emergence of promising techniques to generate radiation in the vicinity of 0.6 THz [13] encourages us to make that choice of frequency, which has been marked out by the black contour line in Fig. 3(a) Fig. 3(b), and the optimized final kinetic energy, read along that line, is reproduced in Fig. 3(c), where an optimal choice of d = 32 m, corresponding to a vacuum core radius of r 1 = 380 m, is evident. In Fig. 3(d), we plot the dispersion curves corresponding to the waveguide with d = 32 m, r 1 = 380 m, to show that at the operating frequency, the TM 01 dispersion curve of our waveguide design is sufficiently linear within the 4.41% intensity FWHM spectral bandwidth. Hence, the electromagnetic fields are well approximated using Eq. (10) and Eq. (12).
The parameters of the final waveguide design are d = 32 m, r 1 = 380 m, v ph = 0.99c, v g = 0.7c, = 5.21/m, T FWHM = 16.7 ps,  2 = 4.54×10 -22 s 2 /m. The 20 mJ pulse yields a |E z0 | of about 0.9 GV/m. The initial parameters of the 1.6 pC, 1 MeV electron bunch with which we will demonstrate the acceleration are  x =  y =  z = 30 m (a 100 fs bunch) ,  x =  y =  z = 0.006, where  x , for instance, denotes the standard deviation of  x . Producing a 1.6 pC, 100 fs electron bunch would be a challenge for typical RF guns, but strides are being made to realize a photocathode RF gun capable of delivering the bunch we have assumed as our input [32]. Although a thorough examination of how performance is impacted by variations in the initial electron bunch properties is beyond the scope of this paper, we expect the results to deteriorate with a larger initial energy spread. 10000 macro-particles, Gaussian-distributed in every dimension of phase space, were employed in the simulation.  Figure 4 shows the evolution of bunch parameters as a function of mean particle position. We see from Fig. 4(a)

Injection point considerations
In our analysis, we have assumed the freedom to inject the electron bunch into any point of the electromagnetic field. According to our computations, the optimum injection point for the electron bunch is a point within the pulse (albeit in its tail). This may be challenging to realize if both the electron bunch and the electromagnetic pulse enter the waveguide from vacuum. The objective of this section is to consider injection of the electron bunch at a point with negligible electric field values and assess the amount by which our predictions would change. The optimum THz waveguide for this case is a waveguide with r 1 = 338 m and d = 33 m.
In addition, v ph = 0.981c,   , k 0 z i = 137.73. We ensure that the electric field's amplitude at the injection point is negligible by making the amplitude 7.4×10 -15 |E z0 |. The evolution of the electron bunch is shown in Fig. 5, where we observe a final kinetic energy of 8.4 MeV (instead of the 9 MeV observed before). The smaller energy gain in this case is partly due to the dispersion and attenuation that the pulse suffers from before the injected bunch begins interacting with the pulse. A final energy close to what is predicted in the previous section should therefore be achievable if the electron bunch and THz pulse can interact before the pulse has travelled too far along the waveguide.

Thermal damage and dielectric breakdown considerations
In this section, we assess the feasibility of the scheme from Sec. 4.1 in terms of its thermal damage and dielectric breakdown prospects. One concern is that the high energy injected into the waveguide and consequent energy dissipation would raise temperature of the copper coating beyond its melting point. Another concern is dielectric breakdown due to the high electric field values in the dielectric.
The energy dG transferred to a differential segment of copper at position z (z = 0 being the start of the waveguide) is related to the associated temperature rise ∆  as From the values in Fig. 6(a), the fact that the melting point of copper is 1084 0 C and also that we have even neglected the conductivity of copper, we can conclude that the metal coating in the designed waveguide withstands the passage of the pulse without melting. Figure 6(b) shows a typical profile of the electromagnetic amplitude of a mode in the transverse direction of the waveguide. The breakdown electric field for diamond has been reported as 10-20 MV/cm, depending on impurities. Reading off the plot we note that the maximum value of the electric field in the dielectric region is about 8 MV/cm. This is close to

Accelerat
In this secti that it is fea bunches as lar All other bun from Sec. 4 is counter-balanced by larger transverse inter-particle spacing. Figure 7(d) shows that due to the larger amount of space charge, the 16 pC expands rather rapidly compared to the 1.6 pCbunch after the pulse has slipped behind the bunch, so a 16 pC-bunch accelerated via this scheme is likely to be useful for a shorter duration after being fully accelerated.

Concurrent phase-limited compression and acceleration of 1.6pC-bunches
In this section, we optimize our waveguide design for simultaneous acceleration and bunch compression. We demonstrate phase-limited (longitudinal) bunch compression of 50 and 62 times for electron bunches of initial kinetic energy 1 MeV and 10 MeV respectively. By "phase-limited" we mean that the maximum compression results do not change substantially when space charge is removed from the simulations.
As in previous sections, we use a 20 mJ, 0.6 THz-centered pulse. For each case (the 1 MeV case and the 10 MeV case), the waveguide and injection conditions are optimized exactly as described in Sec. 4.1, except that in addition to   , z i , r 1 and v ph , we also optimize over pulse duration T FWHM (keeping total energy constant at 20 mJ), for a total of five optimization parameters. The initial conditions of the electron bunch, unless otherwise specified, are the same as those in Sec. 4.1.
To optimize for simultaneous acceleration and compression, the figure-of-merit found to be most useful is the ratio of energy to bunch-length of the electron bunch. Unlike in Sec. 4.1, where we optimized using a single particle, here we optimized using 100 macro-particles and included the effects of space charge. The optimized results are then verified with simulations that use10000 macro-particles.
For the 1 MeV case, our optimized parameters are   , k 0 z i = 13.3, r 1 = 447 m FWHM = 13.1 ps (7.86 cycles). The evolution of the electron bunch parameters under these optimal conditions are presented in Figs. 8(a)-8(c), where we observe a small net acceleration and a phase-limited compression of the electron bunch from 100 fs (30 m) to about 2 fs over an interaction distance of about 18 mm. Note that there is a limited time window during which the electron bunch remains maximally compressed. Conceptually, this is unavoidable due to the presence of space charge which causes the bunch to expand after the bunch has slipped from the THz pulse.
For the 10 MeV case, our optimized parameters are   , k 0 z i = 206r 1 = 597m FWHM = 170.5 ps (102.3-cycle). The evolution of the electron bunch parameters under these optimal conditions are presented in Figs. 8(d)-8(f), where we observe a phaselimited compression of the electron bunch from 100 fs to 1.61 fs over an interaction distance of 42 cm. Although the bunch is compressed by a slightly larger factor than in the 1 MeV case, the much larger interaction distance suggests that the superior strategy to obtain a high energy, compressed bunch is to compress it before acceleration.

Conclusion
The quest to realize an efficient, practical compact accelerator for electron bunches of substantial charge will likely involve a tradeoff between the large wavelengths but low acceleration gradient of RF accelerators, and the high acceleration gradient but small wavelengths available at optical frequencies. The trade-off between acceleration gradient and wavelength, together with the emergence of efficient methods to generate coherent pulses at THz frequencies, make electron acceleration at THz frequencies a promising candidate for the substantial acceleration and compression of pico-Coulomb electron bunches. In this paper, we numerically demonstrated the acceleration of a 1.6 pC electron bunch from a kinetic energy of 1 MeV to one of 10 MeV over an interaction distance of about 20 mm, using a 20 mJ pulse centered at 0.6 THz in a dielectric-loaded metallic waveguide. We have also analyzed the implications of using an arbitrarily distant injection point, as well as the prospects of dielectric breakdown and thermal damage for our optimized design. In addition, observing that prohibitively optimized the bunch compre over an inter interaction di respectively. THz, and enc path to compa d 1 e J, e d d e r n bunches, eteriorates inally, we ation and ed to 2 fs s over an on bunch red at 0.6 ation as a 4192. We nowledges