Search for variability in Newton’s constant using local gravitational acceleration measurements

It is known for a while that the different measurements of Newton’s Gravitational constant ($G$) don’t agree with each other Bethapudi and Desai (2017); Mohr et al. (2012); Gillies (1997). There has also been a claim for a 5.9 year periodicity in $G$ Anderson et al. (2015), although this was rebutted using both Bayesian and frequentist analysis Pitkin (2015); Desai (2016) (see also  Schlamminger et al. (2015)). Some precision gravity experiments have also found tentative deviations from Newtonian gravity Krishak and Desai (2020), which could be connected to some of the unsolved problems in Cosmology such as Dark Energy, Dark Matter, Hubble tension, etc Perivolaropoulos and Kazantzidis (2019). However, we note that a large number of searches for the variation of $G$ have also been carried out using a plethora of astrophysical observations, such as helioseismology, lunar laser ranging, Type 1a Supernovae, binary pulsars, all of which have obtained null results Will (2006); Zyla et al. (2020). From these measurements, the variance in $\frac{\dot{G}}{G}$ is less than $10^{-12}$ Will (2006). At face value, these limits already rule out any putative variations in $G$, such as those claimed in  Anderson et al. (2015). However, any limit on the variation of $G$ from astrophysical or cosmological observation always has some uncertainty and model-dependence in it. Therefore, it behooves us to search for any variation or periodicity in $G$ using independent terrestrial probes.

Motivated by these considerations,  Dai (2021) (D21, hereafter) reconsidered the variability of $G$ using the variance of local Newtonian gravitational acceleration measurements ($\delta g/\bar{g}$). The key premise behind this work, was that the variance of $G$ can be bounded from the gravitational acceleration measurements:

where

Here, $\bar{g}$ refers to the mean of the gravitational acceleration during the observational period.

D21 used the $\frac{\delta g}{\bar{g}}$ data compiled by the International Geodynamics and Earth Tide Service (IGETS), using data from 16 observational stations. They considered two models for the variation of $G$, namely a sinusoidal and a linear model for the variation of the relative acceleration. D21 constructed a $\chi^{2}$ functional based on the residuals between the observed data and various theoretical models. For the sinusoidal model, they looked for a minima at 5.9 years. For the linear model, the value of the slope was constrained from the likelihood distribution. Thereafter, the relative amplitude of the sinusoidal model was constrained to be less than $3\times 10^{-9}$, whereas for the linear model $\dot{G}/G<10^{-10}\rm{yr^{-1}}$. In D21, the observed errors in $\frac{\delta g}{\bar{g}}$ were not considered for the analysis. The reason for this is that the variances of $\frac{\delta g}{\bar{g}}$ contain contributions not only from the acceleration measurement equipment, but also from local environment. These measurement errors could be smaller than the fluctuations from the local environmental effects. Therefore, since D21 were looking for long term trends in the variation of $G$, the measurement errors were ignored.

However, the median and mean fractional errors in $\frac{\delta g}{\bar{g}}$ were about 6% and 64% respectively with 38 data points having fractional errors greater than 100%. Furthermore, ignoring these errors we cannot get an independent assessment of the goodness of fit or robust error estimates on the model parameters. Therefore, it is imperative to redo the analysis in D21, after incorporating the measurement errors in the analysis. This also allows us independently test if the total errors are under-estimated or not.

In this work, we therefore revisit this issue by analyzing the same data (kindly provided to us by Dr. Dai), and by incorporating the errors in $\frac{\delta g}{\bar{g}}$ , while constructing the likelihood. We test for periodicities using two sinusoidal functions as well as a linear model for the variation in $\frac{\delta g}{\bar{g}}$ . We note, however, that the IGETS measurements analyzed in this work are not the only $\frac{\delta g}{\bar{g}}$ measurements available. There also exist gravitational acceleration measurements from several hundred FG5 gravimeters distributed around the world. We shall also analyze these data in a future work.

The manuscript is organized as follows. Our analyses is described in Sect. II and Sect. III, wherein we account for an additional unknown systematic error. We conclude in Sect. IV. We have also made our analysis codes used in this work publicly available.

The data showing the $\frac{\delta g}{\bar{g}}$ measurements, which we use for our fits is the same as in D21, and is shown in Fig. 1. We fit the data to the following four functions. The first function is a constant model. We then consider two sinusoidal functions, with the first one consisting of a mixed sine+cosine, and the second one consisting of sinusoidal with a phase. The last model we consider is a linear model. We enumerate these four functions below:

Note that Eq. 2 and Eq. 4 have also been used by D21. Eq. 3 is a variant of the sinuosoidal model used in D21, with a single sinusoid along with a phase. A similar model is usually used to fit annual modulation in dark matter experiments Krishak et al. (2020); Bhagvati and Desai (2021). We find the best-fit parameters for each of these models using the following likelihood ($L$)

where $y_{i}$, $\sigma_{i}$ denote the $\frac{\delta g}{\bar{g}}$ measurements and their errors, respectively; $y(t)$ represents the model proposed for the variation of $\frac{\delta g}{\bar{g}}$ . For all the models (except the constant model) the sampling of this likelihood has been done using the emcee Foreman-Mackey et al. (2013) MCMC sampler.

For each of the models, we construct the allowed regions in case of multiple parameters and tabulate the best-fit $\chi^{2}$ in Table 1. We now present our results for each of the models.

1. 1.

   Constant model For a constant model, the maximum likelihood solution is the weighted mean. The best-fit constant value is equal to $(-9.89\pm 6.2)\times 10^{-11}$. The minimum $\chi^{2}$ value is shown in Table 1. As we can see that the $\chi^{2}$/dof is greater than 20, implying that the constant model is not a good fit to the data.

2. 2.

   First sinusoidal model - Eq. 2 We now do a fit to Equation 2. This function is almost the same as that used in D21. The main difference being that in D21, a separate $a_{i}$ was used for every station, whereas here we assume $a_{i}$ is the same between all the stations. However, one caveat because of this assumption is that the uncertainty in $\bar{g}$, which is estimated for every station introduces an unknown shift in $a_{i}$. Since this model contains a sinusoidal function, the lower bound on the period should be greater than the sampling time, in order to avoid the overfitting problem which may happen for sinusoidal functions Andrae et al. (2010); Krishak and Desai (2019). The median sampling time is 12.4 days, while the mean sampling time is 23.5 days. Therefore, we choose uniform priors on $\omega$, viz. $\omega\in[0.55,45.92]$ rad/yr corresponding to a period of 11.4 years (half the total duration of the dataset) and 50 days, respectively. This range is broad enough to cover the previously observed peak at 5.9 years. We also choose uniform priors on $a$, $A_{c}$, and $A_{s}$ given by: $a\in[-2\times 10^{-8},2\times 10^{-8}]$, $A_{c}\in[0,2\times 10^{-8}]$, and $A_{s}\in[0,2\times 10^{-8}]$.

   Our results showing the 68%, 90% and 95% credible intervals for the free parameters of this model are shown in Fig. 2. We find that the marginalized credible intervals for $\omega$ show four peaks. These values are tabulated in Table 1. The best-fit values for the other parameters for the value of $\omega$ with the maximum likelihood are $a=(-2.97^{+1.00}_{-2.49})\times 10^{-10}$, $A_{c}=(0.12^{+0.10}_{-1.98})\times 10^{-9}$, and $A_{s}=(2.14^{+2.03}_{-0.15})\times 10^{-9}$. None of the peaks in $\omega$ correspond to a period of 5.9 years. The $\chi^{2}$ per degree of freedom for all these peaks are again greater than 20, indicating that none of the peaks in $\omega$ provide a good fit to the data.

   Similar to D21, we find that $\sqrt{A_{c}^{2}+A_{s}^{2}}$ for the best-fit likelihood is equal to $\sim 2.1\times 10^{-9}$, and hence is about four orders of magnitude more stringent than the periodicity ascertained from previous $G$ measurements Anderson et al. (2015).

3. 3.

   Second sinusoidal model - Eq. 3 We now do a similar fit to the second sinusoidal model listed in Eq. 3, which involves a phase factor. The priors for $a$, $w$, and $A$ are the same as those for $a$, $w$, and $A_{c}$, respectively used in Eq. 2. For $\phi$, we use a uniform prior given by $\phi\in[0,2\pi]$ . The best-fit values at the maximum likelihood obtained for these parameters are: $a=(-3.12^{+0.84}_{-2.32})\times 10^{-10}$, $A=(2.24^{+0.25}_{-0.12})\times 10^{-9}$, and $\phi=(4.30\pm 0.5)$ rad. The credible regions for the parameters of this model can be found in Fig. 3. Once again, we find four peaks for $\omega$, whose values are listed in Table 1. Again, we don’t see a peak at $\omega$, corresponding to the period of 5.9 years. Furthermore, the $\chi^{2}$ per degree of freedom for all the $\omega$ peaks are again greater than 20, indicating that none of the $\omega$ peaks provide an adequate fit to the data. The best-fit value of $A$ is $\sim 2.2\times 10^{-9}$. Therefore, for both the sinusoidal models analyzed, we get a limit, which is about four orders of magnitude more stringent than the amplitude of the sinusoidal variation inferred from $G$ measurements.

4. 4.

   Linear Model - Eq. 4 Finally, we fit the data to the linear model, specified in Eq. 4, which was also used in D21. The priors which we use for $a$ and $b$ are $a\in[-2\times 10^{-8},2\times 10^{-8}]$ and $b\in[-2\times 10^{-8},2\times 10^{-8}]$ /year, respectively. The best-fit credible regions are shown in Fig. 4. The best fit values for $a$ and $b$ are: $a=(-0.23\pm 0.01)\times 10^{-8}$ and $b=(2.03\pm 0.11)\times 10^{-10}$ /year. The best-fit $\chi^{2}$ values can be found in Table 1. Once again, we find that $\chi^{2}$/dof values are greater than 20, indicating that the linear model is also not an adequate fit. The limits which we obtain are not as stringent as those using astrophysical observations, although ours have been obtained using only terrestrial measurements.

   Then, using the same ansatz as D21, we can constrain the linear variation in $G$ using:

   $$\frac{\dot{G}}{G}\lessapprox b+\sigma_{b}=2.14\times 10^{-10}\text{yr.}$$

   This is about twice more stringent than the limit obtained in D21.

Therefore, in summary we find that none of the three time-dependent models (two sinusoidal based functions and a linear model) provide a good fit to the data, once we incorporate the errors in the measurements. However, even the constant model is not a good fit. The only way to get a good fit would be to add an additional systematic uncertainty (or intrinsic scatter) to the observed errors so that the $\chi^{2}$/dof is close to 1, similar to the analysis done in  Desai (2016); Pitkin (2015).

Since the $\chi^{2}$/dof for all the four models considered is much greater than one, one possibility is that the measurement errors are underestimated and there are additional sources of systematic errors. Some examples of these are enumerated in D21, and include tides, polar motion of the Earth, instrumental effects related to the speed of light, gradient height and vertical transfer, etc. Other possibilities include environment effects such as post-glacial rebound Bilker-Koivula et al. (2021) or water table changes Kazama et al. (2012). A detailed modelling of these additional systematics is beyond the scope of this work. Therefore, we redo our search for variability by adding an unknown intrinsic scatter term ($\sigma_{scatter}$) in quadrature to the measurement errors and marginalize over it, similar to various cosmological analyses Pradyumna et al. (2021). The premise is that this intrinsic scatter encapsulates all other unknown systematic errors in the dataset. We use log-uniform priors for $\sigma_{scatter}\in[10^{-9},10^{1}]$. We also carry out model comparison between the constant and the time-varying models using AIC and BIC information theory-based tests Krishak and Desai (2020), where AIC/BIC are defined as follows:

where $k$ and $n$ correspond to the number of free parameters and the number of data points, respectively. The models with smaller AIC/BIC values are the favored ones and the significance can be assessed using qualitative strength of evidence rules Krishak and Desai (2020).

The $\Delta$AIC and the $\Delta$BIC values, calculated with respect to the constant model, along with the $\chi^{2}/DOF$ values are given in Table 2. We now summarize our results for each of these models with the addition of intrinsic scatter.

1. 1.

   Constant Model (with intrinsic scatter): We choose uniform priors for $C\in[-2\times 10^{-8},2\times 10^{-8}]$, where $C$ is the constant value for $\frac{\delta g}{\bar{g}}$ . The best-fit parameters are $C=(1.27\pm 3.20)\times 10^{-10}$ and $\log_{10}(\sigma_{scatter})=(-8.24\pm 0.02)$. The best-fit $\chi^{2}/DOF$ values are close to 1.

2. 2.

   First Sinusoidal Model (with intrinsic scatter): For this analysis also, we includ negative priors on $A_{c}$ and $A_{s}$. Uniform priors are taken for $A_{c}\in[-2\times 10^{-8},2\times 10^{-8}]$, and $A_{s}\in[-2\times 10^{-8},2\times 10^{-8}]$. The priors on $a$ and $\omega$ are the same as those mentioned in Section II. The best-fit values are $a=(1.76^{+3.26}_{-3.10})\times 10^{-10}$, $A_{c}=(-0.59^{+0.97}_{-1.21})\times 10^{-9}$, $A_{s}=(-1.32^{+1.23}_{-1.82})\times 10^{-10}$, and $\log_{10}(\sigma_{scatter})=(-8.25\pm 0.02)$. We see one prominent peak for $\omega$ occurring at $\omega=6.10$ rad/year, and another small peak occurring at $\omega=25.59$ rad/year, if we ignore the peak occurring close to the lower prior on $\omega$, corresponding to a time-period of $11.4$ years. This is contrary to the multiple peaks observed in Figure 2, whose biggest peak occurs at $\omega=24.47$ rad/year. Once again, we do not find any peak corresponding to the period of 5.9 years.

3. 3.

   Second Sinusoidal Model (with intrinsic scatter): The prior for $A$ is the same as the prior used on $A_{c}$, and the priors on $a$, $\omega$, and $\phi$ remain the same as those considered for the analysis in Sect II. We find a peak in $\omega$ at $\omega=6.11$ rad/year, if we ignore the peak corresponding to the $T=11.4$ years (close to the lower bound on $\omega$). This is different from the two peaks observed in the analysis without intrinsic scatter where the biggest peak occurs at $\omega=24.50$ rad/year, as shown in Figure 3. Again, we do not find a peak in the likelihood at 5.9 years. The best-fit values of the other parameters are $a=(1.91\pm 3.06)\times 10^{-10}$, $A=(0.16^{+2.17}_{-2.38})\times 10^{-9}$, $\phi=(3.38^{+1.95}_{-1.72})$, and $\log_{10}(\sigma_{scatter})=(-8.25\pm 0.02)$.

4. 4.

   Linear Model (with intrinsic scatter): The priors on $a$ and $b$ remain the same as those taken previously in Section II. The best-fit values are $a=(-3.28\pm 0.51)\times 10^{-9}$, $b=(3.67\pm 0.45)\times 10^{-10}$, and $\log_{10}(\sigma_{scatter})=(-8.28\pm 0.02)$.

Therefore, the best-fit intrinsic scatter which we obtain for all the four models is $10^{-8}$ and is about 2.8 times the average error in $\frac{\delta g}{\bar{g}}$ measurements. With the addition of intrinsic scatter, the best-fit $\chi^{2}$/dof is now close to one. However, this is a consequence of fitting for an intrinsic scatter in the likelihood. If we compare the constant model to the two sinusoidal models, we find that the constant model is decisively favored over the sinusoidal model (with $\Delta$ AIC/BIC $>10$). However, when we compare the constant model with the linear model, we find that the linear time-varying model is marginally favored over the constant model. Since both $\Delta$AIC and $\Delta$BIC values between the constant and linear model are less than 10, additional data would be necessary to ascertain if this difference persists and is significant.

In a recent work, D21 analyzed the IGETS $\frac{\delta g}{\bar{g}}$ measurements, in order to assess the variability and periodicity in $G$ measurements previously reported in literature Anderson et al. (2015). Another impetus for searching for a variation in $G$ has been to address cosmological questions such as Dark energy, Dark matter, Hubble tension, etc. Although, the limits on variation of $G$ from astrophysical observations Will (2006) already rule out the periodicities in $G$ claimed in  Anderson et al. (2015), by many orders of magnitude, it is important to look for a variability in $G$ using independent terrestrial probes. This was the main motivation for D21 as well as this work.

The main premise for the analysis in D21 is that the $\frac{\delta g}{\bar{g}}$ measurements provide an upper bound to the fractional variation in $G$. For this purpose, D21 analyzed the $\chi^{2}$ functional between the data and measurements, and looked for a minimum at 5.9 years for the sinusoidal functions. They do not find any such periodicity or linear time-dependent variations. The limits which they obtain rule out the claims in  Anderson et al. (2015).

However, D21 did not incorporate the errors in $\frac{\delta g}{\bar{g}}$ , while evaluating the $\chi^{2}$. In this work, we improve upon the work in D21 by considering the measurement errors in all our analyses. We considered four functions to model the putative behavior of $\frac{\delta g}{\bar{g}}$ : a constant model, two sinusoidal models (Eq. 2 and Eq. 3), and a linear model (Eq. 4). For each of these models, we maximized a log-likelihood (Eq. 5) function by incorporating the errors to get the allowed regions on the parameters. These allowed regions can be found in Fig. 2, Fig. 3, and Fig. 4 for the two sinusoidal and linear models. The best-fit values for $\chi^{2}$ for each of these models is tabulated in Table 1. As we see, the $\chi^{2}$/dof for all the models hitherto analyzed in this work are greater than 20, indicating that none of them provide a good fit to the data.

Given these results, we therefore redo our search for a variability in $\frac{\delta g}{\bar{g}}$ by incorporating an unknown intrinsic scatter, while fitting all the four models. The best-fit results from these analyses are tabulated in Table 2. The best-fit intrinsic scatter we obtain is about 2.8 times larger than the errors in $\frac{\delta g}{\bar{g}}$ measurements. We carry out a model comparison between the four models using AIC and BIC. We find that the constant models are decisively favored over the sinusoidal ones. However, we find that the linear model is marginally favored over the constant model. Additional data would be needed to ascertain if the difference becomes significant.

Based on our analysis without intrinsic scatter, we find from our fits to the two sinusoidal models that the relative variance of $G$ for any sinusoidal variation is $\leqslant 2\times 10^{-9}$, which is about four orders of magnitude more stringent than the amplitude of periodic variations inferred from previous $G$ measurements. Using the best-fit linear model, we can constrain $\frac{\dot{G}}{G}<2.14\times 10^{-10}$/year, which is about a factor of two more stringent than the corresponding limit obtained in D21. Therefore, although our limits on variation of $G$ are not as stringent as those obtained from astrophysical observations, they are still slightly better than those in D21, after we incorporate the errors in measurements.

To conclude, we concur with D21 that there is no evidence for a periodicity or a linear variation in the $\frac{\delta g}{\bar{g}}$ measurements, even after incorporating the errors, but without adding an intrinsic scatter. When we add an intrinsic scatter, the linear time evolution model for $G$ is marginally favored over a constant model.

In the spirit of open science analysis, we have made all codes used in this work publicly available, which can be found at the github link https://github.com/srinixbhagvati/Search-for-Variability-in-Newtons-Constant