5D tomographic phase-space reconstruction of particle bunches

We propose a new beam diagnostics method to reconstruct the phase space of charged particle bunches in 5 dimensions, which consist of the horizontal and vertical positions and divergences as well as the time axis. This is achieved by combining a quadrupole-based transverse phase-space tomography with the adjustable streaking angle of a polarizable X-band transverse deflection structure (PolariX TDS). We demonstrate with detailed simulations that the method is able to reconstruct various complex phase-space distributions and that the quality of the reconstruction depends on the number of input projections. This method allows for the identification and visualization of previously unnoticed detailed features in the phase-space distribution, and can thereby be used as a tool towards improving the performance of particle accelerators, or performing more accurate simulation studies.


I. INTRODUCTION
Detailed knowledge of the particle beam properties is required for simulating, optimizing and operating highquality particle accelerators.Ideally, the full phase-space distribution should be known, however, only lower dimensional projections and statistical beam parameters are typically measured.While the full phase-space distribution is available from start-to-end simulation studies, such particle distributions do not always accurately represent the actual beam.This is due to, for example, an incomplete knowledge of the beamline and its instabilities like drifts and jitters or an incomplete knowledge of the initial properties of the beam.
At the same time, imperfections in the beamline or collective effects can impact the phase-space distribution and lead to correlations between the planes, or complex, non-Gaussian phase-space distributions that are not well described by lower dimensional projections or the statistical beam parameters.For example, space-charge forces [1] or coherent-synchrotron radiation [2] can lead to correlations between the longitudinal and transverse planes.Rotated beamline elements or stray magnetic fields on the photocathode and in the gun region can lead to 2dimensional (2D) emittance growth [3,4].For ultra-short bunches in the fs-range the energy gain deviations induced by transverse offsets in a radio frequency (RF) gun and the travel time difference induced by transverse offsets in transverse focusing elements also affect the phase space distribution [5].
To characterize such effects in detail, diagnostic methods to reconstruct high-dimensional phase-space distributions are required.For a low-energy 2.5 MeV H − beam, a measurement of the 6D charge density has been performed [6] in a dedicated beamline with slit masks.As an alternative, tomographic methods can be used which are readily applicable in most existing beamlines to access phase-space distributions in 2D with micrometer transverse resolution [7][8][9][10].Methods exist to measure higher-dimensional distributions, such as the timeresolved transverse phase-space distribution [11], the full 4D transverse phase space [12][13][14][15], as well as the 3D charge density distribution [16][17][18][19].Furthermore, efforts to obtain information on the longitudinal phase space [20] and projections of the 6D phase space [21] are made using machine learning techniques.
In recent years, a collaboration between CERN, DESY, and PSI has successfully developed and tested a novel polarizable X-band transverse deflection structure (PolariX TDS) [19,22,23].In addition to the conventional TDS features, this structure is capable of streaking bunches along any transverse direction, opening up new opportunities for more detailed diagnostic methods, such as a recently experimentally demonstrated 3D charge-density reconstruction [19].Several PolariX TDSs have been installed at DESY and PSI, including two structures at the ARES linear accelerator [24][25][26][27], a facility dedicated to accelerator research and development at DESY.
Here, we present a technique enabled by the PolariX TDS that allows for the tomographic reconstruction of the 5D (x, x , y, y , t) phase space, where x, x and y, y are the transverse position and divergence and t is time.The working principle consists of combining the streaking along various transverse directions with a transverse phase-space tomography based on a quadrupole scan [28,29].The unique features of the 5D phase space tomography method, i.e., resolving complex phase-space distributions and their correlations, are studied through simulations based on the ARES beamline layout.The reconstruction accuracy is characterized for two different test distributions and the influence of the number of rotation and streaking angles on the reconstruction accuracy is investigated.The 5D phase-space information obtained with this method is useful to optimize and improve the beam quality of, e.g., free-electron lasers or beam driven plasma accelerators by detecting previously FIG. 1. Sketch of the ARES beamline layout used for the 5D phase-space tomography simulation studies and working principle of the method.The gap between the third and fourth quadrupole is occupied by a bunch compressor which is not used for the presented study.Two PolariX TDSs are present in the ARES beamline.As mentioned in the text, only the first TDS is used in the presented simulation study.
unnoticed correlations and other features in the distribution.Furthermore, the reconstructed phase-space distribution is useful to benchmark existing simulation codes, or as an input to perform detailed simulation studies of advanced acceleration schemes.

II. WORKING PRINCIPLE
The presented 5D tomography method consists of performing a 4D transverse phase-space tomography for each longitudinal slice of the bunch.Here, the underlying principle of the tomography is to use lower-dimensional projections of an object along different angles to reconstruct its distribution in a higher-dimensional plane.These higher-dimensional distributions are obtained by using a tomographic reconstruction algorithm such as the Filtered Back-projection [30], the ART (Algebraic Reconstruction Technique) [31], the SART (Simultaneous Algebraic Reconstruction Technique) [32], or the MENT (Maximum Entropy Tomography) [33] algorithm.
The transverse tomography of a particle bunch in an accelerator is performed by rotating the transverse phase spaces (x, x ) and (y, y ).This rotation is done by changing the phase advance µ x,y experienced by the particle distribution.It is especially useful to perform this in the normalized phase space, where for linear optics the phase advance µ x,y is equivalent to the rotation angle θ x,y of the phase spaces [34].The conversion from real phase-space coordinates x, y and x , y to normalized phase-space co-ordinates x N , y N and x N , y N is given by: where β x,y and α x,y are the Courant-Snyder parameters [35].To also resolve the correlations between the two transverse planes, their phase advances are controlled simultaneously.
The longitudinal information is obtained by streaking the bunch with the PolariX TDS along various transverse directions.This is required due to the overlap of the longitudinal information with the transverse one in the plane parallel to the streaking direction.Therefore, the full transverse distribution is only available when combining the projections of the bunch under various streaking directions.
The working principle of the 5D phase-space tomography is depicted in Fig. 1.For a fixed combination of transverse rotation angles (θ x , θ y ), the beam is streaked with the TDS, and its projection is recorded on a downstream screen.This is repeated for different streaking angles.These screen images are used to reconstruct the 3D charge-density distribution (x, y, t) of the bunch [16][17][18][19].The transverse profiles are reconstructed at the location of the screen.The longitudinal information is reconstructed at the TDS center.This assumes that the longitudinal distribution does not change between the two locations which is valid as long as the longitudinal displacement of a particle does not exceed the longitudinal resolution of the reconstruction.The reconstruc-tion procedure is repeated for all desired phase advance combinations.Each 3D reconstruction can be regarded as the projection of the (x, x , y, y , t) phase space onto the (x, y, t) plane for different transverse phase advance combinations.To reconstruct the higher-dimensional distribution, the 3D reconstructions are combined.This is done for each longitudinal slice t s individually.Using the reconstructed sliced projections (x, y) ts , the transverse charge-density distribution (x, x , y, y ) can be reconstructed for each slice t s , similar to the reconstruction of the 4D transverse charge density in [12,13].This distribution is reconstructed at a location upstream of the TDS and of all the quadrupoles that are used to scan the phase advance.By combining all longitudinal slices, the 5D charge-density distribution is obtained.
To obtain the optimal beamline settings for the 5D tomography, several factors have to be considered.First, the setup should allow for accurate tuning of the transverse phase advance in both planes in a range of 180°t o obtain a good sampling of the transverse distribution.Second, the longitudinal resolution of the reconstruction should be sufficient to identify all features of interest.This longitudinal resolution R is given by [36] where σ off is the transverse spot size at the screen when the TDS is switched off, p is the average momentum of the bunch, c the speed of light in vacuum, e the elementary charge, f the TDS X-band working frequency, V the peak voltage, and L the drift length between the TDS center and the downstream screen.If transverselongitudinal correlations are present in the distribution, the projected transverse spot size can be larger than the spot size of each individual slice.The longitudinal resolution could therefore be improved by determining the resolution using the transverse spot size of each slice.These spot sizes of the individual slices can be obtained from the screen image of the streaked bunch.Since the spot size along the streaking direction is the one of interest, the screen image streaked perpendicular to the streaking angle of interest has to be considered.With this approach, also a variable longitudinal resolution along the bunch would be possible.The transverse spot size is directly linked to the beta functions at the screen.These need to be matched to a constant value per plane in order to obtain a uniform longitudinal resolution for all transverse phase advance combinations.For a fixed TDS setting, the final longitudinal resolution is given by the largest transverse spot size of all streaking directions for all phase advance combinations.Furthermore, the screen resolution needs to be taken into account to ensure that the pixel size can resolve the features of the projected bunch.
The longitudinal information of the distribution is obtained by converting the screen coordinate in the direction of the streaking angle u * to the longitudinal position along the bunch ∆t according to where S is the shear parameter defined as [11] When measuring the longitudinal bunch information with a TDS, typically, screen images are recorded at both RF zero crossings to account for transverse-longitudinal correlations within the bunch.Here, screen images are recorded only at the first zero crossing to keep a feasible measurement and analysis time.

III. PROOF-OF-PRINCIPLE SIMULATION STUDIES
The present study investigates the capabilities of the 5D tomography method to reconstruct phase-space densities with correlations and complex features.This is done by studying two different distributions in simulation studies and comparing the known input distribution to the reconstruction result.For these studies, the ARES beamline layout shown in Fig. 1 is used, which consists of five quadrupoles and two PolariX TDSs.The first distribution is based on realistic beam parameters at ARES.This distribution is a Gaussian distribution with imprinted correlations in the (x -z), (y -z) and (x -y ) planes.It is depicted in Fig. 2 with the label "original".The second distribution is used to test the ability of the method to reconstruct arbitrary phase-space distributions.It consists of three superimposed Gaussian beams with transverse offsets with respect to each other and longitudinal correlations in the transverse momenta, resulting in a complex phase-space structure.This distribution is depicted in Fig. 3 with the label "original".Both distributions are generated in ocelot [37,38] and feature the same initial Courant-Snyder parameters, which allows for the use of the same quadrupole settings for the two distributions.The second distribution features a larger transverse emittance due to its multi-bunch structure.All the initial beam parameters are listed in Table I, where E is the average energy, σ E the RMS energy spread, Q the charge, σ t the RMS bunch duration, and x,y ( n x,y ) the geometric (normalized) RMS emittance.The five quadrupoles and one of the PolariX TDSs are used to achieve the desired phase advances and streaking angles.To determine the longitudinal resolution of the reconstructed distribution, the maximum unstreaked spot size at the screen is considered.It is obtained by tracking the beam distribution through the beamline for all phase-advance combinations with the TDS switched off.The maximum RMS transverse spot size is then determined by comparing all transverse projections onto all streaking planes.From Eq. 2, a peak voltage V of  3.6 MV is required to achieve a longitudinal resolution of 20 fs, with an X-band TDS frequency of 11.992 GHz, a maximum spot size σ off of 247 µm, and a drift L of 7.21 m.In practice, each of the available PolariX TDSs would be able to streak the beam with up to 20 MV peak voltage.By using both TDSs with 20 MV each, assuming the same unstreaked spot size and the drift length of L = 6.67 m, a resolution of 1.91 fs could be achieved.
In this study, only the first TDS in Fig. 1 is used to achieve the desired 20 fs, and a constant shear parameter for all phase-advance settings and streaking angles is assumed.The streaking angle of this TDS is varied over a range of 180°in 40 steps.The quadrupole strengths are set to scan the transverse phase advances over a range of 180°in 40 steps.The projections of the distributions are recorded on the screen station downstream of the TDS.The simulated screen has a size of 2.02 cm × 2.02 cm to fit ±4 σ of the streaked RMS bunch duration, and 2001 × 2001 pixels, which corresponds to a resolution in the range of the ARES screen stations.Moreover, it ensures a minimum of 4 pixel per transverse RMS spot size and more than 8 pixel per transverse spot size for 97 % of the scanned phase-advance combinations.For the reconstruction a region of interest of 201 pixels is selected in the transverse direction which fits ±4 σ of the transverse RMS beam size.The reconstruction of the transverse distribution is performed at the location of the screen station upstream of the first used quadrupole.
All simulations are performed in ocelot [37,38] using distributions with 4 000 000 particles.Second-order transfer maps are used for all elements except the TDS, where only a first-order transfer map is available.Spacecharge effects are not taken into account because according to the laminarity parameter defined in [1,39], the bunch is emittance dominated throughout the beamline.This is confirmed by simulations including space-charge forces for the phase-advance setting with the maximal laminarity parameter along the beamline.The tracked distribution is compared to the corresponding simulation without space charge for vertical and horizontal streaking.No significant differences between the two cases were observed and therefore space-charge effects are neglected in the presented studies.
The final transverse reconstruction is performed in normalized phase space.For this, the reconstructed sliced (x, y) projections at the screen downstream of the Po-lariX TDS are normalized using the beta functions obtained by propagating the initial Courant-Snyder parameters through the beamline for each phase-advance combination.The final reconstruction, with the transverse distribution obtained at the screen station upstream of the first used quadrupole and the longitudinal information obtained at the TDS center, has a size of 201 bins in all transverse dimensions and 80 bins in the longitudinal plane, which, respectively, results in a resolution of 5 µm/ √ m and 20 fs.All tomographic reconstructions in this study use a Python scikit-image [40] implementation of the SART algorithm [32].This algorithm uses an iterative solver to obtain the reconstruction.A good reconstruction is already obtained in a single iteration.Increasing the number of iterations improves the reconstruction of sharp features in the distribution at the cost of an increased noise level [32].For each tomographic reconstruction in this study, two iterations over the algorithm are performed, where the second iteration uses the first reconstruction as an input for the initial reconstruction estimate.A condition is imposed in the algorithm to ensure charge invariance and strictly positive chargedensity values.

IV. RECONSTRUCTION RESULTS
For the first distribution of a Gaussian beam with imprinted correlations, the reconstructed 5D phase space is shown in Fig. 4 as a projection onto the (x N , y N , t) phase space.The reconstructed Courant-Snyder parameters, emittances, and bunch lengths are listed in Table I.Tomographic methods such as the SART tend to introduce spurious charge density even far off-axis where no charge should be present [19] especially when only a limited number (< 100) of rotation angles are used.Therefore, to analyze the emittances and Courant-Snyder parameters of the reconstructed distributions only values within a ±3 σ range from the bunch center are considered.The bunch duration σ τ is calculated as the standard deviation of the distribution projected onto the time axis.Excellent agreement with relative discrepancies of 5 % is achieved for the beam parameters.The imprinted correlations in the distribution are analyzed in normalized phase space and quantified by the slope of a linear fit which is applied to the corresponding 2D projections of the phase space.The original and reconstructed slopes are listed in Table II and show relative discrepancies be-FIG.4. 3D render of the reconstructed 5D phase space of the Gaussian distribution with correlations.The distribution is projected onto the normalized (x N , y N , t) phase space.The charge density is normalized to the maximum value of the 3D distribution.The transverse (x N , y N ) central slice and the slices ± 200 fs from the center are shown exemplarily.In addition, the 2D projections of the 3D distribution are shown in a blue color coding.The 3D render is done with pyvista [41], a Python implementation of VTK [42].

Plane
Original Reconstruction (x N -y N ) −0.848 −0.814 Additionally, all ten 2D projections of the 5D distribution are shown in Fig. 2 and compared to the original input distribution.The projections are normalized to the maximum value of the projection of the original distribution and the normalized difference between original and reconstructed projection are shown.All ten projections are reconstructed with excellent agreement with a maximum relative discrepancy of 13 %.The largest discrepancy appears in the (x N , y N ) projection where generally a widening of the reconstructed distribution is observed.
The second distribution has a complex phase-space distribution and therefore a larger maximal unstreaked transverse RMS spot size at the screen of 722.07 µm.As a result, some adjustments are required compared to the settings described in Section III.Using the same TDS and quadrupole settings as for the first distribution, a longitudinal resolution of 57.21 fs is achieved, and 27 lon- gitudinal slices are reconstructed.Due to the bigger spot size, a larger region of interest of 431 bins is selected for the reconstruction, which fits ±3 σ of the maximum RMS spot size.All other parameters remain unchanged with respect to the reconstruction of the first distribution.
The reconstructed bunch duration and Courant-Snyder parameters for the multi-structure beam agree well with the original parameters with discrepancies of around 4 %.However, the reconstructed emittances show large deviations of up to 24 %.Possible reasons could be an insufficient number of projection angles and screen resolution to resolve the sharp and thin features of the distribution.Due to the substantial increase in required computational power and memory when simulating many projection angles and a high screen resolution, this is not investigated further.The reconstructed 5D phase space is shown in Fig. 5 as a projection onto the (x N , y N , t) phase space which exhibits a more interesting structure than the previously shown (x N , y N , t) phase space.The substructure consisting of three sub-beams is well reconstructed and clearly visible.All 2D projections of the 5D phase-space distribution are shown in Fig 3 .Also here the complex substructure of the bunch is recovered and the reconstruction of the detailed features can be appreciated.The projections onto the (x N , y N ), (x N , z), and (y N , z) phase spaces are accurately reconstructed and show relative discrepancies below 20 %.The largest relative discrepancies of up to 60 % appear in the (x N , x N ), (y N , y N ), and (x N , y N ) projections where a blurring of the reconstructed phase space is observed.However, this blurring does not prevent a precise qualitative description of the phase space.Therefore, this result demonstrates the method's potential to reconstruct the 5D phase-space distribution of complex distributions.

V. INFLUENCE OF ROTATION ANGLES
The influence of the number of transverse rotation angles, as well as streaking angles, on the reconstruction accuracy is studied in this section.This is done to find an optimal balance between the required preparatory simulation and later experimental measurement time, which is determined by the number of scanned angles, and the final accuracy of the reconstruction.For this, several combinations of rotation angles and streaking angles are compared.The same Gaussian distribution with correlations as studied in Section III and IV is used.The tracking is performed for 60 transverse rotation angles for both transverse planes and 100 streaking angles.The tomographic reconstruction is performed using different sub-datasets to study a variety of rotation and streaking angle combinations without having to repeat the computational expensive tracking which runs approximately 42 h on 6 AMD EPYC 75F3 nodes on a computing cluster.As a measure for the accuracy of the reconstruction, the relative discrepancies η i,j = (σ recon i,j − σ input i,j )/σ input i,j are computed for each longitudinal slice j and transverse 1D projection of the reconstruction, where i = x, x , y, y .The values σ recon i,j are obtained from the standard deviation of a Gaussian fit to the reconstructed distribution, the values σ input i,j are obtained from the standard deviation of the input distribution to be less sensitive to a low number of particles in some slices.The value η i,j of each slice is weighted by the corresponding charge Q j of the slice, and the η i,j values are summed for each transverse plane i individually: Figure 6 shows the discrepancies for the four transverse planes and various combinations of rotation and streaking angles.The reconstruction accuracy improves with increasing number of rotation and streaking angles.Increasing the number of streaking angles from 50 to 100 results in no significant improvement of the reconstruction accuracy which in this case is dominated by the number of transverse rotation angles.To obtain discrepancies of 5 % in all reconstructed planes, ∼ 50 transverse rotation angles and ∼ 50 streaking angles are required.To obtain discrepancies of 10 % in all planes, ∼ 30 transverse rotation angles and ∼ 25 streaking angles are required.In both cases this is limited by the accuracy of the reconstruction in the x plane.This is probably due to the (x − t) correlation in the distribution and the resulting smaller RMS slice divergences.Furthermore, a difference between the x and y planes is observed.Again, this is likely caused by the correlations in these planes which are larger in (x −t) than in (y −t).For a pure Gaussian distribution without correlations, no significant differences between the accuracy of the x, x planes and y, y planes is observed.
The relative discrepancy of the bunch duration is obtained by comparing the standard deviation of the input distribution to the standard deviation of a Gaussian fit of the 1D projection of the reconstructed distribution onto t.It is unaffected by the number of transverse rotation and streaking angles and the average relative discrepancy is 0.273 ± 0.004 %.The independence on the number of rotation angles can be explained by the imposed constraint of the method to preserve the charge of each longitudinal slice once the 3D distribution (x, y, t) is reconstructed.No impact of the number of streaking angles on the accuracy of the bunch duration reconstruction is observed.
The presented study uses the SART algorithm for the tomographic reconstruction.The 5D tomography method itself, however, is independent of the used recon-struction algorithm.Therefore, one possible approach to further reduce the required transverse rotation and streaking angles is to explore other reconstruction algorithms, such as the MENT algorithm [33,43].Recent studies [44,45] using a machine learning approach to reconstruct the phase space show promising results and their application to the 5D tomography method could be explored.

VI. CONCLUSION
The presented tomographic method enables the accurate reconstruction and visualization of the full 5D (x, x , y, y , t) phase-space distribution of accelerated bunches.Simulation studies for a Gaussian distribution with a longitudinally correlated transverse momentum offset show excellent agreement with the original beam parameters and correlations with discrepancies 5 %.For a complex, multi-beam phase-space distribution, all features are clearly reconstructed but are more smeared out.This is expected to stem from the thin and sharp features of the distribution which require a larger number of rotation angles and a small pixel size of the recording screen to be reproduced accurately.In general, an improved reconstruction accuracy is achieved by increasing the number of transverse rotation and streaking angles.To improve the accuracy of the method for small numbers of rotation and streaking angles, different reconstruction algorithms could be explored.Experimental measurements are planned in the near future to demonstrate the practical applicability of the method.

FIG. 2 .FIG. 3 .
FIG.2.All ten 2D projections of the 5D phase-space distribution of the Gaussian distribution with correlations.The top rows show the projection of the original distribution, the middle rows show the projections of the reconstructed phase-space distribution, and the bottom rows show the difference between the two.All projections are displayed in normalized transverse phase space and the density is normalized to the maximum value of the corresponding 2D original projection.

FIG. 5 .
FIG. 5. 3D render of the reconstructed 5D phase space of the multi-beam distribution.The distribution is projected onto the (xN , y N , t) phase space.The charge density is normalized to the maximum value of the 3D distribution.The transverse (xN , y N ) slices at −246 fs, 62 fs and 369 fs are shown exemplarily.In addition, the 2D projections of the 3D distribution are shown in a blue color coding.The 3D render is done with pyvista [41], a Python implementation of VTK [42].

FIG. 6 .
FIG.6.Relative discrepancies, as defined in Eq. 5, between the reconstructed and input distribution of all transverse planes for various numbers of transverse rotation and streaking angles.The dotted gray and black lines show the 5 % and 10 % levels, respectively.

TABLE I .
Beam Parameters of the Original and Reconstructed Distributions