Optimized event generator for strong-field QED simulations within the hi-$\chi$ framework

TL;DR

This paper introduces an optimized event generator for strong-field QED simulations that significantly reduces computational costs by minimizing rate calculations, enhancing performance over previous methods within the hi-$hi$ framework.

Contribution

It presents a novel method that minimizes rate computations per QED event, boosting simulation efficiency by over tenfold within an open-source Python toolkit.

Findings

Performance improved by more than an order of magnitude.

Method reduces computational resource requirements.

Open-source implementation available for collaborative use.

Abstract

Probabilistic generation of photons and electron-positron pairs due to the processes of strong-field quantum electrodynamics (SFQED) is often the most resource-intensive part of the kinetic simulations required in order to model current and future experimental studies at high-intensity laser facilities. To reduce its computational demands one can exploit tabulation of the precomputed rates, time-step sub-cycling, dynamic down-sampling of particle/photon ensembles and other approaches. As the culmination of previous improvements, the method described here provides the opportunity to make the minimal possible number of rate computations per QED event and, therefore, this method can increase performance by more than an order of magnitude. The computational routine is publicly available as a part of the open-source framework hi-$\chi$ designed as a Python-controlled toolbox for collaborative development.

Tables

Table 1. Table 1: Details of approximation of the C 1 ​ ( χ ) subscript 𝐶 1 𝜒 C_{1}(\chi) term of the Compton rate.

$\chi$ interval	$C_1\left(\chi \right)$ trend	$f\left(\chi \right)$	Additional information about approximation of $C_1\left(\chi \right)$
	$\sqrt[3]{\frac{2}{3\chi }}$	$\frac{2}{3\chi }$	Degree $1$ polynomial
	$1$	$\sqrt[6]{\frac{2}{3\chi }}$	Degree $6$ polynomial
	$1$	$\sqrt[3]{\frac{2}{3\chi }}$	Rational polynomial $\frac{P\left(f\left(\chi \right)\right)}{Q\left(f\left(\chi \right)\right)}$, degrees $6$ and $4$, respectively
	$\sqrt[3]{{\left(\frac{2}{3\chi }\right)}^2}$	$\frac{2}{3\chi }$	Degree $6$ polynomial
	$1$	$\frac{2}{3\chi }$	Degree $6$ polynomial
$\frac{2}{3\chi }⩾100$	$1$	$\frac{3\chi }{2}$	Degree $4$ polynomial

Table 2. Table 2: Details of approximation of the C 2 ​ ( χ ) subscript 𝐶 2 𝜒 C_{2}(\chi) term of the Compton rate.

$\chi$ interval	$C_2\left(\chi \right)$ trend	$f\left(\chi \right)$	Additional information about approximation of $C_2\left(\chi \right)$
	$1$	$\sqrt[3]{\frac{2}{3\chi }}$	Degree $3$ polynomial
	$\sqrt[3]{{\left(\frac{2}{3\chi }\right)}^7}$; $\sqrt[3]{{\left(\frac{2}{3\chi }\right)}^{11}}$; $\sqrt[3]{\frac{2}{3\chi }}$; $\sqrt[3]{{\left(\frac{2}{3\chi }\right)}^5}$; $\frac{2}{3\chi }$; ${\left(\frac{2}{3\chi }\right)}^2$	${\left(\frac{2}{3\chi }\right)}^2$	$C_2\left(\chi \right)$ considered as a sum of $6$ functions with individual trends, these functions are approximated independently. For each term the degree of the polynomial is $10$.
	$1$	$\sqrt[4]{\frac{2}{3\chi }}$	Rational polynomial $\frac{P\left(f\left(\chi \right)\right)}{Q\left(f\left(\chi \right)\right)}$, degrees $10$ and $2$, respectively
$\frac{2}{3\chi }⩾100$	$1$	$\frac{3\chi }{2}$	Degree $7$ polynomial

Table 3. Table 3: Details of approximation of the B 1 ​ ( χ γ ) subscript 𝐵 1 subscript 𝜒 𝛾 B_{1}(\chi_{\gamma}) term of the Breit-Wheeler rate.

${\chi }_{\gamma }$ interval	$B_1\left({\chi }_{\gamma }\right)$ trend	$f\left({\chi }_{\gamma }\right)$	Additional information about approximation of $B_1\left({\chi }_{\gamma }\right)$
	$1$	$\sqrt[3]{\frac{8}{3{\chi }_{\gamma }}}$	Degree $4$ polynomial
	$1$	$\sqrt[3]{\frac{8}{3{\chi }_{\gamma }}}$	Degree $8$ polynomial
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[3]{\frac{8}{3{\chi }_{\gamma }}}$	Degree $7$ polynomial
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[3]{\frac{8}{3{\chi }_{\gamma }}}$	Degree $7$ polynomial
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\frac{3{\chi }_{\gamma }}{8}$	Degree $7$ polynomial
$\frac{8}{3{\chi }_{\gamma }}⩾200$	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$1$	Without polynomial

Table 4. Table 4: Details of approximation of the B 2 ​ ( χ γ ) subscript 𝐵 2 subscript 𝜒 𝛾 B_{2}(\chi_{\gamma}) term of the Breit-Wheeler rate.

${\chi }_{\gamma }$ interval	$B_2\left({\chi }_{\gamma }\right)$ trend	$f\left({\chi }_{\gamma }\right)$	Additional information about approximation of $B_2\left({\chi }_{\gamma }\right)$
	${\left(\frac{3{\chi }_{\gamma }}{2}\right)}^4$	$\sqrt[3]{\frac{3{\chi }_{\gamma }}{2}}$	Degree $5$ polynomial
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[3]{\frac{2}{3{\chi }_{\gamma }}}$; $\sqrt[6]{\frac{3{\chi }_{\gamma }}{2}}$	The coefficient $Polynom\left(f\left({\chi }_{\gamma }\right)\right)$ in 17 is computed as the sum of two polynomials of degree 7
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[3]{\frac{2}{3{\chi }_{\gamma }}}$; $\sqrt[4]{\frac{3{\chi }_{\gamma }}{2}}$
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[6]{\frac{2}{3{\chi }_{\gamma }}}$; $\sqrt[4]{\frac{3{\chi }_{\gamma }}{2}}$
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\sqrt[6]{\frac{2}{3{\chi }_{\gamma }}}$	Degree $7$ polynomial
	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\frac{3{\chi }_{\gamma }}{2}$	Degree $7$ polynomial
$\frac{2}{3\chi }⩾50$	$e^{-\frac{8}{3{\chi }_{\gamma }}}$	$\frac{3{\chi }_{\gamma }}{2}$	Degree $7$ polynomial

Equations

$$P(t^{\ast }\le t;t\le \Delta t)=P_{\mathrm{e}}(t^{\ast }\le t).$$

$$P(t^{\ast }\le t;t>\Delta t)$$

$$=1-{\mathrm{e}}^{-\lambda \Delta t}+{\mathrm{e}}^{-\lambda \Delta t}(1-{\mathrm{e}}^{-\lambda (t-\Delta t)})$$

$$=1-{\mathrm{e}}^{-\lambda t}=P_{\mathrm{e}}(t^{\ast }\le t)$$

$$\frac{\mathrm{d}{W}_{\gamma }}{\mathrm{d}\epsilon }=\frac{\sqrt{3}}{2\pi }\frac{\alpha }{{\tau }_C}\frac{\chi }{\gamma }\left[\frac{(1-\epsilon )}{\epsilon }{F}_1(z)+\epsilon F_2(z)\right],z=\frac{2}{3\chi }\frac{\epsilon }{1-\epsilon },$$

$$W_{\gamma }(\chi )=\frac{\sqrt{3}}{2\pi }\frac{\alpha }{{\tau }_C}\frac{\chi }{\gamma }{\int }_0^1\left(\frac{(1-\epsilon )}{\epsilon }{F}_1(z)+\epsilon F_2(z)\right)\mathrm{d}\epsilon$$

$${\mathrm{cdf}}_{\gamma }(\chi ,\epsilon )=\frac{W_{\gamma }(\chi ,\epsilon )}{W_{\gamma }(\chi )},$$

$$\frac{\mathrm{d}W}{\mathrm{d}{\epsilon }_e}=\frac{\sqrt{3}}{2\pi }\frac{\alpha }{{\tau }_C}\frac{m{c}^2{\chi }_{\gamma }}{\hslash \omega }\left[({\epsilon }_e-1){\epsilon }_e{F}_1(z_p)+F_2(z_p)\right],$$

$$W_p({\chi }_{\gamma })=\frac{\sqrt{3}}{2\pi }\frac{\alpha }{{\tau }_C}\frac{m{c}^2{\chi }_{\gamma }}{\hslash \omega }{\int }_0^1\left(({\epsilon }_e-1){\epsilon }_e{F}_1(z_p)+F_2(z_p)\right)\mathrm{d}{\epsilon }_e.$$

$${\mathrm{cdf}}_p({\chi }_{\gamma },{\epsilon }_e)=\frac{W_p({\chi }_{\gamma },{\epsilon }_e)}{W_p({\chi }_{\gamma })},$$

$$\sqrt{\int\big{(}f(\chi,\varepsilon)-f_{0}(\chi,\varepsilon)\big{)}^{2}\,\mathrm{d}\varepsilon\left/\int f_{0}(\chi,\varepsilon)^{2}\,\mathrm{d}\varepsilon\right.},$$

$$\begin{array}[]{rcl}E_{y}&=&A/\sqrt{2}\times\cos(kx)\sin(\omega t),\\ E_{z}&=&A/\sqrt{2}\times\cos(kx)\cos(\omega t),\\ B_{y}&=&-A/\sqrt{2}\times\sin(kx)\sin(\omega t),\\ B_{z}&=&-A/\sqrt{2}\times\sin(kx)\cos(\omega t),\\ \end{array}$$

$$W(\chi )=C_1(\chi )+C_2(\chi ),$$

$$C_i(\chi )\approx Polynom(f(\chi ))\cdot Trend(C_i(\chi ))$$

$$W_{\gamma}(\chi)=\frac{\alpha}{\tau_{C}}\frac{1}{\gamma}\left\{\!\!\begin{array}[]{ccll}\frac{5}{2\sqrt{3}}\chi&\approx&1.44\chi,&\text{if }\chi\ll 1\\[6.0pt] \frac{14\Gamma(\frac{2}{3})}{9\sqrt[3]{3}}\chi^{2/3}&\approx&1.46\chi^{2/3},&\text{if }\chi\gg 1\end{array}\right.,$$

$${\mathrm{cdf}}_{\gamma }(\chi ,\epsilon )=\frac{1}{W(\chi )}\frac{2\cdot 3^{2/3}}{\Gamma (\frac{1}{3})}{\epsilon }^{1/3}{\chi }^{2/3}.$$

$${\mathrm{cdf}}_{\gamma }^{-1}(\chi ,r_2)=D(\chi )\cdot \sqrt[3]{r_2}$$

$$W_p({\chi }_{\gamma })=B_1({\chi }_{\gamma })+B_2({\chi }_{\gamma }),$$

$$B_i({\chi }_{\gamma })\approx Polynom(f({\chi }_{\gamma }))\cdot Trend(B_i({\chi }_{\gamma }))$$

$$W_{p}(\chi_{\gamma})=\frac{\alpha}{\tau_{C}}\frac{mc^{2}}{\hbar\omega}\left\{\!\!\begin{array}[]{ccll}\frac{3\sqrt{3}}{16\sqrt{2}}\chi_{\gamma}\mathrm{e}^{-8/3\chi_{\gamma}}&\approx&0.23\chi_{\gamma}\mathrm{e}^{-8/3\chi_{\gamma}},&\text{if }\chi_{\gamma}\ll 1\\[6.0pt] \frac{20\pi^{2}}{7\sqrt[3]{3}\Gamma(\frac{1}{3})^{4}}\chi_{\gamma}^{2/3}&\approx&0.38\chi_{\gamma}^{2/3},&\text{if }\chi_{\gamma}\gg 1\end{array}\right.,$$

$${\mathrm{cdf}}_p({\chi }_{\gamma },{\epsilon }_e)=\frac{1}{W({\chi }_{\gamma })}\frac{3}{8\sqrt{\pi }}({\chi }_{\gamma }{\epsilon }_e{)}^{3/2}\exp (-\frac{2}{3{\chi }_{\gamma }{\epsilon }_e})$$

$${\mathrm{cdf}}_p^{-1}({\chi }_{\gamma },r_2)=L({\chi }_{\gamma })\cdot \frac{{\chi }_{\gamma }}{LambertW\left[M({\chi }_{\gamma })\cdot {\left(\frac{{\chi }_{\gamma }}{r_2}\right)}^{\frac{2}{3}}\right]},$$

Full text

1 Introduction

Modern high-intensity lasers provide the opportunity to reach electromagnetic field intensities sufficient for invoking processes of strong-field quantum electrodynamics (SFQED), namely non-linear Compton scattering and Breit-Wheeler pair production [1]. These processes trigger a variety of phenomena, for both single-particle and collective dynamics of emerging plasmas, with potential applications as novel particle and radiation sources, as well as for fundamental studies relevant to SFQED, plasma physics and astrophysics [2].

Numerical simulations with probabilistic account for SFQED processes play an important role in both theoretical studies and experiment design for the upcoming experimental programs in this field. A particularly fruitful methodology is based on extending the Particle-In-Cell (PIC) method with a module that accounts for the SFQED processes via probabilistic generation of macro-photons, electrons and positrons that sample the 3D3P distribution by an ensemble of reasonable size [3, 4, 5, 6, 7, 8, 9, 10]. The implementations of this approach, also referred to as QED-PIC codes, include EPOCH [6], PIConGPU [11], VPIC 2.0 [12], OSIRIS [13], Smilei [9], Tristan-MP [14], Calder [15], PICADOR [16], PICSAR-QED [10] to name a few.

Although the QED-PIC method is practical, it can become computationally demanding due to two main reasons: (1) the prolific generation of new particles/photons and (2) the extremely short time step sometimes required in order to resolve the characteristic temporal interval between SFQED events (as compared to the laser wave period and other macroscopic scales of the simulated processes). The former is commonly addressed by resampling the ensemble through either merging/coalescing macro-particles/ photons [17, 18, 19, 20] or by removing some of the macro-particles/photons and redistributing their statistical weight among the remaining ones [21, 3, 22, 23].

The way to handle computational demands caused by the short characteristic time between SFQED events depends on the method used for sampling these events. The use of rejection sampling implies that for each time step and for each photon/particle two random numbers are generated: one that defines the proposed energy of photon/particle to be generated and another that determines whether the proposal is accepted or rejected [3]. This method is particularly simple to implement as it requires the computation of two special functions, more precisely the first and second synchrotron functions [7], found in some standard libraries. The low-energy singularity in the probability distribution for Compton scattering can be resolved via increased sampling of low energy proposals through the method of importance sampling [24]. By setting the requirements on the time step and implementation of time step sub-cycling to meet said requirements, it also becomes possible to meet the general rejection sampling requirements for any simulation setting [7]. For each particle/photon this procedure permits the generation of arbitrarily many events (even local cascades) during each global time step, making the whole computational method versatile for studies of extreme processes [25, 26, 27].

However, the use of rejection sampling has a disadvantage: even if the restrictions on the sub-step duration are mildened, the method requires multiple, and sometimes expensive, computations of the involved special functions per accepted SFQED event. [28] The reason for this is that events need to be rare (per time step) in order to be accurately counted. In the limit of the time step being large there is a unit probability of accepting proposals, resulting in one event per time step. This leads to an artificially uniform generation of events and the exclusion of cases when more than a single event happens within a single time sub-step.

In order to avoid unnecessary computations of the special functions one can use inverse sampling, in which a distribution is sampled through its cumulative distribution, which has already been used in the area of PIC simulations to account for ionization processes (see e.g. [29]). For SFQED processes the use of inverse sampling implies the following procedure [30, 6, 31]. For each particle/photon we generate a uniformly distributed in the unit interval random number $p$ and compute the so-called optical depth $d$ by inverting $p=1–e^{-d}$. If we represent the probability of an event not happening within the time interval $(t,t+t_{1})$ as $P(t,t+t_{1})=1-\exp\left(-\int_{t}^{t+t_{1}}R(\tau)d\tau\right)$, where $R(\tau)$ is the instantaneous event rate over the particle’s trajectory, then the exponent plays the role of a cumulative optical depth. As the particle propagates, the integral in the exponent is computed using first-order Eulerian integration until it reaches the value of $d$, at which point the event is assigned to occur. The energy of the outgoing particles is determined randomly using inverse sampling through pre-computed look-up tables that store the inverse cumulative probability as a two-dimensional function of energy and the quantum nonlinearity parameter $\chi$. A new value of optical depth is then generated in order to determine the moment for the next event to happen, and so on.

In the inverse sampling procedure that was just described one of the most computationally demanding parts is the evaluation of the rate, which is performed once per particle per global time step. This is done in many state-of-the-art QED-PIC codes, including [6, 10]. For many cases where the global time step is set by the requirements of the field solver, this is not a very limiting restriction, as it becomes necessary to reevaluate the rate over this time scale anyways in order to account for the field evolution. However, as was the case for the rejection sampling scheme described earlier, the global time step has to be sufficiently small in order to ensure that events are rare. This means that several evaluations of the rate are necessary for the same reason as described before. Furthermore, choosing the global time step according to the highest rate expected to be observed in a simulation can result in additional computational demands.

While there are many techniques for reducing computational demands, it is difficult to assess performance in an exact and general manner, as performance strongly depends on factors such as problem, method, hardware and various implementation details. We have here described a number of methods that are commonly employed to improve the performance of SFQED computations. However, depending on factors such as those just listed, the performance for any given scenario might naturally vary. Therefore, the development and implementation of a unified methodology for comparing the performance of PIC codes is a very difficult problem. However, we can focus on minimizing the time step sub-cycling. Regardless of the hardware and software environment and implementation features of the various components of the QED-PIC method, reducing the time-step sub-cycling should speed up simulations. That is what this paper focuses on.

In this article we present the implementation of a method that combines the inverse sampling methodology described above with sub-cycling. In doing so, within each global time step we allow multiple SFQED events originated from a single particle/photon, including from secondary particles. After the optical depth is reached, we generate new particles, assign optical depths for each of them and continue the process until the end of the global time step. At this point we re-evaluate the rate for all particles/photons in the ensemble, as the global time step often defines the time scale over which the field can no longer be assumed to remain unchanged. For this reason this procedure implies the minimum number of possible rate evaluations: either one per event or one per global time step, depending on whether the optical depth is shorter or longer than the global time step. Furthermore, our implementation, described in more detail in Section 3, also employs high-precision approximations of the cumulative distribution functions with analytical treatment of limiting cases. The computational accuracy is presented in Section 4 and we perform comprehensive testing in Section 5.

2 Motivation

The QED-PIC event generator can be implemented in several different ways, but one thing in common between all of them is that they sample the events in accordance with the partial rates as described later in equations 3 and 6. The two most common ways of sampling the spectra is to employ either rejection sampling or inverse transform sampling. In the former case, an energy value is first picked at random after which rejection sampling is performed on the partial rate $(\mathrm{d}W/\mathrm{d}E)$ to determine if the event should occur. In the latter case, it is first decided if the event occurs at all using the event rate ($W=\int(\mathrm{d}W/\mathrm{d}E)\,\mathrm{d}E$), then the energy is determined through inverse transform sampling of the partial rate.

There are advantages and disadvantages to each of the two approaches. For example, the latter requires tabulation of the inverse cumulative distribution function, but in both cases the probability of a QED event to occur within one time step must, naturally, be less than unity. More rigorously, the probability of a single event within one time step must in fact be much smaller than unity in order for the number of events over an extended duration of time to be accurate. The reason for this is that we traditionally neglect the occurrence of multiple events within a single time step. The total number of events over a large number of time steps should therefore, on average, be underestimated (not accounting for any quantum interference effects). This discrepancy is minimized by making the events rare. Furthermore, a high event probability will also lead to an apparent, but unphysical, temporal correlation between such events.

This is typically only an issue in regimes of very high $\chi$ and as the simplest way to ensure sufficient accuracy is decreasing the global time step, it can for the most part be ignored. At the same time this issue can not be completely disregarded, as decreasing the global time step can be highly inefficient. Ultimately, an increased temporal resolution is only necessary for particles experiencing high values of $\chi$, and typically only for a limited period of time as the strong fields giving rise to the high $\chi$ are often well localized in both time and space. The general solution is therefore to implement what is called sub-cycling, wherein the time step is reduced only for the particles for which it is required, as determined by their instantaneous value of $\chi$. This can be performed entirely in the particle pusher, without any changes to the remainder of the PIC loop, and is commonly done by splitting the time step into an integer number of sub-steps for the selected particles.

Although only a fraction of the ensemble is typically contributing to the QED processes, the requirements on the time step will affect the efficiency of the entire event generator, as its accuracy relies on a large rejection-acceptance ratio. This in turn implies that the number of rate computations per accepted event is large. An improvement can be obtained by computing when the next emission will occur, instead of if it will occur within the current time step. This minimizes the number of expensive computations per event, even allowing it to approach unity without loss of accuracy under certain conditions. Furthermore, this way a complete separation of time scales can be obtained between the QED-event generator and the classical PIC loop, allowing the global time step to be set through consideration of only the usual PIC constraints, such as the stability of the field solver.

3 Method

3.1 Memorylessness

If we assume that events are discrete and independent then the time between events can be described through the exponential distribution. Because the exponential distribution is memoryless, this can be exploited to create an event generator that is independent of the chosen global time step of the PIC simulation. This is achieved by computing the time until the next event. If the event is determined to occur within the current time step, it can be performed at the determined time through adaptive sub-cycling. If not, the event can be discarded and the simulation pushed to the next time step.

We can easily show that this produces consistent sampling of the distribution function, independently of the choice of the time step $\Delta t$, using the cumulative distribution function (CDF) of the exponential distribution, $P_{\mathrm{e}}(t^{\ast}\leq t)=1-\exp(-\lambda t)$, where $\lambda$ is the event rate and $t^{\ast}$ is the time until the next event. We further assume that $\lambda$ remains unchanged over one time step. For $t\leq\Delta t$, the probability for an event to occur within a current time step is trivial, as it is entirely determined by the exponential distribution.

[TABLE]

With $t>\Delta t$ however, the CDF instead contains two terms: (1) The probability for an event to occur within a concurrent time step ($P_{\mathrm{e}}(t^{\ast}\leq\Delta t)$) plus (2) The probability that the event was rejected in a concurrent time step and accepted in the succeeding time step.

[TABLE]

We thus have that $P(t^{\ast}\leq t)=P_{\mathrm{e}}(t^{\ast}\leq t)$ regardless of the value of $\Delta t$. As such, the choice of time step may, with regard to the event generator, simply be seen as arbitrary synchronization points, allowing us to completely decouple any time step requirement due to particle events from the constraints due to the fields.

This important property leads to two significant improvements in comparison with the general approach [6, 30, 31]. Firstly, we do not need to impose the smallness on the time step to integrate the optical depth. Instead, at each time step we sample events for each particle, and if the event does not occur on the current step, we simply move the particle for the whole time step and discard the event. Otherwise, we move the particle during the fraction of the time step corresponding to the optical depth, to the position where secondary particles are produced. This procedure ensures temporal locality, i.e. all event processing is performed within a single time step and no state (e.g. an optical depth) needs to be saved between time steps. Secondly, after the event we can sample for new events from both initial and secondary particles and move them until all particles move to the end of the time step. A graphical sketch of how this sub-cycling works is shown in Figure 1. Thus, we can account for multiple events during a single time step, which also relax time constraints on a time step in the case of a prolific secondary particles production.

3.2 Determining the next event

For both Compton scattering and Breit-wheeler pair production the event rate is simply determined by the total rate $\lambda=W(\chi)=\int(\mathrm{d}W/\mathrm{d}E)\,\mathrm{d}E$. Knowing that the time between emissions $t^{\ast}$ is governed by the exponential distribution it can, through inverse sampling, be obtained as $t^{\ast}=-\log(r)/W(\chi)$, where $r\sim U(0,1)$ is a uniformly distributed random number. This process is equivalent to determining the optical depth of a particle.

3.3 Photon emission

The partial rate of nonlinear Compton scattering is given by [32, 33, 34, 35]:

[TABLE]

where $\varepsilon=\hbar\omega/\gamma mc^{2}$ is the photon energy normalized to the energy of the radiating particle, $\hbar$ is the Planck constant, $\omega$ is the photon frequency, $\gamma$ is the Lorentz factor, $m$ is the electron mass, $c$ is the speed of light, $\alpha$ is the fine structure constant and $\tau_{C}$ is the Compton time. $F_{1}(x)=x\int_{x}^{\infty}K_{5/3}(x^{\prime})\mathrm{d}x^{\prime}$ and $F_{2}(x)=xK_{2/3}(x)$ are the first and second synchrotron function, respectively.

3.3.1 Event rate

Consider the total rate of photon emission:

[TABLE]

As described above, the main idea of the new method is estimating the time $t^{\ast}$ after which the photon is emitted. In order to achieve this, the Compton rate $W_{\gamma}(\chi)$ is approximated by a piecewise function over $\chi$, the local trend is identified and refined using a polynomial approximation (intervals and local trends are given in A). The tables specifying intervals over $\chi$ and the local trends are provided in the supplementary materials section for a more compact narration. The employment of polynomials up to the $9^{th}$ order allows us to achieve a relative precision of $10^{-3}$, which will be demonstrated below. Polynomial values can be computed using the Horner scheme [36], based on Fused Multiply-Add (FMA) instructions, computing multiplication and addition ($a\cdot b+c$) in one operation with higher accuracy and performance. As a result, these computations are relatively lightweight for modern CPUs. We note that the proposed approach allows for even higher precision, but that would require additional computations and would be significantly slower. We believe that the proposed parameters of approximation provide an optimal compromise between accuracy and computational efficiency, while accuracy remains at least the same as in the prior scheme [7].

3.3.2 Determining the energy (Compton)

The energy can be obtained through inverse transform sampling of the partial rate, requiring the computation of the cumulative distribution function

[TABLE]

where $W_{\gamma}(\chi,\varepsilon)=\int_{0}^{\varepsilon}\frac{\mathrm{d}W}{\mathrm{d}\varepsilon^{\prime}}(\chi,\varepsilon^{\prime})\,\mathrm{d}\varepsilon^{\prime}$.

In order to obtain the required value of $\varepsilon$ the inverse function $\operatorname{cdf}_{\gamma}^{-1}(\chi,\varepsilon)$ (Compton inverse CDF) must be evaluated for $r\sim U(0,1)$, for which the method of tabulation of an inverse function was used. The domain of the inverse function was represented by a mesh as follows. The mesh is uniform over $r_{2}\in\left[0.05,0.93\right]$ (96 points), inverted logarithmic condensing to $1$ over $r_{2}\in\left[0.93,1\right]$ (96 points), and logarithmic over $\chi\in\left[1.5^{-30},1.5^{30}\right]$ (61 points). The inverse function was evaluated in the nodes. In further computations linear interpolation is used between nodes in the main region and asymptotic approximations are used over $r_{2}\in\left[0.0,0.05\right]$. Within this approach an acceptable relative deviation was considered to be below $10^{-2}$ in the specified domain. The typical value was $10^{-3}$. We experimentally found that such a complex distribution of points allows us to approximate the function accurately enough to avoid the problems described in [37]. Implementation details are given in B.

3.3.3 The particle push algorithm with account for QED effects

The constructed approximations allow us to formulate a new particle push algorithm of charged particles during a single global time step within the main PIC cycle with account for QED effects (algorithm 1). The algorithm functions as follows. Until the particle exits the loop (lines 3-14), signifying the end of the current time step, the time of emission of the next photon is continuously updated (lines 4-6). If a photon must be emitted within this time step, the particle is pushed to coordinates corresponding to the time of emission (line 11), after which a photon is emitted (line 14) with the according energy (line 13). Otherwise, the particle is pushed to coordinates corresponding to the end of the time step (line 8). The emitted photons are pushed using an algorithm described in Section 3.4.3. Note that $\chi$, $rate$, and $t^{\ast}$ are updated at the beginning of each time sub-step, while the fields are assumed to be constant throughout the global time step in accordance with the usual assumptions of the PIC method.

3.4 Photon decay into an electron-positron pair

3.4.1 Event rate (Breit-Wheeler)

The partial rate of Breit-Wheeler pair production is given by [33]:

[TABLE]

where $z_{p}=\frac{2}{3\chi_{\gamma}(1-\varepsilon_{e})\varepsilon_{e}}$. This expression shows the probability of decay of a photon into an electron-positron pair with energies $\varepsilon_{e}\hbar\omega$ and $(1-\varepsilon_{e})\hbar\omega$, respectively.

Similar to the way photon emission is handled, the time of photon decay is determined from the total pair production rate:

[TABLE]

By computing the Breit-Wheeler rate $W_{p}(\chi_{\gamma})$ the value of $t^{\ast}$, at which point in time the event occurs, can be determined in accordance with Sec.3.2. For fast computation a piecewise function with identified local trends and refinement using polynomials is used (intervals and local trends are given in C). The demanded precision is the same as in the case of photon emission.

3.4.2 Determining the energy (Breit-Wheeler)

In order to obtain the dimensionless energy of the generated electron $\varepsilon_{e}$ (the energy of the positron is then $1-\varepsilon_{e}$) the same approach as in Sec.3.3.2 is used.

[TABLE]

where $W_{p}(\chi_{\gamma},\varepsilon_{e})=\int_{0}^{\varepsilon_{e}}\frac{\mathrm{d}W}{\mathrm{d}\varepsilon_{e}^{\prime}}(\chi_{\gamma},\varepsilon_{e}^{\prime})\,\mathrm{d}\varepsilon_{e}^{\prime}$.

In order to evaluate the inverse function $\operatorname{cdf}_{p}^{-1}(\chi_{\gamma},\varepsilon_{e})$ (Breit-Wheeler inverse CDF) a mesh logarithmic over $r_{2}\in\left[0.0,0.05\right]$ (64 points), uniform over $r\in\left[0.05,0.5\right]$ (64 points), and logarithmic over $\chi\in\left[1.5^{-30},1.5^{30}\right]$ (61 points) is used. Since the function is symmetric with respect to $r_{2}=0.5$, the values $r\in\left[0.5,1\right]$ are not used and in this case $\operatorname{cdf}_{p}^{-1}(\chi_{\gamma},r)=1-\operatorname{cdf}_{p}^{-1}(\chi_{\gamma},1-r)$. Implementation details are given in D.

3.4.3 Particle push algorithm with account for QED effects (for photons)

The algorithm of particle push for photons with account for QED effects by the end of a global time step (algorithm 2) is similar to the algorithm for charged particles. Unlike the algorithm 1, we do not need a while loop, since the photon is removed during pair production and no further work is required on it. If a photon must decay within this time step, the photon is pushed to coordinates corresponding to the time of decay (line 10), after which an electron-positron pair is generated (line 13) with the according energies (line 12) and the photon is deleted (line 14). Otherwise, the photon is pushed to coordinates corresponding to the end of the time step (line 7). The generated particles are pushed using the algorithm described above in Sec.3.3.3.

4 Validation

We here present the accuracy of the implemented functions by comparing the functions to the true values, which are computed to a high accuracy (of the order of $10^{-15}$) through a combination of piecewise integration and special treatment of integrand divergences.

4.1 Event rates

The implemented Compton rate, $W_{\gamma}(\chi)$, defined in equation 4, and its relative accuracy are presented in Figure 2 and compared with the lowest order asymptotic approximations. The implemented function shows a relative error of less than $2\times 10^{-6}$ across all values of $\chi$. Similarly, the implemented Breit-Wheeler rate, $W_{p}(\chi_{\gamma})$, defined in equation 7, is also presented in Figure 2. The implemented Breit-Wheeler function shows a relative error of less than $10^{-6}$ for all $\chi_{\gamma}\gtrsim 10^{-2}$. For $\chi_{\gamma}\lesssim 10^{-2}$ the relative error increases to around $10^{-3}$, but as the Breit-Wheeler rate (as well as the absolute error) in this domain is less than $10^{-70}\,\text{s}^{-1}$, any such errors are computationally irrelevant.

4.2 Spectra

The particle spectra are defined through the derivative of the CDF with respect to energy. Because of the involvement of the derivative, the numerical discretization and the two-dimensional nature of the CDF, the accuracy of the produced spectra becomes more difficult to evaluate compared to the accuracy of the rate functions. Care must therefore be taken when interpreting the errors of the inverse CDF, and how it relates to the produced spectra.

With these caveats in mind, the inverse CDF and its resulting number spectra (given by the partial rate) are presented in Figure 3 for both the Compton and the Breit-Wheeler processes. Both the CDFs and the spectra show good agreement over a wide range of values of both energy ($\varepsilon$, $\varepsilon_{e}$) and the quantum nonlinearity parameter ($\chi$, $\chi_{\gamma}$). Nevertheless, because the CDF becomes non-invertible as $\varepsilon\rightarrow 1$, the relative error contains a divergence in this limit. However, as was the case for the Breit-Wheeler rate in certain regimes, this generally occurs where the absolute errors of the spectra are unquestionably negligible.

The normalized root mean squared (RMS) errors of both inverse CDFs and spectra are presented in Figure 4 as functions of $\chi$ ($\chi_{\gamma}$). This is computed according to

[TABLE]

where $f$ is the implemented function, and $f_{0}$ the true function. In general, the inverse CDF and its resulting spectra show typical errors on the order of $10^{-2}$-$10^{-3}$ depending on the regime and process.

5 Numerical verification

Here we conclude by testing our implementation on a few benchmark cases and comparing against previously published results.

5.1 Cascade development

In this benchmark we simulate the development of a cascade seeded by a single electron with an initial energy of $\varepsilon_{0}=$100\text{,}\mathrm{GeV}$$, moving in a strong constant magnetic field ($H_{0}=0.2E_{\mathrm{S}}$) oriented perpendicularly to the electron’s direction of motion. In Figure 5 we show the number of electrons and positrons with energy greater than $10^{-3}\varepsilon_{0}$ as a function of time. The results are averaged over $1000$ simulations and are performed for two different time steps, $10^{-3}t_{\mathrm{rad}}$ and $0.5t_{\mathrm{rad}}$, where $t_{\mathrm{rad}}=3.85\times\gamma_{0}^{1/3}(E_{S}/H_{0})^{2/3}\hbar^{2}/(mce^{2})$. The benchmark shows excellent agreement between the two choices of time step, as well as with previously published results [38, 4, 7]. Some deviations can be seen for $t/t_{\mathrm{rad}}>6$, however all discrepancies are shown to lie within 1-2 standard deviations.

5.2 Cascade spectra

Next, we again simulate a cascade seeded by a single electron in a strong constant magnetic field oriented perpendicularly to the electron’s direction of motion, but here with a field strength of $H_{0}=10^{-3}E_{\mathrm{S}}$ and initial electron energy $\gamma_{0}=1000$. The distribution of electrons $\Phi_{-}(\gamma,t)$, positrons $\Phi_{+}(\gamma,t)$ and photons $\Phi_{\gamma}(\gamma,t)$ are presented in Figure 6 as functions of energy and after $1\,\text{fs}$ of cascade development. Our results are averaged over $10^{8}$ simulation runs and compared with the results of [6, 7], showing excellent agreement.

5.3 QED cascade in circularly polarized standing wave

In this section we demonstrate capabilities of the developed method of accounting for QED processes. We consider a model example, which is nevertheless very instructive and widely used for physical insight: a QED cascade in the fields of a standing circularly polarized plane wave [39, 40, 41, 42, 43, 44, 45]. This wave is homogeneous in directions along both electric and magnetic fields, these fields rotate in time and their amplitudes are locally constant. QED processes are more frequent in the vicinity of the electric field antinode where electrons and positrons quickly gain energy after photon emission and photons have a high probability to decay. In this field configuration, when the QED cascade becomes steady-state, the number of electrons, positrons and photons grows exponentially in time, moreover, radiation losses suppress particle escape from the antinode region [44]. Thus, the considered configuration of fields can lead to likely the most frequent (in comparison with other field structures) calls of functions handling QED processes, and the calculations of these functions may take a significant part of computational time.

We performed series of simulations of a QED cascade in the standing circularly-polarized plane wave with different time steps and wave amplitudes. We analyzed the results and compared the computational performance using different event generators: AEG (the adaptive event generator) [7], mAEG (the modified adaptive event generator) [28] and FQED (the fast QED-event generator considered in this paper). For simulations with these event generators we used the same configurations. The simulation box was $64\times 2\times 2$ cells and $2\lambda\times 2\lambda\times 2\lambda$ along $x$, $y$ and $z$, where $\lambda=0.8$ $\mu$m is the laser wavelength. The time step was varied in range from $T/6000$ to $T/50$, where $T\approx 2.67$ fs is the laser period. QED cascade development was considered during $10T$. That period of time is sufficient in order to determine the cascade growth rate.

Fields were set analytically:

[TABLE]

other field components were zero. Amplitude $A$ was considered in the range from $700R_{f}$ to $2100R_{f}$, where $R_{f}=2\pi mc/eT=1.3\times 10^{8}$ G or statV/cm. For each amplitude and each time step we performed six simulations: two simulations for each QED-event generator with different seeds for the random number generator (RNG). Results of simulations, as well as their performance at certain physical parameters, were averaged over these two simulations.

Initially electrons and positrons were uniformly distributed in the whole simulation box and the quantity of each was $5\times 10^{4}$. One macroparticle (aggregation of a number of physical particles in PIC-codes) was equal to one real particle. The boundary conditions for particles were periodic due to the plane structure of the continuous wave. In order to avoid memory shortage as a result of QED cascade development we used a global leveling resampling technique [22], which does not allow the number of macroparticles of any individual type to substantially exceed a manually specified threshold. The chosen resampling method produces minimal numerical artifacts in the particle ensemble[22]. The resampling threshold was set such that resampling was run when the number of macroparticles of any given type exceeded $10^{6}$.

For tests of the developed event generator we used 1 two-CPU node of the Joint Supercomputer Center of RAS with Intel Xeon Gold 6248R. We run tests using 1 MPI process and 96 threads.

First, we analyze physical consistency of simulations performed with different event generators. We consider the QED cascade growth rate $\Gamma$ as the main result to compare. In order to determine $\Gamma$ we calculate the total quantity of electrons in the simulation box as a function of time $N_{e}(t)$ and then $\Gamma=\ln{\left[N_{e}(10T)/N_{e}(8T)\right]}/\left(2T\right)$. As reference values of the growth rate we consider values $\Gamma_{\mathrm{AEG}}$ obtained with the AEG event generator previously used in many studies (for example, see [22, 25, 26, 27, 44]).

Figure 7 (a) shows that $\eta_{\Gamma}=\left(\Gamma_{\mathrm{mAEG,FQED}}-\Gamma_{\mathrm{AEG}}\right)/\Gamma_{\mathrm{AEG}}$ (relative discrepancy of growth rates) is about 1%. Within the amplitude range $A/R_{f}>1000$, where $\Gamma T>1$, the relative discrepancy is minimal ($\lesssim 1\%$). However, within the lower amplitude range $A/R_{f}<1000$ for both the mAEG and the FQED event generators $\eta_{\Gamma}$ is several percent. This increase of $\left|\eta_{\Gamma}\right|$ can be explained by low values of $\Gamma T<1$ and consequently a larger dispersion of results which may demand averaging over more than two performed simulations with different RNG seeds.

To determine the dependence of $\Gamma$ on time step we choose $A=1414R_{f}$ at which $\Gamma T>1$ and the dispersion of results is assumed insignificant. For all three event generators fluctuations of the growth rate were within 1% (average growth rate $\approx 3.215T$), see Fig. 7 (b). The use of the FQED event generator results in minimal fluctuations within 0.6%, while AEG leads to the maximal fluctuations around 1%. The maximal considered time step is $dt=T/50$, since, according to our simulations, for longer time steps the time resolution of particle motion becomes insufficient for the given field amplitude $A=1414R_{f}$.

Besides growth rate, which is calculated based on the amount of particles in the simulation box, we also show [see Fig. 7 (c)] that spatial electron distributions obtained in simulations with different event generators are very close to each other. These distributions are local characteristics and they also confirm that the developed QED-event generators agree nearly perfectly.

Second, we compare the achievable performance of different event generators. The main part of the computations for the considered setup is related to handling of particles as the total duration of each simulation is only slightly longer than the total time spent in the particle loop stage (Fig. 8 (a), (c)). Moreover, a substantial part of computations within the particle loop are devoted to QED processes, especially for $A>1414R_{f}$. An increase of amplitude makes QED processes more probable leading to production of more pairs and photons, substantially increasing both the time spent in the particle loop as well as the overall computational time. However, owing to resampling, which at $A\gtrsim 1000R_{f}$ occurs not only for photons but also for pairs, $t(A)$ is a concave function (see Fig. 8 (a)).

Within the considered amplitude range the FQED event generator is the fastest. Its speedup in relation to AEG increases with amplitude and can exceed the factor of 20 (Fig. 8 (b)). Another event generator, mAEG, demonstrates qualitatively similar behavior, but speedup is around 2-2.5.

For a given wave amplitude $A=1414R_{f}$ we also investigate performance at different time steps [Fig. 8 (b), (d)]. At time steps $dt\ll T/3000$ all three event generators show very close performance and the time of computations scales as $\sim 1/dt$. AEG and mAEG do not need sub-cycling and mainly lines 7-9 of the algorithm 1, and lines 6-8 of the algorithm 2 of FQED are enabled.

With increasing time step the AEG event generator is the first to have $t$ deviate from the $1/dt$ dependence. At $dt>dt_{AEG}\approx T/3000$ it enables sub-cycling more frequently and as a result at $dt>T/150$ the performance of simulations with this event generator becomes nearly independent on time step. The QED-event generator mAEG determines sub-steps more precisely [28]. So, for this event generator the computational time $t$ scales as $1.2\left(dt/T\right)^{-1.08}$ up to $dt_{mAEG}\approx T/1000$ and at $dt\gtrsim T/100$ the slope of the function $t(dt)$ substantially decreases.

The event generator FQED shows the best performance among the considered event generators. The threshold time step at which $t(dt)$ changes its dependence is $dt_{FQED}\approx T/500$, and deviation from this dependence with increasing $dt$ is quite slow. The maximal speedup is achieved at the maximal possible dt: performance of simulations with FQED can be more than an order of magnitude better than that with AEG and up to an order of magnitude better than that with mAEG.

Owing to properties of particle motion [44], the characteristic time of QED processes of the considered problem is roughly $t_{\mathrm{QED}}\sim 1.5T\left(R_{f}/A\right)^{0.645}\approx T/70$. Having compared $dt_{AEG,mAEG,FQED}$ with $t_{QED}$ we conclude that AEG uses excessively short time steps during sub-cycling ($t_{QED}/dt_{AEG}>40$). The more precise estimate of sub-steps for sub-cycling implemented in mAEG ($t_{QED}/dt_{mAEG}\approx 14$) significantly improves performance, but the FQED event generator based on both the optical depth methodology and sub-cycling yields the best performance. The necessity of sub-cycling for this event generator (lines 11-14 of the algorithm 1, and lines 10-14 of the algorithm 2) is evident at time steps closest to $t_{QED}$ ($t_{QED}/dt_{FQED}\approx 7$), and the optical depth methodology allows creating particles or photons with a minimal number of function calls required in order to determine QED rates and generate random numbers.

To sum up, the developed FQED event generator of QED processes presents advantages compared to both mAEG and AEG: it has approximately the same accuracy but speeds up performance by an order of magnitude or more. The greatest benefit of FQED is observed when the QED rates are greatest (field amplitudes $\gtrsim 1000R_{f}$) and the time step is quite long. However, it is important to note that there may be limiting factors (computational or physical) that prohibit the use of long time steps at the considered field amplitudes and thus hinder computational performance. For example, particle motion or plasma oscillations may demand (within the employed numerical schemes) very detailed time resolution, with which the speedup provided by FQED is not so impressive, even for larger field amplitudes.

5.4 QED cascade in two colliding tightly focused laser beams

Although the developed event generator FQED demonstrated exceptional performance in the previous example, it is important to note that the benefit of using this method is problem-dependent. The proportion of the particle ensemble that takes part in QED processes is crucially important, as well as how strong the electromagnetic fields are, and what time resolution is sufficient. These aspects can be interconnected. For example, if we aim to enhance QED processes, laser intensity could be increased by using tighter laser focusing. This narrows the strongest-field region, making it more inhomogeneous. As a result, particles may escape the strong field region more quickly making a smaller portion of the particle ensemble continuously engaged in QED processes. Additionally, stronger field amplitudes may require shorter time steps. As a result, the interplay of different factors can diminish some of the benefits of FQED.

In this section we consider one of the basic schemes of expected experiments devoted to QED cascades using counter-propagating multipetawatt laser pulses (see, for example, [3, 46, 47, 48, 49]). Thanks to remarkable progress in laser science and technology such laser pulses are now almost in operation [50]. In contrast to the example in the previous section here we consider a more realistic numerical setup. Two identical counter-propagating multipetawatt laser beams with focusing angles equal to $45^{\circ}$ irradiate a seed hydrogen-like target in the form of a sphere with radius $1\lambda$ located in the mutual focus of these beams. Electron density of the target for numerical tests is assumed to be $n_{e}=5n_{cr}$ where $n_{cr}=1.38\times 10^{21}$ cm*-3*, the laser beam power is 25 PW (the maximal field amplitude in mutual focus is $2600R_{f}$), the laser pulse duration is 30 fs and the pulse has a $sin^{2}$ envelope.

The simulation box was $264\times 296\times 296$ cells and $8.2\lambda\times 9.1\lambda\times 9.1\lambda$ along $x$, $y$ and $z$, where $\lambda=0.9$ $\mu$m. The time step was $T/64$, where $T\approx 3$ fs. The time of simulations was $41T$. The larger simulation box and the longer time of simulation are necessary for the modeling of tightly focused laser pulses using the Total field – Scattered field technique [51] by setting currents at the closed interface near the PML (perfectly matched layer) [52]. The PML ensures absorbing boundary condition for fields; particles and photons were also absorbed at the boundaries of simulation box. E-mode laser beams propagate along the $x$ axis and the largest component of the electric field is along the $z$ axis. The center of the seed sphere was located in the coordinate origin. Initially this sphere contains $10^{7}$ macroelectrons and hydrogen-like macroions, this quantity of macroparticles is equivalent to a statistical weight of particles of 1100. The threshold for resampling of particles and photons was $10^{7}$ macroparticles. We used the minimal sufficient spatial and time resolution for clearer demonstration of the capabilities of different event generators. Such a rough resolution can be used in practice during the optimization of various laser pulse or seed target parameters, which may require quite many simulations.

For these tests we used 2 two-CPU nodes of the Joint Supercomputer Center of RAS with Intel Xeon Gold 6248R. We run tests using 2 MPI processes with 96 threads per process.

In order to demonstrate agreement between simulations performed using different QED-event generators we show the time evolution of the quantity of electrons $N_{e}$ in the vicinity of the plane $y=0$ (-$0.03\lambda<y<0.03\lambda$) (Fig. 9 (a)). At the beginning of interaction the initial target is deformed and destroyed by the leading edges of laser pulses. As a result, $N_{e}$ decreases between iterations $500<n_{it}<1000$. At this stage laser fields become relativistic but insufficient for fast QED processes. However, excessive sub-cycling within AEG increases the time required to compute the particle loop stage in comparison with the corresponding time in the cases of mAEG and FQED. For analysis this time was summed over each consecutive 100 iterations, we call the summed time $t_{100}$. FQED and mAEG show the same performance at this stage (Fig. 9 (b)).

At the next stage laser fields become strong enough, a QED cascade is triggered and a rapid growth of $N_{e}$ occurs ($1000<n_{it}<1400$). Strong field amplitudes and a not very short time step allow FQED to reveal its remarkable performance at this stage: AEG is up to 6 times slower (mAEG is only slightly slower).

When electron-positron plasma density becomes comparable with the relativistic critical density, the nonlinear self-consistent regime of laser interaction with extremely dense electron-positron plasma takes place ($1400<n_{it}<1600$). Field amplitudes sharply decrease inside dense plasma due to absorption and screening, so the speedup of particle loop computations with FQED and mAEG decreases (Fig. 9 (b)). Between FQED and mAEG the former shows slightly better performance.

The final stage is expansion of plasma due to decreasing of field intensity of incident pulses. Here performance becomes very similar with all three QED-event generators because $N_{e}$ strongly decreases (Fig. 9 (a)). We should note that in the case of mAEG $N_{e}$ shows deviation from that obtained with AEG and FQED. This discrepancy can be explained by several factors. First, $N_{e}(n_{it})$ is a random value, and, when the quantity of particles decreases, the dispersion of this random value increases. The reason of discrepancy of $N_{e}$ at $n_{it}\approx 1000$ (at the end of the deformation and the destruction of the initial spherical target) is similar. Second, resampling leads to different ensembles of macroparticles in simulations with different QED-event generators. As a result, the absorption of macroparticles with different statistical weight at the boundaries of the simulation box additionally increase dispersion of $N_{e}$.

As a result, a remarkable speedup of the particle loop stage with the FQED event generator during a third ($1000<n_{it}<1800$) of all iterations results in a decrease of the total particle loop time from 1952 s in the case of AEG to 786 s. This fact is demonstrated by the cumulative sum of $t_{100}$ in Fig. 9 (c).The mAEG event generator ensures performance (811 s for particle loop stage) very close to that in the case of FQED. The speedup for this problem is much less than that in the previous section and the reasons lie in the properties of interaction. Laser pulses have quite a short duration, so the maximal field amplitude is achieved for a relatively short period of time. The maximal fields are in quite a small region which particles escape relatively quickly. Also the dense generated electron-positron plasma decreases the field amplitude due to absorption and prevents penetration of strong fields inwards. Moreover, besides the particle loop stage there are other computational stages among which current deposition, communications and data output required approximately 900 s, 1100 s and 1100 s, respectively, independently on the type of QED-event generator. Thus the overall computational time demonstrates a less significant difference in performances: overall times were 6912 s, 5211 s and 5049 s with AEG, mAEG and FQED, respectively. Nevertheless, the newly developed QED-event generator can speed up simulations of such a problem, which is not very demanding for computation of QED processes, by 30%, while the speedup of the particle loop stage can achieve the factor of 2.5.

6 Conclusion

In this paper, we proposed a new algorithm for taking into account QED effects in resource-intensive strong-field quantum electrodynamics numerical simulations. Unlike other methods, our scheme allows us to avoid excessive splitting of the time step, which leads to a substantial speed up of simulations. Our implementation is publicly available as part of the open-source hi-$\chi$ project designed as a Python-controlled toolbox for collaborative development. We plan to further optimize the code for better utilization of computational resources of modern CPUs.

7 Authorship contribution statement

VV: Conceptualization, software development, validation, verification, writing draft, editing. JM: Conceptualization, validation, writing draft. AB: Validation, verification, writing draft, editing. AM: Validation, writing draft, editing. EE: Conceptualization, supervision, validation, verification, writing draft, editing. IM: Conceptualization, supervision, methodology, funding acquisition, writing draft, editing.

8 Acknowledgements

VV, AB, AM, EE, and IM acknowledge the use of computational resources provided by the Lobachevsky University and Joint Supercomputer Center of the Russian Academy of Sciences. VV and IM acknowledge the support of the Ministry of Science and Higher Education of the Russian Federation, project FSWR-2023-0034. AB and AM acknowledge the support of the Ministry of Science and Higher Education of the Russian Federation, project FFUF-2024-0030.

The authors would like to thank Arkady Gonoskov for useful discussions.

Appendix A Compton Rate Approximation

The partial rate $\mathrm{d}W/\mathrm{d}\varepsilon$ of nonlinear Compton scattering is given in equation 3. In order to approximate the total rate function $W(\chi)$ we considered it as a sum of two functions:

[TABLE]

where $C_{1}$ and $C_{2}$ correspond to terms in (Eq.3) containing the first and second synchrotron functions, respectively. Each term was approximated independently by means of a local trend approximation with refinement using a polynomial factor in order to approximate the trend-free function:

[TABLE]

This approach factors out the main trend of the initial function, which often results in simpler forms of the approximating polynomial. The main trends of the Compton rate were chosen based on its asymptotic form:

[TABLE]

Basic information about the functions $C_{1}(\chi)$ and $C_{2}(\chi)$ is presented in Tables 1 and 2, respectively. The exact form of the functions, including polynomial coefficients, is available in the open-source code hi-$\chi$ [53].

Appendix B Compton inverse CDF Approximation

As described earlier in Sec.3.3.2, in the main domain linear interpolation based on precomputed and tabulated values on the mesh is employed.

At small values of $\varepsilon$, the cumulative distribution function scales as

[TABLE]

Therefore, the Compton inverse CDF function at the lower bound of $r_{2}$ ($r_{2}<0.05$) is approximated by the following asymptotic form:

[TABLE]

The coefficient $D(\chi)$ is calculated from the conditions of the Compton inverse CDF function continuity at the closest mesh nodes.

Appendix C Breit-Wheeler Rate Approximation

In order to approximate the function $W_{p}(\chi_{\gamma})$ we consider it as a sum of two functions:

[TABLE]

corresponding to terms in (6) containing the first and second synchrotron functions, respectively. Each term is again approximated independently using the same approach as in A:

[TABLE]

The main trends of the Breit-Wheeler rate were chosen based on its asymptotic form:

[TABLE]

Basic information about the functions $B_{1}(\chi_{\gamma})$ and $B_{2}(\chi_{\gamma})$ is presented in Tables 3 and 4, respectively.

Appendix D Breit-Wheeler Inverse CDF Approximation

As described above in Sec.3.4.2, linear interpolation based on precomputed and tabulated values on the mesh is employed in the main domain.

At small values of $\varepsilon_{e}$, the cumulative distribution function scales as follows:

[TABLE]

Therefore, the Breit-Wheeler inverse CDF function at the lower bound of $r_{2}$ ($r_{2}<10^{-5}$) is approximated by the following asymptotic form:

[TABLE]

where the $LambertW$ function [54] is employed. The coefficients $L$ and $M$ are calculated from the conditions of function continuity on the closest mesh nodes.