Terahertz-optical intensity grating for creating high-charge, attosecond electron bunches

Ultrashort electron bunches are useful for applications like ultrafast imaging and coherent radiation production. Currently, however, the shortest achievable bunches, at attosecond time scales, have only been realized in the single or very few electron regime, limited by Coulomb repulsion and electron energy spread. Using ab initio simulations and theoretical analysis, we show that highly-charged bunches are achievable by subjecting relativistic (few MeV-scale) electrons to a superposition of terahertz and optical pulses. Using realistic electron bunches and laser pulse parameters which are within the reach of current compact setups, we provide two detailed examples: one with final bunches of ~1 fC contained within sub-400 as durations, and one with bunches of>25 electrons contained within 20 as durations. Our results reveal a route to achieve such extreme combinations of high charge and attosecond pulse durations with existing technology.

Here, we use ab initio numerical simulations and analytical theory to show that high-charge electron bunches of attosecond-scale durations can be produced by interfering coherent terahertz and optical frequency pulses. We study two regimes of operation: in the first regime, 5 MeV electrons [30] are compressed into bunches of 400 as duration, each containing ∼ 1 fC of charge. In the second regime, 5 MeV electrons are compressed into attobunches of 20 as duration, each containing ∼ 25 electrons. By comparison, theoretical predictions of electron bunch compression using realistic bunches have so far been limited to about 200 as [23] in the single-electron regime. Experimentally, the shortest electron bunches produced to date is about 800 as in the single-electron regime [24], and indirect measurements have revealed durations as short as 260 as [25].
In addition, we obtain fully closed-form expressions for the dynamics of electrons subject to a general combination of counter-propagating pulses. Given a specific initial electron bunch configuration, these analytical tools enable us to predict various key features of our compression scheme, such as the bunch duration at focus (maximum compression), and the final kinetic energy (KE) spread. Our analytical predictions agree well with our ab initio numerical simulation results in regimes where space charge effects are negligible. Our work complements existing theoretical formulations for the behavior of charged particles in counter-propagating electromagnetic fields, which are confined to the subrelativistic regime [21,31].
In the proposed scheme ( Fig. 1a(i)-(iii)), the counter-propagating terahertz and optical pulses interfere to form an intensity grating, which is velocity-matched to the relativistic (few MeV-scale) electrons. The ponderomotive force, which is proportional to the negative intensity gradient, compresses the electrons into a train of attobunches. Similar bunch compression schemes based on intensity gratings have been studied only in the regime where both electromagnetic pulses are at optical/infrared frequencies, where they have been shown to be practical for electron acceleration [32][33][34], and for the compression of non-relativistic, single and few-electron bunches [20,21,23,25]. Here, we show that the idea can be applied to realistic electron bunches [30,[35][36][37] using terahertz pulses of < 100 µJ and optical pulses of < 100 mJ, which can readily be obtained with today's technology [38][39][40][41]. Due to the suppression of space charge effects at relativistic energies [42,43], the resulting attobunches can hold substantially higher charge than current attobunches, which are frequently generated in the non-relativistic, single-electron regime [23][24][25].
We take advantage of the ability of our recently-developed multi-electron ab initio software (see Methods) to exactly model the behavior of electrons in non-paraxial laser fields. This is critical for accuracy since terahertz pulses from compact sources usually operate in the near-single-cycle limit and have beam waists tightly focused down to dimensions close to the central wavelength [39] in order to achieve the required on-axis field strengths. Notably, the electromagnetic wavepackets we use to model both terahertz and optical fields are finite-energy, exact solutions to Maxwell's equations. We also take into account both the near-field and far-field inter-electron interactions, which ensures the accuracy of our model in this multi-electron, attosecond pulse regime.

II. RESULTS
A. High-charge attosecond electron bunches Figure 1 presents results from 2 regimes of our study: (i) a regime where 20 as electron bunch durations are realized with a relatively low charge per bunch and (ii) a regime where 400 as electron bunch durations are realized with fC-scale charge per bunch. The durations of the compressed bunches are stated using full width at half maximum (FWHM) values. In all simulation results presented in this section, the optical (co-propagating) and terahertz (counter-propagating) pulses have central wavelengths of λ 1 = 0.65 µm and λ 2 = 300 µm respectively, are linearly-polarized (LP) in x, and propagate in the ±z direction. The electron bunches have a mean KE of KE = 5 FIG. 1. High-charge, relativistic (5 MeV) attosecond electron pulses formed by a terahertz-optical intensity grating. The scheme we study is shown in a: a(i) A co-propagating optical pulse (blue waveform) and a counter-propagating few-cycle terahertz pulse (green waveform) are incident on a relativistic electron bunch (yellow ellipse) of mean velocity v0. a(ii) The pulses overlap, forming a sub-luminal intensity grating (solid turquoise profile) which co-propagates with the electron bunch. a(iii) After the interaction, the velocity spread imparted by the grating compresses the electrons into a train of attobunches. The heatmaps in b(i)-(ii) show the electron density time-evolution using a centered coordinate system z − z . The electron density spatial distributions at the focal times are shown in c(i)-(ii). In panels b(i) and c(i), electrons with of 10 −5 % relative kinetic energy (KE) spread, and an average of 125 electrons per grating wavelength, interact with a 90.3 mJ optical pulse and a 39.0 µJ terahertz pulse, resulting in electron bunches of 19 as duration (FWHM). In panels b(ii) and c(ii), a 20 fC bunch with an (FWHM) duration of about 16.5 fs and relative KE spread of 0.146% interact with a 6.66 mJ optical pulse and a 16.9 µJ terahertz pulse, resulting in electron bunches of < 400 as (FWHM) containing approximately 1 fC of charge within FWHM. Only the two central grating wavelengths are pictured. b and c are the result of electrodynamic simulations in which non-paraxial electromagnetic fields as well as near-and far-field space charge effects are exactly taken into account.
MeV. The velocity of the intensity grating, v grating , is matched to the mean velocity of the electrons, v 0 , by choosing wavelengths, λ 1 and λ 2 , such that [20,34]: where c is the speed of light in free space. In the lab frame, the grating wavelength is given by where β = β 0ẑ = (v 0 /c)ẑ is the normalized mean velocity of the electron bunch propagating in z. The corresponding Lorentz factor is γ 0 = 1/ 1 − β 2 0 . The mean KE of the electrons is KE = (γ 0 − 1)m e c 2 , where m e is the electron rest mass.
Equation (1) shows the necessity of combining very disparate counterpropagating laser wavelengths where relativistic electrons are concerned: for v gr close to the speed of light, β 0 ∼ 1, λ 2 λ 1 is necessary. The use of relativistic electrons thus takes us into a regime beyond what has been studied for compressing non-relativistic electrons, and gives us an opportunity to leverage the developments of high-intensity terahertz pulses in combination with optical pulses in our scheme.
Figures. 1b(i) and 1c(i) show the electron density distribution obtained by averaging over 300 sets of ab initio simulations using an initial 5 MeV electron bunch containing 1250 electrons (0.2 fC) uniformly distributed across 10 grating wavelengths. The initial energy spread is σ KE / KE ≈ 10 −5 %. Relative energy spread values as low as σ KE / KE = 4 × 10 −4 % have been predicted to be possible for existing RF gun setups [44]. Each resulting attobunch has about 125 electrons contained within each λ gr , and 26 electrons within the FWHM. The non-paraxial optical and terahertz pulses have energies of 90.3 mJ and 39.0 µJ respectively. The optical pulse has a duration of 80 fs (intensity FWHM) and a peak on-axis field strength E 01 ≈ 4.96 × 10 10 V/m. The terahertz pulse has a 1 ps duration (intensity FWHM) and a peak on-axis field strength E 02 ≈ 2.95 × 10 8 V/m. Both laser pulses have the same waist radius, w 0 = 450 µm. During interaction, the bunch has a radius of about 15 µm. The full set of electron bunch and laser pulse parameters is given in the Methods section.
The second scenario, shown in Figs. 1b(ii) and 1c(ii), involves compressing a KE = 5 MeV, 20 fC (total charge), 16.5 fs FWHM duration electron bunch of relative σ KE / KE ≈ 0.146% and 8 µm radius, into a train of sub-400 as duration, fC-scale electron bunches. The electron density heatmap and distribution are averaged over 200 sets of simulation results. The initial electron bunch was modelled after the bunch experimentally demonstrated in [30] (see Methods section). Both pulsed lasers have the same beam waist: w 0 = 200 µm. The optical pulse has a duration of 30 fs (intensity FWHM) and an on-axis peak field strength E 01 ≈ 5 × 10 10 V/m, corresponding to a pulse energy of 6.66 mJ. The terahertz pulse has a duration of 1 ps (intensity FWHM) and an on-axis peak field strength E 02 ≈ 4.18 × 10 8 V/m, corresponding to a pulse energy of 16.9 µJ. Such optical and terahertz pulses are readily achieveable today in a table-top setup [38][39][40][41].
At the focus, we observe the formation of electron bunches with almost 1 fC within a FWHM duration of 367 as, as seen in Fig. 1c(ii). These results show that a combination of terahertz and optical technology can be enabling concepts for the realization of high-charge electron bunches of sub-fs durations.

B. Analytical expressions for charged-particle dynamics in counter-propagating fields
In this section, we present fully closed-form expressions for the behaviour of charged particles subject to a pair of counter-propagating electromagnetic pulses. These expressions have been used to predict verious key properties of our bunch compression scheme -including the focal time (maximum compression), the bunch duration at focus, and the final KE spread -and show excellent agreement with the results of our ab initio simulations in regimes where space charge and non-paraxial laser pulse effects are small (see Fig. 2).
We start from the Newton-Lorentz equations of motion, which describe the dynamics of electrons in arbitrary electromagnetic fields. Treating the counter-propagating laser pulses as plane waves, and considering an electron moving in an arbitrary direction, such that the transverse (x, y-direction) momenta are small compared to the longitudinal (z-direction) momentum, we obtain the normalized electron velocity long after interaction as (see Supplementary Information sections S.1-S.3 for detailed derivations): where the primes on the variables indicate that they are evaluated in the frame moving at normalized velocity β = β 0ẑ . We define this to be the primed frame. For our electron bunch compression scheme, we take β 0 as the mean normalized velocity of the electron bunch being compressed. E 0j and T j respectively refer to the electric field amplitude and pulse duration of the laser pulse labelled by subscript j, where j = 1 (j = 2) refers to the laser pulse which co-propagates (counter-propagates) with the electron bunch. ω = k c is the angular frequency, ∆θ is the relative angle between the polarization vectors of the two laser pulses (which we set to 0 for this work), φ 0 is a phase constant that depends on the carrier envelope phase of each laser pulse, and β z,i is the initial normalized electron velocity. The intensity peaks of the counter-propagating laser pulses overlap at position z = z OL and time t = t OL (in the primed frame), and we define the longitudinal electron position at the time t = t OL to be z OLe in the zero field strength limit. α a is defined as We also obtain the corresponding electron position long after interaction as where When the bunch has vanishing longitudinal velocity spread, i.e., β z,i = 0, the general expression for the focal time, defined as the time between t OL and the electrons reaching maximum compression, t comp , is: Here, K 0 = (2e 2 E 01 E 02 cos ∆θ)/(m e c 2 ). In the special case where we consider the electrons near the center of the intensity grating (z OLe ≈ z OL ) and Eq. (7) reduces to which agrees with the analytical results obtained in [21], modulo an addition factor 2/π which comes from our choice of a Gaussian pulse profile. The overlap of the optical and terahertz pulses results in a finite-length intensity grating in which electrons further from the center of the intensity grating generally experience a weaker compressive force. This effect is taken into account through the exponential factors in Eqs. (3) and (5), as well as through α b .
The results in Fig. 2 show the excellent agreement between our analytical predictions (circles) and numerical results when the laser pulses are modelled as plane waves (crosses). The discrepancy between the plane wave simulations and the exact numerical results using non-paraxial pulses (triangles) shows the importance of taking into account the transverse profiles of the focused optical and terahertz pulses in our simulations. Nevertheless, we also note that these exact results follow the trend predicted by our theory relatively well in the regime considered in Fig. 2. In Fig. 2, the 5 MeV, 10 µm-radius electron bunch was modelled using 3.75 × 10 5 particles, and has a uniform random distribution in z over a length of λ gr . The bunch is normally-distributed in x and y. The momentum spread for all cases was normally-distributed and isotropic in all directions: σ γβx = σ γβy = σ γβz . We used the following initial relative KE spreads: σ KE / KE = 0.02%, 0.06%, 0.10%, and 0.14%. The corresponding momentum spreads were σ γβi = 1.9615 × 10 −3 , 5.8848 × 10 −3 , 9.8075 × 10 −3 , and 1.3731 × 10 −2 in that same order (i ∈ {x, y, z}). All electron bunch and laser pulse parameters are listed in the Methods section. Figures. 2a(i) and 2b(i) show that a larger initial electron bunch KE spread makes it more difficult to compress the bunch unless higher laser field strengths field strengths are used. Our findings in Fig. 2 reveal the importance of low energy spread in realizing attosecond bunches. As seen in Figs. 2a(i), a change in relative initial KE spread from 0.02% to 0.14% can cause the minimum electron bunch durations at the focus to increase by almost an order of magnitude. In the limit where E 01 E 02 is small, we see from Figs. 2a(i)-(iii) that it is possible to obtain fs-scale electron bunches with a very small (practically negligible) change in energy spread, at the cost of a longer focal time. In Fig. 2a(ii), the decrease in the focal time approximately as 1/(E 01 E 02 ) agrees with the trend predicted by Eq. (7).

III. DISCUSSION
In this section, we present an overview of the electron kinetic energies which can be matched using sources of coherent light at various wavelengths. The interest in working with electrons of larger kinetic energies is due to the relativistic suppression of space charge effects, which allows shorter bunch durations to be achieved in this compression scheme. The development of intense, coherent terahertz sources on a table-top scale [38,39,46], as Fig. 3 shows, unlocks a range of electron kinetic energies spanning 4 orders of magnitude (keV to 10 MeV). By contrast, using only wavelengths falling in the optical to near-infrared regime (0.4 µm to 1.4 µm) would limit us to electron kinetic energies of 100 keV or less.
A number of laser-based sources of intense terahertz radiation, suitable for the use in the present compression scheme, have already been reported in the literature. Single-cycle and quasi-single-cycle terahertz radiation centered at 1 to 2 THz with peak field strengths on the order of 1 MV/cm have been achieved using optical rectification of LiNbO 3 using tilted pulse front pumping [38,39,46] and optical rectification of organic crystals with high non-linear constants [40]. In addition, terahertz pulse energies on the order of tens of µJ are routinely produced [45] from these compact sources. Difference frequency generation (DFG) of optical parameteric amplifiers have been used to produced narrow-band, multi-cycle pulses at mid-infrared frequencies (15-30 terahertz) and higher fields strengths of 100 MV/cm [41]. While ultra-broadband terahertz radiation can be produced using plasma ionization [47], the field strengths are typically lower than those achieved using optical rectification. However, they could potentially be used for the compression of low-charge or single-electron bunches with small energy spreads over longer focal distances. Operating regimes for electron pulse compression. The colormap shows the electron bunch kinetic energies which can be matched by a range of co-propagating and counter-propagating wavelengths (λ1 and λ2 respectively). Only the region corresponding to λ1 < λ2 is plotted. The region in which λ1 > λ2 corresponds to a counter-propagating grating. The region bounded by the black dashed lines mark out the terahertz regime [45] for λ2, and the black dotted lines within further segregate this frequency range into bandwidths that are currently attainable through optical rectification of LiNbO3 and organic crystals, difference frequency generation (DFG), and plasma ionization. The yellow star marks the 31 keV, non-relativistic case, studied in [20], whereas the blue star marks the 5 MeV, relativistic case which we study in this paper.
Our findings in Fig. 2, which show stronger compression is enabled by higher field strengths, motivate the development of more powerful compact coherent terahertz sources.
We note that greater flexibility in our choice of wavelength for matching a given electron kinetic energy can be achieved by tilting the counter-propagating pulses [21,23,25,33]. In this case, however, too large a tilt angle will lead to restrictions on the transverse size of the electron bunch. Nevertheless, the concept of tilting laser pulses could be implemented in the terahertz-optical scheme to accomodate an even wider range of electron kinetic energies.

IV. CONCLUSION
We have analyzed the feasibility of using counter-propagating terahertz and optical pulses to compress relativistic electrons into a train of attosecond-duration bunches. Due to the space-charge suppression at few MeV-scale energies, significant amounts of charge can be contained within each attobunch, compared to previously realized attobunches that have only single or very few electrons. Our ab initio simulations take the near and far-field space charge effects into account, and use exact, non-paraxial pulse profiles to model single-cycle tightly-focused terahertz pulses; this is a significant advance over previous numerical studies of similar intensity grating compression schemes, which assumed non-interacting electrons and planar or paraxial electromagnetic waves.
We presented results for attosecond electron bunch compression in two regimes. In the first, the initial electron bunch contains 20 fC of charge and is comparable to the bunches that can be produced by existing few-MeV scale electron sources. In this case, we showed that the electrons can be compressed into smaller bunches of sub-400 as durations (FWHM), each containing up to 1 fC of charge. The second case involved the compression of a low-charge electron cloud into attobunches with sub-20 as durations (FWHM), containing ≈ 26 electrons. Such short-duration bunches could be used, for instance, as sources of high quality coherent radiation through processes like inverse Compton scattering [11], Smith-Purcell radiation [48], and transition radiation [49]. We find that the realization of this scenario depends on having kinetic energy spreads which are extremely low but feasible [44]. Our results indicate that attosecond-scale electron bunches are not inherently limited to the few-to-single-electron regime, which has been the focus of other studies.
We also formulated an analytical theory that agrees well with the numerical results and provides good order-ofmagnitude estimates even for non-paraxial laser pulses. After carefully analyzing the assumptions underlying the theory, we extrapolated it to other values of the mean electron kinetic energy, spanning over four orders of magnitude. We argued that electron bunch compression can be achieved over this huge kinetic energy range via different combinations of terahertz and optical pulse frequencies. With the continued progress of research into generating single and few-cycle, intense terahertz radiation, this scheme could be used to produce a wide range of sub-fs, multi-electron bunches. Possible applications include UEI, as well as the production of high-brightness, coherent radiation from compact and accessible setups.

Numerical integration of the equations of motion
Numerical simulations of the electron-laser interaction were carried out by solving the exact relativistic equation of motion in the lab frame due to the Lorentz force: The electric and magnetic field vectors are denoted as E and B. The E, and B we use are finite-energy, exact, closed-form solutions to Maxwell's equations for non-paraxial pulses. The dimensionless velocity vector is β, and the corresponding Lorentz factor is γ = (1 − | β| 2 ) −1/2 . The charge and rest mass of the simulated charged particle are q and m respectively, and the speed of light in free space is c. The numerical integration was carried out by implementing a fifth-order, adaptive step Runge-Kutta algorithm. The pair-wise space charge interactions between all electrons were computed using exact expressions for both near-and far-fields obtained from the Liénard-Wiechert potentials. Our simulations also took into account radiation reaction via the Landau-Lifshitz formula, but we have verified that radiation reaction has negligible effect in all scenarios studied in the paper. For more details of our algorithm, see the Supplementary Information in [50] and [51].

Exact, finite-energy, pulsed solutions to Maxwell's equations
As the terahertz pulses we use are focused down to dimensions close to the size of the central wavelength and are single-cycle, non-paraxial effects, such as non-vanishing and non-neglible longitudinal fields, arise. In order to accurately incorporate fully non-paraxial pulses, we adopt the method used in [52] to construct exact solutions to Maxwell's equations as opposed to the paraxial wave equation. Expressed in terms of the electric and magnetic Hertz vector potentials, Π and Π * respectively, the field expressions are given by [53]: where Re denotes taking the real part of the complex field expression. The Hertz vector potentials corresponding to a non-paraxial, pulse linearly-polarized in x are [52]: The free space electric permittivity and magnetic permeability are denoted 0 and µ 0 respectively. The speed of light in free space is defined such that 0 µ 0 = 1/c 2 . The phasor, Ψ(x, y, z, t), which incorporates a Poisson power spectrum [54], has the following form: This phasor is valid in the sub-cycle and multi-cycle limits, and does not have unphysical negative frequencies.
We have defined We have denoted some constant phase shift as φ 0 , the central angular frequency as ω 0 , a real-valued, positive Lorentzinvariant controlling the pulse duration as s, the real-valued confocal parameter controlling the beam waist as a, and R = x 2 + y 2 + (z + ia) 2 . The imaginary number is i = √ −1 .

Initial electron phase space distribution
The electron bunches were initialized in the lab frame at transverse and longitudinal waists (minimum transverse and longitudinal extent) where the deviation of each momentum and sptial component from the mean (γβ i − γβ i , where i ∈ {x, y, z}, and r − r where r = (x, y, z)) are uncorrelated. The x and y-coordinates were initialized with a Gaussian distribution with equal standard deviations (SD) σ x = σ y , which we define to be the bunch waist radius, and centered on the propagation axis, x = y = 0. We initialized the momentum components with a Gaussian distribution in γβ x , γβ y , and γβ z . The mean momenta were γβ x = γβ y = 0, and γβ z = γ 0 β 0 , where β 0 corresponds to the bunch central velocity and γ 0 = (1 − β 0 2 ) −1/2 is the corresponding Lorentz factor. For a differential absolute kinetic energy spread ∆γ, the corresponding differential spread in the total momentum is: In order to obtain a dependence of ∆(γβ) on the individual momentum components ∆(γβ i ) for i ∈ {x, y, z}, we take the derivative of (γβ) 2 with respect to γβ: Using the chain rule and expressing the total derivatives in terms of partial derivatives, we get: where (γβ) 2 = (γβ x ) 2 + (γβ y ) 2 + (γβ z ) 2 . Once again, for small spreads away from the mean, non-zero momenta, γβ = γ 0 β 0 , γβ x , γβ y , and γβ z , we get: For the special case we consider, where the mean transverse momenta vanish and γ 0 β 0 ≈ γβ z (and γ 0 β 0 = 0), Eq. (17) simplifies to and the transverse momenta can be initialized independently of ∆(γβ) and ∆(γβ z ). We assume all spreads to be the SD: ∆(γβ) = σ γβ , ∆(γβ i ) = σ γβi for i ∈ {x, y, z}, ∆γ = σ γ . For a given σ γ , we initialized the momenta components and computed the exact σ γ corresponding to the initialized momentum distributions until the initialized values of momentum and well as energy spread were within 0.01% of the desired values. In the main text, the relative KE spread, was defined to be σ γ /(γ 0 − 1) = σ KE / KE .
To simulate transverse and longitudinal focusing of the initial electron bunch, each electron was traced back in time by by some fixed interval in the absence of all forces. This corresponds to the bunch centroid z propagating in the negative z-direction. After this re-tracing procedure, the bunch aquires a negative velocity chirp in all directions, which is representative of a bunch leaving a focusing element. This is the new initial bunch with which each simulation is carried out. We time the interaction between the counter-propagating laser pulses and the electron bunch such that laser intensity pulse peaks coincide with each other and the electron bunch at the same time.

Case 1 parameters: low-charge, low-KE spread bunch
In this section, we state in full the parameters used for the case of a low-charge, low-KE spread bunch. For a 5 MeV bunch, we have: β 0 = 0.9956760 and γ 0 = 10.7649968. The bunch charge simulated was 0.2 fC. In our simulations, each electron is represented by one particle. The electrons were initialized in z following a random uniform distribution spanning 10λ gr . This truncated bunch simulates the central portion of an electron bunch which is long enough for the central bunch region to be approximately unform in charge density (in z). The bunch radius when transverse waist is reached is σ x = σ y = 15 µm. The relative KE SD chosen was σ KE / KE = 10 −5 %, and assuming the bunch centroid propagates in only the +z-direction resulting in vanishing mean momenta in x and y, this corresponds to a longitudinal momentum spread of σ γβz ≈ 9.8072 × 10 −7 . For simplicity, we configured the bunch such that the momentum spread was isotropic: σ γβx = σ γβy = σ γβz . The bunch distribution in each phase space plane was uncorrelated at initialization, as stated in the previous section using the parameters above. The bunch centroid was then propagated backward in time by ∆t = γ 0 (10 ps) ≈ 107.65 ps in the absence of forces.
For both laser pulses, we chose their focal points to coincide with each other, centered on-axis: x f oc = y f oc = 0 m at z f oc ≈ 3.217068 × 10 −2 m. The optical and terahertz pulse intensity peak positions relative to z f oc at initial time t = 0 s were z ≈ −3.23101 × 10 −2 m and z ≈ 3.23101 × 10 −2 m respectively. Both beam waists were w 0 = 450 µm. The optical pulse intensity FWHM duration chosen was 80 fs. The terahertz pulse FWHM duration was 1 ps long. The on-axis peak field strengths of the optical and terahertz pulses at the focus are E 01 = 4.96 × 10 10 V/m and E 02 = 2.95 × 10 8 V/m respectively. Corresponding to these parameters, the optical and terahertz pulses have pulse energies of about 90.3 mJ and 39.0 µJ respectively. In this section, we state all parameters used to create a 20 fC, 7 fs (SD) electron bunch with realistic properties as well as the laser pulse parameters chosen. We adapted the electron bunch parameters from reference [30] and assumed, for simplicity, that the momentum and spatial distributions were normally distributed in x, y, and z. From Fig. 2(c) and (d) in [30], we estimated σ x = σ y ≈ 8 µm, which we have defined to be our radius. The transverse emittance values given were x,rms,n ≈ 16 nm rad and y,rms,n ≈ 20 nm rad. We take the definition of the normalized, rms emittance to be where is the expectation value operator. Assuming the bunch is uncorrelated in phase space at the focus in all 3 phase space planes, the bunch is centered on-axis: x 0 = y 0 = 0, and the bunch centroid motion is confined to z: γβ x,0 = γβ y,0 = 0, we have x,rms,n = σ x σ γβx = 16 nm rad and y,rms,n = σ y σ γβy = 20 nm rad. Using σ x = σ y = 8 µm, we have σ γβx = 2.0 × 10 −3 and σ γβy = 2.5 × 10 −3 . From Fig. 3(b) in [30], we estimate the longitudinal emittance to be z,rms,n ≈ 30 nm rad. Once again assuming that this value corresponds to the longitudinal waist, where σ z = σ τe β 0 c and (SD bunch duration) σ τe = 7 fs, we have: z,rms,n = σ τe β 0 cσ γβz =⇒ σ γβz ≈ 1.4356 × 10 −2 .
This corresponds to a relative KE SD of σ KE / KE ≈ 0.146%. Consistent with our assumption that the z-spatial distribution is also Gaussian, the bunch FWHM corresponding to σ τe = 7 fs is 2 2 log(2)σ τe ≈ 16.5 fs, which is the value stated in the main text. The uncorrelated phase space bunch distribution at initialization was propagated back by a constant time interval of ∆t = γ 0 (2.12 ps) ≈ 22.82 ps to simulate a bunch focusing transversally and longitudinally.
For both laser pulses, we once again chose their focii to coincide with each other and the interaction point, centered on-axis: x f oc = y f oc = 0 m at z f oc ≈ 6.849799 × 10 −3 m. The optical and terahertz pulse intensity peak positions relative to z f oc at initial time t = 0 s are z ≈ −6.879220 × 10 −3 m and z ≈ 6.879220 × 10 −3 m respectively. The beam waists were both chosen to be w 0 = 200 µm. The optical pulse intensity FWHM duration chosen was 30 fs. The terahertz pulse intensity FWHM duration chosen was 1 ps. The on-axis peak field strengths of the optical and terahertz pulses at the focus were E 01 ≈ 5 × 10 10 V/m and E 02 ≈ 4.18 × 10 8 V/m respectively. Corresponding to these parameters, the optical and terahertz pulses have energies of about 6.66 mJ and 16.9 µJ respectively. Parameters used to plot results in Fig. 2 In this section, we list the parameters used in our numerical results plotted in Fig. 2. Our numerical simulations for this cases do not account for space charge effects. We carried out simulations in both the plane wave and non-paraxial limits with 3.75 × 10 5 electrons.
The laser pulses used in our non-paraxial simulations have beam waists of w 0 = 300 µm. The terahertz pulse has a (intensity) FWHM duration of 1 ps. The optical pulse has a (intensity) FWHM duration of 30 fs. The optical pulse energy was kept constant at 60 mJ. For our study, 4 different terahertz pulse energies are used: 0.52 µJ, 2.06 µJ, 4.65 µJ, and 8.26 µJ. These energies correspond to values of the product of peak field strengths, E 01 E 02 , of about (to within 1%) 0.5 × 10 19 V 2 /m 2 , 1.0 × 10 19 V 2 /m 2 , 1.5 × 10 19 V 2 /m 2 , and 2.0 × 10 19 V 2 /m 2 in that order. The foci of both pulses were centered on the propagation axes, x f oc = y f oc = 0 m, and at z f oc ≈ 1.610413 × 10 −2 m. The initial optical and terahertz pulse peak positions relative to z f oc were z ≈ −1.617374 × 10 −2 m and z ≈ 1.617374 × 10 −2 m respectively. The pulse peaks overlapped each other and the electron bunch at the same time. For Fig. 2b(i)-(iii), the results for E 01 E 02 ≈ 1.5 × 10 19 V 2 /m 2 were omitted to keep the other results visible.
The plane wave simulations were carried out using plane laser pulses with a Poisson pulse profile [54], making them valid in both the single-cycle and multi-cycle limits. The laser pulse parameters used were the same as those used for the non-paraxial pulse simulations.