RaDiO: an efficient spatiotemporal radiation diagnostic for particle-in-cell codes

TL;DR

This paper introduces RaDiO, an efficient algorithm for capturing detailed three-dimensional, time-resolved electromagnetic radiation emitted by charged particles in PIC simulations, preserving coherence and spectral accuracy.

Contribution

RaDiO is a novel radiation diagnostic algorithm that accurately models 3D, spatiotemporal electromagnetic fields in PIC codes, maintaining coherence effects and spectral fidelity.

Findings

Recovers theoretical radiation spectra accurately

Efficiently computes full spatiotemporal radiation features

Successfully applied to high harmonic generation in plasma mirrors

Abstract

This work describes a novel radiation algorithm designed to capture the three-dimensional, space-time resolved electromagnetic field structure emitted by large ensembles of charged particles. % in particle-in-cell (PIC) codes. The algorithm retains the full set of degrees of freedom that characterize electromagnetic waves by employing the Li\'enard-Wiechert fields to retrieve radiation emission. Emitted electric and magnetic fields are deposited in a virtual detector using a temporal interpolation scheme. This feature is essential to accurately predict field amplitudes and preserve the continuous character of radiation emission, even though particle dynamics is known only in a discrete set of temporal steps. Our algorithm retains and accurately captures, by design, full spatial and temporal coherence effects. We demonstrate that our numerical approach recovers well known theoretical radiated spectra in standard scenarios of radiation emission. We show that the algorithm is computationally efficient by computing the full spatiotemporal radiation features of High Harmonic Generation through a plasma mirror in a Particle-In-Cell (PIC) simulation.

Equations

$$\begin{array}{rl}\mathbf{E}(\mathbf{x},t_{det})& =e{\left[\frac{\mathbf{n}-\mathbf{\beta }}{{\gamma }_p^2(1-\mathbf{\beta }\cdot \mathbf{n}{)}^3{R}^2}\right]}_{\text{ret}}+\frac{e}{c}{\left[\frac{\mathbf{n}\times [(\mathbf{n}-\mathbf{\beta })\times \dot{\mathbf{\beta }})]}{(1-\mathbf{\beta }\cdot \mathbf{n}{)}^{3}R}\right]}_{\text{ret}},\\ \mathbf{B}(\mathbf{x},t_{det})& =[\mathbf{n}\times \mathbf{E}{]}_{\text{ret}},\end{array}$$

$$t_{det}=t_{ret}+|{\mathbf{r}}_{part}-R_{cell}{\mathbf{n}}_{cell}|/c,$$

$$t_{det}=t_{ret}+R_{cell}/c-{\mathbf{r}}_{part}\cdot {\mathbf{n}}_{cell}/c,$$

$$d{t}_{rad}=dt(1-v_p/c)\simeq dt/\left(2{\gamma }_p^2\right),$$

$$y(t)$$

$$x(t)$$

$$\frac{d^{2}I}{d\omega d\Omega }=\frac{e^2{\omega }^2{\gamma }^2}{3{\pi }^{2}c{\omega }_{\beta }K}{\left(\frac{1}{{\gamma }_p^2}+{\theta }^2\right)}^2\left[\frac{{\theta }^2}{{\gamma }_p^{-2}+{\theta }^2}{K}_{2/3}^2({\rm Y})+K_{1/3}^2({\rm Y})\right],$$

$$\frac{dI}{d\omega }=\sqrt{3}\frac{e^2{\gamma }_p\omega }{c{\omega }_c}{\int }_{\omega /{\omega }_c}^{\infty }{K}_{5/3}(x)dx,{\omega }_c=\frac{3}{2}K{\gamma }_p^2{\omega }_{\beta }$$

Full text

††journal: Computational Physics Communications

1 Introduction

Radiative processes in plasma are ubiquitous in astrophysics [1] and in laboratory settings. In plasma acceleration experiments, for example, they are important to the development of compact light sources [2], commonly employed in probing ultra-fast processes. Radiation emission mechanisms in plasma result from collective effects associated with the self-consistent dynamics of a large number of charged particles in the presence of strong electric and magnetic fields. Ab-initio numerical models, that can capture the motion of single particles, play an important role in this context, not only to validate theoretical advances, but also to predict radiation emission from experiments and in conditions where analytical models are not available.

Among the different numerical techniques, the Particle-in-Cell (PIC) [3] scheme provides a standard model to compute the motion of ensembles of charged particles. In its standard version, the PIC scheme consists in a loop that iteratively computes electric and magnetic fields by solving a discretized version of the full set of Maxwell’s equations in a grid, and then determines the next positions of the charged particles according to the relativistic Lorentz force. PIC codes are thus capable, by design, to retain most classical radiation emission processes.

The resolution required to capture radiation in the PIC algorithm poses quite stringent limitations on the shortest wavelengths that can be captured directly in a simulation, given that increasing the grid resolution will lead to a significant increase in the computational load. Consider a relativistic charged particle, with relativistic factor $\gamma_{p}$, undergoing a periodic motion with period T: The corresponding radiation wavelength, $\lambda_{rad}$, is proportional to $\lambda_{\mathrm{rad}}\propto cT/\gamma_{p}^{2}$ Hence, the spatial resolution required to capture $\lambda_{rad}$ is $\gamma_{p}^{2}$ times higher than the resolution needed to describe the particle trajectory. Furthermore, because of the Courant–Friedrichs–Lewy condition, the required temporal resolution is also $\gamma_{p}^{2}$ higher than standard. This results in an increase of $\gamma_{p}^{4}$ operations per simulation, pushing the limits of current computational capabilities, thereby motivating the development of advanced algorithms to compute radiation emission in PIC codes.

The standard approach to avoid the increased computational load and obtain high-frequency radiation emission from PIC simulations consists in performing additional radiation calculations outside the PIC loop using particle trajectory information obtained with the PIC algorithm. Many simulation codes have been developed over the recent years following this strategy. The code JRAD [4] receives a set of charged particle trajectories in order to compute the radiated spectra from the Fourier transform of the Liénard-Wiechert potentials; PIConGPU [5, 6, 7] follows a similar strategy, but can compute the emitted spectrum as the simulation progresses; the PIC codes OSIRIS [8] and EPOCH employ Monte-Carlo approaches to compute the spectrum of radiation from QED processes at run time (see, e.g. Ref. [9]). These tools have been successfully used to predict the radiation properties of laboratory plasmas (in plasma based accelerators [10]), Quantum Electrodynamics [11] and astrophysical plasmas.

However, the spatiotemporal profile of radiation is also important in fields such as astrophysics, where it can reveal the properties of rotating black holes [1, 12] for example. It can also play an important role in advanced microscopy based on twisted light with helical wavefronts [13]. Furthermore, this approach also provides a natural description of orbital angular momentum of light. To address this, we propose a new algorithm that retrieves the spatiotemporal radiation profile instead. This complementary approach includes built-in spatial and temporal coherence effects that are important to describe unexplored features of radiation emission, such as superradiant emission [14], for example. Our scheme can be used whenever the charged particle motion is well resolved, regardless of whether the spatial or temporal resolution is sufficient to resolve the resulting electromagnetic radiation.

The PIC simulation framework provides a direct and natural application to our present work and we focused on the implementation of this algorithm into the OSIRIS code naming our tool RaDiO, which stands for Radiation Diagnostic for OSIRIS. This diagnostic is composed of two distinct but equally useful counterparts: one implemented as a post-processing tool that uses previously generated trajectories to find the radiation that was emitted along them, and the other implemented as a run-time diagnostic for the PIC code OSIRIS, that uses the simulation data at each time step to compute the radiation.

This paper is structured as follows. In Section 2, we describe the theoretical framework behind radiation emission processes, which lays the groundwork for the development of the algorithm. Section 3 describes the implementation of the algorithm in detail, exploring key aspects like the temporal interpolation scheme. In Section 4, we benchmark our code against theoretical predictions and the results obtained with other radiation codes. Section 5 contains the study of the radiation emitted during the reflection of laser pulses by a plasma mirror. And, finally, Section 6 presents the conclusions.

2 Spatiotemporal electromagnetic field structure

The Fourier transformed Liénard-Wiechert fields [15] are commonly employed to predict the radiation spectra from charged particle trajectories. Here, instead, we calculate the Liénard-Wiechert fields directly, as these formulas provide the emitted electromagnetic fields at a certain position in space-time. The spatiotemporal Electric ($\mathbf{E}$) and Magnetic ($\mathbf{B}$) field structure of the radiation emitted by a charged particle according to the Liénard-Wiechert formulas is given by:

[TABLE]

with $\gamma_{p}=1/\sqrt{1-\beta^{2}}$. In Equation (1), the subscript ret denotes calculations using values at the retarded time, $\mathbf{n}$ is the unit vector oriented from the particle position to the region in space where we are interested in capturing the emitted radiation. The virtual region in space-time where radiation is deposited is henceforth denoted as the detector and will be described in more detail in Section 3. In addition, $\bm{\beta}=\mathbf{v}/c$ and $\dot{\bm{{\beta}}}=\dot{\mathbf{v}}/c$ are respectively, the particle velocity normalized to the speed of light, $c$ and the corresponding acceleration. Here the dot represents the time derivative. The direction of $\bm{\beta}$ and $\dot{\bm{\beta}}$ with respect to the virtual detector and $\mathbf{n}$ are schematically represented in Figure 1. Moreover $e$ is the electron charge and the quantity $R$ is the distance from the particle to the detector. For the purpose of determining the radiated fields, the first term in Equation (1) can be dropped if $R\gamma_{p}^{2}\dot{\beta}/c\gg 1$. This condition is usually satisfied in the far field ($R\gg c/\dot{\beta}$) for sufficietly relativistic particles ($\gamma_{p}\gg 1$). The second term in Equation (1) thus corresponds to emission of propagating electromagnetic waves, describing the so-called acceleration fields.

Equation (1) describes the emitted electric, $\mathbf{E}$, and magnetic, $\mathbf{B}$, fields at a given position, $x$ and time $t$, calculated from quantities obtained at the retarded time $t_{\text{ret}}$. For a given light ray that reaches the detector at a time $t_{\text{det}}$, $t_{\text{ret}}$ is the instant of time when emission has occurred. The time of arrival $t_{\text{det}}$ is given by:

[TABLE]

where $\mathbf{r}_{part}$ is the position of the particle and $R_{cell}\mathbf{n}_{cell}$ is the position of the detector’s cell. In order to enhance computational performance, it is useful and possible to simplify Equation (2) in the far field, which gives  [15]:

[TABLE]

Supplemented by the additional conditions given by Equations (2-3), Equation (1) can thus be used to retrieve the full set of spatiotemporal degrees of freedom of the radiation emitted by accelerated charges. By mapping the emitted radiation at each timestep in the particle trajectory to the corresponding time of arrival at the detector, the actual temporal resolution of the relativistic particle trajectory can be much coarser than the required one to describe the radiated fields.

An estimate of the maximum resolution that can be accurately obtained using Equations (2-3) can be found using the simplified picture shown in Figure 2: The particle located at $x_{0}$ emits a photon 1 at $t=t_{0}$. As the photon travels at c, in the next time step it will have travelled an extra $dt(\text{c}-v_{p})$ than the particle, which emits a second photon at $t=t_{1}$. Considering that a particle emits a photon at every time-step, the time interval between the arrival of two consecutive photons at the detector, provided that they are emitted by a relativistic particle, is given by Equation (4):

[TABLE]

with $dt$ being the temporal distance between emissions and the temporal resolution of the simulation providing for the particle trajectory.

Therefore, we are able to capture radiation with frequencies up to $2\gamma_{p}^{2}$ times larger than the ones used to sample the particle’s motion, as our detector time grid can be as fine as $dt_{\text{det}}=dt/2\gamma_{p}^{2}$. By consequence, the simulation time step can be much larger than the typical period of the emitted radiation. It is also important to note that the resolution in the detector should not be increased indefinitely as resolving time grids finer than $dt/2\gamma_{p}^{2}$ could generate non-physical information. A thorough analysis of these limits can be found in the Supplementary Material. The next section describes our implementation of the radiation algorithm and illustrates the reasons behind the different limits in resolution.

3 Algorithm and Implementation

The calculations of Equation (1) can be fully integrated either into a pre-existing code that computes the trajectories of charged particles (e.g. the PIC scheme) or be used as a post-processing tool that that computes Equation (1) on a set of pre-calculated trajectories. The algorithm consists of two main parts: calculating and obtaining the radiated fields and depositing them in a discretized grid. In this section we discuss the general steps and approach to incorporate the radiation algorithm considering these two components.

3.1 Radiation calculation algorithm

The virtual detector is a key feature of the radiation diagnostic. It is the region of space where radiation is tracked during a given time period. We consider two geometries of the virtual detector, (i) a spherical one [Figure 3 a)], where the grid is defined using spherical coordinates ($\mathbf{e_{\theta}},\mathbf{e_{\phi}},\mathbf{e_{r}}$) and (ii) a cartesian one [Figure 3 b)], where the grid is defined using cartesian coordinates ($\mathbf{e_{x}},\mathbf{e_{y}},\mathbf{e_{z}}$). RaDiO has the capability to compute the radiation in both types of geometries.

In order to track the emitted radiation at each time step of the trajectory we need to evaluate Equation (1) in every cell of the virtual detector. The radiation emitted at each time step of a given trajectory lies on a spherical shell that expands from the position of the particle at the time of emission, $t_{\text{ret}}$, at the speed of light. The intersection of the radiation shell with the detector consists of a circumference, whose radius increases with $t_{\text{det}}$. Figure 4 illustrates this picture, by showing the intersection of the radiation shell with a cartesian detector. The top of Figure 4 shows the detector at three different $t_{\text{det}}$. The bottom of Figure 4 shows the radiation arriving at each one of the highlighted cells as a function of $t_{\text{det}}$, which can be calculated using Equation (2) or Equation (3).

The illustration of Figure 4 suggests a clear approach to track the radiation reaching the detector from the emission of one particle at a given time step $t_{\text{ret}}$: Loop through each spatial cell of the detector and to compute $t_{\text{det}}$ at which radiation arrives. All the required quantities to compute Equation (1) are known or can be easily calculated (see additional details below). This approach features the quality of avoiding loops through the temporal cells in the detector. Thus, radiation computing time becomes independent from the temporal resolution of the detector and the total computing time is proportional to the number of time steps in the PIC simulation, $N_{\text{t\_PIC}}$ multiplied by the number of particles, $N_{\text{part}}$ multiplied by the number of spatial cells in the detector, $N_{\text{sp\_cell}}$.

This approach is summarized in Algorithm 1. It comprises two different loops: one through the particles that emit radiation (denoted as radiative particles) and another through the detector spatial cells. The quantities $t,\,R,\,\mathbf{n},\,\mathbf{\beta}$, $\dot{\beta}$ and $t_{det}$ are required in order to evaluate Equation (1). All of these quantities are either readily available or can be directly calculated from other quantities that are available in the simulation, such as the position of the particle ($\mathbf{x_{part}}$), the momentum of the particle ($\mathbf{p}$) and the time of emission $t$, as well as quantities that are part of the radiation module, such as the position of each detector cell $\mathbf{x_{cell}}$ or the previous velocity of the particle $\beta_{prev}$. These calculations are also shown in Algorithm 1.

Because $t_{\text{det}}$ can be computed at each time of emission, $t_{\text{ret}}$, using Equation (2) or Equation (3), it is in principle possible to conceive a temporally gridless detector. This approach could provide a very accurate description of the radiated fields, particularly if complemented by a post-interpolation scheme with the goal of retaining the continuous nature of radiation emission. Such approach, however, would require storing as many spatial detector arrays as the number of steps in the particle trajectory, for every particle in the simulation ($N_{\text{t\_PIC}}\times N_{\text{part}}\times N_{\text{sp\_cell}}$). High memory consumption would thus be the main limitation of such algorithm. To face this issue, RaDiO deposits radiation in a grid detector with up to 3 dimensions (1 temporal dimension and up to 2 spatial dimensions) with the spatial cells being distributed according to a spherical or cartesian geometry, and uses a temporal interpolation scheme to mimic continuous radiation emission between two consecutive PIC time-steps for every particle in the simulation.

The implementation shown in Algorithm 1 can be applied to both post-processing diagnostics, which calculates the radiation given a set of pre-calculated trajectories, and to run-time diagnostics, in which the radiation calculations are performed at run time during the trajectory calculation. In the latter scenario, the calculation and deposition of the emitted radiation can take place in a sub-step of the particle push loop, created specifically for that purpose. This sub-step comes right after pushing the particles, in such a way that the newly calculated positions and momenta can be used, in conjunction with the corresponding stored values from the previous iteration, to compute the required quantities to determine the radiated fields. In the post-processing version, all required quantities can be readily calculated by considering the positions and momenta from consecutive time-steps.

3.2 Deposition of the radiated fields in a virtual detector

According to Eqs. (2) and (3), each PIC simulation timestep corresponds to a given detector time. In general, consecutive time steps in the trajectory will deposit radiation in non-consecutive detector time cells. A simple prescription that only deposits the radiated fields in the temporal cells that are closest to the predictions given by Eqs. (2) and (3) will therefore generate noisy radiation patterns that are non-physical because particles emit radiation continuously. To re-gain the continuous character of radiation emission, and remove the artificial noise induced by the discretization of the trajectories in time, RaDiO interpolates the fields emitted by each particle between every two consecutive PIC time steps.

The interpolation scheme in RaDiO assumes that particles radiate constant fields between each consecutive PIC timestep. In order to deposit the fields across different temporal cells, we weigh the contribution of each deposition by the time until the next deposition. In fact, the value of the radiation in a time slot is the integral of the radiation in the interval delimited by two consecutive detector time-steps. Incidentally, real-life applications often employ an integrator detector, which takes the information about radiation arriving in-between detector time steps into account. This deposition scheme can be implemented by following Algorithm 2, below.

Each variable in Algorithm 2 is calculated at each PIC time-step and for each particle. Here, slot(…) is a function that returns the index of the slot in the detector’s time-array ($t_{\text{array}}$) where $t_{\text{det}}$ falls, $t_{det}$ is the time of the current deposition and $n_{\text{slot}}$ is the corresponding time-slot position in the detector array. In addition, $n_{\text{itr}}$ is an iterator that runs from $n_{\text{slot,prev}}$, the detector time slot where particle deposited radiation in the previous PIC time-step, until $n_{\text{slot}}$. The quantity $t_{\text{tmp}}$ is an auxiliary variable for the calculation of the time difference between depositions. It runs from, $t_{\text{det,prev}}$, the time of the previous deposition, to $t[n_{\text{itr}}]$, the time for the actual deposition.

Figure 5 shows an example case that clarifies this deposition scheme. Each of these depositions correspond to radiation emitted at a different PIC time step by a single particle. This interpolation can be performed while the simulation is running, as it only requires information about the radiated field in the previous time step. In fact, for the example present in Figure 5 the deposition algorithm would go as follows:

1. At PIC iteration 4, radiation arrives at the detector at $t_{\text{det}}=t^{\prime}_{4}$.

2. $n_{\text{itr}}$ is set to 2, the slot of the previous deposition, at $t^{\prime}_{3}$, $t_{\text{tmp}}$ is set to $t^{\prime}_{3}$, the time of the previous deposition, we enter the loop, the scale factor is calculated: $(t_{3}-t_{\text{tmp}})/dt_{\text{det}}$, with $t_{\text{array}}[n_{\text{itr}}+1]=t_{3}$ and $E(t^{\prime}_{3})(t_{3}-t^{\prime}_{3})/dt_{\text{det}}$ is deposited in the second time slot, $t_{2}$.

3. $n_{\text{itr}}$ is incremented to 3, $t_{\text{tmp}}$ is set to $t_{3}$, the time of the previous deposition, the scale factor is calculated: $(t_{4}-t_{3})/dt_{\text{det}}$, with $t_{\text{array}}[n_{\text{itr}}+1]=t_{4}$ and $E(t^{\prime}_{3})(t_{4}-t_{3})/dt_{\text{det}}$ is deposited in the third time slot, $t_{3}$.

4. $n_{\text{itr}}$ is incremented to 4, $t_{\text{tmp}}$ is set to $t_{4}$, we exit the loop, the scale factor is calculated: $(t_{\text{det}}-t_{4})/dt_{\text{det}}$ and $E(t^{\prime}_{3})(t^{\prime}_{4}-t_{4})/dt_{\text{det}}$ is deposited in the time slot $t_{4}$.

Using this approach, radiation can be computed and deposited using only the information from the current and the previous time steps. This algorithm interpolates radiation coming from a single particle, but can be repeated for all particles in the simulation, as stated in Algorithm 1, in order to capture radiation from all particles.

3.3 Practical example: Helical trajectory

Here we look at a practical example, in which an electron with $\gamma_{p}=57.3$ undergoes an helical motion with amplitude $0.014\,\text{c}/\omega_{p}$ and frequency $\omega_{0}=\omega_{p}$, corresponding to a $K$ parameter of $K=0.8$, $K$ is a trajectory parameter that can be taken as a scaled pitch angle the maximum angle of the particle trajectory, normalized to the Lorentz factor $\gamma_{p}$ and given by $K=\gamma_{p}r_{0}\omega_{0}/c$. The helical motion was described by the PIC algorithm with a temporal resolution of $0.1\;\omega_{p}^{-1}$. Here, $\omega_{p}$, is an arbitrary normalizing frequency. The radiation generated by a particle undergoing such trajectory has a distinctive, expanding spiral spatiotemporal signature. This is shown in Figure 6, which represents the radiated electric field along the $y$ direction deposited onto a spherical 2D detector with an angular aperture of $0.1$ rad placed in the direction of the longitudinal motion of the particles ($x$ axis), with radius $R=10^{6}\;\text{c}/\omega_{p}$. The temporal resolution of the detector was $1.33\times 10^{-5}\;\omega_{p}^{-2}$ and the spatial resolution was $58\;\text{\textmu}\text{rad}$. Figure 6 shows a snapshot of the detector at four different temporal positions. The starting point of the spiral follows the circular motion of the particle in the $y-z$ plane. Between each snapshot the radiation spiral makes two turns, thus the temporal distance between each snapshot is approximately equal to two periods of the emitted radiation. Given the trajectory parameters, the radiation period was expected to be of about $\sim 2\times 10^{-3}\;\omega_{p}^{-1}$, about 10 times smaller than the smaller period that could be resolved using only the PIC algorithm.

4 Benchmarking

In order to benchmark our algorithm, we consider the example of a single relativistic particle emitting synchrotron radiation. Synchrotrons have a magnetic field structure that imposes a sinusoidal trajectory to relativistic electrons that go through the device, thus leading to the emission of high frequency photon beams in the X-UV or X-ray regions of the spectrum. The trajectory of the particle would then be given by:

[TABLE]

where $\beta_{x0}$ is the initial velocity of the particle along the longitudinal $x$ direction, ${\gamma_{x0}=(1-\beta_{x0}^{2})^{-1/2}}$ its longitudinal Lorentz factor, $r_{\beta}$ the amplitude of the sinusoidal trajectory, and $\omega_{\beta}$ its frequency.

As far as we are aware, the only explicit analytical formulas capturing the spatiotemporal radiation profile of synchrotron radiation are found in [16], which gives a semi-analytical model for the emitted field lines. However, direct quantitative comparisons between the visual depiction of field lines and the actual value of the emitted field in a region of space can be difficult (see Supplementary Material for a qualitative comparison). On the other hand, the spectral properties of radiation are well documented [17, 18], so we Fourier transformed the data in the virtual detector with respect to time and compared these spectra to the theoretical predictions. The corresponding intensity spectrum ($I$) with respect to the frequency $\omega$ and solid angle $\Omega$ of the emitted radiation , valid for ultra relativistic particles as an asymptotic limit expression ($\gamma_{p}\gg 1$) and assuming very large number of periods in the trajectory,  [2], is given by:

[TABLE]

where, $\theta$ is the observation angle in the direction perpendicular to the trajectories plane. In addition, $K_{n}$ is the modified Bessel function and $\Upsilon$ is a numerical parameter given by $\Upsilon=\frac{\omega\gamma_{p}}{3\omega_{\beta}K}\left(\gamma_{p}^{-2}+\theta^{2}\right)^{-3/2}$, with $K$ being the aforementioned $K$ parameter.

Equation (7) can be integrated over all angles, returning the frequency spectrum [2]:

[TABLE]

We have benchmarked our algorithm against Equations (7) and (8). The benchmarks were performed using the two dimensional sinusoidal trajectory of a relativistic electron ($\gamma_{p}=50$) with an amplitude of $r_{\beta}=2\,c/\omega_{p}$, $k_{\beta}=0.1\,\omega_{p}/c$ ($K=10$) in the transverse $x-y$ plane and $dt=0.01c/\omega_{p}$, where $\omega_{p}$ is a normalizing frequency. We simulated a line of a spherical detector, placed in the $z-x$ plane, $10^{5}\,c/\omega_{p}$ away from the axis origin with an angular aperture of $0.1$ rad around the $x$ axis. This detector had 512 spatial cells and 131072 temporal cells, resulting in a temporal detector resolution of $2.98\times 10^{-5}c/\omega_{p}$ The results are shown in Figure 7 which features a plot of the detected electric field in the $\mathbf{e_{\phi}}$ direction (perpendicular to the motion plane) for each spatiotemporal cell.

The radiation is composed of several periodically spaced peaks, whose shape can be observed in the lineout [Figure 7 b)]. The short burst nature of the radiation (equivalent to a broad band spectrum), consistent with the large value of the K parameter, is clear from Figure 7. Instead of displaying a purely sinusoidal profile with a single wavelength, the electric field consists of sharply peaked bursts containing many different wavelengths. Moreover, it is possible to observe that consecutive peaks have opposite sign. This is a direct result of the sinusoidal nature of the electron trajectory in which the acceleration $\dot{\bm{\beta}}$ switches sign between peaks. Furthermore, it is possible to note that for higher angles the radiation bursts arrive later, creating the parabola-like structures that can be seen in the upper plot. This delay becomes more significant as the particle approaches the detector’s surface, resulting in a decrease of the curves’ aperture.

This result can also be understood in terms of the spatiotemporal reasoning regarding the estimation for the typical radiation frequency presented in the previous section. However, instead of depicting the emitted radiation parallel to the motion of the particle, we picture them emitted at an angle $\theta$. The temporal distance between the emission and arrival of light ray emitted at a given longitudinal position $x$ is then given by: $c\Delta t_{rad}=\sqrt{R^{2}+x^{2}-2xR\cos{\theta}}$, where all quantities are defined as in Equation (3). This expression shows that the time of arrival increases with $\theta$ and also that it is scaled by the longitudinal position $x$. Thus, as the particle approaches the detector and $x$ grows larger, the parabolic structures left on the detector become tighter.

Figure 7 b), which depicts lineouts of $E_{\phi},$ also shows that the peaks become wider and less intense for larger angles. This is in concordance with the predictions for the spectrum [see Equation (7)], which features a decrease in the number of harmonics for larger angles, resulting in broader and less intense peaks off-axis.

In order to further understand the angular dependent frequency spectra, Figure 8 compares the theoretical result, given by Equation (7), with the simulated result, given by the Fourier transform over time of the field shown in Figure 7 a). The spectrum is symmetric with respect to $\theta=0$. Thus, the upper half of Figure 8 a) ($\theta>0$) shows the simulated results and the bottom half ($\theta<0$) the theory. As expected, the theoretical line, being the assymptotic limit of a continuous harmonic distribution with a very large number of oscilations in the trajectory [2], corresponds to the envelope of the numerical result, showing excellent agreement. This is evident from the lineout of the radiated spectra displayed in Figure 8 b).

The simulated integrated spectrum over all angles, which yields the frequency distribution of the emitted radiation, can be benchmarked against Equation (8). Figure 8 shows excellent agreement between numerical and theoretical results, as the intensity of most peaks matches the expected result with small relative error which rises as frequency increases.

To further confirm the validity of our numerical approach, we benchmarked the frequency integrated spectrum, $dI/d\Omega$ against the spectrum provided by the post-processing spectral code JRad [4], which computes the radiated fields using the spectral version of the Liénard-Wiechert potentials. The results of this comparison, shown in Figure 8 d), are in excellent agreement.

4.1 Coherence tests

Because RaDiO captures the emitted fields in space and in time it can also naturally describe temporal and spacial interference effects. This feature is essential to accurately portrait temporal and spatial coherence, present in superradiant emission scenarios for example. This is an intrinsic feature of our spatiotemporal approach, which allows us to directly obtain the fields radiated by every simulation particle, including interference effects by design.

To test our ability to acurately model temporal and spatial coherence, we ran simulations using two particles with opposite charges and sinusoidal trajectories, similar the one defined in Equations (5) and (6). The two particles, with particle 1 being positively charged and particle 2 being negatively charged underwent this sinusoidal trajectory in perpendicular planes (particle 1 in plane $x-y$ and particle 2 in plane $x-z$) and the detector was the same as the one used in the previous section.

Figure 9 shows the simulated radiated electric field profile as a function of $\theta$ and $t_{\mathrm{det}}$ for three different configurations: one with only particle 1 [Figure 9 (a)], one with only particle 2 [Figure 9 (b)] and other with both particles [Figure 9 (c)]. As the two trajectories lie in different planes, the spatiotemporal signatures of the radiation emitted by each particle are noticeably distinct because the detector plane lies on the plane of the trajectory of particle 2, then being perpendicular to the plane of particle 1. By comparing Figure 9 (a) with Figure 9 (b), we can hence readily identify the radiation coming from each particle in Figure 9 (c).

As both particles have opposite charges, the field on axis for a given particle will have the opposite sign as the on-axis field for the other particle. Thus, the radiation emitted by both particles will interfere destructively on-axis. This happens exactly at $\theta=\pi/2$. Thus, if we look at the time averaged squared field (insets in each panel of Figure 9), we see that although $\langle E^{2}\rangle_{t}$ is maximum at $\theta=\pi/2$ for the simulations with only one of the particles (insets of Figure 9 [a] and Figure 9 [b]), the opposite happens when we capture the fields radiated by both particles (inset of Figure 9 [ c]).

Our algorithm captures coherence effects of the simulation particles by default, but in a PIC code, each particle in the simulation represents a cloud of $N$ real particles with a size close to cell size that follow the same dynamics, this is the so called macroparticle aprroximation. In our code, however, we calculate the radiation emitted by the macroparticles in the simulation as if they were point charges with charge equal to the total charge inside the macroparticle ($Nq$). This is in fact equivalent to assuming that each of the $N$ particles inside the macroparticle radiates coherently. The assumption that they all radiate coherently holds either for all wavelengths if $N=1$, or for wavelengths larger than the cell size if $N\gg 1$. For wavelengths shorter than the cell size, in general, we cannot say it holds, as such an assumption depends on information about particles that are not being simulated. For example, if standard macroparticle approximation is still valid at scales smaller than the cell size, the emitted radiation should be incoherent for wavelengths shorter than the cell size and the result should be corrected with a filter function (see Supplementary Material for a deeper analyisis).

The detailed study of the conditions that allow assuming that each of the $N$ particles inside the macroparticle radiates coherently is out of the scope of this work. It will be up to the user to decide whether it holds or not. If this assumption does not hold, then results given by our code will be correct for wavelengths larger than the cell size, but could be overestimated for wavelengths smaller than the cell size. Nevertheless our code can, in general, accurately predict the qualitative aspects of the emitted radiation for all wavelengths.

5 Example: Radiation from a plasma mirror

When an electromagnetic wave collides with a target such as a metallic surface or an overdense plasma, it is unable to propagate and gets reflected. The process of reflection has long been well understood and thoroughly explained at the macroscopic level by Maxwell’s laws and classical electrodynamics. In the plasma, the phenomenon is commonly explored using a fluid theory approach. Such description predicts the damping of the wave near the surface of the reflective material (it becomes an evanescent wave) and the appearance of a reflected wave. At the electron level, however, the phenomenon is not always trivial, in particular at relativistic laser intensities (with peak normalized vector potential $a_{0}=eA_{0}/(m_{e}c)=1$ ), which lead to High Harmonic Generation (HHG) in plasma mirrors [19, 20]. Several theoretical frameworks have been proposed to describe the underlying mechanisms of HHG, each with different regimes of applicability (see e.g. [21, 22])

PIC simulations are commonly employed to deepen the understanding of the physical processes underlying laser reflection and harmonic generation in plasma mirrors. An accurate description of HHG in standard PIC simulations, for instance, is computationally challenging because spatial and temporal PIC grids need to properly resolve the high harmonics. Thus, to accurately capture high harmonics up to the $10^{th}$ or $100^{th}$ order, PIC simulations require spatiotemporal resolution up to one-two orders of magnitude higher than one required to resolve the fundamental harmonic. The use of RaDiO may thus be computationally advantageous in HHG simulations, as it allows capturing high frequency harmonics without increasing the PIC resolution.

In this section, we present 3D Osiris simulations of an HHG scenario where the laser propagates in the longitudinal $x$ direction and is linearly polarized along the transverse $z$ direction. The laser uses a $\sin^{2}$ temporal profile with 12 full periods ($T_{0}=8\pi~{}\omega_{p}^{-1}$, $\omega_{0}=\omega_{p}/4$) and Gaussian perpendicular profile with spot-size ($W_{0}=2\lambda_{0}$, $\lambda_{0}$ is the central laser wavelength). The plasma mirror consists in an overdense plasma slab with plasma frequency $\omega_{p}$ (and density $n_{p}$, $16$ times larger than the critical density $n_{c}$ for that laser pulse) with thickness 100 $c/\omega_{p}$, much higher than the non-relativistic plasma skin-depth ($l_{s}\sim~{}c/\omega_{p}$ in this case). As the laser gets reflected, we capture the reflected fields both in the PIC grid through Maxwell’s equations and in a virtual detector through RaDiO. We chose to compute the radiation emitted by all plasma electrons located within the plasma cylinder with a radius of three laser spot sizes around the focus. The virtual cartesian detector was located at $x=-160~{}c/\omega_{p}$, ranging from $y=-160~{}c/\omega_{p}$ to $160~{}c/\omega_{p}$, with temporal resolution $dt_{\text{det}}=0.0384~{}\omega_{p}^{-1}$, about five times smaller than the PIC temporal resolution $dt_{\text{PIC}}=0.1792~{}\omega_{p}^{-1}$. The PIC simulation box ranged from $x=-288~{}c/\omega_{p}$ to $x=108~{}c/\omega_{p}$, with a resolution $dx=0.96~{}c/\omega_{p}$ in the longitudinal direction and from $y,z=-160~{}c/\omega_{p}$ to $y,z=160~{}c/\omega_{p}$ with resolution $dy,dz=0.32~{}c/\omega_{p}$ in the transverse direction. This PIC grid is able to resolve 26 points per laser wavelength. Each cell contains 16 simulation particles. A 2-D slice of the setup is shown in Figure 10, the laser propagates from left to right.

We start by capturing the radiation in the absence of HHG, by using a non-relativistic laser intensity, with peak normalized vector potential $a_{0}=0.1$. Figure 11, top, shows the trajectories of a random sample of 512 plasma particles. The zoomed-in region clearly displays the typical figure-8-like motion induced in the plasma particles by the laser pulse. This motion originates the radiation, which is captured both in the PIC grid and in the virtual radiation detector. By comparing the radiation in the detector to the reflected pulse in the PIC grid (Figure 11, bottom), we show that the beam reflection is a direct result of the charged particles’ trajectories induced by the incident beam.

Next to investigate a scenario with strong HHG, we used a high-intensity laser ($a_{0}=4.2$) in a setup similar to the one shown on Figure 10. In this case, we see the clear effect of the increased intensity on the trajectory of the sampled particles (Figure 12, top), with a similar figure-8 motion for the first few laser periods, but with increased amplitude overall and stronger deviation from the standard figure-8 motion.

As a result of this more extreme motion, the reflected laser beam is noticeably different from the incident beam. This is made clear in the comparison shown at the bottom of Figure 12. The differences between incoming and reflected laser pulse electric field profile are due to the existence of high laser harmonics present in the reflected beam. The presence of the high harmonics is also clearly visible in the spectrum of Figure 12. The frequency spectrum shows that the reflected laser captured by RaDiO contains at least 13 harmonics, while the PIC algorithm, which only resolves the plasma relevant scales correctly captures the emission of the first 4 odd harmonics. The PIC grid is able to resolve the original harmonic with 26 points per wavelength, but as the harmonic order increases, past the 7th order only RaDiO’s resolution can capture the signal correctly. In this case the RaDiO frequency spectrum captures frequencies at least 4 times higher than the OSIRIS PIC grid, being able to capture harmonics at least until the 25th order, as expected from the employed laser intensity ($a_{0}=4.2$).

6 Conclusions

The radiation diagnostic for OSIRIS (RaDiO) was successfully implemented, benchmarked and tested in several scenarios, including production runs. While not described here, it should be also be noted that the algorithm was fully parallelized allowing for large simulations. RaDiO is a novel radiation diagnostic that captures the spatiotemporal features of high frequency radiation in PIC codes. A key aspect of our algorithm is the development of a temporal interpolation scheme for depositing radiation. This is essential to preserve the continuous character of radiation emission and to obtain correct values for the amplitude of the radiated fields. The algorithm is general and only requires knowledge about the trajectories of a an arbitrarily large ensemble of charged particles ($>10^{6}$) thus we can apply it to generally enhance the capabilities of any algorithm that predicts the trajectories of charged particles, apart from PIC codes. We described the implementation of RaDiO into OSIRIS and provided benchmarks with well established theoretical models for synchrotron emission. These comparisons showed excellent agreement, therefore adding a high level of confidence to future runs.

We also provided an illustration where we used RaDiO to probe the spatiotemporal features of radiation emitted in the context of laser reflection by a plasma mirror. At lower laser intensities, RaDiO fully recovers the PIC simulation result. This further confirms the validity of RaDiO in a setting where temporal and spatial coherence effects are critical. A simulation at higher laser intensity demonstrated the generation of high harmonics beyond the predictions of the PIC algorithm, showing that RaDiO allows for a complete characterization of the reflected beam along with all the harmonics, without increasing the overall PIC resolution, and effectively demonstrating that RaDiO can be effectively used to predict high frequency radiation from PIC codes [14].

RaDiO is a flexible diagnostic tool that can can be further expanded to include additional features such as higher order interpolation schemes, for example using an advanced particle pusher recently developed [23], the option to compute the electromagnetic field potentials in addition to the electromagnetic fields, or the capability to convert radiation to/from relativistic Lorentz boosted frames. Although this diagnostic does not interact with the particles, it could also be employed together with a QED code that captures radiation reaction and affects the particle’s trajectories and capture radiation compatible with QED effects as long as the emission is purely classical. Because it captures the radiation in space and in time, RaDiO may also be useful in describing the production of spatiotemporally structured beams [24].

Acknowledgments

We acknowledge the Partnership for Advanced Computing in Europe (PRACE) for access to the Leibniz Research Center on SuperMUC and the Barcelona Supercomputing Center on Marenostrum 4. This work was partially supported by the EU Accelerator Research for Innovation for European Science and Society (EU ARIES) under grant agreement no. 738071 (H2020-INFRAIA-2016-1). JV acknowledges the support of FCT (Portugal) Grant No. IF/01635/2015/CP1322/CT0001 and MP acknowledges the support of FCT (Portugal) Grant No. PD/BD/150411/2019.