Time-dependent $G$ in Einstein's equations as an alternative to the cosmological constant

In this work, we investigate cosmologies where the gravitational constant varies in time, with the aim of explaining the accelerated expansion without a cosmological constant. We achieve this by considering a phenomenological extension to general relativity, modifying Einstein's field equations such that $G$ is a function of time, $G(t)$, and we preserve the geometrical consistency (Bianchi identity) together with the usual conservation of energy by introducing a new tensor field to the equations. In order to have concrete expressions to compare with cosmological data, we posit additional properties to this tensor field, in a way that it can be interpreted as a response of spacetime to a variation of $G$. Namely, we require that the energy this tensor represents is nonzero only when there is a time variation of $G$, and its energy depends on the scale factor only because of its coupling to $G$ and the matter and radiation energy densities. Focusing on the accelerated expansion period, we use type Ia supernovae and baryon acoustic oscillation data to determine the best fit of the cosmological parameters as well as the required variation in the gravitational constant. As a result, we find that it is possible to explain the accelerated expansion of the Universe with a variation of $G$ and no cosmological constant. The obtained variation of $G$ stays under 10 \% of its current value in the investigated redshift range and it is consistent with the local observations of $\dot{G}/G$.


I. INTRODUCTION
The concordance model in cosmology, ΛCDM, is extremely successful in being able to explain most of the current cosmological observations with great precision [1][2][3][4][5][6]. On the other hand, this model also has important problems, one of which is its inability to explain the nature of the titular Λ, or the cosmological constant, which remains to be an ad hoc addition to general relativity, employed in order to explain the late stage accelerated expansion of the Universe. While this constant behaves in the same way as a vacuum energy, its value is many orders of magnitude smaller than the estimations of quantum field theory [7]. Finding a cosmological solution with Λ = 0 solves one aspect of this problem, since this solution would be compatible with a renormalization to zero, a more natural value compared to the current estimations of Λ [8].
In this paper, we present an alternative picture, in which the accelerated expansion attributed to the cosmological constant appears naturally as a result of a variation of G in a relativistic model of gravitation. We achieve this by positing a phenomenological time variation of the gravitational constant G in Einstein's field * ekimtaylan@gmail.com equations. As is well-known [9,10], and as we will later demonstrate more clearly, this scenario creates a coupling between the matter and radiation energy density and G, which breaks the energy conservation of matter and radiation. We solve this issue in a general way by adding a new dynamical term in Einstein's equations. Our approach here is similar to Jordan-Brans-Dicke (JBD) theories [11,12], in that we decouple the density evolution of the matter and radiation from the variation in the gravitational constant. However, we determine the gravitational constant phenomenologically from the observations, instead of obtaining it according to first principles (unlike the JBD Lagrangian, which is obtained from Mach's principle). We focus on the accelerated expansion period and show that even a small variation of the gravitational constant can produce a similar effect to a cosmological constant in the present epoch.
In order to obtain specific cosmological models to test against the data, we make additional physical assumptions about the new dynamical term. We consider this tensor to represent only the reaction of the spacetime geometry to any variation of the gravitational constant, such that it can be interpreted as a property of the spacetime itself. Therefore, we require (a) that, apart from its dependence on other terms, the energy represented by this tensor should not be directly affected by the expansion of the space and (b) that only the terms coupled to the evolution of the gravitational constant should be present. After this, we propose a Taylor expansion for the gravitational function G around today (scale factor a 0 = 1), taking only the first few terms.
We then confront this picture with the observations of type Ia supernovae (SNIa) and baryon acoustic oscillations (BAO). While treating SNIa data, we also take into account the effect that a variation of G might have on SNIa intrinsic luminosity, employing the approach widely used in the cosmology literature [13][14][15][16][17], which assumes the supernovae intrinsic luminosity to be proportional to the Chandrasekhar mass. Additionally, we use the results of a more recent analysis by Ref. [18], which provides an opposite luminosity-G relation, in order to see to what extent our results are affected by the physics of SNIa. On the other hand, we do not modify BAO data under the assumption that they will not be considerably affected by a small variation of the gravitational constant. As a result, we show that the variation of G required to fit these observations is small enough to be compatible with the local constraints ofĠ/G [19].
The paper is organized as follows: in Sec. II, we specify the problem that arises when G evolves as a function of time and detail our general solution. In Sec. III, we obtain the cosmological model, together with a specific G function to test against the data. After detailing the data and the methodology used in Sec. IV, we present the results in Sec. V and conclude in Sec. VI.

II. FIELD EQUATIONS WITH VARYING G
In this work, we take a phenomenological approach to modeling the variation of the gravitational constant. More explicitly, instead of determining how G should behave from some starting assumptions (such as Mach's principle), and deriving our equations from an appropriate Lagrangian, we treat it as a free parameter in Einstein's field equations, whose value is to be determined by observations. This approach of letting the data choose the preferred variation of G allows us to be more general in terms of assumptions beyond general relativity. In other words, rather than offering a new theory of gravity, we stick to general relativity and extend it phenomenologically by allowing G to vary in time.
Let us, then, consider general relativity with a timedependent gravitational constant. We start with modified field equations such that This equation is local; namely, it relates the Einstein tensor at a given spacetime event to the energymomentum tensor and the gravitational constant at the same event. In general, this need not be the case, and we find in the literature various nonlocal theories generalizing Einstein's equations by incorporating retardation effects through a susceptibility function [20]. However, any nonlocal approach behave as being quasilocal when the susceptibility is very stitched around zero. In this situation, the response time of spacetime itself is supposed to be very small compared to other characteristic times (here, a cosmological timescale).
In accordance with the cosmological principle, we consider only a time dependence of G in Eq. (1). In this case, the Bianchi identity implies a nonzero covariant derivative for the stress energy tensor: Therefore, without any other prescription, the energymomentum will no longer be conserved. More precisely, if one assumes the standard form for a cosmological perfect fluid, T µν = (ρ + p)u µ u ν + pg µν , where u µ = (1, 0, 0, 0) is the Hubble flow and p and ρ are the pressure and energy densities, respectively, for the usual matter and radiation, one obtainsρ where H =ȧ/a is the Hubble parameter. The source term in the right-hand side creates a coupling between the energy density and G. For instance, one obtains for nonrelativistic matter (p = 0) which implies a dependence of the matter mass to G such that m ∝ G −1 . This means either that the number of baryons is no longer conserved or that the rest energy of one individual particle depends on G. Either of these options leads to questionable conclusions from a particle physics perspective. Therefore, we aim to preserve the usual conservation relation for D µ T µν = 0. From a Lagrangian perspective, this means that we want to keep the usual √ −g coupling of matter and gravity √ −gL matter , as in Ref. [21]: A simple way of decoupling matter conservation from G and satisfying Eq. (5) is to add a new dynamical component S µν to Einstein's equations. Any symmetric rank 2 tensor can be uniquely decomposed as Au µ u ν + 2q (µ u ν) +Bγ µν +π µν , where A and B are two scalar functions, u µ is a vector field such that u µ u µ = −1, q µ is a vector field transverse to u µ (q µ u µ = 0), γ µν = u µ u ν +g µν is the transverse projector, and π µν is a symmetric transverse tensor such that π µν u µ = 0 and π µν γ µν = 0. Taking u µ as the Hubble flow and assuming isotropy and homogeneity, one has finally q µ = 0 (no energy flux with respect to a Hubble observer) and π µν = 0 (no anisotropic pressure). Therefore, we are left with only two scalar functions, which can be rewritten as with Φ(t) and Ψ(t) being arbitrary functions of time. Then, our modified equations read where 8π is put for convenience. The application of the Bianchi identity to Eq. (7), together with (5), then leads to The new component S µν is clearly not conserved if G depends on time. Here, the only equation containing information is the temporal one (ν = 0) due to spatial symmetry. We can use this equation to obtain a relation between Φ and Ψ. Defining an effective equation of state parameter w = Φ/Ψ, one getṡ This is essentially a generalization of the energy conservation Eq. (5) for a component coupled to the variation of G. Equation (9) shows that this component is also coupled to matter and radiation whenĠ = 0. Conversely, the energy densities for matter and radiation, ρ, do not depend on Φ or G directly, but relate to them through the background relation of H, as intended. In addition, since Eq. (9) is a general expression, it is also valid for a constant G, in which case the right-hand side becomes zero and the equation represents the conservation of energy for models of dark energy fluids uncoupled to matter and radiation.
Using the Friedmann-Lemaître-Roberston-Walker metric with these modified Einstein equations, we directly obtain the cosmological equations in the usual manner: where κ is a constant accounting for the spatial curvature of the Universe.

III. SELECTING SPECIFIC MODELS FOR Φ AND Ψ
In order to use these equations for cosmological analyses, we need to specify a function for Φ. This requires making additional assumptions about the nature of this new component. To do this, let us first solve Eq. (9) for Φ. This can be done by defining an auxiliary function ξ(t) satisfying the equationξ/ξ = 3H(1 + w). With this, Eq. (9) becomes We integrate this function with limits from t = 0 to any time t. However, since we have the singularity at a(0) = 0, ρ approaches infinity at the initial point, there is a possibility that the integration will have a similar behavior at the lower boundary, which would make Φ diverge. With this in mind, we treat the lower boundary with some care: Let F be a primitive ofĠρξ such that Then, the condition for Φ(t)ξ(t) to be finite anywhere is which defines a constant we call C 1 . A similar argument exists for the auxiliary function ξ. First, solvingξ/ξ = 3H(1 + w), we have Integrating by parts gives (17) Assuming that the function w is chosen with care, the expression with the exponential will converge. This leads to the condition for ξ(t) to be finite: This defines a second constant C 2 . However, the actual value of C 2 has no importance, since this factor cancels out in Eq. (13) and does not effect the value of Φ.
One particularly simple and interesting case is when w is a constant. This leads to the simple expression where we can see that the scale factor dependence of 1/ξ resembles that of a matter or radiation density for the appropriate values of w. With these, we can express Φ from Eq. (14) as We now define a critical energy density ρ c , such that H 2 0 = 8πG 0 ρ c /3, and Ω = ρ/ρ c = Ω m a −3 +Ω r a −4 , where the subscript 0 refers to the present time. Dividing by H 2 0 , Eq. (10) becomes For clarity, we also define the parts of Φ in Eq. (14) coupled to radiation and matter separately, such that with C = C 1 /G 0 ρ c , and κ is redefined to include the H 2 0 term. There are some important differences between this equation and the usual Friedmann-Lemaître equation in ΛCDM. First, there is the factor G(a)/G 0 in front of the usual terms for the matter and radiation contributions. This, of course, comes from the direct effect of changing G on the gravitational energies of these components. Skipping κ for the moment, we can see three additional terms that come from the energy component of S µν , which we introduced in order to ensure the conservation of energy. The latter of these terms couple to the energy densities and the G variation, while C is a constant which will be discussed in more detail shortly. As we can see from these, the actual variation of G does not need to be very large, since the integral terms f r and f m can generate the more significant portion of the energy contribution as long as there is a nonzero evolution of G.
Finally, κ is the usual curvature term and can be related to other parameters by evaluating Eq. (22) today (a = 1, or z = 0): At this point, we need additional assumptions in order to determine the full expression for Φ. This can be achieved most straightforwardly by choosing a relation for ξ(a). In order to see what this function represents more explicitly, let us focus on Eq. (22). Since the functions f r and f m are nonzero only if G evolves, they can be attributed to an extra energy contribution in the Friedmann-Lemaître equation arising from a variation of G. On the other hand, the constant C produces in the Friedmann-Lemaître equation an energy contribution C/ξ(a) that exists irrespective of whether G evolves or not. As the scale factor dependence of this term is only through ξ(a), the latter function essentially determines how the energy contribution of C changes with the expansion of the Universe. Now, if S µν represents the reaction of spacetime to the varying gravitational constant, it makes sense to expect that the new component Φ should not change because space expands, but rather because G varies. This means that, if G is constant, Φ should also be unchanging, which leads to the choice ξ(a) = Cst. Going back to Eq. (19), this implies w = −1 or, equivalently, Φ = −Ψ.
With ξ being constant, the C/ξ term in Eq. (22) is the same as the cosmological constant Λ in the standard picture. However, in the present case, we have other terms in Eq. (22) that appear when G varies with time, and we may not need this contribution at all. In order to have Φ represent only the response of spacetime to the evolving G, we keep only the terms that depend on G, which means we choose C = 0. By this choice, we get rid of the first cosmological constant problem, namely, the identification of Λ with vacuum energy and the resulting large discrepancy with quantum field theory estimations, and we can test whether complying with cosmological observations is still possible without a cosmological constant. Therefore, we want to see if it is possible to explain the accelerated expansion with the secondary effect of the variation of G instead of the vacuum energy.
For the cosmological tests we approximate G with a power series expansion around a = 1 to represent the series expansion of an unknown function: This expansion around a = 1 is chosen to capture the behaviour of G accurately around the low-redshift regime we investigate. When the high-redshift regime is considered, such a parametrization could be attached to any extrapolation capturing the behavior of the high-redshift regime without altering the current results at low redshift.
Then, with ξ = Cst and replacing G with Eq. (24), f r and f m become, respectively, The terms are written in a way to facilitate the comparison with the usual Ω r and Ω m terms in Eq. (22). This illustrates that Ω r a −4 and Ω m a −3 will dominate over f r and f m as a gets smaller in the past. Of course, these functions f r and f m become zero when G is constant, i.e., b 1 = b 2 = b 3 = 0. With also C = 0, Eq. (22) reduces to the usual Friedman-Lemaître equations for CDM.
In the rest of the paper, we will assume a flat universe, i.e., κ = 0. This allows us to determine one of the b i parameters of the expansion of G in terms of the others, using Eq. (23), namely, What has been done so far is to formulate a phenomenological variation of the gravitational constant within general relativity in a geometrically consistent way, also preserving the usual energy conservation. In this process, we obtained, in the Friedmann-Lemaître equations, a cosmological constant term together with additional terms that couple to the matter and radiation components. Taking this cosmological constant term to be zero, we are left solely with an energy contribution stemming from the coupling of matter and radiation with a variation of G. In the next sections, we will show that this picture is compatible with low-redshift cosmological probes to a great degree and is also able to conform to local constraints on the evolution of the gravitational constant (Ġ/G).

IV. DATA AND METHODOLOGY
In this section, we describe the cosmological probes and the methodology used in this work. We use the usual low-redshift cosmological probes, type Ia supernovae, and baryon acoustic oscillations, together with the χ 2 minimization method to constrain our model parameters. On the other hand, an investigation of the Cosmic Microwave Background measurements is left for a future work, since the calculation of the effects this variable G would have on the CMB is an involved task that deserves a more through treatment. Consequently, we only consider the low-redshift regime when presenting and discussing our results.

A. Type Ia supernovae
We use the SNIa measurements from the SDSS-II/SNLS3 Joint Light-Curve Analysis (JLA) dataset and its covariance matrix provided by Ref. [22]. We obtain the observed distance modulus following the standardization method given by the authors: In this equation, m, X, and C are the observed magnitude in the B-band rest frame and the shape and color standardization parameters for the different SNIa, respectively, provided in the public dataset. The remaining parameters α, β, and M are nuisance parameters, determined together with the cosmological parameters from the fit to the data. The former two are the same for all SNIa, while the latter is the absolute magnitude in the B-band rest frame. Depending on the stellar mass of the host galaxy (M stellar ), it is given by an additional nuisance parameter ∆M : When the gravitational strength changes due to the variation of the gravitational constant, SNIa intrinsic luminosity should also change due to the G dependence of Chandrasekhar's mass, such that L ∝ M Ch ∝ G −3/2 [13][14][15][16][17]. This modifies the observed distance modulus: If gravity was stronger in the past, for instance, supernovae would be dimmer, so their distances would actually be smaller than they appear. The required relation between the distance modulus and G is directly obtained from the definition of the distance modulus. Considering that the absolute magnitude is related to luminosity via the flux, F ∝ L and M = −2.5 log F (10 pc), the distance modulus is given by However, there are also other approaches in the literature. The authors of Ref. [18] present an opposite relation, by using a semianalytical model to predict the SNIa light curves when the gravitational constant changes with the redshift. Their numerical relation is converted to an approximate expression in Ref. [23] as L ∝ G 1.46 . Since there is no consensus on the precise nature of supernovae physics, as a second case we also use the distance modulus derived from this luminosity-G dependence: in order to check the sensitivity of our calculations to the effect of changing G on SNIa luminosities. We compare the SNIa distance modulus to the predictions of our cosmological models using the standard definition with being the comoving distance for a flat space given in natural units, where c = 1.

B. Baryon acoustic oscillations
In this work, we use a variety of isotropic and anisotropic measurements of baryon acoustic oscillations given in the literature. Isotropic observations measure the quantity D V /r d , where r d is the length of the standard ruler and D V relates to cosmology as Anisotropic observations measure two quantities in the transverse and radial directions: In both cases, there is a degeneracy between H 0 and r d , so we calculate them together as a single parameter. In this work, we assume that the variation of G is small enough to not influence the BAO and use the data without modifications.
Further from the peak value, the likelihoods of BAO observables diverge from a Gaussian distribution. In order to take this into account and be more conservative with our estimations, we follow the recipe in Ref. [30] and replace the standard ∆χ 2 G = −2 ln L G likelihood expression for a Gaussian distribution with where S/N stands for the detection significance, in units of σ. We consider a detection significance of 2.4σ for 6dFGS, 2σ for SDSS-MGS, 9σ for BOSS DR12, 4σ for eBOSS DR14, and 5σ for the Ly-α forest.

C. Determination of the parameter constraints
Using a frequentist approach, we obtain the best-fit values for the parameters by minimizing the expression where r pred and r obs are the vectors that include the model prediction and the observations at each redshift, respectively, and C is the covariance matrix of the observations. We add the χ 2 values corresponding to each probe with the assumption that they are statistically independent. To minimize this function, we use PYTHON's iminuit module 1 , an implementation of SEAL Minuit, developed at CERN [31].  Table I shows the χ 2 values of our varying-G model and of the standard flat ΛCDM model. The best-fit values of the parameters are also shown. We remind that b 1 is obtained through Eq. (27). In order to obtain the uncertainty of b 1 we generate 10 6 random sets of parameters from an N -dimensional Gaussian centered at the best fit and with the corresponding covariance matrix. For each one of these sets we derive the value of b 1 and compute the uncertainty from the standard deviation.

V. RESULTS AND DISCUSSION
The reconstruction of the G(z) function is shown in Fig. 1, where the red line is drawn using the best-fit values and the gray lines show sample lines with ∆χ 2 < 1. Again, these lines are obtained by generating random sets of parameter values from an N -dimensional Gaussian centered at the best fit and with the corresponding covariance matrix. The top panel in Fig. 1 shows the G(a) function itself, while the second panel shows the first derivative with respect to the scale factor. In both plots the functions have been normalized with respect to G 0 , the present-day value of G. Constructed in a similar way, Fig. 2 shows the ratio between the different terms in Eq. (22) that drive the expansion at the considered epoch, −Ω m f m and Ω m a −3 G/G 0 . This graph shows that the former term starts to dominate at the late stages, causing the accelerated expansion.
It is clear from Table I that the varying G model has almost the same χ 2 value as the flat ΛCDM model we use for comparison. Therefore, this model is indeed capable of explaining SNIa and BAO observations without the addition of a cosmological constant. While it might seem surprising that such a small variation of the gravitational constant can result in a contribution large enough to supply most of the energy in the Universe, it is apparent from Eq. (22) and Fig. 2 that the major effect of G on the energy balance does not come from the G/G 0 term. This energy instead comes from the last term on the right-hand side of this equation, Ω m f m , which appears as a secondary effect of a time evolution of G. This also explains why the first-order term in the G function can be so small: The contribution of b 1 in Eq. (26) is not any larger compared to the other terms, b 2 and b 3 , inside the brackets. One of the major challenges when predicting a variation of the gravitational constant comes from the highly tight constraints measured from local observations. The review in Ref. [19] compiles the different observations and givesĠ which is obtained from the lunar laser ranging experiment [32], as the tightest constraint on the variation of G.
With the series expansion given by Eq. (24), the quantitẏ can be simply evaluated using the b 1 parameter: where the prime denotes the derivative with respect to the scale factor and H 0 is around 67 − 76 × 10 −12 yr −1 , according to various recent measurements [33]. Looking at Table I, we can see that our results are compatible with this constraint at the one-sigma level. Moreover, since our model is compatible with b 1 = 0 at one sigma, even lower values ofĠ G at z ≈ 0 are consistent with our results.
Bounds on the variation of G for earlier times also exist, such as stellar observations for low redshifts and big bang nucleosynthesis (BBN) and CMB measurements for much higher redshifts. In the former case, the constraints are much less limiting than the Solar System observations [19], and they usually assume a monotonic G evolution (as in Ref. [34], for instance), which is not the case in our calculations. In the case of BBN and CMB, the limits on G evolution concern much higher redshifts than the ones our analysis considers, and, therefore, they are outside the scope of this work.
One new way of falsifying alternative gravity theories is provided by the recent observation of the GW170817 neutron star merger, as it has shown with great accuracy the gravitational wave propagation speed to be equal to the speed of light [35]. Here, we will not provide a full mathematical demonstration but only point out that, since we do not change the geometry of spacetime, the gravitational wave propagation speed remains the same as in standard general relativity. What our approach rather does is analogous to adding background source terms, which does not affect the propagation speed. On the other hand, as shown in Ref. [36], modifications of the gravitational constant may cause the standard siren luminosity distance to differ from its electromagnetic counterpart. As the capabilities of gravitational wave observatories increase, in the future this may potentially be used to probe the history of a possible G evolution, but since the only available observation, GW170817, is from a very low redshift, it does not put an additional constraint to our model at the present.
Turning to the other values in Table I, we see that, for the varying-G model, the best-fit value of Ω r is consistent with zero, implying a very small ρ r at the present epoch, as expected. On the other hand, we obtain Ω m = 0.284± 0.017, which is similar to the usual value for ΛCDM. Therefore, we see that our model does not change the matter content drastically. The best-fit value of the H 0 r d parameter is also in agreement with the standard results.
When we repeat our analysis with the SNIa intrinsic luminosity-G relation from Ref. [18], as discussed in Sec. IV A, we find a χ 2 value slightly larger than the results discussed so far, χ 2 = 698.48, but still perfectly compatible with the value found for ΛCDM, χ 2 = 698.05. The values of the series expansion parameters for G also change somewhat, with b 1 = 0.22 ± 0.14 , Most notably, b 1 becomes compatible with zero at two sigma instead of one. On the whole, the differences are not too drastic considering that the two supernovae luminosity models are completely opposite to each other, which leads to the conclusion that our approach does not depend heavily on the exact nature of supernovae physics. This is not surprising, since the absolute variation we predict of the gravitational constant is small in both cases, and its effect on supernovae should likewise be slight.

VI. CONCLUSION
In this work, we show that, when a phenomenological variation of the gravitational constant is allowed, general relativity can explain the low-redshift accelerated expansion of the Universe without a cosmological constant. When G is taken as a time-dependent function in Einstein's field equations, the enforcement of the Bianchi identity and the usual energy conservation causes a new term to appear. This term represents the coupling between the variation of G and the energy density of matter and radiation, and we determine further properties of it by requiring that it can be interpreted as a reaction of spacetime to the variation of G. We then test the resulting model with SNIa and BAO data.
The comparison with observational data shows that this extra term can cause the late-time accelerated expansion with a deviation less than 10% from the current value of G in the considered redshift range. We show that the required corrections to the gravitational constant are essentially on the second and third order, while the first order turns out to be small, consistent with zero within one standard deviation. This implies that the most strict bounds on a possible variation of the gravitational constant from local observations are also satisfied.
From these results, we observe that in our varying-G model the main driving force behind the accelerated expansion is not the direct effect of the gravitational constant itself but the influence of the additional term that appears because of the time dependence of G. While this term is dominated by the contributions of matter and radiation for higher redshifts, it starts to supply most of the energy in the Universe during the late stages and thereby facilitates the acceleration. As a result, we see that this model can explain the late-stage accelerated expansion of the Universe without a cosmological constant and without requiring a large impact on the small-scale gravitational processes.
Finally, let us mention that varying-G models also have a broader potential in explaining other cosmological tensions. As an example, it has been shown recently that a cosmological Brans-Dicke model with a cosmological constant can alleviate the tension on the Hubble constant H 0 [37]. While this discussion is outside the scope of this paper, it shows the potential of considering varying-G models as interesting alternatives to the standard general relativity.