Dynamical Analysis of the Covarying Coupling Constants in Scalar–Tensor Gravity

Abstract

A scalar–tensor theory of gravity was considered, wherein the gravitational coupling G and the speed of light c were admitted as space–time functions and combined to form the definition of the scalar field $\varphi$. 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 in our gravity model. The effect of the cosmological coupling $\Lambda$ is accommodated within a possible behavior of $\varphi$. We analyzed the dynamics of $\varphi$ in the phase space, thereby showing the existence of an attractor point for reasonable hypotheses on the potential $V\left(\varphi \right)$ and no particular assumption on the Hubble function. The phase space analysis was performed both with the linear stability theory and via the more general Lyapunov method. Either method led 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 arrives at the globally asymptotically stable fixed point and the dynamics of $\varphi$ ceases. This result realized our main motivation: to provide 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 covarying couplings G and c impact the cosmic evolution after the dynamical system settles to equilibrium. The secondary goal of our work was to investigate how this impact occurs. This was performed by constructing the generalized continuity equation in our scalar–tensor model and considering two possible regimes for the varying speed of light—decreasing c and increasing c—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.

1. Introduction

The possibility of the 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 a deep connection with scalar–tensor theories, themselves a popular candidate for viable modified theories of gravity. This paper explored the possible interwoven variation of the couplings G and c (while also accommodating the effects of $\Lambda$) 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 $\alpha$. While G is a dimensionful constant, $\alpha$ is a dimensionless constant. Uzan [1,2] 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 normalized 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 19th 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 the 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 from binary neutron stars [19]. Almost all of them have resulted in constraints on $\dot{G}/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 that could vary spatially and temporally. The 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 a geometrical nature [22]. A more general example within this class is the scalar–tensor theory in the Lyra manifold [23,24]. One of the reasons why scalar–tensor theories of gravity raise so much interest (see [22,25] and the references therein) is their well-known equivalence with some modified gravity theories [26], cf., e.g., [27,28]. Among these possible modifications of gravity, we mention the theories $f\left(R\right)$ [29,30,31], $f\left(R,{\nabla }^{n}R\right)$ [32], $f\left(R,{\square }^{n}R\right)$ [33], and $f\left(T\right)$ [34,35], which have been proposed as attempts to deal with the GR quantization problem [36,37], to realize inflationary models [38], and to address the dark energy problem [39]. References [40,41,42] studied the phenomenology of extended Jordan–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 [43,44,45,46].

The constancy of the speed of light is the foundation of the 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., Reference [57] considered varying G and $\Lambda$; the variation of $G/c^2$ was studied in [58]; G and $\alpha$ changed concomitantly in [59]; covarying c, G, and $\Lambda$ were the subject of [60,61,62]; the set $\{c,G,\alpha \}$ was allowed to vary in [63]; finally, References [64,65] suggested that all the couplings G, c, $\Lambda$, ħ, and $k_B$ must covary.

After this contextualization of the subject of varying physical constants, let us concentrate on the topic to be explored here. Our focus in this paper was on the variation of c and G and the interrelationship of their variation. We will see that this interrelationship yields naturally an effective cosmological constant. Gupta [66] introduced the general ansatz:

$$\frac{\dot{G}}{G}=\sigma \frac{\dot{c}}{c}$$

(1)

inspired by the general constraint appearing in the proposal by Costa et al. [61]. The quantity $\sigma$ is a constant parameter. The value $\sigma =3$ is strongly favored by several phenomenological applications of Equation (1) in cosmology [64,66,67,68] and astrophysics [65,69,70,71]. This preferred value for $\sigma$ is corroborated by other authors, cf., e.g., [63]. For this reason, we set off to investigate the fundamental reasons that might be underlying this fact.

We did this by assuming a modification in the standard theory of gravity: General Relativity (GR). This is justified by the fact that GR does not allow for varying physical constants. In fact, the couplings $G,$ c, and $\Lambda$ showing up in the Einstein field equations:

$$G_{\mu \nu }+\Lambda g_{\mu \nu }=\frac{8\pi G}{c^4}{T}_{\mu \nu }$$

(2)

do not vary in time or space. Recall the definition of the Einstein tensor: $G_{\mu \nu }=R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R$, where $R_{\mu \nu }$ is the Ricci tensor and R is the curvature scalar. We followed the metric signature and Riemann tensor definition of [72], i.e., the Minkowski metric in Cartesian coordinates reads $\left({\eta }_{\mu \nu }\right)=\mathrm{diag}\left(-1,+1,+1,+1\right)$, $R_{\mu \nu }=R_{\phantom{\rho }\mu \rho \nu }^{\rho }$ with the curvature tensor calculated by $R_{\phantom{\rho }\sigma \mu \nu }^{\rho }={\partial }_{\mu }{\mathsf{\Gamma }}_{\phantom{\rho }\nu \sigma }^{\rho }+{\mathsf{\Gamma }}_{\phantom{\rho }\mu \lambda }^{\rho }{\mathsf{\Gamma }}_{\phantom{\rho }\nu \sigma }^{\lambda }-\left(\mu \leftrightarrow \nu \right)$, $R=g^{\mu \nu }{R}_{\mu \nu }$, and ${\mathsf{\Gamma }}_{\phantom{\rho }\mu \nu }^{\rho }=\frac{1}{2}{g}^{\rho \sigma }\left({\partial }_{\nu }{g}_{\sigma \mu }+{\partial }_{\mu }{g}_{\nu \sigma }-{\partial }_{\sigma }{g}_{\mu \nu }\right)$ stands for the Christoffel symbols built from the metric tensor $g_{\mu \nu }$ and its first-order derivatives. Covariant derivatives were built from the Christoffel, e.g., ${\nabla }_{\mu }{V}^{\nu }={\partial }_{\mu }{V}^{\nu }+{\mathsf{\Gamma }}_{\phantom{\nu }\mu \rho }^{\nu }{V}^{\rho }$, for an arbitrary contravariant vector $V^{\nu }$.

Our generalization of GR stems from an action integral that is inspired by the Einstein–Hilbert action:

$$S_g=\frac{1}{16\pi }\int d^{4}x\sqrt{-g}\left[\frac{c^3}{G}\left(R-2\Lambda \right)\right],$$

(3)

but with the couplings $\left\{G,c,\Lambda \right\}$ being spacetime functions. At this point, we emphasize the natural appearance of the couplings $\left\{G,c,\Lambda \right\}$ 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 $\left(R-2\Lambda \right)$ 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 $\left\{G,c,\Lambda \right\}$ justifies why we considered 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., [64,65,66,67,68,71,73].

As usual, the four-volume in the action $S_g$ is $d^{4}x=d{x}^0{d}^3\mathbf{x}$. It is worth stressing 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 such as having to deal with the opened $x^0$-differential $d{x}^0=d\left(ct\right)=cdt+tdc$ are avoided; the speed of light will not appear in the metric components $g_{\mu \nu }$ explicitly. Consequently, the determinant of the metric tensor $g=g\left(x^{\mu }\right)$ 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\to t^{\prime }$ compensates for this fact in each reference frame. This strategy is somewhat different from the procedure of the c-flation framework developed in [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 Equation (1) indicates that the couplings G and c are intertwined. In fact, it could be recast as

$$\frac{\dot{\varphi }}{\varphi }=0,$$

(4)

where

$$\varphi \equiv \frac{c^3}{G}$$

(5)

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 $\dot{\varphi }=\frac{d\varphi }{d{x}^0}$.

The simplest way to satisfy (1) is to admit that $\dot{\varphi }=-\varphi \left[\frac{\dot{G}}{G}-3\frac{\dot{c}}{c}\right]=0$ due to constant couplings, in which case $\dot{G}=0$ and $\dot{c}=0$. This is the way of GR, but it is not the most-general possibility and not the one we adopted here. An earlier attempt to relax the constancy of fundamental couplings in the context of a theory of 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]:

$${\varphi }_{\mathrm{BD}}=\frac{1}{G}.$$

(6)

Our proposal is to extend the scope of the 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 $\varphi$, as previously envisioned by Brans and Dicke (in fact, we geometrized not only the gravitational coupling—through G—but also the causality coupling c). This is true even from a dimensional point of view. In the Brans–Dicke picture, ${\varphi }_{\mathrm{BD}}$ is measured in units of $\left(\mathrm{mass}\right)\times {\left(\mathrm{time}\right)}^2/{\left(\mathrm{length}\right)}^3$. In our scalar–tensor gravity, $\varphi$ has dimensions of $\left(\mathrm{mass}\right)/\left(\mathrm{time}\right)$.

The fact that a varying c enters the scalar field $\varphi$ in Equation (5) may be a cause for 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 [3]. There could be different types of speed of light, viz. the spacetime speed of light $c_{\mathrm{ST}}$ related to Lorentz invariance and causality and the electromagnetic speed of light $c_{\mathrm{EM}}$ built from the couplings in electromagnetism: the (vacuum) electric permittivity ${ϵ}_0$ and the magnetic permeability ${\mu }_0$. Here, we interpreted c as the causality coupling $c_{\mathrm{ST}}$ and ensured 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_{\mathrm{ST}}=c_{\mathrm{EM}}=c_0$ currently.

The main goal of this paper was to investigate if the dynamics of $\varphi$ can naturally lead to $\dot{G}/G=3\left(\dot{c}/c\right)$ as an attractor solution. This would give the fundamental reason for why the value $\sigma =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\left(x^0\right)$ and $c\left(x^0\right)$ is established, a natural question poses itself: “What are the consequences of $G=G\left(c\right)$ for the cosmic evolution?”. Our secondary goal in this paper was 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 the recent accelerated regime). We built our scalar–tensor model in a metric-compatible, four-dimensional, torsion-free, Riemannian manifold. Approaches considering a 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 and derive the equations of motion for the tensor field $g_{\mu \nu }$ and the scalar field $\varphi$, which together build the gravitational field. The field equations are specified for the FLRW metric [74], since our main interest was to study the dynamics of $\varphi$ in the cosmological context. This study is carried out in Section 3, where it is shown that an attractor solution leading to $\sigma =3$ is attainable under reasonable hypotheses for the potential $V\left(\varphi \right)$ and the dominant matter–energy content during the evolution of $\varphi$ (both the linear stability theory and Lyapunov’s method are briefly reviewed in Section 3.1 for the sake of completeness of our presentation and convenience for the reader). This finding comes from the linear stability analysis of our dynamical system—performed in Section 3.2; the result is confirmed and strengthened in Section 3.3, where we used the general Lyapunov’s method to show that the critical point ${\varphi }_{\mathrm{eq}}$ leading to $\sigma =3$ is a globally asymptotically stable fixed point. Section 3.4 gives some details about the system’s attractor point and emphasizes the consequences of our modified gravity. The cosmic evolution after $\varphi$ 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\left(x^0\right)$, 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 $\mathbf{\varphi }={\mathit{c}}^{\mathbf{3}}/\mathit{G}$

The motivations in the previous section led us to describe gravity from a precisely specified form for the action integral. In fact, this form is hinted at by the coefficient of the curvature scalar R in Equation (7), namely $\left(c^3/G\right)$, which was defined as the scalar field $\varphi$. Accordingly, we introduce

$$\begin{array}{cc}\hfill S& =\frac{1}{16\pi }\int d^{4}x\sqrt{-g}\left(\varphi R-\frac{\omega }{\varphi }{\nabla }_{\mu }\varphi {\nabla }^{\mu }\varphi -V\left(\varphi \right)\right)\hfill \\ \hfill & \phantom{=}+\int d^{4}x\sqrt{-g}\left(\frac{1}{c}{\mathcal{L}}_m\right),\hfill \end{array}$$

(7)

where $\omega$ is a dimensionless constant of order one (it is argued that the Brans–Dicke theory recovers general relativity in the limit $\omega \to \infty$, cf., e.g., [74]. However, this is not always the case: Reference [75] presented a counter-example) and $\varphi$ is a scalar field. The constant $\omega$ appears as the coefficient of the kinetic term for the field $\varphi$; the denominator of this term includes $\varphi$ for dimensionality consistency. The cosmological coupling $\Lambda$ was not written explicitly in (7) because it can be included in the potential $V\left(\varphi \right)$. In fact, we will show that a constant $\Lambda$ 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 $\left(\mathrm{energy}\right)\times \left(\mathrm{time}\right)$. This includes inserting the factor $\left(1/c\right)$ explicitly inside the matter term integral containing the matter Lagrangian density ${\mathcal{L}}_m$, itself carrying dimensions of energy density.

The action (7) is formally the same as the Brans–Dicke action. While the Brans–Dicke scalar field is usually interpreted as an effective $G^{-1}$ (see [22]), our $\varphi$ 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 [50,52,53]—for reasons ranging from different modeling to the basic fact that variations of the gravitational coupling were disregarded therein. Additionally, it should be pointed out that (7) does not correspond to the same treatment given in [61]; indeed, the covariant c-flation framework considers G and c as independent fields, while here, these varying couplings entered 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).

The variation of the action in Equation (7) with respect to the metric $g^{\mu \nu }$ yields the BD-like field equations for $g_{\mu \nu }$ (it is a matter of interpretation to label the model explored in this paper as a Brans–Dicke theory. Our model is not Brans–Dicke if one considers as Brans–Dicke a scalar–tensor theory where the scalar field $\varphi$ is related strictly to $1/G$. Our model could be said to be Brans–Dicke if the latter is a scalar–tensor theory for $\varphi$ regardless of its interpretations in terms of the fundamental couplings. We prefer the first interpretation. However, even in the second interpretation, our paper showed that new physics arises from defining $\varphi$ in terms of both G and c: this novel perspective allows for the entangled variation of G and c in a theoretical structure accommodating a consistent variational principle, the covariant conservation of an effective energy momentum tensor, and a fundamental explanation for recent enticing phenomenological results pointing towards $\dot{G}/G=3\left(\dot{c}/c\right)$.):

$$\begin{array}{cc}\hfill R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R+\frac{V}{2\varphi }{g}_{\mu \nu }& =\frac{1}{c}\frac{8\pi }{\varphi }{T}_{\mu \nu }\hfill \\ \hfill & +\frac{\omega }{{\varphi }^2}\left({\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\rho }\varphi {\nabla }_{\rho }\varphi \right)\hfill \\ \hfill & +\frac{1}{\varphi }\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi -g_{\mu \nu }\square \varphi \right),\hfill \end{array}$$

(8)

where

$$T_{\mu \nu }=\frac{-2}{\sqrt{-g}}\frac{\delta \left(\sqrt{-g}{\mathcal{L}}_m\right)}{\delta g^{\mu \nu }}$$

(9)

is the stress–energy tensor of ordinary matter. Notice that $T_{\mu \nu }$ in (9) has the regular dimensions of the energy density. However, it should be underscored that the variation of the matter action:

$$S_m=\int d^{4}x\sqrt{-g}\left(\frac{1}{c}{\mathcal{L}}_m\right)$$

(10)

with respect to $g^{\mu \nu }$ is written as

$$\delta S_m=-\int d^{4}x\sqrt{-g}\left(\frac{1}{c}\right)\left(\frac{1}{2}{T}_{\mu \nu }\right)\delta g^{\mu \nu }$$

(11)

in the context of our scalar–tensor model. Here, the factor $\left(1/c\right)$ must participate in the form of $\delta 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 $\left(1/c\right)$ accompanying ${\mathcal{L}}_m$ (or $T_{\mu \nu }$) brings about one of the interesting features 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 $\varphi$ through its kinetic and potential terms. In Equation (9), no dependence of ${\mathcal{L}}_m$ on the field $\varphi$ is assumed. However, it should be noted that we cannot evade a non-trivial coupling with the (varying) speed of light within the functional derivative entering the definition of $T_{\mu \nu }$. This feature will be dealt with later—see Section 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 its own stress–energy 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 appeared in a different form: the covariant c-flation approach introduces $T_{\mu \nu }$ directly into the action via a Lagrange multiplier.

The variation of S with respect to the scalar field $\varphi$ leads to the Equation Of Motion (EOM) for the scalar part of our scalar–tensor model:

$$\frac{2\omega }{\varphi }\square \varphi +R-\frac{\omega }{{\varphi }^2}{\nabla }^{\rho }\varphi {\nabla }_{\rho }\varphi -\frac{dV}{d\varphi }-16\pi \frac{1}{c^2}\frac{dc}{d\varphi }{\mathcal{L}}_m=0.$$

(12)

This EOM says that “$\varphi$ is a dynamical field; it changes in space and time.” As a consequence of the definition of $\varphi$ in terms of G and c, the latter conclusion yields the covarying character of the gravitational coupling and the causality coupling. We show below that the simultaneous variation of G and c occurs in such a way that the condition $\dot{G}/G=3\left(\dot{c}/c\right)$ 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 $\varphi$ depends only on the cosmic “time” $x^0=ct$ due to the requirements of homogeneity and isotropy. These requirements demand the FLRW line element:

$$d{s}^2=-{\left(d{x}^0\right)}^2+a^2\left(x^0\right)\left[\frac{d{r}^2}{1-k{r}^2}+r^2\left(d{\theta }^2+{sin}^2\theta d{\phi }^2\right)\right].$$

(13)

In this set of coordinates, the scale factor $a\left(x^0\right)$ bears the dimension of length. Conversely, the comoving coordinates $\left\{r,\theta ,\phi \right\}$ 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

$$H=\frac{\dot{a}}{a}$$

(14)

and has the dimension of inverse length.

The matter content appearing in both field Equations (8) and (12) shall be modeled by the perfect fluid stress–energy tensor:

$$T_{\phantom{\mu }\nu }^{\mu }=\mathrm{diag}\left\{-\epsilon ,p,p,p\right\},$$

(15)

where p is the pressure and $\epsilon$ is the energy density (both with dimensions of energy per unit volume).

Due to (15), the 00-component of the $g_{\mu \nu }$ field Equation (8) reads (additional details on how to calculate Equations (16)–(18) are given in Appendix A):

$$H^2=\frac{1}{c}\frac{8\pi }{3\varphi }\epsilon +\frac{\omega }{6}{\left(\frac{\dot{\varphi }}{\varphi }\right)}^2-H\frac{\dot{\varphi }}{\varphi }-\frac{k}{a^2}+\frac{V}{6\varphi }.$$

(16)

This is the (first) Friedmann equation of our scalar–tensor cosmology with covarying G and c through the field $\varphi$.

The second Friedmann equation (or acceleration equation) is calculated by taking the trace of Equation (8), using the curvature scalar built from the metric in (13), utilizing the trace of the stress energy tensor, and Equations (12) and (16). It reads:

$$\begin{array}{cc}\hfill \dot{H}& =-\frac{1}{c}\frac{8\pi }{\left(2\omega +3\right)\varphi }\left[\left(\omega +2\right)\epsilon +\omega p\right]-\frac{\omega }{2}{\left(\frac{\dot{\varphi }}{\varphi }\right)}^2\hfill \\ \hfill & \phantom{=}+2H\frac{\dot{\varphi }}{\varphi }+\frac{k}{a^2}+\frac{1}{2\left(2\omega +3\right)\varphi }\left(\varphi \frac{dV}{d\varphi }-2V+16\pi \frac{\varphi }{c}\frac{dc}{d\varphi }\frac{{\mathcal{L}}_m}{c}\right).\hfill \end{array}$$

(17)

Moreover, in the cosmological context, Equation (12) reduces to

$$\ddot{\varphi }+3H\dot{\varphi }=\frac{1}{2\omega +3}\left[8\pi \frac{1}{c}\left(\epsilon -3p\right)-\varphi \frac{dV}{d\varphi }+2V-16\pi \frac{\varphi }{c}\frac{dc}{d\varphi }\frac{{\mathcal{L}}_m}{c}\right].$$

(18)

This equation determines the dynamics of $\varphi$. It depends on $H=H\left(x^0\right)$, which means that (16)–(18) should be solved simultaneously. The latter also depends on the matter content through the factor $\left(\epsilon -3p\right)$, which demands us to specify the nature of the perfect fluid at play in the cosmological era under scrutiny. Finally, Equation (18) exhibits a dependence on the potential $V\left(\varphi \right)$ and its derivative with respect to the field $\varphi$. Ultimately, we will have to say what is the form that we expect for potential $\varphi$, 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 regard to dynamical systems approach in cosmology, we refer the reader to the interesting papers [79,80,81,82,83,84,85,86,87].

3. The Dynamics of $\mathbf{\varphi }$

The potential $V\left(\varphi \right)$ is assumed to be analytical, but otherwise completely general. It was taken as dominant over the contribution of the term involving the matter Lagrangian; more specifically,

$$\left|16\pi \frac{\varphi }{c}\frac{dc}{d\varphi }\frac{{\mathcal{L}}_m}{c}\right|\ll \left|-\varphi \frac{dV}{d\varphi }+2V\right|.$$

(19)

This condition may be satisfied either for $\left|\frac{\varphi }{c}\frac{dc}{d\varphi }\right|\ll 1$ or $\left|\frac{{\mathcal{L}}_m}{c}\right|\ll \left|-\varphi \frac{dV}{d\varphi }+2V\right|$. We shall consider the first terms of the expansion of $V\left(\varphi \right)$ in a power series:

$$V\left(\varphi \right)=V_0+V_1\varphi +V_2{\varphi }^2,$$

(20)

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 for the potential description around a local minimum of the function. It is 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., [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 [91]. A power series expansion up to second-order can be found in [89]; the interpretation of the different terms that participate in a truncated series such as that in Equation (20) were discussed in [22,92], for instance. The form (20) encodes a few cases of special interest:

• Let $V_1=V_2=0$. This realizes a constant potential that could be rescaled by setting $V_0=0$.
• For $V_0=V_2=0$ and $V_1=2\Lambda$, one recovers the ordinary cosmological constant ($\Lambda =\mathrm{constant}$) already present in GR. The enticing feature here is that $V\left(\varphi \right)=2\Lambda \varphi$ would exhibit all the fundamental couplings in gravity $\left\{G,c,\Lambda \right\}$ when we substitute $\varphi =c^3/G$. (Nevertheless, recall that $\Lambda$ would still be a genuine constant here since it defines $V_1$.)
• If $V_0=V_1=0$ and $V_2=\frac{1}{2}{m}_{\varphi }^2$, the potential accounts for a massive term in $\varphi$’s field equation. In fact, in this case, $V\left(\varphi \right)=\frac{1}{2}{m}_{\varphi }^2\varphi$ and $\frac{d^{2}V}{d{\varphi }^2}=m_{\varphi }^2$, thus corroborating the interpretation. The delicate part is to interpret the role of this mass in $\varphi =c^3/G$. It could be conjectured that the source of $\varphi$’s mass is either a massive photon or a massive graviton.

In view of the above features, we considered the $V\left(\varphi \right)$ in Equation (20) adequate enough to allow for a sufficiently general dynamical analysis of our scalar–tensor theory of covarying physical couplings in gravity.

Substitution of Equation (20) into (18) yields a form for the field equation of $\varphi$ that does not depend on $V_2$. This is not the final possible simplification to this equation. The next step is to consider periods in the cosmic history when

$$\left|\frac{1}{c}\left(\epsilon -3p\right)\right|\ll \left|\varphi \frac{dV}{d\varphi }-2V\right|,$$

(21)

or, more specifically,

$$\left(\epsilon -3p\right)=0.$$

(22)

The above constraints not only guarantee a great simplification of Equation (18), but they are also very reasonable. Indeed, Equation (22) is consistent with a radiation-filled Universe—described by the equation of state $p=\epsilon /3$. By abiding by Equation (22), we assumed that the dynamics of the field $\varphi$ 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 $\varphi$ would have taken place.

In any case, the essential ingredient to the remainder of our analysis is that the condition in Equation (21) is satisfied during all the dynamics of $\varphi$. This means that the combination $\left(\varphi \frac{dV}{d\varphi }-2V\right)$ built from the potential $V\left(\varphi \right)$ should dominate the matter–energy content entering into play via $\frac{1}{c}{T}_{\mu \nu }$. This configuration could be easily realized within varying c frameworks during regimes where $\left(c/c_0\right)\gg 1$, for example in the very early Universe, before the phase transition of the decaying $c\left(x^0\right)$ admitted by some VSL models [1].

Plugging (22) into (18) results in

$$\ddot{\varphi }+3H\dot{\varphi }-\left(\frac{V_1}{2\omega +3}\right)\varphi -\left(\frac{2V_0}{2\omega +3}\right)=0.$$

(23)

A solution to this equation leads to $\varphi =\varphi \left(x^0\right)$, which constrains the dynamics of $\left(c^3/G\right)$. The physical solution to (23) must avoid the point $\varphi =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 $\varphi$ by defining

$$\phi \equiv \varphi +2\frac{V_0}{V_1},$$

(24)

where we admit $V_0\ne 0$ and $V_1$ finite. In terms of this new field variable, Equation (23) reads

$$\ddot{\phi }+3H\dot{\phi }-\frac{V_1}{2\omega +3}\phi =0.$$

(25)

The solution to this equation is not so easy because $H=H\left(x^0\right)$, and this is given by Equation (16). (Even after we solve Equation (25) for $\varphi =\varphi \left(x^0\right)$ using $H=H\left(x^0\right)$, we should have to justify a particular ansatz for $c=c\left(x^0\right)$, which would lead to $G\left(x^0\right)$.)

3.1. Elements of Stability Theory for Dynamical Systems

We shall study Equation (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 given in the following subsections, whose discussions are based on the papers [79,82] and the references therein. Check also [94].

3.1.1. Linear Stability Theory

Let us suppose a physical system can be described by n variables $\mathbf{y}=y^1,y^2,\dots ,y^n$. Each of these variables is a real-valued function, so that $\mathbf{y}\in M\subset {\mathbb{R}}^n$. Now, consider that the dynamics of this physical system can be described by a set of n differential equations of first order:

$$\dot{\mathbf{y}}=\mathbf{f}\left(\mathbf{y}\right),$$

(26)

where the dot indicates the derivative with respect to a parameter representing the evolution of this system, say a time coordinate ($x^0$). Here, $\mathbf{f}\left(\mathbf{y}\right)$ represents a set of n analytical functions (which can be nonlinear) of the variables $\mathbf{y}$. Whenever there is no explicit dependence of $\mathbf{f}$ with respect to $x^0$, the set of equations is called an autonomous system.

If there exists a point ${\mathbf{y}}_0\in M$ such that $\mathbf{f}\left({\mathbf{y}}_0\right)=0$, then

$${\left.\dot{\mathbf{y}}\right|}_{\mathbf{y}={\mathbf{y}}_0}=0.$$

(27)

If there is a time $x_{\left(f\right)}^0$ where $\mathbf{y}\left(x_{\left(f\right)}^0\right)={\mathbf{y}}_0$, then the dynamics of the system ceases and $\mathbf{y}\left(x^0\right)={\mathbf{y}}_0$ for any value of $x^0>x_{\left(f\right)}^0$. 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 neighborhood $U\subset M$ of ${\mathbf{y}}_0$. Roughly speaking, if the trajectories in U diverge from the fixed point, it is considered unstable; if the trajectories converge to ${\mathbf{y}}_0$, then it is stable.

The simplest stability analysis that can be carried out for such a system consists of considering the equations in a small region U where $\mathbf{f}\left(\mathbf{y}\right)$ can be expanded in a Taylor series about ${\mathbf{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

$$\begin{array}{cc}\hfill \dot{\mathbf{y}}& \approx \mathbf{f}\left({\mathbf{y}}_0\right)+\left(y^i-y_0^i\right)·\frac{\partial \mathbf{f}}{\partial y^i}\left({\mathbf{y}}_0\right)\hfill \\ \hfill & =\left[\left(\mathbf{y}-{\mathbf{y}}_0\right)·\mathrm{grad}\right]\mathbf{f}\left({\mathbf{y}}_0\right).\hfill \end{array}$$

(28)

A change of variable, $\mathbf{z}=\mathbf{y}-{\mathbf{y}}_0$, can always be proposed so that the equilibrium point is the origin of the coordinate system in the new variable:

$$\dot{\mathbf{z}}\approx \left(\mathbf{z}·\mathrm{grad}\right)\mathbf{f}\left(0\right),$$

(29)

or in a matrix form:

$$\dot{\mathbf{z}}\approx \mathbf{M}·\mathbf{z},$$

(30)

$$\left(\mathbf{M}\right)={\left(\begin{array}{cccc}\frac{\partial f^1}{\partial z^1}& \frac{\partial f^1}{\partial z^2}& \cdots & \frac{\partial f^1}{\partial z^n}\\ \frac{\partial f^2}{\partial z^1}& \frac{\partial f^2}{\partial z^2}& \cdots & \frac{\partial f^2}{\partial z^n}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{\partial f^n}{\partial z^1}& \frac{\partial f^n}{\partial z^2}& \cdots & \frac{\partial f^n}{\partial z^n}\end{array}\right)}_{{\mathbf{y}}_0}.$$

(31)

$\mathbf{M}$ is called the Jacobian matrix or stability matrix. As long as $det\mathbf{M}\ne 0$, we can look for a new transformation of variables $\mathbf{z}\to \mathbf{Z}$ where the transformed Jacobian matrix becomes diagonal: the elements of the diagonal are actually the eigenvalues, ${\lambda }_i$, of the Jacobian matrix. In this case, the set of differential equations will have the following solution:

$$Z^i=Z_0^i{e}^{{\lambda }_i{x}^0}.$$

(32)

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 of the study of their real parts. Notice that the Jacobian matrix is composed of constant coefficients once $\mathbf{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 $\mathbf{f}$ is valid.

The above considerations shall be employed in Section 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 a dynamical system when the linear stability analysis fails is Lyapunov’s method, presented in the next section.

3.1.2. Lyapunov’s Method

Section 3.1.2 follows closely the treatment given in [79]. This method is more powerful and general than the one in Section 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 as well, something the linear stability theory fails to accomplish. A critical point (fixed point) $\mathbf{y}={\mathbf{y}}_0\in M\subset {\mathbb{R}}^n$ of the dynamical system (26):

$$\dot{\mathbf{y}}=\mathbf{f}\left(\mathbf{y}\right),$$

(33)

is classified as a hyperbolic point if all of the eigenvalues of the stability matrix have a non-zero real part. Otherwise, it is said to be non-hyperbolic.

Lyapunov’s method for checking the stability of a dynamical system consists of finding a Lyapunov function $V\left(\mathbf{y}\right)$, such that:

(i)

$V\left(\mathbf{y}\right)$ is differentiable in a neighborhood U of ${\mathbf{y}}_0$;

(ii)

$V\left(\mathbf{y}\right)>V\left({\mathbf{y}}_0\right)$;

(iii)

$\dot{V}\le 0$, $\forall \mathbf{y}\in U$.

There is no known systematic way of deriving the Lyapunov function $V\left(\mathbf{y}\right)$. If a Lyapunov function exists, the Lyapunov stability theorem establishes that the critical point ${\mathbf{y}}_0$ is a stable fixed point if the requirement $\dot{V}\le 0$ is fulfilled. The critical point ${\mathbf{y}}_0$ is an asymptotically stable fixed point if there is a Lyapunov function $V\left(\mathbf{y}\right)$ satisfying the criterion $\dot{V}<0$. In addition, if ${lim}_{∥\mathbf{y}∥\to \infty }V\left(\mathbf{y}\right)\to \infty$, $\forall \mathbf{y}\in U$, then the critical point ${\mathbf{y}}_0$ is classified as globally stable or globally asymptotically stable, respectively. See [79].

Notice that Requirement (iii) can be cast into the form:

$$\dot{V}\left(y_1,...,y_n\right)=\frac{\partial V}{\partial y_1}{\dot{y}}_1+\cdots +\frac{\partial V}{\partial y_n}{\dot{y}}_n=\mathrm{grad}V·\mathbf{f}\left(\mathbf{y}\right)\le 0.$$

(34)

Equation (33) was used in the last step. The relation in Equation (34) will be important in Section 3.3 below.

3.2. Phase Space Analysis with $H\left(x^0\right)\approx H_{*}$ via the Linear Stability Theory

We avoided the difficulty mentioned below in Equation (25) by assuming, in this subsection, that the Hubble function is constant during the time interval in which the dynamics of the field $\varphi$ occurs. In other words,

$$H\left(x^0\right)\approx H_{*}=\mathrm{constant}.$$

(35)

It should be emphasized that this condition is not as restrictive as it may seem. It does not mean that $H\left(t\right)$ is always constant, but only during the time interval of $\phi$’s evolution. We are not demanding the Universe to be stationary (with $\dot{a}=\mathrm{constant}$); we are essentially assuming that H evolves slowly during the time it takes $\phi$ to reach an eventual attractor point. Again, this is a first-order phase transition scenario common to several models of the VSL [1].

In order to show that the picture described in the previous paragraph is achievable, let us analyze the dynamics of $\phi$ in the phase space resorting to the linear stability theory [79]. Let

$$p_{\phi }\equiv \dot{\phi }$$

(36)

be the momentum associated with $\phi$. Then, Equations (25) and (36) become the following coupled system of equations:

$$\left\{\begin{array}{c}\dot{\phi }=p_{\phi }\hfill \\ {\dot{p}}_{\phi }=-3H_{*}{p}_{\phi }+\frac{V_1}{2\omega +3}\phi \hfill \end{array}\right..$$

(37)

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 for the case in hand. The result stemming from this analysis should not suffer from any pathologies, such as those mentioned, e.g., in [83]—see also [95]. A quick inspection of (37) revels the equilibrium point:

$$p_{\phi ,0}={\phi }_0=0.$$

(38)

The next step is to analyze the stability of the fixed point according to the Lyapunov criterion [94]. The Lyapunov coefficients ${\lambda }_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 with this matrix and read:

$${\lambda }_{\pm }=\frac{3H_{*}}{2}\left(-1\pm \sqrt{1+\frac{V_1}{3}\frac{1}{\left(1+\frac{2}{3}\omega \right)}{\left(\frac{2}{3H_{*}}\right)}^2}\right).$$

(39)

The autonomous system can be considered stable if the real parts of both eigenvalues satisfy [79]

$$\mathrm{Re}\left[{\lambda }_{\pm }\right]⩽0.$$

(40)

The analysis of the direction fields in the phase space allowed 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 criterion. Here, we assumed that

$$H_{*}>0,$$

(41)

i.e., our Universe is expanding.

We split our analysis to cover different possible ranges of values for the parameters $V_1$ and $\omega$.

Case 1: Negative-valued $V_1$

Here, we took $V_1<0$. There are two subcases to be analyzed according to the values of the parameter $\omega$. These subcases are considered separately next.

Subcase 1.1: Negative-valued $\omega$ with $\omega <-\frac{3}{2}$

In the instance $V_1<0$ and $\omega <-\frac{3}{2}$, we may write

$$\frac{\left|V_1\right|}{3}\frac{1}{\left|\frac{2}{3}\left|\omega \right|-1\right|}>0.$$

This condition is important to decide on the behavior of the square root showing up in Equation (39). As a consequence,

$${\lambda }_{+}=\frac{3H_{*}}{2}\left(-1+\sqrt{1+\frac{V_1}{3}\frac{1}{\left(1+\frac{2}{3}\omega \right)}{\left(\frac{2}{3H_{*}}\right)}^2}\right)>0.$$

(42)

We concluded that ${\lambda }_{+}$ is a positive real number, so that $\mathrm{Re}\left[{\lambda }_{+}\right]>0$ is always satisfied, and the associated equilibrium point is unstable according to Equation (40). Since the analysis of ${\lambda }_{+}$ already reveals instability, there is no need to proceed to the study of ${\lambda }_{-}$. The direction fields for the autonomous system (37) in the case where $V_1<0$ and $\omega <-\frac{3}{2}$ are shown in Figure 1.

Figure 1 clearly shows that the trajectories diverge from the equilibrium point $\left({\phi }_0,p_{\phi ,0}\right)$$=\left(0,0\right)$, 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 $\omega$ if a reasonable physical interpretation for the dynamics of $\varphi =c^3/G$ is to be established. We will decide on that later on, after we conclude our stability analysis.

Subcase 1.2: Values of $\omega$ satisfying $\omega >-\frac{3}{2}$

Here, we need both $V_1<0$ and $\omega >-\frac{3}{2}$. Under this requirement, the eigenvalues in (39) can be expressed as

$${\lambda }_{\pm }=\frac{3H_{*}}{2}\left(-1\pm \sqrt{1-\frac{\left|V_1\right|}{3}\frac{1}{\left|1+\frac{2}{3}\omega \right|}{\left(\frac{2}{3H_{*}}\right)}^2}\right).$$

(43)

The relevant part for deciding the nature (real or complex) of the eigenvalues ${\lambda }_{\pm }$ is the second term within the square root. If it is greater than 1, then

$${\lambda }_{\pm }=\frac{3H_{*}}{2}\left(-1\pm i\sqrt{\frac{\left|V_1\right|}{3}\frac{1}{\left|1+\frac{2}{3}\omega \right|}{\left(\frac{2}{3H_{*}}\right)}^2-1}\right).$$

and both eigenvalues ${\lambda }_{\pm }$ become complex numbers with negative real parts. Therefore, the equilibrium point is stable, cf. Equation (40).

On the other hand, if

$$0<\frac{\left|V_1\right|}{3}\frac{1}{\left|1+\frac{2}{3}\omega \right|}{\left(\frac{2}{3H_{*}}\right)}^2⩽1,$$

(44)

the square root in (43) is a real number smaller than one, and both eigenvalues ${\lambda }_{\pm }$ are real and negative. This automatically satisfies the criterion $\mathrm{Re}\left[{\lambda }_{\pm }\right]⩽0$ for stability, leading us to conclude, again, that the equilibrium point $\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$ is stable.

The direction fields in the plane $(\phi ,p_{\phi })$ are analyzed in Figure 2. As we can see in the graph, the trajectories in the phase space converge to the equilibrium point $\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$, 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.

Case 2: Positive-valued $V_1$

The stability analysis of the equilibrium point $\left({\phi }_0,p_{\phi ,0}\right)$ in the case $V_1>0$ also depends on the range of values assumed by $\omega$. The two subcases are again $\omega <-\frac{3}{2}$ and $\omega >-\frac{3}{2}$.

Subcase 2.1: Negative-valued $\omega$ with $\omega <-\frac{3}{2}$

The case $\omega <-\frac{3}{2}$ yields

$${\lambda }_{\pm }=\frac{3H_{*}}{2}\left(-1\pm \sqrt{1-\frac{V_1}{3}\frac{1}{\left|\frac{2}{3}\left|\omega \right|-1\right|}{\left(\frac{2}{3H_{*}}\right)}^2}\right),$$

with $\frac{V_1}{3}\frac{1}{\left|\frac{2}{3}\left|\omega \right|-1\right|}>0$. The analysis here follow exactly the same steps as for Subcase 1.2 above. If $\frac{V_1}{3}\frac{1}{\left|\frac{2}{3}\left|\omega \right|-1\right|}{\left(\frac{2}{3H_{*}}\right)}^2>1$, the eigenvalues are complex numbers with negative real parts, meaning that the equilibrium point is stable. On the other hand, if $\frac{V_1}{3}\frac{1}{\left|\frac{2}{3}\left|\omega \right|-1\right|}{\left(\frac{2}{3H_{*}}\right)}^2⩽1$, the two eigenvalues are real negative numbers, implying that the system has a stable fixed point.

The direction fields for the dynamics of $\phi$ when $V_1>0$ and $\omega <-\frac{3}{2}$ are displayed in Figure 3. The trajectories converge to $\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$; for this reason, it is classified as a stable equilibrium point.

Subcase 2.2: Values of $\omega$ satisfying $\omega >-\frac{3}{2}$

In this situation, $V_1>0$ and $\omega >-\frac{3}{2}$, the eigenvalue ${\lambda }_{+}$ is a real positive number. As a consequence, the stability criterion $\mathrm{Re}\left[{\lambda }_{\pm }\right]⩽0$ is violated irrespective of the behavior of the eigenvalue ${\lambda }_{-}$. Hence, the equilibrium point $\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$ is unstable. This conclusion is confirmed by the direction fields in Figure 4 showing trajectories diverging from the equilibrium point.

3.3. Phase Space Analysis with $H=H\left(x^0\right)$ via Lyapunov’s Method

Presently, we generalized the autonomous system studied in Section 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 [83] (and Ref. [79]). We were, thus, compelled to use other methods of dynamical analysis. In particular, Lyapunov’s method—reviewed in Section 3.1.2—proved to be convenient in the present case.

The generalization of the system in Equation (37) that takes into account $H=H\left(x^0\right)$ is

$$\left\{\begin{array}{c}\dot{\phi }=p_{\phi }\hfill \\ {\dot{p}}_{\phi }=-3H\left(x^0\right)p_{\phi }+\frac{V_1}{2\omega +3}\phi \hfill \end{array}\right..$$

(45)

The equilibrium point is the same as in Equation (38), namely:

$${\mathbf{y}}_0=\left(\begin{array}{c}{\phi }_0\\ p_{\phi ,0}\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right).$$

(46)

The (vector) function $\mathbf{f}\left(\mathbf{y}\right)$ that determines the dynamics reads

$$\mathbf{f}\left(\mathbf{y}\right)=\left(\begin{array}{c}{f}_1\\ f_2\end{array}\right)=\left(\begin{array}{cc}0& 1\\ \frac{V_1}{2\omega +3}& -3H\left(x^0\right)\end{array}\right)\left(\begin{array}{c}\phi \\ p_{\phi }\end{array}\right),$$

(47)

where

$$\mathbf{y}=\left(\begin{array}{c}\phi \\ p_{\phi }\end{array}\right),$$

(48)

as is immediately concluded from the comparison of the system (45) with the definition (33).

Now, we propose the following Lyapunov function $V\left(\mathbf{y}\right)$:

$$V\left(\mathbf{y}\right)=V\left(\phi ,p_{\phi }\right)=a{\phi }^2+b{p}_{\phi }^2,$$

(49)

with a and b constants. These constants will be determined in accordance with the Lyapunov stability theorem momentarily.

The function in Equation (49) immediately satisfies Property (i) in Section 3.1.2. Moreover, notice that

$$V\left({\mathbf{y}}_0\right)\to V\left({\phi }_0,p_{\phi 0}\right)=V\left(0,0\right)=0.$$

(50)

This fact is relevant for checking if Requirement (ii) is fulfilled $\forall \mathbf{y}\in U$.

We can compute the first identity in Equation (34) for the specific case of our dynamical system:

$$\begin{array}{cc}\hfill \dot{V}& =\frac{\partial V}{\partial \phi }\dot{\phi }+\frac{\partial V}{\partial p_{\phi }}{\dot{p}}_{\phi }=2a\phi \dot{\phi }+2b{p}_{\phi }{\dot{p}}_{\phi }\hfill \\ \hfill & =2\left(a+\frac{V_1}{2\omega +3}b\right)p_{\phi }\phi -6bH\left(x^0\right)p_{\phi }^2,\hfill \end{array}$$

(51)

where we made use of Equation (45). The constants a and b are arbitrary; we used this freedom to choose

$$\left(a+\frac{V_1}{2\omega +3}b\right)=0,\mathrm{and}b>0.$$

(52)

Accordingly, Equation (51) yields

$$\dot{V}=-6bH\left(x^0\right)p_{\phi }^2<0.$$

(53)

The choice in (52) led us to satisfy Property (iii) in Section 3.1.2. This is so because the Hubble function is always positive for an expanding Universe. Actually, we have $\dot{V}<0$ since $H\left(x^0\right)>0$ and $p_{\phi }^2>0$, for $p_{\phi }\ne 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 Section 3.1.2 would not be satisfied.

Incidentally, now is the time to attempt to verify Property (ii). Plugging Equation (52) into (49),

$$V\left(\phi ,p_{\phi }\right)=\left(-\frac{V_1}{2\omega +3}{\phi }^2+p_{\phi }^2\right)b.$$

(54)

Since $b>0$ and ${\phi }^2>0$ and $p_{\phi }^2>0$ in a neighborhood of the equilibrium point, the decision on the fulfillment of Requirement (ii) depends entirely on the values of parameters $V_1$ and $\omega$. Essentially, the stability will be achieved if

$$\frac{V_1}{2\omega +3}<0.$$

(55)

We identify four possibilities:

(a)

$V_1<0$ and $\omega <-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, then Condition (55) is not satisfied, and the system is not stable;

(b)

$V_1<0$ and $\omega >-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, then Condition (55) is satisfied, and the system is asymptotically stable;

(c)

$V_1>0$ and $\omega <-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, then the system is asymptotically stable;

(d)

$V_1>0$ and $\omega >-\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$, then the system is not stable.

These possibilities are consistent with those unveiled by the linear stability analysis performed in Section 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$. Figure 5 shows the direction fields for different values of the Hubble function $H\left(x^0\right)$ and fixed values of $\omega$ and $V_1$ under Condition (b).

If we restrict our parameter space $\left(\omega ,V_1\right)$ to satisfy the conditions for stability, i.e., Conditions (b) or (c) above, then we verify that the equilibrium point ${\mathbf{y}}_0=\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$ is globally asymptotically stable, since

$$\underset{\sqrt{{\phi }^2+p_{\phi }^2}\to \infty }{lim}V\left(\phi ,p_{\phi }\right)\to \infty ,$$

cf. the Lyapunov Stability Theorem—enunciated below Point (iii) in Section 3.1.2.

Finally, we plot the parameter space of the pair $\left(\omega ,V_1\right)$. The stability of the system depends solely on the values of $V_1$ and $\omega$, both when $H\approx H_{*}$ (Section 3.2) and when $H=H\left(x^0\right)$ (this section). The shaded region in Figure 6 displays the region of the parameter space where the dynamical system is stable.

3.4. 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,\omega >-\frac{3}{2}$ and the case with $V_1>0,\omega <-\frac{3}{2}$. These different ranges of values for $V_1$ lead to two different interpretations. Let us reconsider Equation (8) and include the potential (20) explicitly:

$$\begin{array}{cc}\hfill R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R+\frac{V_1}{2}{g}_{\mu \nu }& =\frac{1}{c}\frac{8\pi }{\varphi }{T}_{\mu \nu }-\frac{\left(V_0+V_2{\varphi }^2\right)}{2}\frac{1}{\varphi }{g}_{\mu \nu }\hfill \\ \hfill & +\frac{\omega }{{\varphi }^2}\left({\nabla }_{\mu }\varphi {\nabla }_{\nu }\varphi -\frac{1}{2}{g}_{\mu \nu }{\nabla }^{\rho }\varphi {\nabla }_{\rho }\varphi \right)\hfill \\ \hfill & +\frac{1}{\varphi }\left({\nabla }_{\mu }{\nabla }_{\nu }\varphi -g_{\mu \nu }\square \varphi \right).\hfill \end{array}$$

(56)

As mentioned at 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 $\varphi$. It may be suggestive to rename this parameter $V_1\leftrightarrow 2\Lambda$, in such a way that the left-hand side of Equation (56) is formally the same as the geometrical side of Einstein field equations. Hence, the condition $V_1>0$ maps to $\Lambda >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 $\Lambda <0$, thus leading to an AdS-like background. (Of course, some conditions apply for achieving a strict dS/AdS background: vacuum, derivatives of the field $\varphi$ as null, etc.)

In both cases, as long as we have the stability of the equilibrium point, the system will evolve so that the trajectories of field $\varphi$ in the phase space $\left(\phi ,p_{\phi }\right)$ converge to $\left({\phi }_0,p_{\phi ,0}\right)=\left(0,0\right)$. When this fixed point is reached, we have simultaneously

$$\left\{\begin{array}{c}\dot{\varphi }\to {\dot{\varphi }}_{\mathrm{eq}}=0\hfill \\ \varphi \to {\varphi }_{\mathrm{eq}}=-2\frac{V_0}{V_1}=\mathrm{constant}\hfill \end{array}\right..$$

(57)

When the system arrives at the equilibrium point, then

$$\frac{\dot{\varphi }}{\varphi }\to \frac{{\dot{\varphi }}_{\mathrm{eq}}}{{\varphi }_{\mathrm{eq}}}=-{\left[\frac{\dot{G}}{G}-3\frac{\dot{c}}{c}\right]}_{\mathrm{eq}}=0,$$

(58)

as implied by Equations (5) and (57). This means that the condition:

$$\frac{\dot{G}}{G}=3\frac{\dot{c}}{c}$$

(59)

will necessarily be satisfied whenever the system is in equilibrium: Equation (1) holds true precisely for the parameter value $\sigma =3$. This value is exact here: the uncertainty of $\sigma =3$ is theoretically zero; it is a natural consequence of the definition of $\varphi$ demanded by basic dimensional analysis and out of the dynamics of $\varphi$ in the phase space.

Another consequence of this dynamics of $\varphi$ is that the system can start with arbitrary initial conditions. The condition $\frac{\dot{G}}{G}=3\frac{\dot{c}}{c}$ may not even be satisfied for this initial condition. It is just a matter of time until the system evolves and converge to Equation (1).

We also point out the appearance of an effective cosmological constant when the dynamics of $\varphi$ relax to the equilibrium point (regardless of the association $V_1\leftrightarrow 2\Lambda$ suggested above). In fact, the dynamics of the scalar field ceases at the equilibrium point, i.e., ${\nabla }_{\mu }\varphi =0$ with $\varphi \to {\varphi }_{\mathrm{eq}}=\mathrm{constant}$; in this situation, the gravitational field Equation (56) approaches

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R+\left(\frac{V_0+V_1{\varphi }_{\mathrm{eq}}+V_2{\varphi }_{\mathrm{eq}}^2}{2{\varphi }_{\mathrm{eq}}}\right)g_{\mu \nu }=\frac{8\pi }{{\varphi }_{\mathrm{eq}}}{T}_{\mu \nu }^{\left(\mathrm{eff}\right)},$$

(60)

where $T_{\mu \nu }^{\left(\mathrm{eff}\right)}\equiv \frac{1}{c}{T}_{\mu \nu }$ is the conserved stress–energy tensor in our scalar–tensor model. The term in the parenthesis above appears as an effective cosmological constant:

$${\Lambda }_0\equiv \left(\frac{V_0+V_1{\varphi }_{\mathrm{eq}}+V_2{\varphi }_{\mathrm{eq}}^2}{2{\varphi }_{\mathrm{eq}}}\right),$$

(61)

and the field $\varphi$, which converges to ${\varphi }_{\mathrm{eq}}$ in Equation (57), can be set as

$${\varphi }_{\mathrm{eq}}=\frac{c_0^3}{G_0},$$

(62)

while Equation (60) assumes the form of GR’s field equation with the 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, Equation (60) assumes the following form at the present time:

$$R_{\mu \nu }-\frac{1}{2}{g}_{\mu \nu }R+{\Lambda }_0{g}_{\mu \nu }=\frac{8\pi G_0}{c_0^4}{T}_{\mu \nu },$$

(63)

which is the familiar form of the Einstein field equations, cf. Equation (2).

It is imperative that the following is crystal clear: the fact that ${\varphi }_{\mathrm{eq}}=c_0^3/G_0=\mathrm{constant}$ does not mean that G and c will be constants after the equilibrium is reached. All that is required thereafter is $\left(\dot{\varphi }/\varphi \right)=-\left(\dot{G}/G-3\dot{c}/c\right)=0$. Accordingly, after the equilibrium is attained, it could be $c=c_{0}f$ and $G=G_0{f}^3$, where $f=f\left(x^0\right)$ is an arbitrary function of the time coordinate: these time-dependent ansatzs for G and c satisfy both the requirements ${\varphi }_{\mathrm{eq}}=c_0^3/G_0$ and $\left(\dot{G}/G-3\dot{c}/c\right)=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 $\frac{\dot{G}}{G}\ne 3\frac{\dot{c}}{c}$ (by taking into account only gravitational field equations, as we have done here), one would have to consider situations out of the equilibrium for $\varphi$. 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 Equation (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 $\omega >\mathrm{40,000}$ [96], while as a dimensionless parameter, it would be expected to be of order $\omega \approx 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 within sight. The system begins with arbitrary initial conditions where $\varphi =c^3/G$ is out of equilibrium. In this phase, the functional form of G and c is not constrained to any particular behavior being interwoven. There are two regions in the parameter space of the pair $\left(\omega ,V_1\right)$ enabling the system to reach a stable critical point, cf. Figure 6. By restricting our model parameters to these regions, one guarantees that the stable equilibrium is reached, $\varphi \to {\varphi }_{\mathrm{eq}}$, and the couplings G and c are forced to evolve respecting $\left(\dot{G}/G-3\dot{c}/c\right)=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 into 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 in Equation (60)—see Section 4.1. We will show that this continuity equation differs from the conventionally used in FLRW cosmology as it contains a term depending on $\dot{c}$. The ensuing VSL model is studied in Section 4.2 assuming radiation and dust matter content. In both cases, the evolution of $a=a\left(x^0\right)$ 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\left(x^0\right)$.

4. Analysis of the Gravity Field Equation after $\mathbf{\varphi }$ Reaches Equilibrium

After $\varphi$ reaches the equilibrium point, the field equation for the gravitational part of our scalar–tensor model, Equation (60), reduces to

$$G_{\phantom{\rho }\nu }^{\rho }=\frac{8\pi }{{\varphi }_{\mathrm{eq}}}\frac{1}{c}{T}_{\phantom{\rho }\nu }^{\rho }-{\Lambda }_0{\delta }_{\phantom{\rho }\nu }^{\rho },$$

(64)

where ${\varphi }_{\mathrm{eq}}=\mathrm{constant}$. However, $c=c\left(x^0\right)$ on the r.h.s. accompanying $T_{\mu \nu }$. In this section, we investigate non-vacuum cosmological solutions stemming from this new feature.

The 00-component and the $ii$-component of Equation (64), under the FLRW metric in (13), read

$$H^2=\frac{8\pi }{3{\varphi }_{\mathrm{eq}}}\frac{1}{c}\epsilon +\frac{{\Lambda }_0}{3}-\frac{k}{a^2}$$

(65)

and

$$\dot{H}+H^2=-\frac{4\pi }{3{\varphi }_{\mathrm{eq}}}\frac{1}{c}\left(\epsilon +3p\right)+\frac{{\Lambda }_0}{3},$$

(66)

which are the Friedmann equation and the acceleration equation, respectively. (Recall that the definition $H=\dot{a}/a$ leads to $\ddot{a}/a=\dot{H}+H^2$.)

4.1. The Covariance of the Effective Stress–Energy Tensor

By taking the covariant divergence of Equation (64) and by using the Bianchi identity $\left({\nabla }_{\rho }{G}_{\phantom{\rho }\nu }^{\rho }=0\right)$ and the metricity condition, we conclude that

$${\nabla }_{\rho }\left(\frac{1}{c}{T}_{\phantom{\rho }\nu }^{\rho }\right)=0.$$

(67)

This is the extended conservation law in our scalar–tensor gravity and justifies the definition of $T_{\mu \nu }^{\left(\mathrm{eff}\right)}=\frac{1}{c}{T}_{\mu \nu }$, which is actually covariantly conserved. The effective stress–energy tensor $T_{\mu \nu }^{\left(\mathrm{eff}\right)}$ is Noether’s current related to the symmetry of (infinitesimal) general coordinate transformations. When we specify Equation (67) for the FLRW metric, we obtain

$$\dot{\epsilon }+3H\left(\epsilon +p\right)=\left(\frac{\dot{c}}{c}\right)\epsilon .$$

(68)

The term on the r.h.s. disappears for $\dot{c}=0$ so that Equation (68) recovers the standard continuity equation of background cosmology in GR in this particular case.

Equations (65), (66), and (68) form the essential set of equations 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:

• An equation of state $p=p\left(\epsilon \right)$;
• A constitutive equation for $c=c\left(x^0\right).$

The EOS is taken as

$$p=w\epsilon ,\mathrm{where}w=\left\{\begin{array}{cc}0,\hfill & \mathrm{dust}\mathrm{matter}\hfill \\ \frac{1}{3},\hfill & \mathrm{radiation}\hfill \end{array}\right..$$

(69)

(Notice the distinction between w—the EOS parameter—and $\omega$—the scalar–tensor constant parameter appearing in Equation (7) as the coefficient of the kinetic term of $\varphi$.) The continuity Equation (68) then becomes:

$$\dot{\epsilon }+\left[3H\left(1+w\right)-\left(\frac{\dot{c}}{c}\right)\right]\epsilon =0.$$

(70)

The first term in the square brackets is the common contribution in the standard GR FLRW cosmology. The second term in the square brackets is the contribution from our scalar–tensor model. The above equation cannot be solved explicitly without a constitutive equation for $c=c\left(x^0\right)$. That is what we will do next.

4.2. Modeling the Varying Speed of Light

In a few models realizing VSL scenarios, it is assumed that $c\left(x^0\right)$ suffered a first-order phase transition at the beginning 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 motivated 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\left(a\right)$, which would be valid after the dynamics of the field $\varphi$ ceases, does not violate the equilibrium condition $\varphi ={\varphi }_{\mathrm{eq}}=\mathrm{constant}$. For instance, if $c\sim 1/a$, it suffices that $G\sim 1/a^3$ for maintaining $\varphi =\frac{c^3}{G}=\mathrm{constant}$.

4.2.1. Speed of Light Scaling as the Inverse of A

Let us assume that

$$c=c_0\left(\frac{a_0}{a}\right).$$

(71)

That is, the speed of light decreases as the Universe increases its size. The consequence of (71) is [97]

$$\frac{\dot{c}}{c}=-H.$$

(72)

This is a reasonable choice that simplifies the continuity equation. In fact, plugging (72) into Equation (70) enables immediate integration to

$$\epsilon ={\epsilon }_0{\left(\frac{a_0}{a}\right)}^{\left(4+3w\right)}.$$

(73)

In particular, $\epsilon \sim a^{-5}$ in a radiation-dominated era ($w=1/3$). This violates the Stefan–Boltzmann law $\epsilon \sim T^4$ because $a\sim T{}^{-1}$, with T representing the temperature. (Of course, the violation of the Stefan–Boltzmann law by radiation in the model $c\sim a^{-1}$ occurs if we insist on a constant $k_B$.) Curiously enough, dust-like matter ($w=0$) does recover the Stefan–Boltzmann law $\epsilon \sim a^{-4}$ in our scalar–tensor model supplemented by the ansatz $c\sim a^{-1}$.

Substituting (71) and (73) into Equation (65) with $k=0$ and ${\varphi }_{\mathrm{eq}}=\left(c_0^3/G_0\right)$ yields

$$H^2=\frac{{\mathsf{\Omega }}_0}{a_0^2}{\left(\frac{a_0}{a}\right)}^{3\left(1+w\right)}+\frac{{\Lambda }_0}{3}$$

(74)

under the definitions:

$${\epsilon }_{c,0}\equiv \frac{1}{a_0^2}\frac{3c_0^4}{8\pi G_0}\mathrm{and}{\mathsf{\Omega }}_0=\frac{{\epsilon }_0}{{\epsilon }_{c,0}},$$

(75)

where $\left[{\epsilon }_{c,0}\right]=\left(\mathrm{energy}\right)/{\left(\mathrm{length}\right)}^3$ and ${\mathsf{\Omega }}_0$ is dimensionless. Equation (74) integrates to

$$\left(\frac{a}{a_0}\right){\left[1+\sqrt{1+\frac{{\mathsf{\Omega }}_0}{\left(a_0^2\frac{{\Lambda }_0}{3}\right)}{\left(\frac{a_0}{a}\right)}^{3\left(1+w\right)}}\right]}^{2/3\left(1+w\right)}={\left[1+\sqrt{1+\frac{{\mathsf{\Omega }}_0}{\left(a_0^2\frac{{\Lambda }_0}{3}\right)}}\right]}^{2/3\left(1+w\right)}exp\left[\sqrt{\frac{{\Lambda }_0}{3}}\left(x^0-x_0^0\right)\right].$$

(76)

In Equation (76), $x_0^0$ is the value of the “time” coordinate $x^0$ today. Equation (76) is the exact solution for the scale factor.

If ${\Lambda }_0{a}_0^2\gg {\mathsf{\Omega }}_0$ (the cosmological constant dominates over the matter–energy content), the above solution reduces to

$$a\simeq a_{0}exp\left[\sqrt{\frac{{\Lambda }_0}{3}}\left(x^0-x_0^0\right)\right]\left({\Lambda }_0{a}_0^2\gg {\mathsf{\Omega }}_0\right).$$

(77)

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 ${\Lambda }_0{a}_0^2\gg {\mathsf{\Omega }}_0$. Therefore, our scalar–tensor model with $c\sim a^{-1}$ is a promising candidate for fitting the observational data. It remains to be investigated how the values of its parameters $\left({\mathsf{\Omega }}_0,{\Lambda }_0\right)$ 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 Equation (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\sim a^{-1}$ in which the speed of light decreases as the Universe expands. How dependent is this scenario 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\left(a\right)$ that is different from that in Equation (71). We 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].

4.2.2. Speed of Light Scaling as the Exponential of a

Gupta’s ansatz for an increasing speed of light has the form [66]:

$$c=c_{0}exp\left[{\left(\frac{a}{a_0}\right)}^{\alpha }-1\right],$$

(78)

with $\alpha$ a positive constant of order one. Therefore,

$$\frac{\dot{c}}{c}=\alpha {\left(\frac{a}{a_0}\right)}^{\alpha }H.$$

(79)

Inserting (79) into (70),

$$\epsilon ={\epsilon }_0{\left(\frac{a_0}{a}\right)}^{3\left(1+w\right)}\left(\frac{c}{c_0}\right).$$

(80)

For the standard non-varying speed of light picture, $c=c_0$, the equation above recovers the traditional result from the background FLRW cosmology of GR. This means that the Stefan–Boltzmann law $\epsilon \sim a^{-4}$ is recovered in the radiation era $\left(w=1/3\right)$, except for the factor $\left(c/c_0\right)$. This is different from the case in the previous subsection, where $c\sim a^{-1}$ and $\epsilon \sim a^{-4}$ were realized only in the presence of dust-like matter.

Now, Equation (80) can be used in the first Friedmann Equation (65) to yield the scale factor. Indeed, taking $k=0$ and using the definitions of ${\epsilon }_{c,0}$ and ${\mathsf{\Omega }}_0$ in Equation (75) result in

$$H^2=\frac{{\mathsf{\Omega }}_0}{a_0^2}{\left(\frac{a_0}{a}\right)}^{3\left(1+w\right)}+\frac{{\Lambda }_0}{3}.$$

(81)

Surprisingly, this is exactly the same equation for H derived in the previous model—see Equation (74). Therefore, the solution devised in Section 4.2.1 for the model with $c\sim a^{-1}$ is also valid in the present case, namely: Equation (76) leading to the same de Sitter-type limiting case of Equation (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\left(x^0\right)$, just like one admits different functions of $a=a\left(x^0\right)$ 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 $\varphi$ reaches its equilibrium point ${\varphi }_{\mathrm{eq}}$; after all, the whole formalism developed in the present Section 4 assumes the validity of Equation (64), which is the BD-like field equation for the gravitational field after $\varphi$ approaches ${\varphi }_{\mathrm{eq}}$.

It should be emphasized that the model in this subsection (Section 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], and the faint young Sun problem [71].

5. 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 were included in the scalar sector of the model through the field $\varphi =c^3/G$.

The field equations for $g_{\mu \nu }$ and $\varphi$ were built and specified for the homogeneous and isotropic cosmological background. The dynamics of our field $\varphi$ was then analyzed in the phase space under some working hypotheses. We took the curvature parameter of the space sector k as null. We used the first three terms in the series expansion of the potential $V\left(\varphi \right)$. We assumed that the dynamics of $\varphi$ happens in a radiation era. We then studied the evolution of $\varphi$ via linear stability theory (reviewed in Section 3.1) assuming that the dynamical period of $\varphi$ is short enough to ensure that the Hubble function is approximately constant therein (Section 3.2). Finally, the dynamical system analysis was extended to allow for a time-dependent Hubble parameter; this required the use of Lyapunov’s method (Section 3.3).

The phase space analysis showed that two sets of conditions upon the free parameters of the model led to trajectories converging to a globally asymptotically stable equilibrium point. These conditions were: (i) $V_1<0,\omega >-\frac{3}{2}$ and (ii) $V_1>0,\omega <-\frac{3}{2}$, where $V_1$ is the constant coefficient of the linear term in the $V\left(\varphi \right)$ series and $\omega$ is the constant coefficient of the kinetic term of $\varphi$ in the model’s action. The physical interpretation of $V_1$ and $\omega$ is commented on momentarily.

We saw in Section 3.4 that $V_1$ could be understood as a type of cosmological constant during the dynamical phase of $\varphi$; therefore, its positivity (or negativity) could impact this version of $\Lambda$. Moreover, $V_1$ appears in the effective cosmological constant ${\Lambda }_0$, which turns up at the end of the dynamical evolution of $\varphi$, 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 ${\Lambda }_0$ includes the other coefficients of the potential as well: $V_0$ is the constant term in the $V\left(\varphi \right)$ series, and $V_2$ is the coefficient of the quadratic term of V in $\varphi$. Our treatment demands $V_0\ne 0$, but otherwise unconstrained, while $V_2$ is totally unconstrained. The equilibrium value of the field, ${\varphi }_{\mathrm{eq}}=c_0^3/G_0$, also enters the definition of ${\Lambda }_0$, which leads to the interesting conclusion that the cosmological constant could depend on $\left(G_0,c_0\right)$, the (constant) equilibrium values of the gravitational coupling and of the speed of light.

The sign of the parameter $\omega$ is related to the physical nature of the scalar field $\varphi$. Indeed, if $\frac{\omega }{\varphi }<0$, the sign of the kinetic term for $\varphi$ in the action indicates that $\varphi$ is a ghost field. This means that its Hamiltonian is unbounded from below and it would be a source of negative energy. The global sign of the potential $V\left(\varphi \right)$ is also important: together with the sign of the kinetic term, it determines if $\varphi$ could also be an unphysical tachyon field breaking causality. In the vicinity of the equilibrium point, $\varphi >0$, so that $\omega <0$ would characterize the field as a ghost field—see, e.g., page 2 of [98] and the references therein. This all means that we have to impose $\omega >0$ for a physically meaningful behavior of $\varphi$.

The big picture in our development is the following: After evolving to the equilibrium stable point, the field $\varphi$ reaches it equilibrium value ${\varphi }_{\mathrm{eq}}=c_0^3/G_0$ and stays there. Then, the dynamics of $\varphi$ ceases, and two things happen: (i) the couplings G and c that are covarying within $\varphi$ assume the dependence $G/G_0={\left(c/c_0\right)}^3$, where $G_0$ and $c_0$ could be interpreted as the numbers that we use as fundamental constants today; (ii) the gravitational field equation of our scalar–tensor model, involving both $g_{\mu \nu }$ and $\varphi$, degenerates to a field equation that is formally the same as the Einstein field equation of general relativity, but bears a time-dependent factor $\left(1/c\right)$ alongside the ordinary energy momentum tensor. At the present time $c=c_0$, and the field equation recovers exactly the 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 [101] by Faraoni and Franconnet.

When the dynamics of $\varphi$ ends, $\dot{\varphi }=0$, and the covarying G and c are forced to obey the relation $\left(\dot{G}/G\right)=\sigma \left(\dot{c}/c\right)$ with $\sigma =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 $\dot{\varphi }=0$. This fact was pointed out in [65,69,70], in relation to orbital timing, the 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 like those in BD theory, it must be underscored 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 ensuring the general covariance of the gravity equation. This was discussed also in Section 4, where we built the Friedmann equations for the cosmology within our scalar–tensor gravity. These equations were solved for two particular ansatzs of $c=c\left(x^0\right)$: (i) a speed of light that decreases as a increases and (ii) a c that increases as the Universe expands. It was shown that both ansatzs 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 relaxing some of the hypotheses adopted in this manuscript. One possibility is to let the dynamics of $\varphi$ evolve in an era dominated by matter, for which $\left(ϵ-3p\right)\ne 0$ and the coupling of $c=c\left(x^0\right)$ with the matter Lagrangian comes into effect. The eventual interrelation between the matter part of the action and $c=c\left(x^0\right)$ brings forth the possibility of c being an independent scalar field. In this context, one would have to take $G=G\left(x^{\mu }\right)$ and $c=c\left(x^{\mu }\right)$ as scalar fields evolving on their own (not enclosed within the field $\varphi$). To put it another way, in the stability analysis performed in this paper, we have one scalar degree of freedom, incarnated in $\varphi$; 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 $\left\{G,c,\Lambda \right\}$—with the first two forming the field $\varphi$. Other works enlarge the set of possible varying fundamental constants to $\left\{G,c,\Lambda ,ħ,k_B\right\}$ and perform some phenomenological modeling in astrophysics and cosmology [65,102]. Pursuing the fundamentals surrounding an eventual variation of all five couplings from a field theory perspective is something that may be performed in the future.