A regime-dependent constitutive exponent in SPARC rotation curves: phenomenological evidence beyond a universal MOND interpolation

1 Introduction

Galaxy rotation curves remain the cleanest astrophysical probe of the gravitational regime in which Newtonian acceleration falls below ≈10⁻¹⁰ m s⁻². The flatness of these curves at large radii, and the tight baryonic Tully–Fisher relation (BTFR; McGaugh et al., 2000; Lelli et al., 2016), are reproduced by two broad families: cold dark matter haloes (Rubin et al., 1980; Navarro et al., 1996) and modified Newtonian dynamics (MOND; Milgrom, 1983; Famaey and McGaugh, 2012). Both predict the radial acceleration relation (RAR; McGaugh et al., 2016) to comparable accuracy.

A central assumption of MOND-class phenomenology, implicit since Milgrom (1983) and explicit in all subsequent parametric ansätze (Famaey and McGaugh, 2012), is that a single interpolation function relates observed to Newtonian accelerations across the full sample of rotation-supported galaxies. Recent Spitzer Photometry and Accurate Rotation Curves (SPARC) database-based work has tested MOND in several ways: differentiating modified gravity from modified inertia through the global deep-MOND parameter Q (Petersen and Lelli, 2020); detecting the external field effect and linking it to large-scale structure (Chae et al., 2020; Chae et al., 2021); contrasting modified gravity, modified inertia, and dark matter predictions through inner-vs-outer rotation-curve behavior (Chae, 2022); and proposing alternative empirical correction terms (Flynn and Cannaliato, 2025).

What has not been carried out, to our knowledge, is a direct test of the universality assumption at the per-galaxy level within a parametric family that contains the canonical deep-MOND scaling as the p = 1 limit. The present work performs exactly that test. We adopt the one-parameter constitutive family β(g) = (a0/g)^p, a phenomenological ansatz in which the deep-MOND limit p = 1 is one specific value among a continuum, and fit p galaxy-by-galaxy on the 175-galaxy SPARC sample. The universal-p hypothesis is then confronted with the resulting population statistics through six complementary statistical diagnostics (BTFR slope, calibrated synthetic comparison, regime correlation, out-of-sample cross-validation, alternative-binning audit, and Bayesian information criterion (BIC)). We make no claim about the underlying physical mechanism; we aim to determine, empirically and at the population level, whether SPARC rotation curves are consistent with a universal p or require structured variation.

The remainder of this article is organized as follows. Section 2 presents the phenomenological framework: the constitutive family β(g) = (a0/g)^p, the closed-form deep-regime relations it generates, and the regularized interpolation used across the full SPARC dynamic range. Section 3 describes the SPARC data and the fitting procedure. Section 4 reports tests of two universal-p benchmarks (p = 3/4 and p = 1) against the baryonic Tully–Fisher relation, the radial acceleration relation, and galaxy-by-galaxy χ². Section 5 fits p galaxy-by-galaxy and characterizes the distribution of best-fit exponents, including a calibrated comparison against single-p synthetic catalogs. Section 6 partitions the sample by the outer acceleration g_outer, reports the monotonic progression of median p across regime bins, validates the regime structure out-of-sample against seven alternative-binning schemes, and compares all models on an information-criterion basis. Section 7 discusses interpretations, limitations, and tests that would strengthen or falsify the result. Section 8 concludes. Supplementary Material provides the algebraic derivation of the deep-regime law (Supplementary Material S1), inverse calibration tables (Supplementary Material S2), the complete galaxy-by-galaxy fit table (Supplementary Table S11; 175 galaxies), reproducibility code listings (Supplementary Material S6, S8.3), and robustness tests of the BTFR slope to regression method and error model (Supplementary Material S9; revised-estimator values in Supplementary Table S10), of the BIC comparison to effective sample size and nuisance-parameter structure (Berger et al., 1998; Supplementary Material S8.4, S10), and of the population statistics to outlier-rejection cuts (Supplementary Material S11, Supplementary Table S7). Multi-split cross-validation (Supplementary Table S9), the a₀-convention sensitivity of the Section 4.3 count (Supplementary Table S12), the complete per-galaxy fit table of the consolidated pipeline (Supplementary Table S11), the full code listings (Supplementary Material S13), and the data-integrity checksums (Supplementary Material S14) are new in this revision.

2 Theoretical framework

2.1 The constitutive family β(g) = (a₀/g)^p

We consider a one-parameter family of interpolating laws that connects the Newtonian regime g = g_N to a deep-regime power law characterized by an exponent p > 0. The family can be written in a modified Poisson form in which a dimensionless constitutive function β(g) modulates the gravitational coupling between the observed acceleration g = |∇Φ| and the Newtonian baryonic source ρ_b. With the convention β = 1 corresponding to unmodulated Newtonian propagation, and β > 1 enhancing the coupling, the form is given in Equation 1

∇ · β − 1 g ∇ Φ = 4 π G ρ _ b . (1)

This form recovers Newton’s limit β → 1 and modifies the coupling in regimes where the medium has β > 1. The constitutive function admits the power-law expansion of Equation 2, in the dimensionless variable g/a₀,with β large in the deep regime (g ≪ a₀) and β → 1 as g → a₀. The exponent p is a free phenomenological parameter of the family. Two universal-p choices anchor the analysis below: p = 1, which reproduces the standard deep-MOND scaling g = √(a₀ g_N); and p = 3/4, adopted here as an intermediate benchmark to be tested rather than as a derived prediction. The benchmark value p = 3/4 was motivated by prior unpublished theoretical work of one of us (M.A.S.) on regime-dependent gravitational coupling (Simonssini and Ribeiro, in preparation). We adopt it purely as a phenomenological benchmark to be tested empirically against SPARC, and the analysis implies no commitment to the underlying framework. Integration of the modified Poisson form in spherical stationary symmetry yields β − 1 g · g = g _ N , that is , g / a 0 ^ p · g = g _ N, which solves to

β g = a 0 / g ^ p , (2)

g = a 0 ^ p / 1 + p · g _ N ^ 1 / 1 + p . (3)

Equation 3 generalizes the standard deep-MOND scaling (recovered for p = 1, g = √(a₀ g_N)) and yields a baryonic Tully–Fisher slope of 2 (1 + p). For the benchmark p = ¾, this is 7/2 = 3.5, with deep-regime law g = a₀^(3/7) · g_N^(4/7); for the MOND value p = 1, the slope is 4.

The baryonic Tully–Fisher relation (BTFR) is an empirical scaling between the baryonic mass of a rotation-supported galaxy and its characteristic flat rotation velocity, M_bar ∝ V_flat^s (McGaugh et al., 2000; Lelli et al., 2019). Within the constitutive family tested here, the deep-regime law g^(1 + p) = g_N a0^p implies s = 2 + 2p: canonical MOND (p = 1) predicts the classical slope of 4, while p = 3/4 predicts 7/2. The BTFR slope is therefore a direct, population-level discriminant among members of the family, independent of the galaxy-by-galaxy fits of Section 5.

For numerical computation across the full SPARC dynamic range, we use the regularized interpolationwhich recovers the Newtonian limit ν → 1 for y ≫ 1 and the deep-regime limit ν → y^(−p/(1 + p)) for y ≪ 1. For p = ¾, the deep-regime exponent is −3/7; for p = 1, it is −1/2 (standard deep-MOND).

ν _ p y = 1 + y ^ − p ^ 1 / 1 + p , y ≡ g _ N / a 0 , (4)

2.2 Phenomenological hypothesis: is p universal?

The universal-p analysis of Section 2.1 treats p as a single number characterizing the whole galaxy population. This is the standard assumption underlying both the canonical MOND interpolation and its variants: a single function ν(g_N/a0), with a single a0, describes all galaxies. The strength of this assumption is empirical, not theoretical: it has held remarkably well at the level of the radial acceleration relation (McGaugh et al., 2016) but is rarely tested directly at the level of per-galaxy fits within a parametric family.

We therefore relax the universality assumption and ask, as a phenomenological question, whether the SPARC sample is consistent with a single best-fitting p_universal or whether per-galaxy fits prefer a distribution of effective exponents. If the distribution is broader than measurement noise allows under a single-p hypothesis, and if its variation is organized by an identifiable physical regime variable (rather than being random), then the universal-p assumption is observationally insufficient, and the data carry information about a regime-dependent coupling. Conversely, if per-galaxy fits cluster tightly around a single p, the universal-p picture is confirmed. We test both possibilities in Section 5 and Section 6. We treat p = 3/4 as a benchmark central value for comparison with p = 1 (MOND), but make no claim that p = 3/4 is fundamentally privileged; the test is whether SPARC distinguishes between universal-p hypotheses at any value.

3 Data and methods

3.1 The SPARC sample

We use the SPARC database (Lelli et al., 2016): 175 disc galaxies with Spitzer 3.6 μm photometry and HI/Hα rotation curves. For each rotation curve point, the catalog provides the radius R, the observed circular velocity V_obs with uncertainty eV, and the baryonic components V_gas, V_disc and V_bul (the latter two tabulated by the SPARC catalog at unit mass-to-light ratio, precisely so that they can be rescaled by the fitted Υ_disc and fixed Υ_bul of Equation 5; no stage of the analysis adopts unit mass-to-light ratios). The baryonic velocity-squared profile iswith the sign convention preserved for V_gas. The Newtonian baryonic acceleration is g_N(R) = V2_bar(R)/R.

V 2 _ bar R = V 2 _ gas + Υ _ disc · sgn V _ disc · V 2 _ disc + Υ _ bul · sgn V _ bul · V 2 _ bul (5)

3.2 Fitting procedure

Predicted circular velocity is given by Equation 6:with ν_p(y) given by Equation 4, and reference cases obtained by setting p = 3/4 (denoted ν_(3/4) below) and p = 1, the member of the present constitutive family that recovers the standard deep-MOND scaling; we use it as the MOND-class benchmark within this family, not as an exhaustive test of all MOND interpolation functions. For each universal-p choice, the constant a0 is fitted globally to SPARC by minimizing the total velocity-space χ² with Υ_disc refitted per galaxy at each trial a0, via a one-dimensional bounded search, over all rotation curve points with positive baryonic g_N (3,389 points). The stellar disc mass-to-light ratio Υ_disc is fitted per galaxy with a log-normal prior centered on log Υ_disc = log (0.5), σ = 0.1 dex (Lelli et al., 2016); the bulge Υ_bul is fixed at 0.7.

V _ pred R = ν g _ N / a 0 · g _ N R · R (6)

For the galaxy-by-galaxy regime analysis (Section 5 and Section 6), we additionally fit, for each galaxy, the constitutive family β(g) = (a0/g)^p with both Υ_disc and p free, using the same Υ prior. The resulting p_best is the effective constitutive exponent of that galaxy under the family. For the regime analysis, we use a fixed central value a0 = 1.5 × 10⁻¹⁰ m s⁻², intermediate between the universal-p best-fitting values for p = 3/4 and p = 1. Results are insensitive to a0 choices within this range. Because a per-galaxy a0 and a per-galaxy p are partially degenerate within the constitutive family, both rescale the location and sharpness of the Newtonian-to-deep transition. Fixing a0 does not suppress physically meaningful freedom; it is reparametrized into p. We verify in Supplementary Material S4.5 that varying the central a0 across the full plausible range (1.12–1.79 × 10⁻¹⁰ m s⁻²) shifts the median p but leaves the p–g_outer correlation invariant (ρ ≈ +0.47–+0.49). The robustness of the population statistics to the Υ_disc prior width (σ = 0.05–0.20 dex) is documented in Supplementary Material S4.4.

The outer acceleration of each galaxy, used as the regime variable, is defined as g_outer = V2_outer/R_outer, with V_outer being the mean of the last five rotation curve points (or last third of the curve, whichever is larger) and R_outer being the outermost measured radius.

We emphasize the separable structure of the optimization. Conditional on the single global parameter a0, the per-galaxy objectives are complementary: the ∼160-dimensional problem decomposes exactly into independent one-dimensional (universal-p case) or two-dimensional (regime analysis) per-galaxy optimizations, nested inside a one-dimensional search over a0. The global optimum is therefore found to within the tolerance of the one-dimensional searches; no high-dimensional joint minimization is performed or required. This structure has been verified independently: re-implementing the procedure described in this section and public SPARC data reproduces the population statistics of Section 5 (median, dispersion, 16%–84% interval, and the p–g_outer correlation).

3.3 The BTFR sample

For the BTFR, we adopt the reproducible catalog criterion V_flat >0 (Lelli et al., 2016), yielding 135 galaxies; the catalog V_flat is itself defined only where the outer rotation curve satisfies the SPARC flatness algorithm, so the selection embeds the flat-curve requirement. V_flat is the outer mean; M_bar is the baryonic mass integrated to the last measured radius using the per-galaxy fitted Υ_disc.

3.4 Out-of-sample validation protocol

To assess whether the regime structure has predictive content, we apply a pre-specified out-of-sample protocol. The sample is randomly partitioned 50/50 into a train set and a test set. In the train set, we compute the median p_best across three regime bins of x = log10 (g_outer/a0). In the revised main analysis, these are equal-population terciles of x; the original fixed-edge bins (x < −1.0; −1.0 ≤ x < −0.5; −0.5 ≤ x < 0.0) are retained as a robustness check and yield the same qualitative result (Section 6.2). Two galaxies have x ≥ 0 (NGC5005 and UGC11914); under the original fixed-edge bins, they fell outside all bins and were dropped, while under the tercile convention, they belong to the uppermost tercile and are included. Their exclusion does not alter any conclusion. For each test galaxy, we look up the bin median from the training set, fix p to that value, and re-fit only Υ_disc (no per-galaxy p fitting on the test set). The procedure is repeated 100 times with independent random splits.

We compare four candidate models on the test set: (a) p = 1 (universal, standard MOND); (b) p = 3/4 (universal); (c) the best single-p of the training set (median of train p_best); (d) the regime rule (bin median from training, applied by regime of the test galaxy). The test set figure of merit is the median reduced χ² across test galaxies.

We use 50–50 train–test splits because, at this sample size, test set variance dominates the validation error; equal halves maximize the precision of the test set comparison while keeping the training bins populated. Robustness under 70–30 and 60–40 splits is reported in Section 6.2.

3.5 Auditing the regime variable

To test whether the acceleration-regime variable is privileged among physical proxies, we repeat the out-of-sample protocol under alternative binnings: a random-balanced control (same bin sizes as the real binning, but galaxies assigned randomly), and tertile-binnings by log M_bar, log R_outer, log V_outer, f_gas, and log Σ_bar. Each scheme uses 50 random splits and is compared with the real g_outer binning using the Mann–Whitney U test.

4 Universal-p tests: BTFR, RAR, and galaxy-by-galaxy comparison

4.1 Baryonic Tully–Fisher relation

Figure 1 shows the BTFR for the 135 flat-curve galaxies with the two universal-p benchmarks and the headline fit. We adopt an error-in-variables estimator (orthogonal distance regression, ODR; cf. Isobe et al., 1990; Hogg et al., 2010) as the headline, because ordinary least squares (OLS) is not consistent under heteroskedastic errors in both variables:

log ⁡ M bar = 3.47 ± 0.10 · log ⁡ V flat + c ± σ c ODR . (7)

OLS yields 3.28 ± 0.10 (R² = 0.899), and a maximum-likelihood fit with intrinsic scatter (0.21 dex) yields 3.34 ± 0.12 on the same sample: a method systematic of ≈0.19. Across the plausible range of error models (σ_logM = 0.10–0.20 dex; velocity-error floors of 1–4 km s⁻¹), the ODR slope spans 3.40–3.61: consistently within ≈1σ of the p = 3/4 prediction (7/2) and incompatible with the canonical MOND slope of 4 at ≥4σ (statistical; minimum 4.1σ) under every error model tested.

Published BTFR slopes for SPARC depend on estimator, velocity definition, and mass proxy: Lelli et al. (2019) report 3.85 ± 0.09 from an orthogonal maximum-likelihood fit with photometric masses at fixed Υ⋆ = 0.50, while vertical and least-trimmed-squares estimators yield systematically shallower slopes (Lelli et al., 2019; appendix; Schombert et al., 2020), and Lelli et al. (2019) note that systematic uncertainties may drive the slope between 3.5 and 4.0. Our estimate, with dynamically integrated masses and a per-galaxy-fitted Υ_disc, is consistent with this pattern. The comparison relevant to this work is internal: with masses, sample, and estimator held fixed, the observed slope is incompatible with s = 4 and consistent with s = 7/2.

4.2 Radial acceleration relation

On the 3389 SPARC data points with positive baryonic g_N, the best-fitting universal a0 values, which are obtained by the velocity-space procedure of Section 3.2, are

a 0 p = 3 / 4 = 1.79 × 10 − 10 m s − 2 , a 0 MOND , p = 1 = 1.12 × 10 − 10 m s − 2 .

Dispersion of residuals: 0.161 dex (MOND) and 0.162 dex (p = 3/4). The two universal laws reproduce the RAR with statistically indistinguishable quality (Figure 2).

4.3 Galaxy-by-galaxy comparison of p = 3/4 vs. MOND

With Υ_disc per galaxy and a0 fixed at the common regime value a0 = 1.5 × 10⁻¹⁰ m s⁻² for both models, p = 3/4 yields a smaller χ² than MOND in 106 of 175 galaxies (60.6%), binomial p = 0.006. Under Monte Carlo marginalization over distance (σ_D/D = 0.15) and inclination (σ_i = 3°) uncertainties, with the catalog inclinations of each galaxy, the marginalized median preference is 99/175 (range 87–113); 53% of realizations are individually significant at p < 0.05. The count exceeds the coin-flip expectation in 99.5% of realizations (Figure 3). The direct preference of p = 3/4 over MOND on SPARC is statistically significant, and its direction survives the dominant systematic uncertainties of the catalog, with per-realization significance attenuated as expected under a 15% distance scatter.

The preference count is mildly sensitive to the a0 convention: at each model’s own global best-fitting a0 (Section 4.2), the count is 103/175 (binomial p = 0.023), which is still significant. Per-galaxy χ² values for both models are provided in Supplementary Table S11.

4.4 Dense universal-p profile

To sample the universal-p family densely, we profile the total χ² over p ∈ [0.1, 2.0] in steps of 0.05, refitting a0 globally and Υ_disc per galaxy at each p. The profile has a single minimum at p = 0.85 (Figure 4), adjacent to the population median of the galaxy-by-galaxy fits (0.82; Section 5.1), with p = 3/4 preferred over p = 1 (Δχ² = 576) and both lying well above the minimum. No universal exponent, including the profile minimum, reproduces the population distribution (Section 5.2), which is the central point: the dense profile constrains where the family fails least, not where it succeeds.

5 Population distribution of the constitutive exponent

5.1 Galaxy-by-galaxy fits of the constitutive family

The universal-p tests of Section 4 fix p across all galaxies. To test whether SPARC requires a regime-dependent exponent (Section 2.2), we fit the constitutive family β(g) = (a0/g)^p galaxy-by-galaxy with both Υ_disc and p free. Of 175 SPARC galaxies, 160 yield sensible fits (0.05 < p_best <5, χ²_red <50). The robustness of the population statistics under more restrictive outlier cuts (0.3 < p < 2 and 0.4 < p < 1.6) is documented in Supplementary Table S7: the median, the intrinsic dispersion exceeding the synthetic single-p dispersion, and the Spearman correlation with the regime variable are all robust to outlier filtering.

Figure 5a shows the resulting distribution of p_best. The median is 0.82 with intrinsic dispersion σ_p = 0.53; the 16%–84% interval is [0.55, 1.35]. The distribution is broad and asymmetric (skewness +1.6), with 11% of galaxies in a tail at p > 1.5.

5.2 Calibration against single-p synthetic catalogs

Is this dispersion noise, or genuine? We generate synthetic SPARC catalogs with a single true exponent p_true ∈ [0.55, 1.05] in steps of 0.05, preserving the SPARC structure of R, V_bar components, and measurement errors; each catalog is then refitted galaxy-by-galaxy under the constitutive family. The calibration is unbiased: the recovered median p is linear with slope 1.00 ± 0.005 in p_true. The p_true that reproduces the SPARC median is linear in the input (slope 1.00 ± 0.005); the consolidated-pipeline median of 0.82 (Supplementary Material S13) corresponds to p_true ≈0.82.

However, no single p_true reproduces the SPARC distribution shape. The recovered synthetic dispersion is σ_p ≈ 0.15 across all p_true tested, which is three times narrower than the SPARC σ_p = 0.53. Kolmogorov–Smirnov tests reject every single p_true against the real SPARC distribution at p < 10⁻¹¹. The dispersion of effective exponents in SPARC is intrinsic, not measurement noise: galaxies vary in their preferred p.

As a corollary, the population median of SPARC falls between the single-p recoveries and is incompatible with the MOND value. The calibration is unbiased (recovered median linear in p_true with slope 1.00; Supplementary Material S2, Supplementary Table S1): synthetic catalogs with p_true = 1 (MOND) recover a population median of 1.010 ± 0.005, and those with p_true = 0.75 recover 0.755 ± 0.005. The observed SPARC median, 0.82, therefore sits slightly above the p = 0.75 recovery and well below the p = 1 recovery, bracketed between the two single-p expectations and far from the MOND value. The decisive rejection of any single universal p comes not from the median but from the distribution shape and width: every single-p catalog is rejected against the real SPARC distribution by the Kolmogorov–Smirnov test (p < 10⁻¹¹) and recovers a dispersion (σ_p ≈ 0.15) three times narrower than observed (σ_p = 0.53). This is consistent with the empirical hypothesis that the population averages a distribution of regime-dependent exponents rather than reflecting a single universal p.

5.3 Correlation with the outer acceleration

The intrinsic dispersion correlates with the local acceleration regime (Figure 5b). The Spearman correlation of p_best with log10 g_outer is ρ = +0.49 with p = 4 × 10⁻¹¹, robust under 20 split-half realizations (all 20 give significant same-sign correlation in both halves) and 1,000 random permutations of p (0 of 1,000 reach |ρ| ≥ 0.48). Galaxies with low outer acceleration prefer p ≈ 0.5–0.7 (below the MOND value); galaxies with high outer acceleration prefer p ≈ 1.0–1.5 (MOND-like or above). Of the alternative galaxy properties tested, including log M_bar, log R_outer, log V_outer, f_gas, log Σ_bar, and log Σ_disc, only log V_outer (ρ = +0.31, p = 6 × 10⁻⁵) and f_gas (ρ = −0.17, p = 0.03) reach statistical significance, both weaker than g_outer. The correlation is specific to the acceleration-regime variable. This correlation is not an artifact of the fixed a0: repeating the galaxy-by-galaxy fits at a0 = 1.12 × 10⁻¹⁰ and 1.79 × 10⁻¹⁰ m s⁻² changes the population median p (to 0.94 and 0.71, respectively) but leaves the Spearman correlation with log g_outer essentially unchanged at ρ = +0.47–+0.49 (p ∼ 10⁻¹⁰; Supplementary Material S4.5).

6 Regime structure and predictive validation

6.1 Monotonic progression of the median exponent across regime bins

Partitioning the sample by the regime variable x = log10 (g_outer/a0) into three physically-motivated bins, we find:

Very deep (lowest tercile of x; median x̃ = −1.05, that is, g_outer ≈1.3 × 10⁻¹¹ m s⁻²): 54 galaxies, median p̃ = 0.63; Deep (middle tercile; x̃ = −0.77, g_outer ≈2.5 × 10⁻¹¹): 53 galaxies, p̃ = 0.83; Intermediate (upper tercile; x̃ = −0.53, g_outer ≈4.4 × 10⁻¹¹): 53 galaxies, p̃ = 1.11. These (x̃, p̃) pairs constitute the empirical regime-response curve that any successor framework must reproduce. The within-tercile dispersion is itself structured, growing across regimes (σ_p = 0.26, 0.36, 0.68 from very deep to intermediate); it is a second moment of the regime dependence, tabulated per galaxy in Supplementary Table S11, that constrains candidate mechanisms as strongly as the medians.

Kruskal–Wallis rejects equality of the three medians at H = 33.1, p = 6 × 10⁻⁸. Mann–Whitney tests for adjacent terciles: Very deep → Deep, p = 4.5 × 10⁻⁴; Deep → Intermediate, p = 5.1 × 10⁻³. Spearman correlation between tercile index and p: ρ = +0.47, p = 4 × 10⁻¹⁰. The progression is unambiguous (Figure 6).

Robustness check using equal-population quartiles (Q1–Q4; original submission pipeline) yields 0.62, 0.77, 0.91, 1.13, which is the same monotonic structure under finer binning.

6.2 Out-of-sample validation

To test whether the regime structure carries predictive content beyond the data in which it was identified, we apply the protocol in Section 3.4: train regime medians on half the sample and predict the other half. Across 100 random splits: ⟨χ²_red⟩_test (50–50, tercile bins): MOND (p = 1) = 3.92 ± 0.46; single-p global = 3.70 ± 0.56; p = 3/4: 3.54 ± 0.47; regime rule = 3.33 ± 0.40.

The regime rule has a smaller test set median χ²_red than MOND in 97 of 100 splits (99/100 under the fixed-edge binning of the original submission), than single-p global in 84 of 100, and than p = 3/4 in 77 of 100. The result is robust to the splitting ratio: under 70–30 splits, the regime rule beats MOND in 91/100 cases, and under 60–40 splits, the regime rule beats MOND in 93/100 cases. The training set tercile medians are stable across splits (0.63, 0.83, 1.11), with separations between bin medians an order of magnitude larger than train-to-train fluctuation.

6.3 Audit against alternative binnings

The previous result raises the question: does the gain come from the physics of acceleration regimes, or from the mere flexibility of three exponents in any binning? To answer, we repeat the out-of-sample protocol with seven alternative-binning schemes, including a random-balanced control and tertile bins by physical variables (Table 1). Each scheme uses 50 random splits.

TABLE 1. Original fixed-edge audit, retained as a robustness check (values from the original submission pipeline). The revised main validation and its absolute χ²_red values are reported in Section 6.2 and Supplementary Table S9.

Table 1: Out-of-sample test set advantages of each regime-rule binning, averaged over 50 random splits. Δχ²_red is the difference relative to each universal baseline. The acceleration-regime variable (rows A and C) yields the largest advantage on all three baselines and is the only variable that beats p = 3/4 on average. Mann–Whitney p in the last column compares each scheme to row A; the 2 g_outer schemes are statistically indistinguishable from each other (p = 0.41), and every other scheme is rejected at p < 10⁻⁴.

The acceleration-regime variable is qualitatively privileged (Figure 7). Both edge-based and tertile-based binning by g_outer yield the largest advantage across all baselines. Variables that correlate weakly with g_outer (f_gas, Σ_bar, V_outer) yield smaller but still positive advantages. Variables that do not correlate (M_bar, R_outer) drop close to the random-balanced control. Random-balanced binning still beats MOND (Δ = +0.38) because three free exponents have more flexibility than one, but it is rejected at p = 4 × 10⁻⁸ against the real g_outer binning.

6.4 Information-criterion comparison

The cross-validation analysis of Sections 6.2, 6.3 shows that the regime rule generalizes out of sample. It does not, on its own, address whether the gain over universal-p models could result purely from the additional parameters (two more p values beyond the single universal p). To address this directly, we report the Bayesian information criterion (Schwarz, 1978), which explicitly penalizes parameter count:where χ² is the total residual summed over all rotation curve points, k is the number of free parameters in the model, and n is the total number of data points. A smaller BIC indicates a more strongly supported model; ΔBIC = BIC_i − BIC_ref between two models is interpreted on the scale of Kass and Raftery (1995): 2 < |ΔBIC| < 6 is positive evidence, 6 < |ΔBIC| < 10 is strong, and |ΔBIC| > 10 is very strong.

BIC = χ 2 + k · ln n , (8)

We compute BIC across the SPARC sample (160 galaxies with sensible fits; n = 3,215 rotation curve points) for the candidate models, in the single objective convention of the deposited listings (Supplementary Material S13; the per-galaxy Υ-prior term included throughout), with parameter counts including the disc mass-to-light ratio Υ_disc per galaxy and the global parameters of each model. The regime rule is preferred over MOND at ΔBIC = −12,112 and over the universal p = 3/4 at ΔBIC = −6,628.

Two features of the comparison address the hierarchical structure of the problem directly. First, the per-galaxy mass-to-light ratios are nuisance parameters common to every model compared; removing all of them symmetrically from every model leaves every ΔBIC invariant to rounding, because the comparison is differential. Second, the critical control compares the regime rule against the single-p free model: two models with identical nuisance structure differing by exactly two global parameters. The BIC penalty for those parameters is 2·ln(n) ≈ 16.2, against a χ² reduction of 3,636 (ΔBIC = −3,620): more than 99% of the gain is data-driven, ruling out parameter inflation as the source of the preference. In addition, the per-galaxy distribution of evidence (Figure 8a,b) shows the preference is distributed across the population, not an artifact of pooling: with p fixed at the training-bin median (equal per-galaxy parameter counts), the regime rule is favored over MOND in 102 of 160 galaxies (median per-galaxy Δχ² = −3.5; 64 galaxies exceed Δχ² < −10 against 33 in the opposite direction, and the per-galaxy contributions of Supplementary Table S11 sum to the pooled value by construction: all comparisons in this section share the single objective convention of Supplementary Material S13, which includes the Υ-prior term). The per-galaxy free exponent, which is the upper envelope of achievable per-galaxy evidence, penalized by ln(n) for its extra parameter, yields a median ΔBIC of −18 over MOND, with 104 of 175 galaxies beyond the |ΔBIC| > 10 “very strong” threshold of Kass and Raftery (1995); both quantities are tabulated per galaxy in Supplementary Table S11.

6.5 Kinematic robustness of the regime variable

A natural concern is that the p–g_outer correlation could be kinematic in origin: g_outer is constructed from points on the outer rotation curve, whose flatness varies across the sample. We quantify this directly. The outer band has mean radius 15.7 ± 16.1 kpc (median 9.6 kpc) under the last-five-points convention and 14.5 ± 13.6 kpc (median 9.3 kpc) under the last-third convention; by the logarithmic-slope criterion |d ln V/d ln R| ≤ 0.05, 17.7% of outer regions are strictly flat, 61% are gently rising (median outer slope +0.10) and 21% are declining, with a median outer velocity spread ΔV/V of 6.2%. V_outer is operationally equivalent to the catalog V_flat, where the latter is defined (Pearson r = 0.999, N = 135). The construction of g_outer is robust to the band convention (r = 0.97 between conventions; median ratio 1.008). The kinematic covariates below are computed on the same band used for g_outer (last max (5, N/3) points). The confounder exists: log g_outer correlates with outer flatness (Spearman ρ = −0.34, p = 1.4 × 10⁻⁵). The exponent, however, does not: p versus ΔV/V gives ρ = −0.14 (p = 0.09), and p versus the outer slope gives ρ = −0.08 (p = 0.31). Decisively, the partial correlation of p with log g_outer controlling for both kinematic covariates is ρ = +0.50 (p = 2 × 10⁻¹¹), which is unchanged from the raw correlation, while the reverse control (p versus flatness at fixed g_outer) vanishes (ρ = +0.04, p = 0.66). The regime structure tracks the acceleration, not the shape or construction of the outer curve. Representative rotation curve fits across the three acceleration regimes are shown in Figure 9.

7 Discussion

7.1 Interpretation of the regime-dependent exponent

The structured failure of universality documented here speaks most directly to the family of modified-gravity proposals to which MOND belongs. Within that family, our results disfavor any theory committed to a single universal interpolation function, while remaining compatible with, and requiring, proposals in which the effective constitutive exponent is non-universal or regime-dependent. Frameworks of this type exist in the literature, ranging from extended MOND variants with environment- or acceleration-dependent interpolation to proposals that derive the constitutive exponent from an underlying mechanism, including emergent-gravity scenarios (Verlinde, 2017); a first-principles framework aimed at deriving the benchmark value p = 3/4, and at confronting the regime progression reported here, is developed separately (Simonssini and Ribeiro, in prep.). On the modified-gravity reading, the regime structure reported here is the primary empirical target any such theory must reproduce.

The same phenomenology admits a standard ΛCDM interpretation in which no modification of gravity is required: the diversity of effective exponents then reflects the diversity of dark matter halo profiles and their coupling to baryonic feedback (e.g., Di Cintio et al., 2014; Katz et al., 2017). On this reading, p is a summary statistic of halo structure rather than of a constitutive law, and the correlation with outer acceleration traces the well-known dependence of halo response on surface density. Our analysis does not discriminate between these two readings; it constrains the shape any successful account must take.

Whatever the ontology, three features require explanation: (i) the rejection of a universal exponent at high significance; (ii) the monotonic progression of the median exponent across acceleration regimes; and (iii) the survival of the p–g_outer correlation under kinematic controls (Section 6.5). These are the portable, theory-independent results of this work.

7.2 Limitations

Several caveats apply. (i) The mass-to-light ratio prior σ(log Υ_disc) = 0.1 dex follows standard SPARC practice but is not derived from stellar-population synthesis; tighter or looser priors will affect per-galaxy χ² differences. (ii) The marginalization assumes uniform σ_D/D = 0.15 and σ_i = 3°, typical of SPARC, but individual distances range in precision from 3% (Cepheid) to 30% (Hubble-flow). (iii) The regime bins are coarse: three bins on a continuous regime variable. The choice of edges does not drive the result; tertile-based binning provides the same progression, but a finer-resolution analysis would benefit from a larger sample. (iv) The galaxy-by-galaxy fitting fixes a0 at a central value (1.5 × 10⁻¹⁰ m s⁻²) for the constitutive family; this is a fitting convention, not a derived quantity. (v) The choice of g_outer as proxy for the local acceleration regime is physically natural but not unique; testing finer regime variables (e.g., gradient of g_N, environmental density) is left for future work.

7.3 Comparison with previous tests of modified gravity

MOND’s BTFR slope prediction of 4 has long been considered one of its tightest (McGaugh, 2012; Famaey and McGaugh, 2012). Observed SPARC slopes have ranged from 3.5 to 4.0 depending on selection, Υ treatment, and regression method (Lelli et al., 2016; Lelli et al., 2019). Our headline value 3.47 ± 0.10 comes from orthogonal distance regression, adopted to account for errors in both variables; OLS and maximum-likelihood alternatives on the same sample are reported in Section 4.1 and Supplementary Table S10. The dependence of the slope on the regression method (OLS, ODR, bisector) and on the error model is documented in Supplementary Table S10 (the original method-and-cut grid, including bisector regression, is retained in Part I of the Supplementary Material). Under every method-cut combination tested, the observed slope sits between the p = 3/4 prediction (3.5) and the canonical MOND prediction (4); the qualitative conclusion that p = 3/4 fits the BTFR more closely than p = 1 is robust to these methodological choices.

The galaxy-by-galaxy preference of 106/175 we report is, on its own, modest: a 61% preference is statistically significant under a binomial test but not overwhelming as an individual measure. The case for a regime-dependent exponent rests instead on the joint pattern across multiple complementary diagnostics: universal-p preference at the BTFR level, calibrated synthetic comparisons rejecting single-p hypotheses, monotonic regime progression, out-of-sample validation, comparative audit against alternative binnings, and information-criterion comparison. The cross-validation and BIC results are individually striking and jointly blind the analysis against the most natural critique: BIC penalizes in-sample complexity, cross-validation penalizes out-of-sample overfit, and both criteria independently select the regime rule against universal-p alternatives. No single test is conclusive; the pattern across tests is what we report.

7.4 Tests that would strengthen or falsify the result

Within the constitutive family itself, Equation 4 supplies a clean forthcoming discriminant: the deep-regime logarithmic slope of the radial acceleration relation tends to 1/(1 + p), yielding 1/2 for canonical MOND (p = 1) and 4/7 ≈ 0.571 for the p = 3/4 benchmark. Current SPARC data do not reach deep enough to separate these; weak-lensing and ultra-diffuse-galaxy measurements pushing below g_N ≈ 10⁻¹² m s⁻² would decide between family members independently of rotation curve fitting.

Three classes of test could decisively strengthen or falsify the present results. (i) Ultra-diffuse galaxies, with g_outer well below the SPARC range, would extend the regime variable beyond the deepest bin probed here. Under the extrapolation of the regime trend reported in Section 6.1, such systems would be expected to display p values systematically below 0.5; the available kinematic samples (van Dokkum et al., 2018; van Dokkum et al., 2019) are small but extensible. (ii) Galaxy clusters, where MOND is known to require residual dark matter, would test whether the regime-dependent constitutive structure extrapolates beyond galactic scales. (iii) Wide binary stars in the low-acceleration regime (Chae, 2023) probe deep-field gravity with stellar data and would discriminate between universal-p and regime-dependent scenarios under different priors. A fourth class of test is provided by galaxies discovered after the present analysis is published: the regime rule of Section 6 is now specified and can be applied to such samples without further adjustment, providing genuine predictive validation.

7.5 Failure mode of the universal-p hypothesis

The central result is not that one fixed exponent replaces MOND. Rather, the failure mode of the universal-p hypothesis is structured. When galaxies are fitted individually under the constitutive family, the recovered exponent does not scatter randomly around a universal value: it drifts monotonically with the outer-acceleration regime, correlates specifically with g_outer, and produces a population dispersion that no single-p generative model reproduces. This distinction is important. Generic per-galaxy freedom would be expected to produce noise-like scatter centered on the universal best value; instead, the residual freedom in p is physically ordered. The cross-validation results (Section 6.2) and the BIC comparison against the universal-p free model (Section 6.4) indicate that the ordered structure is not merely a consequence of adding parameters.

7.6 Implications beyond the galactic scale

The SPARC sample probes the gravitational regime of disc galaxy outskirts, where Newtonian accelerations fall to ≈10⁻¹¹–10⁻⁹ m s⁻². The structured failure of universality reported here has implications across a wider range of scales, all of which lie beyond the present analysis but provide concrete tests of the framework.

Galaxy clusters. Clusters present the principal known stress test for MOND: even with the standard interpolation, residual unseen mass remains required at the centers of clusters (Famaey and McGaugh, 2012). A regime-dependent exponent of the form detected here would, under naive extrapolation of the monotonic regime progression of Section 6, predict effective exponents p > 1 in the high-acceleration cluster cores, which strengthens rather than weakens the residual-mass discrepancy. Conversely, in the cluster outskirts where g_outer drops to galaxy-comparable values, the empirical regime-dependent extrapolation predicts p ≈ 0.6–0.8, consistent with continuity into the galactic regime. This generates a falsifiable prediction: if the regime structure is fundamental, the effective MOND interpolation reconstructed from cluster outskirt kinematics should match the low-regime p inferred from SPARC. The available dynamical data on rich clusters provide the basis for such a test; a dedicated re-analysis under the constitutive family is straightforward and is left for future work.

Weak gravitational lensing. Galaxy–galaxy lensing measurements of the effective gravitational potential extend the acceleration regime well below the SPARC range, into g ≲ 10⁻¹² m s⁻². A regime-dependent exponent predicts a specific deviation from the universal MOND lensing signal in this deep regime: under extrapolation of the regime progression, the effective p in the deepest accessible bin should sit below the lowest SPARC-tercile median of 0.63. Modern wide-field weak-lensing surveys provide the data volume required to discriminate at this level.

8 Conclusion

To our knowledge, this work provides the first statistical test of the universality assumption underlying MOND-class interpolations at the per-galaxy level, within a parametric family that contains the canonical deep-MOND scaling as the p = 1 limit, validated across six complementary tests. We tested whether the SPARC galaxy sample requires a regime-dependent extension of the MOND interpolation framework. Within the phenomenological family β(g) = (a0/g)^p, the universal value p = 3/4 yields a baryonic Tully–Fisher slope of 3.5, substantially closer to the observed 3.47 ± 0.10 (ODR) than the canonical MOND value of 4, consistent with the p = 3/4 prediction within statistical uncertainty. Galaxy-by-galaxy, p = 3/4 is preferred over the MOND value p = 1 in 106 of 175 galaxies (binomial p = 0.006; the direction of the preference survives marginalization over distance and inclination uncertainties, with per-realization significance attenuated as expected).

The two universal-p choices reproduce the radial acceleration relation with statistically indistinguishable quality (rms 0.161 vs. 0.162 dex on 3,389 points).

Beyond the universal-p tests, fitting the constitutive family galaxy-by-galaxy reveals that the effective exponent is not constant. The population median is 0.82 with intrinsic dispersion σ = 0.53, which is three times broader than any single-p synthetic catalog can reproduce. The dispersion is structured: p correlates with the outer-acceleration regime at ρ = +0.49 (p < 10⁻⁹), and across three regime bins, the median p increases monotonically 0.63 → 0.83 → 1.11 (Kruskal–Wallis p = 6 × 10⁻⁸). A regime rule calibrated on half the sample predicts the held-out half with lower median χ²_red than MOND in 97 of 100 splits (99/100 under fixed-edge bins). Among seven alternative binning methods tested, including random control and tertile bins by five other physical variables, the acceleration-regime variable is qualitatively privileged at p < 10⁻⁴. An information-criterion comparison yields ΔBIC = −12,112 in favor of the regime rule over canonical MOND, with the critical control against the universal-p free model (ΔBIC = −3,620) ruling out “extra parameters” as the explanation: the BIC complexity penalty is dominated by the χ² reduction by roughly two orders of magnitude.

Taken together, these results indicate that the universal-p assumption underlying canonical MOND-style interpolations is observationally insufficient on SPARC. The data support a regime-dependent effective coupling, organized specifically by the local acceleration scale. We do not claim to identify the underlying physical mechanism: candidate interpretations span at least five families in the existing literature: curvature in the true MOND interpolation function that no single power-law captures, modified-inertia formulations of MOND, the external field effect, emergent-gravity scenarios, and baryonic-physics modifications to dark matter halo profiles that produce acceleration-dependent effective interpolations. The empirical phenomenon is the structured failure of universality itself; any of these mechanisms must reproduce it. Discriminating between them requires data and tests beyond the present sample, including cluster dynamics, weak gravitational lensing, ultra-diffuse galaxies, and the cosmological viability of any modified-gravity proposal that reproduces the regime structure (Section 7.6).

Whether the regime structure detected here reflects a fundamental property of the gravitational interaction, an emergent feature of baryon–dark matter interplay, or curvature in the true MOND interpolation function, the structured failure of universality is itself a testable prediction, one that any successor framework must reproduce. The present work locates the empirical phenomenon; the next step is to determine which class of physics generates it.

Generative AI statement

The author(s) declared that generative AI was used in the creation of this manuscript. Claude (Anthropic) was used for assistance with manuscript formatting and Portuguese-to-English language revision only. All scientific content, analyses, theoretical derivations, and conclusions are the responsibility of the authors. No generative AI was used in the design or execution of the analyses, or in the production of figures, tables, or scientific text content.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fspas.2026.1890948/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Jirong Mao, Chinese Academy of Sciences (CAS), China

Reviewed by: Yunchong Wang, Stanford University, United States

Arkady Bolotin, Ben-Gurion University of the Negev, Israel