Study of Phase Reconstruction Techniques applied to Smith-Purcell Radiation Measurements

On a particle accelerator the longitudinal profiles of a particle bunch can not easily be measured. Several indirect measurement techniques have been established relying on the measurement of the spectrum of radiation emitted by the bunch either when it crosses a different material [1] or when it passes near a different material [2, 3]. This emitted spectrum encode the longitudinal profile through the relation:

where $I(\lambda)$ is the emitted intensity as a function of the wavelength $\lambda$. $I_{1}(\lambda)$ is the intensity of the signal emitted by a single particle and $F(\lambda)$ is a form factor that encodes the longitudinal and transverse shape of the particle bunch. Recovering the longitudinal profile requires to invert this equation however this is not straightforward as the information about the phase of the form factor can not be measured and therefore is not available.

A phase reconstruction algorithm must therefore be used to recover this phase. Several methods exist (see for example [4]). In this article we describe how we implemented two of these methods and compared their performances.

When it is only possible to measure the amplitude of the complex signal, it is necessary to recover the phase of the available data. We assume that the function of the longitudinal beam density is analytical. For an analytic function this is easier because the real and imaginary part are not completely independent. The Kramers-Kronig relations [4] helps restore the imaginary part of an analytic function $\varepsilon(\omega)$ from its real part and vice versa.

To recover the phase from the amplitude, the function should be written as: $log(F(\omega))=log(\rho(\omega))+i\Theta(\omega)$ with $\rho(\omega)$ its amplitude and $\Theta(\omega)$ its phase. The Kramers-Kronig relations can then be applied as follows:

The basis of this relationship are the Cauchy-Riemann conditions (analyticity of function).

In some cases this phase can also be obtained simply by using the Hilbert transform of the spectrum:

As the Hilbert transform (H) is related to the Fourier transform (${\cal F}$):

the calculation of phase can use an optimised FFT code and is therefore much faster than calculating the Kramers-Kronig’s integral. We implemented in Matlab these two different phase reconstruction methods. The Hilbert transform method has the advantage of being directly implemented in Matlab, allowing a much faster computing.

To test the performance of these methods we created a small Monte-Carlo program that randomly simulates profiles (${\cal G}(x)$) made of the combination of 5 gaussians according to the formula ${\cal G}(x)=\sum_{i=1}^{5}A_{i}\exp{\frac{-(\frac{x}{mX}-\mu_{i})^{2}}{2\sigma^{2}_{i}}}$ where $mX=2^{16}$ and $A_{i}$, $\mu_{i}$ and $\sigma_{i}$ are random numbers with $x\in[1;mX]$, $A_{i}\in[0;1]$, $\mu_{i}\in 0.5+[-11.44;+11.44]\times 10^{-9}$ and $\sigma_{i}\in[3;9]\times 10^{-9}$ . The values of these ranges have been chosen to generate profiles that are not disconnected (that is profiles whose intensity drops to almost zero between two peaks) without being perfect gaussian. We checked that our conclusions are valid across this range.

Using this formula we generated 1000 profiles, then took the absolute value of their Fourier transform ${\cal F}=\|\mbox{FFT}\left({\cal G}\right)\|$ and sampled at a limited number of frequency points ($F_{i}={\cal F}(\omega_{i})$) as would be done with a real experiment in which the number of measurement points is limited (limited number of detectors or limited number of scanning steps).

To estimate the performance of the reconstruction several estimators are available. We choose to use the $\chi^{2}$, defined as follow:

where $O_{i}$ is the observed value , $E_{i}$ is the expected (simulated) value, $\omega_{i}=1/\sqrt{O_{i}+E_{i}}$ is the weight of the point, N is the number of points.

However two very similar profiles but with a slight offset, will give a worse $\chi^{2}$ than a profile with oscillations (see figure 1). This can be partly mitigated (in the case of horizontal offset) by offsetting one profile with respect to the other until the $\chi^{2}$ is minimized.

Also we decided to look at the FWHM which was generalized as FWXM where $X\in[0.1;0.9]$ is the fraction of the maximum value at which the full width of the reconstructed profile was calculated (with this definition the standard FWHM is noted FW0.5M). We created an estimator $\Delta_{FWXM}$ defined as follow:

where $FWXM_{\mbox{orig}}$ and $FWXM_{\mbox{reco}}$ are the FWXM of the original and reconstructed profiles respectively.

Here two profiles that are similar but slightly offset (in position or amplitude) will nevertheless return good values of this estimator despite returning a rather large $\chi^{2}$.

To ensure that the choice of the parameters $\sigma_{i}$ and $\mu_{i}$ for the simulations does not bias significantly the results, their value has been varied and this is shown in figure 2.

Different distributions have been used for the frequencies $\omega_{i}$: linear, logarithmic, triple-sine. In most sampling schemes we used 33 frequencies to make it comparable with the Triple-sine distribution used in [5]. These sampling schemes are defined as follow:

• •

  Triple-sine This sampling matches that of the E-203 experiment at FACET [5]. Eleven detectors are located every $10^{o}$ around the interaction point and 3 different sets of wavelengths are used, giving the following distribution:

  $$\frac{c}{\omega_{i}}=l_{n}(1-cos(\Theta_{i}))$$

  (7)

  with $l_{n}=50,250,1500\mu m$ and $\Theta_{i}$ varying between $40^{o}$ and $140^{o}$ by steps of $10^{o}$.

• •

  Linear sampling Here sampling points are distributed uniformly. The first and last points of sampling are the first ($\omega_{0}$) and last ($\omega_{f}$) points used in the Triple-sine sampling. The following formula gives the sampling frequencies:

  $$\omega_{i}=\omega_{0}+(\omega_{f}-\omega_{0})/32\times[0:1:32].$$

  (8)

• •

  Logarithmic sampling. Here sampling points are distributed logarithmically:

  $$\omega_{0}*exp(log(\omega_{f}/\omega_{0})\times[0:1:32]/32).$$

  (9)

  For this sampling also the first and last points are the same as in Triple-sine sampling.

The study of the sampling is important, as it shows the best position of the detectors and also how to optimize the system. Linearly sampled spectrum gives the best result as shown in figure 3.

However, the linear sampling may not be practical to realize in a the real world. One needs to take into account the spatial size of the detectors (about 10 degrees in the case of E-203) and there is also a limit on the start and end points of detectors location (35-145 degrees for E-203). So linear sampling at a wide range of frequencies is impossible with this number of points. An investigation of how many linear sampling points can be used for a given angle difference between detectors shows that such physical constraints reduce strongly the number of detectors that can be used. Figure 4 shows examples of detector positions. The position of the red points is calculated using formula 7 and blue are possible detector positions which don’t break the minimum detector distance (MDD) given on top of each plot.

Figure 5 shows a comparison of the performances achieved with such positioning for different MDD. In each case the triple sine sampling (Ts) is better than the linear sampling (Ls) and close from the maximum number of detectors with linear sampling (Lsmx). So the Ts configuration is favored and will be used in the rest of this paper. The comparison between Ts1, Ts5 and Ts10 shows that reconstruction performances are limited by the MDD.

The choice of 33 frequencies for the sampling of the spectrum was made to match the current layout used on E-203. However it is important to check if there is an optimum value. To perform this check we used the same simulations and the same simulated spectrum but sampled with 3 to 140 points. The effect of changing the sampling frequencies on the $\chi^{2}$ is shown in figure 6. This study uses Triple sine sampling with 1000 profiles for each point and both reconstruction method.

It can be seen at figure 6 that beyond about 33 sampling points the gain on the reconstructed $\chi^{2}$ is marginal.

After applying the sampling procedure the data need to be interpolated and extrapolated to have a larger number of points in the spectrum. Interpolation is done using Piecewise Cubic Hermite Interpolating Polynomial (PCHIP) [8], as suggested in [7]. The interpolation function must satisfy the following criteria: it must conserve the slope at the two endpoints (to have a continuous derivative) and respects monotonicity. PCHIP interpolation has been chosen as it matches these requirements.

For low frequency extrapolation two methods have been investigated: Gaussian or Taylorian.

In the Gaussian method, we defined the extrapolation as follow:

Where $\rho_{HF}(\omega)$ is the extrapolated spectrum at low frequency and the constants A, B, and C were chosen from the following conditions:

• •

  $\rho_{LF}(0)=1$

• •

  $\rho_{LF}(\omega_{0})=\rho(\omega_{0})$

• •

  $\rho_{LF}^{\prime}(\omega_{0})=\rho^{\prime}(\omega_{0})$

The extrapolation relies in the fact that according to the central limit theorem in the time space the expected profile is Gaussian-like and in the frequency space it will also be Gaussian.

The other extrapolation method is based on Taylor expansion with the following definition:

Approximation to the 4th order gives the following low frequency (LF) extrapolation:

Conditions for the constants A, B and C are the same. Comparison of different LF extrapolation can be found in figure 7 and the performances of these methods in figure 8.

In the rest of this paper we used the Gaussian method.

Several high frequency (HF) extrapolation methods were also tested. The most common [7, 10] is :

where $\rho_{HF}(\omega)$ is the extrapolated spectrum at high frequency and $A=\rho_{f}\omega_{f}^{4}$, where $\rho_{f}$ is spectrum value of final point $\omega_{f}$.

The second method uses the same consideration as in Lai and Sievers [9]: Assuming that the bunch size is finite with two end points at $z=0$ and at $z=\sigma_{z}$ then the longitudinal charge distribution ($S$) must follow $S(0)=S(\sigma_{z})=0$. An integration by parts gives :

The first term vanishes because of the boundary conditions, so for large $\omega$, $F(\omega)$ is proportional to $\omega^{-2}$ and two conditions have to be matched :

• •

  $\rho_{HF}(\omega_{fmax})=\rho(\omega_{fmax})$

• •

  $\rho_{HF}^{\prime}(\omega_{fmax})=\rho^{\prime}(\omega_{fmax})$

where $\omega_{fmax}$ is the last sampled point of the spectrum. To satisfy the boundary condition two constants are needed, giving a two-terms extrapolation :

or extrapolation with degree of frequency as free parameter:

where the A and B – coefficients which are calculated from the last data samples and the boundary conditions as follow:

• •

  $B={\rho_{HF}^{\prime}(\omega_{fmax})\omega_{fmax}}/{\rho(\omega_{fmax})}$

• •

  $A={\rho(\omega_{fmax})}/{\omega_{fmax}^{B}}$

The requirement of finite bunch size requires $B\leq 2$, so in the case where the fit gives $B>-2$ we use $B=-2$. Two other extrapolation methods also have been investigated:

• •

  $\rho_{HF}(\omega_{f})=0$ for $\omega_{f}>\omega_{fmax}$

• •

  $\rho_{HF}(\omega_{f})=\rho_{real}(\omega_{f})$ for $\omega_{f}>\omega_{fmax}$ where $\rho_{real}$ is the real spectrum.

These HF extrapolation methods are compared in figure 9 and 10.

Thus, by virtue of the above arguments and simulations, it’s naturally to choose the high-frequency extrapolation by power function.

After applying extrapolation and interpolation, the spectrum recovery is completed. Then we used different reconstruction techniques to reconstruct the original profile. For each reconstruction method some profiles are very well reconstructed whereas some other are not so well reconstructed. Examples of well reconstructed profiles are shown in figure 11 and examples of poorly reconstructed profile are shown in figure 12.

The $\Delta_{FWXM}$ and $\chi^{2}$ distribution of the 1000 simulations which were made and then reconstructed using the Hilbert transform method and Kramers-Kornig reconstruction are shown in figure 13. There is a good concordance in FWHM between two methods indicating that they are both good at finding the bunch length. However, the Hilbert method gives lower $\chi^{2}$ indicating that this method is better at reconstruction of the bunch profile.

The fact that the phase recovery method based on the Kramers-Kronig relation gives a worst $\chi^{2}$ than the method based on Hilbert relation has been investigated. It is caused by the presence of negative components in the tails of the profiles. Figure 14 highlights this issue for one of the profiles.

Figure 15 shows the different FWXH values for different values of X. This shows that at different height of the profiles the quality of reconstruction varies: there is a better agreement in the tails (X=10%) than at the top of the profile (X=90%). Figure 16 shows the modulus of the difference between the original and reconstructed profiles. One can see oscillations in the difference between the original and reconstructed profile.

While doing this work we also became aware of the discussion in [6] where it is argued that these reconstruction method have more difficulties with lorentzian profiles than gaussian profiles. Therefore we simulated 1000 Lorenzian profiles and performed a similar study. This is shown in figure 17. Although the $\chi^{2}$ is slightly worse in that case than in the case of gaussian profiles there still a good agreement between the original and reconstructed profiles.

In our discussion so far we considered only the ideal case where no noise is added to the measured spectrum. However in a real experiment a noise component has to be added to the measured spectrum. This noise was added as follow :

where $O_{i}$ is the observed value, $O_{i}^{\prime}$ is the observed value with noise, $n_{i}$ is a random number between 0 and 1 (all numbers between 0 and 1 being equiprobable), and $N_{max}$ is the maximum noise for that simulation (depending on the case this can be 5%, 10%, 20%, 30%, 40% or 50%). This study was done using linear sampling with 33 samples and 1000 simulated profiles for each noise value. Figure 18 shows how the $\chi^{2}$ is modified when this noise component is added.

We performed extensive simulation to estimate the performance of two phase recovery methods in the case of multi-gaussian and Lorenzian profiles. In both cases we found that when the sampling frequencies are chosen correctly we obtained a good agreement between the original and reconstructed profiles (in most cases $\Delta_{FWXM}<10\%$; $\chi^{2}\sim 10^{-6}$). This confirms that such methods are suitable to reconstruct the longitudinal profiles measured at particle accelerators using radiative methods.

The authors are grateful for the funding received from the French ANR (contract ANR-12-JS05-0003-01), the PICS (CNRS) "Development of the instrumentation for accelerator experiments, beam monitoring and other applications" and Research Grant #F58/380-2013 (project F58/04) from the State Fund for Fundamental Researches of Ukraine in the frame of the State key laboratory of high energy physics and the IDEATE International Associated Laboratory (LIA).