# Gauss-Bonnet Dark Energy and the Speed of Gravitational Waves

José Jaime Terente Díaz  Departamento de Física Teórica, Universidad
Complutense de Madrid, E-28040 Madrid, Spain  Konstantinos Dimopoulos
Consortium for Fundamental Physics, Physics Department, Lancaster University,
Lancaster LA1 4YB, UK  Mindaugas Karčiauskas  Departamento de Física Teórica,
Universidad Complutense de Madrid, E-28040 Madrid, Spain  Antonio Racioppi
National Institute of Chemical Physics and Biophysics, Rävala 10, 10143
Tallinn, Estonia

###### Abstract

Gauss-Bonnet Dark Energy has been a popular model to explain the accelerated
expansion of the Universe. Quite generically it also predicts the speed of
gravitational waves $c_{\mathrm{GW}}$ to be different from the speed of light.
This fact alone led some authors to exclude such models in view of the new
tight observational constraints on $c_{\mathrm{GW}}$. However, the behaviour
of $c_{\mathrm{GW}}$ depends on the choice of the Gauss-Bonnet (GB) coupling
function. It is possible to construct models where $c_{\mathrm{GW}}$ is always
equal to the speed of light. More generally, $c_{\mathrm{GW}}$ is a time
dependent function with instances where both speeds coincide. Nevertheless, we
observe that the bound on $c_{\mathrm{GW}}$ excludes scenarios where the GB
term directly affects the expansion of the Universe, even if the constraint on
the variation of the coupling function does not appear to be strong. We
perform the dynamical systems analysis to see if the expansion of the Universe
could be affected indirectly by modulating the behaviour of the scalar field,
which modulates the GB coupling. It is shown that either the bounds on
$c_{\mathrm{GW}}$ are violated by many orders of magnitude, or it might be
very difficult to find models that are consistent with other cosmological
observations.

##  I Introduction

The detection of gravitational waves (GW) [1, 2, 3] opens a new window to
observe and measure the Universe. Most directly, it enables testing General
Relativity (GR) in regimes that were not accessible before and constrain
possible modifications of the laws of gravity. They also provide new ways to
test Dark Energy (DE) models. Many of such models rely on gravity
modifications and therefore are subject to such constraints.

A very clear demonstration is provided in Ref. [4]. A lucky coincidence of
being able to detect GW emitted by the merger of two neutron stars as well as
the electromagnetic counterpart of this event made it possible to put very
stringent constraints on the speed of GW, $c_{\mathrm{GW}}$. The delay between
arrival times of GW and $\gamma$-rays led to the bound

$\displaystyle\left|\alpha_{\mathrm{T}}\right|$ $\displaystyle<$
$\displaystyle 10^{-15}\,,$ (1)

where $\alpha_{\mathrm{T}}$ parametrises the deviation of $c_{\mathrm{GW}}$
from the speed of light

$\displaystyle\alpha_{\mathrm{T}}$ $\displaystyle\equiv$ $\displaystyle
c_{\mathrm{GW}}^{2}-1\>,$ (2)

in natural units, where $c=\hbar=1$.

Many classes of modified gravity theories predict $\alpha_{\mathrm{T}}\neq 0$.
The constraints on $\alpha_{\mathrm{T}}$ in Ref. [4] excluded a lot of well
motivated and otherwise attractive models and considerably narrowed down the
space of available modifications [5, 6].

Among the excluded models – it is claimed in Ref. [5] – is the Gauss-Bonnet
Dark Energy (GBDE) one. This model has many attractive features. The Gauss-
Bonnet term itself is a unique combination of curvature terms squared

$\mathcal{G}\equiv
R^{2}-4R_{\mu\nu}R^{\mu\nu}+R_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}\,,$ (3)

where $R$, $R_{\mu\nu}$ and $R_{\hphantom{\rho}\sigma\mu\nu}^{\rho}$ are the
Ricci scalar, tensor and Riemann tensor respectively. Nevertheless, this
combination leads to metric tensor equations of motion that are second order.
The Gauss-Bonnet (GB) term is quite ubiquitous in actions of low-energy
effective string theory, be it at tree or one loop level [7, 8, 9, 10, 11, 12,
13]. The corresponding modification can be written as $\xi\mathcal{G}$ term in
the Lagrangian, where $\xi$ is the GB coupling. If the latter is a constant,
the GB term is a surface term and can be integrated out (although it can still
be important for other aspects of the theory, such as regularization [14,
15]). However, on quite generic grounds, one might expect that the GB term
also couples to scalar fields of the theory, such as moduli or dilaton fields,
making $\xi$ field dependent.

The possibility of explaining DE with the GB term was first investigated in
Refs. [16, 17, 18, 19, 20]. One of the attractive features of such models is
that they provide the means to safely cross the phantom divide, that is, enter
the regime where the DE equation of state is $w<-1$, without instabilities.
The best fit value of $w$ is smaller than $-1$ [21]. If such $w$ is associated
with a scalar field, it leads to many instabilities and such a model is likely
excluded by observations [22]. On the other hand, modifications of gravity in
GBDE in some parameter range allow for $w<-1$. In this model $w$ is time
dependent; it briefly dips below $-1$ before settling on $w=-1$ [18]. Hence,
it can accommodate this low value without leading to contradictions.

Quite generically GBDE predicts $\alpha_{\mathrm{T}}\neq 0$. This fact alone
led the authors of Ref. [5] to claim that GBDE is ruled out. Here we would
like to point out that $c_{\mathrm{GW}}$, predicted by GBDE, is not a
constant. Moreover, the constraint in eq. (1) is an upper bound which is
applicable only at the very latest stages of the evolution of the Universe.
Hence, to assess the implications of these constraints for GBDE we need to
study it more carefully.

In this work we use the dynamical systems analysis to look for viable models
of GBDE and compute the evolution of $\alpha_{\mathrm{T}}$ parameter. The
crucial quantity in such models is the GB coupling $\xi\left(\phi\right)$. It
determines the dynamics of the universe as well as the evolution of the
$\alpha_{\mathrm{T}}$ parameter. Applying the bound in eq. (1) to the
variation of $\xi\left(\phi\right)$ and the rate of its variation, we find
that the constraints appear weak. Nevertheless, the bound in eq. (1) prevents
the GB term from affecting the expansion of the Universe directly. The
remaining possibility is for this term to affect the expansion indirectly, by
modifying the behaviour of the scalar field. To investigate this issue we
apply the dynamical systems analysis. We also apply this analysis to the case
where $\alpha_{\mathrm{T}}=0$ by construction, which is allowed by the model.

In Section II we introduce the model, derive dynamical equations and show the
bounds on $\xi\left(\phi\right)$ that follow from eq. (1). In Section III we
assume an exponential potential $V\left(\phi\right)$ and write the dynamical
equations in terms of dimensionless variables, which are used in the following
sections. The dynamical systems analysis is applied to models with the
exponential GB coupling $\xi\left(\phi\right)$ in Section IV and it is applied
to the linear function $\xi\left(\phi\right)$ in Section V. The case of
$\alpha_{\mathrm{T}}=0$ is studied in Section VI.

##  II Scalar-Gauss-Bonnet Dark Energy and Constraints on $c_{\mathrm{GW}}$

We start with the scalar-Gauss-Bonnet action

$\displaystyle S$ $\displaystyle=$
$\displaystyle\int\mathrm{d}^{4}x\sqrt{-g}\left[\frac{1}{2}m_{\mathrm{Pl}}^{2}R+\xi\left(\phi\right)\mathcal{G}-\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V\left(\phi\right)+\mathcal{L}_{\mathrm{m}}\right]\,,$
(4)

where $\mathcal{G}$ is defined in eq. (3). For brevity we address to the above
action as the Gauss-Bonnet (GB) action in this work.
$\mathcal{L}_{\mathrm{m}}$ is the Lagrangian of the matter sector. If $\xi$ is
constant, the GB term is a total derivative and does not affect the dynamics
of the system. We assume the background spacetime to be homogeneous, isotropic
and flat, described by the FRW metric
$g_{\mu\nu}=\mathrm{diag}\left[-1,a^{2}\left(t\right),a^{2}\left(t\right),a^{2}\left(t\right)\right]$,
where $t$ is the cosmic time and $a$ is the scale factor.

In this model the speed of tensor mode propagation is determined by the rate
of change of the coupling function $\xi\left(\phi\right)$ [23, 24, 25, 26]
(see also Refs. [27, 28])

$\displaystyle\alpha_{\mathrm{T}}$ $\displaystyle=$
$\displaystyle\frac{8\left(\ddot{\xi}-\dot{\xi}H\right)}{m_{\mathrm{Pl}}^{2}+8\dot{\xi}H}\,.$
(5)

It is clear from this expression that the constraint in eq. (1) can be
satisfied if one of the following two conditions is fulfilled. The first
option is to choose the coupling function $\xi\left(\phi\right)$ such that

$\displaystyle\ddot{\xi}$ $\displaystyle=$ $\displaystyle H\dot{\xi}\,.$ (6)

This choice is discussed in refs. [26, 29] (other related references can be
found in these articles) in the context of inflation. However, there is
another, more generic possibility. We notice in eq. (5) that $\alpha_{T}$ is
suppressed by the Planck mass. Therefore, as long as the conditions
$\left|\ddot{\xi}\right|/m_{\mathrm{Pl}}^{2},\>H\left|\dot{\xi}\right|/m_{\mathrm{Pl}}^{2}<10^{-15}$
are satisfied, the GB action in eq. (4) is compatible with the constraints on
the speed of gravitational waves. We can write these conditions in a more
useful way

$\displaystyle\frac{\left|\ddot{\xi}\right|}{H^{2}},\,\frac{\left|\dot{\xi}\right|}{H}$
$\displaystyle<$ $\displaystyle
10^{-15}\left(\frac{m_{\mathrm{Pl}}}{H}\right)^{2}\,,$ (7)

which emphasises the change of the coupling function and the rate of change of
this coupling over one Hubble time. If this condition is to be imposed on
inflation models, where typically $H\sim 10^{-5}m_{\mathrm{Pl}}$, this bound
can be tight. In that case, to limit $\alpha_{T}$ within the allowed range, it
is better to look for models that satisfy the condition in eq. (6). However,
the constraints from the observations of GRB170817A apply not to the early
Universe, but the late one, when the Hubble parameter is more than fifty
orders of magnitude smaller than during inflation. Indeed, plugging
$H_{0}^{2}/m_{\mathrm{Pl}}^{2}\sim 10^{-120}$ into eq. (7) we find

$\displaystyle\frac{\left|\ddot{\xi}\right|}{H_{0}^{2}},\,\frac{\left|\dot{\xi}\right|}{H_{0}}$
$\displaystyle<$ $\displaystyle 10^{105}\,.$ (8)

That is, $\xi$ and $\dot{\xi}$ need to vary by more than 100 orders of
magnitude, $\Delta\xi\>,\Delta\dot{\xi}<10^{105}$, over the age of the
Universe to violate the bound. This appears to demonstrate that the constraint
on $c_{\mathrm{GW}}$ is exceptionally weak and might give hope that GBDE
models remain viable. Unfortunately, as we show in this work, at least for
simple functions $\xi\left(\phi\right)$, this turns out not to be the case.

To understand the implications of eq. (1) for Gauss-Bonnet Dark Energy (GBDE)
models better, let us first write the homogeneous dynamical equations in the
FRW background as [16, 30]

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{\rho_{\phi}+\rho_{\mathrm{m}}}{3\left(m_{\mathrm{Pl}}^{2}+8H\dot{\xi}\right)}\>,$
(9) $\displaystyle\dot{H}$ $\displaystyle=$
$\displaystyle-\frac{\rho_{\phi}+P_{\phi}+\rho_{\mathrm{m}}+P_{\mathrm{m}}+8H^{2}\left(\ddot{\xi}-\dot{\xi}H\right)}{2\left(m_{\mathrm{Pl}}^{2}+8H\dot{\xi}\right)}\>,$
(10)

where $\rho_{\phi}$ and $P_{\phi}$ are the energy and pressure densities of
the homogeneous scalar field $\phi$ respectively. They are defined to be

$\displaystyle\rho_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\dot{\phi}^{2}+V\left(\phi\right)\>,$ (11)
$\displaystyle P_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\dot{\phi}^{2}-V\left(\phi\right)\>.$ (12)

Similarly $\rho_{\mathrm{m}}$ and $P_{\mathrm{m}}$ are the energy and pressure
densities of the matter field.

The acceleration of spatial slices can be parametrised using the Hubble flow
parameter

$\displaystyle\epsilon_{H}$ $\displaystyle\equiv$
$\displaystyle-\frac{\dot{H}}{H^{2}}\>.$ (13)

Alternatively, it is common to use the deceleration parameter
$q\equiv\epsilon_{H}-1$ for this purpose. The spatial slices expand in an
accelerating fashion if $\epsilon_{H}<1$ ($q<0$).

Plugging eqs. (9) and (10) into eq. (13) we can write

$\displaystyle\epsilon_{H}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\left(1+\frac{P_{\phi}+P_{\mathrm{m}}}{\rho_{\phi}+\rho_{\mathrm{m}}}\right)+\frac{1}{2}\alpha_{T}\>,$
(14)

where we also made use of eq. (5). At the present epoch $\epsilon_{H}\simeq
0.5$. Thus, in view of eq. (1) we see that the last term must be negligible.
This rules out any direct effect of the Gauss-Bonnet term to the expansion of
the Universe. Neglecting that last term, we arrive at the expression which can
also be obtained in a typical, General Relativistic quintessence models [31,
32].

But even if observations exclude the scenario where the GB term affects the
expansion of the Universe directly, there remains a possibility that it does
so indirectly, by modifying the behaviour of the scalar field $\phi$. As we
will see next, such a possibility is also excluded, at least for an
exponential potential $V\left(\phi\right)$.

##  III The Dynamical System

Equations (9) and (10) can be supplemented with dynamical equations governing
the evolution of the $\phi$ field and $\rho_{\mathrm{m}}$

$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}$ $\displaystyle=$
$\displaystyle 24\xi_{,\phi}\left(\dot{H}+H^{2}\right)H^{2}\>,$ (15)
$\displaystyle\dot{\rho}_{\mathrm{m}}+3H\rho_{\mathrm{m}}\left(1+w_{\mathrm{m}}\right)$
$\displaystyle=$ $\displaystyle 0\>,$ (16)

where $w_{\mathrm{m}}\equiv P_{\mathrm{m}}/\rho_{\mathrm{m}}$ is the
barotropic parameter of the matter component. And we assume a matter fluid
with $0\leq w_{\mathrm{m}}<1$.

To analyse the generic behaviour of this dynamical system, it is convenient to
normalise the dynamical degrees of freedom and write them in a dimensionless
form, such as

$x\equiv\frac{\phi^{\prime}}{\sqrt{6}m_{\mathrm{Pl}}}\,,\quad
y\equiv\frac{\sqrt{V}}{\sqrt{3}m_{\mathrm{Pl}}H}\>,\quad u\equiv
4\sqrt{6}\frac{H^{2}\xi_{,\phi}}{m_{\mathrm{Pl}}}\,,\text{ and
}z\equiv\frac{\sqrt{\rho_{\mathrm{m}}}}{\sqrt{3}m_{\mathrm{Pl}}H}\,.$ (17)

The prime in the definition of $x$ and the equations below denotes the
derivatives with respect to the e-fold number

$\displaystyle N$ $\displaystyle\equiv$ $\displaystyle\ln a\,,$ (18)

where we normalised $a$ such that $a=1$ today.

The definitions in eq. (17) are particularly useful if the scalar field
potential is an exponential function

$\displaystyle V$ $\displaystyle=$ $\displaystyle
V_{0}\mathrm{e}^{-\lambda\phi/m_{\mathrm{Pl}}}\>,$ (19)

where we take $\lambda>0$ to be a constant. We will always use the above
ansatz in this work. In that case the dynamical equations are self-similar and
the explicit dependence on the Hubble parameter drops out of those equations.
In particular, eqs. (15) and (16) can be written as

$\displaystyle x^{\prime}$ $\displaystyle=$
$\displaystyle\left(\epsilon_{H}-3\right)x+\sqrt{\frac{3}{2}}\lambda
y^{2}+u\left(1-\epsilon_{H}\right)\>,$ (20) $\displaystyle y^{\prime}$
$\displaystyle=$ $\displaystyle\left(\epsilon_{H}-\sqrt{\frac{3}{2}}\lambda
x\right)y\>,$ (21) $\displaystyle z^{\prime}$ $\displaystyle=$
$\displaystyle\left[\epsilon_{H}-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\right]z\>,$
(22)

where

$\displaystyle\epsilon_{H}$ $\displaystyle=$
$\displaystyle\left[3x^{2}+\frac{3}{2}\left(1+w_{\mathrm{m}}\right)z^{2}+\left(ux\right)^{\prime}-ux\right]\frac{1}{1+ux}$
(23)

and the constraint in eq. (9) becomes

$\displaystyle 1$ $\displaystyle=$ $\displaystyle x^{2}+y^{2}+z^{2}-2ux\>.$
(24)

A few comments about these equations are in order. First, notice the important
difference from analogous equations of models in General Relativity (see e.g.
[33]). In those models the dimensionless variables $y$ and $z$ are constrained
within the range $\left[0;1\right]$ and $x\in\left[-1;1\right]$. In the case
of the GB model such a restriction does not apply. Due to _a priori_
undetermined sign of $u$, maximum values of $\left|x\right|$, $y$ and $z$ are
not limited to $1$. Second, it might appear that eq. (23) is inconsistent
because it diverges as $ux\rightarrow-1$. But this impression is wrong because
the constraint equation (24) forbids such values.

When doing dynamical analysis of this system, it is convenient to use the
equation for $u$ too. Taking the derivative of the expression in eq. (17) we
find

$\displaystyle u^{\prime}$ $\displaystyle=$
$\displaystyle-2\epsilon_{H}u+24H^{2}\xi_{,\phi\phi}x\>.$ (25)

##  IV The Exponential Gauss-Bonnet Coupling

###  IV.1 Dynamics

To understand the qualitative behaviour of this dynamical system, we find its
fixed points and investigate their stability. The fixed points are defined as
points, or regions, in the phase space where
$x^{\prime}=y^{\prime}=z^{\prime}=0$ is satisfied. First, we study the case of
the exponential GB function, given by

$\xi=\xi_{0}\mathrm{e}^{\kappa\phi/m_{\mathrm{Pl}}}\,,$ (26)

which allows us to write eq. (25) as

$\displaystyle u^{\prime}$ $\displaystyle=$ $\displaystyle\left(\sqrt{6}\kappa
x-2\epsilon_{H}\right)u\>.$ (27)

The computation details for finding fixed points are provided in the Appendix
and the results are summarised in Table 1.

| $x_{\mathrm{c}}$ | $y_{\mathrm{c}}$ | $u_{\mathrm{c}}$ | $z_{\mathrm{c}}$ | Expansion rate $\epsilon_{H\mathrm{c}}$ | If $\kappa$ is | $\epsilon_{H\mathrm{c}}u_{\mathrm{c}}$
---|---|---|---|---|---|---|---
M | 0 | 0 | 0 | 1 | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$ | any $\kappa$ | 0
K± | $\pm 1$ | 0 | 0 | 0 | $3$
I | $\frac{\lambda}{\sqrt{6}}$ | $\sqrt{1-\frac{\lambda^{2}}{6}}$ | 0 | 0 | $\frac{1}{2}\lambda^{2}$
Sc | $\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\lambda}$ | $\frac{\sqrt{\frac{3}{2}\left(1-w_{\mathrm{m}}^{2}\right)}}{\lambda}$ | 0 | $\sqrt{1-\frac{3\left(1+w_{\mathrm{m}}\right)}{\lambda^{2}}}$ | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$
dS | 0 | $1$ | $-\sqrt{\frac{3}{2}}\lambda$ | 0 | 0
S2 | $\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\kappa}$ | 0 | $\sqrt{\frac{3}{2}}\frac{3\left(w_{\mathrm{m}}^{2}-1\right)}{\kappa\left(3w_{\mathrm{m}}+1\right)}$ | $\sqrt{1+\frac{3\left(3w_{\mathrm{m}}-7\right)\left(1+w_{\mathrm{m}}\right)^{2}}{2\kappa^{2}\left(3w_{\mathrm{m}}+1\right)}}$ | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$ | $\kappa\neq\lambda$ | $\neq 0$
G | $\beta$ | 0 | $\frac{\beta^{2}-1}{2\beta}$ | 0 | $\frac{5\beta^{2}+1}{\beta^{2}+1}=\sqrt{\frac{3}{2}}\kappa\beta$
S3 | $\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\lambda}$ | $\frac{\sqrt{\frac{3}{2}\left(1-w_{\mathrm{m}}^{2}\right)+\frac{\lambda}{\sqrt{6}}\left(1+3w_{\mathrm{m}}\right)u_{\mathrm{c}}}}{\lambda}$ | $u_{\mathrm{c}}$ | $\sqrt{1-\frac{3\left(1+w_{\mathrm{m}}\right)-\frac{\lambda}{\sqrt{6}}\left(5+3w_{\mathrm{m}}\right)u_{\mathrm{c}}}{\lambda^{2}}}$ | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$ | $\kappa=\lambda$
IV | $x_{\mathrm{c}}$ | $\sqrt{1-x_{\mathrm{c}}^{2}+2x_{\mathrm{c}}\cdot\frac{3x_{\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda}{1+\sqrt{\frac{3}{2}}\lambda x_{\mathrm{c}}}}$ | $\frac{3x_{\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda}{1+\sqrt{\frac{3}{2}}\lambda x_{\mathrm{c}}}$ | 0 | $\sqrt{\frac{3}{2}}\lambda x_{\mathrm{c}}$
Table 1: The values of $x$, $y$, $u$ and $z$ parameters at the fixed points
and curves. The conditions for the existence of these fixed points and curves
can be derived from the fact that all $x_{\mathrm{c}}$, $y_{\mathrm{c}}$,
$u_{\mathrm{c}}$ and $z_{\mathrm{c}}$ values must be real. Parameters
$\lambda$, $\kappa$ and $\beta$ are defined in eqs. (19), (26) and (75)
respectively. Both $\lambda$ and $\kappa$ are strictly positive.

Looking at this table, notice that all the fixed points with
$u_{\mathrm{c}}=0$ coincide with the ones analysed in ref. [33], as expected.
In this reference the authors analyse a quintessence model within the theory
of GR and the exponential scalar field potential. In other words, all the
fixed points that are present in a similar setup in GR, they are also present
in GB models. However, even if the fixed points coincide, the presence of the
GB term might change their stability, as we will show below.

The fixed point M in Table 1 corresponds to the case where the scalar field is
diluted and only the matter field remains. The K± fixed points correspond to
kination, where the universe is dominated by the kinetic energy of the scalar
field. In the case of I and $\lambda<\sqrt{2}$ this fixed point represents the
power law inflation [34]. In the scaling fixed point (Sc) the evolution of the
scalar field adjusts to mimic the behaviour of the matter field. Therefore,
the expansion rate of the Universe is given by
$\epsilon_{H}=\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$.

The GB term introduces two more scaling solutions: the fixed point S2, which
exists for $\kappa\neq\lambda$, and the fixed curve S3, which exists if
$\xi\left(\phi\right)V\left(\phi\right)=\mathrm{constant}$. For our purpose,
the most interesting new fixed point is the de Sitter one (dS), where
$\dot{H}=0$. This fixed point is very robust, and exists for large variety of
$\xi\left(\phi\right)$ functions.

Some discussion of DE models with $\kappa=\lambda$ and various solutions were
provided in Ref. [16]. Some dynamical analysis with $\kappa\neq\lambda$ was
performed in Ref. [18] (see also [35]). Here we modify and extend the analysis
to make it more generic. In this section, we take $\kappa\neq\lambda$.

Often, the stability of fixed points can be determined by taking a linear
perturbation of equations (20)–(25) around those points. In the case of an
exponential GB function in eq. (26) those linear equations can be written as

$\displaystyle\delta x^{\prime}$ $\displaystyle=$
$\displaystyle\left(\epsilon_{H\mathrm{c}}-3\right)\delta x+\sqrt{6}\lambda
y_{\mathrm{c}}\delta y+\left(1-\epsilon_{H\mathrm{c}}\right)\delta
u+\left(x_{\mathrm{c}}-u_{\mathrm{c}}\right)\delta\epsilon_{H}\,,$ (28)
$\displaystyle\delta y^{\prime}$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{3}{2}}\lambda y_{\mathrm{c}}\delta
x+\left(\epsilon_{H\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda
x_{\mathrm{c}}\right)\delta y+y_{\mathrm{c}}\delta\epsilon_{H}\,,$ (29)
$\displaystyle\delta u^{\prime}$ $\displaystyle=$ $\displaystyle\sqrt{6}\kappa
u_{\mathrm{c}}\delta x+\left(\sqrt{6}\kappa
x_{\mathrm{c}}-2\epsilon_{H\mathrm{c}}\right)\delta
u-2u_{\mathrm{c}}\delta\epsilon_{H}\,,$ (30) $\displaystyle\delta z^{\prime}$
$\displaystyle=$
$\displaystyle\left[\epsilon_{H\mathrm{c}}-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\right]\delta
z+z_{\mathrm{c}}\delta\epsilon_{H}\,,$ (31)

where

$\displaystyle\delta\epsilon_{H}$ $\displaystyle=$
$\displaystyle\frac{u_{\mathrm{c}}\delta x^{\prime}+x_{\mathrm{c}}\delta
u^{\prime}+\left[6x_{\mathrm{c}}-\left(\epsilon_{H\mathrm{c}}+1\right)u_{\mathrm{c}}\right]\delta
x-x_{\mathrm{c}}\left(1+\epsilon_{H\mathrm{c}}\right)\delta
u+2\epsilon_{H\mathrm{c}}z_{\mathrm{c}}\delta
z}{1+u_{\mathrm{c}}x_{\mathrm{c}}}$ (32)

is the linearised eq. (23). The constraint equation fixes the dynamics onto
the three-dimensional hypersurface in the four-dimensional phase space. The
linearised version of that equation is given by

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(x_{\mathrm{c}}-u_{\mathrm{c}}\right)\delta
x+y_{\mathrm{c}}\delta y+z_{\mathrm{c}}\delta z-x_{\mathrm{c}}\delta u\>.$
(33)

We next compute the eigenvalues of the system of equations (28)-(31) and
determine their stability. The eigenvalues at the fixed point M are

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle\epsilon_{H\mathrm{c}}-3\>,$ (34) $\displaystyle m_{2}$
$\displaystyle=$ $\displaystyle\epsilon_{H\mathrm{c}}\>,$ (35) $\displaystyle
m_{3}$ $\displaystyle=$ $\displaystyle-2\epsilon_{H\mathrm{c}}\>,$ (36)

where $\epsilon_{H\mathrm{c}}=\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$ is the
value of the Hubble flow parameter at M. We can see that this fixed point is
always a saddle point for the range of $w_{\mathrm{m}}$ values that we
consider.

The eigenvalues at the kination fixed points (K±) are

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\left(1-w_{\mathrm{m}}\right)\>,$ (37) $\displaystyle
m_{2}$ $\displaystyle=$ $\displaystyle-6\pm\sqrt{6}\kappa\>,$ (38)
$\displaystyle m_{3}$ $\displaystyle=$ $\displaystyle
3\mp\sqrt{\frac{3}{2}}\lambda\>,$ (39)

where the upper sign corresponds to the point K+. We can see that the
eigenvalue $m_{1}$, which corresponds to the eigenvector
$\boldsymbol{v}=\left(0,0,0,1\right)$111We order eigenvector components as
$\boldsymbol{v}=\left(v_{x},v_{y},v_{u},v_{z}\right)$., is always positive.
Hence, these two fixed points are never stable.

The eigenvalues at the scaling fixed point Sc are

$\displaystyle m_{1}$ $\displaystyle=$ $\displaystyle
2\left(\frac{\kappa}{\lambda}-1\right)\epsilon_{H\mathrm{c}}\>,$ (40)
$\displaystyle m_{\pm}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\left(3-\epsilon_{H\mathrm{c}}\right)\left[1\pm\sqrt{1-\frac{8\epsilon_{H\mathrm{c}}}{3-\epsilon_{H\mathrm{c}}}\cdot\left(1-2\frac{\epsilon_{H\mathrm{c}}}{\lambda^{2}}\right)}\right]\>.$
(41)

As can be seen from Table 1, this fixed point exists ($z_{\mathrm{c}}^{2}\geq
0$) only if $\lambda^{2}>2\epsilon_{H\mathrm{c}}$. Such a condition makes the
real part of $m_{\pm}$ always negative. Therefore the stability of this fixed
point is determined solely by the sign of $m_{1}$. That is, the scaling fixed
point Sc is a saddle point for $\kappa>\lambda$.

The eigenvalues at the de Sitter fixed point dS are

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\>,$ (42)
$\displaystyle m_{\pm}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left[1\pm\sqrt{1+\frac{8\lambda^{2}}{3\left(2+3\lambda^{2}\right)}\left(1-\frac{\kappa}{\lambda}\right)}\right]\>.$
(43)

Notice that the condition for the stability of this fixed point is exactly
opposite from the one required by the scaling fixed point Sc: for
$\kappa>\lambda$ the scaling fixed point is a saddle and the de Sitter one is
the attractor. None of the interesting solutions pass through the fixed points
G or S2 so we don’t analyse their stability here.

To visualise the behaviour of the system we integrate numerically a set of
trajectories and show the phase portraits in Figure 1. In all of those
simulations we start with a negligible GB contribution, $u_{0}=-10^{-25}$ to
be precise. If the initial value of $u_{0}$ is too large, the phase portrait
changes drastically. However, the requirement for the negligible GB
contribution at the initial stages of the evolution is consistent with the
required scenario to explain Dark Energy. All the trajectories in the plots in
Figure 1 start from $y_{0}=10^{-3}$ and move towards the dS attractor at
$y=1$.

The physically interesting trajectories are those that start close to the K±
or M fixed points. The former set corresponds to the kination initial
conditions and the latter ones corresponds to a universe where matter
dominates initially. In both cases most trajectories are first attracted to
the scaling fixed point Sc. But because this fixed point is a saddle point for
$\kappa/\lambda>1$, eventually all the trajectories are repelled and move to
the de Sitter attractor dS.

This represents the desirable sequence of events: initially the universe is
dominated by the kinetic energy of the scalar field, which is the case for
quintessential inflation models, or the matter component, which is often
encountered in the quintessence models. Next, the system moves into the
scaling fixed point, and for a long time the Universe evolves with an
effective equation of state that of the matter component.

In the case of GR models and the exponential potential $V\left(\phi\right)$,
the scaling fixed point Sc is an attractor [33]. That is, all the trajectories
converge onto this point and remain there. This is problematic, because the
universe is not accelerating at Sc in contrast to observations. The GB term,
on the other hand, converts this point into a saddle one and provides an
escape route. The scalar field eventually can come to dominate and cause the
universe to expand in an accelerated fashion in a new de Sitter attractor
point dS.

Figure 1: The phase portrait of the dynamical system given by eqs. (20), (21)
and (27), with the matter component represented by a pressureless fluid,
$w_{\mathrm{m}}=0$. All trajectories start with very small values of $y$ and
$u$. The first two columns show a projection of the 3 dimensional phase space.
The full 3D portrait is displayed in the last column. Red regions in the
second column mark forbidden regions of the phase space. For those values the
constraint in eq. (24) cannot be satisfied. The upper plot corresponds to
$\lambda=1.3\times\sqrt{3}$ and the lower one to $\lambda=5\sqrt{3}$. In both
cases $\kappa/\lambda=3/2$.

In Figure 2 we show the time evolution of density parameters and the effective
equation of state of the “dark fluid”. Density parameters are defined as

$\displaystyle\Omega_{\mathrm{m}}$ $\displaystyle\equiv$
$\displaystyle\frac{\rho_{\mathrm{m}}}{3m_{\mathrm{Pl}}^{2}H^{2}}=z^{2}\>,$
(44) $\displaystyle\Omega_{\phi}$ $\displaystyle\equiv$
$\displaystyle\frac{\rho_{\phi}}{3m_{\mathrm{Pl}}^{2}H^{2}}=x^{2}+y^{2}\>,$
(45) $\displaystyle\Omega_{\mathrm{GB}}$ $\displaystyle\equiv$
$\displaystyle-\frac{8H\dot{\xi}}{m_{\mathrm{Pl}}^{2}}=-2ux\>.$ (46)

In terms of these variables the constraint in eq. (24) (the Friedmann
equation) can be written as

$1=\Omega_{\phi}+\Omega_{\mathrm{m}}+\Omega_{\mathrm{GB}}\>.$ (47)

As $\Omega_{\mathrm{GB}}$ can be positive as well as negative, the value of
each $\Omega$ parameter is not bounded by $\left|\Omega\right|\leq 1$.

In order to apply the observational bounds on the equation of state of Dark
Energy we define a new parameter $w_{\mathrm{f}}$. It can be interpreted as
the equation of state of an effective “dark fluid” that causes the accelerated
expansion. To do that eq. (14) can be written as

$\displaystyle\epsilon_{H}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\frac{\rho_{\phi}\left(1+w_{\mathrm{f}}\right)+\rho_{\mathrm{m}}\left(1+w_{\mathrm{m}}\right)}{\rho_{\phi}+\rho_{\mathrm{m}}}\>,$
(48)

where

$\displaystyle w_{\mathrm{f}}$ $\displaystyle\equiv$
$\displaystyle\frac{P_{\phi}+8H^{2}\left(\ddot{\xi}-\dot{\xi}H\right)}{\rho_{\phi}}$
(49)

and $\rho_{\phi}$ and $P_{\phi}$ are the energy and pressure densities of the
scalar field defined in eqs. (11) and (12). In terms of dimensionless
variables in eq. (17) the last expression can be also written as

$\displaystyle w_{\mathrm{f}}$ $\displaystyle=$
$\displaystyle-1+\frac{\frac{2}{3}\epsilon_{H}\left(1+2ux\right)-z^{2}\left(1+w_{\mathrm{m}}\right)}{x^{2}+y^{2}}\>,$
(50)

where $\epsilon_{H}$ is meant to be substituted with eq. (23).

We next run a large number of numerical simulations of eqs. (20), (21) and
(27) varying $\lambda$, $\kappa$ parameters as well as the initial conditions
$x_{0}$ (but always with $u_{0}=-10^{-25}$ and $w_{\mathrm{m}}=0$) and select
those models which have regions where $\Omega_{\mathrm{m}}=0.3147\pm 0.0074$
and $w_{\mathrm{f}}=-0.957\pm 0.08$ [21] are satisfied for the same value of
$N$. The scale factor is normalised such that $a=1$ ($N=0$) at that moment.

In Figure 2 we show two of such models. On the L.H.S. column we can see the
phase portraits, where these models (red curves) are drawn from and on the
R.H.S. column we find the time evolution of the density parameters and
$w_{\mathrm{f}}$. Both models have an initial period of kination, which
quickly gives way to the matter domination . The model in the upper plot, has
a long period of the scaling behaviour before GB energy takes over.
Eventually, all the models settle down at the de Sitter attractor point, where
$\Omega_{\phi}=1$ and all other $\Omega$’s vanish.

In the lower panel of Figure 2 we can also notice a quite generic feature of
GB Dark Energy models, namely, that for a brief period of time the effective
equation of state of the dark fluid can drop below $-1$.

Figure 2: The time evolution of the density and the effective equation of
state for two selected trajectories (denoted by red curves on the left
column). They are selected such that $\Omega_{\mathrm{m}}$ and
$w_{\mathrm{f}}$ are consistent with present day observations.

###  IV.2 The Speed of Gravitational Waves

Obviously, to select a realistic model of DE the consistency with
observational constraints on $\Omega_{\mathrm{m}}$ and $w_{\mathrm{f}}$
parameters is a necessary but not sufficient condition. There are some other
requirements that a viable model of cosmology must satisfy. Among those
requirements, especially in the case of GB model, is a negligibly small
deviation of the speed of gravitational waves $c_{\mathrm{GW}}$ from the speed
of light.

As it was pointed out in [5], generically in scalar-GB models
$c_{\mathrm{GW}}\neq 1$. Due to the tight observational constraints on
$c_{\mathrm{GW}}$ (see eq. (1)) it was deemed that GBDE models are excluded.
However, such constraints do not fix $c_{\mathrm{GW}}=1$, they only place
upper bounds on the deviation from $1$, albeit very strong ones. And, as one
could naively conclude from eq. (8), that bound is not very constraining for
GB models of DE. Moreover, in GB models, $c_{\mathrm{GW}}$ is not a constant,
but varies with time. But the constraint in eq. (1) applies only for a short
period over the history of the Universe.

To investigate this issue let us first write eq. (5) in terms of dimensionless
variables, defined in eq. (17)

$\displaystyle\alpha_{T}$ $\displaystyle=$
$\displaystyle\frac{\left(ux\right)^{\prime}+\left(\epsilon_{H}-1\right)ux}{ux+1/2}=\frac{\Omega_{\mathrm{GB}}^{\prime}+\left(\epsilon_{H}-1\right)\Omega_{\mathrm{GB}}}{\Omega_{\mathrm{GB}}-1}\>,$
(51)

where $\Omega_{\mathrm{GB}}$ is defined in eq. (46). It might appear that
$\alpha_{T}$ diverges at $ux=-1/2$, or equivalently at
$\Omega_{\mathrm{GB}}=-1$. But, as can be seen from the constraint equation
(24), such a value is not allowed.

We can immediately notice from eq. (51) that $\alpha_{T}$ vanishes at the two
fixed points, Sc and dS. These fixed points are the most interesting ones.
Unfortunately, if the GB model is to be a good model of our Universe, the
current stage of the evolution cannot be represented by any of these two fixed
points. Instead, we should find ourselves somewhere on the trajectory between
Sc and dS, as it is also demonstrated in Figure 2.

To see how $\alpha_{T}$ evolves with time we selected a large number of
numerical solutions that satisfy the $\Omega_{\mathrm{m}}$ and
$w_{\mathrm{f}}$ bounds discussed above. Some of those solutions are shown in
Figure 3. For clarity we can draw only a few of them. However, the ones
displayed in Figure 3 are representative of the whole set. First thing we
notice is that the maximum value of $\alpha_{T}$ is always very close to
$N=0$, exactly where the observational bounds in eq. (1) apply. Moreover that
value is always $\left|\alpha_{T}\right|\sim\mathcal{O}\left(1\right)$, which
clearly falls out of the allowed range. Hence, with high confidence we can
conclude that GBDE models with an exponential coupling constant
$\xi\left(\phi\right)$ are excluded by the observational constraint on the
speed of gravitational waves.

Figure 3: The deviation of the speed of gravitational waves from the speed of
light (see eq. (2)) as a function of time. $N=0$ represents the moment where
the model predicts $\Omega_{\mathrm{m}}$ and $w_{\mathrm{f}}$ values
consistent with present day observations. This is also the moment where the
bound in eq. (1) applies. Obviously all models violate this bound by many
orders of magnitude. The approximate excluded region is marked by a red
vertical stripe.

To get some insight why $\alpha_{T}$ is so large at this particular moment, we
can write eq. (51) in an alternative form

$\displaystyle\alpha_{T}$ $\displaystyle=$ $\displaystyle
2\left(\epsilon_{H}-3\right)+3\left(1-w_{\mathrm{m}}\right)\frac{\Omega_{\mathrm{m}}}{1-\Omega_{\mathrm{GB}}}+6\frac{y^{2}}{1-\Omega_{\mathrm{GB}}}\>.$
(52)

This equality can be obtained either by using the dynamical equations
(20)–(25) to eliminate the time derivative from eq. (51), or directly from eq.
(14). In any case, it is important to notice that this expression is very
generic: it is valid for any potential $V\left(\phi\right)$ and any GB
coupling $\xi\left(\phi\right)$, not necessarily exponential. At $N=0$
observations require $\epsilon_{H}\sim 0.5$. Hence, the absolute value of the
first term is of order $\sim 1$. Observations also require $w_{\mathrm{m}}=0$
and $\Omega_{\mathrm{m}}\sim 0.3$. We can use the latter in the constraint in
eq. (47) to write $\Omega_{\phi}+\Omega_{\mathrm{GB}}\sim 0.7$. Moreover,
between the scaling fixed point Sc and the de Sitter one dS, the GB density
parameter $\Omega_{\mathrm{GB}}>0$. Therefore, $0<\Omega_{\mathrm{GB}}<0.7$
and the second term in eq. (52) must be of order $\mathcal{O}\left(0.1\right)$
to $\mathcal{O}\left(1\right)$. Finally $y^{2}\leq\Omega_{\phi}<0.7$ and the
last term of order $\mathcal{O}\left(1\right)$ at the maximum. Barring precise
cancellations, this leads to the conclusion that generically
$\left|\alpha_{T}\right|\sim\mathcal{O}\left(0.1\right)-\mathcal{O}\left(1\right)$.

222The observation that $\Omega_{\mathrm{m}}\sim 0.1$ quite generically leads
to $\left|\alpha_{T}\right|\gtrsim 0.1$ raises another question. As it follows
from eq. (8), such a large value of $\left|\alpha_{T}\right|$ implies a very
large change in the GB coupling. But one has to remember that the GB term is
only the lowest order term in the series of low-energy effective string theory
corrections [9, 12, 13]. Since $\xi$ varies so much, one must wonder if it is
consistent to neglect the higher order corrections over the whole range of the
evolution.

##  V The Linear Gauss-Bonnet Coupling

As we have seen above, models with an exponential GB coupling could
potentially provide a reasonable history of the Universe and explain DE.
Unfortunately, such models predict the speed of GW in the current Universe
that violates the allowed values by many orders of magnitude. In this section
we investigate another option: the linear GB coupling

$\displaystyle\xi$ $\displaystyle\propto$ $\displaystyle\phi\>.$ (53)

As it is shown in the Appendix, such a dynamical system has the same fixed
points as the $\epsilon_{H\mathrm{c}}u_{\mathrm{c}}=0$ subset in Table 1.

To study the stability of those fixed points, we can linearise eqs. (20)–(25)
and use $\xi_{,\phi\phi}=0$. Equivalently, we can use the results in section
IV by setting $\kappa=0$. This gives the eigenvalues at the Sc point

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle-3\left(1+w_{\mathrm{m}}\right)\>,$ (54) $\displaystyle m_{\pm}$
$\displaystyle=$
$\displaystyle-\frac{1}{2}\left(3-\epsilon_{H\mathrm{c}}\right)\left[1\pm\sqrt{1-\frac{8\epsilon_{H\mathrm{c}}}{3-\epsilon_{H\mathrm{c}}}\cdot\left(1-2\frac{\epsilon_{H\mathrm{c}}}{\lambda^{2}}\right)}\right]\>.$
(55)

As discussed in section IV the real part of $m_{\pm}$ eigenvalues are always
negative. We can also see from the first equation above that $m_{1}$ is
negative too. Hence, this fixed point is always an attractor.

On the other hand, the eigenvalues at the de Sitter fixed point dS are

$\displaystyle m_{1}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\>,$ (56)
$\displaystyle m_{\pm}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left[1\pm\sqrt{1+\frac{8\lambda^{2}}{3\left(2+3\lambda^{2}\right)}}\right]\>.$
(57)

It is clear that $m_{1},\>m_{+}<0$ and $m_{-}>0$. Hence, this fixed point is a
saddle.

What the above result shows is that a realistic cosmological scenario is
impossible with the linear GB coupling. The scaling fixed point is an
attractor and there are no solutions which display a long matter dominated
period followed by an accelerated expansion. This is also demonstrated in
Figure 4.

One might wish to generalise the analysis presented here to higher order
polynomial functions $\xi\left(\phi\right)$. Unfortunately, such a functional
form of $\xi\left(\phi\right)$ is not amenable to the presented methods of
analysis, as the equations become non-self-similar. However, we expect the
results to be qualitatively similar to the ones presented in this section.

Figure 4: The phase space of the linear model, with $\xi\left(\phi\right)$
defined in eq. (53). The notation is the same as in Figure 1. All three rows
correspond to the same model with $\lambda=5\sqrt{3}$ and $w_{\mathrm{m}}=0$,
but different initial conditions of $u$. In the first row $u_{0}=-10^{-25}$.
As we can see the evolution remains within the $u=0$ plane. To demonstrate
that the scaling fixed point Sc is an absolute attractor we also run
simulations with large initial values $\left|u_{0}\right|$. In the second row
we can see that some trajectories are initially attracted towards the de
Sitter fixed point dS, but since for a linear $\xi\left(\phi\right)$ this
point is a saddle, eventually all trajectories converge onto Sc. Notice that
in the last column we plot $-u$ on the vertical axis of the first two plots
and $u$ of the third plot.

##  VI The Case of $\alpha_{T}=0$

As can be seen from the expression of $\alpha_{T}$ in eq. (5) this parameter
can be made to vanish if one arranges for the GB coupling in such a way that
its time evolution obeys the condition in eq. (6). Unfortunately, this
condition does not provide us with a functional form of
$\xi\left(\phi\right)$. Nevertheless, it provides with enough information to
investigate the relevant aspects of such a dynamical system.

First, notice that eq. (6), written in terms the dimensionless variables in
eq. (17), becomes

$\displaystyle 24H^{2}\xi_{,\phi\phi}x^{2}$ $\displaystyle=$ $\displaystyle
u\left[\left(\epsilon_{H}+1\right)x-x^{\prime}\right]\,.$ (58)

It provides an additional constraint that can be used to eliminate the
$\xi_{,\phi\phi}$ term in eq. (72), at least for fixed points with
$x_{\mathrm{c}}\neq 0$. Otherwise the LHS of the above equation vanishes in
any case. The fixed points are computed in the Appendix and summarised in
Table 2. As one would expect, all the fixed points with
$u_{\mathrm{c}}\epsilon_{H\mathrm{c}}=0$ are the same as in Table 1.

To determine the stability of fixed points we again linearise eqs. (20)–(25)
and also eq. (58). At the de Sitter fixed point dS this linear system reduces
to

$\displaystyle\delta x^{\prime}$ $\displaystyle=$ $\displaystyle\delta x\,,$
(59) $\displaystyle\delta u^{\prime}$ $\displaystyle=$
$\displaystyle\left(4+3\lambda^{2}\right)\delta x\,,$ (60)
$\displaystyle\delta z^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\delta z\,,$ (61)

where we used that the linearised eq. (58) at dS implies the constraint
$\delta x^{\prime}=\delta x$. Combined with the linearised eqs. (20) and (24)
we find $\delta u=\left(4+3\lambda^{2}\right)\delta x$. It is easy to compute
the eigenvalues, which are

$\displaystyle m_{1}$ $\displaystyle=$ $\displaystyle 0\>,$ (62)
$\displaystyle m_{2}$ $\displaystyle=$ $\displaystyle 1\>,$ (63)
$\displaystyle m_{3}$ $\displaystyle=$
$\displaystyle-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\>.$ (64)

We can see that the last two eigenvalues have opposite signs. Hence, in this
case dS is a saddle fixed point, no longer an attractor.

In the neighbourhood of the scaling fixed point Sc, the linear equation for
$\delta u$ takes the form

$\displaystyle\delta u^{\prime}$ $\displaystyle=$
$\displaystyle\left(1-\epsilon_{H\mathrm{c}}\right)\delta u\>,$ (65)

where $1-\epsilon_{H\mathrm{c}}=-\frac{1}{2}\left(1+3w_{\mathrm{m}}\right)$ is
always negative within the range $0\leq w_{\mathrm{m}}<1$. Therefore, in the
$u$ direction of the phase space this fixed point is attractive and there are
no trajectories that flow from Sc to dS.

In summary, the condition $c_{\mathrm{GW}}=1$ implies the scaling fixed point
Sc to be an attractor, just as in models of GR [33] and dS becomes a saddle
point. This leads to the conclusion that there are no solutions which
reproduce a long, matter-like domination period and asymptotically approach
the de Sitter solution, which is required to reproduce the evolution of the
Universe.

| $x_{\mathrm{c}}$ | $y_{\mathrm{c}}$ | $u_{\mathrm{c}}$ | $z_{\mathrm{c}}$ | Expansion rate $\epsilon_{H\mathrm{c}}$
---|---|---|---|---|---
M | 0 | 0 | 0 | 1 | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$
K± | $\pm 1$ | 0 | 0 | 0 | $3$
I | $\frac{\lambda}{\sqrt{6}}$ | $\sqrt{1-\frac{\lambda^{2}}{6}}$ | 0 | 0 | $\frac{1}{2}\lambda^{2}$
Sc | $\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\lambda}$ | $\frac{\sqrt{\frac{3}{2}\left(1-w_{\mathrm{m}}^{2}\right)}}{\lambda}$ | 0 | $\sqrt{1-\frac{3\left(1+w_{\mathrm{m}}\right)}{\lambda^{2}}}$ | $\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$
dS | 0 | $1$ | $-\sqrt{\frac{3}{2}}\lambda$ | 0 | 0
$\overline{\text{{IV}}}$ | $\sqrt{\frac{2}{3}}\frac{1}{\lambda}$ | $\frac{2}{\sqrt{3}\lambda}$ | $\sqrt{\frac{3}{2}}\frac{1}{\lambda}\left(1-\frac{\lambda^{2}}{2}\right)$ | 0 | 1
Table 2: Fixed points of the dynamical system that satisfies $\alpha_{T}=0$.

##  VII Summary and Conclusions

In this work we investigate Gauss-Bonnet Dark Energy (GBDE) models.
Generically such models predict the speed of gravitational waves different
from the speed of light. In view of the tight observational constraints on
such deviations, denoted by $\alpha_{\mathrm{T}}$ (see eq. (1)), such models
are considered to be disfavoured. However, the deviation is time dependent and
the bound, although tight, is an upper bound, which is only applicable for
very late Universe. Hence, before excluding GBDE models we need to perform a
more detailed analysis. Moreover, if the bound in eq. (1) is expressed in
terms of the variation of the GB coupling function $\xi\left(\phi\right)$, it
might appear to be a weak bound, as shown in eq. (8).

To see if GBDE models can indeed survive the observational constraints on
$\alpha_{\mathrm{T}}$ we perform the dynamical systems analysis. We assume
that the scalar field has an exponential potential and find that such a
dynamical system quite generically has the scaling and de Sitter fixed points
(among others), denoted by Sc and dS respectively in this work.

At Sc the scalar field adjusts in such a way that it mimics the behaviour of
the background matter component. In particular, the equation of state of the
scalar field is the same as that of the matter component. In the case of the
exponential GB coupling, depending on the magnitude of the exponent, Sc is not
stable: it is a saddle point. This is in contrast to the General Relativistic
quintessence models with an exponential potential [33]. In the latter setup
the scaling fixed point is an attractor, which makes it impossible to use as
an explanation for the accelerated expansion of the Universe. In GBDE Sc can
have an unstable direction which links to dS, the latter being an attractor.
Using numerical solutions we show that if the universe starts with a very
small GB term and it is either kination or matter dominated, initially all
solutions evolve towards the Sc fixed point. They linger in the neighbourhood
of Sc for a long period of time and eventually change its course towards dS.
This behaviour is beneficial to modelling quintessential inflation, because it
allows, without extreme fine-tuning, to bridge the enormous energy density gap
(110 orders of magnitude) between infaltion and dark energy. We find numerous
models that follow this scenario and can predict observationally allowed
values for the matter energy density and DE equation of state.

Moreover, at Sc and dS fixed points the speed of gravitational waves is
exactly the same as that of the speed of light, i.e. $\alpha_{\mathrm{T}}=0$.
Unfortunately, if this model is to represent the evolution of the actual
Universe, we cannot be living either on the scaling or de Sitter fixed points,
but somewhere in between. However, as we show in Figure 3,
$\left|\alpha_{\mathrm{T}}\right|$ changes from $0$ to $\sim 1$ in between
these fixed points. Moreover, as we argue below eq. (52), if
$\Omega_{\mathrm{m}}\sim 0.1$, which corresponds to the current value, quite
generically one expects
$\left|\alpha_{\mathrm{T}}\right|\sim\mathcal{O}\left(0.1\right)-\mathcal{O}\left(1\right)$,
which is ruled out by many orders of magnitude. We conclude that the bound in
eq. (1) makes GBDE with an exponential $\xi\left(\phi\right)$ function
unviable.

In view of the above conclusion, we investigated other choices of
$\xi\left(\phi\right)$. First, we demonstrate that for a linear function
$\xi\left(\phi\right)$ the de Sitter fixed point is a saddle point, and the
Scaling fixed point becomes an attractor. This makes it impossible to find any
viable solution that would be consistent with the evolution of the Universe.
We cannot apply our method to study more generic functions
$\xi\left(\phi\right)$, for example monomials. This is because dynamical
equations lose their self-similar character. However, it is reasonable to
expect that the cases studied above are quite representative, and the
behaviour of models with a more general $\xi\left(\phi\right)$ function is
qualitatively very similar.

The advantage of the above described analysis is that we have an explicit
$\xi\left(\phi\right)$ function. But as we showed, this does not provide a
viable cosmological solution. Another hope to make GBDE models conform to
observational constraints is to impose the condition $\alpha_{\mathrm{T}}=0$,
as can be seen in eq. (5). In this case we lose the benefit of having an
explicit functional form of $\xi\left(\phi\right)$, but we gain an additional
constraint equation (58). As it is demonstrated in Section VI, the stability
of fixed points in this case is very similar to the linear model. That is,
there are no solutions which provide a long period of matter domination
followed by an accelerated expansion.

In summary, we find that a GBDE model with an exponential scalar field
potential and an exponential GB coupling function could provide a realistic
model of DE. However, the recent bounds on the speed of GW rules out this
possibility by many orders of magnitude. If, on the other hand, we look for
models that do satisfy $\alpha_{\mathrm{T}}=0$, then it is impossible to find
a scenario consistent with other cosmological observations.

The negative conclusions reached in this work apply to the metric formulation
of the Gauss-Bonnet model. But we know that some modified gravity models that
violate $\alpha_{\mathrm{T}}$ bound in the metric formulation, become viable
again in the Palatini formalism [36]. One can hope that a similar modification
could save the GBDE model too.333Analogously, if the GB term is coming from
the breaking of the Weyl symmetry, also a Weyl term should be added to the
action (see for instance Refs. [37, 38]).  We intend to study this possibility
in future publications.

###### Acknowledgements.

M.K. is supported by the María Zambrano grant, provided by the Ministry of
Universities from the Next Generation funds of the European Union. This work
is also partially supported by the MICINN (Spain) projects
PID2019-107394GB-I00/AEI/10.13039/501100011033 (AEI/FEDER, UE). K.D. was
supported, in part, by the Lancaster-Manchester-Sheffield Consortium for
Fundamental Physics under STFC grant: ST/T001038/1.A.R. is supported by the
Estonian Research Council grant PRG1055. For the purpose of open access, the
authors have applied a Creative Commons Attribution (CC BY) license to any
Author Accepted Manuscript version arising.

##  Appendix A The Fixed Points

In this section we derive fixed points of the dynamical system in eqs.
(20)–(24) that are summarised in Table 1. Such points, or higher dimension
structures are regions of the phase space where
$x^{\prime}=y^{\prime}=z^{\prime}=0$. We denote the constant values of $x$,
$y$, $z$ at these regions by $x_{\mathrm{c}}$, $y_{\mathrm{c}}$,
$z_{\mathrm{c}}$. _A priori_ we do not impose the condition $u^{\prime}=0$ at
fixed points, but it follows from the equations. In order to see that we can
take the derivative of the constraint equation with respect to $N$. At a fixed
point this gives

$\displaystyle\left.u^{\prime}\right|_{x=x_{\mathrm{c}}}x_{\mathrm{c}}$
$\displaystyle=$ $\displaystyle 0\,.$ (66)

This equation allows for $u^{\prime}\neq 0$ if $x_{\mathrm{c}}=0$. However, if
we plug $x^{\prime}=x_{\mathrm{c}}=0$ into eq. (20) and (23) we find

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\sqrt{\frac{3}{2}}\lambda
y_{\mathrm{c}}^{2}+u\left(1-\epsilon_{H\mathrm{c}}\right)\>,$ (67)

where, in this case, $\epsilon_{H\mathrm{c}}$ is the value of $\epsilon_{H}$
in eq. (23) at $x^{\prime}=y^{\prime}=z^{\prime}=x_{\mathrm{c}}=0$. As it is
clear from the above equation, all quantities on the RHS are constant but $u$.
Hence, $u=u_{\mathrm{c}}$ must be also a constant even for $x_{\mathrm{c}}=0$.
In summary, the algebraic equations to find fixed point values are given by

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(\epsilon_{H\mathrm{c}}-3\right)x_{\mathrm{c}}+\sqrt{\frac{3}{2}}\lambda
y_{\mathrm{c}}^{2}+u_{\mathrm{c}}\left(1-\epsilon_{H\mathrm{c}}\right)\>,$
(68) $\displaystyle 0$ $\displaystyle=$
$\displaystyle\left(\epsilon_{H\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda
x_{\mathrm{c}}\right)y_{\mathrm{c}}\>,$ (69) $\displaystyle 0$
$\displaystyle=$
$\displaystyle\left[\epsilon_{H\mathrm{c}}-\frac{3}{2}\left(1+w_{\mathrm{m}}\right)\right]z_{\mathrm{c}}\>,$
(70)

where

$\displaystyle\epsilon_{H\mathrm{c}}$ $\displaystyle=$
$\displaystyle\left[3x_{\mathrm{c}}^{2}+\frac{3}{2}\left(1+w_{\mathrm{m}}\right)z_{\mathrm{c}}^{2}-u_{\mathrm{c}}x_{\mathrm{c}}\right]\frac{1}{1+u_{\mathrm{c}}x_{\mathrm{c}}}$
(71)

and all the values have to satisfy the constraint in eq. (24). There are four
possible solutions of this system:

A

$\left(x,y,u\right)=\left(x_{\mathrm{c}},0,\frac{3\left(w_{\mathrm{m}}-1\right)}{3w_{\mathrm{m}}+1}x_{\mathrm{c}}\right)$,
which is valid for any value of $x_{\mathrm{c}}$. The expansion rate in this
case is given by
$\epsilon_{H\mathrm{c}}=\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$. We can see
that this is the scaling solution, where the scalar field adjusts to the
equation of state of the matter component.

B

$\left(x,y,u\right)=\left(x_{\mathrm{c}},0,\frac{x_{\mathrm{c}}^{2}-1}{2x_{\mathrm{c}}}\right)$
with
$\epsilon_{H\mathrm{c}}=\frac{5x_{\mathrm{c}}^{2}+1}{x_{\mathrm{c}}^{2}+1}$

C

$\left(x,y,u\right)=\left(x_{\mathrm{c}},\sqrt{1-x_{\mathrm{c}}^{2}+2x_{\mathrm{c}}\frac{3x_{\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda}{1+\sqrt{\frac{3}{2}}\lambda
x_{\mathrm{c}}}},\frac{3x_{\mathrm{c}}-\sqrt{\frac{3}{2}}\lambda}{1+\sqrt{\frac{3}{2}}\lambda
x_{\mathrm{c}}}\right)$ with $\epsilon_{H\mathrm{c}}=\sqrt{\frac{3}{2}}\lambda
x_{\mathrm{c}}$

D

$\left(x,y,u\right)=\left(\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\lambda},\frac{\sqrt{\frac{3}{2}\left(1-w_{\mathrm{m}}^{2}\right)+\frac{\lambda}{\sqrt{6}}\left(1+3w_{\mathrm{m}}\right)u_{\mathrm{c}}}}{\lambda},u_{\mathrm{c}}\right)$
with $\epsilon_{H\mathrm{c}}=\frac{3}{2}\left(1+w_{\mathrm{m}}\right)$. This
is a second scaling solution, which is valid for any (allowed by the
constraint) value of $u_{\mathrm{c}}$.

As we can see, these “fixed points” are actually curves in the three
dimensional phase space. In the case of A to C, the curves are parametrised by
the $x_{\mathrm{c}}$ value, which can vary within the range where $y$, $u$ and
$z$ remain real. In the case of D, the value of $x$ is fixed to
$x_{\mathrm{c}}=\sqrt{\frac{3}{2}}\frac{1+w_{\mathrm{m}}}{\lambda}$, but $u$
is free to vary within the similarly defined range.

We can infer more about the structure of fixed points if we use the definition
$u$ in eq. (17). Taking the derivative of that expression we get eq. (25),
which we rewrite it here for convenience:

$\displaystyle u^{\prime}$ $\displaystyle=$
$\displaystyle-2\epsilon_{H}u+24H^{2}\xi_{,\phi\phi}x\>.$ (72)

Hence, at fixed points we get the relation

$\displaystyle 12H^{2}\xi_{,\phi\phi}x_{\mathrm{c}}$ $\displaystyle=$
$\displaystyle\epsilon_{H\mathrm{c}}u_{\mathrm{c}}\>.$ (73)

If we take an exponential GB function, given in eq. (26), the above equation
becomes

$\displaystyle\sqrt{\frac{3}{2}}\kappa x_{\mathrm{c}}u_{\mathrm{c}}$
$\displaystyle=$ $\displaystyle\epsilon_{H\mathrm{c}}u_{\mathrm{c}}\>.$ (74)

Note, that we did not cancel $u_{\mathrm{c}}$ factors, because
$u_{\mathrm{c}}=0$ is an allowed solution. Plugging various values, that are
consistent with this equation, into the system A–D, we obtain fixed points
that are summarised in Table 1. In that table $x_{\mathrm{c}}=\beta$ at the
fixed point G. This is a solution of the cubic equation, which is given by

$\displaystyle\beta$ $\displaystyle\equiv$
$\displaystyle\frac{1}{3\sqrt{3}\kappa}\left(5\sqrt{2}+\frac{9\kappa^{2}-50}{\alpha^{1/3}}-\alpha^{1/3}\right)\>,$
(75)

where

$\displaystyle\alpha$ $\displaystyle\equiv$
$\displaystyle\sqrt{2}\left(27\kappa^{2}-250\right)+\sqrt{2\left(27\kappa^{2}-250\right)^{2}+\left(9\kappa^{2}-50\right)^{3}}$
(76)

and $\kappa$ is an exponent of $\xi\left(\phi\right)$.

It is interesting to note that all fixed points with
$u_{\mathrm{c}}\epsilon_{H\mathrm{c}}\neq 0$ (S2, G, S3 and IV) exist only if
the GB function $\xi\left(\phi\right)$ is of the form

$\displaystyle\xi\left(\phi\right)$ $\displaystyle=$ $\displaystyle
c_{1}\frac{\phi}{m_{\mathrm{Pl}}}+c_{2}\mathrm{e}^{\kappa\frac{\phi}{m_{\mathrm{Pl}}}}\>,$
(77)

where $c_{2}\neq 0$ and
$\kappa=\sqrt{2/3}\cdot\epsilon_{H\mathrm{c}}/x_{\mathrm{c}}$. Or, more
precisely, it is sufficient that the GB function asymptotically approaches
this solution, as the trajectory in the phase space gets closer to those fixed
points. This result can be obtained by integrating eq. (73) and using the fact
that for $x_{\mathrm{c}}\epsilon_{H\mathrm{c}}\neq 0$ we have

$\displaystyle H$ $\displaystyle=$ $\displaystyle
H_{0}\mathrm{e}^{-\frac{\epsilon_{H\mathrm{c}}}{\sqrt{6}x_{\mathrm{c}}}\frac{\phi}{m_{\mathrm{Pl}}}}\>.$
(78)

For the linear $\xi\left(\phi\right)$, i.e. for $\xi_{,\phi\phi}=0$, it
follows from eq. (73) that only fixed points with
$\epsilon_{H\mathrm{c}}u_{\mathrm{c}}=0$ are present. Their values are
summarised in Table 1.

In the case of $\alpha_{T}=0$ constraint, eq. (58) at a fixed point can be
written as

$\displaystyle
24\left.H^{2}\xi_{,\phi\phi}\right|_{\mathrm{c}}x_{\mathrm{c}}^{2}$
$\displaystyle=$
$\displaystyle\left(\epsilon_{H\mathrm{c}}+1\right)x_{\mathrm{c}}u_{\mathrm{c}}\>.$
(79)

For fixed points with $x_{\mathrm{c}}=0$ this equation vanishes identically.
Looking at the system A–D, we see that only two such fixed points exist:
$\left(x,y,u\right)=\left(0,0,0\right)$ and
$\left(x,y,u\right)=\left(0,1,-\sqrt{\frac{3}{2}}\lambda\right)$. These
correspond to points M and dS in Table 1. For $x_{\mathrm{c}}\neq 0$ we can
equate eq. (73) with (79) and get

$\displaystyle 2\epsilon_{H\mathrm{c}}u_{\mathrm{c}}$ $\displaystyle=$
$\displaystyle\left(\epsilon_{H\mathrm{c}}+1\right)u_{\mathrm{c}}\>.$ (80)

For $u_{\mathrm{c}}=0$ the equation vanishes identically and we obtain the
same fixed points as $u_{\mathrm{c}}=0$ points in Table 1. On the other hand,
there is only one fixed point that satisfies $u_{\mathrm{c}}\neq 0$. That
fixed point must have $\epsilon_{H\mathrm{c}}=1$. Looking at the system A–D,
we can see that only one such point is allowed, which is the relation C with
$\left(x,y,u\right)=\left(\sqrt{\frac{2}{3}}\frac{1}{\lambda},\frac{2}{\sqrt{3}\lambda},\sqrt{\frac{3}{2}}\frac{1}{\lambda}\left(1-\frac{\lambda^{2}}{2}\right)\right)$.
All these fixed points are summarised in Table 2.

## References

  * Abbott  _et al._ [2016a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 061102 (2016a), arXiv:1602.03837 [gr-qc] .
  * Abbott  _et al._ [2016b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 116, 241103 (2016b), arXiv:1606.04855 [gr-qc] .
  * Abbott  _et al._ [2017a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott _et al._ (LIGO Scientific, Virgo, Fermi-GBM, INTEGRAL), Astrophys. J. Lett. 848, L13 (2017b), arXiv:1710.05834 [astro-ph.HE] .
  * Ezquiaga and Zumalacárregui [2017] J. M. Ezquiaga and M. Zumalacárregui, Phys. Rev. Lett. 119, 251304 (2017), arXiv:1710.05901 [astro-ph.CO] .
  * Amendola  _et al._ [2018] L. Amendola, I. Sawicki, M. Kunz, and I. D. Saltas, JCAP 08, 030 (2018), arXiv:1712.08623 [astro-ph.CO] .
  * Callan  _et al._ [1969] C. G. Callan, Jr., S. R. Coleman, J. Wess, and B. Zumino, Phys. Rev. 177, 2247 (1969).
  * Gross and Sloan [1987] D. J. Gross and J. H. Sloan, Nucl. Phys. B 291, 41 (1987).
  * Boulware and Deser [1985] D. G. Boulware and S. Deser, Phys. Rev. Lett. 55, 2656 (1985).
  * Antoniadis, Gava, and Narain [1992] I. Antoniadis, E. Gava, and K. S. Narain, Nucl. Phys. B 383, 93 (1992), arXiv:hep-th/9204030 .
  * Antoniadis, Rizos, and Tamvakis [1994] I. Antoniadis, J. Rizos, and K. Tamvakis, Nucl. Phys. B 415, 497 (1994), arXiv:hep-th/9305025 .
  * Bento and Bertolami [1996] M. C. Bento and O. Bertolami, Phys. Lett. B 368, 198 (1996), arXiv:gr-qc/9503057 .
  * Sami  _et al._ [2005] M. Sami, A. Toporensky, P. V. Tretjakov, and S. Tsujikawa, Phys. Lett. B 619, 193 (2005), arXiv:hep-th/0504154 .
  * Aros  _et al._ [2000] R. Aros, M. Contreras, R. Olea, R. Troncoso, and J. Zanelli, Phys. Rev. Lett. 84, 1647 (2000), arXiv:gr-qc/9909015 .
  * Olea [2005] R. Olea, JHEP 06, 023 (2005), arXiv:hep-th/0504233 .
  * Nojiri, Odintsov, and Sasaki [2005] S. Nojiri, S. D. Odintsov, and M. Sasaki, Phys. Rev. D 71, 123509 (2005), arXiv:hep-th/0504052 .
  * Amendola, Charmousis, and Davis [2006] L. Amendola, C. Charmousis, and S. C. Davis, JCAP 12, 020 (2006), arXiv:hep-th/0506137 .
  * Koivisto and Mota [2007a] T. Koivisto and D. F. Mota, Phys. Rev. D 75, 023518 (2007a), arXiv:hep-th/0609155 .
  * Koivisto and Mota [2007b] T. Koivisto and D. F. Mota, Phys. Lett. B 644, 104 (2007b), arXiv:astro-ph/0606078 .
  * van de Bruck  _et al._ [2017] C. van de Bruck, K. Dimopoulos, C. Longden, and C. Owen,  (2017), arXiv:1707.06839 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Cline, Jeon, and Moore [2004] J. M. Cline, S. Jeon, and G. D. Moore, Phys. Rev. D70, 043543 (2004), arXiv:hep-ph/0311312 .
  * Noh and Hwang [2001] H. Noh and J.-c. Hwang, Phys. Lett. B 515, 231 (2001), arXiv:astro-ph/0107069 .
  * Hwang and Noh [2005] J.-c. Hwang and H. Noh, Phys. Rev. D 71, 063536 (2005), arXiv:gr-qc/0412126 .
  * Hwang and Noh [2002] J.-c. Hwang and H. Noh, Phys. Rev. D 66, 084009 (2002), arXiv:hep-th/0206100 .
  * Odintsov and Oikonomou [2019] S. D. Odintsov and V. K. Oikonomou, Phys. Lett. B 797, 134874 (2019), arXiv:1908.07555 [gr-qc] .
  * Kawai, Sakagami, and Soda [1998] S. Kawai, M.-a. Sakagami, and J. Soda, Phys. Lett. B 437, 284 (1998), arXiv:gr-qc/9802033 .
  * Kawai and Soda [1999] S. Kawai and J. Soda, Phys. Lett. B 460, 41 (1999), arXiv:gr-qc/9903017 .
  * Oikonomou [2022] V. K. Oikonomou, Astropart. Phys. 141, 102718 (2022), arXiv:2204.06304 [gr-qc] .
  * Heydari-Fard, Razmi, and Yousefi [2016] M. Heydari-Fard, H. Razmi, and M. Yousefi, Int. J. Mod. Phys. D 26, 1750008 (2016).
  * Ratra and Peebles [1988] B. Ratra and P. J. E. Peebles, Phys. Rev. D37, 3406 (1988).
  * Wetterich [1988] C. Wetterich, Nucl. Phys. B 302, 668 (1988), arXiv:1711.03844 [hep-th] .
  * Copeland, Liddle, and Wands [1998] E. J. Copeland, A. R. Liddle, and D. Wands, Phys.Rev. D57, 4686 (1998), arXiv:gr-qc/9711068 [gr-qc] .
  * Lucchin and Matarrese [1985] F. Lucchin and S. Matarrese, Phys.Rev. D32, 1316 (1985).
  * Tsujikawa and Sami [2007] S. Tsujikawa and M. Sami, JCAP 01, 006 (2007), arXiv:hep-th/0608178 .
  * Kubota  _et al._ [2021] M. Kubota, K.-Y. Oda, K. Shimada, and M. Yamaguchi, JCAP 03, 006 (2021), arXiv:2010.07867 [hep-th] .
  * Corianò, Maglio, and Theofilopoulos [2022] C. Corianò, M. M. Maglio, and D. Theofilopoulos, Eur. Phys. J. C 82, 1121 (2022), arXiv:2203.04213 [hep-th] .
  * Fernandes  _et al._ [2022] P. G. S. Fernandes, P. Carrilho, T. Clifton, and D. J. Mulryne, Class. Quant. Grav. 39, 063001 (2022), arXiv:2202.13908 [gr-qc] .