Analytic calculation of the electric field of a coherent THz pulse

Accelerators can serve as sources for intense THz radiation by producing sub-ps long electron bunches, which generate synchrotron radiation coherently in the THz regime. In this paper, we present an analytic method to calculate the electric field pulse emitted by an electron bunch with arbitrary bunch shape. Our method is applied to a Gaussian bunch as well as a complex bunch profile resulting from a nonlinear beam dynamics simulation.


I. INTRODUCTION
Coherent radiation is produced whenever the emitted wavelength is comparable or larger than the source size.Its intensity increases with the number of emitters squared, as opposed to the incoherent radiation, which increases only linearly.Since we take an electron bunch as a source, consisting of several hundred million electrons, coherent radiation from electron bunches is enhanced over incoherent radiation by this factor.For coherent radiation to be in the THz regime, the electron bunch must be shorter than 1 ps.Progress in accelerator science has enabled the creation of coherent THz radiation in storage rings [1,2].However, further decreasing the bunch size, and, therefore, reaching frequencies above 1 THz, is limited by instabilities.These instabilities are not present in linear accelerators.The linear accelerator FLUTE [3], which is currently under construction at Karlsruhe Institute of Technology, allows the production of bunches with bunch lengths ranging from 100 fs to a few fs and charges from a few nC to pC, respectively.
Because the coherence condition depends on both the bunch shape and the frequency, the spectrum of coherent radiation contains information on the longitudinal bunch profile.Conversely, due to the large amount of electrons, even small substructures in a bunch can lead to a significant production of coherent radiation.Thus, coherent radiation offers a tool for the noninvasive study of bunch profiles and beam dynamics [4].
Besides the spectrum of the emitted THz radiation, the temporal profile as well as the peak value of the corresponding electric field pulse are important [5].To design an accelerator that optimizes these parameters requires the calculation of the electric field pulse from the bunch profile.
For this, one needs to sum the electric fields at a given frequency emitted by each electron and then integrate over all frequencies.Since the number of electrons in a bunch is large, the sum can be replaced by an integral over the bunch profile.The resulting double integral could, in principle, be performed numerically for an arbitrary bunch profile.However, it turns out that the numerical integration converges poorly and we resort to analytic methods in this paper.For Gaussian profiles the pulse shape of a coherent THz pulse has been calculated analytically in [6,7].In this paper, we present an analytic method that applies to arbitrary bunch shapes.The method is based on first doing a cubic spline interpolation of an arbitrary bunch profile.Subsequently, the integrals over the bunch profile and frequency are computed analytically and expressed in terms of the coefficients of the spline interpolation.
The paper is organized as follows: In Sec.II we derive how the electric field pulse can be calculated from a spectral integral.The field is then calculated analytically for a Gaussian electron bunch profile in Sec.II A. Next, in Sec.II B, we derive our main result: an expression for the electric field of an arbitrary bunch shape based on a cubic spline interpolation.The two methods are compared in Sec.III A for a Gaussian bunch.Finally, we employ the spline interpolation in Sec.III B to compute the electric field of a simulated complex bunch profile.Section IV discusses the results.

II. DERIVATION
The electric field pulse coherently emitted by an electron bunch is obtained by summing over the individual fields of all electrons [8,9].Here, we assume that the field of a single electron, as detected by an observer at distance R and time t, is given by a plane wave propagating in the positive z-direction [7] E ω ðR; tÞ ¼ ReðE 0 ðωÞe −iωðt−R=cÞ−iϕ Þ; ( where c denotes the speed of light.The details of the radiation mechanism result in a frequency dependence of the spectral amplitude E 0 ðωÞ.In the following, the phase ϕ is assumed to be constant for all electrons.For two electrons, the field is the superposition of two plane waves, but the wave of the second electron has an additional phase factor of iωΔτ, which takes the relative displacement Δτ of the two electrons into account.Summing over all N electrons we obtain In the second line, we took the continuum limit and replaced the sum over all electrons by an integral over the normalized charge distribution ϱðτÞ.Notice that the integral is the Fourier transform F of the charge distribution and ϱðωÞ ≡ F ½ϱðtÞ in the following.Finally, we integrate over the whole spectrum and the result reads In the following, we set R to zero in Eq. ( 3) for convenience, as it is merely a constant shift in observation time t.To solve Eq. ( 3) we must first calculate the Fourier transform ϱðωÞ and then solve the spectral integral.
Our main focus is on coherent synchrotron radiation (CSR) as the radiation mechanism.As the spectral amplitude E 0 ðωÞ we take the angle integrated synchrotron spectrum [10,11] with e and ϵ 0 denoting the electric charge and vacuum permittivity, respectively.In the second line, we used the low frequency approximation, valid for ω ≪ ω c .Finally, we substituted the definition of the critical frequency ω c ≡ 3γ 3 c=2ρ in terms of the bending radius ρ.Notice that the final result is independent of the energy γ.Also notice that the use of Eq. ( 4) in the spectral integral of Eq. ( 3) would lead to a UV divergence, if it were not for the highfrequency cutoff provided by the Fourier transformed charge distribution ρðωÞ.The latter goes to zero for frequencies larger than the inverse bunch length σ b .The use of the low frequency spectrum is, therefore, justified as long as 1 ≪ σ b ω c .If this condition is not satisfied, i.e., for very short bunches, the upper limit of integration in Eq. ( 3) must be replaced by ω c .The cutoff in this case is provided by ω c because the synchrotron spectrum decays as ffiffiffiffiffiffiffiffiffiffiffi ω=ω c p expð−ω=ω c Þ for frequencies larger than ω c [10].For the remainder of the article, we assume 1 ≪ σ b ω c and use Eq. ( 4).
The main problem is then to compute for a given charge distribution.The coherent electric field pulse of a bunch with charge Q is then given by

A. Gaussian bunch
In this section we consider a Gaussian bunch of width σ b , that is, the profile is given by The Fourier transform is well known and reads explicitly Inserting Eq. ( 7) into Eq.( 5) and splitting into real and imaginary parts leads to These integrals were solved using Eqs.(3.952,7-8) of [12] together with the property expðxÞ 1 F 1 ðb − a; b; −xÞ ¼ Figure 1 shows a plot of the electric field for a Gaussian bunch with σ b ¼ 50 fs for different values of the phase ϕ.For ϕ ¼ 0 we get a so-called half-cycle pulse which is symmetric and centered at t ¼ 0, corresponding to the center of the charge distribution.Notice that the pulse shape resembles the underlying Gaussian charge profile.For 0 < ϕ < 180°a minimum develops at t > 0 (see dashed curve in Fig. 1), leading to a single-cycle pulse, which is antisymmetric at ϕ ¼ 90°.The pulse shape for ϕ ¼ 180°is the negative of the one for ϕ ¼ 0. Phases with 180°< ϕ < 360°have a minimum at t < 0 (see dotted curve in Fig. 1).However, the integral over all time is zero in all cases [13].This implies that the larger the minimum, the smaller the maximum field becomes, as can be seen in Fig. 1.
An arbitrary bunch profile can be written as a sum of Gaussians.The electric field, being linear in ρðωÞ, is then a superposition of fields given by Eq. ( 9).Using a multi-Gaussian fit is possible, but becomes tedious in practice.

B. Arbitrary bunch shape
We proceed in three steps to analytically compute Eq. ( 5) for an arbitrary bunch shape.The first step is to use a spline interpolation for the bunch profile.Second, we compute the Fourier-transform ϱðωÞ of the interpolated profile.
The integral in Eq. ( 5) can then be calculated analytically in the third step.

Spline interpolation
We assume that the bunch profile is given at N þ 1 discrete, not necessarily equidistant, points ðt i ; ϱ i Þ, with i ¼ 0…N.The spline interpolation connects these points by cubic polynomials, such that the resulting curve is continuous up to the second derivative.The cubic polynomial ϱ i ðtÞ that interpolates between the points ðt i−1 ; ϱ i−1 Þ and ðt i ; ϱ i Þ is given by where Δt i ≡ t − t i .The Heaviside functions ΘðtÞ ensure that the polynomial "lives" only between t i−1 and t i and is zero outside this interval.For each of the N intervals ½t i−1 ; t i the coefficients c 3 to c 0 need to be determined such that In terms of the second derivatives at the data points Of course, Eq. (11c) can be simplified, but it will be more convenient in this form later.The p i themselves can be computed very efficiently numerically [15], given the free boundary values p 0 and p N .Here, we use the boundary values for natural splines, that is

Fourier transform
With the coefficients in Eq. ( 10) expressed in terms of the data points, we can now proceed to compute the Fourier transform ϱ i ðωÞ.From Eq. (10) we see that we have to compute integrals of the form Z c n;i t n ½ΘðΔt i−1 Þ − ΘðΔt i Þe iωt dt;  9).The parameters are Q ¼ 100 pC, R ¼ 1 m, ρ ¼ 1.1 m, and σ b ¼ 50 fs.For ϕ ¼ 0 (continuous) we get a symmetric socalled half-cycle pulse whereas ϕ ¼ 65°(dashed) and ϕ ¼ 325°( dotted) yield single-cycle pulses with minima at t > 0 and t < 0, respectively.
ANALYTIC CALCULATION OF THE ELECTRIC FIELD OF … Phys.Rev. ST Accel.Beams 17, 050701 (2014) where we used Eq.(17.23,6) of [12] for the Fourier transform of the Heaviside function ΘðtÞ and δðωÞ denotes the Dirac δ-function.Results for other powers of n can be obtained by applying F ½t n fðtÞ ¼ ð−iÞ n f ðnÞ ðωÞ; to Eq. ( 13), where f ðnÞ denotes the nth derivative of a function f.The contribution from the δ-function vanishes because the exponential functions cancel for ω ¼ 0.
The same is true for the higher order derivatives of the δ-function appearing for n > 0 [14].We then have The sums with odd powers of ω can further be simplified by inserting the values of the c n;i from Eqs. (11) [14].For the sum proportional to 1=ω 3 we obtain In the last step, we used the boundary condition for natural splines, Eq. (12).Similarly, the sum proportional to 1=ω reduces to Here, we assumed that the charge distribution approaches zero at its endpoints.The Fourier transform in Eq. ( 14) then simplifies to
The second problem requires more effort.Its root is the IR behavior of the integral.Using a Taylor expansion around ω ¼ 0 we find the bunch compressor is shown in Fig. 3a [17], together with a fit of the profile to four Gaussians (continuous line).The simulations take into account the self-interaction of the bunch with its own radiation field in the compressor.This self-interaction gives rise to the complex substructure.Depending on the detector resolution of an experiment, the measured bunch profile will be more smooth.Here, we are mainly interested in benchmarking our methods for a complex bunch.The electric field pulse was calculated from this using [17] ρ ¼ 1.108 m, R ¼ 1 m, and ϕ ¼ 0°. Figure 3b shows the electric field pulse using Eq. ( 9) for the multi-Gaussian fit (continuous curve) and Eq. ( 19) for the spline interpolation (points).Again, the two methods agree.Notice that the fields from the two methods closely follow their respective charge profiles for ϕ ¼ 0°.Since the multi-Gaussian fit does not "see" the substructure at around −30 fs in Fig. 3a, its electric field does not reproduce it in Fig. 3b either.Similarly, since the multi-Gaussian fit correctly finds the position of the maximum current at about 53 fs but underestimates its magnitude by about 7%,  9) (continuous) as well as the spline interpolation and Eq.(19) (points).Notice how, for the phase ϕ ¼ 0°used here, each pulse shape follows its underlying temporal charge profile.

FIG. 3 .
FIG. 3. (a) Simulated longitudinal charge profile of a 100 pC bunch at FLUTE.Points are simulated data, while the continuous line is a fit of the data to a sum of four Gaussians.(b) Calculated electric field E based on the multi-Gaussian fit and Eq.(9) (continuous) as well as the spline interpolation and Eq.(19) (points).Notice how, for the phase ϕ ¼ 0°used here, each pulse shape follows its underlying temporal charge profile.