Comment on "Test of the Gravitational Force Law on Cosmological Scales Using the Kinematic Sunyaev--Zeldovich Effect''

Abstract: Gallardo et al. [Phys. Rev. Lett. \textbf{136}, 151002 (2026)] constrain the index of the force law $g\propto r^{-n}$ with pairwise kinematic Sunyaev--Zeldovich measurements, $n=2.1\pm0.3$, and exclude the deep-MOND value $n=1$ at $3.3σ$. The exclusion does not follow from their calculation. Their templates for $n\neq2$ [their Eq.~(8)] raise the pair-averaged Newtonian acceleration to the power $n/2$, whereas the pair average of the pointwise law [their Eq.~(30) for deep MOND] is the average of the transformed field; the step that identifies the two attributes to the three-point function a correction that is of second order and not small. In a Gaussian linear field with linearly biased tracers, Gaussian integration by parts shows that the correctly averaged template is exactly proportional to the Newtonian one, with a constant that the analysis absorbs into the mean optical depth. Gaussian realizations with peak-selected tracers find the same proportionality to within 2\% (0.5\% for the fiducial selection), and refitting the data of the Letter with the resulting template gives $χ^2=20.1$--$20.2$, against $20.1$ for the Newtonian template and $33.8$ for Eq.~(8): in every controlled case examined the data do not discriminate between $n=1$ and $n=2$, and the $3.3σ$ exclusion, with the $10σ$ forecast that rescales it, rests on the order in which the average and the power were taken. What a self-consistent MOND cosmology predicts for this observable is a separate question, which the Letter does not address.

• The paper challenges the interpretation of Gallardo et al. and claims the exclusion of the deep-MOND value $n=1$ could be due to a methodological error, not the actual physical phenomena.
• The paper's author, Alberto Melchiorri, argues that the templates used for inference in the pairwise response do not correctly reflect the pointwise force law used in the paper.
• The original analysis replaced an average of the transformed acceleration with the square root of the averaged acceleration, leading to significant implications for the inferred force law and making the proposed $3.3$ σ exclusion of $n=1$ invalid.

The comment challenges the interpretation of the pairwise kinematic Sunyaev–Zeldovich (kSZ) analysis by Gallardo et al., which reported a constraint on a phenomenological force-law index, $n=2.1\pm0.3$, and claimed a $3.3\sigma$ exclusion of the deep-MOND value $n=1$ (Gallardo et al., 15 Apr 2026). Melchiorri does not dispute the kSZ measurement or the numerical fit obtained with the published pipeline. Instead, he argues that the non-Newtonian templates used to infer the force-law exponent do not represent the pairwise response generated by the pointwise force law that they are intended to test (2610.34437).

The central distinction between averaging and nonlinear transformation

The disputed force law is written pointwise as a nonlinear transformation of the Newtonian acceleration,

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$

For a point source, this produces $g\propto r^{-n}$. In the deep-MOND case, $n=1$, the transformation becomes

$$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$

The observable, however, is not the acceleration at an isolated point. It is a tracer-pair average of the line-of-sight acceleration difference at separation $r$,

$$\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$$

According to the comment, Gallardo et al. instead construct their generalized templates by first computing the pair-averaged Newtonian acceleration and then raising its magnitude to the power $n/2$:

$3.3\sigma$0

These operations coincide for $3.3\sigma$1, because the transformation is linear. They do not generally coincide for $3.3\sigma$2. The distinction is the standard noncommutativity of nonlinear functions with ensemble averages, but it is decisive here because the radial dependence of the published $3.3\sigma$3 templates is precisely what drives the claimed statistical discrimination.

The comment reproduces the published fit using the public data products: $3.3\sigma$4 from 15 separation bins, with $3.3\sigma$5 for the Newtonian template and $3.3\sigma$6 for the published $3.3\sigma$7 template. Thus, the criticism is directed at the theoretical identification of the template, not at the numerical implementation or the underlying kSZ measurement.

The invalid treatment of the nonlinear correction

The problematic step appears in the treatment of the MOND square-root relation. The original analysis replaces an average of the transformed acceleration, schematically $3.3\sigma$8, with the square root of the averaged acceleration, $3.3\sigma$9, and assigns the neglected term to the three-point correlation function.

Melchiorri shows that this expansion is not controlled in the relevant regime. Expanding around the mean gives a leading correction proportional to the variance,

$n=1$0

not to the three-point function. Consequently, the correction is present even for a Gaussian field, for which the connected three-point function vanishes. The argument used in the Newtonian linear-theory derivation therefore cannot be transferred to a nonlinear square-root mapping.

There is also a scale-separation problem. The acceleration at a tracer is generated by the full surrounding density field, whereas the conditional mean acceleration induced by a partner at separations of $n=1$1–$n=1$2 is only a small shift. In the realizations considered by the comment, the ratio of the tracer acceleration rms to this conditional mean is approximately $n=1$3 at $n=1$4, $n=1$5 at $n=1$6, $n=1$7 at $n=1$8, and $n=1$9 at $$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$0. The expansion parameter is therefore not small at any of the relevant separations.

The implication is direct: the radial shape used to exclude $$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$1 is not obtained by consistently averaging the pointwise MOND acceleration law. It is generated by applying the nonlinear power to a conditional mean that is much smaller than the underlying local acceleration fluctuations.

Exact result for Gaussian fields with linearly biased tracers

The comment then analyzes a controlled model that closely matches the ingredients entering the published template. The density field is Gaussian and linear, the Newtonian acceleration is a linear functional of that field, and tracers are represented by linear bias weights,

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$2

The pair average is

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$3

For the nonlinear vector function

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$4

Gaussian integration by parts yields an exact relation between correlations involving $$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$5 and correlations involving the original acceleration $$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$6. Isotropy reduces the tensor derivative to a scalar coefficient,

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$7

with

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$8

The pair-averaged transformed acceleration therefore satisfies

$$\mathbf{F}(\mathbf{g})\propto \mathbf{g}|\mathbf{g}|^{n/2-1}.$$9

where $g\propto r^{-n}$0 is independent of separation. The result is exact within the stated Gaussian, linearly biased model for the considered range of $g\propto r^{-n}$1.

This is the key analytical finding. The correctly averaged nonlinear force law has the same radial shape as the Newtonian pairwise acceleration. Its dependence on $g\propto r^{-n}$2 appears only through a constant amplitude. Since the kSZ fit has a free amplitude degenerate with the mean optical depth $g\propto r^{-n}$3, that constant is not observable in the analysis. Within this model, the pairwise statistic consequently has no shape information with which to distinguish $g\propto r^{-n}$4 from $g\propto r^{-n}$5.

The result is stronger than a numerical indication that the published template is inaccurate: under the controlled assumptions, it proves that the entire $g\propto r^{-n}$6-dependent radial shape in the published construction is an artifact of taking the nonlinear power after the pair average.

Numerical validation beyond linear bias

The analytical result relies on Gaussianity and linear tracer bias, so the comment tests it numerically in more restrictive and more realistic tracer selections. Gaussian realizations of the Planck 2018 linear density field at $g\propto r^{-n}$7 are generated in periodic boxes of $g\propto r^{-n}$8 and $g\propto r^{-n}$9, with smoothing scales of $n=1$0, $n=1$1, and $n=1$2. Tracers are selected either through linear-bias weights or as density peaks above several thresholds. The MOND transformation is applied pointwise to the acceleration before pair differences and pair averaging are computed.

For linearly biased tracers, the ratio

$n=1$3

is constant to $n=1$4 peak-to-peak in the large box and $n=1$5 in the smaller box over $n=1$6–$n=1$7. The measured mean agrees with the analytical proportionality constant to $n=1$8.

Peak selection violates the assumptions of the exact derivation, yet the ratio remains nearly constant. For the fiducial selection, using a $n=1$9 box, $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$0 smoothing, and peaks above $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$1, the mean ratio agrees with the analytical constant evaluated at the tracers to $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$2. Across the tested selections and smoothing scales, deviations remain below approximately $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$3.

Figure 1: The correctly pair-averaged deep-MOND-to-Newtonian ratio remains nearly constant with separation, whereas the published $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$4 template predicts a variation by a factor of approximately $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$5.

This numerical behavior directly contrasts with the published $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$6 template. The latter predicts a radial variation by a factor of $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$7 over the same separation interval. The simulations therefore distinguish the two constructions sharply: the pointwise-transformed and then pair-averaged acceleration is nearly proportional to the Newtonian template, while the power of the pair-averaged Newtonian acceleration produces a strong and unsupported scale dependence.

Consequences for the reported significance

Melchiorri refits the kSZ data using the numerically measured ratio between the correctly averaged MOND and Newtonian accelerations. The resulting goodness of fit is effectively unchanged from the Newtonian case:

Template or selection	Minimum $$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$8
Newtonian template	$$\mathbf{F}(\mathbf{g})=\sqrt{a_0}\frac{\mathbf{g}}{\sqrt{|\mathbf{g}|}}.$$9
Linear-bias corrected template	$r$0–$r$1
Peak-selected corrected templates	$r$2–$r$3
Published $r$4 template	$r$5

Including the covariance of the simulated ratio changes the minimum $r$6 by less than $r$70.04$r$80.07$in the large-box tests.</p> <p>The reported$r$9 exclusion of $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$0 therefore derives entirely from the radial shape of the published phenomenological template. Once the pair average is performed in the order required by the pointwise force law, the deep-MOND and Newtonian templates fit the data equally well in every controlled case examined. The claimed $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$1 future forecast is consequently not a forecast for a correctly specified MOND template; it is a rescaling of a significance generated by the same disputed template construction.

The comment carefully preserves a narrower interpretation of the published fit. The constraint $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$2 can remain valid as a constraint on the exponent of the phenomenological family

$\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$3

but it cannot be interpreted as a measurement of the exponent in the underlying pointwise force law. This distinction separates an empirical fit to a chosen template family from a physical test of gravity.

Scope and unresolved physical issues

The analytical and numerical results do not constitute a complete cosmological prediction for MOND. They establish the inconsistency of the template prescription under the assumptions examined and show that the published exclusion does not follow from the stated force law. Outside Gaussian fields and the tested tracer selections, exact proportionality is not guaranteed, although the numerical departures found here are small.

More fundamentally, the pointwise relation used in the original analysis is not itself a complete MOND calculation for a generic cosmological density field. The vanishing of the solenoidal component in the relation between Newtonian and MOND fields is guaranteed by special symmetries such as planar, cylindrical, or spherical symmetry, not by considering a finite number of sources or by imposing small primordial non-Gaussianity. A realistic halo acceleration is sourced by the full three-dimensional environment.

In addition, the observed pairwise velocity is a dynamical quantity. It requires solving the modified field equation and integrating the resulting acceleration through the cosmological growth history with specified initial conditions. Applying a modified acceleration relation at one epoch does not supply that velocity prediction. The comment also notes that the alternative prescription considered by Gallardo et al. is explicitly not the prediction of a MOND theory for any $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$4.

The principal open question is therefore specific: what pairwise kSZ signal follows from a fully specified relativistic or nonrelativistic modified-gravity theory, including its cosmological field equation, nonlinear environmental response, and growth history? The paper under discussion does not answer this question, and neither does the analysis it critiques.

Conclusion

Melchiorri’s comment identifies a noncommutativity error at the foundation of the claimed force-law constraint. Gallardo et al. apply the nonlinear transformation to a pair-averaged Newtonian acceleration, whereas the stated force law requires transforming the local acceleration before performing the pair average. In a Gaussian linear field with linearly biased tracers, the correctly ordered calculation is exactly proportional to the Newtonian result; simulations with peak-selected tracers preserve this proportionality to within roughly $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$5. Refitting the data then yields $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$6 for both the Newtonian and deep-MOND constructions, rather than $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$7 for the published $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$8 template. The reported $\langle X_n\rangle_{12}(r) = \left\langle \left[\mathbf{F}(\mathbf{g}_2)-\mathbf{F}(\mathbf{g}_1)\right] \cdot\mathbf{e}_{12} \right\rangle_{12}.$9 exclusion, and its projected $n/2$0 extension, therefore constrain a phenomenological template rather than the index of the underlying force law.