First Simulations for the EuAPS Betatron Radiation Source

1. Introduction

The EuAPS (EuPRAXIA Advanced Photon Sources) project [1], within the broader EuPRAXIA project [2], will provide a source of soft X-rays produced by betatron oscillations of a self-injected electron beam in a plasma bubble created in the Laser Wakefield Acceleration (LWFA) configuration. This process is commonly referred to as betatron radiation emission. 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 [3][4][5], while a complete analytical model is lacking for the general case, which includes non-uniform acceleration and dephasing. Furthermore, in EuAPS, the beam undergoing betatron oscillations will be self-injected through a highly nonlinear wave breaking process, resulting in a strong variability of the beam’s initial conditions 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 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, which allows neglecting the radiation recoil experienced by the electrons. The electron dynamics have been calculated using the numerical code FBPIC [6] [7]. For the radiation, a dedicated code has been developed to numerically calculate the delayed fields produced by the superposition of electron trajectories on a detector [8]. The code, whose structure will be detailed in the next section, has been developed in a forward approach to facilitate its future integration into FBPIC or Smilei [9] PIC codes, for direct radiation calculation. In the last section the results of the analysis of a representative shot will be presented, where the versatility of the code has allowed highlighting the presence of a collective transverse oscillation of the beam and decomposing the spectral analysis into different particle chunks which could display some 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 (due to the ion channel) and a quadratic energy variation (driven by the dephasing effect), showing good agreement in terms of energy range with the most coherently oscillating chunk.

2. Materials and Methods

The structure of the python 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. The presented code was developed in python.

2.1. 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 $\gamma$ factor plasma particles have to be discarded. Some of these beam particles may appear in the simulation window at one point in time and then disappear at a later moment. 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}\times 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_{fin}$ as they are the final data format. For vector quantities such as position and momentum, the arrays will be of size $[N_{timestep}\times N_{IDs}\times 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_{fin}$ will be performed using temporary arrays $A_{temp}$ of smaller dimensions $[N_{timestep}\times 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}\ge N_{maxID}$.
• Creation of temporary arrays $A_{temp}$ with dimensions $[N_{timestep}\times N_{maxID}]$ or $[N_{timestep}\times N_{maxID}\times 3]$.
• 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 leftmost block in Figure 1, green/yellow array matrix)
• Creation of final arrays $A_{fin}$ filled with NaN values, with dimensions $[N_{timestep}\times N_{IDs}]$ or $[N_{timestep}\times N_{IDs}\times 3]$ (see leftmost block in Figure 1, blue array matrix).
• ID alignment: a single loop with index (j) over $U_{IDs}$. A counter array $A_c$ with size $N_{timestep}$ is incremented only at positions i where ${\left\{I{D}_{temp}\right\}}^{i,A_c}={\left\{U_{IDs}\right\}}^j$. The same indices for which this condition is satisfied are used to extract the variables from $A_{temp}$ and insert them into $A_{fin}$. A visual representation of the loop is presented in Figure 1, going from left box to the right one.

In this way, the array $A_{fin}$, 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 green in figure), the algorithm assigns the corresponding values to a single column in $A_{fin}$ (highlighted blue in figure). Proceeding in this orderly fashion through a simple loop over unique IDs reduces computation time compared to using classical search algorithms provided by the numpy library.

2.2. Spectrum calculation

Using positions, moments, and fields obtained from the PIC simulations, structured as indicated in the previous section, it is possible to perform parallel calculation of the electric field signals produced by each trajectory on a detector at any given spatial position. The well-known expression for calculating the retarded fields takes the following form [8]:

$$\overrightarrow{E}(\overrightarrow{r},t)=e{\left[\frac{\overrightarrow{n}-\overrightarrow{\beta }}{{\gamma }^2{(1-\overrightarrow{\beta }·\overrightarrow{n})}^3{R}^2}\right]}_{ret}+\frac{e}{c}{\left[\frac{\overrightarrow{n}\times ((\overrightarrow{n}-\overrightarrow{\beta })\times \dot{\overrightarrow{\beta }})}{{(1-\overrightarrow{\beta }·\overrightarrow{n})}^{3}R}\right]}_{ret}$$

(1)

where e is the charge of the electron, c is the speed of light, $\gamma$ is the Lorentz factor, R is the distance between the trajectory and the detector, $\overrightarrow{n}$ is the unit vector pointing from the trajectory to the detector, and $\overrightarrow{\beta }$ and $\dot{\overrightarrow{\beta }}$ 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

$$\overrightarrow{A}(\overrightarrow{r},t)={\left(\frac{c}{4\pi }\right)}^{\frac{1}{2}}{\left[R\overrightarrow{E}\right]}_{ret}$$

(2)

and performing its Fourier transform on the detector

$$\overrightarrow{A}(\overrightarrow{r},\omega )=\frac{1}{\sqrt{2\pi }}\int \overrightarrow{A}(\overrightarrow{r},t)e^{-i\omega t}dt$$

(3)

the differential intensity over frequency $\omega$ and solid angle $\Omega$ is obtained as follows

$$\frac{d^{2}I}{d\Omega d\omega }=2{\left|\overrightarrow{A}(\overrightarrow{r},\omega )\right|}^2$$

(4)

where the dimensions of the last expression is [erg s]. In Eq.1 and Eq.2, the subscript $ret$ indicates that the field at point $\overrightarrow{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 $(\overrightarrow{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 2 a). In the case of relativistic particles, interpolating the trajectories is a delicate and potentially risky operation, as even small variations in quantities like $\overrightarrow{\beta }$ and $\dot{\overrightarrow{\beta }}$ can lead to huge discrepancies in energy or even result in unphysical situations (like superluminal velocities). The "forward" approach used in the development of 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 trajectory and detector at each instant, the fields calculated with Eq.1 will be associated with the time $t_d=t_t+R/c$ (Figure 2 b).

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 2 b, 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}\ne 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):

1.

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 done 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.

2.

Direct FT Calculation (DFC): In this faster approach, we exploit the principle of 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.

3. 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-10keV$. Therefore, a validation of the correct functioning of the code and the determination of its applicability limits is mandatory. As a validation approach, a well-known and analytically tractable case has been chosen: the undulator emission [5][10]. In the following first section, the numerically computed differential intensity in frequency and solid angle $d^{2}I/d\Omega d\omega$ 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.

3.1. Validation

It is possible to demonstrate that the differential intensity produced along the axis by a particle during propagation in an undulator takes on the following appearance (for further details, refer to Chapter 2 of [10]):

$$\frac{d^{2}I}{d\Omega d\omega }=\frac{e^2{N}^2{\gamma }^2}{c}\sum _{m=1}^{\infty }{F}_m\left(\xi \right)sin{c}^2({\nu }_m/2)$$

(5)

where $\xi$ is a variable defined from the undulator force, $F_m\left(\xi \right)$ is a function composed of Bessel functions of the first kind, ${\nu }_m$ is a linear variable of frequency $\omega$, N is the number of undulator periods and m is the harmonic number. Comparisons between the intensities calculated numerically using Eqs. 1-4 and Eq. 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/\gamma$, 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/\gamma$. Let’s assume that the transverse coordinate x as a function of the longitudinal coordinate z is approximately

$$x\left(z\right)=\frac{K}{k_U\gamma }cos\left(k_{U}z\right)$$

(6)

with K undulator strength and $k_U$ undulator wavevector. Then, the angle of the tangent to the trajectory will be given by

$$\theta =\frac{dx}{dz}=-\frac{K}{\gamma }sin\left(k_{U}z\right)\approx -\frac{K{k}_U}{\gamma }z$$

(7)

Calculating from Eq. 7, the $\Delta \theta$ as a function of the $\Delta z$ of the trajectory and equating it to $1/\gamma$, it is possible to determine the $\Delta t$ of the trajectory generating a visible signal on the detector by substituting the velocity $v=\beta c$:

$$\Delta z=\frac{1}{K{k}_U}\Delta t\approx \frac{1}{K{k}_U\beta c}$$

(8)

It has been verified that by choosing a timestep sufficiently smaller than the one indicated by the condition given by Eq.8 (e.g. $\Delta t_{sim}=0.1\Delta t$), the oscillation of the field on the detector is correctly reconstructed, and there is a coincidence between the theoretical and numerically calculated spectra. For larger $\Delta 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.

3.2. 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. The case involves a plasma density $n_p=2.44\times {10}^{18}c{m}^{-3}$ and a laser pulse energy $E_l=2.5J$. Other relevant laser pulse parameters, common to all simulations, include the wavelength ${\lambda }_l=800nm$, temporal length ${\tau }_l=25fs$, and beam waist $w_l=15\mu m$, assuming a perfectly Gaussian transverse spot. The laser pulse is polarized along the $\widehat{x}$ axis. In this instance, 10000 macroparticles were self-injected, and 600 timesteps of the PIC simulation were selected, requiring a computation time of $30min$ for a total of 100 different detector positions. Figure 5 visually represents the analysis of the presented case. In panel $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 $\mp \sigma$ (standard deviation) is shown in dashed lines. A clear collective oscillation of the beam along the $\widehat{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 resulting field intensity and photon flux.

In panel $5B$, the intensity of total radiation (integrated over all selected frequencies corresponding to the range ${10}^{-1}-{10}^{2}keV$) is represented in 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 what was mentioned about the trajectories in panel $5A$. In panel $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 panel $5A$. The subfigure on the top shows the average $\gamma$ 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 $\widehat{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 $\gamma$ reported in the subfigure. It shows a 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 panel $5D$, the spectra emitted on-axis by each of the three isolated chunks in panel $5C$ are shown, along with the numerical spectrum calculated from the analytical trajectory (orange). In the top-right inset, the total spectrum integrated over solid angle is presented. Moving from light blue to dark blue, three effects can be observed:

• Decrease in the cutoff energy.
• Decrease in the extension of the low-energy background (left).
• Increase in peak intensity.

The first effect can be attributed to the progressive reduction in the trajectories average amplitude: experiencing less acceleration (recall that transverse force is linear with distance from 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).

4. 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 Fourier transformed to retrieve the spectrum. The forward delayed fields calculation approach, presented in subSection 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 ad it’s currently used on a daily base in a parametric scan of laser-plasma parameters process. A systematic evaluation of 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.