Revealing the dynamics of ultrarelativistic non-equilibrium many-electron systems with phase space tomography

The description of physical processes with many-particle systems is a key approach to the modeling of numerous physical systems. For example in storage rings, where ultrarelativistic particles are agglomerated in dense bunches, the modeling and measurement of their phase-space distribution is of paramount importance: at any time the phase-space distribution not only determines the complete space-time evolution but also provides fundamental performance characteristics for storage ring operation. Here, we demonstrate a non-destructive tomographic imaging technique for the 2D longitudinal phase-space distribution of ultrarelativistic electron bunches. For this purpose, we utilize a unique setup, which streams turn-by-turn near-field measurements of bunch profiles at MHz repetition rates. To demonstrate the feasibility of our method, we induce a non-equilibrium state and show that the phase-space distribution microstructuring as well as the phase-space distribution dynamics can be observed in great detail. Our approach offers a pathway to control ultrashort bunches and supports, as one example, the development of compact accelerators with low energy footprints.


I. INTRODUCTION
The modeling of systems composed of a large number of interacting particles are of utmost importance in physics [1] with applications reaching from molecular dynamics simulations [2], the description of traffic dynamics [3] up to the treatment of quantized many-particle systems [4].One of the idealized model systems is the many-electron system.In applications, many-electron systems can be affected by relativistic effects, quantized emission of photons and coherence, if they are spatially compressed.For example, relativistic free electrons are utilized in particle accelerators for high-energy physics experiments and the generation of synchrotron radiation.In the latter case, the fact that accelerated ultrarelativistic electrons produce broadband emission from the microwave to the hard x-ray range is harnessed in evermore advanced synchrotron light sources throughout the word for widespread applications in science, industry and medicine.For such applications, accelerator-based light sources need to provide stable emission.This is achieved by storage rings, where the particles circulate in pulsed structures consisting of electron bunches.If sufficiently compressed, these bunches emit coherent synchrotron radiation (CSR) at wavelengths corresponding to the Fourier transform of their longitudinal density profile [5][6][7].Compared to incoherent radiation, CSR can be several orders of magnitude higher in intensity, which makes it particularly interesting for user applications.Hence, to provide an ever-increasing brilliance at frequencies in the THz range and above [8], scientists aim to increase the electron density by compressing the bunch as short as possible.However, there is a trade-off to this optimization: above a certain density threshold the collective radiation field exerts significant forces on the particles leading to a buildup and self-amplification of substructures on the bunch profile.This formation of substructures of an electron bunch with the associated dynamics is called microbunching instability (MBI).As a consequence, irregular bursts of CSR are emitted, which can at short wavelengths be even more intense than the CSR during stable bunch formation at an electron density below the stability threshold [9].The MBI is the main limiting factor for the generation of stable and intense CSR in storage rings run in short bunch mode.The MBI is furthermore considered as a natural process occurring during the emission of solar flares [10].During the MBI the electron bunch might be regarded as a non-equilibrium thermodynamic system with a steady flow of energy exhibiting rich structural and dynamic selforganized patterns similar to other non-equilibrium thermodynamic systems, such as the Rayleigh-Bérnard convection, Turing instability or the Belousov-Zhabotinsky reaction [11].Here, a deep understanding of the physical laws determining the development of the electron bunch density bears the potential to further control or stabilize the occurrence of the CSR bursts, which then could provide a bright THz radiation source for user applications [12].The details of the dynamics during the MBI are, however, still unclear and subject of intense research efforts [13][14][15][16][17][18][19][20][21][22][23][24][25].With respect to accelerator diagnostics, measuring the FIG. 1.
The setup to measure charge density profiles of electron bunches in a storage ring with single-shot electrooptical sampling of the near-field is shown (see [32][33][34]).The closed metallic beam pipe of the accelerator is here indicated as a transparent cylinder.
phase space distribution (PSD) of the electron bunch is key to a complete characterization of its physical state.While on the theoretical side, phase space dynamics can be simulated with great detail (for example by numerically solving the corresponding Vlasov-Fokker-Planck equation [17,20]), for non-destructive experiments with electron bunches the PSD remains an elusive quantity.Early ideas by Hancock et.al. [26] for the reconstruction of the longitudinal PSD of proton bunches and other concepts for heavy ion colliders [27], linear accelerators [28,29] as well as theoretical considerations [30], did not mature into a diagnostic tool for electron storage rings, where challenges from the high revolution frequencies of rings and fading out of information about substructure in the far-field have to be overcome.Hence, many interpretations of experimental results can only be taken as indirect conclusions about the state of the PSD of the electron bunches.In this paper, we describe how this diagnostics gap can be closed using single-shot electro-optical (EO) sampling of the electron bunch near-field in combination with turnby-turn measurements at MHz repetition rates using an ultra-fast line-array detection system [31].The longitudinal PSD not only determines the emission spectrum of the bunch, but also provides insights in other fundamental properties such as bunch length, energy spread or, more sophisticated, relativistic intra-bunch interactions.Throughout the paper we restrict our attention to the longitudinal PSD.We therefore consider the spatial coordinate in the direction of the bunch motion and approximate the bunch as one-dimensional line charge density.

II. RECONSTRUCTION METHOD
The working principle of our setup is visualized in Fig. 1 (see [34] for a detailed description): to measure the charge density of a relativistic bunch we overlap broadband linearly polarized chirped laser pulses from a fiber laser with the Coulomb field of the electron bunch in a gallium phosphide EO crystal.The timing is adjusted by a phase-locked loop synchronizing the laser repetition rate to the revolution frequency of the storage ring.Due to the relativistic motion of the bunch, the Coulomb field of an individual relativistic electron is highly compressed perpendicular to its motion in a radial pancake-like structure.Consequently, the time-dependent Coulomb field of the electron bunch leaked into the EO crystal is proportional to the charge density profile.Due to the chirp in the linearly polarized beam, the charge density profile is imprinted on the laser pulse as a wavelength dependent change of the polarization.In the next stage, this information is transformed into a wavelength dependent intensity profile by a polarizing cube beam splitter.In this process, the wavelength encodes the time-dependence.The profile is finally decoded by a grating and read out by an ultra-fast line array camera [31] in the spectral domain.By comparing the intensity profile of the laser pulse to a reference measurement, taken without the presence of electrons, we can directly deduce the longitudinal charge profile of the electron bunch [34].From the bunch profile measurements, we can reconstruct the PSD using the following approximation: during the time of half a synchrotron oscillation period (in our case ∼ 61 µs), we assume that the change of the microstructures is small.Then the resulting dynamics can be described by a well-known rotation of the PSD [19], which is a result of the phase focusing of the electron bunch to the radio frequency phase of the storage ring [36].In other words, we approximate the dynamics of the PSD by a rigid rotation with a rotation period of ∼ 122 µs.In Fig. 2 a, we show a typical configuration of the PSD taken from a simulation of the electron bunch dynamics based on the Vlasov-Fokker-Planck equation using Inovesa, a simulation program developed at the Karlsruhe Institute of Technology (KIT) [20].With respect to the phase space, an EO near-field measurement can be interpreted as a projection of PSD onto the axis related to a generalized coordinate (in our case we measure the arrival time).During consecutive measurements, performed on a turn-by-turn basis by our setup, the PSD of the electron bunch rotates, so that the experiment can be interpreted as a tomographic measurement of the electron bunch PSD.If we concatenate consecutive measurements taken during half a rotation period we obtain, in this interpretation, the Radon transform of the PSD resulting in the so-called sinogram (see Fig. 2 a bottom).The sinogram can be related to the 2D Fourier transform of the original image (or PSD) via the Fourier slice theorem [37].The reconstruction of the PSD is then performed with a filtered back-projection [37], a well-known procedure widely used for example in medicine for computed tomography scans.The approximation of the PSD dynamics as a rigid rota-  [20].We used a standard ramp filter to avoid the blurring effect of the back-projection [35] in combination with a high frequency cut-off at a relative frequency of 0.1 to reduce the exaggeration of edges.Left side: comparison between the original and reconstructed PSD.The "original" PSD is evaluated at the 90 • frame, while for the reconstruction the complete sinogram has to be taken into the account.Right side: several comparisons during different times of a bursting cycle.Even in situations with prominent substructures, where bunch self-interaction is expected to be strong, the reconstruction shows a very good agreement with the original PSD, even in fine details.
tion is an idealization made to motivate the reconstruction algorithm.This situation might be different especially above the stability threshold, where the substructures actually form because of deviations from the rigid rotation model [7].Hence, to validate our approach for the condition present during the experiments we use Inovesa simulations.The simulations provide a PSD during every time step, so that we can calculate the EO sinogram (which is a simple projection) and compare the corresponding reconstructed PSD to the original simulated PSD.In Fig. 2 b, we show that even though the PSD dynamics can be complicated during the MBI due to collective effects of the microstructures, the filtered back-projection can not only reconstruct the PSD on a coarse grained level, but also resolve subtleties of the microstructures (we marked a few of them with white FIG. 3. a) Experimental data: three revolution plots of 6000 consecutive measurements (corresponding to 18 synchrotron oscillations) of the electron bunch charge density taken from a longer data set with 500,000 measurements (∼ 0.18 s).For each revolution plot, the measurements are concatenated from left to right.The value for ∆t is determined by the pixel number of the line array, calibration measurements [34] and the mean bunch profile position.b) Corresponding reconstructions of the PSD for different time intervals.The numbers refer to the time intervals marked in the revolution plots on the left side, respectively.Upper side: reconstructions for different subtraction percentages of the mean bunch profile are compared for the first revolution plot.The dynamic cycle of the microbunching instability is repeated over and over again, as it can be seen by comparing the rows of the reconstructed PSDs for the different revolution plots.To provide the optimal contrast, we normalized the data for every image independently.

arrows).
To demonstrate this principle experimentally, we measured profiles of a single electron bunch during the MBI at the Karlsruhe Research Accelerator (KARA).During the experiments, the bunch contained roughly 2•10 9 electrons and was accelerated to a highly relativistic motion with a Lorentz factor γ ≈ 2500.To induce the MBI, we compressed the longitudinal bunch size with the low-α c mode of the storage ring to a few picoseconds.We then observed the development of the bunch profile at a repetition rate of 2.72 MHz corresponding to the revolution time of a bunch in the storage ring.In Fig. 3a, we show three different sections each consisting of 6000 successive profile measurements stacked from left to right.The sections are part of a larger data set with 500,000 turn-by-turn measurements.The representation of the sections in Fig. 3a is also referred to as "revolution plot" (see [34]).In each displayed revolution plot, the noise-induced [38] coherent synchrotron oscillation appears as a clearly visible center of mass oscillation of the bunch profile.We utilize this effect to determine the synchrotron oscillation period, which is needed to extract the sinograms for PSD reconstruction: if we consider a horizontal cross section in the revolution plots, the charge density fluctuates as a consequence of the synchrotron oscillation.Hence, the frequency of the synchrotron oscillation can be conveniently determined from the data set using a Fourier transform.To provide a precise estimation we used all bunch profiles of the data set and averaged the resulting spectrum along the horizontal cross sections.In the case shown here, the synchrotron frequency is estimated to be 8.28 kHz corresponding to 328 turns.Therefore, roughly 18 synchrotron oscillations are displayed for each revolution plot in Fig. 3 a.From the determination of the synchrotron oscillation frequency, we conclude that a PSD reconstruction can be performed from a sinogram consisting of about 328/2=164 consecutive measurements (the angle resolution is about 1.1 • ), i.e. a sinogram corresponds to a revo-FIG.4. The Euclidian distance between a reference PSD and each PSD from the subsequent data is shown (in blue).The PSD configurations at the minima and maxima as well as the reference PSD are displayed as insets.The Euclidian distance follows an oscillatory pattern as expected from a complete rigid rotation of a PSD structure with 7-fold symmetry.
lution plot with exact 164 measurements.In Fig. 3 b, we show reconstructions of the PSD during different times.
For the calculations we used a filtered back-projection with a ramp filter [35].Compared to the simulations, we have a bunch profile sampling rate that is approximately two times lower, so that we doubled the high frequency cut-off for the filtered back-projection from 0.1 to 0.2 of the highest frequency.The time intervals used for every reconstruction are indicated in each revolution plot in Fig. 3 a as white horizontal lines.Each of these short intervals on its own can be regarded as a sinogram of the PSD.In general, the observed microstructuring of the bunch profile is small compared to the overall profile.For the top revolution plot in Fig. 3, we show how the contrast of the PSD reconstructions can be enhanced: we calculate the mean bunch profile for each revolution plot and subtract a scaled version from every bunch profile before applying the filtered back-projection.Such a procedure might introduce artifacts and should be handled carefully.A comparison for different scaling factors is visible at top Fig. 3 b.We note that a 40 % subtraction can improve the contrast significantly without introducing additional features (deduced from a careful comparison with the PSD at 0 % subtraction) as it might be the case when subtracting 100 % of the mean bunch profile.At this point, we emphasize that the measured PSDs in Fig. 3 a are remarkably similar to the simulated ones in Fig. 2 b.In the next step, we use the measured PSDs, to analyze the development and the structure of the electron bunch distribution in phase-space.We first analyze the PSD reconstruction from the revolution plot in Fig. 3 a at the top: at the beginning of the measurement, substructures due to the MBI are not visible in the revolution plot at time interval 1.Hence, the PSD, shown in the first column of Fig. 3 b, corresponds approximately to a 2D Gaussian-like distribution.In the further development, diagonal stripes appear in the revolution plot.Here, the reconstruction of the PSD in Fig. 3 b at column 2 shows a star-like structuring with a symmetry center roughly located at the center of mass position.Afterwards, the structure of the PSD during time intervals 3 and 4 undergoes a complex dynamic development and the stripe pattern becomes irregular.Finally, during the time intervals 5 and 6, the stripe pattern becomes less visible and the PSD images correspond to Gaussian distributions overlaid with weak substructures.This relaxation during which the substructures disappear takes rather long in comparison to the fast appearance of the substructures in the beginning.The overall process shown here retraces a typical excitation and relaxation pattern during the MBI (see [39]).This is underlined by the reproduced dynamics demonstrated by the reconstructions for the two other revolution plots in Fig. 3 a.While this shows how the overall dynamics can be easily analyzed with the reconstructed PSDs, it is also possible to extract information about the dynamical development of structural features of the PSD: as an example, we calculated all PSDs during one complete synchrotron oscillation period starting at 18.74 ms of the measured data set (around the second step of the dynamic cycle in the first revolution plot in Fig. 3), where a star-like shape with a 7-fold symmetry is observed.In the next step, we investigate the correlation of every image with the reference frame obtained from the first 164 revolutions using a simple Euclidean distance [40,41].The result, shown in Fig. 4, features a well-defined oscillatory pattern.For the interpretation of the results it is important to note, that the reconstructed PSDs from time intervals up to half a synchrotron period after the reference frame are statistically dependent.This is due to the overlapping data used for the reconstructions (i.e. a sliding window is used for the turn-by-turn reconstructions).Therefore, the correlation should be investigated for more than half a synchrotron oscillation period, such as in Fig. 4, where the oscillatory pattern is recorded for a complete synchrotron oscillation period.While a distinct reconstruction provides structural information, the correlation of the PSD over time adds dynamic information, which is different from a selfcorrelation of a single PSD measurement.Here, the clear occurrence of seven minima indicates that the star-like structures of the PSD plots (in Fig. 4) have a 7-fold symmetry and might rotate with a rigid 360 • motion during one complete synchrotron oscillation period, similar to the model in [42].This situation might change for other times, different dynamic regimes of the MBI and different experimental settings of the storage ring.Since we did not enhance the contrast of the PSD reconstruction by mean bunch profile subtraction for this correlation study, this method demonstrates how situations with even chal-lenging signal-to-noise ratios can be handled.Finally, the Euclidean distance shows a slow but significant up and down bending, particularly visible at the progression of the minima positions.We attribute this effect to a mismatch between the rotation center and the center of gravity of the PSD (rather than an asymmetry of the PSD itself), which is consistent with the observation of the synchrotron oscillation in the revolution plots [38] (even when the substructures are faded out, see especially Fig. 3).

III. SUMMARY AND OUTLOOK
In this publication, we demonstrated a robust method to reconstruct the longitudinal PSD from turn-by-turn measurements of the longitudinal bunch profile using a filtered back-projection.In doing so we close an important diagnostics gap for electron storage rings.Our approach is based on the approximation of the electron bunch dynamics by a rigid rotation of the PSD during half a synchrotron oscillation period.From comparison to simulations, we provide empirical evidence, that the filtered back-projection performs very well, even under conditions with strong intra-bunch interactions.Applying the PSD reconstruction to experimental data, we are able to observe a typical dynamic cycle in a reproducible manner.
For future research, it might be interesting to analyze the filtered back-projection performance with respect to the PSD dynamics governed by the Vlasov-Fokker-Planck equation more deeply.Moreover, considering that the filtered back-projection is originally not intended to be used for PSD tomography of electron bunches, it might be possible to construct a specialized algorithm deduced from the Vlasov-Fokker-Planck equation relaxing the rigid rotation assumption.Finally, our approach offers a glimpse into the physics of equilibrium/non-equilibrium ultrarelativistic systems.It therefore supports the development of advanced accelerator concepts with low energy footprints, where coherent emission from compressed, high-density electron beams plays an import role [43].

A. Experimental parameters and setup
The measurements were enabled through an synchronized interplay between three sophisticated instruments: an accelerator providing strongly self-interacting relativistic electron bunches, a single-shot EO sampling experiment using a synchronized laser system and a fast line array camera system measuring at MHz repetition rates.The experiments were conducted at the Karlsruhe Research Accelerator (KARA), which was operated in a dedicated short bunch mode [44] to compress the bunch length to a few picoseconds.The machine parameter settings, which were similar to those in reference [25], are displayed in Table I.During the measurements, the beam energy was set to 1.3 GeV using a radio frequency voltage of 799 kV.The data of Fig. 3 was taken at a bunch current of 0.85 mA in single bunch operation mode at the so-called sawtooth bursting regime (see [25]).
For the EO measurement, we used chirped laser pulses from a custom-build regenerative amplifier [33] with a central wavelength of 1050 nm, an optical power of a few mW and a spectral width of about 80 nm.The ytterbium-doped fiber oscillator of the laser system operates at 62.5 MHz and is actively synchronized to the 500 MHz of the master oscillator.A pulse-picker is then used to reduce the repetition rate of the laser system (see also [34] for further details).Afterwards, the spectrum of the laser pulses was resolved by a grating and measured with a KIT-developed ultrafast spectrometer KALYPSO [31] (KArlsruhe Linear ar-raY detector for MHz-rePetition rate SpectrOscopy) using a 256-pixel back-illuminated silicon line array.The repetition rate of the line array readout was 2.72 MHz, so that that the electron bunch profile was detected on a turn-by-turn basis.

B. Details of the post-measurement data processing and the back-projection parameters
To calculate the bunch profile, we first determine the background signal of the line array by blocking the laser beam.In the next step, we record the signal of the laser pulses propagation through the EO (GaP) crystal without the overlap with the electron bunch.To do so, the phase of the laser synchronization system is changed, so that the laser pulse arrives before the electron bunch at the EO crystal.We subtract the background signal from the measurements with and without overlap and obtain the modulated and unmodulated signal.The result is shown in Fig. 5 a.For further calculations, we determined these offsets by averaging the intensity for all measurements in the region indicated by the intervals in Fig. 5 a. Afterwards, the offsets were removed from every data set.Figure 5 a displays in red and blue the corrected values for the modulated and unmodulated signals.Especially the kink is reduced significantly for these data sets.With the corrected data, we calculated the bunch profiles according to [34] by the division between the modulated and the unmodulated signal.In the next step we determined which pixel should be considered for further analysis.Clearly such a decision depends on the signal-to-noise ratio for the calculated charge density of the bunch profiles.Here we used the following criterion: we consider all pixels for which half of the maximum signal modulation can be distinguished from a background signal with an accuracy of 90 %.We assume that the maximum signal is about 0.2 (determined from average bunch profiles).As described in [34], let σ ρ,i be the standard deviation for a hypothetical repeated, single-shot measurement for pixel i, then we consider all pixel for which: is fulfilled.Figure 5 b displays the dependence of σ ρ on the pixel number.From this procedure the pixels numbers between 48 and 201 are selected.While (1) provides an objective and reproducible measure, which data should be considered for further evaluation, the reconstruction of the PSD needs further refinement.
For the filtered back-projection it is assumed that the PSD rotates around a center point, which is located in the middle of the sinogram.Thus, we have to determine at which pixel the center point of rotation is located and then further cut the considered data range symmetrically around this center point.In order to calculate this center point, we average all bunch profiles in a section for which the reconstruction is considered (see Fig. 6).We then fit a normal distribution to the bunch profile considering pixels from 48 to 114, because later the profile data is overlapped the wake field of the electron bunch [45], which results in a distortion of the line shape.Figure 6 visualizes this calculation for the first section in Fig. 3a.The mean of the normal distribution is closest to pixel 102.For the other sections in Fig. 3b, we obtain similar fits and with pixel 102 and pixel 101 similar corresponding means.Finally, we chose the largest interval for which the center rotation point is in the middle and which is part of the interval from pixel 48 and 201.Thus, for the first two revolution plots, we used the pixels between pixel 48 and 156 and for the last revolution plot we used the pixels between 48 and 154.In general, these small differences have no influence on the data interpretation at the current signal-to-noise ratio, i.e. it would be also possible to use the same interval for all reconstructions.

C. Validation of the back-projection method with simulations
In general, our motivating assumption that the dynamics of an electron bunch can be approximated by a rigid rotation is a strong simplification.In fact, static deformation [46,47] as well as the dynamic formation of substructures [48] is caused by nonlinearities in the potential.
Thus, to validate the reconstruction empirically under the conditions present during the experiment, we simulate the PSD dynamics by numerically solving the Vlasov-Fokker-Planck equation using the Inovesa software program.Our group already demonstrated that the software can precisely predict the bunch current depending bursting behavior of electron bunches over long time scales [20].The PSDs shown in Fig. 2 were simulated using Inovesa v1.1.0(doi:10.5281/zenodo.3466767).This version implements the phase noise of the radio frequency cavities to excite the synchrotron oscillation, which is typically not part of simulations [49].The resulting potential is modified by the wake potential due to self-interaction with the CSR of the bunch and calculated according to the parallel plates shielding model [50].The simulations were carried out, closely to the experimental conditions, with a bunch current of 0.85 mA at the sawtooth bursting regime.The list of physical parameters is displayed in Table II.

FIG. 2
FIG. 2. a) Reconstruction principle of the PSD distribution with electro-optical sampling.Top: we assume a rigid rotation.Shown are the PSD configurations at different rotation angles and the final reconstruction.Bottom: the EO measurement (projection on the time axis) of the different PSD configurations can be interpreted as a Radon transform.Each PSD configuration is represented by a single column in the corresponding sinogram.The Radon transform is connected to the 2D Fourier transform of the original image via the Fourier slicing theorem.The original PSD can be reconstructed with a filtered back-projection.b) Validation of the reconstruction algorithm with simulated data obtained from the Vlasov-Fokker-Planck solver Inovesa[20].We used a standard ramp filter to avoid the blurring effect of the back-projection[35] in combination with a high frequency cut-off at a relative frequency of 0.1 to reduce the exaggeration of edges.Left side: comparison between the original and reconstructed PSD.The "original" PSD is evaluated at the 90 • frame, while for the reconstruction the complete sinogram has to be taken into the account.Right side: several comparisons during different times of a bursting cycle.Even in situations with prominent substructures, where bunch self-interaction is expected to be strong, the reconstruction shows a very good agreement with the original PSD, even in fine details.

FIG. 5
FIG. 5. a) Removal of the Gotthard chip dependent offset from the measured data.b) Definition of the region of interest.
. and B.K. were supported by BMBF contract No. 05K16VKA.E.B. and G.N. acknowledge support by the Helmholtz Presidents strategic fund IVF "Plasma accelerators".The authors thank N. Hiller, B. Hrer, T. Boltz, M. Brosi and P. Schreiber for fruitful discussions.We further thank J. Schmid and O. Speck for their technical support.

TABLE I .
Machine parameters during the experiment.

TABLE II .
Parameter settings used for the Inovesa simulation of the PSD dynamics.