Hints for a gravitational constant transition in Tully-Fisher data

We use an up to date compilation of Tully-Fisher data to search for transitions in the evolution of the Tully-Fisher relation. Using an up to date data compilation, we find hints at $\approx 3\sigma$ level for a transition at critical distances $D_c \simeq 9 Mpc$ and $D_c \simeq 17 Mpc$. We split the full sample in two subsamples according to the measured galaxy distance with respect to a splitting distance $D_c$ and identify the likelihood of the best fit slope and intercept of one sample with respect to the best fit corresponding values of the other sample. For $D_c \simeq 9 Mpc$ and $D_c \simeq 17 Mpc$ we find a tension between the two subsamples at a level of $\Delta \chi^2>17\; (3.5\sigma)$. Using a Monte-Carlo simulation we demonstrate that this result is robust with respect to random statistical and systematic variations of the galactic distances. If the tension is interpreted as due to a gravitational strength transition, it would imply a shift of the effective gravitational constant to lower values for distances larger than $D_c$ by $\frac{\Delta G}{G}\simeq -0.1$. Such a shift is of the anticipated sign and magnitude but at somewhat lower distance (redshift) than the gravitational transition recently proposed to address the Hubble and growth tensions ($\frac{\Delta G}{G}\simeq -0.1$ at transition redshift $z_t\lesssim 0.01$ ($D_c\lesssim 40 Mpc$)).


Introduction
The Tully-Fisher relation (TFR) [1] has been proposed as an empirical relation that connects the intrinsic optical luminosity of spiral galaxies with their observed maximum velocity v rot in the rotation curve as follows: (1) where s 4 is the slope in a logarithmic plot of (1), and A is a constant (log(A) is the zero point or intercept). The constants s and A appear to depend very weakly on galaxy properties, including the mass to light ratio, the observed surface brightness, the galactic profiles, H I gas content, size, etc. [2]. They clearly also depend, however, on the fundamental properties of gravitational interactions as demonstrated below.
The baryonic Tully-Fisher relation (BTFR) is similar to Equation (1) but connects the galaxy's total baryonic mass (the sum of mass in stars and H I gas) M B with the rotation velocity as follows: where A B 50M km −4 s 4 [3]. This allows to include gas-rich dwarf galaxies that appear in groups and have stellar masses below 10 9 M .
A simple heuristic analytical derivation for the BTFR can be obtained as follows [4]. Consider a star in a circular orbit of radius R around a galactic mass M rotating with velocity v. Then, the following holds: where G eff is the effective Newton's constant involved in gravitational interactions and S the surface density S ≡ M/R 2 , which is expected to be constant [5]. From Equations (2) and (3), the following is anticipated: Therefore, the BTFR can, in principle, probe both galaxy formation dynamics (through, e.g., S) and possible fundamental constant dynamics (through G eff ). An interesting feature of the BTFR is that despite the above heuristic derivation, it appears to be robust, even in cases when the galaxy sample includes low S and/or varying S galaxies [6,7]. In fact, no other parameter appears to be significant in the BTFR. The BTFR has been shown to have lower scatter [2,8,9] than the classic stellar TFR and also to be applicable for galaxies with stellar masses lower than 10 9 M . It is also more robust than the classic TFR [10][11][12][13] since the parameters A B (intercept) and s (slope) are very weakly dependent on galactic properties, such as size and surface brightness [2].
The low scatter of the BTFR and its robustness make it useful as a distance indicator for the measurement of the Hubble constant H 0 . A calibration of the BTFR using Cepheid and TRGB distances leads to a value of H 0 = 75 ± 3.8 km s −1 Mpc −1 [14].
This value of H 0 is consistent with local measurements of H 0 , using SnIa calibrated with Cepheids (H 0 = 73.2 ± 1.3 km s −1 Mpc −1 ) [15], but is in tension with the value of H 0 obtained using the early time sound horizon standard ruler calibrated using the CMB anisotropy spectrum in the context of the standard ΛCDM model (H 0 = 67.36 ± 0.54 km s −1 Mpc −1 ) [16]. The tension between the CMB and Cepheid calibrators is at a level larger than 4σ and constitutes a major problem for modern cosmologies (for a recent review and approaches see Refs. [17][18][19][20][21]).
The Hubble tension may also be viewed as an inconsistency between the value of the standardized SnIa absolute magnitude M calibrated using Cepheids in the redshift range 0 < z < 0.01 (distance ladder calibration) and the corresponding M value calibrated using the recombination sound horizon (inverse distance ladder calibration) for 0.01 < z < z rec . Thus, a recently proposed class of approaches to the resolution of the Hubble tension involves a transition [22,23] of the standardized intrinsic SnIa luminosity L and absolute magnitude M at a redshift z t 0.01 from M = (−19.24 ± 0.037) mag for z < z t (as implied by Cepheid calibration) to M = (−19.4 ± 0.027) mag for z > z t (as implied by CMB calibration of the sound horizon at decoupling) [24]. Such a transition may occur due to a transition in the strength of the gravitational interactions G eff , which modifies the SnIa intrinsic luminosity L by changing the value of the Chandrasekhar mass. The simplest assumption leads to [25,26], even though corrections may be required to the above simplistic approach [27].
The weak evolution and scatter of the BTFR can be used as a probe of galaxy formation models as well as a probe of possible transitions of fundamental properties of gravitational dynamics since the zero point constant A B is inversely proportional to the square of the gravitational constant G. Previous studies investigating the evolution of the best-fit zero point log A B and slope s of the BTFR have found a mildly high z evolution of the zero point from z 0.9 to z 2.3 [28], which was attributed to the galactic evolution inducing a lower gas fraction at low redshifts after comparing with the corresponding evolution of the stellar TFR (STFR), which ignores the contribution of gas in the galactic masses.
Ref. [28] and other similar studies assumed a fixed strength of fundamental gravitational interactions and made no attempt to search for sharp features in the evolution of the zero point. In addition, they focused on the comparison of high redshift with low redshift effects without searching for possible transitions within the low z spiral galaxy data. Such transitions, if present, would be washed out and hidden from these studies, due to averaging effects. In the present analysis, we search for transition effects in the BTFR at z 0.01 (distances D 40 Mpc), which may be due to either astrophysical mechanisms or to a rapid transition in the strength of the gravitational interactions G eff , due to fundamental physics.
In many modified gravity theories, including scalar tensor theories, the strength of gravitational interactions G eff measured in Cavendish-type experiments measuring force F between masses (F = G eff m 1 m 2 r 2 ), is distinct from the Planck mass corresponding to G N that determines the cosmological background expansion rate (H 2 = 8πG N 3 ρ tot ). For example, in scalar tensor theories involving a scalar field φ and a non-minimal coupling F(φ) of the scalar field to the Ricci scalar in the Lagrangian, the gravitational interaction strength is as follows [29]: while the Planck mass related G N is as follows: Most current astrophysical and cosmological constraints on Newton's constant constrain the time derivative of G eff at specific times, assume a smooth power-law evolution of G eff , or constrain changes of the Planck mass-related G N instead of G eff (CMB and nucleosynthesis constraints [30]). Therefore, these studies are less sensitive in the detection of rapid transitions of G eff at low z.
The current constraints on the evolution of G eff and G N are summarized in Table 1, where we review the experimental constraints from local and cosmological time scales on the time variation of the gravitational constant. The methods are based on very diverse physics, and the resulting upper bounds differ by several orders of magnitude. Most constraints are obtained from systems in which gravity is non-negligible, such as the motion of the bodies of the solar system, and the astrophysical and cosmological systems. They are mainly related in the comparison of a gravitational time scale, e.g., period of orbits, with a non-gravitational time scale. One can distinguish between two types of constraints, from observations on cosmological scales and on local (inner galactic or astrophysical) scales. The strongest constraints to date come from lunar ranging experiments.
In the first column of Table 1, we list the used method. The second column contains the upper bound ∆G G max of the fractional change of G during the corresponding timescale. Most of these bounds assume a smooth evolution of G. In the third column, we present the upper bound on the normalized time derivative Ġ G max . The fourth column is an approximate time scale over which each experiment is averaging each variation, and the fifth column refers to the corresponding study where the bound appears. Entries with a star ( * ) indicate constraints on G N , while the rest of the constraints refer to the gravitational interaction constant G eff .
In the present analysis, we search for a transition of the BTFR best-fit parameter values (intercept and slope) between data subsamples at low and high distances. We consider sample dividing distances D c ∈ [2, 60] Mpc, using a robust BTFR dataset [12,[31][32][33], which consists of 118 carefully selected BTFR datapoints, providing distance, rotation velocity baryonic mass (D, V f , M B ) as well as other observables with their 1σ errorbars. We focus on the gravitational strength Newton constant G eff and address the following questions: • Are there hints for a transition in the evolution of the BTFR? • What constraints can be imposed on a possible G eff transition, using BTFR data? • Are these constraints consistent with the level of G eff required to address the Hubble tension?
The structure of this paper is the following: In the next section, we describe the datasets involved in our analysis and present the method used to identify transitions in the evolution of the BTFR at low z. We also show the results of our analysis. In Section 3, we summarize, present our conclusions and discuss possible implications and extensions of our analysis.

Search for Transitions in the Evolution of the BTFR
The logarithmic form of the BTFR (Equation (2)) is as follows: and a similar form for the TFR. Due to Equation (4), the intercept b ≡ logA B depends on both the galaxy formation mechanisms through the surface density S and on the strength of gravitational interactions through G eff . A controversial issue in the literature is the type of possible evolution of the slope and intercept of the TFR and the BTFR. Most studies have searched for possible evolution in high redshifts (redshift range z ∈ [0, 3]) with controversial results. For example, several studies found no statistically significant evolution of the intercept of the TFR up to redshifts of z ∼ 1.7 [52][53][54][55][56][57][58], while other studies found a negative evolution of the intercept up to redshift z 3 [59][60][61][62][63][64][65][66]. Similar controversial results in high z appeared for the BTFR, where [60] found no significant evolution of the intercept since z 0.6, while [64] found a positive evolution of the intercept between low-z galaxies and a z 2 sample. In addition, cosmological simulations of disc galaxy formation based on cosmological N-body/hydrodynamical simulations have indicated no evolution of the TFR based on stellar masses in the range z ∈ [0.1] [67], indicating also that any observed evolution of the TFR is an artifact of the luminosity evolution.
These studies have focused mainly on comparing high-z with low-z samples, making no attempt to scan low redshift samples for abrupt transitions of the intercept and slope. Such transitions would be hard to explain in the context of known galaxy formation mechanisms but are well motivated in the context of fundamental gravitational constant transitions, which may be used to address the Hubble tension [22,23]. Thus, in this section, we attempt to fill this gap in the literature.
We consider the BTFR dataset shown in Appendix A based on the data from [12,[31][32][33] of the flat rotation velocity of galaxies vs. the baryonic mass (stars plus gas) consisting of 118 datapoints, shown in Table A1. The sample is restricted to those objects for which both quantities are measured to better than 20% accuracy and includes galaxies in the approximate distance range D ∈ [1,130] Mpc. This is a robust low z dataset (z < 0.1) with low scatter showing no evolution of velocity residuals as a function of the central surface density of the stellar disks.
Our analysis is distinct from previous studies in two aspects: • We use an exclusively low z sample to search for BTFR evolution. • We focus on a particular type of evolution: sharp transitions of the intercept and slope.
In this context, we use the dataset shown in Table A1 of Appendix A [12,[31][32][33], consisting of the distance D, the logarithm of the baryonic mass logM B and the logarithm of the asymptotically flat rotation velocity logv rot of 118 galaxies along with 1σ errors. We fix a critical distance D c and split this sample in two subsamples Σ 1 (galaxies with D < D c ) and Σ 2 (galaxies with D > D c ). For each subsample, we use the maximum likelihood method [68] and perform a linear fit to the data setting y i = log(M B ) i , x i = log(v rot ) i , while the parameters to fit are the slope s and the intercept b of Equation (7). Thus, for each sample j (j = 0, 1, 2 with j = 0 corresponding to the full sample and j = 1, 2 corresponding to the two subsamples Σ 1 and Σ 2 ), we minimize the following: with respect to the slope s j and intercept b j . We fix the scatter to σ s = 0.077, obtained by demanding that ,min is the minimized value of χ 2 for the full sample and N 0 is the number of datapoints of the full sample. We thus find the best fit values of the parameters s j and b j , (j = 0, 1, 2) and also construct the 1σ − 3σ likelihood contours in the s − b parameter space for each sample (full, Σ 1 and Σ 2 ) for a given value of D c . We then evaluate the ∆χ 2 kl (D c ) of the best fit of each subsample k, best fit with respect to the likelihood contours of the other subsample l. Using these values, we also evaluate the σ-distances (d σ,kl (D c ) and d σ,lk (D c )) and conservatively define the minimum of these σ-distances as follows: For example, for the σ-distance of the best fit of Σ 1 with respect to the likelihood contours of Σ 2 , we have the following: and d σ, 12 is obtained as a solution of the following equation [68]: where Q −1 is the inverse regularized incomplete Gamma function, M is the number of parameters to fit (M = 2 in our case) and Erf is the error function.     Figure 3 alongside with their 1σ errors. The minimum ∆χ 2 between the best fits of the two samples is also shown. The corresponding σ-tension in parenthesis is obtained in the context of two free parameters from Equation (11). Notice that, even though the parameter values appear to be consistent, the value of ∆χ 2 between the subsamples reveals the tension at D c = 9 Mpc and D c = 17 Mpc. The typical qualitative feature of d σ (D c ) corresponding to the real sample disappears if we homogenize the sample by randomizing both the velocities and the galactic masses, using the measured values of the velocities and the estimated values of the galactic masses in the context of the best-fit BTFR. In order to construct such a homogenized BTFR sample from the real sample, we use the following steps: • We assign to each galaxy a randomly chosen distance obtained from a Gaussian distribution with mean equal to the measured distance and standard deviation equal to the 1σ error of the measured distance. • We assign to each galaxy a randomly chosen logv rot obtained from a Gaussian distribution with mean equal to the measured logv rot and standard deviation equal to the 1σ error of the measured logv rot .
• For each galaxy, we use the random logv rot obtained in the previous step to calculate the corresponding BTFR logM B , using the best-fit slope and intercept of the real full dataset (first row of Table 2). We then obtain a random logM B for each galaxy from a Gaussian distribution with mean equal to the BTFR calculated logM B and standard deviation equal to the 1σ error of the measured logM B . • We repeat the above process 100 times, thereby generating 100 homogeneous Monte Carlo samples (HMCS) based on the SPARC dataset. • For each HMCS, we find the σ distances d σ (D c ) and for each D c , we find the mean σ distance and its standard deviation over the 100 HMCS. We thus construct the orange region in Figure 2    . The best-fit contours of the slope and intercept for the entire dataset, as well for 3 different cases of split distance (D c ). The red contours correspond to the dataset with galaxies that have a distance below D c , whereas the cyan contours correspond to galaxies with distances above D c .
Clearly, the forms of d σ (D c ) generated from the homogenized Monte Carlo samples have the expected property to be confined mainly between 0σ and 2σ in contrast to the real measured sample, where d σ (D c ) extends up to 4σ or more. Thus, the real dataset is statistically distinct from a homogeneous BTFR dataset.
The two maxima of d σ are more clearly illustrated in Figure 3, where the likelihood contours are shown in the parameter space s (slope)-b (intercept) for the full sample (upper left panel) and for three pairs of subsamples Σ i , including those corresponding to the peaks shown in Figure 1 (D c = 17 and D c = 9). For both d σ maxima, the tension between the two best-fit points is mainly due to the different intercepts, while the values of the slope are very similar for the two subsamples. In contrast, for D c = 40 Mpc, where the σ distance is much lower (about 1σ, lower right panel), both the slope and the intercept differ significantly in magnitude but the statistical significance of this difference is low. Notice that the use of different statistics, such as the 1σ range of the best-fit intercept and slope shown in Table 2, or the level of likelihood contour overlap in Figure 3 would not reveal the tension between far and nearby subsamples. In contrast, the σ-distance statistic demonstrates the effect and the Monte Carlo results of Figure 2 verify the fact that such a large σ-distance would be rare in the context of a homogeneous sample.
The statistical significance of the different Tully-Fisher properties between near and far galaxies, which abruptly disappears for dividing distance D c 20 Mpc, could be an unlikely statistical fluctuation, a hint for systematics in the Tully-Fisher data 1 , an indication for an abrupt change in the galaxy evolution or a hint for a transition in the values of fundamental constants and, in particular, the strength of gravitational interactions G eff . The best-fit values of the intercept and the slope for the cases shown in Figure 3 are displayed in Table 2 along with their 1σ errors.    The best-fit logM B − logv rot lines corresponding to Equation (7) for the near-far galactic subsamples are shown in Figure 4, superimposed with the datapoints (red/blue correspond to near/far galaxies). The full dataset corresponds to the upper-left panel. The difference between the two lines for D c = 9 Mpc and D c = 17 Mpc is evident, even though their slopes are very similar. The statistical significance of this difference disappears for larger values of the splitting distance (e.g., D c = 40 Mpc), even though the slopes of the two lines become significantly different in this case.  Figure 5. The distances to galaxies beyond 20 Mpc are determined using the Hubble flow with H 0 = 73 km/sec Mpc, and thus, there is no effect of their peculiar velocities. Galaxies closer than about D 20 Mpc are clearly not in the Hubble flow and their redshift is affected significantly by their in-falling peculiar velocities, which tend to reduce their cosmological redshifts. The detected transitions at about 9 Mpc and 17 Mpc correspond to cosmological redshifts of z 0.005, which is lower than the transition redshift required for the resolution of the Hubble tension (z t 0.07 is the upper redshift of SnIa-Cepheid host galaxies).
In the context of the above-described analysis, we have ignored the possible systematic uncertainties induced on the estimated baryonic masses M B , due to systematic uncertainties in the measurement of galactic distances. In particular, different sub-samples of galaxies in the SPARC database are affected by different systematic uncertainties. The SPARC sample includes galaxies with both direct and indirect distance measurements. Thus, the identified mismatch of the Tully-Fisher parameters between low-and highdistance subsamples could, in principle, be due to such a systematic uncertainty of the galactic baryonic masses of Hubble flow galaxies. In order to examine this possibility, we have constructed new Monte Carlo samples where we not only vary randomly the distances but also add a fixed shift of ∆logM B along the vertical axis (mass) for all the datapoints where the mass is estimated using the Hubble flow with H 0 = 73 km/s/Mpc. The distances of these points are calculated using the Hubble flow, assuming H 0 = 73 km s −1 Mpc −1 , and correcting for Virgo-centric infall. We have considered four cases of systematic shifts (fixed values of ∆logM B ):−0.1 dex, −0.05 dex, +0.05 dex and +0.1 dex. The results for the σ-distance ranges in terms of the splitting distance D c for each one of the above four cases are shown in Figure  6. The corresponding likelihood contours for the subsamples corresponding to D c = 9 Mpc (maximum mismatch) are shown in Figure 7. Clearly, the mismatch features at D c = 9 Mpc and D c = 17 Mpc remain in all four cases that explore this type of systematic uncertainty. In particular, the 9 Mpc peak height varies from about 4σ for ∆logM B = 0.1 dex to about 3σ for ∆logM B = −0.1 dex. We thus conclude that this type of systematic uncertainty is unable to wash out the mismatch effect we have identified.
If the intercepts' transitions are interpreted as being due to a transition in G eff , we can use Equation (4) along with the observed intercept transition amplitude shown in Table 2 to identify the magnitude and sign of the corresponding G eff transition. The intercept transition at D c = 17Mpc indicated in Table 2 corresponds to the following: Since A B is found to be higher at larger distances (early times), G eff should be lower, due to Equation (4). The corresponding fractional change in G eff is easily obtained by differentiating the logarithmic form of Equation (4) as follows: This sign (weaker gravity at early times) and magnitude of the G eff transition is consistent with the gravitational transition required for the resolution of the Hubble and growth tensions in the context of the mechanism of Ref. [22].  Figure 6. The red contours correspond to the dataset with galaxies that have distance below 9 Mpc, whereas the cyan contours correspond to galaxies with distances above 9 Mpc. The σ distance between the two best fits varies between 3σ and 4σ.

Conclusions -Discussion
We used a specific statistic on a robust dataset of 118 Tully-Fisher datapoints to demonstrate the existence of evidence for a transition in the evolution of BTFR. This evidence was verified by a wide range of Monte Carlo simulations that compare the real dataset with corresponding homogenized datasets constructed using the BTFR. It indicates a transition of the best-fit values of BTFR parameters, which is small in magnitude but appears at a level of statistical significance of more than 3σ. It corresponds to a transition of the intercept of the BTFR at a distance of D c 9 Mpc and/or at D c 17 Mpc (about 80 million years ago or less). Such a transition could be interpreted as a systematic effect or as a transition of the effective Newton constant with a 10% lower value at early times, with the transition taking place about 80 million years ago or less. The amplitude and sign of the gravitational transition are consistent with a recently proposed mechanism for the resolution of the Hubble and growth tensions [22,23]. However, the time of the transition is about 60 million years later than the time suggested by the above mechanism (100-150 million years ago corresponding to D c 30-40 Mpc and z 0.007-0.01).
The effect shown in our analysis could be attributed to causes other than a gravitational transition. One such possible cause would be the presence of systematic errors affecting the estimate of galactic masses or rotation velocities for particular distance ranges. Even if this is the case, it is important to point out these inhomogeneities, which may require further analysis to identify their origin. Alternatively, if the causes of the detected mismatch are physical, they could also be due to variation of conventional galaxy formation mechanisms, which may involve other types of modifications of gravitational physics (e.g., effects of MOND gravity). The BTFR is an observationally tight empirical correlation and has therefore been used as a test of various modified gravity models (Refs. [71][72][73] offer comprehensive reviews on the cosmological implications of such models), including modified Newtonian dynamics (MOND) [74,75] and Grumiller modified gravity [76].These models have been shown to be consistent with BTFR for specific values of their acceleration parameters. The BTFR has also been used as a test of the properties of Cold Dark Matter and galaxy formation mechanisms in the context of ΛCDM [77,78].
An interesting effect in the direction of the one observed in our analysis was also reported in Ref. [79]. There, the authors found a transition of the Cepheid magnitude behavior in the range of 10-20 Mpc, which could explain the Hubble tension (see Figure 4 of Ref. [79]). The authors claimed that this transition is probably due to dust property variation, but there is currently a debate on the actual cause of this mismatch.
An important extension of this analysis is the search for similar transition signals and constraints in other types of astrophysical and geophysical-climatological data of Earth paleontology. For example, a wide range of solar system anomalies were discussed in Ref. [80], which could be revisited in the context of the gravitational transition hypothesis. Of particular interest, for example, is the 'Faint young Sun paradox' [81], which involves an inconsistency between geological findings and solar models about the temperature of the Earth about 4 billion years ago. Another interesting extension of this study would be the use of alternative methods for the identification of transition-like features in the data, e.g., the use of a Bayesian analysis tool, such as the internal robustness described in Refs. [82,83].
Alternatively, other astrophysical relations that involve gravitational physics, such as the Faber-Jackson relation between intrinsic luminosity and velocity dispersion of elliptical galaxies or the Cepheid star period-luminosity relation, could also be screened for similar types of transitions as in the case of BTFR. For example, the question to address in the Cepheid case would be the following: 'What constraints can be imposed on a transition-type evolution of the absolute magnitude (M v )-period (P) relation of Population I Cepheid stars?' This relation may be written as follows: where s = −2.43 ± 0.12 and b = −4.05 ± 0.02 [84,85].
In conclusion, the low z gravitational transition hypothesis is weakly constrained in the context of current studies but it could lead to the resolution of important cosmological tensions of the standard ΛCDM model. We have demonstrated the existence of hints for such a transition in the evolution of the Tully-Fisher relation.
Author Contributions: LP contributed in the conceptualization, the methodology, the writing, as well as the general supervision of the project. GA contributed in the formal data analysis and the writing. IA contributed in the data curation and the literature investigation. All authors have read and agreed to the published version of the manuscript.  A possible source of systematics is the Malmquist bias, which would imply that the detected more distant galaxies are also more massive and may, therefore, display different slopes and intercepts in different mass bins [69,70].