Add differential cross-section in linear Breit Wheeler process

[aeriforme]

Overview

Until now, electron-positron pairs produced by the linear Breit-Wheeler module were emitted isotropically in the center-of-momentum (CM) frame. This PR samples the polar angle of the products from the Breit-Wheeler differential cross section, keeping the azimuthal angle uniform. The total cross section, and hence the pair
production rate, is unchanged.

Kinematics

For two colliding macrophotons with lab momenta $\boldsymbol{p}_1$, $\boldsymbol{p}_2$ colliding at an angle $\psi$, the energy of each photon in the CM frame follows from the Lorentz invariance of the total four-momentum norm,

$$ (E^\star)^2 = \dfrac{1}{2} \ c^2 |\boldsymbol{p}_1||\boldsymbol{p}_2|(1 - \cos\psi) . $$

By symmetry, each produced lepton carries energy $E^*$ in the CM frame, so its Lorentz factor, velocity and momentum are

$$ \gamma^\star = \dfrac{E^\star}{m_e c^2}, \qquad \beta = \sqrt{1 -( \gamma^\star)^{-2}}, \qquad |\boldsymbol{p}^\star| = m_e c \sqrt{(\gamma^\star)^2 - 1} . $$

Differential cross section

Define $x = \cos\theta$ the polar angle of the outgoing leptons with respect to the collision axis in the CM frame. The differential cross-section is (Berestetskii, Lifshitz & Pitaevskii, Quantum Electrodynamics, paragraph 86; see also Ribeyre et al., Plasma Phys. Control. Fusion 60, 104001, 2018):

$$ \dfrac{d\sigma}{d x} = \dfrac{\pi r_e^2}{2} \ \beta \ (1-\beta^2) \ f(\beta, x), \qquad f(\beta,x) = \dfrac{1 + 2\beta^2 - 2\beta^4 - 2\beta^2(1-\beta^2)x^2 - \beta^4 x^4} {\left(1 - \beta^2 x^2\right)^2} . $$

This is integrated analytically,

$$ F(x) = \int f \ dx = -\dfrac{(1-\beta^2)^2 \ x}{1-\beta^2 x^2} - x + \dfrac{3-\beta^4}{\beta} \ \mathrm{atanh}(\beta x), $$

and satisfies

$$ \beta\left[F(1) - F(-1)\right] = (3-\beta^4) \ \ln \ \frac{1+\beta}{1-\beta} + 2\beta(\beta^2-2), $$

which is the total cross section already implemented in LinearBreitWheelerCrossSection, i.e. $\int_{-1}^{1} (d\sigma/dx), dx = \sigma_{\rm tot}$ (see the corresponding code). This consistency is asserted in the new test.

Sampling

The polar angle is drawn by inverse-transform sampling of the exact cumulative distribution. To resolve the strong forward/backward peaking at high energy, the sampling is performed in the transformed variable

$$ y = \mathrm{atanh}(\beta \cos\theta) \in [-y_{\max}, \ y_{\max}], \qquad y_{\max} = \mathrm{atanh} \ \beta = \mathrm{arccosh} \gamma^\star, $$

in which $1 - \beta^2 x^2 = \mathrm{sech}^2 y$ and the (unnormalized) density and its integral read

$$ \frac{d\sigma}{dy} \propto g(\beta,y) = \frac{1}{\beta}\left[1 + \tanh^2 y + 2(1-\beta^2)\left(1 - (1-\beta^2)\cosh^2 y\right)\right], $$

$$ \tilde{F}(y) = \frac{1}{\beta}\left[-(1-\beta^2)^2 \tanh y \ \cosh^2 y - \tanh y + \left(3-\beta^4\right) y \right] \equiv F \left(\tanh y/\beta\right). $$

Since $g$ is even, $\tilde F$ is odd and the CDF is simply

$$ \mathrm{CDF}(y) = \frac{1}{2} + \frac{1}{2} \ \frac{\tilde F(y)}{\tilde F(y_{\max})} . $$

For $u \sim \mathcal{U}[0,1)$, the equation $\mathrm{CDF}(y) - u = 0$ is solved with a safeguarded Newton-Raphson iteration (analytic derivative $g/(2\tilde F(y_{\max}))$, bisection fallback whenever the Newton step leaves the min-max range, at most 20 iterations, tolerance $32 \ \epsilon$). Finally

$$ \cos\theta = \tanh(y)/\beta, \qquad \varphi \sim \mathcal{U}[0, 2\pi) . $$

Numerical robustness at high energy

For $\gamma^\star \gg 1$ the distribution collapses towards $|\cos\theta| \to 1$ and a naive evaluation loses all precision. This is avoided by computing $1-\beta^2 = (\gamma^\star)^{-2}$ directly (no subtractive cancellation), $\beta = \tanh(\mathrm{arccosh} \ \gamma^\star)$, and by sampling in $y$ rather than in $\cos\theta$.

Reconstruction of the momenta

The polar angle is measured from the collision axis, defined as the direction of the first photon in the CM frame. With the CM velocity $\boldsymbol{v}_c = c \ (\boldsymbol{p}_1+\boldsymbol{p}_2)/(|\boldsymbol{p}_1|+|\boldsymbol{p}_2|)$ and $\gamma_c = (1-v_c^2/c^2)^{-1/2}$,

$$ \boldsymbol{p}^\star_1 = \boldsymbol{p}_1 + \boldsymbol{v}_c \left[\frac{\gamma_c-1}{v_c^2} \ ( \boldsymbol{v}_c\cdot\boldsymbol{p}_1) - \frac{\gamma_c |\boldsymbol{p}_1|}{c}\right], \qquad \hat{\boldsymbol{n}} = \frac{\boldsymbol{p}^\star_1}{|\boldsymbol{p}^\star_1|} . $$

An orthonormal triad $(\boldsymbol{e}_1, \boldsymbol{e}_2, \hat{\boldsymbol{n}})$ is built around $\hat{\boldsymbol{n}}$ and the momentum of the first product in the CM frame is

$$ \boldsymbol{p}^\star = |\boldsymbol{p}^\star| \left(\sin\theta \ \cos\varphi \ \boldsymbol{e}_1 + \sin\theta \ \sin\varphi \ \boldsymbol{e}_2 + \cos\theta \ \hat{\boldsymbol{n}}\right), $$

which is then boosted back to the lab frame with the inverse transformation (Eq. 13 of F. Perez et al., Phys. Plasmas 19, 083104, 2012). The momentum of the second product follows from total momentum conservation, as before. Note that $d\sigma/d\cos\theta$ is even in $\cos\theta$, so the result is independent of which of the two photons defines the collision axis.

Tests

Two new tests, test_3d_linear_breit_wheeler_angular and test_3d_linear_breit_wheeler_angular_relativistic, collide two monoenergetic photon populations head-on along $x$ with $u_x = \pm 2.8$ and $u_x = \pm 10^6$ (in units of
$m_e c$), so that the CM frame coincides with the lab frame. analysis_angular.py checks:

1. total energy and momentum conservation from the reduced diagnostics;
2. the lepton Lorentz factor, $\gamma_{\rm lepton} = u_\gamma$;
3. that the weighted Kolmogorov-Smirnov distance between the sampled $\cos\theta$ and the analytical CDF is below $10^{-2}$;
4. that a 20-bin weighted histogram in $y$ agrees with the exact CDF differences to within 15%;
5. the consistency of $\tilde F$ with the closed-form angular integral.

[EZoni]

Here's a list of further comments that I propose we review and address, if necessary, in follow-up PRs: review_6141.md.

[EZoni]

Here's a list of further comments that I propose we review and address, if necessary, in follow-up PRs: review_6141.md.

Now tracked in #7241.