PICSAR-QED: a Monte Carlo module to simulate Strong-Field Quantum Electrodynamics in Particle-In-Cell codes for exascale architectures

Physical scenarios where the electromagnetic fields are so strong that Quantum ElectroDynamics (QED) plays a substantial role are one of the frontiers of contemporary plasma physics research. Investigating those scenarios requires state-of-the-art Particle-In-Cell (PIC) codes able to run on top high-performance computing machines and, at the same time, able to simulate strong-field QED processes. This work presents the PICSAR-QED library, an open-source, portable implementation of a Monte Carlo module designed to provide modern PIC codes with the capability to simulate such processes, and optimized for high-performance computing. Detailed tests and benchmarks are carried out to validate the physical models in PICSAR-QED, to study how numerical parameters affect such models, and to demonstrate its capability to run on different architectures (CPUs and GPUs). Its integration with WarpX, a state-of-the-art PIC code designed to deliver scalable performance on upcoming exascale supercomputers, is also discussed and validated against results from the existing literature.


Introduction
One of the frontiers of modern physics research deals with physical scenarios where the electromagnetic fields are so strong that Quantum ElectroDynamics (QED) plays a substantial role, specifically in the so-called strong-field regime (sf-QED). These scenarios range from the interaction region of an ultra-intense laser pulse with a plasma [1][2][3][4][5][6][7][8][9][10] to extreme astrophysical objects, such as pulsar magnetospheres [11,12], black-holes [13][14][15][16], or gamma-ray bursts [17]. Numerical modeling is essential to gain insights into these scenarios and to assist the experimental investigation of ultra-intense laser-matter interaction. Since a kinetic description of the plasma is usually required, Particle-In-Cell (PIC) codes [18][19][20] are often the numerical tool of choice. Moreover, PIC simulations can include the most relevant sf-QED effects in these regimes [21][22][23][24][25], such as the emission of high-energy photons via the inverse Compton process [26] (also known as nonlinear synchrotron emission) and the decay of a high-energy photon into an electron-positron pair via the nonlinear Breit-Wheeler pair production process [26][27][28][29][30] (see Fig. 1 for a scheme showing the core algorithms of a PIC code including sf-QED effects).
Realistic simulations often demand a substantial amount of computational resources, to the point that the most ambitious numerical campaigns can only be performed on the most powerful High-Performance Computing (HPC) facilities [31]. Most of those machines offload a conspicuous fraction of their calculations to specialized hardware [32] (i.e., Graphics Processing Units, GPUs) or make use of CPUs specifically designed for HPC needs [33]. Few PIC codes in use in the plasma physics community can efficiently take advantage of those machines. Moreover, only a fraction of these codes is distributed as free and open-source software. Requiring the capability to simulate at least the most relevant QED processes further restricts the choice. OSIRIS [34], Picador [35], VPIC 2.0 [36], and PIConGPU [37] are popular codes able to take advantage of modern, GPU-based supercomputers. However, only PIConGPU and VPIC 2.0 are available as free and open-source software. Moreover, while OSIRIS and Picador have very comprehensive QED modules, PIConGPU doesn't implement nonlinear Breit-Wheeler pair production yet, and VPIC 2.0 doesn't implement QED processes. Smilei [38], EPOCH [20], and Tristan-MP v2 [39] are well-known opensource, massively parallel PIC codes with comprehensive QED modules, but they are not currently designed to take advantage of GPU-based supercomputers. Calder [23,40] is also well-known for its very comprehensive QED modules [25,41], but it is not Figure 1: The figure shows the essential sub-steps of a timestep of a PIC code. A gold border marks the sub-steps related to QED effects. In the core PIC algorithm there are two main actors: 1) the plasma, which is simulated with a collection of macroparticles (each representing several real particles) and 2) the electromagnetic field, which is simulated on a grid. In each timestep, charged particles move according to the fields interpolated from the grid. The motion of these charged particles generates currents on the grid, which are a source for the evolution of the electromagnetic field. In order to simulate sf-QED two additional sub-steps are added. Particles undergoing QED effects have an additional property, called optical depth. This quantity is evolved at each timestep. When the optical depth becomes negative, an event such as Breit-Wheeler pair production or nonlinear synchrotron emission occurs. See [20,24] and section 3 for more details.
open-source nor optimized for large-scale GPU-based HPC machines. This paper presents PICSAR-QED [42], a module part of the PICSAR library [43], which has been coupled with the WarpX [44,45] PIC code in order to simulate sf-QED processes relevant for extreme plasma physics scenarios. WarpX is an opensource code developed within the framework of the Exascale Computing Project [46] and designed to provide scalable performance on upcoming exascale supercomputers. PICSAR-QED implements the methods needed to model sf-QED processes in PIC codes, and WarpX+PICSAR-QED is -to the best of the authors' knowledge -the first open-source PIC code able to simulate sf-QED processes on large-scale GPU-based supercomputers. PICSAR-QED implements primitives designed to be portable across different architectures (CPUs and GPUs). It's conceived to be easily included in existing projects, and is released as a standalone open-source project, contributing a carefully validated module to the open-source PIC ecosystem, of which other PIC codes not currently implementing sf-QED could take advantage. In this work, we first review the physical models implemented in PICSAR-QED (section 2). Then, we review the numerical methods to model sf-QED processes in PIC codes, we discuss the specific implementation choices made in PICSAR-QED, and we present a detailed validation of the methods provided by the library (section 3). In section 4 we show performance benchmarks on different architectures. Finally, we describe in section 5 the integration of PICSAR-QED in WarpX and present benchmarks with existing results from the literature.

Physical processes implemented in PICSAR-QED
The extreme plasma physics scenarios mentioned in the introduction are characterized by electromagnetic fields so strong that relevant QED processes pertain to the strongfield regime of QED. Strong-field here refers to the electromagnetic fields being comparable to the QED critical field E s , also known as Schwinger field [47][48][49]: where m e is the electron mass, c is the speed of light, q e is the elementary charge and is the reduced Plank constant. Such a tremendous field is far beyond current experimental capabilities, being roughly three orders of magnitude higher than the strongest electric fields available on Earth [50]. However, it can be approached in the reference frame of ultra-relativistic particles. Indeed, in this case, the actual parameter of interest for sf-QED is the parameter χ, defined as the ratio between the electric field in the reference frame of the particle and E s . χ is Lorentz-invariant and is called quantum parameter.
For an electron or a positron with 4-momentum p µ propagating in a region where the electromagnetic field tensor is F µν , χ is defined as follows [29] : where E is the electric field, B is the magnetic field, v is the velocity of the particle, and γ is its Lorentz factor. We attain the so-called full quantum regime of sf-QED when χ > 1, while strong-field QED effects rapidly vanish for χ < 1. For high-energy photons with 4-momentum p µ , an analogous Lorentz-invariant χ γ parameter can be defined: where γ γ is the photon energy normalized to m e c 2 and c is a velocity vector with a magnitude equal to the speed of light. The threshold of the full quantum regime of sf-QED is χ γ ∼ 1 for photons as well. The physical models used to include the most relevant sf-QED processes in PIC codesnonlinear Breit-Wheeler pair production and inverse Compton emission -are well known from the literature [21, 22, 24-28, 30, 51] and are briefly reviewed below. Schwinger pair production, which is a particularly extreme physical process where the electromagnetic field is strong enough to generate electron-positron pairs from the quantum fluctuations of the vacuum, is also briefly discussed.

Nonlinear Breit-Wheeler pair production
Nonlinear Breit-Wheeler pair production is the decay of a high-energy photon propagating in a strong field into an electron-positron pair. The differential pair production cross section for a photon with quantum parameter χ γ reads as follows [23,27]: where χ + and χ − are respectively the quantum parameter of the emitted positron and of the emitted electron, K α are the modified Bessel functions of the second kind of order α, α f is the fine structure constant, and Since particles in very intense electromagnetic fields are usually ultra-relativistic (i.e., γ 1), the so-called ultra-relativistic approximation can be used. In this approximation, the 3-momenta of the product particles are aligned with that of the highenergy photon. Moreover, since |v| ≈ c, the square root terms in the expressions for χ and χ γ become almost identical. Finally, since γ ≈ |p|/(m e c), using total momentum conservation we can write: This result is used to replace χ + with χ γ − χ − , so that the differential cross section and the total cross section can be rewritten as follows: where It is noteworthy that very good and simple asymptotic approximations exist for T (χ γ ) [27]: Equation (8) allows determining the probability of a photon to decay into an electronpositron pair within a given time interval. Indeed, if the timestep of the simulation is ∆t, this probability is p decay = 1 − exp(−∆tdN/dt). In a Monte Carlo approach, a random number r decay from a uniform probability distribution between zero and one can be drawn and compared with p decay : if r decay < p decay Breit-Wheeler pair production occurs for the given photon (or, in a PIC simulation, a macro-photon, that is a numerical particle representing several real photons). In practice, however, a different approach is typically followed [22], in order to avoid a random number extraction per particle at each iteration. In this approach, each macro-photon has a randomly initialized quantity, called optical depth, which is reduced at each time-step according to the total cross section. As soon as this quantity reaches zero, Breit-Wheeler pair production occurs (see section 3 for a more in-depth discussion). The energy of the generated particles can be determined using Eq. (7), by calculating the cumulative probability distribution with respect to χ − : The quantum parameter of the electron χ − can be sampled by solving P (χ γ , χ − ) = r, where r is a random number drawn from a uniform probability distribution in the range [0, 1]. The quantum parameter of the positron is then simply χ + = χ γ − χ − . In the ultra-relativistic limit, determining the energy and momenta of the generated particles is straightforward if the quantum parameters χ ± are known. Indeed, their kinetic energy can be calculated as K ± = (γ γ − 2)χ ± /χ γ .

Inverse Compton photon emission
The inverse Compton photon emission is the emission of a high-energy photon from a charged particle (e.g., an electron or a positron) propagating in a strong electromagnetic field. PICSAR-QED implements the model described in [22], which is summarized here for completeness (it is worth noting that a slightly different notation is adopted: in [22] χ is replaced with η and χ γ is replaced with 2χ). The differential cross-section for the inverse Compton scattering process reads as follows: where Y (χ, χ γ ) = 2 3 Equation (12) can be re-written as where we introduced ξ = χ γ /χ . As for Breit-Wheeler pair production, the ultrarelativistic approximation applies. Therefore, since the photon is emitted within a cone of amplitude α ≈ 1/γ, for γ 1 the photon can be safely considered to be emitted along the direction of the momentum of the emitting particle. Within the ultra-relativistic approximation, it is also trivial to show that ξ < 1. The total cross section is obtained by integrating Eq. (14) over χ γ from 0 up to χ: where Equation (17) allows determining the probability of an electron or a positron to emit a high-energy photon via inverse Compton emission, with a procedure identical to that described for Breit-Wheeler pair production. As for Breit-Wheeler pair production, the quantum parameter χ γ of the generated photon is determined using the cumulative probability distribution: Once χ γ is known, the energy of the generated photons can be determined trivially using the ultra-relativistic approximation: γ γ = (γ − 1)ξ. Finally, the kinetic energy of the emitting particle must be reduced by γ γ m e c 2 .

Schwinger pair production
Schwinger pair production is the generation of electron-positron pairs from the fluctuations of the quantum vacuum in the presence of a sufficiently strong electromagnetic field. An expression for the Schwinger pair production rate per unit volume can be found in [52]: where = E/E s and η = H/E s . E and H are given by where F and G are the invariants of the electromagnetic field and are equal to Electron-positron pairs generated via the Schwinger pair production process can be initialized at rest. In principle, the electromagnetic field should lose an amount of energy equal to 2m e c 2 when pairs are created via the Schwinger process. However, since the field loses significantly more energy while accelerating these particles to relativistic velocities immediately after their creation, the small energy loss due to the rest-mass energy of the pair can be safely disregarded. Schwinger pair production is implemented in PICSAR-QED. However, since the implementation is relatively simple with respect to the other sf-QED process (the pair production rate is not very expensive to compute and product particles are initially at rest), it will not be further discussed.

Numerical implementation
The total cross-sections for Breit-Wheeler pair production and inverse Compton emission have quite complex expressions, featuring special functions and multiple integrals, which would be too computationally expensive to evaluate for each particle at each time step. Indeed, in the standard PIC algorithm, the number of operations per particle per timestep is relatively small if compared with what would be required to compute QED cross-sections. Therefore, their evaluation at runtime would largely dominate the simulation time, unacceptably slowing down the simulation. For this reason, as documented in the literature [22,23,25], the standard approach is to reformulate the total cross-sections as a product between simple numerical factors and a numerically expensive function, which is pre-computed and stored in a one-dimensional lookup table.
For instance, the right hand side of Eq. (8) is the product of a constant (α f m e c 2 / ), the normalized photon energy γ γ (which is a simple function of the photon momentum), the quantum parameter χ γ (which is a simple function of the photon momentum and of the electromagnetic field), and the function T (χ γ ), which contains all the other terms of the total cross section (see Eq. (9)). Similarly, for the inverse Compton emission total cross-section, all the numerically expensive terms can be absorbed into the G(χ) function (see Eq. (15) and Eq. (17)). The cumulative probability distributions -required to determine the properties of the product particles -are also unpractical to compute at runtime. Therefore, Eq. (11) and (18) are pre-computed over a finite set of parameters and the result is stored in two-dimensional lookup tables. As mentioned in section 2, in principle a random number per particle at each timestep should be drawn in order to determine if a sf-QED process occurs (two underlying assumptions are that the QED cross sections do not vary significantly over one timestep and that the probability of a QED process to occur during a timestep is significantly smaller than one). However, generating pseudo-random numbers can have a significant numerical cost, depending on the algorithm. Therefore, as documented in the literature [22], the preferred approach is to assign a quantity τ , called optical depth, to each particle which may undergo a sf-QED process. τ is extracted from an exponential probability distribution p(x) = exp(−x), and at each iteration it is updated as τ n+1 = τ n − dN/dt∆t, where dN/dt is the total cross-section of either Breit-Wheeler pair production or inverse Compton photon emission. This second approach is equivalent to the former, but with a reduced computational cost. Moreover, from a numerical point of view, the simpler loop on the particles to update the optical depth offers more opportunities for the compiler to optimize the code (e.g., exploiting Single Instruction on Multiple Data parallelization on CPU architectures). This section describes the specific implementation choices made for PICSAR-QED. In particular, we provide details on how lookup tables are calculated and how interpolation within these tables is performed. We also assess how the precision of the lookup tables (number of points and use of single or double precision) affects the accuracy of the results. This is particularly important from the perspective of a user. Indeed, although the general idea of the method to implement sf-QED processes in PIC codes has already been described in the literature, to the best of the authors' knowledge, detailed guidelines on how to choose the parameters of the lookup tables have never been published. At the end of the section, we finally discuss specific choices aimed at achieving portability across multiple architectures. Before delving into the implementation choices made for PICSAR-QED, it is important to clarify that methods to compute the lookup tables are provided for CPU architectures only. This is due to the fact that they require special functions not yet implemented for GPUs (e.g., Bessel functions of fractional order), and rely on the CPU-only library Boost for sophisticated quadrature methods, such as tanh-sinh [53,54]. Computing the lookup tables typically requires only few tens of seconds on a multi-core CPU, so, in principle, they could be generated at the beginning of each simulation. In practice, it is often more convenient to store them on disk and load them whenever needed (lookup tables typically require only few megabytes of storage).

Nonlinear Breit-Wheeler pair production: implementation choices and benchmarks
Two lookup tables are needed for Breit-Wheeler pair production: a one-dimensional table for T (χ γ ) and a two-dimensional table for the cumulative probability distribution P (χ γ , χ). For T (χ γ ), PICSAR-QED library adopts a solution very similar to that described for the Smilei PIC code [23,38]. The T (χ γ ) table is generated between a minimum value χ γ,min and a maximum value χ γ,max , with N χγ points logarithmically distributed between the extrema (actually, ln T is stored in the table). The choice of a logarithmic scale for χ γ allows spanning several orders of magnitude with a limited number of points, following the strategy proposed in [22]. Outside the extrema of the table, we use the approximations in Eq. (10). In practice, this is not a significant issue, provided that χ min 0.1 and χ max 1000, since the asymptotic limit in Eq. (10) is a very good approximation (besides, at χ γ ∼ 0.1 Breit-Wheeler pair production rapidly becomes negligible). On the other hand, in the range χ min < χ γ < χ max , we perform an interpolation. Specifically, in order to calculate T (χ * γ ), we first individuate two contiguous tabulated values χ γ,n and χ γ,n+1 such that χ γ,n ≤ χ * γ < χ γ,n+1 . Then, we compute T (χ * γ ) as : where T γ,n and T γ,n+1 are the tabulated values corresponding to χ γ,n and χ γ,n+1 . Figure 2 provides detailed results on how different choices of N χγ and performing all the calculations in single or double precision affects the accuracy of the table. As expected, we find that increasing the number of table points reduces the error. In order to achieve an error below few percents for χ γ > 0.1, tables must be calculated with at least 128 points. Calculating the tables and performing the interpolation in single or in double precision does not seem to affect the final error significantly. For P (χ γ , χ), PICSAR-QED adopts a significantly more complex strategy, which partially differs with respect to implementations described elsewhere. First of all, we consider P (χ γ , ξ), where ξ = χ/χ γ , so that 0 ≤ ξ ≤ 1. Moreover, we can exploit the symmetry P (χ γ , ξ) = 1 − P (χ γ , 1 − ξ) to store the table only in the range 0 ≤ ξ ≤ 0.5. P (χ γ , ξ) is then generated in the range χ γ,min ≤ χ γ ≤ χ γ,max , with N χγ points logarithmically distributed between the extrema, and in the range 0 ≤ ξ ≤ 0.5, with N ξ linearly spaced points. If χ γ < χ γ,min or χ γ > χ γ,max , we replace χ γ with either χ γ,min or χ γ,max . This means that a user must choose those extrema in such a way that the whole χ γ range relevant for a given application is included in the table. When a photon with χ γ decays into an electron-positron pair via Nonlinear Breit-Wheeler pair production, a random number r is extracted from a uniform distribution between 0 and 1, which is used to calculate the quantum parameters of the generated particles. If r ≤ 0.5, χ is the quantum parameter of the generated electron. Otherwise χ is the quantum parameter of the generated positron. In this second case, we replace r → (1 − r), so as to enforce 0 ≤ r ≤ 0.5. At this point, as in the previous case, we individuate the two contiguous tabulated values χ γ,n and χ γ,n+1 such that χ γ,n ≤ χ * γ < χ γ,n+1 . We can now define N ξ new values P m : where P n,m is the table value corresponding to χ γ,n and ξ m . By performing a binary search, we can find m * such that: We can finally calculate χ with a second linear interpolation:

Inverse Compton photon emission: implementation choices and benchmarks
As for nonlinear Breit-Wheeler, two lookup tables are needed for inverse Compton photon emission: a one-dimensional table for G(χ) and a two-dimensional table for the cumulative probability distribution P (χ, χ γ ). For G(χ), the PICSAR-QED library adopts a solution very similar to that described for the Smilei PIC code [23,38]. The G(χ) table is generated between a minimum value χ min and a maximum value χ max ,  is indeed proportional to the number of product particles generated with χ = ξχ γ . with N χ points logarithmically distributed between the extrema (actually, ln G is stored in the table). The choice of a logarithmic scale for χ allows spanning several orders of magnitude with a limited number of points, following the strategy proposed in [22]. Outside the extrema of the table we use either the first or the last value stored in the table, while within this range we perform an interpolation (which means that a user must select those extrema in order to cover all the χ range of interest). Specifically, in order to calculate G(χ * ), we first individuate the two contiguous tabulated values χ n and χ n+1 such that χ n ≤ χ * < χ n+1 . Then, we compute G(χ * ) as : where G n and G n+1 are the tabulated values corresponding to χ n and χ n+1 . Figure 4 provides detailed results on how different choices of N χ and performing all the calculations in single or double precision affects the accuracy of the table. As for Breit-Wheeler pair production, increasing the number of table points reduces the error and calculating the tables and performing the interpolation in single or in double precision doesn't seem to affect the final error significantly. However, in this case, a resolution as low as 32 points is already enough to reduce the error below the percent level across the whole χ range considered here. As far as P (χ, χ γ ) is of concern, as for Breit-Wheeler pair production, PICSAR-QED adopts a significantly more complex strategy, which partially differs with respect to implementations described elsewhere. First of all, we consider P (χ, ξ), where ξ = χ γ /χ, so that 0 ≤ ξ ≤ 1. P (χ, ξ) is then generated in the range χ min ≤ χ ≤ χ max , with N χ points logarithmically distributed between the extrema, and in the range ξ min ≤ ξ ≤ 1, with N ξ logarithmically distributed points. ξ min must be low enough that photons below the threshold contribute negligibly to the total energy loss via inverse Compton. If χ < χ min or χ > χ max , we replace χ with either χ min or χ max , so that the extrema must be selected in order to include all the χ range of interest. In the table we actually store ln P instead of P . When an electron or a positron with quantum parameter χ emits a high-energy photon via inverse Compton process, a random number r is extracted from a uniform distribution between 0 and 1. At this point, as in the previous case, we individuate the two contiguous tabulated values χ n and χ n+1 such that χ n ≤ χ * < χ n+1 . We can now define: ln P m = ln P n,m + (ln χ γ − ln χ γ,n ) ln P n+1,m − ln P n,m ln χ γ,n+1 − ln χ γ,n where P n,m is the table value corresponding to χ n and ξ m . By performing a binary search, we can find m * such that: We can finally calculate χ γ with a second linear interpolation:

Portability across different architectures
From a technical point of view, PICSAR-QED is a C++14, header-only library, designed to integrate easily into other projects, and to provide methods able to run efficiently on different computing architectures. This section describes how these goals are achieved.
Support of different unit systems -Internal calculations in PICSAR-QED are performed adopting Heaviside-Lorentz units with 1 MeV chosen as the reference energy. However, the interface of PICSAR-QED also supports SI units and normalized units where the speed of light, the elementary charge, and the electron mass are equal to one and either a reference length or a reference frequency is used (in these cases the value of the reference quantity must be provided). In practice, since the library is written in C++, this is achieved via templates, in order to avoid code duplication. The choice of the units is performed at compilation time, in order to avoid overheads at runtime. Since lookup tables are adimensional, once generated they can be used with any choice of unit system. Support for single and double precision -We provide simple and double precision versions of each method (in practice, since the library is written in C++ this is also achieved via templates). This flexibility is crucial, especially for methods that should be used at runtime, since on several architectures running in single precision leads to large speedups. Moreover, it allows integrating PICSAR-QED with codes using single precision, avoiding expensive runtime floating point casts. Lookup tables can be computed either in single or double precision, with the former case being significantly faster due to relaxed tolerance required for numerical integration. PICSAR-QED also provides an option to compute the lookup tables in double precision and store them in Avoid dependencies on specific pseudo-random number generators -PICSAR-QED is a Monte Carlo module. Therefore, pseudo-random numbers are needed at runtime. Since pseudo-random number generators have different interfaces in different libraries and performance portability frameworks, we decided not to force the use of a specific pseudo-random number generator, nor to include a pseudo-random number generator in PICSAR-QED. Our design specifications requires that a random number (uniformly distributed between zero and one) is passed to each runtime function requiring randomness. This gives complete freedom to the users on how to generate such random numbers.
Compatibility with different architectures -PICSAR-QED provides a collection of methods that can be divided into two categories: runtime methods (which are actually needed during a Particle-In-Cell simulation) and lookup table generation methods, which are needed only to generate lookup tables for later use. Only runtime functions need to be portable on different architectures, while methods for lookup-table generation need only to run on CPUs (moreover, their compilation for some architectures, namely GPUs, is not currently possible). In order to achieve portability across different architectures, all the runtime methods are pre-pended with some macros, whose values must be set appropriately to compile the code for CPUs or GPUs, or to use performance portability frameworks like Kokkos [55] and AMReX [56], as explained in detail in appendix A. While Kokkos is primarily designed as a performance portability framework, AMReX is actually a library designed to support massively parallel block-structured adaptive mesh refinement (AMR) applications, but it also offers features enabling performance portability of the applications built on top of it, like the PIC code WarpX.
Another key concept enabling portability concerns the data structures. In this regard, PICSAR-QED provides containers, such as those used internally for the lookup-tables, which must be initialized with methods running on CPUs, while at the same time being available in GPU kernels, if the library is compiled for those architectures. Achieving this may require some effort from the user, but the amount of effort is minimal for Kokkos and AMReX, as well as for for programming models like CUDA [57], as shown in appendix B.

Performance Benchmarks on different architectures
Since PICSAR-QED can be compiled for different architectures and integrated with different performance portability frameworks, we have carried out extensive performance benchmarks of the four most important kernels of the library: • Breit-Wheeler optical depth evolution • Breit-Wheeler pair production • Inverse Compton optical depth evolution

• Inverse Compton photon emission
Those benchmarks were carried out with a test program using CUDA on an NVIDIA Quadro GV100 GPU, with a test program using OpenMP on a dual-socket machine with Intel Xeon Gold 6152 CPUs and on an AMD EPYC 7302 CPU, and with a test program using the Kokkos library on all the aforementioned architectures. In some selected cases, the effect of changing the number of threads and enabling non IEEEcompliant aggressive floating-point optimizations ("fast math") was tested as well. We also performed some initial benchmarks on the Fujitsu A64FX CPU, which demonstrate that PICSAR-QED can be used on this architecture, but the results are too preliminary to be included in a fair benchmark. In all cases the benchmark was carried out with 10 8 particles, each one having ten real components: the three components of the momentum, the six components of the electromagnetic field, and the optical depth. GeV /c, 1 GeV /c] and the optical depth from an exponential distribution. 1D lookup tables were generated with 256 points, while 2D lookup tables were generated with 256 × 256 points. χ min and χ γ,min were chosen to be, respectively, 0.001 and 0.02. χ max and χ γ,max were both 10 3 , while ξ min was 10 −12 . Each kernel was tested in double precision and in single precision for all the 10 8 particles. Appendix C provides details on how the code was compiled in each case. Table 1 reports the results of these benchmarks.
The results show remarkable performances with the NVIDIA GV100 GPU, and a near perfect halving of the time spent in each kernel passing from double precision to single precision (as for many architectures, with NVIDIA GV100 the theoretical maximum FLOPS in single precision is twice the maximum FLOPS attainable in double precision).
On the other hand, we rarely observe such halving in the case of CPUs. However, this can be explained as follows. With CPUs, the performance doubling passing from double precision to single precision occurs normally only if a kernel can take full advantage of SIMD vectorization. For SIMD vectorization to occur, a kernel must respect relatively stringent requirements. Pair production and photon emission kernels are too complex to be fully vectorized, especially due to the branching conditions in the binary search step included in those kernels. Optical depth evolution kernels are significantly easier to vectorize, although, as shown in the table, this requires to relax floating point arithmetic (relaxing strict IEEE compliance for floating point operations allows most compilers to perform more extensive optimizations, such as reordering of commutative math operations and using approximate versions of some functions). We also observed that using Kokkos introduces an overhead, which is small (even negligible in some cases) for the optical depth evolution kernels. However, we observed a significant overhead in the case of the pair production kernel and of the photon emission kernel. Investigation of this issue is currently in progress. In any case, although the overhead for some kernels is substantial, in a complete PIC simulation this is not likely to be a significant   issue. Indeed, while the optical evolution kernels are executed at each timestep for each particle, pair production and photon emission kernels are relevant only for the few of them undergoing a QED process in a given timestep (for the simulation to be reliable, the timestep must be chosen small enough for pair production and photon emission to occur significantly less than every timestep for each particle [22]). A fairer comparison between different computing architectures requires to take into account their different power consumption (i.e., 250 W for NVIDIA Quadro GV100, 280 W for a pair of Intel Xeon Gold 6152, and 155 W for AMD EPYC 7302). Therefore, table 2 reports, for a selection of cases, the total energy required to apply each kernel to all the particles. We note that, without using Kokkos, performing a test on a NVIDIA Quadro GV100 GPU requires one order of magnitude less energy than performing the same test on CPU.

Integration of PICSAR-QED with WarpX PIC code
In order to integrate PICSAR-QED with WarpX, we followed standard procedures described in the literature [24]. We first added photon macro-particles. Then, for particles involved in QED processes, we added a new real component: the optical depth. We then modified the particle evolution routines in order to update this component for particles participating in QED processes. Only for particles whose optical depth reaches zero, we call either a routine to generate a pair, if the particle is a photon, or a routine to emit a photon, if the particle is an electron or a positron. In order to avoid the generation of too many low-energy photon particles, WarpX allows setting a threshold for keeping such particles after creation. If the photon energy is below such threshold, the emitting particle is slowed down as if the photon were emitted, but the photon particle is discarded immediately after creation. By default, the threshold is set to 2m e c 2 , because particles with lower energy will almost certainly not decay into a pair.
Since simulating QED effects makes sense only when QED cross-sections are sufficiently high, WarpX offers the possibility to select a minimum quantum parameter threshold below which Breit-Wheeler pair production is ignored and photon emission is simulated with a classical model [58,59]. When WarpX is compiled for GPUs, all the routines dealing with QED processes (optical depth initialization and evolution, pair production and photon emission) are carried out in GPU kernels, with the exception of the initial lookup tables generation. Lookup tables generation typically requires a few minutes, and can benefit from multi-core parallelization (which can bring the time required to compute the tables down to few tens of seconds). The user can either generate tables by setting the parameters in WarpX input file or load a lookup table generated beforehand.

Benchmarks with existing literature
In order to test WarpX+PICSAR-QED, we reproduced the three test cases mentioned in [22]. In this paper by Ridgers et al., the authors describe their implementation of a QED module to simulate QED processes in laser-plasma interaction. They test this module (which is built-in in the EPOCH PIC code) in the following three physical scenarios: (i) an electron population with initial Lorentz factor γ 0 = 1000 propagating in a perpendicular, static, magnetic field B 0 = 10 −3 E s /c. This case corresponds to an initial quantum parameter χ = 1. The duration of the simulation is 1 fs.
(ii) an electron population with initial Lorentz factor γ 0 = 4120 counter-propagating with a circularly polarized plane wave with E 0 = 1.22 × 10 −4 E s . This case corresponds to an initial quantum parameter χ = 1. The duration of the simulation is 3 fs.
(iii) an electron population with initial Lorentz factor γ 0 = 1000 propagating in a perpendicular, static, magnetic field B 0 = 9 × 10 −3 E s /c. This case corresponds Figure 6: The plots show particle energy spectra obtained with WarpX+PICSAR-QED for the 3 simulation cases described in [22]. For each case the simulation was performed multiple times, in single and double precision and using lookup tables of different sizes. Only simulations performed in single precision are shown. In each plot, the black line is the expected result, obtained by directly solving the equations for the evolution of the distribution functions of particles and photons, as described in [22,60].
to an initial quantum parameter χ = 9. The duration of the simulation is 0.1 fs.
We simulated these three cases with WarpX+PICSAR-QED, using the code in single precision and with lookup tables resolutions ranging from 32 points to 256 points in each dimension. Fig. 6 reports the energy spectra of the particles at the end of the simulation, for the three cases and for the 4 lookup table resolutions tested in the benchmark. Our results closely reproduce those published in the paper by Ridgers et al., which confirms the correctness of the QED modules in WarpX. Concerning the results reported in Fig. 6, we repeated the tests also in double precision, observing only minor differences. Running the code in single precision (and therefore using lookup tables calculated in single precision), however, has an effect on the photon energy spectrum at very low energy, as shown in fig. 7. The spectra are very similar, but in single precision the low energy part of the spectrum is truncated. This is simply due to the smaller range of real numbers that single precision can represent, and is a minor concern, since photons with Figure 7: Photon energy spectra obtained with WarpX+PICSAR-QED for case 1 described in [22]. The plot shows results obtained in single and double precision, with lookup tables having 256 points in each direction. The vertical line marks the 2m e c 2 energy threshold below which photons are deleted by default right after creation.
such a low energy have a negligible impact on the total radiative losses. Moreover, those photons are orders of magnitude below 2m e c 2 , which means that they cannot decay into electron-positron pairs via Breit-Wheeler pair production.

Conclusions
This paper presents PICSAR-QED, a Monte Carlo module to simulate strong-field QED processes in PIC codes designed for next-generation supercomputers. PICSAR-QED is conceived to be easily included in other processes and to be portable to different architectures. The outcome of several benchmarks, which demonstrate that the library produces reliable results and can run efficiently on CPUs and GPUs, is reported. Moreover, a detailed investigation on how the parameters used to generate the lookup tables affect the simulation results is provided. The main outcomes of this investigation are that (i) performing the calculations of the QED routines in single precision does not affect the accuracy of the results, and (ii) lookup tables with 128 or 256 points in each direction are sufficient to ensure negligible interpolation errors. Finally, the integration of PICSAR-QED with a state-of-the-art PIC code, WarpX, and the validation of such integration are discussed. WarpX+PICSAR-QED has been already used for production simulations of laser-plasma interaction at QED relevant intensities [61].

Appendix A: portability across different architectures
In PICSAR-QED, all the runtime methods are pre-pended with two macros, PXRMP_GPU_QUALIFIER and PXRMP_FORCE_INLINE, as exemplified here for the function calculating the Schwinger pair production rate: By setting those two macros according to table 3, runtime functions can be either: • compiled directly for CPUs.
• compiled directly for NVIDIA GPUs using CUDA.
• made compatible with the Kokkos performance portability library (and then compiled for CPUs or GPUs).
• made compatible with the AMReX library (and then compiled for CPUs or GPUs).  Table 3: Reference table on how to set PXRMP_GPU_QUALIFIER and PXRMP_FORCE_INLINE to compile PICSAR-QED methods for CPUs and GPUs and/or to make them compatible with AMReX or Kokkos. * for most compilers this must be inline __attribute__((always_inline)), which is the default value in PICSAR-QED.

Appendix B: details on the data structures used in PICSAR-QED
This appendix provides some details on the underlying data structures used for the lookup tables by PICSAR-QED. Lookup tables use several 1D vectors to store table data (even 2D lookup tables use 1D vectors as their underlying data structure). Those vectors must be manipulated on CPU when the lookup table is generated or loaded from disk, but they must be accessible within GPU kernels if PICSAR-QED is compiled for GPU architectures. Moreover, a machinery to seamlessly pass data between the CPU and the GPU is needed. AMReX already provides such a vector type: amrex::Gpu::DeviceVector (with the additional prescription of calling a GPU synchronization method after table initialization). In pure CUDA, a very thin wrapper around the device_vector provided by CUDA thrust library is enough to fulfill these requirements: PICSAR-QED calls the method pxr_sync() whenever needed, if it detects that a type having such a method is used as the underlying vector type. On CPUs, lookup tables can be passed directly as a constant reference to runtime functions. On GPUs, an additional step is needed: each table can export a "table  view", which is a lookup table of the same kind internally using non-owning pointers to  the raw table data. This "table view" can be safely passed by copy to GPU kernels. This more complex procedure also works on CPUs. Therefore, if libraries such as AMReX or Kokkos are used, the same code can be compiled for both GPUs and CPUs.
Appendix C: technical details on PICSAR-QED benchmarks Table 4 reports several details (compiler versions, compilation options . . . ) on the configurations used for the benchmarks discussed in section 4. To ensure reproducibility of our results, we also provide the git commit number corresponding to the version of the software that we have used in this work: • PICSAR-QED [42] : c16b642e3dcf860480dd1dd21cefa3874f395773 • WarpX [62] : b83f2949a1ac2eed003e991e9653b8427716bf14 • AMReX [63] : b15b1cf8d282cbb2c0d0bc0c7b049a79375ea63c The code used for the performance benchmarks will be made available on the repository of PICSAR-QED.