First Simulations for the EuAPS Betatron Radiation Source: A Dedicated Radiation Calculation Code

: X-ray production through betatron radiation emission from electron bunches is a valuable resource for several research fields. The EuAPS (EuPRAXIA Advanced Photon Sources) project, within the framework of EuPRAXIA, aims to provide 1–10 keV photons (X-rays), developing a compact plasma-based system designed to exploit self-injection processes that occur in the highly nonlinear laser-plasma interaction (LWFA) to drive electron betatron oscillations. Since the emitted radiation spectrum, intensity, angular divergence, and possible coherence strongly depend on the properties of the self-injected beam, accurate preliminary simulations of the process are necessary to evaluate the optimal diagnostic device specifications and to provide an initial estimate of the source’s performance. A dedicated tool for these tasks has been developed; electron trajectories from particle-in-cell (PIC) simulations are currently undergoing numerical analysis through the calculation of retarded fields and spectra for various plasma and laser parameter combinations. The implemented forward approach evaluation of the fields could allow for the integration of the presented scheme into already existing PIC codes. The spectrum calculation is thus performed in detector time, giving a linear complex exponential phase; this feature allows for a semi-analitical Fourier transform evaluation. The code structure and some trajectories analysis results are presented.


Introduction
The EuAPS (EuPRAXIA Advanced Photon Sources) project [1], within the broader EuPRAXIA project [2], will provide an X-ray source produced by the betatron oscillations of a self-injected electron beam in a plasma bubble created in the Laser Wakefield Acceleration (LWFA) configuration.A gas jet produces a gas distribution consisting of a helium-neon mixture.A laser pulse with energy O(1) J and duration O(10) f s is focused on the gas, ionizing it in a highly nonlinear regime and in causing the self-injection phenomena of electrons [3].The transverse focusing force generated by the plasma blowout region determines the formation of electron betatronic orbits, leading to radiation emission.This process is commonly referred to as betatron radiation emission [4][5][6].Although it shares some similarities with an undulator, the properties of the emitted radiation are generally different, primarily due to the dependence of the oscillation strength on the distance from the axis and to the beam energy variation.The characteristics of the betatron radiation spectrum are well-described in the case of zero longitudinal acceleration and low oscillation strength [7,8], while a complete analytical model is lacking for the general case, which includes non-uniform acceleration and dephasing.Furthermore, in EuAPS, since the beam undergoing betatron oscillations will be self-injected through a highly nonlinear wave-breaking process, it will be subject to initial conditions with strong variability, depending on parameters such as plasma density and the spatio-temporal characteristics of the laser pulse.It is evident that detailed simulations of the process are necessary in order to assess the range of energy and angular divergences of the produced radiation and enable realistic considerations of the experimental setup required for measurements and radiation transfer to users.The dynamics of the electrons and the emitted radiation have been calculated separately due to the significant difference in energy scales; electron energies are of the order of O( 102 ) MeV, while single-photon energies are O(1-10) keV.This allows for neglecting the radiation recoil effect experienced by the electrons while emitting photons, justifying the separated treatment.The electron dynamics have been calculated using the numerical code FBPIC [9,10].For the radiation, a dedicated code has been developed to numerically calculate the delayed fields produced on a detector by the superposition of self-injected beam electrons along their trajectories [11].This paper will be mainly dedicated to the radiation calculation code.Its structure will be presented in the Materials and Methods section; the initial PIC data-sorting algorithm will be shown in detail, together with the delayed fields "forward" calculation, which is meant to facilitate the code's future integration into FBPIC or Smilei [12] PIC codes to achieve a direct radiation calculation.In the Results section, an analysis of a representative LWFA shot will be shown, where the versatility of the code allows us to highlight the presence of a collective transverse oscillation of the beam and to decompose the spectral analysis into different particle chunks which could display different amounts of coherence.The obtained results have been compared with the numerical spectrum calculated from the analytical trajectory of an ideal electron subjected to a linear focusing force and a quadratic energy variation [4,13], showing good agreement in terms of energy range with the most coherently oscillating chunk.

Materials and Methods
The structure of the code for computing the forwarded fields, and the corresponding spectrum produced by PIC-computed electron trajectories, will now be presented.The first subsection will provide some comments on the necessary data structuring, while the following two will briefly outline two of the possible approaches for calculating the radiation spectrum.We emphasize that both the data structuring and radiation calculation are developed for offline analysis but are structured in a way that allows for an easy extension to an online calculation directly integrated into a PIC code.The sorting algorithm presented in the next section is designed to operate concurrently on existing PIC data and iteratively on particle IDs.As it stands, the algorithm can easily be applied to an online adaptive sorting of particle trajectories to optimize a spectrum calculation by Fourier transform.The presented code was developed in python.

Data Structuring
To effectively work with numerical trajectories computed using PIC simulations, a significant effort is needed to structure the data.The IDs corresponding to the beam particles need to be selected, while the other low γ factor plasma particles have to be discarded.The PIC simulation window co-moves with the driving laser pulse; plasma electrons enter the moving window, while some of them get trapped in the plasma wake and later leave it.This dynamic is such that some of these beam particles may appear in the simulation window at one point in time and then disappear at a later moment, resulting in uneven trajectory lengths.To address this issue, we employ a particle selection script to process the raw output data from the PIC simulation.This script allows us to choose a specific set of particle IDs at a given point in time within the simulation, based on their position and momentum values.Once we have identified these selected IDs, we collect data related to these particles for all preceding and subsequent time steps of the simulation.We want to stress again that the number of beam particles may neither remain constant nor follow a consistent trend; the total number of processed particles will be, in general, larger than the maximum instantaneous "size" of the beam.Furthermore, even when the number of IDs remains the same from one timestep to the next, their order may not necessarily be preserved.The proposed approach involves inserting the data into at least two-dimensional arrays with dimensions [N timestep × N IDs ], where N timestep is the number of timesteps in the PIC simulation, and N IDs is the number of unique particle IDs:; these arrays will be called A f in as they are the final data format.For vector quantities such as position and momentum, the arrays will be of size [N timestep × N IDs × 3].In timesteps where a particular ID is missing, the value NaN (Not a Number) has been chosen for insertion, since it can be more easily identified compared to zero (indeed, zero could be the actual value of the variable for that specific ID at that particular timestep).The progressive filling of A f in will be performed using temporary arrays A temp of smaller dimensions [N timestep × N maxID ].
Here, N maxID represents the maximum instantaneous size of the self-injected beam in terms of particle number.This step will facilitate an orderly extraction of variables.
The data structuring takes place in the following steps: • First loop over the timesteps to identify the maximum number of concurrent IDs N maxID in a single timestep (maximum instantaneous "size" of the self-injected beam).

•
Identification of the array of unique IDs U IDs , with size N IDs ≥ N maxID .This condition is given by the fact that the total number of beam IDs can be greater then the maximum instantaneous number of IDs in the beam, due to the window entering/leaving dynamics of the self-injected electrons.

• Creation of temporary arrays
Second loop over the timesteps, sorting the variables by ID values, and inserting them into the A temp arrays.Note that, in a timestep (row of A temp ), the number of present IDs may generally be < N maxID .In these cases, the variables sorted by ID will be accumulated in the leftmost columns of A temp , while NaNs will be left in the free cells (see the leftmost block in Figure 1, green/yellow array matrix).In this way, the array A f in , which contains trajectories and other quantities ordered by time and separated by particle ID, is efficiently populated.This is achieved by advancing in an ordered manner along a "battlefront-like" line A c , which adapts to the non-uniformity of particle IDs' original distribution (see orange boxed cells in green/yellow array matrices in Figure 1, evolving from left to right) where A c finds the "right" ID cells in A temp (highlighted in green in thw figure), and the algorithm assigns the corresponding values to a single column in A f in (highlighted in blue in the figure).Proceeding in this orderly fashion through a simple loop over unique IDs reduces the computation time compared to using the classical search algorithms provided by the numpy library.

Spectrum Calculation
Using positions, moments, and fields obtained from the PIC simulations, structured as indicated in the previous section, it is possible to perform a parallel calculation of the electric field signals produced by each trajectory on a detector at any given spatial position.Collective effects like beam bunching due to emitted radiation can be ignored thanks to the wide energy spread and the high ondulator strength > 1 typical of the self injected beams; the superposition of many different radiation frequencies averages out the net effect.Moreover, the aforementioned separation of the beam and radiation descriptions avoids any possibility of simulating such a process with the presented scheme.The well-known expression for calculating the retarded fields takes the following form [11]: where e is the charge of the electron, c is the speed of light, γ is the Lorentz factor, R is the distance between the trajectory and the detector, ⃗ n is the unit vector pointing from the trajectory to the detector, and ⃗ β and ⃗ β are the velocity and acceleration of the particle normalized to c.Note that the given expression is valid in CGS units.Then, defining the quantity and performing its Fourier transform on the detector the differential intensity over frequency ω and solid angle Ω is obtained as follows: where the dimensions of the last expression is [erg s].In Equations ( 1) and ( 2), the subscript ret indicates that the field at point⃗ r at time t must be evaluated by finding an appropriate value of R that provides an intersection between the back-propagating light cone from point (⃗ r, t) and the trajectory of the particle.The values of the variables necessary for the field calculation should then be evaluated at time t − R/c, taking into account the signal propagation at velocity c.This conventional "backward" approach allows us to select a grid of times t d at which to evaluate the field on the detector and requires interpolating the trajectory at the points indicated by the reasoning just described, thus deriving a grid of times on the trajectory t t (Figure 2a).In the case of relativistic particles, interpolating the trajectories is a delicate and potentially risky operation, as even small variations in quantities like ⃗ β and ⃗ β can have a huge impact on energy or even result in unphysical situations (like superluminal velocities).It would be safer to interpolate momenta and then retrieve speed and acceleration, but the interpolation itself can have a big impact on calculation performance.The "forward" approach used in our code avoids interpolating the trajectories (preserving causality) and consists of evaluating the trajectory quantities only at the instants provided by the PIC simulation, thus choosing a grid of times on the trajectory t t and deriving a grid of times on the detector t d .Given the distance R between the trajectory and detector at each instant, the fields calculated with Equation (1) will be associated with the time t d = t t + R/c (Figure 2b).One of the advantages of the forward approach is its natural adaptability to the evolution of a numerical simulation.In the developments following this work, integrating this scheme directly into a PIC code will enable the reconstruction of the spectrum using complete information about trajectory evolution.Besides providing improved convergence of the Fourier integrals required for calculating spectral components, this operation will significantly reduce the amount of data to be saved during the simulation.The critical aspect of this approach concerns the temporal alignment of signals produced on the detector by different particles.Indeed, from Figure 2b, it is evident that adding a second trajectory that shares the same evenly spaced time grid (t t,2 = t t,1 ) would result in different arrival times of the signal on the detector compared to those produced by the first trajectory (t d,2 ̸ = t d,1 ) due to the particles' different positions.Furthermore, changing the position of the detector may also alter the order in which the elements of t d,2 and t d,1 alternate.Thus, two approaches to spectrum calculation are outlined, each with different computational complexities and different levels of extracted information (Figure 3): Total Field Calculation (TFC): To proceed in this manner, it is necessary to interpolate the fields produced by each particle on the detector.This is carried out using a specifically designed algorithm where • All the fields on the detector are combined into a single array, which is then sorted based on t d,all ; • Concurrently, a container of size t d,all is filled by interpolating the field of each particle i at all the missing time points t d,all−i ; • Finally, the ordered fields from the first step are directly summed into the container, providing the temporal profile of the total field.
Then, to compute the spectrum, it is possible to perform a Fourier Transform (FT) directly on the total field on the detector, using t d,all as the array of times.

Direct FT Calculation (DFC):
In this faster approach, we exploit the principle of the superposition of fields to decompose a single FT of the total field into a sum of FTs on the fields of individual particles.The obtained spectrum is identical to the TFC case, but the information about the temporal profile of the field is lost.For TFC, the effect of sorting and summing an increasing number of arrays composed of a growing number of elements results in an overall quadratic scaling.Conversely, for DFC, the expected and measured scaling is linear.

Results
The code presented in the previous section is currently in use and will be extensively employed in the next year to consolidate the numerical research phase, aiming to identify the most advantageous conditions for obtaining high photon fluxes in the energy range of 1-10 keV.Therefore, a validation of the correct functioning of the code and the determination of its applicability limits are mandatory.As a validation approach, a well-known and analytically tractable case has been chosen: the planar undulator emission [6,14].In the following first section, the numerically computed differential intensity in frequency and solid angle d 2 I/dΩdω from analytical undulator particle trajectories will be compared with the analytical expressions for the intensity itself.The following second section will present some analysis results of trajectories and betatron radiation emitted by self-injected particles in an LWFA process.

Validation
It is possible to demonstrate that the differential intensity produced along the axis by a particle during mid-plane orbits in a planar undulator takes on the following appearance (for further details, refer to Chapter 2 of [14]): where ξ is a variable defined from the undulator force, F m (ξ) is a function composed of Bessel functions of the first kind, ν m is a linear variable of frequency ω, N is the number of undulator periods, and m is the harmonic number.Comparisons between the intensities calculated numerically from particle trajectories using Equations ( 1)-( 4) and with the analytical expression given by Equation ( 5) are presented in Figure 4; with a temporal sampling rate sufficient to accurately reproduce the field profile, the results exhibit an excellent level of agreement.To estimate the minimum trajectory sampling rate, the radiation on the axis was considered, as it is the scenario where the "searchlight-like" emission determines the shortest signal time window on the detector.Assuming that the appreciable radiation is concentrated within a cone with an angle of the order of 1/γ, and that the temporal duration of the signal on the detector corresponds to one complete passage of the emission cone over the detector's position, the corresponding angle covered by the tangent to the trajectory will also be of the order of 1/γ.Let us assume that the transverse coordinate x as a function of the longitudinal coordinate z is approximately with K undulator strength and k U undulator wavevector.Then, for small angles, the angle of the tangent to the trajectory will be given by Calculating from Equation (7), the ∆θ as a function of the ∆z of the trajectory and equating it to 1/γ, it is possible to determine the ∆t of the trajectory generating a visible signal on the detector by substituting the velocity v = βc: It has been verified that, by choosing a timestep small enough compared to the one indicated by the condition given by Equation ( 8) (e.g., ∆t sim = 0.1∆t), the oscillation of the field on the detector is correctly reconstructed, and there is an agreement between the theoretical and numerically calculated spectra.For larger ∆t sim , the numerically calculated spectrum is generally higher than the analytical one, while still maintaining harmonic alignment.A more detailed investigation of how the numerical spectrum deviates from the analytical one as a function of undersampling level will be pursued in the future.If a scaling law is identified, it would be possible to perform preliminary spectrum calculations using undersampled trajectories and rapidly assess the order of magnitude of the number of produced photons for the choosen plasma-laser parameters configuration.

Betatron Radiation Analysis: A Peculiar Case
As a demonstration of the performance of the presented code, we showcase some results from the analysis of betatron radiation emitted by a self-injected electron beam during an LWFA shot simulated with FBPIC.Since the employed setup is an ideal one (bigaussian laser, flat top plasma density), these results must be considered as a benchmark for future more realistic plasma simulations.The case involves a plasma density n p = 2.44 × 10 18 cm −3 and a laser pulse energy E l = 2.5 J.Other relevant laser pulse parameters, common to all simulations, include the wavelength λ l = 800 nm, temporal length τ l = 25 f s, and beam waist w l = 15 µm, assuming a perfectly Gaussian transverse spot.The laser pulse is polarized along the x axis.In this instance, 10,000 macroparticles were self-injected, and 600 timesteps of the PIC simulation were selected, requiring a computation time of 30 min for a total of 100 different detector positions.Figure 5 visually represents the analysis of the presented case.In Figure 5A, the trajectories x − z of the selected self-injected beam are shown in gray, the median of the x coordinate (y coordinate) as a function of z is depicted in full orange (blue), and the corresponding mean ∓σ (standard deviation) is shown with dashed lines.A clear collective oscillation of the beam along the x axis is observed, which coincides with the laser's polarization.Whether this behavior is an isolated occurrence or a systematic feature, and whether it is a physical characteristic or a PIC numerical artifact, are aspects currently under evaluation.Nonetheless, if this characteristic proves to be genuine and controllable, it could provide significant advantages in terms of the resulting field intensity and photon flux.
In Figure 5B, the intensity of total radiation (integrated over all selected frequencies corresponding to the range 10 −1 -10 2 keV) is represented as a color-coded variable on the angular spot centered along the laser propagation axis; the differences in intensity and radiation divergence along the two axes confirm that which was mentioned about the trajectories in Figure 5A.In Figure 5C, a more detailed analysis of the beam was performed by dividing the particles into three chunks (represented by three shades of blue) based on their similarity with the oscillating median in Figure 5A.The subfigure on the top shows the average γ of the particles in each chunk as a function of z.The core of the beam (dark blue) exhibits an initial clear collective oscillation along the x axis (on the left), which gradually fades due to a dephasing process and reappears partially in the final part (on the right).In orange, an analytical trajectory is overlaid, calculated assuming a parabolic energy profile obtained from a fit of the average γ reported in the subfigure.This shows good agreement with the initial and final sections of the trajectories, demonstrating that the core moves approximately like a single ideal electron.It can be demonstrated that, in the case of parabolic longitudinal acceleration, the analytical trajectory for an electron subjected to a linear focusing transverse force is approximated by a linear combination of Legendre polynomials.A dedicated work focused on this type of analysis, aimed at obtaining an approximate analytical expression of the spectrum, is currently under development.Finally, in Figure 5D, the spectra emitted on-axis by each of the three isolated chunks in Figure 5C are shown, along with the numerical spectrum calculated from the analytical trajectory (orange).In the top-right inset, the total spectrum integrated over the solid angle is presented.Moving from light blue to dark blue, three effects can be observed:

•
A decrease in the cutoff energy.• A decrease in the extension of the low-energy background (left).

•
An increase in the peak intensity.The first effect can be attributed to the progressive reduction in the trajectories' average amplitude; when experiencing less acceleration (recall that the transverse force is linear with the distance from the axis), the related emitted radiation will be less energetic.The second and third effects seem to be related to the increased coherence and energy as one moves towards the core; a more coherent oscillation will enhance total intensity, while the concurrent reduction of noise will give a neater spectrum profile.Finally, although the (scaled for comparison) spectrum associated with the analytical trajectory (orange) shows different characteristics due to the evident simplification of the system, the energy extension remains consistent, particularly with that of the spectrum associated with the core trajectories (dark blue).

Conclusions
The code presented, specifically developed for calculating betatron radiation spectra within the EuAPS project, is suited for the calculation of radiation emission from any set of particle trajectories whose energy greatly exceeds that of the radiation itself (i.e., radiation recoil is not included).So far, trajectories obtained from PIC calculations are employed to compute delayed fields radiation on a detector.These fields are then Fouriertransformed to retrieve the spectrum.The forward delayed fields calculation approach, presented in Section 2.2, enables direct integration into any PIC (Particle-in-Cell) code, easing data management and extracting all possible trajectories' information for spectrum calculation.The code was validated through analytical undulator spectra comparison and it is currently used on a daily basis in the parametric scanning of laser-plasma parameter processes.A systematic evaluation of the calculated intensity scaling for different levels of undersampling may also enable a quick and approximate usage based on a few trajectory points, providing an order of magnitude estimation of the number of produced photons.

Figure 1 .
Figure 1.Three steps of the alignment loop for variables in sorted arrays.A c advances only when and where there is a match between the current ID and the local IDs.The variables in A f in are updated (light blue) with the corresponding values extracted from A temp (light green).Note that the number of columns in A temp and A f in is not generally the same.
A f in filled with NaN values, with dimensions [N timestep × N IDs ] or [N timestep × N IDs × 3] (see leftmost block in Figure 1, blue array matrix).• Initialization of a counter array A c with size N timestep set to all zeros; this object will keep track of the timesteps where the jth ID variables are found during the alignment process.• ID alignment: a single loop with index (j) over U IDs .A c is incremented only at positions i where {ID temp } i,A c = {U IDs } j .The same indices for which this condition is satisfied are used to extract the variables from A temp and insert them into A f in .A visual representation of the loop is presented in Figure 1, going from the left box to the right one.

Figure 2 .
Space-time diagrams for delayed field usage patterns.The detector's position is shown in dashed orange, the particle's trajectory in solid green, and the light cones of the signals in continuous black at 45 • angles.(a) Backward approach, with evenly spaced times on the detector, requiring interpolation along the trajectory.Given the highly relativistic regime, in some situations this could lead to causality violation.(b) Forward approach, with evenly spaced times along the trajectory, resulting in uneven times on the detector.No interpolation is needed.Causality is always preserved.

Figure 3 .
Schemes for total spectrum calculation through Fourier transform.(a) Visual schematization of the TFC and DFC approaches.Broken lines of different colors represent the fields calculated by different particles.The shaded yellow line represents the total field at the detector.(b) Performance analysis of the two approaches as a function of the number of particles analyzed, based on the convergence of the mean squared error with respect to the polynomial degree of the runtime fit.

Figure 4 .
Comparison between numerical spectra (full lines) and analytical spectra (dashed lines) along with the evaluation of the theoretical position of the first harmonic (vertical dashed lines) for an electron in an undulator.Evaluations were performed for several values of the undulator strength K, and the on-axis spectra for K = 0.3 and K = 3 are shown in panels (a) and (b), respectively.

Figure 5 .
Figure 5. Analysis of betatron radiation spectrum emitted by a PIC-simulated LWFA self-injected beam.(A) Trajectories x − z of the beam, with medians and means ∓σ along the two transverse axes ( x and ŷ).(B) Radiation intensity on angular spot.(C) Division of beam particles into three chunks based on their similarity level with the median along the x axis; superimposed: analytical trajectory of a single electron with quadratically approximated acceleration (orange).(D) Spectra of the three chunks from figure (C) along with the (scaled) numerical spectrum of the analytical trajectory (orange) and the total spectrum integrated over the solid angle (top right).