🔬 physics

Monte Carlo sampling of first-order QED processes in laser and pulsar plasmas

This paper introduces a novel Monte Carlo sampling method for first-order QED processes in laser and pulsar plasmas that replaces traditional numerical tabulation and root-finding with closed-form Padé approximants, enabling efficient and accurate event generation directly within radiative particle-in-cell codes.

Original authors: Verneri Sarjomaa, Joonas Nättilä

Published 2026-09-25

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Verneri Sarjomaa, Joonas Nättilä

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the most extreme corners of the universe, where gravity crushes matter into dense stars or where human-made lasers concentrate light to unimaginable intensities, the rules of everyday physics begin to fray. Here, the vacuum of space is not empty but teeming with activity, governed by a branch of physics known as quantum electrodynamics. This field describes how light and matter interact when electromagnetic fields are so strong that they can spontaneously create particles out of nothing. To understand these environments, scientists rely on computer simulations that act as virtual laboratories. These simulations must track individual particles, such as electrons and photons, as they race through intense fields, occasionally colliding or transforming into new particles. Because these events happen randomly and unpredictably, the computers use a method called Monte Carlo sampling, which is essentially a way of rolling digital dice to decide when an event occurs and how the energy is shared between the resulting particles.

For decades, running these simulations has been a slow and cumbersome task. To decide the outcome of a particle interaction, the computer must consult massive tables of pre-calculated numbers. These tables act like a dictionary, storing the complex probabilities of every possible event. When the simulation needs an answer, it looks up the closest match in the table and guesses the exact value by interpolating between the stored points. This process of looking up, guessing, and repeating is computationally expensive, slowing down the simulation and limiting how much of the universe researchers can model. The bottleneck is not the physics itself, but the way the computer accesses the data needed to perform the calculations.

A team of researchers at the University of Helsinki has found a way to bypass this bottleneck entirely. Instead of relying on massive tables of numbers, they have developed a set of mathematical shortcuts that describe these quantum processes using simple, standard functions. In their new approach, the computer does not need to search for a value or guess between two points. Instead, it can calculate the exact outcome of a particle interaction directly, using a formula that yields a precise answer in a single step. The researchers focused on two specific high-energy processes: synchrotron radiation, where an electron emits a photon as it spirals through a magnetic field, and the nonlinear Breit-Wheeler process, where a high-energy photon splits into an electron and a positron.

The core of their discovery lies in replacing the complex, tabulated probabilities with elegant approximations. In the old method, determining the energy of a new particle required solving a difficult equation that had no simple solution, forcing the computer to iterate through guesses until it found the right answer. The new method replaces this struggle with a direct calculation. The researchers constructed their approximations using a specific type of mathematical ratio that is flexible enough to match the complex behavior of the real physics but simple enough to be solved instantly. When the computer needs to know how much energy a new particle will carry, it simply plugs the current conditions into their formula and gets the answer immediately.

The results of this new method are remarkably accurate. When the researchers tested their formulas against the exact, traditional calculations, the difference was less than one percent across the entire range of conditions relevant to both laser experiments and astrophysical environments. In some cases, the error was even smaller, dropping to a fraction of a percent. This level of precision means that the new method does not sacrifice accuracy for speed. The researchers verified this by running millions of simulated events and comparing the distribution of the resulting particles to the theoretical predictions. The patterns produced by their new formulas matched the expected physics perfectly, confirming that the shortcuts do not distort the reality they are trying to model.

This advancement removes the need for the heavy data tables that have long slowed down these simulations. Because the new formulas require no storage space for pre-calculated data and no time-consuming lookups, they can be inserted directly into the most advanced simulation codes used by physicists today. This allows researchers to simulate complex scenarios, such as the behavior of plasma around neutron stars or the dynamics of matter in ultra-intense laser beams, with a level of detail and speed that was previously difficult to achieve. The method is robust enough to handle the extreme conditions found in nature, where quantum effects dominate, and it works equally well for the controlled environments of laboratory lasers.

The researchers also noted that their approach is general enough to be applied to other similar processes in quantum physics. While they focused on the two most common first-order interactions, the same technique could be adapted to describe other ways particles interact in strong fields. By turning a complex, iterative search into a direct calculation, this work simplifies the computational landscape of high-energy physics. It allows scientists to focus less on the mechanics of the simulation and more on the physical phenomena they are trying to understand, opening the door to more detailed and expansive models of the universe's most energetic events.

Technical Summary: Monte Carlo sampling of first-order QED processes in laser and pulsar plasmas

Problem Statement
Simulations of strong-field quantum electrodynamics (QED) in high-intensity laser experiments and astrophysical compact-object magnetospheres rely on Monte Carlo (MC) sampling to model stochastic processes like synchrotron radiation and the nonlinear Breit–Wheeler pair production. Standard MC algorithms require two random numbers per event: one to determine the time of occurrence (based on the total process rate) and another to determine the energy partition between produced particles (based on the cumulative probability distribution).

Currently, these simulations rely on pre-calculated lookup tables containing the total rates and cumulative probabilities, which involve integrals of modified Bessel functions. Evaluating these integrals at every time step is computationally prohibitive, necessitating tabulation and numerical inversion (often via interpolation or root-finding algorithms). This approach introduces overhead, requires memory for storage, and complicates the integration of QED modules into radiative particle-in-cell (PIC) codes.

Methodology
The authors propose replacing the tabulated functions and numerical inversion procedures with closed-form approximations constructed entirely from elementary functions. The methodology focuses on two specific first-order QED processes:

Synchrotron Radiation: An electron emitting a photon.
Nonlinear Breit–Wheeler Process: A photon decaying into an electron-positron pair.

The approach involves two main approximation steps:

Auxiliary Functions ( $T(\chi)$ ): The total rate of the process is governed by an auxiliary function $T(\chi)$ of the quantum parameter $\chi$ . The authors approximate the double Bessel integrals defining $T(\chi)$ using products of elementary functions (exponentials and polynomials) designed to match the known asymptotic limits at low and high $\chi$ .
Cumulative Probabilities ( $p(\chi, \chi_1)$ ): The probability distribution for the energy fraction $r = \chi_1/\chi$ of the produced particle is approximated using fourth-order Padé approximants. This specific functional form is chosen because it is flexible enough to track the exact probability across the full range of quantum parameters and, crucially, is analytically invertible.

By setting the cumulative probability equal to a random number $\zeta$ , the inversion of the Padé approximant reduces to solving a quartic equation. The authors provide explicit closed-form solutions for the roots of these equations, allowing the quantum parameter of the produced particle ( $\chi_1$ ) to be calculated directly without iterative root finding or table lookups.

The parameters for these approximations are determined by minimizing a weighted least-squares error against numerically computed exact values. The weighting scheme prioritizes regions of high particle production probability.

Key Contributions

Closed-Form Sampling: The paper derives explicit analytical formulas for the sampled quantum parameters ( $\chi_\gamma$ for synchrotron radiation and $\chi_e$ for Breit–Wheeler) that depend only on the incoming particle's quantum parameter and a random number.
Elimination of Tabulation: The proposed method removes the need for lookup tables, interpolation, and numerical root-finding algorithms in MC sampling of these specific QED processes.
Symmetrization Strategy: For the Breit–Wheeler process, the authors utilize the symmetry of the electron-positron pair production (interchangeability of the two particles) to restrict the approximation to the range $r \in [0.5, 1)$ , applying manual symmetrization for $r < 0.5$ .
Implementation Details: The paper provides specific coefficients for the Padé approximants and auxiliary function corrections, along with the explicit quartic solutions required for implementation.

Results
The authors validate the approximations through two independent tests:

Spectral Derivative Test: They compare the derivative of the approximated cumulative probability (which yields the differential spectrum) against the exact theoretical spectrum. The results show deviations of approximately 1% or less. Regions with larger relative errors correspond to spectral regions where the probability density is vanishingly small, keeping the absolute error negligible.
Direct Sampling Test: They generate histograms of $N = 9 \times 10^4$ sampled events using the closed-form formulas and compare them to the theoretical distributions. The sampled spectra agree with the exact theory at the percent level.

Specific accuracy metrics include:

Synchrotron Radiation ( $T_S$ ): Relative error stays below 0.7% for $\chi_e \in (10^{-3}, 10^2)$ and reaches 1.6% near $\chi_e = 5 \times 10^3$ .
Breit–Wheeler ( $T_{BW}$ ): Relative error is better than $3.3 \times 10^{-4}$ for $\chi_\gamma \gtrsim 1$ and remains below $5.5 \times 10^{-3}$ near $\chi_\gamma \approx 0.2$ .
Cumulative Probabilities: The sampled distributions for both processes match theoretical predictions within 1%.

Significance and Claims
The paper claims that these approximations enable the direct insertion of strong-field QED sampling into radiative particle-in-cell codes without the computational overhead of table lookups or iterative inversion. The authors emphasize that the method is general and can be extended to other first-order QED processes (such as photon absorption and pair annihilation) or different background fields (such as pure magnetic fields) with appropriate modifications to the quantum parameter definition.

The authors remain modest regarding performance claims, noting that they do not report specific speedup factors because the computational cost depends heavily on the host code, hardware, and the specific tabulation schemes being replaced. They also acknowledge that while single-process errors are small, the accumulation of these errors over a full simulation involving cascading processes (e.g., electrons emitting photons that create pairs) has not yet been quantified. Finally, they note a mathematical divergence at $\zeta = 1$ (where the produced particle takes all the energy), which is unphysical; they recommend clipping random numbers slightly away from unity to avoid floating-point errors.

More like this