Study of Phase Reconstruction Techniques applied to Smith-Purcell Radiation Measurements
Nicolas Delerue, Joanna Barros, Oleg Bezshyyko, Vitalii Khodnevych
aa r X i v : . [ phy s i c s . acc - ph ] D ec STUDY OF PHASE RECONSTRUCTION TECHNIQUES APPLIED TOSMITH-PURCELL RADIATION MEASUREMENTS
N. Delerue ∗ , J. Barros, LAL, Orsay, France PARIS SACLAYO. Bezshyyko, V. Khodnevych, Taras Shevchenko National University of Kyiv, Ukraine Abstract
Measurements of coherent radiation at accelerators typ-ically give the absolute value of the beam profile Fouriertransform but not its phase. Phase reconstruction tech-niques such as Hilbert transform or Kramers Kronig recon-struction are used to recover such phase. We report a studyof the performances of these methods and how to optimizethe reconstructed profiles.
LONGITUDINAL BUNCH PROFILEMEASUREMENT AT PARTICLEACCELERATORS
On a particle accelerator the longitudinal profiles of a par-ticle bunch can not easily be measured. Several indirectmeasurement techniques have been established relying onthe measurement of the spectrum of radiation emitted bythe bunch either when it crosses a different material [1] orwhen it passes near a different material [2, 3]. This emittedspectrum encode the longitudinal profile through the rela-tion: I ( λ ) = I ( λ )( N + | F ( λ ) | N ) (1)where I ( λ ) is the emitted intensity as a function of thewavelength λ . I ( λ ) is the intensity of the signal emittedby a single particle and F ( λ ) is a form factor that encodesthe longitudinal and transverse shape of the particle bunch.Recovering the longitudinal profile requires to invert thisequation however this is not straightforward as the informa-tion about the phase of the form factor can not be measuredand therefore is not available.A phase reconstruction algorithm must therefore be usedto recover this phase. Several methods exist (see for exam-ple [4]). In this article we describe how we implementedtwo of these methods and compared their performances. RECONSTRUCTION METHODS
When it is only possible to measure the amplitude of thecomplex signal, it is necessary to recover the phase of theavailable data. We assume that the function of the longitu-dinal beam density is analytical. For an analytic functionthis is easier because the real and imaginary part are notcompletely independent. The Kramers-Kronig relations [4]helps restore the imaginary part of an analytic function ε ( ω ) from its real part and vice versa.To recover the phase from the amplitude, the functionshould be written as: log ( F ( ω )) = log ( ρ ( ω )) + i Θ( ω ) ∗ [email protected] with ρ ( ω ) its amplitude and Θ( ω ) its phase. The Kramers-Kronig relations can then be applied as follows: Θ( ω ) = 2 ω π P Z + ∞ ln ( ρ ( ω )) ω − ω dω (2)The basis of this relationship are the Cauchy-Riemann con-ditions (analyticity of function).In some cases this phase can also be obtained simply byusing the Hilbert transform of the spectrum: Θ( ω ) = − π P Z + ∞−∞ ln ( ρ ( ω )) ω − ω dω. (3)As the Hilbert transform ( H ) is related to the Fourier trans-form ( F ): F ( H ( u ))( ω ) = ( − isgn ( ω )) F ( u )( ω ) , (4)the calculation of phase can use an optimised FFT codeand is therefore much faster than calculating the Kramers-Kronig’s integral. We implemented in Matlab these two dif-ferent phase reconstruction methods. The Hilbert transformmethod has the advantage of being directly implemented inMatlab, allowing a much faster computing. DESCRIPTION OF THE SIMULATIONS
To test the performance of these methods we created asmall Monte-Carlo program that randomly simulates pro-files ( G ( x ) ) made of the combination of 5 gaussians accord-ing to the formula G ( x ) = P i =1 A i exp − ( xmX − µ i ) σ i where mX = 2 and A i , µ i and σ i are random numbers with x ∈ [1; mX ] , A i ∈ [0; 1] , µ i ∈ . − .
44; +11 . × − and σ i ∈ [3; 9] × − . The values of these rangeshave been chosen to generate profiles that are not discon-nected (that is profiles whose intensity drops to almost zerobetween two peaks) without being perfect gaussian. Wechecked that our conclusions are valid across this range.Using this formula we generated 1000 profiles, thentook the absolute value of their Fourier transform F = k FFT ( G ) k and sampled at a limited number of frequencypoints ( F i = F ( ω i ) ) as would be done with a real experi-ment in which the number of measurement points is limited(limited number of detectors or limited number of scanningsteps).To estimate the performance of the reconstruction severalestimators are available. We choose to use the χ , definedas follow: χ = X i ω i ( O i − E i ) /N, (5)here O i is the observed value , E i is the expected (simu-lated) value, ω i = 1 / √ 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 χ than a profile with oscillations (see fig-ure 1). This can be partly mitigated (in the case of horizon-tal offset) by offsetting one profile with respect to the otheruntil the χ is minimized.Also we decided to look at the FWHM which was gen-eralized as FWXM where X ∈ [0 .
1; 0 . is the fraction ofthe maximum value at which the full width of the recon-structed profile was calculated (with this definition the stan-dard FWHM is noted FW0.5M). We created an estimator ∆ F W XM defined as follow: ∆ F W XM = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) F W XM orig − F W XM reco
F W XM orig (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) X ∈{ . . . . . } (6)where F W XM orig and
F W XM reco are the FWXM ofthe original and reconstructed profiles respectively.Here two profiles that are similar but slightly offset (inposition or amplitude) will nevertheless return good valuesof this estimator despite returning a rather large χ . N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalSine noiseOffset
Figure 1: Example of profiles giving very different χ despite being relatively similar. χ sine noise = 3 . × − , χ offset = 7 . × − ;For profile with sine noise: FW0.1M=0.0241,FW0.2M=0.044 FWHM=0.0621 FW0.8M=0.1849FW0.9M=0.3619. As FWXM calculated from top ofprofile, for all profiles FWXM=0.To ensure that the choice of the parameters σ i and µ i forthe simulations does not bias significantly the results, theirvalue has been varied and this is shown in figure 2.Different distributions have been used for the frequen-cies ω i : linear, logarithmic, triple-sine. In most samplingschemes we used 33 frequencies to make it comparable withthe Triple-sine distribution used in [5]. These samplingschemes are defined as follow:• Triple-sine
This sampling matches that of the E-203 ex-periment at FACET [5]. Eleven detectors are located −7 −6 −5 µ / µ init χ Kramers−KronigHilbert10 µ / µ init ∆ F W H M Kramers−KronigHilbert10 −1 −8 −7 −6 −5 −4 σ / σ init χ Kramers−KronigHilbert10 −1 σ / σ init ∆ F W H M Kramers−KronigHilbert
Figure 2: The constraints on the parameters σ i (top) and µ i (bottom) due to effect of scaling in terms of the χ and Delta
FWHM ratio. For each point 1000 simulations weremade.very o around the interaction point and 3 differentsets of wavelengths are used, giving the following dis-tribution: cω i = l n (1 − cos (Θ i )) (7)with l n = 50 , , µm and Θ i varying between o and o by steps of o .• Linear sampling
Here sampling points are distributeduniformly. The first and last points of sampling are thefirst ( ω ) and last ( ω f ) points used in the Triple-sinesampling. The following formula gives the samplingfrequencies: ω i = ω + ( ω f − ω ) / × [0 : 1 : 32] . (8)• Logarithmic sampling . Here sampling points are dis-tributed logarithmically: ω ∗ exp ( log ( ω f /ω ) × [0 : 1 : 32] / . (9)For this sampling also the first and last points are thesame as in Triple-sine sampling.The study of the sampling is important, as it shows the bestposition of the detectors and also how to optimize the sys-tem. Linearly sampled spectrum gives the best result asshown in figure 3. Triple sine Linear Logarithmic00.20.40.60.811.21.4 x 10 −6 χ HilbertKramers−KronigTriple sine Linear Logarithmic00.010.020.030.040.050.060.070.08 ∆ F W H M HilbertKramers−Kronig
Figure 3: Comparison of different samplings with χ cri-terium (top) and ∆ F W HM (bottom).However, the linear sampling may not be practical to re-alize in a the real world. One needs to take into account the spatial size of the detectors (about 10 degrees in the case ofE-203) and there is also a limit on the start and end pointsof detectors location (35-145 degrees for E-203). So linearsampling at a wide range of frequencies is impossible withthis number of points. An investigation of how many linearsampling points can be used for a given angle difference be-tween detectors shows that such physical constraints reducestrongly the number of detectors that can be used. Figure 4shows examples of detector positions. The position of thered points is calculated using formula 7 and blue are possi-ble detector positions which don’t break the minimum de-tector distance (MDD) given on top of each plot.
MDD=5 o Triple sine (Ts)Linear (Ls)MaxLinear(Lsmx)MDD=10 o Triple sine (Ts)Linear (Ls)MaxLinear(Lsmx)
Figure 4: Detector position for linear sampling for a mini-mal detector distance of o (top) and o (bottom) .Figure 5 shows a comparison of the performancesachieved with such positioning for different MDD. In eachcase the triple sine sampling (Ts) is better than the linearsampling (Ls) and close from the maximum number of de-tectors with linear sampling (Lsmx). So the Ts configura-tion is favored and will be used in the rest of this paper. Thecomparison between Ts1, Ts5 and Ts10 shows that recon-struction performances are limited by the MDD.The choice of 33 frequencies for the sampling of the spec-trum was made to match the current layout used on E-203.However it is important to check if there is an optimumvalue. To perform this check we used the same simulationsand the same simulated spectrum but sampled with 3 to 140points. The effect of changing the sampling frequencies onthe χ is shown in figure 6. This study uses Triple sine s1 Lsmx1 Ts1 Ls5 Lsmx5 Ts5 Ls10 Lsmx10 Ts1000.020.040.060.080.10.12 ∆ F W H M Comparison of different sampling methods, Hilbert methodLs1 Lsmx1 Ts1 Ls5 Lsmx5 Ts5 Ls10 Lsmx10 Ts100123456 x 10 −7 χ Comparison of different sampling methods, Hilbert method
Figure 5: Comparison of different sampling with number ofMDD with χ criterium (top) and ∆ F W HM (bottom). Ls islinear sampling with o , o , MDD and Ts is Triple sinesampling; mx mean the reconstruction use the maximumnumber of detectors (blue and red dots in figure 4).sampling with 1000 profiles for each point and both recon-struction method.It can be seen at figure 6 that beyond about 33 samplingpoints the gain on the reconstructed χ is marginal.After applying the sampling procedure the data need tobe interpolated and extrapolated to have a larger number ofpoints in the spectrum. Interpolation is done using Piece-wise Cubic Hermite Interpolating Polynomial (PCHIP) [8],as suggested in [7]. The interpolation function must sat-isfy the following criteria: it must conserve the slope at thetwo endpoints (to have a continuous derivative) and respectsmonotonicity. PCHIP interpolation has been chosen as itmatches these requirements.For low frequency extrapolation two methods have beeninvestigated: Gaussian or Taylorian.In the Gaussian method, we defined the extrapolation asfollow: ρ LF ( ω ) = Ae − ( ω − B ) / C (10)Where ρ HF ( ω ) is the extrapolated spectrum at low fre-quency and the constants A, B, and C were chosen fromthe following conditions:• ρ LF (0) = 1 • ρ LF ( ω ) = ρ ( ω ) • ρ ′ LF ( ω ) = ρ ′ ( ω ) The extrapolation relies in the fact that according to thecentral limit theorem in the time space the expected profile −6 Number of detectors χ Kramers−KronigHilbert0 50 100 150 20000.020.040.060.080.10.12 Number of detectors ∆ F W H M Kramers−KronigHilbert
Figure 6: Effect of the sampling frequencies on the χ (top)and ∆ F W HM (bottom).is Gaussian-like and in the frequency space it will also beGaussian.The other extrapolation method is based on Taylor expan-sion with the following definition: F ( ω ) = Z ∞ dtS ( t ) e − i ( ωt ) = Z ∞ dtS ( t ) ∞ X k =0 ( − iωt ) k k ! == ∞ X k =0 (cid:18) ( − iω ) k k ! Z ∞ dtS ( t ) t k (cid:19) = ∞ X k =0 (cid:18) ( − iω ) k k ! < t k > (cid:19) (11)Approximation to the 4th order gives the following lowfrequency (LF) extrapolation: ρ LF = | F ( ω ) | = p A + Bω + Cω (12)Conditions for the constants A, B and C are the same.Comparison of different LF extrapolation can be found infigure 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 werealso tested. The most common [7, 10] is : ρ HF ( ω ) = Aω − , (13)where ρ HF ( ω ) is the extrapolated spectrum at high fre-quency and A = ρ f ω f , where ρ f is spectrum value of finalpoint ω f . −3 Time, [ps] N o r m a li z ed i n t en s i t y , a . u . Profile N2
Real profileGaussTeylorReal LF spectrum 10 −1 I n t en s i t y , a . u . Profile N2
Real spectrumSampled pointsGaussTeylorReal LF spectrum
Figure 7: Comparison of different LF extrapolation: exam-ple of spectrum (top) and profile (bottom) and histogramwith mean χ for each method (bottom). Gaussian and Tay-lorian methods are described in the text. "Real LF spec-trum" means that the real LF spectrum is used. For thissimulation the Hilbert method of phase recovery and Aω B high frequency extrapolation were used.The second method uses the same consideration as in Laiand Sievers [9]: Assuming that the bunch size is finite withtwo end points at z = 0 and at z = σ z then the longitudinalcharge distribution ( S ) must follow S (0) = S ( σ z ) = 0 . Anintegration by parts gives : F ( ω ) = Z ∞ dzS ( z ) e i ( ωc ) z == S ( z ) i ωc e i ( ωc ) z (cid:12)(cid:12)(cid:12) σ z − S ′ ( z ) (cid:16) i ωc (cid:17) e i ( ωc ) z (cid:12)(cid:12)(cid:12) σ z + . . . (14)The first term vanishes because of the boundary conditions,so for large ω , F ( ω ) is proportional to ω − and two condi-tions have to be matched :• ρ HF ( ω fmax ) = ρ ( ω fmax ) • ρ ′ HF ( ω fmax ) = ρ ′ ( ω fmax ) where ω fmax is the last sampled point of the spectrum. Tosatisfy the boundary condition two constants are needed,giving a two-terms extrapolation : ρ HF ( ω ) = Aω − + Bω − (15) Gauss Teylor Real0.05750.05750.05750.05750.05750.05750.05750.05750.0575 ∆ F W H M Gauss Teylor Real3.2313.23153.2323.23253.2333.23353.234 x 10 −7 χ Figure 8: Comparison of different LF extrapolation: his-togram with mean χ for each method (top) and ∆ F W HM (bottom).or extrapolation with degree of frequency as free parameter: ρ HF ( ω ) = Aω B (16)where the A and B – coefficients which are calculatedfrom the last data samples and the boundary conditions asfollow:• B = ρ ′ HF ( ω fmax ) ω fmax /ρ ( ω fmax ) • A = ρ ( ω fmax ) /ω Bfmax
The requirement of finite bunch size requires B ≤ , soin the case where the fit gives B > − we use B = − .Two other extrapolation methods also have been investi-gated:• ρ HF ( ω f ) = 0 for ω f > ω fmax • ρ HF ( ω f ) = ρ real ( ω f ) for ω f > ω fmax where ρ real is the real spectrum.These HF extrapolation methods are compared in figure 9and 10.Thus, by virtue of the above arguments and simulations,it’s naturally to choose the high-frequency extrapolation bypower function. STUDY OF THE RECONSTRUCTIONPERFORMANCE
After applying extrapolation and interpolation, the spec-trum recovery is completed. Then we used different re-construction techniques to reconstruct the original profile. N o r m a li z ed i n t en s i t y , a . u . Profile N34
Real profileAw −2 +Bw −3 ZeroAw B Aw −4 Real HF spectrum
20 30 40 50 60 7010 −4 −3 −2 Frequency, [THz] I n t en s i t y , a . u . Profile N34
Real profileSampled spectrumAw −2 +Bw −3 ZeroAw B Aw −4 Real HF spectrum
Figure 9: Comparison of different HF extrapolations : ex-ample of spectrum (top) and profile (bottom). For these sim-ulations the Hilbert reconstruction method of phase recov-ery and Gaussian LF extrapolations were used.
Zero Aw^−4Aw^−2+Bw^−3 Real Aw^B012345 x 10 −7 χ GaussiansZero Aw^−4Aw^−2+Bw^−3 Real Aw^B00.020.040.060.080.1 ∆ F W H M Gaussians
Figure 10: Comparison of different HF extrapolation forGaussian : histogram with mean χ (top) and ∆ F W HM (bottom). For each reconstruction method some profiles are very wellreconstructed whereas some other are not so well recon-structed. Examples of well reconstructed profiles are shownin figure 11 and examples of poorly reconstructed profile areshown in figure 12. −3 Profile N541 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.8 40.9 41 41.100.0020.0040.0060.0080.010.012 Profile N658 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.8 40.9 41 41.100.0020.0040.0060.0080.010.012 Profile N914 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig
Figure 11: Examples of well reconstructed profile. Theoriginal profile is in blue and the profiles reconstructed withthe Hilbert transform and the full Kramers-Kronig proce-dures are in red and black respectively.The ∆ F W XM and χ distribution of the 1000 simula-tions which were made and then reconstructed using theHilbert transform method and Kramers-Kornig reconstruc-tion are shown in figure 13. There is a good concordancein FWHM between two methods indicating that they areboth good at finding the bunch length. However, the Hilbertmethod gives lower χ indicating that this method is betterat reconstruction of the bunch profile.The fact that the phase recovery method based on theKramers-Kronig relation gives a worst χ than the method −3 Profile N227 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.8 40.9 41 41.100.0020.0040.0060.0080.01 Profile N231 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.8 40.9 41 41.100.0020.0040.0060.0080.01 Profile N667 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig
Figure 12: Example of poorly reconstructed profile. Theoriginal profile is in blue and the profiles reconstructed withthe Hilbert transform and the full Kramers-Kronig proce-dures are in red and black respectively.based on Hilbert relation has been investigated. It is causedby the presence of negative components in the tails of theprofiles. Figure 14 highlights this issue for one of the pro-files. Figure 15 shows the different FWXH values for dif-ferent values of X. This shows that at different height of theprofiles the quality of reconstruction varies: there is a betteragreement in the tails (X=10%) than at the top of the profile(X=90%). Figure 16 shows the modulus of the differencebetween the original and reconstructed profiles. One cansee oscillations in the difference between the original andreconstructed profile.While doing this work we also became aware of the dis-cussion in [6] where it is argued that these reconstructionmethod have more difficulties with lorentzian profiles thangaussian profiles. Therefore we simulated 1000 Lorenzian ∆ FWHM N u m be r o f p r o f il e s ∆ FWHM distribution
HilbertKramers−Kronig0 0.5 1 1.5 2 2.5 3x 10 −6 χ N u m be r o f p r o f il e s χ distribution HilbertKramers−Kronig
Figure 13: ∆ F W HM (top) and χ (bottom) distributionof 1000 simulations reconstructed using the Hilbert trans-form method (black line) and Kramers-Kronig reconstruc-tion method (red line). −3 Profile N541 N o r m a li z ed i n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig 40 40.2 40.4 40.6 40.8−2024x 10 −5 Area 140.95 40.96 40.97 40.98 40.999.49.6x 10 −3 Area 241 41.2 41.4 41.6 41.8 42−4−2024x 10 −5 Area 3
Figure 14: Example of reconstructed profile with zoomson the peak and tails. One can see that the profile recon-structed using the Kramers-Kronig method has a negativecomponent. This will dominate the final χ and explainswhy the χ obtained by this method is higher as shown infigure 13. W0.1MFW0.2M FWHM FW0.8MFW0.9M00.050.10.150.20.25 ∆ F W X M HilbertKramers−Kronig
Figure 15: ∆ F W XM for 1000 profiles with both methods. I n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.8 40.9 41 41.1024 x 10 −3 I n t en s i t y , a . u . Time, [ps]
Hilbert abs. diff.Kramers−Kronig abs. diff.40.7 40.75 40.8 40.85 40.9 40.95 4100.0050.01 Profile N914 I n t en s i t y , a . u . Time, [ps]
OriginalHilbertKramers−Kronig40.7 40.75 40.8 40.85 40.9 40.95 41024 x 10 −4 I n t en s i t y , a . u . Time, [ps]
Hilbert abs. diff.Kramers−Kronig abs. diff.
Figure 16: Original and reconstructed profile and their dif-ference for bad profile (top) and good profile (bottom).profiles and performed a similar study. This is shown in fig-ure 17. Although the χ is slightly worse in that case thanin the case of gaussian profiles there still a good agreementbetween the original and reconstructed profiles.In our discussion so far we considered only the ideal casewhere no noise is added to the measured spectrum. How-ever in a real experiment a noise component has to be addedto the measured spectrum. This noise was added as follow: O ′ i = O i × [1 + ( n i N max )] (17)where O i is the observed value, O ′ i is the observed valuewith noise, n i is a random number between 0 and 1 (all num-bers between 0 and 1 being equiprobable), and N max is themaximum noise for that simulation (depending on the casethis can be 5%, 10%, 20%, 30%, 40% or 50%). This study −6 χ N u m be r o f p r o f il e s Kramers−KronigHilbert0.5 1 1.5 2 2.5020406080100120 ∆ FWHM N u m be r o f p r o f il e s Kramers−KronigHilbert
Figure 17: Distribution of the χ in the case of a Lorenziandistribution.was done using linear sampling with 33 samples and 1000simulated profiles for each noise value. Figure 18 showshow the χ is modified when this noise component is added. DISCUSSION
We performed extensive simulation to estimate the per-formance of two phase recovery methods in the case ofmulti-gaussian and Lorenzian profiles. In both cases wefound that when the sampling frequencies are chosen cor-rectly we obtained a good agreement between the originaland reconstructed profiles (in most cases ∆ F W XM < ; χ ∼ − ). This confirms that such methods are suitableto reconstruct the longitudinal profiles measured at particleaccelerators using radiative methods.The authors are grateful for the funding received from theFrench ANR (contract ANR-12-JS05-0003-01), the PICS(CNRS) "Development of the instrumentation for accelera-tor experiments, beam monitoring and other applications"and Research Grant
10 20 30 40 5000.511.522.533.5 x 10 −5 Percent of noise in amplitude χ Kramers−KronigHilbert0 10 20 30 40 5000.020.040.060.08 Percent of noise in amplitude ∆ F W H M Kramers−KronigHilbert
Figure 18: Mean χ and ∆ F W XM as function of noise am-plitude.
EFERENCES [1] M.-A Tordeux and J. Papadacci. A new OTR based beamemittance monitor for the linac of LURE. EPAC 2000.[2] A. Cianchi, et al. Non-intercepting diagnostic for high bright-ness electron beams using optical diffraction radiation in-terference (odri).
Journal of Physics: Conference Series ,357(1):012019, 2012.[3] V. Blackmore et al. First measurements of the longitudi-nal bunch profile of a 28.5 GeV beam using coherent Smith-Purcell radiation.
Phys. Rev. ST Accel. Beams , 12:032803,Mar 2009.[4] O. Grimm and P. Schmüser. Principles of longitudinalbeam diagnostics with coherent radiation.
TESLA FEL note ,page 03, 2006.[5] L. Andrews et al. Reconstruction of the time profile of20.35 GeV, subpicosecond long electron bunches by meansof Coherent Smith-Purcell radiation.