Temporal diagnostics of femtosecond electron bunches with complex structures using sparsity-based algorithm

Femtosecond electron bunches with complex temporal structures play a crucial role in THz generation, free-electron lasers and plasma wakefield accelerators. The ultrashort electron pulse duration can be reconstructed from the coherent transition radiation (CTR) spectrum based on prior knowledge. A weighted greedy sparse phase retrieval (WGESPAR) algorithm is developed to reduce the ambiguities for reconstructing the distribution of the beam current. This algorithm achieves better performance than iterative algorithms, especially for truncated noisy spectra of multibunch structures. Based on the WGESPAR algorithm, the complex temporal structures of femtosecond electron bunches generated from laser wakefield accelerators can be successfully reconstructed.


I. INTRODUCTION
Laser plasma accelerators [1] and bunch compression techniques [2] have a proven capability for the generation of few-femtosecond(fs) electron beams for plasma wakefield accelerators, ultrafast electron microscopy, compact ultrafast x-ray sources such as free-electron lasers and Thomson scattering sources [3][4][5][6][7][8]. The temporal diagnostics of few-femtosecond electron bunch profiles hence have become a crucial issue in enhancing the performance of these facilities.
Methods of picosecond and sub-picosecond electron bunch length measurement, such as deflecting cavities, electro-optic techniques and frequency-domain methods [9][10][11][12], have been well developed for conventional accelerators. However, the temporal characteristics of few-femtosecond electron bunches cannot directly utilize these methods due to their limited resolution. Because the amplitude of the transverse deflecting voltage is one of the critical constraining factors, high powered lasers have been introduced for the desired sub-fs resolution [13][14][15][16]. For electro-optic techniques, the resolution is strongly influenced by the crystal response function [17], and the achievable temporal resolution has been demonstrated down to a few tens of femtoseconds. Frequency-domain methods have the potential to achieve sub-femtosecond resolution by scanning the autocorrelation curve of coherent transition radiation (CTR), diffraction radiation or coherent synchrotron radiation [18][19][20].
To measure the duration of few-femtosecond electron bunches, the method of CTR has frequently been adopted. CTR is generated when electron bunches propagate through inhomogeneous media, and the longitudinal bunch profile can be reconstructed from the attained CTR spectrum by Kramer-Kronig relations [21] or the fitting of Gaussian shape bunches [1,[22][23][24]. Conventionally, Kramer-Kronig relations work well in accelerator physics, but extrapolation to the unknown frequency * jfhua@tsinghua.edu.cn ranges based on assumptions on the bunch shape is often required, which brings unavoidable uncertainty for complex temporal structures [25]. The Bubblewrap algorithm has been developed for the temporal characterization of femtosecond electron bunches, and a two-bunch structure can be reconstructed [26,27]. However, the retrieved phase would be sensitive to noise and have ambiguities, which consequently lead to longer bunch lengths or even inaccurate bunch shapes. Recently, a greedy sparse phase retrieval (GESPAR) algorithm has been applied in coherent diffraction imaging and achieves good performance with prior knowledge [28]. The GESPAR algorithm can be heuristically applied in bunch duration reconstruction from the CTR spectrum.
In this paper, a WGESPAR algorithm is developed to reconstruct the temporal profiles of few-femtosecond electron bunches with complex structures. In Sec. II, we describe the theory of CTR diagnostics and conduct a Monte Carlo simulation for CTR generated by ultrashort electron bunches passing through a radiator. In Sec. III, a WGESPAR algorithm is developed to be specifically applied in the reconstruction of complex longitudinal bunch structures from the CTR spectrum; the algorithm's performance is compared with the Gerchberg-Saxton (GS) algorithm. The effectiveness of the WGES-PAR algorithm is confirmed by the synthetic simulation of OSIRIS Particle-in-Cell (PIC) code and Monte Carlo code.

II. THEORY AND MONTE CARLO SIMULATION OF CTR SPECTRUM
Transition radiation is emitted when charged particles pass through inhomogeneous media in the general case [29]. In the transverse dimensionless case, the CTR energy spectrum density dI dωdΩ observed in the direction of θ is given by Here, N denotes the total number of electrons and ρ(ω) , where σ t is the electron bunch duration. dIe dωdΩ means the spatial energy distribution of the CTR generated by a single electron [30]: This is called the Ginzburg-Frank formula. Here, β and β denote the normalized velocity of the electrons and their image charges, respectively, andn is the unit vector pointed toward the observer. dIe dωdΩ is not correlated with the optical frequency; hence, the longitudinal profile is only related to the form factor ρ(ω).
Since the form factor of CTR is dependent on the distribution of electron bunches, the forward CTR radiation will inevitably be altered as the electron bunches pass through the radiator. Here, we use Al foil as the radiator. We conduct a series of Monte Carlo simulations on Geant4 [21,31] to estimate the actual CTR spectrum considering with the effect of the Al foil and bunch transverse distribution. In the simulation, CTR can be precisely calculated by transporting each electron separately through an Al foil and subsequently superposing the field of all electrons on the collection plane. Here, the collection plane is set to be 646 mm from the radiator, and the collection angle is 15 mrad. In Fig. 1, 100 MeV electrons are emitted randomly from a point source with a divergence of 2 mrad (rms) and a duration of 2 fs (rms). After passing through a 76 um Al foil located 1 cm from the electron source, the arrival time, position and momentum of each electron on the rear surface of the foil are recorded. From the calculated CTR spectrum, we can find that the spatial distribution of electron bunches is nearly invariable, and thus, the CTR spectrum shape is close to the optimal case ( Fig. 1), validating this method for characterizing the longitudinal structure of electron beams from the form factor. However, the simulated spectrum amplitude is nearly half the theoretical calculation under the influence of electron divergence and the limited collection angle. According to the simulation, a 10 pC electron bunch could generate more than 300 nJ CTR in the wavelength range of 1-9 um, which provides practical guidance for the implementation of beam temporal profile measurements such as the design of spectrometer layouts and the choice of mid-infrared detector. From the Ginzburg-Frank formula, the opening angle of the CTR radiation cone can be defined as θ ∼ 1 γ , where γ stands for the Lorentz factor. Since the observation angle is limited by the spectrometer aperture, the collection efficiency will be influenced by the electron energy and the polar angle between the velocity direction and the beamline. Considering the polar angle and electron energy, we quantify the collection efficiency as the ratio of the collected CTR energy in the preset collection angle to the total radiated CTR energy shown in Fig. 2. To quantify the current profile reconstructed from the CTR spectrum, the effective charge can be regarded as the initial electron charge multiplied by the square root of the collection efficiency ratio. Combining the concept of effective charge with phase retrieval algorithms, the shape of the beam current successfully reconstructed from the obtained CTR spectrum should be the same as the current contributed by the total effective charge.

III. ALGORITHMS FOR LONGITUDINAL BUNCH PROFILE RECONSTRUCTION
The reconstruction of the beam current from the CTR spectrum can be regarded as one-dimensional phase retrieval problems, which have been studied for decades. The reconstruction process can be summarized as seeking a current profile ρ(t), the Fourier transformρ(ν) = +∞ −∞ ρ(t) exp (−2πiνt)dt of which has the same ampli-tude as the measured CTR spectrum |ĝ(ν)|, and the missing phase information is expected to be retrieved from the CTR spectrum modulus. For physical considerations, the function ρ(t) must be nonnegative and constrained in a finite support.
Several algorithms, such as the Kramer-Kronig relations [21] and iterative algorithms, including the GS algorithm and Hybrid Input-Output Algorithm (HIO) [32], which have been applied in astronomy, Xray diffraction imaging and other applications, have been developed to solve the phase retrieval problem. The emergent Bubblewrap algorithm is another iterative algorithm based on the GS and HIO algorithms. Iterative algorithms can reconstruct longitudinal bunch profiles with sufficient frequency domain information, but they may lead to longer bunch tails if low-frequency information is missed. Recently, it has been proved that one-dimensional phase retrieval problems still have intrinsic "zero-flipping" ambiguities even with the assumption of finite support and non-negativity [33,34]. With the development of compressed sensing in the past few years, sparsity-based algorithms [35] have been widely applied in the field of signal processing because the reconstruction is robust even with highly noisy spectra. Accordingly, optimization-based algorithms, such as semidefinite programming, the GESPAR algorithm and the sparse Fienup algorithm [36], have been proposed to solve phase retrieval problems with a sparsity priori. The GESPAR algorithm has been recently adapted in coherent diffraction imaging and achieved good performance with sparsity as prior knowledge [28]. We develop a weighted GESPAR (WGESPAR) algorithm for electron bunch profile reconstruction, and we compare its performance with the GS algorithm.
A. One-dimensional phase retrieval algorithms The extensively used GS algorithm and HIO algorithm are error-reducing algorithms used to obtain a better approximation at each iteration. Suppose that the Fourier transform of the bunch profile ρ(t) is given bŷ ρ(ν) = |ĝ(ν)| e −φ(ν) . Here, ρ(t) is a real, non-negative signal with compact support, and |ĝ(ν)| is the square root of the band-limited spectrum measured in the experiment. A random function ρ 0 (t) can act as the initial beam current profile. The k-th trial solution ρ k (t) is Fourier transformed, yielding |ρ k (ν)| e −φ k (ν) . To satisfy the Fourier-domain constraints, a better estimation ofρ(ν) is given byρ k (ν) = |ĝ(ν)| e −φ k (ν) . Then, the resultingρ k (ν) is inverse Fourier transformed, yielding the function ρ k (t). ρ k (t) is then modified to satisfy the timedomain constraints, yielding ρ k+1 (t), which is regarded as the input value for the next iteration. Normally, a relaxation factor r is introduced to improve performance. The relaxation iterative format is called an HIO algorithm. When r = −1, the HIO algorithm simplifies to the GS algorithm. In iterative algorithms, such as GS and HIO, overfitting usually occurs, leading to unphysical results when some information is missing and the spectrum data are highly noisy. From a different perspective, the GESPAR algorithm considers the phase retrieval as a nonlinear least-squares optimization problem, therein minimizing the error E = ρ 2 − |ĝ| 2 between the square of the discrete Fourier transform of the reconstructed signalŝ ρ 2 and the measured spectrum |ĝ| 2 in a confined set of signals. Here, bothρ andĝ are vectors. Substitutinĝ ρ = F ρ as the discrete fourier transform(DFT) of the time-domain signal ρ by multiplying it with the DFT matrix F , the objective optimization function can be written Here F i is the ith row of the DFT matrix F . The non-negative signal vector ρ can be approximated by a projection on a specific base set, for example, a base set of Gaussian bunches. The approximated signal takes the form ρ = Dx, where D is the dictionary of the preset base set and x means the projection of the signal ρ onto that base. We make an assumption that the electron current profile can be approximated by a linear combination of less than s bases defined by the dictionary D. With those assumptions, the approximation of the signal ρ can be solved uniquely by minimizing the error function E.
For practical spectrum measurements, the noise introduced by the detector should be considered. According to the relation dω = − 2πc λ 2 dλ, the spectrum data acquired in the experiment will be characterized by different noise levels, which will influence the reconstruction results if all the spectrum data are given the same weight. To solve this problem, the optimization function is rewrit-ten as the weighted square of the error between the reconstructed spectrum and the approximated real spectrum where w i is the normalized weight introduced as an evaluation of the importance of each data point. Here, , where σ i stands for the standard variance of spectrum data i. Finally, we obtain the mathematical formulation of the phase retrieval problem subject to the non-negativity, finite support and sparsity constraints: Here, x 0 stands for the number of non-zero elements of x. Note that the dictionary is chosen in the support so that the support constraint can be automatically satisfied. For example, if the dictionary is chosen as a matrix based on Gaussian pulses with finite support, the reconstructed signal can only be the linear combination of those Gaussian signals such that it also satisfies the support priority.
The flow chart of the WGESPAR algorithm is shown in Fig. 3, which is different from its original form [26] considering the non-negative condition and the weighted optimization function. Let n be the index array of the non-zero elements of the vector x, where x n means the non-zero element vector of x such that the zero norm of x is equal to the length of the index vector n: x 0 = |n|. The index array has a random initialization value n 0 with a length of s. At the k-th iteration, the optimization problem arg min is solved by the non-linear least-square trust-region algorithm [37] and yields x nk and E n k with index array n k . Then, we swap the index of M components of x nk with the smallest absolute value with the P negative maximums of ∇E(x). After this swap, the index array becomes n k , and E n k is yielded by solving the optimization function again. If the optimization function value E n k is smaller than E n k , the index vector n k is reserved as the non-zero element index of the next loop n k+1 until the value E n k of each swap is higher than the original E n k , and we accept Dx nk as the reconstructed beam current.

B. Comparison of WGESPAR and GS algorithms
To compare the WGESPAR algorithm with the GS algorithm, we generate pulses with random positions and durations in a finite support. The spectra of those test signals are sampled with a frequency step of 10 THz deduced from the resolution of a commercial mid-infrared detector. We apply the GS and WGESPAR algorithms to rebuild the phase of the truncated spectra within the range of 50 THz-500 THz. In the GS algorithm, we do not apply restrictions in the frequency range below 50 THz. The base of the WGESPAR algorithm is composed of Gaussian pulses. There are two ranges of pulse widths in that base, one from 1 fs to 4 fs with 0.5 fs increments and one from 0.5 fs to 1 fs with 0.1 fs increments. The positions of all Gaussian pulses are mapped uniformly on the support interval, and the position interval is set to be 0.5 fs. The dictionary is generated corresponding to the Gaussian base. During the reconstruction, the WGES-PAR algorithm is run several times until the error is less than a preset threshold or the run times are larger than a preset number N . We set M = 2, P = 30, N = 5. M and P depend on the sparsity and complexity of the selected base.
For the GS algorithm, the retrieved pulse length has been mentioned to be elongated because of a lack of information at longer wavelengths [26]. From Fig. 4, we can find the same phenomenon. In contrast, with the assumption of bunch series whose pulse width of every single bunch is less than 5 fs, bunch trains with a support length of 40 fs can be well reconstructed by WGESPAR even without frequency information below 50 THz.
In the experiment, the influence of noise should be seriously considered for bunch profile reconstruction. Typical MIR pyroelectric cameras usually have a noise level of a few nJ, which is only two orders less than the expected CTR energy. To test the stability of these two algorithms, 10% gaussian noise is added to the virtual "detector" pixel. When the frequency goes to zero, the noise of the spectrum goes to infinity because of the relationship dω = − 2πc λ 2 dλ. Fig. 4(b,d) shows the temporal distribution rebuilt from the noisy spectrum. We can find that the GS algorithm is more sensitive to noise and fails to reconstruct the signal, even though it can work for the ideal case, whereas the WGESPAR algorithm tends to show a better tolerance to noise. This is because iterative algorithms conform to all noisy data of the spectrum, whereas optimization methods only find the signal in the preset group whose spectra are closest to the measured spectra.
Electron beams with specific current profiles have special applications in plasma wakefield accelerators and THz generators [3,4]. In Fig. 5, WGESPAR shows good performance for special bunch profile reconstruction such as triangular and multi-bunch spikes superposed on long bunch backgrounds. The reconstructed beam current maintains main structural features such as the slope of the bunch and the bunch train intervals. If other prior-ies (such as an approximation of the bunch shape) are provided, the reconstruction result will be improved.

C. Beam current reconstruction for combined simulation of PIC and Monte Carlo
In the previous section, we demonstrate that the WGESPAR algorithm can effectively reconstruct complex current structures by reducing the number of feasible solutions. To test the WGESPAR algorithm for the reconstruction of real beam currents, we have performed 2D PIC simulations for downramp injected electron beams using OSIRIS. A 35 fs (rms), 800 nm laser is focused at the front edge of a plasma using a spot with a diameter of 10 um (rms), and the peak plasma density is set as 1.74 × 10 19 cm −3 , with a sharp transition at 100 um. Strong downramp injection of electrons with a twopulse current structure have been observed in the simulation. The injected quasi-monoenergetic electron bunch can be accelerated to 100 MeV with an energy spread of 40%( Fig. 6(a)). Using the accelerated electron bunches as the radiation source, the process of CTR generation is simulated by the Monte Carlo method discussed in Sec. II. Due to the limited collection aperture for CTR collection in experiment, electrons with energy higher than 20 MeV are considered as the CTR source. The longitudinal electron bunch profiles are reconstructed from spectra simulated by the Monte Carlo code with a sampling rate of 5 THz from 35 THz to 350 THz.
We use the WGESPAR algorithm for phase retrieval with a support of [−30, 30] fs and the non-negativity constraint. A Gaussian base with pulse width varying from 1 fs to 4 fs, with an increment of 0.5 fs, is selected. We set the WGESPAR algorithm parameters as follows: M = 2, P = 30 and N = 5. The sparsity parameter varies from 2 to 10, and the result with the lowest sparsity whose spectrum deviation from the measured spectrum is within a preset threshold is chosen. Fig. 6 shows a typical reconstruction result.
We have previously discussed the dependence of the CTR spectrum on the electron energy and the polar angle. The effective charge is introduced by taking these two factors into account. In the sense of effective charge, the WGESPAR algorithm can successfully reconstruct the main features of two peak beam currents (Fig. 6), where the triangular shape of the sub-bunch can also be clearly recognized. The relative deviations of the current amplitude, beam intervals and bunch length of each subbunch are less than 5%. A plateau after the main bunch is missing because of the lack of information under 35 THz.

IV. CONCLUSIONS
It has been proved that fs beam current profiles can be reconstructed using a CTR spectrum with sufficient bandwidth. However, the achieved spectra in the experiment are usually truncated due to the experimental layouts, which can lead to the inaccurate reconstruction of complex longitudinal bunch structures. An optimizationbased algorithm, WGESPAR, has been developed for the temporal characterization of electron bunches even with truncated noisy spectra. We find that the WGES-PAR algorithm exhibits better tolerance than iterative algorithms to noisy spectra with a sparsity assumption. Using the synthetic simulation with Particle-in-Cell and Monte Carlo codes, the WGESPAR algorithm can be successfully applied to measure ultrashort electron bunches with complex temporal structures generated in laserplasma accelerators, thereby paving the way for applications in ultrafast science with fs electron beams.