Address review comments in #6141

[EZoni]

@aeriforme and I are opening this blank issue to keep track of a number of review comments in #6141 that we would like to review and address if necessary, in follow-up PRs.

1. Correctness / robustness (high priority)

LinearBreitWheelerUtil.H

• Near-threshold (s → 1) is unguarded. gamma_star = sqrt(s) can be marginally below 1 from round-off, making std::acosh(gamma_star) return NaN; and β → 0 makes every /beta in the integral and derivative a 0/0. Clamp with gamma_star = amrex::max(gamma_star, 1) and add an explicit small-β fallback (isotropic sampling, which is the correct β → 0 limit) below some threshold (e.g. β < 1e-3 or s - 1 < few·eps). This matters in single precision, where s - 1 loses all significance near threshold.

• The Lorentz-boost path (vc_sq > 0) is not covered by either new test. Both new inputs use head-on, equal-energy photons along x, so v_CM = 0, gc = 1, and the entire new boost-of-photon-1 + triad code is exercised only in the trivial n ∝ p1_lab branch. Add at least one test with (a) unequal photon energies and (b) a collision axis not aligned with a coordinate axis, so both the boost and the |nz| ≥ 0.9 branch are covered.

• Possible cancellation when computing n in a strongly boosted frame. bfact = (gc-1)/vc_sq·(v·p1) - gc·|p1|/c is a difference of large near-equal terms when photon 1 is strongly blueshifted in the lab frame; p1_CM can lose most of its significant digits. Evaluate the accuracy for large gc (single precision especially) and consider a reformulation, e.g. using the back-to-back relation between the two photons in the CM frame rather than boosting photon 1 alone.

• No convergence guarantee after the 20-iteration loop. The safeguarded Newton is fine in practice (the pdf in y varies only by ∼2× over the interval), but there is no post-loop check. Add an AMREX_ASSERT in debug builds on the final |G|, and document the worst-case error. Note that pure bisection from ±acosh(1e6) ≈ ±14.5 would need ∼50 iterations, so the loop count is only safe because Newton is expected to dominate.

• std::numeric_limits<ParticleReal>::min() is used as the "zero boost" test. This is inherited from the previous code, but it is now guarding a much larger code block; vc_sq below the smallest normal number is essentially "exactly zero". Consider a physically scaled threshold instead (relative to c²), and make sure the two vc_sq > min() tests stay in sync (they are now duplicated in two places).

• Missing includes. New uses of std::tanh, std::cosh, std::sinh, std::acosh, std::sqrt, std::numeric_limits in both LinearBreitWheelerUtil.H and LinearBreitWheelerCrossSection.H; add <cmath> and <limits> explicitly rather than relying on transitive AMReX includes.

• Which product gets the sampled angle. The angle is measured from photon 1 and assigned to the first product; this is unbiased only because dσ/dcosθ is even in cosθ. Add a one-line comment stating this, so a future non-symmetric extension (e.g. polarized photons) does not silently inherit a bias from the arbitrary photon-1/photon-2 ordering.

2. Simplification / performance (medium priority)

LinearBreitWheelerUtil.H

• g_star = sqrt(1 + p_star_sq/(me²c²)) is identical to gamma_star = sqrt(s). Reuse gamma_star and delete p_star_sq entirely (it has no other remaining use).

  Reuse gamma_star, remove p_star_sq in LinearBreitWheelerUtil.H #7244

• s is computed twice (powi<2>(2*E_ratio) appears in both p_star_sq and s). Compute once. Also make the amrex::Math::powi qualification consistent (one call is qualified, the neighbouring pre-existing one is not).

• Each Newton iteration calls LBWIntegralOfDifferentialCrossSectionTransformed and LBWDifferentialCrossSectionTransformed, which each recompute tanh(y) and cosh(y). Merge into a single value_and_derivative helper to halve the transcendental count in the hot loop (up to 20 iterations per pair-production event, on GPU).

• beta = std::tanh(std::acosh(gamma_star)) can be std::sqrt(1 - one_minus_beta2) (or sqrt(s-1)/sqrt(s)); check whether the transcendental form is actually needed for precision, and comment if so.

• Use amrex::Math::sincos(phi, ...) instead of separate std::cos/std::sin.

• LBWDifferentialCrossSectionTransformed computes sech2_y and then divides by it; use cosh_y*cosh_y directly and drop the division.

• Magic numbers: 20 (max iterations), 32 (epsilon multiplier), 0.9 (axis-selection threshold) should be named constexpr values with a short justification.

• px_star, py_star, pz_star are declared uninitialized far above their assignment; declare them const at the point of use.

• Rename I1 to something explicit (integral_at_ymax), and u to random_cdf_value.

• Consider factoring "sample CM-frame direction + boost to lab" into a shared utility, since the fusion / Compton / BW modules now duplicate a large part of this boost code.

• Check whether removing #include "Utils/ParticleUtils.H" breaks any other symbol used in this header (only RandomizeVelocity appears to have been used).

• Measure GPU register pressure / kernel occupancy before and after; this function grew substantially and is inside a particle kernel.

LinearBreitWheelerCrossSection.H

• Scalar arguments are taken by const amrex::ParticleReal&; pass by value (matches AMReX device-code practice, avoids reference-through-registers).

• Missing blank line before #endif at the end of the file.

• The two new functions are sampling helpers rather than cross sections; either rename the file or add a comment explaining why they live here.

3. Tests (medium priority)

analysis_angular.py

• The Python re-implementation of I(y) duplicates the C++ formula and can silently drift out of sync. Either add a small C++ unit test for LBWIntegralOfDifferentialCrossSectionTransformed / LBWDifferentialCrossSectionTransformed, or add a cross-reference comment in both files. The integral_matches assertion tests the analysis script's own algebra, not the simulation — move it to a unit test.

• Thresholds ks_distance < 0.01 and max_rel_err < 0.15 are hard-coded. 1.36/sqrt(N) for the sample sizes here is close to 0.01, so the KS check may be borderline/flaky. Derive the threshold from the effective sample size (weights are non-uniform, so use N_eff = (Σw)²/Σw²) and document the margin.

• Only electrons are checked. Add: positron distribution, electron/positron back-to-back-ness in the CM frame, and azimuthal uniformity of atan2(uz, uy) — the latter exercises the new phi sampling and the orthonormal-basis code, which is currently untested.

• rel_err = |hist - expected| / expected will divide by zero if any bin has zero expected counts; add a guard or a minimum-count mask.

• The theory curve uses np.gradient(cdf_fine, y_fine); use the analytic derivative instead (it already exists in C++, and mirrors what is actually sampled).

• Plotting is unconditional and adds a matplotlib dependency plus PNG output on every CI run; consider gating it behind a --plot flag.

• load_reduced_diagnostic / get_reduced_column duplicate helpers that exist elsewhere in WarpX analysis scripts; move to a shared module under Tools/ if a suitable location exists.

• The docstring hard-codes ux = +/-2.8 and +/-1.e6; it already reads them from the inputs, so drop the hard-coded values or mark them as "as configured in the shipped inputs".

Input files

• amr.n_cell = 8, amr.max_grid_size = 8 produces a single box while the test requests 2 MPI ranks — one rank idles and the domain decomposition is not exercised. Use max_grid_size = 4.

  Not true. I see this grids summary in both tests:

  Grids Summary:
     Level 0   2 grids  512 cells  100 % of domain
               smallest grid: 8 x 8 x 4  biggest grid: 8 x 8 x 4
  

• positron.do_not_deposit = 1 but electron has no do_not_deposit; with algo.maxwell_solver = none neither species needs deposition. Make them consistent.

• Add a comment in inputs_test_3d_linear_breit_wheeler_angular_relativistic explaining why density = 1e45 and ux = 1e6 (statistics for the strongly forward-peaked regime), so it is not mistaken for a physical setup.

• Consider adding a third case near threshold (ux slightly above 1) once item 1 is fixed — that is the regime most likely to produce NaNs.

Checksums

• Key ordering was changed inconsistently between the two regenerated benchmarks (lev=0 moved after electron in many_photons, moved to the top in two_photons), suggesting they were produced with different tooling. Regenerate all four with the standard WarpX checksum tool so ordering is uniform.

• Confirm that the existing two_photons / many_photons analysis scripts do not assert isotropy anywhere (only the checksums appear in this diff); if they do, update them.

• Newton iteration counts can differ between CPU/GPU libm, which could perturb these checksums; note the expected tolerance margin, and re-check if these tests are ever enabled on GPU CI.

4. Documentation (low priority)

• Docs/source/usage/parameters.rst: state that the polar angle is measured from the collision axis in the CM frame, and that the azimuthal angle is uniform for unpolarized photons.

• The following sentence, "The implementation follows the same numerical algorithm as that of fusion reactions", is now only true for the pairing/probability part, not the product kinematics. Qualify it.

• Check Docs/source/theory/multiphysics/collisions.rst (or equivalent) for a statement that BW emission is isotropic and update it; add the explicit dσ/d(cosθ) formula and a note on the y = atanh(β cosθ) change of variable used for sampling.

• LinearBreitWheelerUtil.H doxygen: parameter units say "normalized momentum ... (in m.s^-1)" — pre-existing error, worth fixing while touching the block.