Adversarial review: Beane threshold to Auger cubic-anisotropy bound
Date: 2026-09-09
Role: hostile physics reviewer
Target: lab/2026-09-08-auger-anisotropy/derive_bound.py
Context: reports/threads/2026-09-09-beane-scaling.md
Scratch code written: lab/2026-09-08-auger-anisotropy/adversary2_checks.py
Verdict
The claimed numerical bound b^-1 > 5.8e10 GeV is not a defensible physics bound. It is a conditional sensitivity estimate for a pure-proton, GZK-dominated, weakly-deflected universe with a one-parameter propagation response. Under those assumptions it is not obviously arithmetically wrong. Those assumptions are exactly the problem.
The two damaging failures are composition and the identification of Auger's spectral suppression with proton GZK attenuation. Either can erase the bound. The continuous-energy-loss objection, which I expected to be fatal, is subtler: stochastic 20% photopion losses are not an absorber, but if the observed factor-two flux loss really is due to that same interaction channel, the first-order response still scales with ln 2. The derivation is physically mislabeled, not automatically numerically killed.
Baseline used below: Argus response coefficient A = 0.74499 (b|p|)^2, 95% scan sensitivity A_95 = 0.3581, giving b|p| < 0.6933 and b^-1 > 5.77e10 GeV at |p| = 40 EeV.
Top 5 Objections Ranked By Damage
| rank |
objection |
damage to b^-1 > 5.8e10 GeV |
| 1 |
Composition invalidates the proton Delta-resonance chain. Auger Xmax/source fits favor mixed composition and decreasing proton flux with energy. For a nucleus, the relevant process is photodisintegration, not gamma_CMB + p -> Delta; the Beane eq. 18 proton threshold no longer maps to the observed event energy. |
Meaningless unless proton fraction and nuclear response are modeled. Optimistic nucleon-scaling gives pure He 0.25x, N 0.071x, Si 0.036x, Fe 0.018x: pure Fe would be ~1.0e9 GeV, a 56x weaker bound. Mixed response scales roughly as sqrt(sum f_A/A^2) before magnetic washout. |
| 2 |
D/lambda = ln 2 assumes the Auger suppression is proton GZK. Modern Auger fits allow low maximum rigidity / source cutoff plus mixed composition; the "disappointing model" explicitly removes CMB photopion GZK. |
If only fraction f_GZK of the log suppression is proton GZK, coefficient scales by f_GZK and bound scales by sqrt(f_GZK): f=0.5 -> 4.1e10 GeV, f=0.25 -> 2.9e10, f=0.1 -> 1.8e10, f=0 -> no bound. |
| 3 |
Magnetic deflection averages the lattice angle along the path. Eq. 18 uses instantaneous proton momentum relative to lattice axes. The proton does not propagate on a straight 100 Mpc line in general. |
For protons in B=1 nG, l_c=1 Mpc, D=100 Mpc, E=40 EeV, RMS deflection is about 6.3 deg, path-averaged l=4 damping 0.97, bound weaker 1.02x only. For 20-30 deg proton deflections seen in some IGM simulations, damping 0.73-0.54, bound weaker 1.17-1.36x. For ~63 deg (10 nG toy field), damping 0.165, bound weaker 2.46x. For nuclei, theta proportional to Z, so the signal is effectively washed out. |
| 4 |
The published 35.8% number is a null scan sensitivity, not a proper signal upper limit. Comparing a scan-max null quantile to a deterministic amplitude is approximate. |
Small. A known-orientation two-sided 95% null scale is 0.217, which would strengthen b^-1 by 1.29x; but a crude known-orientation Neyman inversion using observed T=0.2158 gives A_95 ~= 0.398, weakening the bound by 0.95x. Net damage is at most tens of percent, not the fatal issue. For unknown orientation, using the scan-max statistic is the right conservative comparison. |
| 5 |
The expansion is used outside its advertised regime, but the first omitted lattice-dispersion term is probably not the dominant error. Beane says ` |
p |
1. Steps 3-4: Amplification And Attenuation
1.1 Sign check: Argus has the sign right
For an absorber proxy F = exp(-tau), tau = D R/c. If the photon threshold increases, the interaction rate decreases:
+d omega/omega > 0 -> dR/R = -alpha d omega/omega -> dF/F = -d tau = +tau alpha d omega/omega
So higher threshold means longer attenuation length and more surviving flux. The sign in derive_bound.py is correct.
1.2 The Wien amplification is order-right but not the exact photopion kernel
Argus computes alpha = -d ln n(>omega)/d ln omega = 14.78 at E = 40.7 EeV, with n(>omega) proportional to int_x^infty u^2/(exp(u)-1) du, x = omega/kT = 16.66.
For an isotropic photon bath and a narrow resonance, the standard energy-loss/rate kernel is closer to
R proportional to int_{epsilon'_thr}^infty d epsilon' sigma(epsilon') ...
int_{epsilon'/(2 gamma)}^infty d epsilon n(epsilon)/epsilon^2
not simply the number density above threshold. In the narrow-resonance limit the inner blackbody integral is int_x^infty du/(exp(u)-1) = -ln(1-exp(-x)), whose logarithmic derivative at x=16.66 is essentially 16.66, not 14.78. That would strengthen the response by 16.66/14.78 = 1.13 and strengthen b^-1 by sqrt(1.13) = 1.06.
This does not break the bound; it shows the coefficient is a hand-built approximation. A real calculation needs SOPHIA/CRPropa/SimProp-style photopion kernels and inelasticities, not a bare CMB photon-count tail. Sources: UHECR propagation review by Allard, arXiv:astro-ph/0501124, emphasizes energy losses and magnetic fields; CRPropa 3.2 is explicitly a propagation framework for high-energy particles, arXiv:2208.00107; SimProp is a Monte Carlo UHECR propagation code, arXiv:1602.01239 and arXiv:1705.03729.
1.3 Continuous energy loss: not fatal if, and only if, the suppression is GZK
The physical objection is real: photopion production degrades a proton by an inelasticity of order kappa ~ 0.2; it does not delete the proton. The absorber label lambda = 1/R is wrong if read literally.
However, the first-order anisotropy response survives in a simple stochastic-loss model if the observed flux suppression is caused by the same interaction channel. For an integral flux above a threshold with source differential index gamma, Poisson mean interaction count mu, and fractional energy retained per interaction (1-kappa),
F(>E) / F0(>E) = sum_n Pois(n;mu) (1-kappa)^{n(gamma-1)}
= exp{ mu[(1-kappa)^{gamma-1} - 1] }.
For Auger's above-ankle gamma ~= 2.6 and kappa = 0.2, a factor-two suppression implies mu = ln 2 / [1 - 0.8^1.6] = 2.31 interactions, not 0.693. But the logarithmic response to a fractional rate change is still
d ln F = -ln(2) dmu/mu.
So the "20% loss, not absorption" objection changes the physical interpretation of D/lambda; it does not by itself change the coefficient if the measured factor-two suppression is genuinely proton GZK. Pile-up below the cut is already what the integral-flux expression accounts for.
Damage factor: ~1.0 under the pure-proton/GZK assumption. The derivation survives this specific attack more than I expected.
1.4 The fatal propagation assumption is not absorption; it is causation
Argus reads D/lambda = ln 2 from Auger arXiv:1002.1975, which reports the spectrum above the ankle is E^-2.6 followed by flux suppression and factor two at log10(E/eV)=19.61. That is a measured spectral fact. It is not a measured proton-GZK optical depth.
Modern Auger source/composition fits explicitly reproduce the data with mixed composition, hard spectra, and low rigidity cutoffs. Auger arXiv:2211.02857 says the data above the ankle are well reproduced by sources emitting mixed composition with a hard spectrum and a low rigidity cutoff. Auger arXiv:1705.03807 similarly reports fits with relatively low maximum injection energies and hard spectral indices. Berezinsky et al.'s "disappointing model" arXiv:0907.5194 is sharper: with average energy per nucleon below 2-4 EeV, CMB pion photoproduction is absent and "GZK cutoff in the spectrum does not exist."
If the suppression is source maximum rigidity rather than propagation, a threshold anisotropy in CMB photopion production produces no flux anisotropy at 40 EeV. The response coefficient is multiplied by the fraction of the observed log suppression actually due to proton GZK losses, f_GZK:
A = f_GZK * 0.74499 (b|p|)^2
b^-1 limit = 5.77e10 GeV * sqrt(f_GZK)
Numerically: f_GZK=0.5 -> 4.08e10 GeV; 0.25 -> 2.88e10; 0.10 -> 1.82e10; 0 -> no bound. This is one of the two largest failures.
1.5 Magnetic deflection and path averaging
For turbulent extragalactic fields, a standard scaling is
theta_rms ~= 0.8 deg * Z * (B/nG) * sqrt(D/10 Mpc) * sqrt(l_c/1 Mpc) * (100 EeV/E).
At E=40 EeV, D=100 Mpc, l_c=1 Mpc, B=1 nG, this gives theta_rms ~= 6.3 deg for a proton. An l-th multipole convolved with small-angle diffusion is damped approximately as exp[-l(l+1) theta^2/4]; for l=4 this is exp[-5 theta^2]. For signal accumulated along the trajectory and compared to the arrival direction, a simple path average gives
W(theta) = [1 - exp(-5 theta^2)] / [5 theta^2].
Computed damping and bound weakening:
| total proton RMS deflection |
l=4 path damping W |
bound weaker by 1/sqrt(W) |
| 3 deg |
0.993 |
1.00x |
| 6.3 deg |
0.970 |
1.02x |
| 10 deg |
0.928 |
1.04x |
| 19 deg |
0.769 |
1.14x |
| 30 deg |
0.544 |
1.36x |
| 63 deg |
0.165 |
2.46x |
Sokolenko et al. arXiv:2101.07207 find feedback-generated intergalactic magnetic bubbles can deflect about half of 10^20 eV protons by >=1 deg and some by 20-30 deg. At 40 EeV those deflections are larger by 2.5x at fixed rigidity unless the path/field environment differs. For protons in quiet void fields, deflection is not fatal. For structured fields or nuclei, the l=4 signal is heavily washed out.
2. Steps 1-2: Harmonic Normalisation
I could not break this. With standard Condon-Shortley, orthonormal complex spherical harmonics,
Y_4^0 = 3/(16 sqrt(pi)) (35 z^4 - 30 z^2 + 3)
Y_4^4 + Y_4^-4 = 2 * (3/16) sqrt(35/(2 pi)) (x^4 - 6x^2 y^2 + y^4)
Using sin^4(theta) cos(4 phi) = x^4 - 6x^2 y^2 + y^4,
Y_4^0 + sqrt(5/14)(Y_4^4 + Y_4^-4)
= 15/(4 sqrt(pi)) [x^4 + y^4 + z^4 - 3/5].
Since Argus defines
K4 = [x^4 + y^4 + z^4 - 3/5] / sqrt(16/525),
the conversion factor is
Ytilde = [15/(4 sqrt(pi)) sqrt(16/525)] K4 = 0.369348781845 K4.
Then
(sqrt(pi)/9) * 0.369348781845 = 0.072739296745.
This matches derive_bound.py. The convention caveat is only this: if Beane's Y_4^{+4}+Y_4^{-4} notation meant real tesseral harmonics already normalized differently, the coefficient would change by sqrt(2)-type factors. But their displayed combination is consistent with the standard complex-harmonic sum, and the Cartesian identity is exact.
Damage factor: 1.0. This step survives.
3. Steps 5-6: Statistical Inversion
The user's worry was: Argus compares a 95% upper limit on a scanned max |S| statistic to a predicted amplitude at a specific orientation. Is that wrong?
Mostly no. The lattice orientation is unknown, and Argus's search statistic is exactly max |S(R)| over orientations. If the prediction is "there exists some lattice orientation," then a scan-calibrated limit is the conservative object to use. It is not too strong relative to a single-orientation limit; it is weaker.
Quantitatively, at E >= 40 EeV, Argus has N=87, per-orientation noise sigma ~= 1.0324/sqrt(87) = 0.1107. A known-orientation two-sided 95% null scale would be 1.96 sigma = 0.2169, while the scan-null 95th percentile is 0.3581. The scan penalty is therefore 0.3581/0.2169 = 1.65 in amplitude. Because b^-1 proportional to sqrt(A_response/A_limit), a known-orientation bound would be stronger by sqrt(1.65) = 1.29: 7.4e10 GeV instead of 5.8e10 GeV.
The real statistical criticism is different: the 35.8% number is the 95th percentile of the null scan maximum. That is a sensitivity scale, not a strict Neyman/Feldman-Cousins upper limit on a positive signal amplitude. Using the observed scan value T_obs=0.2158 and a crude known-orientation Gaussian inversion gives A_95 ~= 0.398, not 0.358; that would weaken b^-1 by sqrt(0.358/0.398)=0.95. Full scan-signal injection would be needed for a publication-grade limit.
Damage factor: ~0.95x for proper-CL pedantry; 1.29x stronger if orientation were known. This is not a fatal objection.
4. Expansion Parameter
Beane eq. 18 is explicitly for |p| << 1/b. Argus inverts to b|p| < 0.693, so the derived limit sits at the edge of the expansion's domain.
I quantified the next term for the simple massless lattice dispersion, not the full proton/Delta threshold. For
sinh^2(bE/2) = sum_i sin^2(bp_i/2)
the expansion of E/p has an l=4 leading anisotropy proportional to q^2 (sum n_i^4), q=bp. Projecting the q^4 terms back onto the same l=4 cubic harmonic gives
next/leading ~= 0.0979 q^2.
So the projected omitted l=4 piece is 1.6% at q=0.4 and 4.8% at q=0.7. Exact massless dispersion values give additional scale context:
q=bp |
axis E/p |
diagonal E/p |
axis - diagonal |
| 0.4 |
0.986926 |
0.991273 |
-0.004347 |
| 0.7 |
0.961459 |
0.974220 |
-0.012761 |
| 1.0 |
0.925468 |
0.950036 |
-0.024567 |
This narrow calculation says the formal lattice sine expansion is not exploding at q=0.7; the next l=4 correction would move the bound only a few percent if the proton threshold inherited the same coefficients. But that is not enough. Eq. 18 is a threshold formula for composite proton/Delta kinematics under assumptions Beane themselves flag. At q~0.7, a full threshold derivation from the exact lattice dispersion for proton, Delta, and photon is required. Without it, quote the result as order-of-magnitude, not a bound.
Damage factor: formal next l=4 term ~0.98-1.03x on b^-1; uncontrolled threshold-systematic could be order unity.
5. Composition
This is the most damaging physics objection. Argus assumes every E >= 40 EeV event is a proton undergoing gamma_CMB + p -> Delta. Auger composition data do not support that as a safe assumption.
Auger Xmax analysis arXiv:1701.06812 reports the data are better described by a mix of light, intermediate, and heavy primaries, and that the proton flux decreases with energy. Auger source fits arXiv:2211.02857 reproduce spectrum and composition above the ankle with mixed composition, hard spectra, and low rigidity cutoff. The 2017 Auger astrophysical interpretation arXiv:1705.03807 reaches the same broad conclusion: relatively low maximum injection energies and hard spectra. Allard arXiv:0906.3156 emphasizes that spectrum alone cannot distinguish source-composition models; Xmax is needed to break the degeneracy.
For nuclei, the dominant energy-loss process is photodisintegration on CMB/EBL photons, not the Delta resonance of a free proton. The relevant Lorentz factor is E/A, and the relevant cross sections are giant dipole resonance / quasi-deuteron / nuclear photopion channels, not Beane's eq. 18. CRPropa/SimProp uncertainty studies explicitly treat proton, nitrogen, and iron scenarios and photodisintegration choices as major ingredients; see arXiv:1508.01824.
An optimistic ceiling if one pretended the bound applies to bound nucleons independently:
response_A proportional to (b E/A)^2
b^-1 limit_A = b^-1_proton / A
Numerically:
| primary |
A |
optimistic b^-1 |
factor vs proton |
| p |
1 |
5.77e10 GeV |
1.000 |
| He |
4 |
1.44e10 GeV |
0.250 |
| N |
14 |
4.12e9 GeV |
0.071 |
| Si |
28 |
2.06e9 GeV |
0.036 |
| Fe |
56 |
1.03e9 GeV |
0.018 |
For a mixture, before magnetic smearing, the response coefficient would scale approximately as sum_A f_A/A^2, so the b^-1 bound scales as sqrt(sum_A f_A/A^2). Example: 25% p + 75% N gives sqrt(0.25 + 0.75/196) = 0.504, so 2.9e10 GeV. If the proton fraction above 40 EeV is small, the claimed bound rapidly becomes irrelevant. If the events are mostly intermediate/heavy nuclei and the physics is photodisintegration, the Beane proton-Delta conversion is simply the wrong calculation.
Damage factor: 0 to 0.5x for plausible source-cutoff/mixed-composition cases; 1/56 in the pure-iron optimistic ceiling; more honestly, no bound without a nuclear lattice-threshold derivation.
Sources
- Beane, Davoudi, Savage, "Constraints on the Universe as a Numerical Simulation," arXiv:1210.1847, https://arxiv.org/abs/1210.1847
- Pierre Auger Collaboration, "Measurement of the energy spectrum of cosmic rays above 10^18 eV using the Pierre Auger Observatory," arXiv:1002.1975, https://arxiv.org/abs/1002.1975
- Allard, "Ultra-high energy cosmic ray propagation in the Universe," arXiv:astro-ph/0501124, https://arxiv.org/html/astro-ph/0501124
- Allard, "Propagation of extragalactic ultra-high energy cosmic-ray nuclei," arXiv:0906.3156, https://arxiv.org/abs/0906.3156
- Berezinsky et al., "Ultra High Energy Cosmic Rays: The disappointing model," arXiv:0907.5194, https://arxiv.org/abs/0907.5194
- Pierre Auger Collaboration, "Measurements of the depth of maximum... composition implications," arXiv:1701.06812, https://arxiv.org/abs/1701.06812
- Pierre Auger Collaboration, "Astrophysical interpretation of Pierre Auger Observatory measurements...," arXiv:1705.03807, https://arxiv.org/abs/1705.03807
- Pierre Auger Collaboration, "Constraining the sources of ultra-high-energy cosmic rays across and above the ankle...," arXiv:2211.02857, https://arxiv.org/abs/2211.02857
- Alves Batista et al., "CRPropa 3.2," arXiv:2208.00107, https://arxiv.org/abs/2208.00107
- Aloisio et al., "SimProp v2r3," arXiv:1602.01239, https://arxiv.org/abs/1602.01239
- Aloisio et al., "SimProp v2r4," arXiv:1705.03729, https://arxiv.org/abs/1705.03729
- di Matteo et al., "Effects of uncertainties in simulations of extragalactic UHECR propagation, using CRPropa and SimProp," arXiv:1508.01824, https://arxiv.org/abs/1508.01824
- Sokolenko et al., "Ultra-high energy cosmic rays deflection by the Intergalactic Magnetic Field," arXiv:2101.07207, https://arxiv.org/abs/2101.07207
View exactly as delivered (raw text)
# Adversarial review: Beane threshold to Auger cubic-anisotropy bound
Date: 2026-09-09
Role: hostile physics reviewer
Target: `lab/2026-09-08-auger-anisotropy/derive_bound.py`
Context: `reports/threads/2026-09-09-beane-scaling.md`
Scratch code written: `lab/2026-09-08-auger-anisotropy/adversary2_checks.py`
## Verdict
The claimed numerical bound `b^-1 > 5.8e10 GeV` is not a defensible physics bound. It is a conditional sensitivity estimate for a pure-proton, GZK-dominated, weakly-deflected universe with a one-parameter propagation response. Under those assumptions it is not obviously arithmetically wrong. Those assumptions are exactly the problem.
The two damaging failures are composition and the identification of Auger's spectral suppression with proton GZK attenuation. Either can erase the bound. The continuous-energy-loss objection, which I expected to be fatal, is subtler: stochastic 20% photopion losses are not an absorber, but if the observed factor-two flux loss really is due to that same interaction channel, the first-order response still scales with `ln 2`. The derivation is physically mislabeled, not automatically numerically killed.
Baseline used below: Argus response coefficient `A = 0.74499 (b|p|)^2`, 95% scan sensitivity `A_95 = 0.3581`, giving `b|p| < 0.6933` and `b^-1 > 5.77e10 GeV` at `|p| = 40 EeV`.
## Top 5 Objections Ranked By Damage
| rank | objection | damage to `b^-1 > 5.8e10 GeV` |
|---|---|---:|
| 1 | **Composition invalidates the proton Delta-resonance chain.** Auger Xmax/source fits favor mixed composition and decreasing proton flux with energy. For a nucleus, the relevant process is photodisintegration, not `gamma_CMB + p -> Delta`; the Beane eq. 18 proton threshold no longer maps to the observed event energy. | Meaningless unless proton fraction and nuclear response are modeled. Optimistic nucleon-scaling gives pure He `0.25x`, N `0.071x`, Si `0.036x`, Fe `0.018x`: pure Fe would be `~1.0e9 GeV`, a 56x weaker bound. Mixed response scales roughly as `sqrt(sum f_A/A^2)` before magnetic washout. |
| 2 | **`D/lambda = ln 2` assumes the Auger suppression is proton GZK.** Modern Auger fits allow low maximum rigidity / source cutoff plus mixed composition; the "disappointing model" explicitly removes CMB photopion GZK. | If only fraction `f_GZK` of the log suppression is proton GZK, coefficient scales by `f_GZK` and bound scales by `sqrt(f_GZK)`: `f=0.5 -> 4.1e10 GeV`, `f=0.25 -> 2.9e10`, `f=0.1 -> 1.8e10`, `f=0 -> no bound`. |
| 3 | **Magnetic deflection averages the lattice angle along the path.** Eq. 18 uses instantaneous proton momentum relative to lattice axes. The proton does not propagate on a straight 100 Mpc line in general. | For protons in `B=1 nG`, `l_c=1 Mpc`, `D=100 Mpc`, `E=40 EeV`, RMS deflection is about `6.3 deg`, path-averaged l=4 damping `0.97`, bound weaker `1.02x` only. For `20-30 deg` proton deflections seen in some IGM simulations, damping `0.73-0.54`, bound weaker `1.17-1.36x`. For `~63 deg` (`10 nG` toy field), damping `0.165`, bound weaker `2.46x`. For nuclei, `theta proportional to Z`, so the signal is effectively washed out. |
| 4 | **The published 35.8% number is a null scan sensitivity, not a proper signal upper limit.** Comparing a scan-max null quantile to a deterministic amplitude is approximate. | Small. A known-orientation two-sided 95% null scale is `0.217`, which would strengthen `b^-1` by `1.29x`; but a crude known-orientation Neyman inversion using observed `T=0.2158` gives `A_95 ~= 0.398`, weakening the bound by `0.95x`. Net damage is at most tens of percent, not the fatal issue. For unknown orientation, using the scan-max statistic is the right conservative comparison. |
| 5 | **The expansion is used outside its advertised regime, but the first omitted lattice-dispersion term is probably not the dominant error.** Beane says `|p| << 1/b`; Argus's bound allows `b|p| = 0.69`. | In a massless lattice-dispersion expansion, the projected next l=4 term is `~0.0979 (bp)^2` times the leading l=4 term: `1.6%` at `bp=0.4`, `4.8%` at `bp=0.7`. That alone would move `b^-1` only `~2-3%`. The real problem is not the formal sine expansion; it is that the proton/Delta threshold and propagation response are not recalculated at `bp ~ 0.7`. Treat as uncontrolled order-unity systematic until derived. |
## 1. Steps 3-4: Amplification And Attenuation
### 1.1 Sign check: Argus has the sign right
For an absorber proxy `F = exp(-tau)`, `tau = D R/c`. If the photon threshold increases, the interaction rate decreases:
```text
+d omega/omega > 0 -> dR/R = -alpha d omega/omega -> dF/F = -d tau = +tau alpha d omega/omega
```
So higher threshold means longer attenuation length and more surviving flux. The sign in `derive_bound.py` is correct.
### 1.2 The Wien amplification is order-right but not the exact photopion kernel
Argus computes `alpha = -d ln n(>omega)/d ln omega = 14.78` at `E = 40.7 EeV`, with `n(>omega) proportional to int_x^infty u^2/(exp(u)-1) du`, `x = omega/kT = 16.66`.
For an isotropic photon bath and a narrow resonance, the standard energy-loss/rate kernel is closer to
```text
R proportional to int_{epsilon'_thr}^infty d epsilon' sigma(epsilon') ...
int_{epsilon'/(2 gamma)}^infty d epsilon n(epsilon)/epsilon^2
```
not simply the number density above threshold. In the narrow-resonance limit the inner blackbody integral is `int_x^infty du/(exp(u)-1) = -ln(1-exp(-x))`, whose logarithmic derivative at `x=16.66` is essentially `16.66`, not `14.78`. That would strengthen the response by `16.66/14.78 = 1.13` and strengthen `b^-1` by `sqrt(1.13) = 1.06`.
This does not break the bound; it shows the coefficient is a hand-built approximation. A real calculation needs SOPHIA/CRPropa/SimProp-style photopion kernels and inelasticities, not a bare CMB photon-count tail. Sources: UHECR propagation review by Allard, arXiv:astro-ph/0501124, emphasizes energy losses and magnetic fields; CRPropa 3.2 is explicitly a propagation framework for high-energy particles, arXiv:2208.00107; SimProp is a Monte Carlo UHECR propagation code, arXiv:1602.01239 and arXiv:1705.03729.
### 1.3 Continuous energy loss: not fatal if, and only if, the suppression is GZK
The physical objection is real: photopion production degrades a proton by an inelasticity of order `kappa ~ 0.2`; it does not delete the proton. The absorber label `lambda = 1/R` is wrong if read literally.
However, the first-order anisotropy response survives in a simple stochastic-loss model if the observed flux suppression is caused by the same interaction channel. For an integral flux above a threshold with source differential index `gamma`, Poisson mean interaction count `mu`, and fractional energy retained per interaction `(1-kappa)`,
```text
F(>E) / F0(>E) = sum_n Pois(n;mu) (1-kappa)^{n(gamma-1)}
= exp{ mu[(1-kappa)^{gamma-1} - 1] }.
```
For Auger's above-ankle `gamma ~= 2.6` and `kappa = 0.2`, a factor-two suppression implies `mu = ln 2 / [1 - 0.8^1.6] = 2.31` interactions, not `0.693`. But the logarithmic response to a fractional rate change is still
```text
d ln F = -ln(2) dmu/mu.
```
So the "20% loss, not absorption" objection changes the physical interpretation of `D/lambda`; it does **not** by itself change the coefficient if the measured factor-two suppression is genuinely proton GZK. Pile-up below the cut is already what the integral-flux expression accounts for.
Damage factor: `~1.0` under the pure-proton/GZK assumption. The derivation survives this specific attack more than I expected.
### 1.4 The fatal propagation assumption is not absorption; it is causation
Argus reads `D/lambda = ln 2` from Auger arXiv:1002.1975, which reports the spectrum above the ankle is `E^-2.6` followed by flux suppression and factor two at `log10(E/eV)=19.61`. That is a measured spectral fact. It is not a measured proton-GZK optical depth.
Modern Auger source/composition fits explicitly reproduce the data with mixed composition, hard spectra, and low rigidity cutoffs. Auger arXiv:2211.02857 says the data above the ankle are well reproduced by sources emitting mixed composition with a hard spectrum and a low rigidity cutoff. Auger arXiv:1705.03807 similarly reports fits with relatively low maximum injection energies and hard spectral indices. Berezinsky et al.'s "disappointing model" arXiv:0907.5194 is sharper: with average energy per nucleon below `2-4 EeV`, CMB pion photoproduction is absent and "GZK cutoff in the spectrum does not exist."
If the suppression is source maximum rigidity rather than propagation, a threshold anisotropy in CMB photopion production produces no flux anisotropy at 40 EeV. The response coefficient is multiplied by the fraction of the observed log suppression actually due to proton GZK losses, `f_GZK`:
```text
A = f_GZK * 0.74499 (b|p|)^2
b^-1 limit = 5.77e10 GeV * sqrt(f_GZK)
```
Numerically: `f_GZK=0.5 -> 4.08e10 GeV`; `0.25 -> 2.88e10`; `0.10 -> 1.82e10`; `0 -> no bound`. This is one of the two largest failures.
### 1.5 Magnetic deflection and path averaging
For turbulent extragalactic fields, a standard scaling is
```text
theta_rms ~= 0.8 deg * Z * (B/nG) * sqrt(D/10 Mpc) * sqrt(l_c/1 Mpc) * (100 EeV/E).
```
At `E=40 EeV`, `D=100 Mpc`, `l_c=1 Mpc`, `B=1 nG`, this gives `theta_rms ~= 6.3 deg` for a proton. An l-th multipole convolved with small-angle diffusion is damped approximately as `exp[-l(l+1) theta^2/4]`; for l=4 this is `exp[-5 theta^2]`. For signal accumulated along the trajectory and compared to the arrival direction, a simple path average gives
```text
W(theta) = [1 - exp(-5 theta^2)] / [5 theta^2].
```
Computed damping and bound weakening:
| total proton RMS deflection | l=4 path damping `W` | bound weaker by `1/sqrt(W)` |
|---:|---:|---:|
| 3 deg | 0.993 | 1.00x |
| 6.3 deg | 0.970 | 1.02x |
| 10 deg | 0.928 | 1.04x |
| 19 deg | 0.769 | 1.14x |
| 30 deg | 0.544 | 1.36x |
| 63 deg | 0.165 | 2.46x |
Sokolenko et al. arXiv:2101.07207 find feedback-generated intergalactic magnetic bubbles can deflect about half of `10^20 eV` protons by `>=1 deg` and some by `20-30 deg`. At 40 EeV those deflections are larger by `2.5x` at fixed rigidity unless the path/field environment differs. For protons in quiet void fields, deflection is not fatal. For structured fields or nuclei, the l=4 signal is heavily washed out.
## 2. Steps 1-2: Harmonic Normalisation
I could not break this. With standard Condon-Shortley, orthonormal complex spherical harmonics,
```text
Y_4^0 = 3/(16 sqrt(pi)) (35 z^4 - 30 z^2 + 3)
Y_4^4 + Y_4^-4 = 2 * (3/16) sqrt(35/(2 pi)) (x^4 - 6x^2 y^2 + y^4)
```
Using `sin^4(theta) cos(4 phi) = x^4 - 6x^2 y^2 + y^4`,
```text
Y_4^0 + sqrt(5/14)(Y_4^4 + Y_4^-4)
= 15/(4 sqrt(pi)) [x^4 + y^4 + z^4 - 3/5].
```
Since Argus defines
```text
K4 = [x^4 + y^4 + z^4 - 3/5] / sqrt(16/525),
```
the conversion factor is
```text
Ytilde = [15/(4 sqrt(pi)) sqrt(16/525)] K4 = 0.369348781845 K4.
```
Then
```text
(sqrt(pi)/9) * 0.369348781845 = 0.072739296745.
```
This matches `derive_bound.py`. The convention caveat is only this: if Beane's `Y_4^{+4}+Y_4^{-4}` notation meant real tesseral harmonics already normalized differently, the coefficient would change by `sqrt(2)`-type factors. But their displayed combination is consistent with the standard complex-harmonic sum, and the Cartesian identity is exact.
Damage factor: `1.0`. This step survives.
## 3. Steps 5-6: Statistical Inversion
The user's worry was: Argus compares a 95% upper limit on a scanned `max |S|` statistic to a predicted amplitude at a specific orientation. Is that wrong?
Mostly no. The lattice orientation is unknown, and Argus's search statistic is exactly `max |S(R)|` over orientations. If the prediction is "there exists some lattice orientation," then a scan-calibrated limit is the conservative object to use. It is not too strong relative to a single-orientation limit; it is weaker.
Quantitatively, at `E >= 40 EeV`, Argus has `N=87`, per-orientation noise `sigma ~= 1.0324/sqrt(87) = 0.1107`. A known-orientation two-sided 95% null scale would be `1.96 sigma = 0.2169`, while the scan-null 95th percentile is `0.3581`. The scan penalty is therefore `0.3581/0.2169 = 1.65` in amplitude. Because `b^-1 proportional to sqrt(A_response/A_limit)`, a known-orientation bound would be stronger by `sqrt(1.65) = 1.29`: `7.4e10 GeV` instead of `5.8e10 GeV`.
The real statistical criticism is different: the 35.8% number is the 95th percentile of the null scan maximum. That is a sensitivity scale, not a strict Neyman/Feldman-Cousins upper limit on a positive signal amplitude. Using the observed scan value `T_obs=0.2158` and a crude known-orientation Gaussian inversion gives `A_95 ~= 0.398`, not `0.358`; that would weaken `b^-1` by `sqrt(0.358/0.398)=0.95`. Full scan-signal injection would be needed for a publication-grade limit.
Damage factor: `~0.95x` for proper-CL pedantry; `1.29x stronger` if orientation were known. This is not a fatal objection.
## 4. Expansion Parameter
Beane eq. 18 is explicitly for `|p| << 1/b`. Argus inverts to `b|p| < 0.693`, so the derived limit sits at the edge of the expansion's domain.
I quantified the next term for the simple massless lattice dispersion, not the full proton/Delta threshold. For
```text
sinh^2(bE/2) = sum_i sin^2(bp_i/2)
```
the expansion of `E/p` has an l=4 leading anisotropy proportional to `q^2 (sum n_i^4)`, `q=bp`. Projecting the `q^4` terms back onto the same l=4 cubic harmonic gives
```text
next/leading ~= 0.0979 q^2.
```
So the projected omitted l=4 piece is `1.6%` at `q=0.4` and `4.8%` at `q=0.7`. Exact massless dispersion values give additional scale context:
| `q=bp` | axis `E/p` | diagonal `E/p` | axis - diagonal |
|---:|---:|---:|---:|
| 0.4 | 0.986926 | 0.991273 | -0.004347 |
| 0.7 | 0.961459 | 0.974220 | -0.012761 |
| 1.0 | 0.925468 | 0.950036 | -0.024567 |
This narrow calculation says the formal lattice sine expansion is not exploding at `q=0.7`; the next l=4 correction would move the bound only a few percent if the proton threshold inherited the same coefficients. But that is not enough. Eq. 18 is a threshold formula for composite proton/Delta kinematics under assumptions Beane themselves flag. At `q~0.7`, a full threshold derivation from the exact lattice dispersion for proton, Delta, and photon is required. Without it, quote the result as order-of-magnitude, not a bound.
Damage factor: formal next l=4 term `~0.98-1.03x` on `b^-1`; uncontrolled threshold-systematic could be order unity.
## 5. Composition
This is the most damaging physics objection. Argus assumes every `E >= 40 EeV` event is a proton undergoing `gamma_CMB + p -> Delta`. Auger composition data do not support that as a safe assumption.
Auger Xmax analysis arXiv:1701.06812 reports the data are better described by a mix of light, intermediate, and heavy primaries, and that the proton flux decreases with energy. Auger source fits arXiv:2211.02857 reproduce spectrum and composition above the ankle with mixed composition, hard spectra, and low rigidity cutoff. The 2017 Auger astrophysical interpretation arXiv:1705.03807 reaches the same broad conclusion: relatively low maximum injection energies and hard spectra. Allard arXiv:0906.3156 emphasizes that spectrum alone cannot distinguish source-composition models; `Xmax` is needed to break the degeneracy.
For nuclei, the dominant energy-loss process is photodisintegration on CMB/EBL photons, not the Delta resonance of a free proton. The relevant Lorentz factor is `E/A`, and the relevant cross sections are giant dipole resonance / quasi-deuteron / nuclear photopion channels, not Beane's eq. 18. CRPropa/SimProp uncertainty studies explicitly treat proton, nitrogen, and iron scenarios and photodisintegration choices as major ingredients; see arXiv:1508.01824.
An optimistic ceiling if one pretended the bound applies to bound nucleons independently:
```text
response_A proportional to (b E/A)^2
b^-1 limit_A = b^-1_proton / A
```
Numerically:
| primary | A | optimistic `b^-1` | factor vs proton |
|---|---:|---:|---:|
| p | 1 | `5.77e10 GeV` | 1.000 |
| He | 4 | `1.44e10 GeV` | 0.250 |
| N | 14 | `4.12e9 GeV` | 0.071 |
| Si | 28 | `2.06e9 GeV` | 0.036 |
| Fe | 56 | `1.03e9 GeV` | 0.018 |
For a mixture, before magnetic smearing, the response coefficient would scale approximately as `sum_A f_A/A^2`, so the `b^-1` bound scales as `sqrt(sum_A f_A/A^2)`. Example: `25% p + 75% N` gives `sqrt(0.25 + 0.75/196) = 0.504`, so `2.9e10 GeV`. If the proton fraction above 40 EeV is small, the claimed bound rapidly becomes irrelevant. If the events are mostly intermediate/heavy nuclei and the physics is photodisintegration, the Beane proton-Delta conversion is simply the wrong calculation.
Damage factor: `0` to `0.5x` for plausible source-cutoff/mixed-composition cases; `1/56` in the pure-iron optimistic ceiling; more honestly, no bound without a nuclear lattice-threshold derivation.
## Sources
- Beane, Davoudi, Savage, "Constraints on the Universe as a Numerical Simulation," arXiv:1210.1847, https://arxiv.org/abs/1210.1847
- Pierre Auger Collaboration, "Measurement of the energy spectrum of cosmic rays above 10^18 eV using the Pierre Auger Observatory," arXiv:1002.1975, https://arxiv.org/abs/1002.1975
- Allard, "Ultra-high energy cosmic ray propagation in the Universe," arXiv:astro-ph/0501124, https://arxiv.org/html/astro-ph/0501124
- Allard, "Propagation of extragalactic ultra-high energy cosmic-ray nuclei," arXiv:0906.3156, https://arxiv.org/abs/0906.3156
- Berezinsky et al., "Ultra High Energy Cosmic Rays: The disappointing model," arXiv:0907.5194, https://arxiv.org/abs/0907.5194
- Pierre Auger Collaboration, "Measurements of the depth of maximum... composition implications," arXiv:1701.06812, https://arxiv.org/abs/1701.06812
- Pierre Auger Collaboration, "Astrophysical interpretation of Pierre Auger Observatory measurements...," arXiv:1705.03807, https://arxiv.org/abs/1705.03807
- Pierre Auger Collaboration, "Constraining the sources of ultra-high-energy cosmic rays across and above the ankle...," arXiv:2211.02857, https://arxiv.org/abs/2211.02857
- Alves Batista et al., "CRPropa 3.2," arXiv:2208.00107, https://arxiv.org/abs/2208.00107
- Aloisio et al., "SimProp v2r3," arXiv:1602.01239, https://arxiv.org/abs/1602.01239
- Aloisio et al., "SimProp v2r4," arXiv:1705.03729, https://arxiv.org/abs/1705.03729
- di Matteo et al., "Effects of uncertainties in simulations of extragalactic UHECR propagation, using CRPropa and SimProp," arXiv:1508.01824, https://arxiv.org/abs/1508.01824
- Sokolenko et al., "Ultra-high energy cosmic rays deflection by the Intergalactic Magnetic Field," arXiv:2101.07207, https://arxiv.org/abs/2101.07207