Electron bunch profile reconstruction based on phase-constrained iterative algorithm

The phase retrieval problem occurs in a number of areas in physics and is the subject of continuing investigation. The one-dimensional case, e.g., the reconstruction of the temporal profile of a charged particle bunch, is particularly challenging and important for particle accelerators. Accurate knowledge of the longitudinal (time) profile of the bunch is important in the context of linear colliders, wakefield accelerators and for the next generation of light sources, including x-ray SASE FELs. Frequently applied methods, e.g., minimal phase retrieval or other iterative algorithms, are reliable if the Blaschke phase contribution is negligible. This, however, is neither known a priori nor can it be assumed to apply to an arbitrary bunch profile. We present a novel approach which gives reproducible, most-probable and stable reconstructions for bunch profiles (both artificial and experimental) that would otherwise remain unresolved by the existing techniques.


I. INTRODUCTION
The problem of retrieving the phase of a signal from a measurement of its power spectrum alone is well known and has been under investigation for a few decades [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Measurements in the frequency domain (power spectrum) are necessary in cases where a direct measurement in the time domain is not possible, for instance because of the limited response time of the available detectors. Most detectors, however, are capable of measuring the power of the signal but the consequent loss of information about the phase is a fundamental problem. In accelerator physics the longitudinal bunch profile is a parameter that becomes progressively more difficult to determine for picosecond (ps) and sub-ps long bunches and micromodulated bunches. These, however, are the bunches that are expected from the next generation of high-gradient particle accelerators (based on laser or beam-driven wakefield acceleration). Interest in this problem was boosted by the realization that future light sources can be driven by fs-electron pulses and that knowledge of the bunch's longitudinal profile is essential for the operation of the accelerator and for the control of the light source stability.
The nondestructive techniques which are used to determine a bunch profile are based on the analysis of the amplitude of the coherent radiation spectrum generated by the bunch and the application of a phase recovery technique [1,5,6]. Some information about the missing phase can be obtained from the so-called "minimal" phase [θ m ðωÞ]. The phase θ m ðωÞ can be calculated by the Kramers-Kronig (KK) method under the assumption that the signal, i.e., the radiation spectrum ρðωÞe iθ m ðωÞ , is a holomorphic function: Published by the American Physical Society under the terms of the Creative Commons Attribution 3.0 License. Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
The quantity ρðωÞ in (1) is the magnitude of the Fourier transform of the signal, which can be derived from the spectral yield. The minimal phase provides a lower limit for the estimate of the missing phase. Although there are situations where the minimal phase is adequate, this is not always the case because of the possible existence of additional contributions, known as the Blaschke phase [10][11][12] which cannot be known from the measured spectral intensity. In addition, the calculation of the minimal phase [θ m ðωÞ] [Eq. (1)] assumes knowledge of ρðωÞ over the whole of the frequency domain. Since this is not possible, it is always necessary to extrapolate beyond the measurement range and to interpolate within it. Alternatively, iterative methods [13][14][15] have been used for the determination of the missing phase; these can give stable solutions for the 2D and 3D cases [16] but it is known [13][14][15][16] that the problem remains unresolved in the 1D case, such as the determination of the time profile of ultrashort bunches of charged particles. This is of particular interest to us because of its importance in relativistic beam diagnostics and free electron lasers [1][2][3][4][5][6][7]. Present and future particle accelerators will produce fs-long bunches whose shape is unlikely to be a simple Gaussian [17][18][19][20] and, given the importance of the temporal profile of the bunch for the understanding of beam-beam interactions and in FEL lasing, any advance in the phase retrieval problem would be highly desirable.
In this paper we present a novel iterative algorithm for the retrieval of the missing phase that gives reproducible, most-probable and stable solutions to the problem of the reconstruction of the time profile of a bunch from a measurement of the spectral yield of a radiative process induced by this bunch, e.g., by Smith-Purcell radiation. We start with a brief overview of the new algorithm and we explain how it differs from other iterative techniques. We then discuss the identification of the most probable solution and the convergence to this solution. We demonstrate the advantages and the limitations of the new algorithm by means of examples of profiles that are usually unresolved by the existing algorithms. We note that the algorithm is applicable to any 1D problem in which: (a) the temporal and spectral amplitudes are linked by a Fourier transform; (b) the temporal amplitude is known to be nonzero only FIG. 2. Reconstruction of a synthetic function (black dotted line) generated by the superposition of two Gaussian functions with widths (σ i ) and amplitudes (A i ) given by The reconstructions are based on: (a) a non-PCI algorithm (the solid and dashed lines show the profiles that are reconstructed when using two different randomly generated sets of initial phases); (b) the PCI (solid red line) and minimal phase (solid green line) algorithms. The dotted black line indicates the original function. The inset to (a) illustrates the absolute integral errors (4) for both realizations shown in (a), illustrating the variation of the convergence trajectory with different initial phases. Fig. 2  within a known time window (τ) and is assumed to be zero outside this window; (c) the value (F) of the integral of the time profile over this window is known (for a chargenormalized distribution F ¼ 1). Finally, we apply the new algorithm to the analysis of spectral yield measurements from FACET (SLAC) [1,[5][6][7]21].

II. THE PHASE-CONSTRAINED ITERATIVE (PCI) ALGORITHM
The new algorithm is based on a combination of the minimal phase θ m and an iterative procedure (with repeated iterations between the frequency and time domains). The Kramers-Kronig (KK) method is used to calculate the minimal phase. We assume that the power spectrum is known over the whole frequency domain by measurement and by interpolation/extrapolation. Since the minimal phase θ m provides a lower limit for the estimate of the phase at each iterative step, we call the new algorithm a phase-constrained iterative (PCI) algorithm. The time window (τ) is usually set at about 5-10 times the expected bunch length which is known independently from numerical simulations and from measurements of the beam transport through the accelerator. In addition to the minimal phase, the integral value F is also used as a limiting condition and in order to monitor the convergence of the procedure to a solution (as discussed below). These are the essential features of the new algorithm. We note that displacements along the time axis, mirror imaging of the temporal profile or, equivalently, the sign of the phase cannot be determined unambiguously but for our specific problem, displacements along the time axis have no physical significance.
The operation of the PCI algorithm [13][14][15] is shown in schematic form in Fig. 1 and summarized by (2): where ϒ is the set of points that violate the constrains on the target function fðtÞ while f n ðtÞ is the estimate of fðtÞ after the nth iterative step and θ n the corresponding phase. The primed quantities are the estimate of fðtÞ after the FT −1 from the frequency domain, i.e., before it has been tested against the time-domain constraints (the abbreviations FT and FT −1 denote the forward and inverse Fourier transforms, respectively); θ m is the minimal phase (1). The iterative part of the PCI algorithm is similar to the Gerchberg-Saxton and hybrid input-output (HIO) group of algorithms and incorporates the concept of the adaptive change of the free parameter β and of the width of the support of the function fðtÞ during each iteration step. The term "support" in the context of this paper is the set of points that satisfies the condition fðtÞ > 0. The quantity β is a "feedback" or "mixing" parameter which changes progressively during the iterations in order to adjust the support as discussed in [13][14][15][16]. The details of the iteration steps are as follows.
(1) The measured amplitude spectrum (interpolated and extrapolated as needed) and the minimal phase θ m are used for the initial inverse transform from the frequency into the time domain in order to derive the zeroth order approximation f 0 0 ðtÞ to the unknown profile fðtÞ. This is the starting point of the algorithm. Any values of that function inside the defined window (τ) which do not satisfy the constraints f 0 n ðtÞ ∉ ϒ [e.g., f 0 n ðtÞ < 0 in the case of electron bunch profiling] are corrected in accordance with a specific procedure. If, for example, f 0 n ðtÞ < 0 then the f nþ1 (t) value could be set to 0. An alternative approach would be to "mix" the new value with the previous estimate, f nþ1 ðtÞ ¼ f n ðtÞ − βf 0 n ðtÞ. In this case the parameter β is set initially at around 1 and must be optimized to avoid algorithm stagnation. In our case we have used the concept of the adaptive change of β and of the width of the support function during each iteration step. We note that the understanding of the influence of β on an iterative algorithm performance is still rather superficial and that it requires further investigation.
(3) If θ < θ m then θ is set equal to θ m . The values of the moduli are discarded.
where ε 1;2 are the set of predefined acceptable errors, appropriate for a specific experiment. The expressions (3) and (4) are the relative and absolute integral errors, respectively, and ε 1;2 ≪ 1. The first is an indication that the algorithm is converging to a solution and the second that this solution is likely to be the correct one. In our experience, it is advisable to run the PCI algorithm twice in order to understand the convergence behavior of the algorithm: the first time with as a large number of steps as practicable in order to identify with some confidence the iteration number for which both (3) and (4) are minimized.

III. COMPARISON OF THE PCI WITH THE ALTERNATIVE TECHNIQUES
To demonstrate the performance of the new algorithm and to compare it with other iterative techniques as well as with the Kramers-Kronig method, we have reconstructed a number of "synthetic" time profiles generated by the superposition of two Gaussian functions of varying widths (σ i ) and amplitudes ( [2][3][4][5] and the Lorentz function f L ¼ 1 π α ðt−t 0 Þ 2 þα 2 (Fig. 6); the Gaussians are centered at t i and the Lorentzian at t 0 . The reconstructions have been compared with those derived from non-PCI algorithms that belong to the group known as "projection algorithms" (e.g., Bubblewrap, Shrinkwrap, etc.), using as criteria, conditions (3) and (4) as well as the error estimate FIG. 5. As in Fig. 2 , where x j and x 0 j are the values of the reconstructed and original functions, respectively. The values of ϵ are shown on the figures. This is only possible for the synthetic profiles (where the original function is known) while in a real experiment this quantity would be unknown. Moreover, this parameter is not particularly sensitive to rather obvious differences in profile shape (see Figs. 3 and 4). In Figs. 2-6 the original profile is indicated by the dotted black lines. Note that all profiles in Figs. 2-6 have been superimposed and normalized to unit amplitude for display purposes only. Figure 2(a) shows the original function (black dotted line) and the solutions (blue dashed and red solid lines) obtained by running the same conventional (i.e., non-PCI) iterative algorithm twice, but each time with different random seed phases. As seen in Fig. 2(a), this algorithm generates different solutions for each run. This observation is also reflected in the behavior of condition (4) as a function of iteration number, shown in the inset to Fig. 2(a). Figure 2(b) shows the profile reconstructions obtained by using the minimal phase and the PCI algorithm (solid green and red lines, respectively), together with the original function (dotted black line). The minimal phase approach generated a perfect result in this case but it is unknown a priori if the minimal phase is sufficient. Figures 3(a), 5(a) and 6(a) show the profiles reconstructed by means of the PCI (solid red line) and a non-PCI (dashed blue line) iterative algorithm, known as "Bubblewrap" [15], as well as the original function (dotted black line). Figures 3(b), 5(b) and 6(b) show the profile reconstruction based solely on minimal phase (solid green line) together with the original profile. The case of Fig. 3 demonstrates the strong deviation of the profile reconstructed using the non-PCI algorithm from the original profile, while the PCI procedure delivers an acceptable approximation to the given function. Figure 4 refers to the same profile as that of Fig. 3 and is used to demonstrate the stability of the solutions provided by the PCI algorithm. The behavior of conditions (3) and (4)  and 180, as marked by circles. The profiles recovered by stopping the PCI at the above iteration numbers are nearly identical and are shown in Fig. 4(b). The small differences between these profiles are due to the finite accuracy of the numerical analysis. Figure 5 illustrates a case where all algorithms are capable of generating a solution that is close to the original function but, as mentioned earlier, unlike other algorithms the PCI algorithm gives consistently unique results which satisfy conditions (3) and (4) simultaneously.
One of the most difficult functions for reconstruction is the Lorentz function (Fig. 6).This is also a case where both the non-PCI and the PCI algorithms produce satisfactory approximations to the original function. The inset of the figure is a comparison of condition (4) for the two algorithms which shows that at iteration n ¼ 13 both algorithms have minima (δ n ¼ 0.063 and 0.066 for the non-PCI and PCI algorithms respectively). This is also a case where the minimal phase approach fails. Figure 7 shows the temporal profile reconstruction of an electron bunch studied at FACET [1,[5][6][7]21]. The dashdotted line represents the profile derived by means of the minimal phase technique while the solid and dashed lines show the profile reconstruction using the PCI and non-PCI algorithms, respectively. In this case the non-PCI algorithm failed to reproduce the profile which was known from independent measurements to be a single bunch [22]. A comparison between the minimal phase and the PCI reconstructions shows no significant differences in the shape of the bunch or in its width. Nevertheless, the PCI profile is smoother, without large oscillations or unphysical (negative) values, which make it more realistic.

IV. SUMMARY AND CONCLUSIONS
We have developed a new phase retrieval algorithm based on repeated iterations between the time and frequency domains, with constraints applied in each domain. The two novel features are: (a) the use of the minimal phase as a lower limit for the missing phase and (b) the use of conditions (3) and (4) as criteria for the convergence of the iteration to the most probable and correct solution. We have applied the algorithm to the reconstruction of a number of complex synthetic profiles, some of which are shown in Figs. 2-6 and, also, to a set of experimental data (Fig. 7). In all of the cases that we have studied, the new algorithm (PCI) gave at least as good and usually better results than the non-PCI algorithm. In certain cases, e.g., Figs. 2 and 5, the minimal phase method would have been adequate for an accurate reconstruction but in both cases the PCI algorithm provided answers that were close to those derived from the Kramers-Kronig method. Since the existence or otherwise, of a significant Blaschke phase contribution cannot be deduced experimentally or known a priori, the PCI algorithm provides an acceptable and reliable reconstruction which does not depend on any initial assumptions. It is worth noting that the PCI algorithm has no predefined number of iterations and the reconstruction is assumed to be complete when both conditions (3) and (4) are satisfied simultaneously. One clear and significant advantage of the PCI algorithm over other similar procedures is the reproducibility of the results: there are no random seed phases at the start of the iterations and the answer returned by the PCI will be the same, every time it is run, provided that the conditions are satisfied simultaneously. We believe that the technique presented in this paper should be of interest to other experiments involving 1D reconstructions based on spectral yield measurements.