Dynamical analysis of the covarying coupling constants in scalar-tensor gravity

A scalar-tensor theory of gravity is considered wherein the gravitational coupling $G$ and the speed of light $c$ are admitted as space-time functions and combine to form the definition of the scalar field $\phi$. The varying $c$ participates in the definition of the variation of the matter part of the action; it is related to the effective stress-energy tensor which is a result of the requirement of symmetry under general coordinate transformations. The effect of the cosmological coupling $\Lambda$ is accommodated within a possible behaviour of $\phi$. We analyze the dynamics of $\phi$ in the phase space, thereby showing the existence of an attractor point for reasonable hypotheses on the potential $V(\phi)$ and no particular assumption on the Hubble function. The phase space analysis is performed both with the linear stability theory and via the more general Lyapunov's method. Either method lead to the conclusion that the condition $\dot{G}/G=\sigma\left(\dot{c}/c\right)$ where $\sigma=3$ must hold for the rest of cosmic evolution after the system gets to the globally asymptotically stable fixed point and the dynamics of $\phi$ ceases. This result provides a physical foundation for the phenomenological model admitting $\left(G/G_{0}\right)=\left(c/c_{0}\right)^{3}$ used recently to interpret cosmological and astrophysical data. The thus co-varying couplings $G$ and $c$ impact the cosmic evolution after the dynamical system settles to equilibrium. This impact is investigated by constructing the generalized continuity equation in our scalar-tensor model and considering two possible regimes for the varying speed of light -- decreasing $c(a)$ and increasing $c(a)$ -- while solving our modified Friedmann equations. The solutions to the latter equations make room for radiation- and matter-dominated eras that progress to a dark-energy-type of accelerated expansion.


Introduction
The possibility of variation of fundamental physical constants has been long explored as a means of solving some issues in the description of nature, especially in astrophysics and cosmology. In spite of being a controversial idea, it has deep connection with scalar-tensor theories, themselves a popular candidate for viable modified theories of gravity. This paper explores the possible interweaved variation of the couplings G and c (while also accommodating the effects of Λ) in a Brans-Dicke-like model. Our motivations will be stated and the construction of our setup will begin during this introductory section, which starts off with a brief perspective on the subject of varying fundamental constants.
The two most studied constants for their potential variation on cosmological time scales are the gravitational constant G and the fine structure constant α. While G is a dimensionful constant, α is a dimensionless constant. Uzan [1,2] has critically reviewed the subject of the variation of fundamental constants and discussed it extensively. Some authors (Ellis and Uzan [3], Duff [4,5]) strongly argue against dimensionful varying fundamental constants. In an attempt to avoid this criticism we normalize the time derivative of the physical coupling by the value of the associated quantity.
Variable fundamental constant theories can be traced back to at least the late nineteenth century, e.g., [6], and the first half of the last century, e.g., [7,8]. However, they gained importance after Dirac [9,10] suggested potential variation of the gravitational constant G based on his large number hypothesis. Many observational methods were suggested to constrain an eventual variation of G including neutron star masses and ages [11]; CMB anisotropies, e.g., [12]; big-bang nucleosynthesis abundances, e.g., [13]; asteroseismology [14]; lunar laser ranging, e.g., [15]; the evolution of planetary orbits, e.g., [16]; binary pulsars, e.g., [17]; supernovae type-Ia (SNeIa) luminosity evolution, e.g., [18]; and gravitational-wave observations on binary neutron stars [19]. Almost all of them have resulted in the constraints onĠ/G well below that predicted by Dirac. There has been significant development on the theoretical side as well. Building on the work of Jordan [20], Brans and Dicke [21] developed a scalar-tensor theory of gravitation wherein 1/G was raised to the status of a scalar field which could vary spatially and temporally. Brans-Dicke (BD) theory may be seen as one representative of the class of scalar-tensor theories in which gravitation manifests itself by both the metric tensor and a scalar field of geometrical nature [22]. A more general example within this class is the scalar-tensor theory in Lyra manifold [23,24]. One of the reasons why scalar-tensor theories of gravity raise so much interest (see [25] and [22] and references therein) is their well know equivalence with some modified gravity theories [26], cf. e.g., [27,28]. Among these possible modifications of gravity we mention the theories f (R) [29,30,31], f (R, ∇ n R) [32], f (R, n R) [33], f (T ) [34,35], which have been proposed as attempts to deal with the GR quantization problem [36,37], to realize inflationary models [38], to address the dark energy problem [39]. Refs. [40,41,42] studied the phenomenology of extended Jordar-Brans-Dicke theories and their implication for cosmological datasets. Moreover, significant astrophysical consequences can be found in the literature within the realm of modified theories of gravity including those in Refs. [43,44,45,46].
The constancy of the speed of light is the foundation of special and general relativity theories, and arguably its variation is the most contentious issue in physics. However, even Einstein considered its possible variation [47]. There are several theories of the variation of the speed light, e.g., those by Dicke [48], Petit [49], and Moffatt [50,51]. Some of these proposals break Lorentz invariance, e.g., [52,53,54], others produce locally invariant theories [55,56].
Attempts have also been made to consider the simultaneous variation of two or more constants, e.g., Ref. [57] considered varying G and Λ; variation of G/c 2 was studied in [58]; G and α were the changing concomitantly in Ref. [59]; co-varying c, G and Λ were the subject of Refs. [60,61,62]; the set {c, G, α} was allowed to vary in [63]; and, finally, Refs. [64,65] suggested that all the couplings G, c, Λ, , and k B must co-vary.
After this contextualization on the subject of varying physical constants, let us concentrate on the topic to be explored here. Our focus in this paper is on the variation of c and G and the interrelationship of their variation. We will see that this interrelationship yields naturally to an effective cosmological constant. Gupta [66] introduced the general ansatzĠ inspired by the general constraint appearing in the proposal by Costa et al. [61]. The quantity σ is a constant parameter. The value σ = 3 is strongly favored by several phenomenological applications of Eq. (1) in cosmology [64,66,67,68] and astrophysics [65,69,70,71]. This preferred value for σ is corroborated by other authors, cf. e.g., Ref. [63]. For this reason, we set off to investigate the fundamental reasons that might be underlying this fact. We do that by assuming a modification in the standard theory of gravity: General Relativity (GR). This is justified from the fact that GR does not allow for varying physical constants. In fact, the couplings G, c, and Λ showing up in Einstein field equations do not vary in time or space. Recall the definition of the Einstein tensor: G µν = R µν − 1 2 g µν R, where R µν is the Ricci tensor and R is the curvature scalar. We follow the metric signature and Riemann tensor definition of Ref. [72], i.e. the Minkowski metric in Cartesian coordinates reads (η µν ) = diag (−1, +1, +1, +1), R µν = R ρ µρν with the curvature tensor calculated by R ρ σµν = ∂ µ Γ ρ νσ + Γ ρ µλ Γ λ νσ − (µ ↔ ν), R = g µν R µν , and Γ ρ µν = 1 2 g ρσ (∂ ν g σµ + ∂ µ g νσ − ∂ σ g µν ) stands for the Christoffel symbols built from the metric tensor g µν and its first-order derivatives. Covariant derivatives are built from the Christoffel, e.g., ∇ µ V ν = ∂ µ V ν + Γ ν µρ V ρ for an arbitrary contravariant vector V ν .
Our generalization of GR stems from an action integral that is inspired in Einstein-Hilbert action, but with the couplings {G, c, Λ} being spacetime functions. At this point, we emphasize the natural appearance of the couplings {G, c, Λ} explicitly in the action integral that is supposed to describe gravity. We feel the urge to underscore this fact for two reasons. The first is: in the ordinary approach of GR, both G and c are only multiplicative constants of the kernel (R − 2Λ) and, as such, can be taken outside the integral with impunity. In fact, people even go further to use units where c = 1, the so called natural units. These practices are very dangerous here because both G and c are space-time dependent functions; accordingly, they are directly affected by the variation process soon to be carried out. Second, the fact that (3) displays the set {G, c, Λ} justifies why we consider only these three couplings, instead of bringing about other fundamental constants as well, such as the Planck constant and the Boltzmann constant k B . Others works considered this enlarged set from a phenomenological perspective-see e.g., Refs. [64,65,66,67,68,71,73]. As usual, the four-volume in the action S g is d 4 x = dx 0 d 3 x. It pays out to stress that the time coordinate x 0 = ct encompasses the speed of light c, which is allowed to vary in our context. Therefore, it is paramount to work with x 0 throughout. In this way, covariance is guaranteed; caveats like having to deal the opened x 0 -differential dx 0 = d (ct) = cdt + tdc are avoided; and the speed of light will not appear in the metric components g µν explicitly.
Consequently, the determinant of the metric tensor g = g (x µ ) will depend explicitly on x 0 , but not on c. The time component x 0 hides c in a convenient form: if c changes, then a time reparameterization t → t compensates for this fact in each reference frame. This strategy is somewhat different from the procedure of the c-flation framework developed in Ref. [61]. Here, x 0 has dimensions of length, which explains the power cube for the speed of light in (3): otherwise the action would not bear the correct dimension of angular momentum (energy × time).
Gupta's proposal in Eq. (1) indicates that the couplings G and c are intertwined. In fact, it could be recast as: is a scalar field comprising the gravitational coupling G and the speed of light c. The dot on top of the quantity means a derivative with respect to x 0 , for instanceφ = dφ dx 0 . The simplest way to satisfy (1) is to admit thatφ = −φ Ġ G − 3ċ c = 0 due to constant couplings, in which caseĠ = 0 andċ = 0. This is the way of GR, but it is not the most general possibility and not the one we will be adopting here. An earlier attempt to relax the constancy of fundamental couplings in the context of a theory for gravity was the Brans-Dicke (BD) scalar-tensor theory [21]. Aiming to fully realize the Machian principle in a geometrical theory of gravitation, Brans and Dicke proposed a varying gravitational coupling through the scalar field [22] φ BD = 1 G .
Our proposal is to extend the scope of BD field to include the spacetime coupling c as part of its very definition. The presence of c changes the nature of the scalar field φ as previously envisioned by Brans and Dicke. 1 This is true even from a dimensional point of view. In the Brans-Dicke picture, φ BD is measured in units of (mass) × (time) 2 / (length) 3 . In our scalar-tensor gravity, φ has dimensions of (mass) / (time). The fact that a varying c enters the scalar field φ in Eq. (5) may be a cause of concern to some. A common criticism upon varying speed of light proposals is that the Maxwell electrodynamics should not be tempered with. This is a wise thought, however one that can be bypassed by arguments such as that by Ellis and Uzan in Ref. [3]. There could be different types of speed of light, viz. the spacetime speed of light c ST related to Lorentz invariance and causality and the electromagnetic speed of light c EM built from the couplings in electromagnetism: the (vacuum) electric permittivity 0 and the magnetic permeability µ 0 . Here we interpret c as the causality coupling c ST and ensure that Lorentz invariance (of Maxwell's theory) is not broken locally by convenient reparameterizations of t at each time slice of the space-time manifold. Moreover, our approach deals with cosmology which makes room for a variation of c with respect to the cosmic time in such a way that c ST = c EM = c 0 nowadays.
The main goal of this paper is to investigate if the dynamics of φ can naturally leadĠ/G = 3 (ċ/c) as an attractor solution. This would give the fundamental reason for why the value σ = 3 is preferred in several apparently different contexts, such as cosmology [66], solar astrophysics [71], and solar system kinematics [70]. After the eventual relation between G x 0 and c x 0 is established, a natural question poses itself: "What are the consequences of G = G (c) for the cosmic evolution?". Our secondary goal in this paper is to answer this question by investigating cosmological solutions attributable to different epochs in the Universe's thermal history (e.g., radiation dominated era, dust-matter era, and recent accelerated regime). We will build our scalar-tensor model in a metric compatible, four-dimensional, torsion-free, Riemannian manifold. Approaches considering richer geometrical structure [24], higher-order derivatives of curvature-based objects [28,32] or torsion-based invariants [35] are possible, but these cases shall be explored elsewhere.
The remainder of the paper is organized as follows. In Section 2 we state the action of our scalar-tensor model, derive the equations of motion for the tensor field g µν and the scalar field φ, which together build the gravitational field. The field equations are specified for the FLRW metric [74], since our main interest will be to study the dynamics of φ in the cosmological context. This study is carried out in Section 3, where it is shown that an attractor solution leading to σ = 3 is attainable under reasonable hypotheses for the potential V (φ) and the dominant matterenergy content during the evolution of φ. 2 This finding comes from the linear stability analysis of our dynamical system-performed in Subsection 3.2; the result is confirmed and strengthened in Subsection 3.3 where we use the general Lyapunov's method to show that the critical point φ eq leading to σ = 3 is a globally asymptotically stable fixed point. Subsection 3.4 gives some details about the system's attractor point and emphasizes the consequences for our modified gravity. The cosmic evolution after φ reaches equilibrium is analyzed in Section 4 wherein it is shown that our scalar-tensor model accommodates radiation-and matter-dominated eras evolving naturally to a de Sitter-type accelerated expansion, possibly recognized as dark energy. This is true for the two parameterizations we adopted for the function c = c x 0 , the first assuming a speed of light that decreases as the universe expands while the second admits c increasing with the scale factor. Our final comments appear in Section 5.
2 Modified gravity with φ = c 3 /G The motivations in the previous section lead us to describe gravity from a precisely specified form for the action integral. In fact, this form is hinted by the coefficient of the curvature scalar R in Eq. (7), namely c 3 /G , which was defined as the scalar field φ. Accordingly, we introduce where ω is a dimensionless constant of order one 3 and φ is a scalar field. The constant ω appears as the coefficient of the kinetic term for the field φ; the denominator of this term includes φ for dimensionality consistency. The cosmological coupling Λ was not written explicitly in (7) because it can be included in the potential V (φ). In fact, we will show that a constant Λ can be realized in our framework. The reader can easily check that all terms in (7) were crafted in order to give S the correct dimension of (energy) × (time). This includes inserting the factor (1/c) explicitly inside the matter term integral containing the matter Lagrangian density L m , itself carrying dimensions of energy density. The action (7) is formally the same as Brans-Dicke action. While Brans-Dicke scalar field is usually interpreted as an effective G −1 (see Ref. [22]), our φ includes both G and c. Notice, moreover, that (7) is not the same as the actions in the Varying Speed of Light (VSL) proposals-such as those in Refs. [50,52,53]-for reasons ranging from different modelling to the basic fact that variations of the gravitational coupling are disregarded therein. Additionally, it should be pointed out that (7) does not correspond to the same treatment given in Ref. [61]; indeed, the covariant c-flation framework considers G and c as independent fields while here these varying couplings enter the form of a single scalar field; on top of that, the very form of the action integral is different both in its gravitational part and in the kinetic term for the scalar field(s).
Variation of the action in Eq. (7) with respect to the metric g µν yields the BD-like Field Equations for g µν : 4 where is the stress-energy tensor of ordinary matter. Notice that T µν in (9) has the regular dimensions of energy density. However, it should be underscored that the variation of the matter action with respect to g µν is written as in the context of our scalar-tensor model. Here the factor (1/c) must participate the form of δS m within the integral sign. This is not the case in non-varying speed of light scenarios (including the pure BD theory). The factor (1/c) accompanying L m (or T µν ) brings about one of the interesting feature in our scalar-tensor model: large-valued varying-c setups suppress the ordinary matter contribution to the action; in which case, the dynamics is mostly determined by φ through its kinetic and potential terms. In Eq. (9), no dependence of L m on the field φ is assumed. However, it should be noted that we can not evade a non-trivial coupling with the (varying) speed of light within the functional derivative entering the definition of T µν . This feature will be dealt with later-see Subsection 4.1; it is related to the definition of an effective stress-energy tensor including c. According to Noether's theorem [76] to every symmetry corresponds a conserved current; in theories of gravity, the conserved current is the own stressenergy tensor and the symmetry is the invariance under general coordinate transformations (realized in practice by infinitesimal translations involving Killing vectors [77,78]). In the paper [61] this delicate point surrounding the underlying symmetry and the associated energy-momentum tensor appears in a different form: the covariant c-flation approach introduces T µν directly in the action via a Lagrange multiplier. Variation of S with respect to the scalar field φ leads to the Equation of Motion (EOM) for the scalar part of our scalar-tensor model: This EOM tells that "φ is a dynamical field; it changes in space and time." As a consequence of the definition of φ in terms of G and c, the latter conclusion yields the co-varying character of the gravitational coupling and the causality coupling. We will show bellow that the simultaneous variation of G and c occurs in such a way that the conditionĠ/G = 3 (ċ/c) is satisfied after the system evolves to the equilibrium stable point in the phase space. In the following we specify our BD-like theory (in the Jordan frame) for cosmology-see e.g., [22]. The main hypothesis is that φ depends only on the cosmic "time" x 0 = ct due to the requirements of homogeneity and isotropy. These requirements demand the FLRW line element: In this set of coordinates the scale factor a x 0 bears the dimension of length. Conversely, the comoving coordinates {r, θ, ϕ} are dimensionless quantities. The curvature of the space section is determined by the parameter k = −1, 0, +1 as usual [74]. The Hubble function H is defined as and has the dimension of inverse length.
The matter content appearing in both field equations (8) and (12) shall be modelled by the perfect fluid stressenergy tensor, where p is the pressure and ε is the energy density (both with dimensions of energy per unit volume). Due to (15), the 00-component of the g µν field equations (8) reads: 5 This is the (first) Friedmann equation of our scalar-tensor cosmology with covarying G and c through the field φ.
The second Friedmann equation (or acceleration equation) is calculated by taking the trace of Eq. (8), using the curvature scalar built from the metric in (13), utilizing the trace of the stress energy tensor, and the Eqs. (12) and (16). It reads:Ḣ Moreover, in the cosmological context, Eq. (12) reduces to: This equation determines the dynamics of φ. It depends on H = H x 0 , which means that (16), (17) and (18) should be solved simultaneously. The latter also depends on the matter content through the factor (ε − 3p), which demands us to specify the nature of the perfect fluid at play in the cosmological era under scrutiny. Finally, Eq. (18) exhibits a dependence on the potential V (φ) and its derivative w.r.t. the field φ. Ultimately, we will have to say what is the form that we expect for potential φ which, in turn, is written to include the physical couplings G and c. This challenge will be addressed next: In the following section we perform the phase space analysis of our model. With regards to dynamical systems approach in cosmology, we refer the reader to the interesting papers [79,80,81,82,83,84,85,86,87].

The dynamics of φ
The potential V (φ) is assumed to be analytical but otherwise completely general. It will be taken as dominant over the contribution of the term involving the matter Lagrangian; more specifically: This condition may be satisfied either for φ c dc dφ We shall consider the first terms of the expansion of V (φ) in a power series: where V i (i = 0, 1, 2) are the first three constant coefficients in the truncated series. This type of expansion up to second order in the power series is customary of the potential description around a local minimum of the function. It is a usual practice in classical mechanics and quantum mechanics to approximate the potential in the vicinity of a local minimum for that of a harmonic oscillator-e.g. Ref. [88]. The usefulness of a power series expansion of the potential up to its first-order terms is clearly noticeable in the inflationary context-e.g., [89]-with the quadratic term being a useful example to illustrate the mechanism of a slowly-rolling scalar field [90]. The power series expansion of V has also applications in classical field theories, such as in Ref. [91]. A power series expansion up to second-order can be found in Ref. [89]; the interpretation of the different terms that participate in a truncated series such as that in Eq. (20) are discussed in [22,92], for instance. The form (20) encodes a few cases of special interest: This realizes a constant potential that could be rescaled by setting V 0 = 0.
2. For V 0 = V 2 = 0 and V 1 = 2Λ one recovers the ordinary cosmological constant (Λ = constant) already present in GR. The enticing feature here is that V (φ) = 2Λφ would exhibit all the fundamental couplings in gravity {G, c, Λ} when we substitute φ = c 3 /G. (Nevertheless, recall that Λ would still be a genuine constant here since it defines V 1 .) 3. If V 0 = V 1 = 0 and V 2 = 1 2 m 2 φ , the potential accounts for a massive term in φ's field equation. In fact, in this case V (φ) = 1 2 m 2 φ φ and d 2 V dφ 2 = m 2 φ thus corroborating the interpretation. The delicate part is to interpret the role of this mass in φ = c 3 /G. It could be conjectured that the source of φ's mass is either a massive photon or a massive graviton.
In view of the above features, we consider the V (φ) in Eq. (20) adequate enough to allow for a sufficiently general dynamical analysis of our scalar-tensor theory of covarying physical couplings in gravity.
Substitution of Eq. (20) into (18) yields a form for the field equation of φ that does not depend on V 2 . This is not the final possible simplification to this equation. The next step further is to consider periods in the cosmic history when or, more specifically, The above constraints not only guarantees a great simplification of Eq. (18) but it is also very reasonable. Indeed, Eq. (22) is consistent with a radiation filled universe-described by the equation of state p = ε/3. By abiding to Eq. (22) we assume that the dynamics of the field φ takes place in this era.
In standard cosmology, the radiation era spans from the period after inflation and reheating all the way until the matter dominated phase. It includes major events such as electroweak unification, quark-hadron transition, neutrino decoupling, and nucleosynthesis [93]. It is a key time period in the cosmic thermal history during which the dynamics of φ would have taken place.
In any case, the essential ingredient to the remainder of our analysis is that the condition in Eq. (21) is satisfied during all the dynamics of φ. This means that the combination φ dV dφ − 2V built from the potential V (φ) should dominate the matter-energy content entering the play via 1 c T µν . This configuration could be easily realized within varying-c frameworks during regimes where (c/c 0 ) 1; for example, in the very early universe, before the phase transition of the decaying c x 0 admitted by some VSL models [1].
Plugging (22) into (18) results in:φ A solution to this equation leads to φ = φ x 0 , which constrains the dynamics of c 3 /G . The physical solution to (23) must avoid the point φ = 0 since the notion of a vanishing speed of light or of an infinitely large gravitational coupling is meaningless. In order to realize this physical requirement, we displace the origin of the scalar field φ by defining: where we admit V 0 = 0 and V 1 finite. In terms of this new field variable, Eq. (23) reads: The solution to this equation is not so easy because H = H x 0 and is given by Eq. (16). 6

Elements of stability theory for dynamical systems
We shall study Eq. (25) from the point of view of dynamical systems. For this reason, it is worth briefly reviewing the elements of standard phase space analysis most commonly used in modern gravity. This is done in the following subsections, whose discussions are based on the papers [79], [82], and references therein. Check also Ref. [94].

Linear stability theory
Let us suppose a physical system can be described by n variables y = y 1 , y 2 , . . . , y n . Each of these variables are real-valued functions, so that y ∈ M ⊂ R n . Now consider that the dynamics of this physical system can be described by a set of n differential equations of first order:ẏ where dot indicates derivative with respect to a parameter representing the evolution of this system, say a time coordinate (x 0 ). Here, f (y) represents a set of n analytical functions (which can be nonlinear) of the variables y. Whenever there is no explicit dependence of f with respect to x 0 , the set of equations is called an autonomous system.
If there exists a point y 0 ∈ M such that f (y 0 ) = 0, theṅ If there is a time x 0 (f ) where y x 0 (f ) = y 0 , then the dynamics of the system ceases and y x 0 = y 0 for any value of x 0 > x 0 (f ) . This point is called a fixed point, critical point or equilibrium point. Whenever a system has a fixed point, it is interesting to analyze the trajectories in a neighbourhood U ⊂ M of y 0 . Roughly speaking, if the trajectories in U diverge from the fixed point, it is considered unstable; if the trajectories converge to y 0 , then it is stable.
The simplest stability analysis that can be carried out for such system consists in considering the equations in a small region U where f (y) can be expanded in a Taylor series about y 0 and approximated to first order. The advantage of this method is that it allows for an analytical solution for the approximated system (the drawback is that the approximation by a Taylor series is not always a good approximation for a nonlinear system). In U , we have:ẏ A change of variable, z = y −y 0 , can always be proposed so that the equilibrium point is the origin of the coordinate system in the new variable:ż or in a matrix form: M is called the Jacobian matrix or stability matrix. As long as det M = 0, we can look for a new transformation of variables z → Z where the transformed Jacobian matrix becomes diagonal: the elements of the diagonal are actually the eigenvalues, λ i , of the Jacobian matrix. In this case, the set of differential equations will have the following solution: If the real part of at least one eigenvalue is positive, then the trajectories will diverge from the fixed point and the system is unstable; if the real part of all eigenvalues are negative, then the trajectories converge to the equilibrium point and the system is stable. If the real part of one of the eigenvalues is zero, then the system can be either stable or unstable-in this case, each system has to be considered separately.
The eigenvalues of the Jacobian matrix are called Lyapunov coefficients and the stability analysis essentially consists in the study of their real parts. Notice that the Jacobian matrix is composed of constant coefficients once f is not explicitly x 0 -dependent. The stability that is characterized by the Lyapunov coefficients is restricted to a region where the linear approximation for f is valid.
The above considerations shall be employed in Subsection 3.2.
When the components of the stability matrix are not constant or when the linear approximation of the Taylor series is not a good approximation, other methods have to be applied. For instance, a powerful method to establish the stability of dynamical system when the linear stability analysis fails is the Lyapunov's method presented in the next section.

Lyapunov's method 7
This method is more powerful and general than the one in Subsection 3.2 in the sense that it does not rely on linear stability: It has the capability of investigating both local and global instability. In addition, Lyapunov's method applies to non-hyperbolic points too, something the linear stability theory fails to accomplish. A critical point (fixed point) y = y 0 ∈ M ⊂ R n of the dynamical system (26), is classified as a hyperbolic point if all of the eigenvalues of the stability matrix have non-zero real part. Otherwise, it is said to be non-hyperbolic. Lyapunov's method for checking stability of a dynamical system consists in finding a Lyapunov function V (y), such that: There is no known systematic way of deriving the Lyapunov function V (y). If a Lyapunov function exists, the Lyapunov Stability Theorem establishes that the critical point y 0 is a stable fixed point if the requirementV ≤ 0 is fulfilled. The critical point y 0 is an asymptotically stable fixed point if there is a Lyapunov function V (y) satisfying the criteriumV < 0. In addition, if lim y →∞ V (y) → ∞, ∀y ∈ U , then the critical point y 0 is classified as globally stable or globally asymptotically stable, respectively. See Ref. [79].
Eq. (33) was used in the last step. The relation in Eq. (34) will be important in Subsection 3.3 below.

Phase-space analysis with H (x 0 ) ≈ H * via the linear stability theory
We avoid the difficulty mentioned below Eq. (25) by assuming, in this subsection, that the Hubble function is constant during the time interval in which the dynamics of the field φ occurs. In other words, It should be emphasized that this condition is not as restrictive as it may seem. It does not mean that H (t) is always constant, but only during the time interval of ϕ's evolution. We are not demanding the universe to be stationary (withȧ = constant); we are essentially assuming that H evolves slowly during the time it takes ϕ to reach an eventual attractor point. Again, this is a first-order phase-transition scenario common to several models of VSL [1].
In order to show that the picture described in the previous paragraph is achievable, let us analyze the dynamics of ϕ in the phase space resorting to the linear stability theory [79]. Let be the momentum associated to ϕ. Then, Eqs. (25) and (36) become the following coupled system of equations: This pair of differential equations has the structure of an autonomous system [94]. The dynamical system (37) is already linear; consequently, the linear stability theory is not just an approximation around the critical point(s) but rather an exact description, completely adequate to the case in hand. The result stemming from this analysis should not suffer from any pathologies, such as those mentioned e.g. in Ref. [83] -see also [95]. A quick inspection of (37) revels the equilibrium point: The next step is to analyze the stability of the fixed point according to the Lyapunov criterium [94]. The Lyapunov coefficients λ i are the eigenvalues of the Jacobian matrix (stability matrix) related to the autonomous system [79]. They are found by solving the characteristic equation associated to this matrix, and read: The autonomous system can be considered stable if the real part of both eigenvalues satisfy [79] Re [λ ± ] 0.
The analysis of the direction fields in the phase space allows us to double check if stability is attained (in an enlarged region surrounding the equilibrium point [82]) and if it is consistent with the Lyapunov criterium. Here we assume that i.e. our universe is expanding. We split our analysis to cover different possible ranges of values for the parameters V 1 and ω.

Case 1: Negative valued V 1
Here we take V 1 < 0. There are two subcases to be analyzed according to the values of the parameter ω. These subcases are considered separately next. In the instance V 1 < 0 and ω < − 3 2 , we may write: This condition is important to decide on the behaviour of the square-root showing up in Eq. (39). As a consequence, We conclude that λ + is a positive real number, so that Re [λ + ] > 0 is always satisfied, and the associated equilibrium point is unstable according to Eq. (40). Since the analysis of λ + already reveals instability, there is no need to proceed to the study of λ − . The direction fields for the autonomous system (37) in the case where V 1 < 0 and ω < − 3 2 are shown in Fig. 1. Fig. 1 clearly shows that the trajectories diverge from the equilibrium point (ϕ 0 , p ϕ,0 ) = (0, 0), thus indicating its instability in the context of the subcase under scrutiny here. This conclusion will guide us towards constraining the interval of values allowed for the parameters V 1 and ω if a reasonable physical interpretation for the dynamics of φ = c 3 /G is to be established. We will decide on that later on, after we conclude our stability analysis. Here we need both V 1 < 0 and ω > − 3 2 . Under this requirement, the eigenvalues in (39) can be expressed as: The relevant part for deciding the nature (real or complex) of the eigenvalues λ ± is the second term within the square root. If it is greater than 1, then and both eigenvalues λ ± become complex numbers with negative real parts. Therefore, the equilibrium point is stable, cf. Eq. (40).
On the other hand, if the square root in (43) is a real number smaller than one, and both eigenvalues λ ± are real and negative. This automatically satisfy the criterium Re [λ ± ] 0 for stability, leading us to conclude, again, that the equilibrium point (ϕ 0 , p ϕ,0 ) = (0, 0), is stable.
The direction fields in the plane (ϕ, p ϕ ) are analyzed in Fig. 2. As we can see in the graph, the trajectories in the phase space converge to the equilibrium point (ϕ 0 , p ϕ,0 ) = (0, 0), indicating its stability. The physical interpretation of this key feature will be explored momentarily. Before that, we need to address the case where V 1 is positive.
The stability analysis of the equilibrium point (ϕ 0 , p ϕ,0 ) in the case V 1 > 0 also depends of the range of values assumed by ω. The two subcases are again ω < − 3 2 and ω > − 3 2 . The case ω < − 3 2 yields: The direction fields for the dynamics of ϕ when V 1 > 0 and ω < − 3 2 are displayed in Fig. 3. The trajectories converge to (ϕ 0 , p ϕ,0 ) = (0, 0); for this reason, it is classified as a stable equilibrium point. In this situation, V 1 > 0 and ω > − 3 2 , the eigenvalue λ + is a real positive number. As a consequence, the stability criterium Re [λ ± ] 0 is violated irrespective of the behaviour of the eigenvalue λ − . Hence, the equilibrium point (ϕ 0 , p ϕ,0 ) = (0, 0) is unstable. This conclusion is confirmed by the direction fields in Fig. 4 showing trajectories diverging from the equilibrium point.

Phase-space analysis with H = H (x 0 ) via Lyapunov's method
Presently, we generalize the autonomous system studied in Subsection 3.2 by allowing the Hubble function to be x 0 -dependent. In this case, the linear stability analysis may be insufficient, or even misleading, as per Ref. [83] (and Ref. [79]). We are thus compelled to use other methods of dynamical analysis. In particular, the Lyapunov's method-reviewed in Subsection 3.1.2-proves to be convenient in the present case.
The generalization of the system in Eq. (37) that takes into account H = H x 0 is: The equilibrium point is the same as in Eq. (38), namely: The (vector) function f (y) that determines the dynamics reads: where as it is immediately concluded from the comparison of the system (45) with the definition (33). Now we propose the following Lyapunov function V (y): with a and b constants. These constants will be determined in accordance with the Lyapunov Stability Theorem momentarily. The function in Eq. (49) immediately satisfies the property (i) in Subsection 3.1.2. Moreover, notice that This fact is relevant for checking if requirement (ii) is fulfilled ∀y ∈ U . We can compute the first identity in Eq. (34) for the specific case of our dynamical system: where we have made use of Eq. (45). The constant a and b are arbitrary; we use this freedom to choose: Accordingly, Eq. (51) yields: The choice in (52) lead us to satisfy property (iii) in Subsection 3.1.2. This is so because the Hubble function is always positive for an expanding Universe. Actually, we haveV < 0 since H x 0 > 0 and p 2 ϕ > 0, for p ϕ = 0. The equality in (iii) would only be achieved if b = 0, but this choice would lead to V = 0 in which case property (ii) in Subsection 3.1.2 would not be satisfied.
Incidentally, it is the time to attempt to verify property (ii). Plugging Eq. (52) into (49): Since b > 0 and ϕ 2 > 0 and p 2 ϕ > 0 in a neighbourhood of the equilibrium point, the decision on the fulfilment of requirement (ii) depends entirely on the values of parameters V 1 and ω. Essentially, the stability will be achieved if We identify four possibilities: (a) V 1 < 0 and ω < − 3 /2 then condition (55) is not satisfied, and the system is not stable; (b) V 1 < 0 and ω > − 3 /2 then condition (55) is satisfied, and the system is asymptotically stable; (c) V 1 > 0 and ω < − 3 /2 then the system is asymptotically stable; (d) V 1 > 0 and ω > − 3 /2 then the system is not stable.
These possibilities are consistent with those unveiled by the linear stability analysis performed in Subsection 3.2.
Here, however, the nature of the equilibrium point is kept even when the Hubble parameter is allowed to change as a function of x 0 . If we restrict our parameter space (ω, V 1 ) to satisfy the conditions for stability, i.e. conditions (b) or (c) above, then we verify that the equilibrium point y 0 = (ϕ 0 , p ϕ,0 ) = (0, 0) is globally asymptotically stable, since

Discussion: the meaning of stability
The results that have physical interest are those where the equilibrium point is stable, i.e. the case where V 1 < 0, ω > − 3 2 and the case with V 1 > 0, ω < − 3 2 . These different ranges of values for V 1 lead to two different interpretations. Let us reconsider Eq. (8) and include the potential (20) explicitly: As mentioned in the beginning of Section 3, we see that the parameter V 1 plays the role of a cosmological constant even in the dynamical phase of φ. It may be suggestive to rename this parameter V 1 ↔ 2Λ, in such a way that the left-hand side of Eq. (56) is formally the same as the geometrical side of Einstein field equations. Hence, the condition V 1 > 0 maps to Λ > 0 which indicates that the background solution to (56) is a de Sitter-like spacetime [74]. Conversely, the condition V 1 < 0 would correspond to Λ < 0 thus leading to AdS-like background. (Of course, some conditions apply for achieving a strict dS/AdS background: vacuum, derivatives of the field φ null, etc.).
In both cases, as long as we have stability of the equilibrium point, the system will evolve so that the trajectories of field φ in the phase space (ϕ, p ϕ ) converge to (ϕ 0 , p ϕ,0 ) = (0, 0). When this fixed point is reached, we have simultaneously: When the system gets to the equilibrium point, theṅ as implied from Eqs. (57) and (5). This means that the conditioṅ will necessarily be satisfied whenever the system is in equilibrium: Eq. (1) holds true precisely for the parameter value σ = 3. This value is exact here: the uncertainty of σ = 3 is theoretically zero, it is a natural consequence of the definition φ demanded by basic dimensional analysis, and out of the dynamics of φ in the phase space. Another consequence of this dynamics of φ is that the system can start with arbitrary initial conditions. The conditionĠ G = 3ċ c may not even be satisfied for this initial condition. It is just a matter of time until the system evolves and converge to Eq. (1).
We also point out the appearance of an effective cosmological constant when the dynamics of φ relax to the equilibrium point (regardless of the association V 1 ↔ 2Λ suggested above). In fact, the dynamics of the scalar field ceases at the equilibrium point, i.e. ∇ µ φ = 0 with φ → φ eq = constant; at this situation, the gravitational field equation (56) approaches: where T (eff) µν ≡ 1 c T µν is the conserved stress-energy tensor in our scalar-tensor model. The term in the parenthesis above appears as an effective cosmological constant, and the field φ, which converges to φ eq in Eq. (57), can be set as and Eq. (60) assumes the form of GR's field equation with cosmological constant. The quantities c 0 and G 0 are constant values of the couplings, which might as well be considered as their present-day values in cosmological terms. In fact, due to the last two equations, Eq. (60) assumes the following form at the present time: which is the familiar form of Einstein field equations, cf. Eq. (2). It is imperative that the following is cristal clear: the fact that φ eq = c 3 0 /G 0 = constant does not mean that G and c will be constants after the equilibrium is reached. All that is required thereafter is φ /φ = − Ġ /G − 3ċ/c = 0.
Accordingly, after the equilibrium is attained, it could be c = c 0 f and G = G 0 f 3 , where f = f x 0 is an arbitrary function of the time coordinate: these time-dependent ansatze for G and c satisfy both the requirements φ eq = c 3 0 /G 0 and Ġ /G − 3ċ/c = 0.
In the face of the comments above, one concludes that, if an astrophysical event takes place after the equilibrium condition is attained, it is essentially impossible to identify the dynamics of G separately from the dynamics of c. In order to identify the eventual dynamics of G and c withĠ G = 3ċ c (by taking into account only gravitational fields equations, as we have done here), one would have to consider situations out of the equilibrium for φ. If we are currently in the equilibrium condition, then one would have to take into consideration events in the past history of the universe. How far in that past one should go is still an open question, hinging on the amount of time that the system takes to converge to the equilibrium point. In any case, the condition in Eq. (21) must be preserved for the dynamical analysis in Section 3 to hold. As an example of what is stated above, we mention the solar system tests constraining the Brans-Dicke parameter ω > 40, 000 [96], while as a dimensionless parameter, it would be expected to be of order ω ≈ 1. In these tests, the Parameterized Post-Newtonian approximation is applied assuming that the dynamics of the scalar field plays a role in the solar system evolution. In our case, the dynamics ceased in the radiation era, way before the solar system was formed. From this perspective, it would be meaningless to try and discard our model based on BD constraints.
From the discussion above, the global picture of the dynamics in our scalar-tensor model should be at sight. The system begins with arbitrary initial conditions where φ = c 3 /G is out of equilibrium. In this phase, the functional form of G and c is not constrained to any particular behaviour interweaved. There are two regions in the parameter space of the pair (ω, V 1 ) enabling the system to reach an stable critical point, cf. Fig. 6. By restricting our model parameters to these regions, one guarantees that the stable equilibrium is reached, φ → φ eq , and the couplings G and c are forced to evolve respecting Ġ /G − 3ċ/c = 0 henceforth. The evolution of the universe from this point onward should be impacted by the fact that G scales as c 3 . In the next section we launch ourselves to the task of determining the background cosmological evolution of our scalar-tensor model after the equilibrium. Our final goal is to solve for the scale factor a as a function of the time-coordinate x 0 . Along the way, we build the continuity equation from the conservation of the effective stress-energy tensor mentioned below Eq. (60)-see Subsection 4.1. We will show that this continuity equation differs from the conventionally used in FLRW cosmology for it contains a term depending onċ. The ensuing VSL model is studied in Subsection 4.2 assuming radiation and dust-matter content. In both cases, the evolution of a = a x 0 tends to an accelerated de Sitter-type solution. As it happens, the latter is true for two classes of VSL models: those with a decreasing speed of light and those with an increasing c = c x 0 .

Analysis of the gravity field equation after φ reaches equilibrium
After φ reaches the equilibrium point, the field equation for the gravitational part of our scalar-tensor model, Eq. (60), reduces to, i.e.
where φ eq = constant. However, c = c x 0 on the r.h.s. accompanying T µν . In this section we investigate non-vacuum cosmological solutions stemming from this new feature. The 00-component and the ii-component of Eq. (64), under the FLRW metric in (13), read: andḢ

The covariance of the effective stress-energy tensor
By taking covariant divergence of Eq. (64) and by using Bianchi identity (∇ ρ G ρ ν = 0) and the metricity condition we conclude that: This is the extended conservation law in our scalar-tensor gravity and justifies the definition of T (eff) µν = 1 c T µν which is actually covariantly conserved. The effective stress-energy tensor T (eff) µν is the Noether's current related to the symmetry of (infinitesimal) general coordinate transformations. When we specify Eq. (67) for the FLRW metric, we get:ε The term on the r.h.s disappears forċ = 0 so that Eq. (68) recovers the standard continuity equation of background cosmology in GR in this particular case. Eqs. (65), (66) and (68) form the essential set of equation of cosmology in our scalar-tensor model allowing for the variation of c. Our proposal is to use the first and the third equations in the set above to determine the cosmic evolution. For that goal, two elements are necessary: 1. An equation of state p = p (ε);

A constitutive equation for
The EOS is taken as (Notice the distinction between w-the EOS parameter-and ω-the scalar-tensor constant parameter appearing in Eq. (7) as the coefficient of the kinetic term of φ.) The continuity equation (68) then becomes: The first term in the square brackets is the common contribution in standard GR FLRW cosmology. The second term in the square brackets is the contribution from our scalar-tensor model. The above equation can not be solved explicitly without a constitutive equation for c = c x 0 . That is what we will do next.

Modelling the varying speed of light
In a few models realizing VSL scenarios, it is assumed that c x 0 suffered a first order phase transition in the begining of the thermal history of the universe-see e.g., [52]. According to this image, the speed of light would have decreased its magnitude as the scale factor increased. This motivates the proposal in the next subsection. Before moving to the next subsection, it is probably worth emphasizing that the assumption of a particular ansatz for c = c (a), which would be valid after the dynamics of the field φ ceases, does not violate the equilibrium condition φ = φ eq = constant. For instance, if c ∼ 1/a it suffices that G ∼ 1/a 3 for maintaining φ = c 3 G = constant.

Speed of light scaling as the inverse of a
Let us assume that c = c 0 a 0 a .
That is: the speed of light decreases as the universe increases its size. The consequence of (71) is [97]: This is a reasonable choice that simplifies the continuity equation. In fact, plugging (72) into Eq. (70) enables immediate integration to: In particular, ε ∼ a −5 in a radiation dominated era (w = 1/3). This violates the Stefan-Boltzmann law ε ∼ T 4 because a ∼ T −1 , with T representing the temperature. (Of course, the violation of the Stefan-Boltzmann law by radiation in the model c ∼ a −1 occurs if we insist in a constant k B .) Curiously enough, dust-like matter (w = 0) does recover the Stefan-Boltzmann law ε ∼ a −4 in the our scalar-tensor model supplemented by the ansatz c ∼ a −1 . Substituting (71) and (73) into Eq. (65) with k = 0 and φ eq = c 3 0 /G 0 yields: under the definitions where [ε c,0 ] = (energy) / (length) 3 and Ω 0 is dimensionless.
In Eq. (76), x 0 0 is the value of the "time"-coordinate x 0 today. Eq. (76) is the exact solution for the scale factor. If Λ 0 a 2 0 Ω 0 (the cosmological constant dominates over the matter-energy content), the above solution reduces to: This is a de Sitter-type accelerated expansion and could describe dark energy. This is true regardless of the value of the parameter w, which means that both a radiation era and a matter dominated universe accommodate a de Sitter acceleration in the regime Λ 0 a 2 0 Ω 0 . Therefore, our scalar-tensor model with c ∼ a −1 is a promising candidate for fitting observational data. It remains to be investigated how the values of its parameters (Ω 0 , Λ 0 ) would be constrained by the data. Regardless, from the theoretical point of view, one could imagine a cosmic evolution from a radiation-dominated era to a matter-dominated universe to a dark-energy-dominated cosmos. This last step would come as a natural evolution of the scale factor in Eq. (76).
The question that poses itself at this point is the following. The picture for the universe's evolution in the previous paragraph was obtained from our scalar-tensor model and the particular ansatz c ∼ a −1 in which the speed of light decreases as the universe expands. How dependent this scenario is on the particular choice of a decreasing c? To put it another way: Could a different scenario with an increasing speed of light accommodate radiation-dominated and matter-dominated eras allowing for a subsequent dark-energy-type evolution?
In order to answer this question, we have to choose an ansatz for c = c (a) that is different from that in Eq. (71). We will do this in the next section by assuming an ansatz introduced by one of us (Gupta) in a series of papers [62,64,65,66,67,68,69,70,71].

Speed of light scaling as the exponential of a
Gupta's ansatz for an increasing speed of light has the form [66]: with α a positive constant of order one. Therefore, Inserting (79) into (70): For the standard non-varying speed of light picture, c = c 0 , the equation above recovers the traditional result from background FLRW cosmology of GR. This means that the Stefan-Boltzmann law ε ∼ a −4 is recovered in the radiation era (w = 1/3) except for the factor (c/c 0 ). This is different from the case in the previous subsection, where c ∼ a −1 and ε ∼ a −4 was realized only in the presence of dust-like matter. Now, Eq. (80) can be used in the first Friedmann equation (65) to yield the scale factor. Indeed, taking k = 0 and using the definitions of ε c,0 and Ω 0 in Eq. (75) results in: Surprisingly, this is exactly the same equation for H derived in the previous model-see Eq. (74). Therefore, the solution devised in Subsection 4.2.1 for the model with c ∼ a −1 is also valid in the present case, namely: Eq. (76) leading to the same de Sitter-type limiting case of Eq. (77). Now we are able to answer the question raised at the end of the last section. A dark-energy phase is not a privilege of decreasing speed of light scenarios. In our scalar-tensor model, an increasing speed of light also accommodates an accelerated expansion. Therefore there is enough freedom to imagine multiple eras of different kinematics for c = c x 0 , just like one admits different functions of a = a x 0 for distinct eras of radiation or matter domination. The universe might as well choose a decreasing speed of light in the beginning of its evolution later evolving to a configuration where c increases as x 0 increases.
By "beginning of its evolution" in the previous sentence we mean the period immediately after φ reaches its equilibrium point φ eq ; after all, the whole formalism developed in the present Section 4 assumes the validity of Eq. (64), which is the BD-like field equation for the gravitational field after φ approaches φ eq . It should be emphasized that the model in this subsection (Subsection 4.2.2) was indeed tested against a plethora of astrophysical and cosmological data with successes, including SNe Ia data [64], BAO and CMB data [66], BBN data [67], quasar data [68], gravitational lensing [69], orbital mechanics [70], the faint young-Sun problem [71].

Concluding remarks
In this paper, we studied a scalar-tensor model for gravity in which both the gravitational coupling G and the (causality) speed of light c are included in the scalar sector of the model through the field φ = c 3 /G.
The field equations for g µν and φ are built and specified for the homogeneous and isotropic cosmological background. The dynamics of our field φ is then analyzed in the phase space under some working hypotheses. We take the curvature parameter of the space sector k null. We use the first three terms in the series expansion of the potential V (φ). We assume that the dynamics of φ happens in a radiation era. We then study the evolution of φ via linear stability theory (reviewed in Subsection 3.1) assuming that the dynamical period of φ is short enough to assure that Hubble function is approximately constant therein (Subsection 3.2). Finally, the dynamical system analysis is extended to allow for time-dependent Hubble parameter; this requires the use of the Lyapunov's method (Subsection 3.3).
The phase space analysis shows that two sets of conditions upon the free parameters of the model lead to trajectories converging to a globally asymptotically stable equilibrium point. These conditions are: where V 1 is the constant coefficient of the linear term in the V (φ) series and ω is the constant coefficient of the kinetic term of φ in the model's action. The physical interpretation of V 1 and ω are commented on momentarily.
We have seen in Subsection 3.4 that V 1 could be understood as a type of cosmological constant during the dynamical phase of φ; therefore, its positivity (or negativity) could impact this version of Λ. Moreover, V 1 appears in the effective cosmological constant Λ 0 which turns up at the end of the dynamical evolution of φ, when the trajectories of the phase space converge toward the equilibrium point and stay there for the rest of the cosmic history. The definition of Λ 0 includes the other coefficients of the potential too: V 0 is the constant term in the V (φ) series, and V 2 is the coefficient of the quadratic term of V in φ. Our treatment demands V 0 = 0 but otherwise unconstrained while V 2 is totally unconstrained. The equilibrium value of the field, φ eq = c 3 0 /G 0 , also enters the definition of Λ 0 which leads to the interesting conclusion that the cosmological constant could depend on (G 0 , c 0 ), the (constant) equilibrium values of the gravitational coupling and of the speed of light.
The sign of the parameter ω is related to the physical nature of the scalar field φ. Indeed, if ω φ < 0 the sign of the kinetic term for φ in the action indicates that φ is a ghost field. This means that its Hamiltonian is unbounded from below and it would be source of negative energy. The global sign of the potential V (φ) is also important: together with the sign of the kinetic term, it determines if φ could also be an unphysical tachyon field breaking causality. In the vicinity of the equilibrium point, φ > 0, so that ω < 0 would characterize the field as a ghost field-see e.g., page 2 of [98] and references therein. This all means that we have to impose ω > 0 for a physically meaningful behaviour of φ.
The big picture in our development is: after evolving to the equilibrium stable point, the field φ reaches it equilibrium value φ eq = c 3 0 /G 0 and stays there. Then, the dynamics of φ ceases and two things happen: (i) the couplings G and c that are co-varying within φ assume the dependence G/G 0 = (c/c 0 ) 3 , where G 0 and c 0 could be interpreted as the numbers that we use as fundamental constants today; and (ii) the gravitational field equation of our scalar-tensor model, involving both g µν and φ, degenerates to a field equation that is formally the same as Einstein field equation of general relativity, but bears a time-dependent factor (1/c) alongside the ordinary energy momentum tensor. At the present time c = c 0 and the field equation recovers exactly Einstein equation of GR. That could explain why GR is so successful in describing local and low-redshift data but should be generalized in a larger scope of application. The relaxation of our model toward a field equation that is formally the same as that of GR is consistent with the works by Damour and Nordtvedt [99,100]-see also Ref. [101] by Faraoni and Franconnet. When the dynamics of φ ends,φ = 0, and the co-varying G and c are forced to obey the relation Ġ /G = σ (ċ/c) with σ = 3. This number is exact as far as the theoretical prediction is concerned (under the hypotheses we have assumed, cf. the beginning of Section 3). Our conclusion is that at least some astrophysical phenomena would be unchanged if G and c vary concurrently while respectingφ = 0. This fact has been pointed out in Refs. [70,69,65], in relation with orbital timing, strong gravitational lensing of SNe Ia, and core-collapse supernova events.
Our scalar-field model may be seen as a version of Brans-Dicke allowing for a varying speed of light (alongside the varying G). Although the field equations in Section 2 look alike those in BD theory, it must be underscored the fact that the varying c appears explicitly in the term containing the energy momentum tensor. This forces the definition of an effective stress-energy tensor, which is covariantly conserved thus assuring general covariance of the gravity equation. This was discussed also in Section 4, where we build the Friedmann equations for the cosmology within our scalar-tensor gravity. These equations are solved for two particular ansatze of c = c x 0 : (i) a speed of light that decreases as a increases, and (ii) a c that increases as the universe expands. It was shown that these both ansatze allow for an accelerated expansion consistent with the effect of dark energy. The fact that we were able to fully resolve the field equation for our scalar-tensor gravity with varying coupling constants shows that the generality of our model does not spoil practicality. The following logical step in our research in the future will be to refine the details of the picture by constraining the models with observational data.
Future perspectives include to relax some of the hypotheses adopted in this manuscript. One possibility is to let the dynamics of φ evolve in an era dominated by matter, for which ( − 3p) = 0 and the coupling of c = c x 0 with the matter Lagrangian comes into effect. The eventual interrelation between the matter part of the action and c = c x 0 brings forth the possibility of c being an independent scalar field. In this context, one would have to take G = G (x µ ) and c = c (x µ ) as scalar fields evolving on their own (not enclosed within the field φ). To put it another way, in the stability analysis performed in this paper we have one scalar degree of freedom, incarnated in φ; this is true even as both G and c are allowed to vary. In a future paper, we will consider two independent degrees of freedom, namely G and c, with their own separate kinetic terms and potentials (with or without a mutual interplay).
In this paper, the varying physical couplings were restricted to the set {G, c, Λ}-with the first two forming the field φ. Other works enlarge the set of possible varying fundamental constants to {G, c, Λ, , k B } and perform some phenomenological modelling in astrophysics and cosmology [65,102]. Pursuing the fundamentals surrounding an eventual variation of all these five couplings from a field theory perspective is something that may be done in the future.
Author Contributions: All authors have equally contributed to all parts of the paper. All authors have read and agreed to the published version of the manuscript.
Funding: Macronix Research Corporation: 2020-23; National Council for Scientific and Technological Development-CNPq: 309984/2020-3. This is Eq. (17) in the main part of the paper.