# Quintessence in the Weyl-Gauss-Bonnet Model

José Jaime Terente Díaz  Konstantinos Dimopoulos  Mindaugas Karčiauskas
Antonio Racioppi

###### Abstract

Quintessence models have been widely examined in the context of scalar-Gauss-
Bonnet gravity, a subclass of Horndeski’s theory, and were proposed as viable
candidates for Dark Energy. However, the relatively recent observational
constraints on the speed of gravitational waves $c_{\textrm{GW}}$ have
resulted in many of those models being ruled out because they predict
$c_{\textrm{GW}}\neq c$ generally. While these were formulated in the metric
formalism of gravity, it was found later that some Horndeski models could be
rescued in the Palatini formalism, where the connection is independent of the
metric and the underlying geometry no longer corresponds to the pseudo-
Riemannian one. Motivated by this and the relation between scalar-Gauss-Bonnet
gravity and Horndeski’s theory, we put forward a new quintessence model with
the scalar-Gauss-Bonnet action but in Weyl geometry. We find the fixed points
of the dynamical system under some assumptions and determine their stability
via linear analysis. Although the past evolution of the Universe as we know it
is correctly reproduced, the constraints on $c_{\textrm{GW}}$ are shown to be
grossly violated for the coupling function under consideration. The case of
$c_{\textrm{GW}}=c$ is regarded also, but no evolution consistent with other
cosmological observations is obtained.

##  1 Introduction

Late Dark Energy (DE) models have become quite popular following the
observational discovery of the current accelerated expansion of the Universe
from Type Ia supernovae [1, 2]. _Planck_ Collaboration in particular finds
that the equation of state (EoS) parameter of this DE is consistent with a
cosmological constant (CC) [3], arguably the simplest candidate for DE. This
CC is interpreted as vacuum energy in particle physics [4, 5, 6], but no
convincing scenario has been found where the actual observed energy scale of
DE can be naturally accounted for by vacuum energy ascribed to a particle
physics model [7, 8].

_Planck_ measurements also leave the door open to a dynamical scalar field
(named ‘quintessence’ [9, 10]), with changing DE EoS parameter. This latter
possibility allows for a more sensible mechanism to explain DE [11, 12, 13]
and paves the way for a unified description (called ‘quintessential inflation’
[14, 15, 16, 17]) of the late accelerated expansion era and early inflation.
The inflationary phase in the early evolution of the Universe serves to
explain the primordial origin of the Large Scale Structures in the Universe
while resolving the horizon and flatness problems of the Hot Big Bang model
[18, 19, 20].

Many of those quintessence DE models rely on modifications of gravity at large
and fall within the scalar-tensor theories chiefly [21, 25, 26, 24, 22, 23].
Consequently, they have been (and still are) under close examination after the
detection of gravitational waves (GWs) made by the Advanced LIGO and Virgo
detectors [27, 28, 29] (see Refs. [30, 31, 32, 33, 34, 35]). The multi-
messenger observation of a binary neutron star system’s merger [36] placed
stringent bounds on the relative difference between the speed of propagation
of GWs, $c_{\textrm{GW}}$, and the speed of light in vacuum, $c$ [37].
Defining the parameter

$\alpha_{T}\equiv c_{\textrm{GW}}^{2}-1~{},$ (1.1)

in units where $c=1$, those bounds yield

$|\alpha_{T}|<10^{-15},$ (1.2)

given the delay between arrival times of the GWs and $\gamma$-rays emitted by
the system’s merger.

Horndeski’s theory [38] is an example of theory of gravity that has been
strongly constrained by Eq. (1.2).111Regarding this bound, its consistency
across the whole spectrum of GW frequencies has been discussed within the
context of Horndeski’s theory in Ref. [39]. This theory is the most general
one in four dimensions leading to second-order field equations of the metric
tensor $g_{\mu\nu}(x)$ and a scalar field $\phi(x)$. It can be shown that, in
order for GWs to propagate at the speed of light ($\alpha_{T}=0$), only the
standard, field-dependent conformal coupling to gravity $G(\phi)R$, where $R$
is the Ricci scalar, is allowed. This reduces significantly the number of
subclasses of Horndeski’s theory conducive to modelling quintessence [40, 31].

Among those subclasses affected by the constraint (1.2), the scalar-Gauss-
Bonnet (SGB) model, the scalar-tensor generalisation of Einstein-Gauss-Bonnet
(EGB) gravity [41], is a well-known case [40] in cosmology. This model
consists of a field-dependent coupling $\xi(\phi)$ to the Gauss-Bonnet (GB)
term

$\mathcal{G}\equiv
R^{2}-4R^{\mu\nu}R_{\mu\nu}+R^{\alpha\beta\mu\nu}R_{\alpha\beta\mu\nu}~{},$
(1.3)

where $R_{\mu\nu}$ and $\tensor{R}{{}^{\alpha}_{\beta\mu\nu}}$ are the Ricci
and Riemann tensors, respectively, in addition to the Einstein-Hilbert action
of General Relativity (GR). The upper bound on $|\alpha_{T}|$ is satisfied if
$\xi$ is constant. This is not surprising as the GB term is a topological term
in four dimensions and does not contribute to the field equations of the
metric [41]. Hence, one is left with Einstein’s theory of GR, which predicts
that GWs propagate at the speed of light. Another possibility is to impose the
following condition:

$\ddot{\xi}=H\dot{\xi}~{},$ (1.4)

which restricts the functional form of the coupling $\xi(t)$, such that
$\dot{\xi}(t)\propto a(t)$. $H(t)\equiv\dot{a}(t)/a(t)$ is the Hubble
parameter and $a(t)$ the scale factor. Overdots denote time derivatives and a
homogeneous, isotropic and spatially-flat spacetime, described by the flat
Friedmann-Lemaître-Robertson-Walker (FLRW) metric
$g_{\mu\nu}=\textrm{diag}\left[-1,a^{2}(t),a^{2}(t),a^{2}(t)\right]$, was
assumed. The consequences of the scaling of $\xi$ after Eq. (1.4) have been
studied in the literature in the context of cosmic inflation [42, 44, 43],
despite the fact that the neutron star system’s merger is a low-redshift
event, significantly lower than the estimated redshift value at the onset of
DE domination [45] and at a time much later than inflation. For this latter
reason, we considered Gauss-Bonnet Dark Energy (GBDE) models with
$\alpha_{T}=0$ in Ref. [46]. Unfortunately, we found that one cannot make
$c_{\textrm{GW}}=1$ if the density parameters of matter and DE are such that
$\Omega_{\textrm{M}}\sim\Omega_{\textrm{DE}}$.

Leaving those two possibilities aside, there are works in which the effects of
the constraint on $|\alpha_{T}|$ are examined and illustrated [48, 47], the
first reference including a case of kinetic coupling to curvature besides the
coupling to the GB term. In Ref. [46], those effects were explored in the GBDE
models mentioned above. We came to the realisation that the bound does not
seem very constraining for those models. Contrary to this intuition however,
it was found that the constraint is easily violated at present
time.222Modifications to the SGB model including a coupling between the matter
sector and the quintessence field or a non-minimal coupling of the scalar
field to gravity seem to satisfy the constraint (1.2) for very specific
initial conditions set in radiation domination and $\xi(\phi)\propto\phi^{2}$.
See results reported in Refs. [49, 50].

Additionally, there has been an intensive research on Horndeski models in the
framework of _the Palatini formalism_ [51, 52, 53, 54, 55].333Recent
quintessential inflation models have been proposed in this formalism in $f(R)$
and $f(\phi,R)$ theories [56, 57]. In the case of $f(\phi)R$ gravity, an
example of quintessence scalar field couplings with metric and torsion can be
found in Ref. [58], where the growth of linear matter perturbations is
studied. A striking result of that investigation has been the realisation that
Horndeski models in this formalism, with a conformal coupling depending on the
kinetic term
$X\equiv-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ (that is,
$G(\phi,X)R$), do allow GWs propagating at the speed of light [53]. This is in
contrast to what the metric one prescribes. Similar attempts were made to
address this question for more complicated Horndeski models [54, 55], but the
connection field equations turn out to be _differential_ rather than
algebraic, so the connection $\Gamma^{\alpha}_{\mu\nu}$ becomes a propagating
degree of freedom of the theory. This entails relevant new aspects in the
physical content of the theory [51]. Even the fact of whether the SGB model
remains a subclass of Horndeski theory in the Palatini formalism is an
unexplored issue in the literature to the best of our knowledge.

Motivated by that relation between the SGB model and Horndeski’s theory in the
metric formalism, and the results of the referenced works involving the speed
of GWs, we propose a quintessence model with the SGB action assuming Weyl
geometry, and determine the behaviour of the $\alpha_{T}$ parameter. A
previous analysis of EGB gravity with a Weyl connection in higher dimensions
was carried out in Ref. [59], but the present work considers a field-dependent
coupling function $\xi(\phi)$ instead. On top of this, we find that the very
same homogeneous and tensor perturbation equations derived here apply to
another connection which, as opposed to the Weyl connection, has non-vanishing
torsion and is metric compatible. The role the projective transformations play
to achieve this is discussed in Sec. C of the Appendix.

Once the equations are derived, we perform a dynamical systems analysis
assuming an exponential potential, in line with Ref. [46]. This potential is
very common in quintessence scenarios [60, 61, 57, 62, 11, 12]. For an
exponential coupling function $\xi(\phi)$, the fixed points are calculated,
those being relevant from the cosmological standpoint. We analyse their
stability and compute the $\alpha_{T}$ parameter along trajectories that
reproduce the observationally constrained values of the density parameter of
matter and the effective EoS parameter of the DE fluid at the present time.
While in the SGB model, Eq. (1.4) must be satisfied in order for $\alpha_{T}$
to vanish, here a different equation is obtained and its stability around the
de Sitter and scaling regimes is examined.

Natural units for which $M_{\textrm{Pl}}\equiv(8\pi G)^{-1/2}$,
$M_{\textrm{Pl}}=2.43\times 10^{18}\,$GeV being the reduced Planck mass, have
been assumed in this work. $G$ is Newton’s gravitational constant.

##  2 Theoretical Framework

###  2.1 Quadratic Gravity in the Palatini Formalism

In the metric formalism of gravity, the connection is set to be

$\mathring{\Gamma}^{\alpha}_{\mu\nu}\equiv\frac{1}{2}g^{\alpha\beta}\left(\partial_{\mu}g_{\beta\nu}+\partial_{\nu}g_{\beta\mu}-\partial_{\beta}g_{\mu\nu}\right),$
(2.1)

which is the well-known ‘Levi-Civita’ (LC) connection [63].444Rings
(overcircles) are used to denote geometric quantities defined with respect to
the LC connection and its first derivatives. We shall adopt this convention
hereon, unless otherwise stated. It is torsion free and metric compatible.
Consequently

$\displaystyle\mathring{\Gamma}^{\alpha}_{\mu\nu}=\mathring{\Gamma}^{\alpha}_{\nu\mu}~{},$
(2.2) $\displaystyle\mathring{\nabla}_{\alpha}g_{\mu\nu}=0~{},$ (2.3)

respectively, where $\mathring{\nabla}_{\mu}$ is the covariant derivative
induced by the LC connection.

In the Palatini formalism of gravity however, the connection
$\Gamma^{\alpha}_{\mu\nu}(x)$ does not depend on the metric tensor
$g_{\mu\nu}(x)$. The dynamics of $\Gamma^{\alpha}_{\mu\nu}$ is governed by
appropriate field equations. While $\Gamma^{\alpha}_{\mu\nu}$ does not
necessarily fulfill Eqs. (2.2) and (2.3), one can relate it to the LC
connection by

$\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}\equiv\Gamma^{\alpha}_{\mu\nu}-\mathring{\Gamma}^{\alpha}_{\mu\nu}~{}.$
(2.4)

Since $\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}$ is the difference of two
connections, it transforms as a tensor under coordinate transformations. This
tensor is called the ‘distortion tensor’ [64]. By defining the torsion and
non-metricity tensors as

$\displaystyle\tensor{T}{{}_{\mu\nu}^{\alpha}}\equiv\Gamma^{\alpha}_{\left[\mu\nu\right]}=\Gamma^{\alpha}_{\mu\nu}-\Gamma^{\alpha}_{\nu\mu}~{},$
(2.5) $\displaystyle Q_{\alpha\mu\nu}\equiv\nabla_{\alpha}g_{\mu\nu}~{},$
(2.6)

respectively, we find the following equations:

$\displaystyle\tensor{T}{{}_{\mu\nu}^{\alpha}}=\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}-\tensor{\kappa}{{}_{\nu\mu}^{\alpha}},$
(2.7) $\displaystyle
Q_{\alpha\mu\nu}=-\kappa_{\alpha\mu\nu}-\kappa_{\alpha\nu\mu}~{},$ (2.8)

where indices are raised and lowered with respect to the metric tensor
$g_{\mu\nu}$. We see that $\tensor{T}{{}_{\mu\nu}^{\alpha}}$ is antisymmetric
in its first two indices, while $Q_{\alpha\mu\nu}$ is symmetric in the last
two.

In this work, we consider $\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}$ of the form

$\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}=\delta_{\mu}^{\alpha}A_{\nu}+\delta_{\nu}^{\alpha}A_{\mu}-g_{\mu\nu}A^{\alpha}.$
(2.9)

$A_{\mu}$ in the above expression is a vector field that we call ‘Weyl vector’
and $\Gamma^{\alpha}_{\mu\nu}$ with that distortion tensor is known as the
‘Weyl connection’ [59]. Then, the corresponding torsion and non-metricity
tensors, using Eqs. (2.7) and (2.8), are respectively,
$\tensor{T}{{}^{\alpha}_{\mu\nu}}=0$ and

$Q_{\alpha\mu\nu}=-2g_{\mu\nu}A_{\alpha}~{}.$ (2.10)

Weyl connection then has non-vanishing non-metricity tensor. A more general
connection is regarded in App. A.

Given the GB term in the metric formalism $\mathring{\mathcal{G}}$

$\mathring{\mathcal{G}}\equiv\mathring{R}^{2}-4\mathring{R}^{\mu\nu}\mathring{R}_{\mu\nu}+\mathring{R}^{\alpha\beta\mu\nu}\mathring{R}_{\alpha\beta\mu\nu}~{},$
(2.11)

one might regard a similar collection of quadratic curvature scalars as the GB
term in the Palatini formalism. However, given the lack of symmetries in the
latter formalism (see App. A), one must take into account that the most
general combination of quadratic terms is more complicated. In fact, such a
combination can be written as [65]

$\displaystyle\mathcal{G}=\alpha
R^{2}+R^{\mu\nu}\left(\beta_{1}R_{\mu\nu}+\beta_{2}R_{\nu\mu}\right)+\tilde{R}^{\mu\nu}\left(\beta_{3}R_{\mu\nu}+\beta_{4}R_{\nu\mu}+\beta_{5}\tilde{R}_{\mu\nu}+\beta_{6}\tilde{R}_{\nu\mu}\right)+\bar{R}^{\mu\nu}\left(\beta_{7}R_{\mu\nu}+\right.$
$\displaystyle\left.+\beta_{8}\tilde{R}_{\mu\nu}+\beta_{9}\bar{R}_{\mu\nu}\right)+R^{\alpha\beta\mu\nu}\left(\gamma_{1}R_{\alpha\beta\mu\nu}+\gamma_{2}R_{\beta\alpha\mu\nu}+\gamma_{3}R_{\mu\nu\alpha\beta}+\gamma_{4}R_{\alpha\mu\nu\beta}+\gamma_{5}R_{\mu\alpha\nu\beta}+\right.$
$\displaystyle\left.+\gamma_{6}R_{\nu\beta\mu\alpha}+\gamma_{7}R_{\beta\nu\mu\alpha}\right),$
(2.12)

where $\tilde{R}_{\mu\nu}$ and $\bar{R}_{\mu\nu}$ are the co-Ricci tensor and
the homothetic curvature tensor, defined in Eqs. (A.5) and (A.6)
respectively.555We refer the reader to Ref. [59], where ‘co-Ricci’ and
‘homothetic’ are used to name the respective tensors as well. In the above,
the constant coefficients $\alpha$, $\left\\{\beta_{i}\right\\}_{i=1,...,6}$
and $\left\\{\gamma_{i}\right\\}_{i=1,...,7}$ are arbitrary. In order to
recover Eq. (2.11) when $\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}=0$ (i.e.
$\Gamma^{\alpha}_{\mu\nu}=\mathring{\Gamma}^{\alpha}_{\mu\nu}$), we set
$(\alpha,\beta,\gamma)=(1,-4,1)$, where

$\displaystyle\beta\equiv\sum^{6}_{i=1}\beta_{i}~{},$ (2.13)
$\displaystyle\gamma\equiv\gamma_{1}-\gamma_{2}+\gamma_{3}~{}.$ (2.14)

The rest of $\gamma$ coefficients are chosen such that
$\gamma_{4}=\gamma_{5}=\gamma_{6}=\gamma_{7}$ [65]. Since
$\mathring{\bar{R}}=0$, there is no need to fix $\beta_{7}$, $\beta_{8}$ and
$\beta_{9}$.

With that parametrisation, $\mathcal{G}$ becomes the Lagrangian of quadratic
terms that we use in the action of the SGB model in the Palatini formalism,
which reads

$S=\int\textrm{d}^{4}x\sqrt{-g}\left[\frac{M^{2}_{\textrm{Pl}}}{2}R-\frac{1}{2}\xi(\phi)\mathcal{G}-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-V(\phi)\right].$
(2.15)

Also, $R=\mathring{R}+\kappa$ (see Eq. (A.9) in App. A). $V(\phi)$ is the
scalar field potential and $\xi(\phi)$ the GB coupling function.

###  2.2 Weyl-Gauss-Bonnet Action

Using the tensors defined in Eqs. (A.2)-(A.7), their symmetry properties that
were shown in App. A, and the scalar quantity $\kappa$ defined in Eq. (A.9),
we can write $\mathcal{G}$ as

$\displaystyle\mathcal{G}=\mathring{\mathcal{G}}+2\mathring{R}\kappa-8\mathring{R}^{\mu\nu}\tilde{\kappa}_{\mu\nu}+2\mathring{R}^{\alpha\beta\mu\nu}\bar{\kappa}_{\alpha\beta\mu\nu}+\kappa^{2}-4\tilde{\kappa}^{\mu\nu}\tilde{\kappa}_{\mu\nu}+\bar{\kappa}^{\alpha\beta\mu\nu}\bar{\kappa}_{\alpha\beta\mu\nu}-\frac{1}{4}\Upsilon\bar{\kappa}^{\mu\nu}\bar{\kappa}_{\mu\nu}~{},$
(2.16)

where $\mathring{\mathcal{G}}$ is defined in Eq. (2.11) and $\Upsilon$ is
given by

$\displaystyle\Upsilon\equiv
1-\beta_{1}+\beta_{2}-2\beta_{7}-4\beta_{9}-2\gamma_{1}+\gamma_{4}~{}.$ (2.17)

As one can see, $\Upsilon$ depends on the arbitrary coefficients, whereas the
rest of the terms in Eq. (2.16) are oblivious to them.

In terms of the Weyl vector and its derivatives, $\mathcal{G}$ reads

$\displaystyle\mathcal{G}=\mathring{\mathcal{G}}+8\left[\left(\mathring{G}^{\mu\nu}-\mathring{\nabla}^{\nu}A^{\mu}\right)\mathring{\nabla}_{\mu}A_{\nu}-\left(\mathring{R}^{\mu\nu}-2\mathring{\nabla}^{\mu}A^{\nu}\right)A_{\mu}A_{\nu}+\left(\mathring{\nabla}_{\sigma}A^{\sigma}+A_{\sigma}A^{\sigma}\right)\mathring{\nabla}_{\rho}A^{\rho}-\right.$
$\displaystyle\left.-\Upsilon\mathring{\nabla}^{\mu}A^{\nu}\left(\mathring{\nabla}_{\mu}A_{\nu}-\mathring{\nabla}_{\nu}A_{\mu}\right)\right].$
(2.18)

Plugging this $\mathcal{G}$ into Eq. (2.15), we obtain the action of what we
call the ‘Weyl-Gauss-Bonnet’ (WGB) model,

$\displaystyle
S_{\textrm{WGB}}=\frac{1}{2}\int\textrm{d}^{4}x\sqrt{-g}\left[M^{2}_{\textrm{Pl}}\mathring{R}-\xi(\phi)\mathring{\mathcal{G}}\right]-4\int\textrm{d}^{4}x\sqrt{-g}\left\\{\frac{3}{4}M^{2}_{\textrm{Pl}}A_{\sigma}A^{\sigma}-\left(\mathring{G}^{\mu\nu}-\mathring{\nabla}^{\mu}A^{\nu}\right)A_{\mu}\partial_{\nu}\xi-\right.$
$\displaystyle\left.-\left(\mathring{\nabla}_{\sigma}A^{\sigma}+A_{\sigma}A^{\sigma}\right)\partial_{\rho}\xi
A^{\rho}-\xi(\phi)\Upsilon\mathring{\nabla}^{\mu}A^{\nu}\left(\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}\right)\right\\}-\int\textrm{d}^{4}x\sqrt{-g}\left[\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+V(\phi)\right]+$
$\displaystyle+\int\textrm{d}^{4}x\sqrt{-g}\mathcal{L}_{\textrm{M}}(g_{\mu\nu},\Psi)~{},$
(2.19)

where we have integrated by parts and used some well-known properties of the
curvature tensors in the metric formalism [63]. We drop the rings from now on
as it is clear that every geometric quantity is defined with respect to the LC
connection and its derivatives once the variables depending on the distortion
tensor have been rewritten in terms of the Weyl vector and its covariant
derivative. Notice also that we have included the matter action
$S_{\textrm{M}}$,

$S_{\textrm{M}}[g_{\mu\nu},\Psi]=\int\textrm{d}^{4}x\sqrt{-g}\mathcal{L}_{\textrm{M}}(g_{\mu\nu},\Psi)~{},$
(2.20)

given our interest in the late Universe. ‘$\Psi$’ is used to denote
collectively the matter fields and $\mathcal{L}_{\textrm{M}}$ is the
corresponding Lagrangian. As we see, the matter fields are minimally coupled
to gravity and the action $S_{\textrm{M}}$ is _not_ a functional of the
connection, only of the metric tensor. This choice is made for simplicity,
although one might explore the possibility of matter fields coupled to the
connection. In that case, besides the energy-momentum tensor, one obtains the
so-called ‘hypermomentum tensor’ from the variation of the action with respect
to the connection [66].

###  2.3 Homogeneous Field Equations

We now present the homogeneous equations of the WGB model. The field equation
of the vector field $A_{\mu}$, derived from Eq. (2.19), is (remember that
rings are dropped for simplicity)

$\displaystyle\frac{3}{2}\left(M^{2}_{\textrm{Pl}}-\frac{4}{3}\partial_{\nu}\xi
A^{\nu}\right)A_{\mu}=\left(\tensor{G}{{}_{\mu}^{\nu}}-2\nabla_{\mu}A^{\nu}\right)\partial_{\nu}\xi+\left(2\nabla_{\sigma}A^{\sigma}+A_{\sigma}A^{\sigma}\right)\partial_{\mu}\xi+\frac{1}{2}\Upsilon\nabla_{\nu}\left(\xi\tensor{\bar{\kappa}}{{}_{\mu}^{\nu}}\right).$
(2.21)

A vector field of the form $A_{\mu}(t)\equiv(-\Phi(t),\mathbf{0})$ is
considered. That is, we take the spatial components of the vector field to
vanish at the background level. This is required to keep the expansion of the
Universe isotropic (see App. B). Such a choice implies
$\bar{\kappa}_{\mu\nu}=0$ in the background. Consequently, $\Upsilon$ does not
play any role in the dynamics of $\Phi$. In particular, the homogeneous part
of Eq. (2.21) reads

$\Phi=2\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}\left(H-\Phi\right)^{2}.$ (2.22)

For $\dot{\xi}\neq 0$, one can write this equation as

$\Phi^{2}-2\Phi\left(H+\frac{M^{2}_{\textrm{Pl}}}{4\dot{\xi}}\right)+H^{2}=0~{}.$
(2.23)

It can be noticed from Eq. (2.22) that $\Phi=0$ for $\dot{\xi}=0$ (constant
$\xi$), and the LC connection is recovered at the background level. Another
possibility that leads to $\Phi=0$ is $H=\Phi$ (see Eq. (2.22)), but this
implies $H=0$, which is not a valid solution when we regard a time dependent
scale factor.

Eq. (2.23) is an algebraic equation that can be readily solved for
$\dot{\xi}\neq 0$

$\Phi=H\left[1+\frac{1}{4\mu}\left(1\pm\sqrt{1+8\mu}\right)\right],$ (2.24)

where $\mu$ is defined as

$\mu\equiv\frac{\dot{\xi}H}{M_{\textrm{Pl}}^{2}}~{}.$ (2.25)

Notice that $\mu$ cannot be smaller than $-1/8$, otherwise $\Phi$ becomes
complex. Since we shall consider expanding universes ($H>0$), this bound has
to be taken into account if $\dot{\xi}$ turns negative.

The _Weyl parameter_ $W\equiv H-\Phi$ is introduced, such that

$\frac{W}{H}=-\frac{1}{4\mu}\left(1-\sqrt{1+8\mu}\right),$ (2.26)

where the positive root in Eq. (2.24) has been discarded. The reason for this
is that if one takes $\dot{\xi}=0$ in Eq. (2.22), then $\Phi=0$. But only the
negative root in Eq. (2.24) yields $\Phi\simeq 0$ (i.e. $W\simeq H$) when
$|\mu|\ll 1$. On the other hand, $|\mu|\gg 1$ in Eq. (2.26) leads to

$\frac{W}{H}\simeq\frac{1}{\sqrt{2\mu}}~{},$ (2.27)

and hence $W\ll H$. For $W$ to be real, $\mu>0$ when $|\mu|\gg 1$ as was
pointed out above in relation to $\Phi$.

Later we will make use of the inverse relation to Eq. (2.26). Using Eq.
(2.22), we can write the $H(W,\dot{\xi})$ function as

$H=W+\Phi=W\left(1+2\frac{\dot{\xi}W}{M^{2}_{\textrm{Pl}}}\right).$ (2.28)

The homogeneous field equations of the metric and the scalar field are given
by (see App. B)

$\displaystyle
H^{2}=\frac{1}{3M^{2}_{\textrm{Pl}}}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\rho_{\textrm{M}}+12\dot{\xi}W^{3}\right]+\Phi^{2},$
(2.29)
$\displaystyle\dot{H}=-\frac{1}{2M^{2}_{\textrm{Pl}}}\left[\dot{\phi}^{2}+\rho_{\textrm{M}}+P_{\textrm{M}}+4\dot{\xi}W^{3}-4\left(\dot{\xi}W^{2}\right)^{\bullet}-12\dot{\xi}\Phi
W^{2}\right]-3\Phi^{2},$ (2.30)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}=-12\xi_{,\phi}W^{2}\left(\dot{W}+HW\right).$
(2.31)

$\rho_{\textrm{M}}$ and $P_{\textrm{M}}$ are the homogeneous energy density
and isotropic pressure of the matter fields, respectively. The energy density
satisfies the following continuity equation:

$\dot{\rho}_{\textrm{M}}=-3H\rho_{\textrm{M}}\left(1+w_{\textrm{M}}\right),$
(2.32)

where $w_{\textrm{M}}\equiv P_{\textrm{M}}/\rho_{\textrm{M}}$ is the EoS
parameter of the background fluid under consideration, which we will assume to
lie in the range $0\leq w_{\textrm{M}}<1$. This range includes the values
corresponding to pressureless matter ($w_{\textrm{M}}=w_{m}=0$) and radiation
($w_{\textrm{M}}=w_{r}=1/3$).

If $\Phi=0$, Eqs. (2.29)-(2.31) are the same as those of the SGB model [44]

$\displaystyle
H^{2}=\frac{1}{3M^{2}_{\textrm{Pl}}}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\rho_{\textrm{M}}+12\dot{\xi}H^{3}\right],$
(2.33)
$\displaystyle\dot{H}=-\frac{1}{2M^{2}_{\textrm{Pl}}}\left[\dot{\phi}^{2}+\rho_{\textrm{M}}+P_{\textrm{M}}+4\dot{\xi}H^{3}-4\left(\dot{\xi}H^{2}\right)^{\bullet}\right],$
(2.34)
$\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}=-12\xi_{,\phi}H^{2}\left(\dot{H}+H^{2}\right).$
(2.35)

This is reasonable as we have not made use of the background equation of
$\Phi$, Eq. (2.22), yet. Now, using Eq. (2.22), Eqs. (2.29)-(2.31) take on the
form

$\displaystyle
W^{2}=\frac{1}{3M^{2}_{\textrm{Pl}}}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\rho_{\textrm{M}}\right],$
(2.36) $\displaystyle\dot{W}+\Phi
W=-\frac{1}{2M^{2}_{\textrm{Pl}}}\left(\dot{\phi}^{2}+\rho_{\textrm{M}}+P_{\textrm{M}}\right),$
(2.37)
$\displaystyle\ddot{\phi}+3W\dot{\phi}+V_{,\phi}=-12\xi_{,\phi}W^{2}\left(W^{2}-\rho_{\textrm{M}}\frac{1+w_{\textrm{M}}}{2M^{2}_{\textrm{Pl}}}\right).$
(2.38)

For $\Phi=0$, these equations reduce to

$\displaystyle
H^{2}=\frac{1}{3M^{2}_{\textrm{Pl}}}\left[\frac{1}{2}\dot{\phi}^{2}+V(\phi)+\rho_{\textrm{M}}\right],$
(2.39)
$\displaystyle\dot{H}=-\frac{1}{2M^{2}_{\textrm{Pl}}}\left(\dot{\phi}^{2}+\rho_{\textrm{M}}+P_{\textrm{M}}\right),$
(2.40) $\displaystyle\ddot{\phi}+3H\dot{\phi}+V_{,\phi}=0~{}.$ (2.41)

We recover the equations of GR above because $\Phi=0$ implies $\dot{\xi}=0$
(unless $H=0$) by virtue of Eq. (2.22).

###  2.4 Equations for Tensor Perturbations

We derived the homogeneous equations and solved the only non-trivial equation
(2.22) that governs the evolution of $A_{\mu}$. Now, we perturb about the flat
FLRW metric and consider tensor perturbations

$\textrm{d}s^{2}=-\textrm{d}t^{2}+a^{2}(t)\left[\delta_{ij}+h_{ij}(t,\mathbf{x})\right]\textrm{d}x^{i}\textrm{d}x^{j},$
(2.42)

where $h_{ij}(t,\mathbf{x})$ are the traceless-transverse tensor components,
such that $\tensor{h}{{}^{i}_{i}}=0=\partial_{i}\tensor{h}{{}^{i}_{j}}$.
Perturbing the action (2.19) up to quadratic terms, we derive the equation of
tensor perturbations666We assume there is no contribution from the matter
action in the form of an anisotropic stress tensor.

$\left(M^{2}_{\textrm{Pl}}-4\dot{\xi}W\right)\ddot{h_{ij}}+\left[3H\left(M^{2}_{\textrm{Pl}}-4\dot{\xi}W\right)-4\left(\ddot{\xi}W+\dot{\xi}\dot{W}\right)\right]\dot{h_{ij}}-\left[M^{2}_{\textrm{Pl}}-4\left(\ddot{\xi}+\dot{\xi}\Phi\right)\right]a^{-2}\partial_{k}\partial^{k}h_{ij}=0~{},$
(2.43)

where $\partial_{k}\partial^{k}\equiv\delta^{kl}\partial_{k}\partial_{l}$.
_Assuming that_ $4\dot{\xi}W\neq M^{2}_{\textrm{Pl}}$, the equation may be
written as

$\ddot{h_{ij}}+\left(3+\alpha_{M}\right)H\dot{h_{ij}}-\left(1+\alpha_{T}\right)a^{-2}\partial_{k}\partial^{k}h_{ij}=0~{},$
(2.44)

where $\alpha_{M}$ and $\alpha_{T}$ are given by

$\displaystyle\alpha_{M}\equiv-\frac{4}{H}\frac{\ddot{\xi}W+\dot{\xi}\dot{W}}{M^{2}_{\textrm{Pl}}-4\dot{\xi}W}~{},$
(2.45)
$\displaystyle\alpha_{T}\equiv-4\frac{\ddot{\xi}-\left(H-2\Phi\right)\dot{\xi}}{M^{2}_{\textrm{Pl}}-4\dot{\xi}W}~{}.$
(2.46)

The speed of GWs predicted by the WGB model is then given by
$c_{\textrm{GW}}^{2}=1+\alpha_{T}$, with $\alpha_{T}$ of Eq. (2.46). Using Eq.
(2.22), we may write the $\alpha_{T}$ parameter entirely in terms of $W$,
$\dot{\xi}$ and $\ddot{\xi}$

$\alpha_{T}=-\frac{4}{M^{2}_{\textrm{Pl}}}\frac{M^{2}_{\textrm{Pl}}\left(\ddot{\xi}-\dot{\xi}W\right)+2\dot{\xi}^{2}W^{2}}{M^{2}_{\textrm{Pl}}-4\dot{\xi}W}~{}.$
(2.47)

##  3 Dynamical Systems Analysis

###  3.1 Dynamical System

Given the homogeneous field equations (2.36)-(2.38), we analyse the generic
behaviour of the dynamical system. We find the corresponding fixed points and
determine their stability keeping in mind the known history of the Universe.
To that end, the following dimensionless variables are defined:

$\displaystyle x\equiv\frac{\dot{\phi}}{\sqrt{6}WM_{\textrm{Pl}}}~{},~{}\
~{}y\equiv\frac{1}{WM_{\textrm{Pl}}}\sqrt{\frac{V}{3}}~{},~{}\
~{}z\equiv\frac{1}{WM_{\textrm{Pl}}}\sqrt{\frac{\rho_{\mathrm{M}}}{3}}~{},$
(3.1)

such that Eq. (2.36), which is the analogue of the Friedmann equation in GR,
can be written as

$1=x^{2}+y^{2}+z^{2}.$ (3.2)

It is noticed that, in contrast to the definitions of $x$, $y$ and $z$ in Ref.
[46], these are given in terms of the Weyl parameter, $W$, not the Hubble one.
$W$ and $H$ are related by Eq. (2.28). Using the dimensionless variables, this
relation can be written as

$H=W\left(1+ux\right),$ (3.3)

where another dimensionless variable,

$u\equiv 2\sqrt{6}\frac{\xi_{,\phi}W^{2}}{M_{\textrm{Pl}}}~{},$ (3.4)

was defined.

A second difference between those definitions and the ones discussed in Ref.
[46] is that $u$ is not constrained by Eq. (3.2). The reason for this
independence is related to the fact that the modifications of gravity are
included in the definitions of $x$, $y$ and $z$ in Eq. (3.1) via the Weyl
parameter. In Ref. [46], on the other hand, the effects of modified gravity
are accounted for solely by the $u$ variable (see Eq. (17) in that paper).

With the definitions of the dimensionless variables above, we can write the
density parameter of the background fluid $\Omega_{\textrm{M}}$ as

$\Omega_{\textrm{M}}\equiv\frac{\rho_{\textrm{M}}}{3H^{2}M^{2}_{\textrm{Pl}}}=\frac{1-x^{2}-y^{2}}{(1+ux)^{2}}~{},$
(3.5)

where Eqs. (3.2) and (3.3) were used.

Using those very same equations, we define the density parameter of the scalar
field $\phi$ too and write it in terms of the dimensionless variables (3.1)
and (3.4)

$\Omega_{\phi}\equiv\frac{\frac{1}{2}\dot{\phi}^{2}+V}{3H^{2}M_{\textrm{Pl}}^{2}}=\frac{x^{2}+y^{2}}{(1+ux)^{2}}~{}.$
(3.6)

Finally, the ‘effective’ density parameter associated with $\xi$ is

$\Omega_{\xi}\equiv\frac{(2+ux)ux}{(1+ux)^{2}}~{},$ (3.7)

such that $1=\Omega_{\phi}+\Omega_{\textrm{M}}+\Omega_{\xi}$.

Assuming an exponential potential

$V(\phi)=V_{0}e^{-\lambda\phi/M_{\textrm{Pl}}},$ (3.8)

such that $V_{,\phi}=-\frac{\lambda}{M_{\textrm{Pl}}}V$ and $\lambda>0$ is a
non-zero constant, the dynamical equations become self-similar given the
definitions in Eqs. (3.1) and (3.4). This class of potentials gives rise to
scaling solutions and they are ubiquitous in quintessence models (see e.g.
Ref. [25], where the exponential potential is the large field limit of the one
proposed therein). Then, taking the derivatives of $x$, $y$, $z$ and $u$, and
using Eq. (2.31), the scalar field equation, the continuity equation (2.32),
and Eq. (3.3), we arrive at

$\displaystyle
x^{\prime}=-\frac{1}{1+ux}\left[(3x+u)(1+ux)-\sqrt{\frac{3}{2}}\lambda
y^{2}+(x+u)\frac{\dot{W}}{W^{2}}\right],$ (3.9) $\displaystyle
y^{\prime}=-\frac{y}{1+ux}\left(\frac{\dot{W}}{W^{2}}+\sqrt{\frac{3}{2}}\lambda
x\right),$ (3.10) $\displaystyle
z^{\prime}=-z\left[\frac{\dot{W}}{W^{2}}\frac{1}{1+ux}+\frac{3}{2}\left(1+w_{\mathrm{M}}\right)\right],$
(3.11) $\displaystyle
u^{\prime}=\frac{2}{1+ux}\left(u\frac{\dot{W}}{W^{2}}+6\xi_{,\phi\phi}W^{2}x\right).$
(3.12)

$\dot{W}/W^{2}$ can be expressed in terms of the dimensionless variables using
Eq. (2.37). Plugging Eq. (2.22) into Eq. (2.37) we find

$-\frac{\dot{W}}{W^{2}}=x\left(3x+u\right)+\frac{3}{2}z^{2}\left(1+w_{\textrm{M}}\right).$
(3.13)

Primes denote derivatives with respect to the number of elapsing e-folds,
which we choose to be $\textrm{d}N\equiv\textrm{d}\ln a=H\textrm{d}t$.

$\dot{W}/W^{2}$ can be related to the Hubble flow parameter $\epsilon_{H}$,
which is defined as [67]

$\epsilon_{H}\equiv-\frac{\dot{H}}{H^{2}}~{}.$ (3.14)

$\epsilon_{H}$ parametrises the acceleration of spatial slices. Using Eq.
(2.28) we find

$\epsilon_{H}=-\frac{1}{1+ux}\left[\frac{\dot{W}}{W^{2}}+\frac{(ux)^{\bullet}}{H}\right].$
(3.15)

Notice that no assumption regarding $\xi(\phi)$ has been made so far. All the
previous equations are valid for any coupling function, including the
exponential one

$\xi(\phi)=\xi_{0}e^{\kappa\phi/M_{\textrm{Pl}}},$ (3.16)

such that $\xi_{,\phi\phi}=\frac{\kappa}{M_{\textrm{Pl}}}\xi_{,\phi}$, where
$\kappa>0$ is a non-zero constant.777The strength of the exponential $\kappa$
in Eq. (3.16) should not be confused with the scalar $\kappa$ defined in Eq.
(A.9), discussed in Sec. 2.2 and App. A. Assuming such a coupling function,
Eq. (3.12) is a self-similar equation too and the explicit dependence on $W$
drops out of the equation888For a linear coupling function
$\xi(\phi)\propto\phi$, the explicit dependence on the Weyl parameter $W$ also
disappears from Eq. (3.12). However, as we will argue later on (see footnote
12), the linear coupling case can be discarded on the basis of our knowledge
of the past history of the Universe.

$u^{\prime}=\frac{2u}{1+ux}\left(\frac{\dot{W}}{W^{2}}+\sqrt{\frac{3}{2}}\kappa
x\right).$ (3.17)

###  3.2 Fixed Points

Fixed points of the dynamical system of equations (3.9), (3.10) and (3.17) are
denoted by $x_{c}$, $y_{c}$ and $u_{c}$, respectively, and correspond to
regions of the phase space where $x^{\prime}=y^{\prime}=u^{\prime}=0$. Among
those, we can list the following ones:

  * •

M: $(x_{c},y_{c},u_{c})=(0,0,0)$ and
$\left.\epsilon_{H}\right|_{c}=\frac{3}{2}\left(1+w_{\textrm{M}}\right)$;

  * •

K$\pm$: $(x_{c},y_{c},u_{c})=(\pm 1,0,0)$ and
$\left.\epsilon_{H}\right|_{c}=3$;

  * •

I:
$(x_{c},y_{c},u_{c})=\left(\frac{\lambda}{\sqrt{6}},\sqrt{1-\frac{\lambda^{2}}{6}},0\right)$
and $\left.\epsilon_{H}\right|_{c}=\frac{\lambda^{2}}{2}$;

  * •

ScI:
$(x_{c},y_{c},u_{c})=\left(\sqrt{\frac{3}{2}\frac{\left(1+w_{\textrm{M}}\right)^{2}}{\lambda^{2}}},\sqrt{\frac{3}{2}\frac{1-w^{2}_{\textrm{M}}}{\lambda^{2}}},0\right)$
and $\left.\epsilon_{H}\right|_{c}=\frac{3}{2}\left(1+w_{\textrm{M}}\right)$;

which are fixed points in GR too because $W_{c}=H_{c}$ at those fixed points
(see Eq. (3.3)). $z_{c}$ is determined from the constraint equation (3.2), so
we do not really need to impose $z^{\prime}=0$; it follows from the constraint
equation.999Notice from Eq. (3.11) that, in order for $z^{\prime}=0$, either
$z_{c}=0$ or
$\left.\epsilon_{H}\right|_{c}=\frac{3}{2}\left(1+w_{\textrm{M}}\right)$ (see
Eq. (3.15)). Both are satisfied (although not simultaneously) by the fixed
points above and by the ones we discuss later. Besides this, we remark that,
unlike $z$, $u^{\prime}=0$ is imposed as $x^{\prime}=y^{\prime}=0$ are
satisfied for $u=u(N)$ and $x_{c}=\pm 1$ and $y_{c}=0$, for example. As was
argued already, this is due to the fact that the modifications of gravity are
included in all the definitions of the dimensionless variables, Eqs. (3.1) and
(3.4), whilst in our previous work, Ref. [46], these were only represented by
$u$, which is defined in Eq. (17) of that paper.

M represents the regime of the background fluid domination. K$\pm$ stands for
the period of kination (dominance of the kinetic energy of the scalar field
$\dot{\phi}^{2}/2$) and I corresponds to a regime of power law inflation (when
$\lambda<\sqrt{2}$, so that $\left.\epsilon_{H}\right|_{\textrm{{I}}}<1$)
where matter fields do not contribute ($z_{\textrm{{I}}}=0$; see Eq. (3.2)).
Notice that this latter fixed point exists if $\lambda<\sqrt{6}$, so
$y_{\textrm{{I}}}$ is real (remember that $\lambda>0$). Lastly, we have the
scaling regime (ScI). In this case, the EoS parameter of the scalar field,

$w_{\phi}\equiv\frac{\frac{1}{2}\dot{\phi}^{2}-V(\phi)}{\frac{1}{2}\dot{\phi}^{2}+V(\phi)}~{},$
(3.18)

becomes equal to that of the background fluid, $w_{\textrm{M}}$. This fixed
point yields

$x^{2}_{\textrm{{ScI}}}+y^{2}_{\textrm{{ScI}}}=\frac{3}{\lambda^{2}}\left(1+w_{\textrm{M}}\right),$
(3.19)

and, by virtue of Eq. (3.2), we obtain the following existence condition:

$\lambda>\sqrt{3(1+w_{\textrm{M}})}\geq\sqrt{3}~{}.$ (3.20)

Therefore, we observe that I and ScI coexist whenever
$\sqrt{6}>\lambda>\sqrt{3}$. If I corresponds to power law inflation, then
both fixed points cannot coexist because $\lambda<\sqrt{2}$ is required, as
was pointed out already.

In addition to the fixed points listed above, we find five other points:

  * •

dS: $(x_{c},y_{c},u_{c})=\left(0,1,\sqrt{\frac{3}{2}}\lambda\right)$ and
$\left.\epsilon_{H}\right|_{c}=0$;

  * •

ScII$\pm$:
$(x_{c},y_{c},u_{c})=\left(\frac{\kappa}{\sqrt{6}}\frac{1+3w_{\textrm{M}}}{3(1-w^{2}_{\textrm{M}})}\Delta_{\pm}(\kappa,w_{\textrm{M}}),0,\frac{\kappa}{\sqrt{6}}\frac{1}{1+w_{\textrm{M}}}\Delta_{\pm}(\kappa,w_{\textrm{M}})\right)$,
where

$\Delta_{\pm}(\kappa,w_{\textrm{M}})\equiv
1\pm\sqrt{1-\frac{18}{\kappa^{2}}\frac{(1+w_{\textrm{M}})(1-w^{2}_{\textrm{M}})}{1+3w_{\textrm{M}}}}~{},$
(3.21)

and $\left.\epsilon_{H}\right|_{c}=\frac{3}{2}\left(1+w_{\textrm{M}}\right)$;

  * •

K̂$\pm$: $(x_{c},y_{c},u_{c})=\left(\pm 1,0,\sqrt{\frac{3}{2}}\kappa\mp
3\right)$ and
$\left.\epsilon_{H}\right|_{c}=\frac{\mp\frac{3\kappa}{2\sqrt{6}}}{1\mp\frac{3\kappa}{2\sqrt{6}}}$.

The first fixed point corresponds to the de Sitter expansion, where the total
energy density consists only of the CC. This fixed point exists in GBDE models
even for $\lambda\neq 0$. It is also found in the SGB model (see Refs. [46,
68]) because $x_{\textrm{{dS}}}=0$ and $W_{\textrm{{dS}}}=H_{\textrm{{dS}}}$
(see Eq. (3.3)).

Besides the de Sitter fixed point, we have a second scaling regime which
consists of two fixed points (ScII$\pm$), depending on the root sign in Eq.
(3.21). These fixed points exist if

$\kappa^{2}\geq
18\frac{(1+w_{\textrm{M}})(1-w^{2}_{\textrm{M}})}{1+3w_{\textrm{M}}}~{},$
(3.22)

so that $x_{\textrm{{ScII}$\pm$}}$ and $u_{\textrm{{ScII}$\pm$}}$ are real.
Then $1\leq\Delta_{+}<2$ and $0<\Delta_{-}\leq 1$. Also, Eq. (3.2) demands

$x^{2}_{\textrm{{ScII}$\pm$}}=\frac{\kappa^{2}}{54}\frac{(1+3w_{\textrm{M}})^{2}}{(1-w^{2}_{\textrm{M}})^{2}}\Delta_{\pm}^{2}<1~{}.$
(3.23)

Combining it with Eq. (3.22) we obtain

$\Delta_{\pm}^{2}<\frac{3(1-w_{\textrm{M}})}{1+3w_{\textrm{M}}}~{}.$ (3.24)

Given that $1\leq\Delta_{+}^{2}<4$ and $0<\Delta_{-}^{2}\leq 1$, we have
$0\leq w_{\textrm{M}}<\frac{1}{3}$ for $\Delta_{+}$ and $0\leq
w_{\textrm{M}}<1$ for $\Delta_{-}$ (we remind the reader that the whole range
of values of $w_{\textrm{M}}$ is assumed to be $0\leq w_{\textrm{M}}<1$; see
below Eq. (2.32)). Therefore, the scaling fixed point ScII$+$ in particular
does not exist if the background fluid under consideration corresponds to
radiation for example, where $w_{\textrm{M}}=w_{r}=1/3$. The scaling regime is
still represented by ScII$-$ in that case though.

Finally, the fourth and fifth fixed points listed above have vanishing
potential energy and the energy density of the background fluid is zero. We
denote them by ‘K̂$\pm$’ and call the corresponding regime ‘pseudo-kination’
since
$3W^{2}_{\textrm{{\^{K}}}\pm}M^{2}_{\textrm{Pl}}=\dot{\phi}_{\textrm{{\^{K}}}\pm}^{2}/2$,
although $W_{\textrm{{\^{K}}}\pm}\neq H_{\textrm{{\^{K}}}\pm}$. Irrespective
of the sign, if $\kappa\gg 1$, then
$\left.\epsilon_{H}\right|_{\textrm{{\^{K}}}\pm}\rightarrow 1$. This value of
the Hubble flow parameter, for large $\kappa$, leads into a universe that does
not accelerate.

For $x_{c}=+1$ (i.e. K̂$+$) and $\kappa<2\sqrt{\frac{2}{3}}$,
$\left.\epsilon_{H}\right|_{\textrm{{\^{K}}}+}<0$ (remember that $\kappa$ is
defined to be positive; see below Eq. (3.16)). This implies that the Universe
accelerates faster than during de Sitter, where it accelerates
exponentially.101010In order to see why it accelerates faster than
exponential, we solve $\epsilon_{H}=-\alpha$, where $\alpha>0$ is constant
($\epsilon_{H}$ is constant at the fixed points K̂$\pm$, indeed)
$H(t)=\frac{H_{0}}{1-\alpha H_{0}t}~{},$ (3.25) where $H_{0}=H(t=0)$. Solving
for $a(t)$ we obtain  $a(t)=\frac{a_{0}}{(1-\alpha H_{0}t)^{1/\alpha}}~{},$
(3.26) such that $a_{0}=a(t=0)$. Assuming an expanding universe, $H_{0}>0$ and
$t<t_{s}\equiv 1/(\alpha H_{0})$. Then, $a(t)$ diverges (approaches
$a(t_{s})$) faster than if $a(t)\propto e^{H_{0}t}$.

Lastly, choosing $\lambda=\kappa$ allows for all $x_{c}$, $y_{c}$ and $u_{c}$
being non-zero at once as can be seen from Eqs. (3.10) and (3.17). The same
occurs in the SGB model (see Ref. [46]). The corresponding fixed ‘curves’ are

  * •

ScIII:
$(x,y,u)=\left(\frac{\frac{3}{2}(1+w_{\textrm{M}})}{\sqrt{\frac{3}{2}}\lambda-\frac{3}{2}u_{c}(1+w_{\textrm{M}})},\sqrt{\frac{\sqrt{\frac{3}{2}}\lambda
u_{c}+\frac{3}{2}(1+w_{\textrm{M}})\left(\frac{3}{2}(1-w_{\textrm{M}})-\sqrt{\frac{3}{2}}\lambda
u_{c}+\frac{1}{2}(1+3w_{\textrm{M}})u_{c}^{2}\right)}{\left(\sqrt{\frac{3}{2}}\lambda-\frac{3}{2}u_{c}(1+w_{\textrm{M}})\right)^{2}}},u_{c}\right)$
and $\epsilon_{H}|_{c}=\frac{3}{2}(1+w_{\textrm{M}})$;

  * •

$\widehat{\textrm{{IV}}}$:
$(x,y,u)=\left(x_{c},\sqrt{1-x_{c}^{2}},-3x_{c}+\sqrt{\frac{3}{2}}\lambda\right)$
and $\epsilon_{H}|_{c}=\frac{\sqrt{\frac{3}{2}}\lambda
x_{c}}{1-3x_{c}^{2}+\sqrt{\frac{3}{2}}\lambda x_{c}}$.

We then have a third scaling solution, ScIII, and a regime of zero energy
density of the background fluid, corresponding to the fixed ‘curve’
$\widehat{\textrm{{IV}}}$. This curve is denoted that way following the
convention in Ref. [46] (see Table I in that paper), although neither of those
two fixed ‘curves’ are found in the SGB model strictly because $W_{c}\neq
H_{c}$, so carets are used on the latter to differentiate it from the fixed
curve in the SGB model (similarly to what was done with the pseudo-kination
fixed points we just discussed).

Those two regimes (especially the scaling one) could be of interest, although
in this work we shall consider $\kappa\neq\lambda$. We leave the stability
analysis of these two fixed curves and possible solutions for future works.

###  3.3 Stability

To perform the stability analysis of the fixed points shown above, we perturb
Eqs. (3.9), (3.10) and (3.17) linearly such that

$\displaystyle\left(1+u_{c}x_{c}\right)\delta
x^{\prime}=-\left[3-9x_{c}^{2}-2u_{c}x_{c}-\frac{3}{2}(1+w_{\textrm{M}})\left(1-3x^{2}_{c}-y_{c}^{2}-2u_{c}x_{c}\right)\right]\delta
x+$
$\displaystyle+y_{c}\left[\sqrt{6}\lambda-3(1+w_{\textrm{M}})(x_{c}+u_{c})\right]\delta
y-\left[1-x_{c}^{2}-\frac{3}{2}\left(1+w_{\textrm{M}}\right)\left(1-x_{c}^{2}-y_{c}^{2}\right)\right]\delta
u~{},$ (3.27) $\displaystyle(1+u_{c}x_{c})\delta
y^{\prime}=-y_{c}\left[\sqrt{\frac{3}{2}}\lambda-
u_{c}-3x_{c}(1-w_{\textrm{M}})\right]\delta
x-\left[\left(\sqrt{\frac{3}{2}}\lambda-3x_{c}-u_{c}\right)x_{c}-\right.$
$\displaystyle\left.-\frac{3}{2}(1+w_{\textrm{M}})\left(1-x_{c}^{2}-3y_{c}^{2}\right)\right]\delta
y+y_{c}x_{c}\delta u~{},$ (3.28) $\displaystyle\left(1+u_{c}x_{c}\right)\delta
u^{\prime}=2u_{c}\left[\sqrt{\frac{3}{2}}\kappa-
u_{c}-3x_{c}\left(1-w_{\textrm{M}}\right)\right]\delta
x+6u_{c}y_{c}\left(1+w_{\textrm{M}}\right)\delta y+$
$\displaystyle+2\left[\left(\sqrt{\frac{3}{2}}\kappa-3x_{c}-2u_{c}\right)x_{c}-\frac{3}{2}\left(1+w_{\textrm{M}}\right)\left(1-x_{c}^{2}-y_{c}^{2}\right)\right]\delta
u~{},$ (3.29)

where the constraint equation (3.2) was used to replace $z^{2}$ in Eq. (3.13).

We begin our stability analysis with the fixed points M and K$\pm$, given the
relevance of the latter in particular in early stages in models of
quintessential inflation. In the case of M, which corresponds to the stage of
the background fluid domination, as was commented already, the eigenvalues,
which we denote by ‘$m$’, are the following ones:

$\displaystyle m_{1}=-\frac{3}{2}\left(1-w_{\textrm{M}}\right),$ (3.30)
$\displaystyle m_{2}=\frac{3}{2}\left(1+w_{\textrm{M}}\right),$ (3.31)
$\displaystyle m_{3}=-3\left(1+w_{\textrm{M}}\right).$ (3.32)

One of the eigenvalues, $m_{2}$, which corresponds to the eigenvector
$\mathbf{v}=(0,1,0)$,111111We order eigenvector components as
$\mathbf{v}=(v_{x},v_{y},v_{u})$. is non-negative, whereas $m_{1},m_{3}<0$. In
view of this, M is a saddle point and is unstable in the $y$ direction.

Regarding K$\pm$, the eigenvalues are given by

$\displaystyle m_{1}=3\left(1-w_{\textrm{M}}\right),$ (3.33) $\displaystyle
m_{2}=3\mp\sqrt{\frac{3}{2}}\lambda~{},$ (3.34) $\displaystyle
m_{3}=-6\pm\sqrt{6}\kappa~{},$ (3.35)

where the upper sign corresponds to the fixed point K$+$. $m_{1}$,
corresponding to the eigenvector $\mathbf{v}=(1,0,0)$, is always positive,
irrespective of $\lambda$ and $\kappa$. Thus, the kination period is unstable
in the $x$ direction. The stability of these fixed points, M and K$\pm$, that
is shown here, resembles the one in the SGB model [68, 46].

We now analyse the stability of the fixed points ScI and dS, which stand out
amongst the rest of the fixed points given our interest in quintessence
cosmology. For ScI, the eigenvalues can be easily calculated to be

$\displaystyle
m_{1}=3(1+w_{\textrm{M}})\left(\frac{\kappa}{\lambda}-1\right),$ (3.36)
$\displaystyle
m_{2,3}=-\frac{3}{4}(1-w_{\textrm{M}})\left\\{1\mp\sqrt{1-8\frac{1+w_{\textrm{M}}}{1-w_{\textrm{M}}}\left[1-\frac{3}{\lambda^{2}}(1+w_{\textrm{M}})\right]}\right\\}.$
(3.37)

Since $\lambda>\sqrt{3(1+w_{\textrm{M}})}$ (see Eq. (3.19)), $m_{2,3}$ are
complex eigenvalues whose real part is always negative. On the other hand, if
$\kappa>\lambda$, $m_{1}$ will be positive, and the fixed point ScI will be a
saddle point. This same feature is shared with the SGB model [46, 68], and it
is an important one because the scaling regime, in which the expansion of the
Universe does not accelerate in contrast to observations [69], is brought to
an end.121212For a linear coupling function $\xi(\phi)\propto\phi$, we have
that $\kappa=0$ because $\xi_{,\phi\phi}=0$. This implies that $m_{1}<0$ in
Eq. (3.36), irrespective of $\lambda$. Then, the scaling fixed point is
stable. Moreover, we see that dS is unstable (see Eq. (3.39), where
$m_{2}>0$). The stability of the scaling fixed point ScI contradicts
observations indicating that the Universe undergoes a phase of current
accelerating expansion, which cannot correspond to the scaling regime with
$\epsilon_{H}=\frac{3}{2}(1+w_{\textrm{M}})$ as commented already. Therefore,
the linear coupling can be discarded. Also, $\kappa=0$ implies that either
$u_{c}=0$ or $\left.\epsilon_{H}\right|_{c}=0$ (see Eqs. (3.17) and (3.15)).
This implies that no other scaling regime with $u_{c}\neq 0$ exists. The same
conclusions can be drawn in the SGB model (see Ref. [46]).

The Universe can enter the de Sitter phase depending on the stability of this
fixed point. The eigenvalues are

$\displaystyle m_{1}=-3\left(1+w_{\textrm{M}}\right),$ (3.38) $\displaystyle
m_{2,3}=-\frac{3}{2}\left[1\mp\sqrt{1-\frac{4}{3}\lambda^{2}\left(\frac{\kappa}{\lambda}-1\right)}\right].$
(3.39)

As long as $\kappa>\lambda$, all the eigenvalues will have negative real
parts, and dS will be a stable fixed point (an attractor). The condition
$\kappa>\lambda$ was sufficient in order for ScI to be unstable (see Eq.
(3.36)), and therefore, the Universe may transition from ScI to dS. This
reproduces the evolution of the Universe as we know it: from radiation/matter
domination to an eventual exponentially expanding universe. The fact that dS
may be a stable fixed point is a well-known result in the SGB model as well
[68, 46].

For completeness, we include the stability analysis of I before getting to the
other two fixed points with $u_{c}\neq 0$. In this case, the eigenvalues read

$\displaystyle m_{1}=\lambda^{2}\left(\frac{\kappa}{\lambda}-1\right),$ (3.40)
$\displaystyle m_{2}=\frac{1}{2}\left(\lambda^{2}-6\right),$ (3.41)
$\displaystyle m_{3}=\lambda^{2}-3\left(1+w_{\textrm{M}}\right).$ (3.42)

Since $\lambda^{2}<6$ in order for I to exist (so that $y_{\textrm{{I}}}$ is
real, as was commented earlier), we have that $m_{2}<0$. Given the range of
values of $w_{\textrm{M}}$, $m_{3}$ may be non-negative ($\lambda^{2}$ would
have to be smaller than $3(1+w_{\textrm{M}})$ otherwise). Considering power
law inflation though, such that $\lambda^{2}<2$, we see that $m_{3}$ is
definitely negative. Regardless of $w_{\textrm{M}}$ and $\lambda$, $m_{1}$
will be positive if $\kappa>\lambda$. Then, I will be unstable if
$\kappa>\lambda$.

The fixed point K̂$+$ could lead to a phase of faster than exponential
acceleration, $\left.\epsilon_{H}\right|_{\textrm{{\^{K}}}+}<0$, if
$\kappa<2\sqrt{\frac{2}{3}}$ (remember that we are assuming that $\kappa$ does
not vanish; see Eq. (3.16)). We study the stability of the two fixed points
K̂$\pm$. The eigenvalues are

$\displaystyle
m_{1}=-\frac{2}{2\mp\sqrt{\frac{3}{2}}\kappa}\left[\pm\sqrt{\frac{3}{2}}\kappa+\frac{3}{2}\left(1+w_{\textrm{M}}\right)\left(2\mp\sqrt{\frac{3}{2}}\kappa\right)\right],$
(3.43) $\displaystyle
m_{2}=\frac{\mp\sqrt{\frac{3}{2}}\lambda}{2\mp\sqrt{\frac{3}{2}}\kappa}\left(\frac{\kappa}{\lambda}-1\right),$
(3.44) $\displaystyle
m_{3}=-\frac{2}{2\mp\sqrt{\frac{3}{2}}\kappa}\left(3\mp\sqrt{\frac{3}{2}}\kappa\right),$
(3.45)

the upper sign corresponding to K̂$+$. In the case of K̂$-$, the fixed point
will be unstable whenever $\kappa>\lambda$ (because $m_{2}>0$), regardless of
the magnitude of $\kappa$ (given that $m_{1}$ and $m_{3}$ are always negative
for K̂$-$). For K̂$+$ and $\kappa<2\sqrt{\frac{2}{3}}$, $m_{1}$ and $m_{3}$
are both negative too. If we assume that $\kappa>\lambda$, then $m_{2}$,
corresponding to the eigenvector $\mathbf{v}=(0,1,0)$, will be negative and
the phase of super-acceleration will be _stable_.131313‘Super-acceleration’ in
relation to $\epsilon_{H}<0$ is used, for example, in Refs. [70, 71]. If not,
it will be unstable in the $y$ direction. Irrespective of the sign of $x_{c}$,
if $\kappa\gg 1$, but $\lambda>\kappa$, the asymptotic phase of zero
acceleration, corresponding to
$\left.\epsilon_{H}\right|_{\textrm{{\^{K}}}\pm}=1$, will be stable as well.

Notice that, for that super-accelerated phase (see Eq. (3.3)),

$\left.\frac{H}{W}\right|_{\textrm{{\^{K}}}+}=1+u_{\textrm{{\^{K}}}+}x_{\textrm{{\^{K}}}+}=1+\left(\sqrt{\frac{3}{2}}\kappa-3\right)$
(3.46)

is negative (remember that $\kappa<2\sqrt{\frac{2}{3}}$). For the rest of the
fixed points treated above, the ratio $H/W$ is $1$ (positive in the case of
ScII$\pm$ as one can check from the values of $x$ and $u$ at those fixed
points; see before Eq. (3.21)) because $u_{c}=0$. The same is true for K̂$-$
given that

$\left.\frac{H}{W}\right|_{\textrm{{\^{K}}}-}=1+u_{\textrm{{\^{K}}}-}x_{\textrm{{\^{K}}}-}=1-\left(\sqrt{\frac{3}{2}}\kappa+3\right)$
(3.47)

is always negative for $\kappa>0$.

Finally, for the second scaling solution, ScII$\pm$, the eigenvalues are

$\displaystyle
m_{1}=\frac{3}{2}\left(1+w_{\textrm{M}}\right)\left(1-\frac{\lambda}{\kappa}\right),$
(3.48) $\displaystyle
m_{2,3}=-\frac{3}{4}\left(1-w_{\textrm{M}}\right)\left\\{1\mp\sqrt{1+8\frac{1+w_{\textrm{M}}}{1-w_{\textrm{M}}}\left[1-\frac{1-3w_{\textrm{M}}}{3(1-w_{\textrm{M}})}\Delta_{\pm}-\frac{12}{\kappa^{2}}\frac{\left(1+w_{\textrm{M}}\right)^{2}}{1+3w_{\textrm{M}}}\right]}\right\\}.$
(3.49)

We see that $m_{1}$, which corresponds to the eigenvector
$\mathbf{v}=(0,1,0)$, is positive (and hence neither of the scaling fixed
points ScII$\pm$ is stable) whenever $\kappa>\lambda$ (same condition as that
of the former scaling solution; see Eq. (3.36)). The instability is in the $y$
direction. On the other hand, $m_{2,3}$ do not depend on the magnitude of
$\lambda$.

From the existence condition of the fixed points (see Eq. (3.22)) we find

$\frac{1-3w_{\textrm{M}}}{3(1-w_{\textrm{M}})}\Delta_{\pm}+\frac{12}{\kappa^{2}}\frac{(1+w_{\textrm{M}})^{2}}{1+3w_{\textrm{M}}}\leq\frac{(1-3w_{\textrm{M}})\Delta_{\pm}+2}{3(1-w_{\textrm{M}})}~{}.$
(3.50)

Also, since $1\leq\Delta_{+}<2$ and $0\leq w_{\textrm{M}}<\frac{1}{3}$ for the
fixed point ScII$+$ (see Eq. (3.24)), then

$\frac{1-3w_{\textrm{M}}}{3(1-w_{\textrm{M}})}\Delta_{+}+\frac{12}{\kappa^{2}}\frac{(1+w_{\textrm{M}})^{2}}{1+3w_{\textrm{M}}}<\frac{2}{3}\frac{2-3w_{\textrm{M}}}{1-w_{\textrm{M}}}\leq\frac{4}{3}~{}.$
(3.51)

In the case of ScII$-$, we assume that $0\leq w_{\textrm{M}}<\frac{1}{3}$ as
well, although $w_{\textrm{M}}$ was constrained to be smaller than $1$. Since
$0<\Delta_{-}\leq 1$, we have

$\frac{1-3w_{\textrm{M}}}{3(1-w_{\textrm{M}})}\Delta_{-}+\frac{12}{\kappa^{2}}\frac{(1+w_{\textrm{M}})^{2}}{1+3w_{\textrm{M}}}\leq
1~{}.$ (3.52)

Thus only in the case of ScII$-$ will $m_{2,3}$ always be real (again,
assuming $0\leq w_{\textrm{M}}<\frac{1}{3}$). For complex $m_{2,3}$, the real
part of both eigenvalues is negative and ScII$+$ is stable if
$\lambda>\kappa$. For the rest of the possibilities, that depends on the value
of $w_{\textrm{M}}$ under consideration.

|  | Fixed point: $(x_{c},y_{c},u_{c})$ | Stability
---|---|---|---
M | $(0,0,0)$ | Unstable
K$\pm$ | $(\pm 1,0,0)$ | Unstable
I | $\left(\frac{\lambda}{\sqrt{6}},\sqrt{1-\frac{\lambda^{2}}{6}},0\right)$ | Unstable whenever $\kappa>\lambda$
ScI | $\left(\sqrt{\frac{3}{2}\frac{(1+w_{\textrm{M}})^{2}}{\lambda^{2}}},\sqrt{\frac{3}{2}\frac{1-w^{2}_{\textrm{M}}}{\lambda^{2}}},0\right)$ | Unstable whenever $\kappa>\lambda$
dS | $\left(0,1,\sqrt{\frac{3}{2}}\lambda\right)$ | Stable unless $\kappa<\lambda$
ScII$\pm$ | $\left(\frac{\kappa}{\sqrt{6}}\frac{1+3w_{\textrm{M}}}{3(1-w^{2}_{\textrm{M}})}\Delta_{\pm},0,\frac{\kappa}{\sqrt{6}}\frac{1}{1+w_{\textrm{M}}}\Delta_{\pm}\right)$ | Unstable whenever $\kappa>\lambda$
K̂$\pm$ | K̂$+$ | $\left(+1,0,\sqrt{\frac{3}{2}}\kappa-3\right)$ |  Unstable whenever $\kappa<\lambda$ and $\kappa<2\sqrt{\frac{2}{3}}$
K̂$-$ | $\left(-1,0,\sqrt{\frac{3}{2}}\kappa+3\right)$ | Unstable whenever $\kappa>\lambda$
Table 1: Summary of the fixed points considered in this work and their
stability according to $\kappa$ and $\lambda$ for an exponential coupling
function (see Eq. (3.16) and remember that $\kappa\neq\lambda$ is assumed). We
determine the stability of each fixed point assuming an EoS parameter of the
background fluid in the range $0\leq w_{\textrm{M}}<1$.
$\Delta_{\pm}(\kappa,w_{\textrm{M}})$ is defined in Eq. (3.21).

A summary of all the fixed points we studied in this section, and their
stability, is provided in Table 1.

###  3.4 Numerical Simulations

In view of the fixed points shown in Sec. 3.2 and their stability, we are
interested in the trajectories that begin close to kination, pass near the
scaling regime ScI and approach the de Sitter fixed point, depicting the
evolution of the Universe as suggested by the observations. To that end, Eqs.
(3.9), (3.10) and (3.17) are solved numerically, taking the constraint
equation (3.2) into account. We select those trajectories in the phase space
that satisfy $\Omega_{\textrm{M}}=0.3147\pm 0.0074$ and
$w_{\textrm{f}}=-0.957\pm 0.080$ [3] at $N=0$, which corresponds to the
present time. To define $w_{\textrm{f}}$, we use the relation [46]

$\epsilon_{H}=\frac{3}{2}\frac{\rho_{\phi}(1+w_{\textrm{f}})+\rho_{\textrm{M}}(1+w_{\textrm{M}})}{\rho_{\phi}+\rho_{\textrm{M}}}~{}.$
(3.53)

$w_{\textrm{f}}$ can be interpreted as the EoS parameter of the DE fluid
causing the accelerated expansion of the Universe [69]. Then

$w_{\textrm{f}}=-1+\frac{1}{x^{2}+y^{2}}\left[\frac{2}{3}\epsilon_{H}-\left(1-x^{2}-y^{2}\right)\left(1+w_{\textrm{M}}\right)\right],$
(3.54)

where $\epsilon_{H}$ is given in Eq. (3.15). The background fluid EoS
parameter is chosen to be that of pressureless matter,
$w_{\textrm{M}}=w_{m}=0$ (therefore $\Omega_{\textrm{M}}=\Omega_{m}$).

As initial values, we set $u_{0}=10^{-23}$ and $y_{0}=10^{-3}$. Such a small
initial value of $u$ is imposed to ensure a negligible initial contribution
from the GB coupling, given that the trajectories in phase space begin near
the fixed points K$\pm$ and M, in accordance with the evolution prescribed by
the quintessence models. The initial value of $x$ is chosen differently within
the range $(-1,1)$. All the trajectories obtained from the simulations end up
at the de Sitter fixed point $y_{\textrm{{dS}}}=1$ (these trajectories
correspond to the different values assigned to $x_{0}$). Since $\delta
u\propto\exp\left[3\left(\frac{\kappa}{\lambda}-1\right)N\right]$ during the
scaling regime (see Eq. (3.36), where $w_{\textrm{M}}=w_{m}=0$),
$\kappa/\lambda$ can be chosen to be of order $\mathcal{O}(1)$, so it takes
longer for the Universe to abandon the scaling regime.

We see in Figs. 1a and 1b the results of some simulations for
$\kappa/\lambda=1.4$ (so the ratio $\kappa/\lambda$ is of order
$\mathcal{O}(1)$, as argued above) depicting the evolution of interest as
explained at the beginning of this section. The $3$D phase portrait is shown
in Fig. 1e. In Figs. 1a and 1e, the M fixed point has been represented by a
square. The kination fixed points K$\pm$ correspond to the left and right
leaning triangles, depending on the sign. The de Sitter fixed point is given
by the circle instead, and the scaling regime corresponds to a star onto which
the trajectories converge. It can be noticed that, before the de Sitter phase
is reached, the density parameter of the scalar field (see Eq. (3.6)) becomes
larger than $1$ given that $\Omega_{\xi}$ (see Eq. (3.7)) turns negative. We
obtained this in many other simulations, with different values of $\lambda$
and $\kappa$. Likewise, the growth of $w_{\textrm{f}}$ reaches the value of
the EoS parameter of pressureless matter and surpasses it for a brief moment
of time. There is an equally short period where it takes on phantom values
before settling down at $w_{\textrm{{dS}}}=-1$, which is the expected value
during de Sitter. Note that slightly phantom DE is preferred by the
observational data [3].

We observe another feature shared with other simulations: the density
parameter $\Omega_{\xi}$ is always the dominant contribution to DE at the
present time. This occurs unless the scalar field energy density takes over
during the scaling regime, which is discarded given that the field is supposed
to be light and may alter the successful predictions of Big Bang
Nucleosynthesis (BBN).

The Universe may enter the second scaling regime ScII$\pm$ instead, where
$y_{\textrm{{ScII}$\pm$}}=0$, starting from kination. If $\Omega_{\textrm{M}}$
is larger than the sum of the rest of the density parameters,
$\Omega_{\phi}+\Omega_{\xi}$, during the scaling regime, then
$\left.\Omega_{\textrm{M}}\right|_{\textrm{{ScII}$\pm$}}>0.5$.141414One might
argue that this is not necessarily true if $\Omega_{\xi}$ is negative.
However, this is discarded in this situation given that $\Omega_{\xi}<0$
implies $-2<ux<0$ (see Eq. (3.7)). At ScII$\pm$, we obtain
$u_{\textrm{{ScII}$\pm$}}x_{\textrm{{ScII}$\pm$}}=\frac{1+3w_{\textrm{M}}}{3(1-w_{\textrm{M}})}u_{\textrm{{ScII}$\pm$}}^{2}~{},$
(3.55) and because $0\leq w_{\textrm{M}}<1$,
$\left.\Omega_{\xi}\right|_{\textrm{{ScII}$\pm$}}$ is never negative. For
pressureless matter ($w_{\textrm{M}}=w_{m}=0$),
$\left.\Omega_{m}\right|_{\textrm{{ScII}$\pm$}}$ reads

$\left.\Omega_{m}\right|_{\textrm{{ScII}$\pm$}}=\frac{54\left(54-\kappa^{2}\Delta_{\pm}^{2}\right)}{\left(54+3\kappa^{2}\Delta_{\pm}^{2}\right)^{2}}~{}.$
(3.56)

Then $\left.\Omega_{m}\right|_{\textrm{{ScII}$\pm$}}>0.5$ implies
$\kappa\Delta_{\pm}<\sqrt{6}$. The fixed point ScII$+$ can be discarded
because $\kappa>3\sqrt{2}$ in order for the fixed point to exist, and
$1\leq\Delta_{+}$ (see Eq. (3.22) and discussion below), which means that
$\kappa\leq\kappa\Delta_{+}<\sqrt{6}$. Regarding the fixed point ScII$-$, we
have (see Eq. (3.21))

$\kappa\Delta_{-}=\kappa-\sqrt{\kappa^{2}-18}<\sqrt{6}~{},$ (3.57)

and hence $\kappa>2\sqrt{6}$ for $\Omega_{m}$ to be larger than
$\Omega_{\phi}+\Omega_{\xi}$ during the scaling regime, in agreement with the
BBN predictions.

(a) Phase portrait for $\lambda=4.47\sqrt{\frac{2}{3}}$ and
$\kappa/\lambda=1.4$. (b) $\Omega$’s and $w_{\textrm{f}}$ drawn from the red
trajectory.

(c) Phase portrait for $\lambda=0.72\sqrt{3}$ and $\kappa/\lambda=4.41$. (d)
$\Omega$’s and $w_{\textrm{f}}$ drawn from the red trajectory.

(e) $3$D phase portrait for $\lambda=4.47\sqrt{\frac{2}{3}}$ and
$\kappa/\lambda=1.4$. (f) $3$D phase portrait for $\lambda=0.72\sqrt{3}$ and
$\kappa/\lambda=4.41$.

Figure 1: Phase portraits and density parameters and $w_{\textrm{f}}$ for some
values of the parameters $\lambda$ and $\kappa$. In Fig. 1a, all the
trajectories begin either close to K$\pm$ or near M. They all converge at ScI
and transition to dS. In Fig. 1c, the red and cyan trajectories in particular
pass near the scaling regime ScII$-$ (triangle pointing downwards) before
reaching dS. The rest converge at I instead (the diamond in between K$+$ and M
in the plot). The stripe indicates the forbidden region explicitly due to the
constraint equation (3.2) and, the triangle pointing upwards, the scaling
fixed point ScII$+$. Same phase portraits in $3$D are shown in Figs. 1e and
1f.

If the fixed points ScII$\pm$ are unstable only in the $y$ direction (see
discussion below Eq. (3.49)), then $m_{2,3}$ should be complex, given that the
real part in that case is always negative (see Eq. (3.49)). Unfortunately,
this cannot happen in the case of ScII$-$ (see Eq. (3.52) and discussion
below), and only at this fixed point can $\Omega_{m}$ dominate over
$\Omega_{\phi}+\Omega_{\xi}$.

An example of trajectories going through the second scaling solution
(corresponding to the fixed point ScII$-$, which is being represented by a
triangle pointing downwards) are the red and cyan ones in Fig. 1c, where
$\kappa/\lambda=4.41$ and $u_{0}=10^{-5}$. The purple line no longer passes
near the scaling regime, but approaches the I fixed point, represented by the
diamond in the plot. The corresponding $3$D phase portrait is shown in Fig.
1f. In Fig. 1c we show the density parameters and $w_{\textrm{f}}$ associated
with the red trajectory, which are compatible with the constraints on
$\Omega_{m}$ and $w_{\textrm{f}}$. Unfortunately, the scaling regime does not
last long enough for this particular case, although we see that $\Omega_{m}$
takes on the value that we expect in view of Eq. (3.56) and the corresponding
value of $\kappa$. Also, in line with what we explained earlier, the fixed
point ScII$+$ falls within the forbidden region because $\kappa>3\sqrt{2}$.
Such a fixed point is represented by a triangle pointing upwards. The
forbidden region corresponds to the pale red coloured stripe, which excludes
values of $x>1$ that violate the constraint equation (3.2).

###  3.5 The Speed of Gravitational Waves

In this last section we examine the evolution of the $\alpha_{T}$ parameter.
We first write it in terms of the dimensionless variables defined in Eqs.
(3.1) and (3.4), and of the Hubble flow parameter $\epsilon_{H}$ (see Eqs.
(2.47) and (3.15))

$\alpha_{T}=-\frac{2}{1-2ux}\left\\{(1+2ux)(ux)^{\prime}-ux\left[1-ux-\epsilon_{H}\left(1+ux\right)\right]\right\\}.$
(3.58)

It can be readily seen that

$\left.\alpha_{T}\right|_{\textrm{{M}}}=\left.\alpha_{T}\right|_{\textrm{{K}}\pm}=\left.\alpha_{T}\right|_{\textrm{{I}}}=\left.\alpha_{T}\right|_{\textrm{{ScI}}}=\left.\alpha_{T}\right|_{\textrm{{dS}}}=0~{}.$
(3.59)

For the rest of the fixed points, the $\alpha_{T}$ parameter is given by

$\displaystyle\left.\alpha_{T}\right|_{\textrm{{ScII}$\pm$}}=-\frac{18(1+w_{\textrm{M}})(1-w^{2}_{\textrm{M}})-2\kappa^{2}(1+3w_{\textrm{M}})\Delta_{\pm}}{27(1+w_{\textrm{M}})(1-w_{\textrm{M}}^{2})-2\kappa^{2}(1+3w_{\textrm{M}})\Delta_{\pm}}\left[2-\frac{\kappa^{2}}{18}\frac{(1+3w_{\textrm{M}})(5+3w_{\textrm{M}})}{(1+w_{\textrm{M}})(1-w^{2}_{\textrm{M}})}\Delta_{\pm}\right],$
(3.60)
$\displaystyle\left.\alpha_{T}\right|_{\textrm{{\^{K}}}\pm}=-\frac{(6\mp\sqrt{6}\kappa)(4\mp\sqrt{6}\kappa)}{7\mp\sqrt{6}\kappa}~{}.$
(3.61)

$\left.\alpha_{T}\right|_{\textrm{{\^{K}}}+}$ in particular is zero if
$\kappa=3\sqrt{\frac{2}{3}}$, which does not satisfy the bound
$\kappa<2\sqrt{\frac{2}{3}}$ and falls out of the super-accelerated phase
region. The other possibility is $\kappa=2\sqrt{\frac{2}{3}}$, but in that
case $u_{c}x_{c}=-1$ and $H_{c}=0$ (see Eq. (3.3)). Since $\kappa>0$,
$\left.\alpha_{T}\right|_{\textrm{{\^{K}}}-}$ never vanishes exactly.

Regarding $\left.\alpha_{T}\right|_{\textrm{{ScII}$\pm$}}$, in the case of the
fixed point ScII$+$, the $\alpha_{T}$ parameter cannot be zero due to the
existence conditions of the fixed point derived below Eq. (3.21) and the range
of values of $w_{\textrm{M}}$ under consideration. One might have
$\left.\alpha_{T}\right|_{\textrm{{ScII}$-$}}=0$ though. However, given
$w_{\textrm{M}}=w_{m}=0$, $\left.\alpha_{T}\right|_{\textrm{{ScII}$-$}}=0$ if
$\kappa^{2}\Delta_{-}(\kappa,0)=36/5$ or $\kappa^{2}\Delta_{-}(\kappa,0)=9$
(see Eq. (3.60)). We do not obtain solutions for $\kappa$ in those cases.
Thus, $\alpha_{T}$ never vanishes at ScII$-$ when $w_{\textrm{M}}=0$.

Although the value of $|\alpha_{T}|$ at the fixed points is useful to have an
indication of its evolution, the current state of the Universe cannot be
described by any of those fixed points. Moreover, $|\alpha_{T}|$ does not have
to be exactly zero, instead it must satisfy the bound in Eq. (1.2). Then, we
show the evolution of $\alpha_{T}$ connecting the different fixed points and
make sure it does not surpass the upper bound from the observational
constraint.

Figure 2:  $|\alpha_{T}|$ predicted by the WGB model for some values of the
parameters $\lambda$ and $\kappa$ and different initial conditions of $x$. The
stripe in the smallest graph on the left indicates the excluded region
starting at redshift $z=0.00980\pm 0.00079$ (see Ref. [45]). We see that the
bound $|\alpha_{T}|<10^{-15}$ is grossly violated for the cases under
consideration. These cases have been selected given that they did satisfy the
constraints on $\Omega_{\textrm{M}}$ and $w_{\textrm{f}}$.

That bound is applicable in the relatively recent history of the Universe. To
be precise, $|\alpha_{T}|$ is constrained from redshift151515If one consults
Ref. [72] the redshift value, $z=0.0099$, is the geocentric redshift of the
host galaxy. The cosmological redshift by contrast, which is the one shown
here and what we need in this case, is given in Ref. [45]. $z=0.00980\pm
0.00079$ until today ($z=0$).161616Do not confuse ‘$z$’ here denoting the
redshift with the dimensionless variable defined in Eq. (3.1) and related with
the energy density of the background fluid. Selecting models which reproduce
trajectories that go through the scaling solution ScI, before entering de
Sitter, with different values of the parameters $\kappa$ and $\lambda$, and
different initial condition $x_{0}$, we plot the evolution of the
$|\alpha_{T}|$ parameter over the number of e-folds $N$ (see Fig. 2). In that
plot, we zoom in on the excluded region, marked by a vertical, pale red
coloured stripe, and we observe that the bound is violated by many orders of
magnitude, similarly to the SGB model [46]. The cases shown in the plot were
chosen given that they did satisfy the constraints on $\Omega_{m}$ and
$w_{\textrm{f}}$ at the present time.

If we examine Eq. (3.58), we can understand why the bound is so strongly
violated. We notice that for $ux$ of order
$\mathcal{O}(0.1)$-$\mathcal{O}(1)$, then $\Omega_{\xi}$ is of order
$\mathcal{O}(0.1)$ (see Eq. (3.7)). We observed numerically that
$\Omega_{\xi}$ is always the density parameter that dominates over the rest by
the time $\Omega_{m}$ and $w_{\textrm{f}}$ satisfy the constraints. In
addition to this, we remark that it grows at that time. Then $ux\rightarrow
1/2$ (i.e. $\Omega_{\xi}\rightarrow 0.6$) while the numerator in Eq. (3.58)
remains of order $\mathcal{O}(0.1)$-$\mathcal{O}(1)$. This causes
$|\alpha_{T}|$ to increase near the present time, although it does not diverge
because $(ux)^{\prime}$ increases eventually as well, even though
$\Omega_{\xi}$ remains constant. This can be seen by simply taking the
derivative of Eq. (3.7)

$\Omega_{\xi}^{\prime}=\frac{2(1-\Omega_{\xi})}{1+ux}(ux)^{\prime}.$ (3.62)

When $\Omega_{\xi}\rightarrow 1$, $\Omega_{\xi}^{\prime}\rightarrow 0$.

The fact that $|\alpha_{T}|$ does not remain below the observational bound
$10^{-15}$ while we approach the present time indicates that the exponential
coupling function in the WGB model must be discarded. Similar results are
obtained in the SGB model [46].

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

In this section, instead of specifying the functional form of $\xi(\phi)$ and
consider examples other than the exponential case, we impose the condition
$\alpha_{T}=0$. In the SGB model, such a condition translates into Eq. (1.4).
In the WGB model, however, $\alpha_{T}=0$ leads to the following equation (see
Eq. (2.47)):

$\ddot{\xi}=\dot{\xi}W\left(1-2\frac{\dot{\xi}W}{M^{2}_{\textrm{Pl}}}\right).$
(3.63)

Using Eq. (2.26), Eq. (3.63) becomes

$\ddot{\xi}=-\frac{\dot{\xi}H}{2\mu}\left(1+2\mu-\sqrt{1+8\mu}\right).$ (3.64)

Now, the equation of $\dot{H}$ in the SGB model when $\ddot{\xi}=H\dot{\xi}$
reads (see Eq. (2.34))

$-2\dot{H}\left(M^{2}_{\textrm{Pl}}-4\dot{\xi}H\right)=\dot{\phi}^{2}+\rho_{\textrm{M}}\left(1+w_{\textrm{M}}\right).$
(3.65)

This implies that in de Sitter, $\dot{\phi}=0$ (so $\dot{H}=0$). However, this
value of $\dot{\phi}$ is problematic when we integrate Eq. (1.4). Indeed, the
integral on the LHS of the equation

$\int_{0}^{\dot{\xi}_{1}}\textrm{d}\ln(|\dot{\xi}|)=\int_{a_{\textrm{{dS}}}}^{a_{1}}\textrm{d}\ln
a~{},$ (3.66)

does not converge. In the WGB model however, since we have $\Phi$ in addition
to $H$, there is a chance that $\dot{H}=0$ but still $\dot{\phi}\neq 0$. To
verify this, we must write an equation analogous to Eq. (3.65). Bearing Eq.
(3.64) in mind, we write Eq. (2.30) as

$-\frac{2M^{2}_{\textrm{Pl}}\dot{H}}{\sqrt{1+8\mu}}=\dot{\phi}^{2}+\rho_{\textrm{M}}(1+w_{\textrm{M}})-4\ddot{\xi}H^{2}\frac{(1-\sqrt{1+8\mu})^{2}}{16\mu^{2}\sqrt{1+8\mu}}-4\dot{\xi}H^{3}\frac{(1-\sqrt{1+8\mu})^{3}}{64\mu^{3}}~{},$
(3.67)

where Eq. (2.26) was used. Plugging Eq. (3.64) into Eq. (3.67), we obtain

$-\frac{2M^{2}_{\textrm{Pl}}\dot{H}}{\sqrt{1+8\mu}}=\dot{\phi}^{2}+\rho_{\textrm{M}}\left(1+w_{\textrm{M}}\right)+3H^{2}M^{2}_{\textrm{Pl}}\frac{\left(1-\sqrt{1+8\mu}\right)^{2}}{16\mu^{2}\sqrt{1+8\mu}}\left(1+4\mu-\sqrt{1+8\mu}\right).$
(3.68)

We see that the last term vanishes when $|\mu|\ll 1$, recovering Eq (3.65) for
$|\dot{\xi}H|\ll M_{\textrm{Pl}}^{2}$ ($\mu$ was defined in Eq. (2.25)).
Unfortunately, each term on the RHS is non-negative meaning that $\dot{H}=0$
implies $\dot{\phi}=0$ and hence $\dot{\xi}=0$, same as in the SGB model.

Knowing this, we proceed to analyse the stability of the de Sitter and scaling
fixed points, dS and ScI respectively, when the new condition $\alpha_{T}=0$
is imposed, in line with was done in Ref. [46]. Firstly, we write Eq. (3.63)
as

$(ux)^{\prime}=\frac{ux}{1+ux}\left[1-\left(3x+2u\right)x-\frac{3}{2}\left(1+w_{\textrm{M}}\right)\left(1-x^{2}-y^{2}\right)\right],$
(3.69)

where we used Eq. (3.13), the constraint equation (3.2) and Eqs. (3.4) and
(3.12), given that we do not assume a coupling function form in particular.
This is equivalent to $\alpha_{T}=0$ in Eq. (3.58). We however make it clear
that we use Eq. (3.12) to write it that way.

We perturb Eq. (3.69) around the de Sitter fixed point (which is still a fixed
point of the system despite Eq. (3.69)), that can only correspond to $x=0$ and
$z=0$ as we just argued

$\delta x^{\prime}=\delta x~{}.$ (3.70)

This equation implies that $\delta x\propto e^{N}$. The rest of the equations
read171717Notice that, although we have not assumed a specific coupling
function $\xi(\phi)$ in this section,
$u_{\textrm{{dS}}}=\sqrt{\frac{3}{2}}\lambda$, as in the case of the
exponential function. This can be seen from Eq. (3.9) with
$x_{\textrm{{dS}}}=0$ and $y_{\textrm{{dS}}}=1$.

$\displaystyle\delta x^{\prime}=-3\delta
x-\sqrt{\frac{3}{2}}\lambda\left(1+3w_{\textrm{M}}\right)\delta y-\delta
u~{},$ (3.71) $\displaystyle\delta y^{\prime}=-3(1+w_{\textrm{M}})\delta
y~{}.$ (3.72)

The second equation gives $\delta y\propto e^{-3(1+w_{\textrm{M}})N}$. Using
Eq. (3.70), Eq. (3.71) gives $\delta u$. Since $\delta y$ dies out but $\delta
x$ grows, $\delta u$ will increase too. Then we see that the de Sitter fixed
point dS is unstable, as in the SGB model [46].

Regarding the scaling regime ScI, since $u_{\textrm{{ScI}}}=0$, this is still
a fixed point despite Eq. (3.69). Linearising Eq. (3.69), we obtain

$\delta u^{\prime}=-\frac{1}{2}\left(1+3w_{\textrm{M}}\right)\delta u~{}.$
(3.73)

Hence $\delta u\propto e^{-\frac{1}{2}(1+3w_{\textrm{M}})N}$, and given the
eigenvalues in Eq. (3.37) (which remain the same despite $\alpha_{T}=0$), we
have that ScI becomes stable, a result similar to the one obtained in the SGB
model in Ref. [46] and in quintessence models in GR. When it comes to the
second scaling regime, where $y_{\textrm{{ScII}$\pm$}}=0$, we can see that
they are not even fixed points of the new system of equations that includes
Eq. (3.69).

Consequently, the WGB model with $\alpha_{T}=0$ cannot reproduce a long
scaling regime with small $u$ followed by a period of accelerated expansion,
and hence can be discarded.

##  4 Summary and Conclusions

Gauss-Bonnet (GB) gravity with $\xi(\phi)\mathcal{G}$ (the scalar-Gauss-Bonnet
(SGB) model), where $\mathcal{G}$ is the GB term, has become a popular dark
energy (DE) model and it is known to predict the propagation speed of
gravitational waves (GWs) $c_{\textrm{GW}}\neq 1$ (in units where $c=1$). For
simple coupling functions $\xi(\phi)$ which reproduce the past evolution of
the Universe as we know it, the SGB model does not satisfy the observational
constraints on $|\alpha_{T}|\equiv|c_{\textrm{GW}}^{2}-1|$ from GRB170817A
[37] when applied to the late Universe in particular [46]. The same occurs
with Horndeski’s theory [38]. The model is actually a subclass of that theory
and one can recover the SGB action for a non-trivial choice of the functions
of the theory [40].

In spite of those well-known results, it was shown recently that Horndeski’s
theory where $G(\phi,X)R$ with
$X\equiv-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi$ does
predict GWs propagating at the speed of light in vacuum when formulated in the
Palatini formalism [53]. All the previous claims of $c_{\textrm{GW}}\neq 1$
were made assuming the metric case instead, where the connection is the Levi-
Civita (LC) one, which fully depends on the metric and its first derivative.
In the Palatini formalism however, the connection is independent of the metric
and obeys some field equations, meaning that the geometry differs from the
pseudo-Riemannian one. Given this, in the present work we considered the SGB
action but in Weyl geometry where the connection (the Weyl connection) has
zero torsion but it is not metric compatible. We named this model Weyl-Gauss-
Bonnet. A similar action was analysed in Ref. [59] but ours considers the
field-dependent coupling to the GB term and restricts to four dimensions.
Thus, the goal in our work was to investigate whether the SGB gravity, as a DE
model, when formulated in a different geometry, fulfilled the constraint (1.2)
and yet displayed a period of matter domination followed by an accelerated
expanding phase of the Universe in accordance with other constraints on the
density and equation of state (EoS) parameters.

We assumed a spatially-flat FLRW metric and a background Weyl vector with
temporal component $\Phi(t)$ and no spatial components breaking spatial
isotropy. The resulting action once the LC and non-Riemannian parts were
written separately was given in Eq. (2.19). The $\alpha_{T}$ parameter was
calculated in Sec. 2.4. In Sec. 3, we performed a dynamical systems analysis
assuming an exponential potential, which is common in quintessence scenarios.
The dimensionless variables in Eqs. (3.1) and (3.4) were formulated in terms
of the newly defined Weyl parameter $W\equiv H-\Phi$, thereby including the
modifications to Einstein’s gravity in contrast to the metric formalism where
these are accounted for by $u$ alone. Fixed points of the dynamical system
were found in Sec. 3.2 under the assumption that $\xi(\phi)$ is an exponential
function (see Eq. (3.16)). This choice makes all the equations self-similar so
the explicit dependence on $W$ disappears. Another possibility was the linear
GB coupling (because $\xi_{,\phi\phi}=0$; see Eq. (3.12)), but it was argued
that this did not lead to the expected evolution of the Universe, as was
noticed in the metric formalism as well [46]. A stability analysis of the
fixed points was carried out in Sec. 3.3.

Among the fixed points calculated in this work, the scaling ScI (where $u=0$)
and de Sitter dS ones had been found in the SGB model of the metric formalism
too [46], and were shown to satisfy the same stability conditions in both the
WGB and SGB models. While the scaling fixed point is unstable if
$\kappa>\lambda$, dS becomes an attractor and the Universe is predicted to
abandon a period of matter domination and enter a regime of accelerated
expansion, where $\kappa$ and $\lambda$ are the strength of the exponentials
involved in the GB coupling and the scalar potential, as shown in Eqs. (3.16)
and (3.8) respectively. If $\kappa<\lambda$ however, the roles are switched
over and de Sitter becomes unstable (a saddle point) whereas ScI is stable. On
top of these and other known fixed points, we found a second scaling regime
consisting of two fixed points ScII$\pm$ (where $y=0$) and a regime of
‘pseudo-kination’ where the Weyl parameter, $W_{c}$, satisfied
$3W_{c}^{2}M^{2}_{\textrm{Pl}}=\dot{\phi}^{2}_{c}/2$ but $W_{c}\neq H_{c}$.
This regime is associated with two fixed points K̂$\pm$ as well. In the case
of ScII$+$, it was shown not to exist if the EoS parameter of the background
fluid fell within the range $\frac{1}{3}\leq w_{\textrm{M}}<1$, thereby
excluding background radiation but including pressureless matter
$w_{\textrm{M}}=0$. Regarding the pseudo-kination regime, we found that K̂$+$
led to a universe that accelerated faster than during de Sitter, where it
accelerates exponentially, if $\kappa<2\sqrt{\frac{2}{3}}$. Also, regardless
of the sign of $x_{\textrm{{\^{K}}$\pm$}}$, if $\kappa\gg 1$, then the
Universe was predicted not to accelerate. Moreover, both situations, the
super-accelerated and uniform expansion phases, were shown to be stable if
$\kappa>\lambda$ and $\kappa<\lambda$, respectively. The scaling regime
ScII$\pm$ was determined to be unstable whenever $\kappa>\lambda$, same as in
the case of ScI. The instability in the former case was found to be in the $y$
direction of the phase space. We found also that ScII$+$ is always stable only
if $\lambda>\kappa$ because two of the three eigenvalues, which do not depend
on $\lambda$, are complex with negative real parts. The summary of the fixed
points and their stability was provided in Table 1.

In Sec. 3.4 we performed some numerical simulations assuming pressureless
background matter $w_{\textrm{M}}=w_{m}=0$ and checked trajectories crossing
both scaling regimes. When it comes to those approaching ScI, we plotted the
2D phase portrait in Fig. 1a and the 3D one in Fig. 1e for some values of
$\lambda$ and $\kappa$. These were chosen among other simulations given that
they satisfied the constraints on $\Omega_{m}$ and $w_{\textrm{f}}$ [3] at
present time ($N=0$), where $w_{\rm f}$ is the EoS parameter of the DE fluid
responsible for the current accelerated expansion of the Universe. The
evolution of the density parameters and $w_{\textrm{f}}$ for one of the
trajectories in the phase portrait was displayed in Fig. 1b. The system’s
behaviour is reproduced as expected: starting from a period of kination, the
Universe enters a long period of background fluid domination which gives way
to the present epoch of DE domination. An eventual period of de Sitter is
predicted after a transient epoch where the effective density parameter
associated with the GB coupling $\xi(\phi)$, $\Omega_{\xi}$, dominates and
$w_{\textrm{f}}$ takes on phantom values (i.e. $w_{\textrm{f}}<-1$). We
pointed out that a common feature of all the simulations was that, during DE
domination, $\Omega_{\xi}$ was always dominant over the density parameter of
the scalar field $\Omega_{\phi}$.

A similar but much shorter period of matter domination was reproduced in the
case of the second scaling solution ScII$-$ (see Fig. 1d). The corresponding
2D phase diagram was shown in Fig. 1c and the 3D phase space in Fig. 1f. In
order for $\Omega_{m}$ to dominate over $\Omega_{\phi}+\Omega_{\xi}$ during
the scaling regime (so BBN predictions are not affected), it was demonstrated
that $\kappa>2\sqrt{6}$. As to ScII$+$, the bound translated into
$\kappa<\sqrt{6}$ and was at odds with the existence condition of that fixed
point. Indeed, we saw that ScII$+$ was located in the excluded region of the
phase portrait 1c (the coloured stripe).

The $\alpha_{T}$ parameter was written in terms of the dimensionless
variables, defined in Eqs. (3.1) and (3.4), their derivatives, and the Hubble
flow parameter $\epsilon_{H}$ in Eq. (3.58). We noticed that it is always non-
vanishing for K̂$-$, and for K̂$+$ only if one considers the super-accelerated
regime. Also, $\alpha_{T}$ never vanishes at ScII$-$ if $w_{\textrm{M}}=0$.
However, it was argued that the value of $|\alpha_{T}|$ at the fixed points
was not enough because the current Universe had to be somewhere in between the
scaling and de Sitter fixed points in the phase space. We plotted in Fig. 2
the evolution of $|\alpha_{T}|$ for different values of the constants
$\lambda$ and $\kappa$ and initial conditions $x_{0}$. Unfortunately, the
constraint in Eq. (1.2) was found to be grossly violated at present time, when
$\Omega_{m}$ and $w_{\textrm{f}}$ satisfy the aforementioned constraints. It
is to that moment that Eq. (1.2) applies approximately given that the binary
neutron star system’s merger occurred in the relatively recent past history of
the Universe. This was indicated in an inset plot, where the constrained
region corresponds to the coloured narrow stripe.

Despite the negative results, similar to the ones we obtained in [46] for the
SGB model in the metric formalism, the analysis performed in this work did not
consider coupling functions other than the exponential and linear cases (so
the equations are self-similar), which might lead to different conclusions. We
then regarded the case of $\alpha_{T}=0$, without making any assumption on the
coupling function, and analysed the stability of the scaling and de Sitter
fixed points in particular given their importance in late time cosmology. We
found that the latter was only possible if $\dot{\phi}=0$, same as in the SGB
model, despite the presence of the new variable $\Phi(t)$ in the equations.
The subsequent analysis of the stability indicated that, whereas the de Sitter
regime became unstable, the scaling one was stable, contradicting the known
evolution of the Universe. Consequently, when $\alpha_{T}=0$ for a non-
constant coupling function $\xi(\phi)$, we encountered the same stability of
ScI and dS as that in Ref. [46].

Assuming that the Weyl vector has vanishing spatial components, one can apply
the very same results of this work to a connection compatible with the metric
but with non-zero torsion. We showed this in App. C. In App. A we introduced
such a connection and the geometric structures used in this manuscript
briefly. The background equations used in Sec. 2 and those including the
homogeneous spatial components of the Weyl vector were detailed in App. B.

As mentioned above, we have regarded only exponential and linear coupling
functions in this work, given that the equations were self-similar and the
dynamical systems analysis could be easily done. On top of this, the case of
$\kappa=\lambda$, which was analytically considered in Ref. [73] in the metric
case, was not tackled here (because, again, it was not amenable to the
presented methods of analysis). It would be interesting to apply an analysis
similar to that of the referenced paper, although two fixed ‘curves’ were
calculated in Sec. 3.2 already, one of them corresponding to a third scaling
solution of cosmological interest. Also, despite that the homogeneous spatial
components of the Weyl vector were fully neglected under the assumption of
spatial isotropy, one might regard the presence of such components so that
they are relevant in the dynamical analysis but subdominant enough not to
generate a large scale anisotropy. During inflation, this and even the fully
dominant presence of the vector field were shown to be consistent with
observations in Ref. [74]. We leave these questions for future publications.

With this work, we hope to have brought attention to the fact that SGB gravity
could be rescued in the Palatini formalism with the aim of addressing the late
time evolution of the Universe, where the stringent observational bounds on
the speed of propagation of GWs are an unavoidable test for scalar-tensor
theories at large.

## Acknowledgments

We are grateful to T. Koivisto for useful comments while this work was being
completed. K.D. is supported (in part) by the STFC consolidated grant:
ST/X000621/1. 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 and PID2022-139841NB-I00
(AEI/FEDER, UE). 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 Palatini Formalism and Curvature Tensors

In this section of the Appendix, we turn our attention to various geometric
structures in the Palatini formalism. The Riemann tensor is defined only with
respect to the connection and its first derivatives as [63]

$\tensor{R}{{}^{\alpha}_{\beta\mu\nu}}\equiv\partial_{\mu}\Gamma^{\alpha}_{\nu\beta}-\partial_{\nu}\Gamma^{\alpha}_{\mu\beta}+\Gamma^{\alpha}_{\mu\lambda}\Gamma^{\lambda}_{\nu\beta}-\Gamma^{\alpha}_{\nu\lambda}\Gamma^{\lambda}_{\mu\beta}~{}.$
(A.1)

In the metric formalism, the above expression can be written in terms of the
second derivatives of the metric by virtue of Eq. (2.1). In the Palatini
formalism, on the other hand, the Riemann tensor does not depend on the
derivatives of the metric tensor. However, in that case, we can write Eq.
(A.1) as the sum of the LC part, and the distortion tensor (defined in Eq.
(2.4)) and its derivatives

$\tensor{R}{{}^{\alpha}_{\beta\mu\nu}}=\tensor{\mathring{R}}{{}^{\alpha}_{\beta\mu\nu}}+\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}~{},$
(A.2)

where we introduced a tensor $\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}$,
given by

$\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}\equiv\mathring{\nabla}_{\mu}\tensor{\kappa}{{}_{\nu\beta}^{\alpha}}-\mathring{\nabla}_{\nu}\tensor{\kappa}{{}_{\mu\beta}^{\alpha}}+\tensor{\kappa}{{}_{\mu\lambda}^{\alpha}}\tensor{\kappa}{{}_{\nu\beta}^{\lambda}}-\tensor{\kappa}{{}_{\nu\lambda}^{\alpha}}\tensor{\kappa}{{}_{\mu\beta}^{\lambda}}.$
(A.3)

An important second-rank tensor that is related to the Riemann one is the
Ricci tensor, which is defined as

$R_{\mu\nu}\equiv\tensor{R}{{}^{\sigma}_{\mu\sigma\nu}}=\mathring{R}_{\mu\nu}+\kappa_{\mu\nu}~{},$
(A.4)

where $\kappa_{\mu\nu}\equiv\tensor{\kappa}{{}^{\sigma}_{\mu\sigma\nu}}$.
While $\mathring{R}_{\mu\nu}$ is symmetric, $R_{\mu\nu}$ does not possess any
symmetry _a priori_.

Another related tensor, which in the metric formalism coincides with the Ricci
tensor, is the co-Ricci tensor

$\tilde{R}_{\mu\nu}\equiv
g^{\alpha\beta}R_{\mu\alpha\nu\beta}=\mathring{R}_{\mu\nu}+\tilde{\kappa}_{\mu\nu}~{},$
(A.5)

where $\tilde{\kappa}_{\mu\nu}\equiv
g^{\alpha\beta}\kappa_{\mu\alpha\nu\beta}$.

The last second-rank tensor that can be obtained from the Riemann tensor is
the ‘homothetic curvature tensor’, $\bar{R}_{\mu\nu}$. In the Palatini
formalism, the Riemann tensor is not generically antisymmetric in its first
two indices. The only symmetry it has is the antisymmetry in the last two
indices, which can be inferred from Eq. (A.1). Therefore, $\bar{R}_{\mu\nu}$
is defined as

$\bar{R}_{\mu\nu}\equiv\tensor{R}{{}^{\sigma}_{\sigma\mu\nu}}=\bar{\kappa}_{\mu\nu}~{},$
(A.6)

where

$\bar{\kappa}_{\mu\nu}\equiv\tensor{\kappa}{{}^{\sigma}_{\sigma\mu\nu}}=\mathring{\nabla}_{\mu}\tensor{\kappa}{{}_{\nu\sigma}^{\sigma}}-\mathring{\nabla}_{\nu}\tensor{\kappa}{{}_{\mu\sigma}^{\sigma}}.$
(A.7)

This tensor is fully antisymmetric. Consequently, the Ricci scalar

$R\equiv
g^{\mu\nu}R_{\mu\nu}=\tensor{\tilde{R}}{{}^{\sigma}_{\sigma}}=\mathring{R}+\kappa~{},$
(A.8)

is the only non-vanishing scalar one can extract from the Riemann tensor. In
the above expression, the scalar $\kappa$ is defined as

$\kappa\equiv
g^{\mu\nu}\kappa_{\mu\nu}=\tensor{\tilde{\kappa}}{{}^{\sigma}_{\sigma}}~{}.$
(A.9)

The Weyl connection considered in Eq. (2.9) is a concrete example of a more
general one with distortion tensor

$\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}=c_{1}\delta_{\mu}^{\alpha}A_{\nu}+c_{2}\delta_{\nu}^{\alpha}A_{\mu}-c_{3}g_{\mu\nu}A^{\alpha},$
(A.10)

where $c_{1}$, $c_{2}$ and $c_{3}$ are constants. The torsion and non-
metricity tensors read (see Eqs. (2.7) and (2.8), respectively)

$\displaystyle\tensor{T}{{}_{\mu\nu}^{\alpha}}=\left(c_{1}-c_{2}\right)\left(\delta_{\mu}^{\alpha}A_{\nu}-\delta_{\nu}^{\alpha}A_{\mu}\right),$
(A.11) $\displaystyle
Q_{\alpha\mu\nu}=-2c_{2}g_{\mu\nu}A_{\alpha}-\left(c_{1}-c_{3}\right)\left(g_{\mu\alpha}A_{\nu}+g_{\nu\alpha}A_{\mu}\right).$
(A.12)

We see that for $(c_{1},c_{2},c_{3})=(1,1,1)$, the torsion tensor vanishes
although the non-metricity tensor is not zero, leading to Eq. (2.10). Besides,
the choice $(c_{1},c_{2},c_{3})=(1,0,1)$ leads to a connection compatible with
the metric (i.e. with vanishing non-metricity tensor), but with non-zero
torsion. A connection with non-zero torsion and non-metricity tensors
corresponds to $(c_{1},c_{2},c_{3})=(0,1,0)$. This latter connection is well-
known for being the most general solution to the connection field equations
derived from the Einstein-Hilbert action in the Palatini formalism [75, 76].

In the case of the Weyl connection, the contributions from the distortion
tensor to the Ricci, co-Ricci and homothetic curvature tensors read

$\displaystyle\kappa_{\mu\nu}=\mathring{\nabla}_{\mu}A_{\nu}-3\mathring{\nabla}_{\nu}A_{\mu}+2A_{\mu}A_{\nu}-g_{\mu\nu}\left(\mathring{\nabla}_{\sigma}A^{\sigma}+2A_{\sigma}A^{\sigma}\right),$
(A.13)
$\displaystyle\tilde{\kappa}_{\mu\nu}=-\mathring{\nabla}_{\mu}A_{\nu}-\mathring{\nabla}_{\nu}A_{\mu}+2A_{\mu}A_{\nu}-g_{\mu\nu}\left(\mathring{\nabla}_{\sigma}A^{\sigma}+2A_{\sigma}A^{\sigma}\right),$
(A.14)
$\displaystyle\bar{\kappa}_{\mu\nu}=4\left(\mathring{\nabla}_{\mu}A_{\nu}-\mathring{\nabla}_{\nu}A_{\mu}\right),$
(A.15)

respectively. One can easily show that the relation among the three tensors
above can be written as

$\kappa_{\mu\nu}=\tilde{\kappa}_{\mu\nu}+\frac{1}{2}\bar{\kappa}_{\mu\nu}~{}.$
(A.16)

Notice that $\tilde{\kappa}_{\mu\nu}$ (and thus $\tilde{R}_{\mu\nu}$) is
symmetric in this particular case. Moreover, $\tilde{\kappa}_{\mu\nu}$ is the
symmetric part of $\kappa_{\mu\nu}$; i.e.
$\kappa_{(\mu\nu)}=\tilde{\kappa}_{\mu\nu}$. The antisymmetric part, on the
other hand, is given by $\kappa_{[\mu\nu]}=\frac{1}{2}\bar{\kappa}_{\mu\nu}$.

The scalar $\kappa$, defined in Eq. (A.9), becomes

$\kappa=-6\left(\mathring{\nabla}_{\sigma}A^{\sigma}+A_{\sigma}A^{\sigma}\right).$
(A.17)

Regarding $\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}$, we have the following
symmetry relations:

$\displaystyle\kappa_{\alpha\beta\mu\nu}=\bar{\kappa}_{\alpha\beta\mu\nu}+\frac{1}{4}g_{\alpha\beta}\bar{\kappa}_{\mu\nu}~{},$
(A.18)
$\displaystyle\kappa_{\alpha\beta\mu\nu}-\kappa_{\mu\nu\alpha\beta}=\frac{1}{4}\left(g_{\alpha\beta}\bar{\kappa}_{\mu\nu}-g_{\mu\nu}\bar{\kappa}_{\alpha\beta}-g_{\beta\nu}\bar{\kappa}_{\mu\alpha}+g_{\mu\alpha}\bar{\kappa}_{\beta\nu}-g_{\mu\beta}\bar{\kappa}_{\alpha\nu}+g_{\alpha\nu}\bar{\kappa}_{\mu\beta}\right),$
(A.19)

where $\bar{\kappa}_{\alpha\beta\mu\nu}\equiv\kappa_{[\alpha\beta]\mu\nu}$.

Among the various options discussed below Eq. (A.12), the second connection
with torsion is addressed in App. C.

##  Appendix B Homogeneous Spatial Vector

For completeness, we derive the homogeneous equations when the Weyl vector has
non-zero homogeneous spatial components. This induces background anisotropy,
as expected. To that end, we vary the action with respect to the inverse
metric first.

The variation of the scalar-Gauss-Bonnet part of the action in Eq. (2.19) is

$\displaystyle\delta
S_{\textrm{SGB}}=\frac{M^{2}_{\textrm{Pl}}}{2}\int\textrm{d}^{4}x\sqrt{-g}\left\\{G_{\mu\nu}-\frac{4}{M^{2}_{\textrm{Pl}}}\left(\tensor{R}{{}_{\mu}^{\alpha}{}_{\nu}^{\beta}}+2\tensor{S}{{}^{\alpha}_{\mu}}\delta_{\nu}^{\beta}-g_{\mu\nu}G^{\alpha\beta}-R_{\mu\nu}g^{\alpha\beta}\right)\nabla_{\alpha}\partial_{\beta}\xi-\right.$
$\displaystyle\left.-\frac{1}{M^{2}_{\textrm{Pl}}}\left[\partial_{\mu}\phi\partial_{\nu}\phi-\frac{1}{2}g_{\mu\nu}\left(g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi+2V\right)\right]\right\\}\delta
g^{\mu\nu},$ (B.1)

where $S_{\mu\nu}$ is the traceless tensor defined as

$S_{\mu\nu}\equiv R_{\mu\nu}-\frac{1}{4}g_{\mu\nu}R~{}.$ (B.2)

We remind the reader we dropped the rings (overcircles) here because, as was
argued below Eq. (2.19), every geometric quantity is defined with respect to
the LC connection once the LC and non-Riemannian parts (the latter depending
on the Weyl vector and its derivatives) have been written separately. Then

$\displaystyle\delta
S_{\textrm{SGB}}=\frac{M^{2}_{\textrm{Pl}}}{2}\int\textrm{d}^{4}x\sqrt{-g}\left\\{\left[3H^{2}-12\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}H^{3}-\frac{1}{M^{2}_{\textrm{Pl}}}\rho_{\phi}\right]\delta
g^{00}-\left[2\dot{H}+3H^{2}-4\frac{\ddot{\xi}}{M^{2}_{\textrm{Pl}}}H^{2}-\right.\right.$
$\displaystyle\left.\left.-8\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}H^{3}-8\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}H\dot{H}+\frac{1}{M^{2}_{\textrm{Pl}}}P_{\phi}\right]g_{ij}\delta
g^{ij}\right\\},$ (B.3)

where the spatially-flat FLRW metric was assumed, and $\rho_{\phi}$ and
$P_{\phi}$ are defined as the energy density and pressure of the canonical
scalar field $\phi$

$\displaystyle\rho_{\phi}\equiv\frac{1}{2}\dot{\phi}^{2}+V(\phi)~{},$ (B.4)
$\displaystyle P_{\phi}\equiv\frac{1}{2}\dot{\phi}^{2}-V(\phi)~{},$ (B.5)

respectively.

Regarding the action of the vector field $A_{\mu}$, which we shall call
‘$S_{A}$’, such that $S_{\textrm{WGB}}=S_{\textrm{SGB}}+S_{A}+S_{\textrm{M}}$
(see Eq. (2.20) where $S_{\textrm{M}}$ was introduced) and

$\displaystyle
S_{A}=-4\int\textrm{d}^{4}x\sqrt{-g}\left[\frac{3}{4}M^{2}_{\textrm{Pl}}A_{\sigma}A^{\sigma}-\left(G^{\mu\nu}-\nabla^{\mu}A^{\nu}\right)A_{\mu}\partial_{\nu}\xi-\left(\nabla_{\sigma}A^{\sigma}+A_{\sigma}A^{\sigma}\right)\partial_{\rho}\xi
A^{\rho}-\right.$ $\displaystyle\left.-\frac{1}{2}\xi\Upsilon
F^{\mu\nu}F_{\mu\nu}\right],$ (B.6)

where we defined

$\displaystyle
F_{\mu\nu}\equiv\nabla_{\mu}A_{\nu}-\nabla_{\nu}A_{\mu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}~{},$
(B.7)

i.e. $F_{\mu\nu}=\frac{1}{4}\bar{\kappa}_{\mu\nu}$, we have

$\displaystyle\delta
S_{A}=-4\int\textrm{d}^{4}x\sqrt{-g}\left\\{\left[\frac{3}{8}M^{2}_{\textrm{Pl}}\Phi^{2}-\frac{1}{2}\dot{\xi}Wa^{-2}A_{i}A^{i}-\frac{9}{2}\dot{\xi}H\Phi
W-\frac{3}{2}\dot{\xi}\Phi^{3}-\frac{1}{2}\xi\Upsilon
a^{-2}\dot{A_{i}}\dot{A^{i}}+\right.\right.$
$\displaystyle\left.\left.+\frac{3}{8}M^{2}_{\textrm{Pl}}a^{-2}A_{i}A^{i}\right]\delta
g^{00}-\frac{3}{2}M^{2}_{\textrm{Pl}}\left[\Phi-2\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}\left(W^{2}-\frac{1}{3}a^{-2}A_{j}A^{j}\right)\right]A_{i}\delta
g^{0i}+\left[\frac{1}{2}\left(\frac{3}{2}M^{2}_{\textrm{Pl}}+\right.\right.\right.$
$\displaystyle\left.\left.\left.+\ddot{\xi}+3H\dot{\xi}-2\dot{\xi}\Phi\right)A_{i}A_{j}+\dot{\xi}\dot{A_{i}}A_{j}+\xi\Upsilon\dot{A_{i}}\dot{A_{j}}+\left(\ddot{\xi}H\Phi-\frac{1}{2}\ddot{\xi}\Phi^{2}+\frac{3}{8}M^{2}_{\textrm{Pl}}\Phi^{2}-\right.\right.\right.$
$\displaystyle\left.\left.\left.-\frac{3}{8}M^{2}_{\textrm{Pl}}a^{-2}A_{k}A^{k}+\dot{\xi}H\dot{\Phi}+\dot{\xi}\dot{H}\Phi-\dot{\xi}\Phi\dot{\Phi}-\frac{1}{2}\dot{\xi}\Phi^{3}+\frac{3}{2}\dot{\xi}H^{2}\Phi-\frac{1}{2}\dot{\xi}Wa^{-2}A_{k}A^{k}-\right.\right.\right.$
$\displaystyle\left.\left.\left.-\frac{1}{2}\xi\Upsilon
a^{-2}\dot{A_{k}}\dot{A^{k}}\right)g_{ij}\right]\delta g^{ij}\right\\},$ (B.8)

after some tedious calculations. $A_{i}(t)$ are the homogeneous spatial
components of the Weyl vector; i.e.

$A_{\mu}(t)\equiv\left(-\Phi(t),A_{i}(t)\right).$ (B.9)

Spatial indices are raised and lowered with respect to $\delta_{ij}$.

The off-diagonal part in Eq. (B.8) vanishes because of the $\Phi$ equation
which, in the presence of the homogeneous spatial component, is

$\Phi=2\frac{\dot{\xi}}{M^{2}_{\textrm{Pl}}}\left(W^{2}-\frac{1}{3}a^{-2}A_{j}A^{j}\right).$
(B.10)

The equation of $A_{i}$, which can be obtained from Eq. (2.21), reads

$\Upsilon\left[\xi\ddot{A_{i}}+\left(\dot{\xi}+\xi
H\right)\dot{A_{i}}\right]=\left(\frac{3}{4}M^{2}_{\mathrm{Pl}}+\dot{\xi}W\right)A_{i}~{}.$
(B.11)

The $M^{2}_{\textrm{Pl}}G_{ij}=T^{\textrm{eff}}_{ij}$ metric field equations,
$T^{\textrm{eff}}_{\mu\nu}$ being the effective energy-momentum tensor, can be
read off from Eqs. (B.3) and (B.8) (including the matter fields)

$\displaystyle-M^{2}_{\textrm{Pl}}\left(2\dot{H}+3H^{2}\right)g_{ij}=-\left(4\ddot{\xi}H^{2}+8\dot{\xi}H^{3}+8\dot{\xi}H\dot{H}-P_{\phi}-P_{\textrm{M}}\right)g_{ij}+4\left(\frac{3}{2}M^{2}_{\textrm{Pl}}+\ddot{\xi}+3H\dot{\xi}-\right.$
$\displaystyle\left.-2\Phi\dot{\xi}\right)A_{i}A_{j}+8\dot{\xi}\dot{A}_{\left(i\right.}A_{\left.j\right)}+8\xi\Upsilon\dot{A_{i}}\dot{A_{j}}+8\left(\ddot{\xi}H\Phi-\frac{1}{2}\ddot{\xi}\Phi^{2}+\frac{3}{8}M^{2}_{\textrm{Pl}}\Phi^{2}-\frac{3}{8}M^{2}_{\textrm{Pl}}a^{-2}A_{k}A^{k}+\dot{\xi}H\dot{\Phi}+\right.$
$\displaystyle\left.+\dot{\xi}\dot{H}\Phi-\dot{\xi}\Phi\dot{\Phi}-\frac{1}{2}\dot{\xi}\Phi^{3}+\frac{3}{2}\dot{\xi}H^{2}\Phi-\frac{1}{2}\dot{\xi}Wa^{-2}A_{k}A^{k}-\frac{1}{2}\xi\Upsilon
a^{-2}\dot{A_{k}}\dot{A^{k}}\right)g_{ij}~{}.$ (B.12)

Assuming that the homogeneous spatial vector lies along the $z$ direction
without the loss of generality, $\mathbf{A}(t)=(0,0,A_{z}(t))$, we see that
$T^{\textrm{eff}}_{xx}=T^{\textrm{eff}}_{yy}\neq T^{\textrm{eff}}_{zz}$, which
means spatial anisotropy (and the off-diagonal components vanish).

Finally, $M^{2}_{\textrm{Pl}}G_{00}=T^{\textrm{eff}}_{00}$ reads

$3H^{2}M^{2}_{\textrm{Pl}}=\rho_{\phi}+\rho_{\textrm{M}}+3\left(\Phi^{2}+a^{-2}A_{i}A^{i}\right)M^{2}_{\textrm{Pl}}-4\xi\Upsilon
a^{-2}\dot{A_{i}}\dot{A^{i}}+12\dot{\xi}W\left(W^{2}-\frac{1}{3}a^{-2}A_{i}A^{i}\right).$
(B.13)

##  Appendix C Projective Transformations and Torsion

In this section of the Appendix, the homogeneous equations (2.36)-(2.38) and
the tensor perturbation equation (2.43) are shown to be valid for a connection
that is metric compatible but has non-zero torsion. For this purpose, we
review the so-called ‘projective transformations’, which are given by (see
Ref. [77])

$\displaystyle\Gamma^{\alpha}_{\mu\nu}\rightarrow\Gamma^{\alpha}_{\mu\nu}+B_{\mu}\delta_{\nu}^{\alpha}~{},$
(C.1)

where $B_{\mu}$ is an arbitrary vector. This amounts to

$\displaystyle\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}\rightarrow\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}+B_{\mu}\delta_{\nu}^{\alpha}~{}.$
(C.2)

The $\kappa$ tensors defined in App. A transform as

$\displaystyle\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}\rightarrow\tensor{\kappa}{{}^{\alpha}_{\beta\mu\nu}}+\delta_{\beta}^{\alpha}\mathring{F}_{\mu\nu}~{},$
(C.3)
$\displaystyle\kappa_{\mu\nu}\rightarrow\kappa_{\mu\nu}+\mathring{F}_{\mu\nu}~{},$
(C.4)
$\displaystyle\tilde{\kappa}_{\mu\nu}\rightarrow\tilde{\kappa}_{\mu\nu}-\mathring{F}_{\mu\nu}~{},$
(C.5)
$\displaystyle\bar{\kappa}_{\mu\nu}\rightarrow\bar{\kappa}_{\mu\nu}+4\mathring{F}_{\mu\nu}~{},$
(C.6)

where
$\mathring{F}_{\mu\nu}\equiv\mathring{\nabla}_{\mu}B_{\nu}-\mathring{\nabla}_{\nu}B_{\mu}$.
The scalar $\kappa$ is invariant under these transformations because the
tensor $\mathring{F}_{\mu\nu}$ is antisymmetric. Notice also that the tensor
$\bar{\kappa}_{\alpha\beta\mu\nu}$, which was defined as the
antisymmetrisation of the first two indices of $\kappa_{\alpha\beta\mu\nu}$
below Eq. (A.19), is invariant as well.

Now, let us apply this transformation to the Weyl connection. The
$\mathcal{G}$ in Eq. (2.16) becomes

$\displaystyle\mathcal{G}\rightarrow\mathcal{G}-4\mathring{F}^{\mu\nu}\mathring{F}_{\mu\nu}-4\Upsilon\left(\frac{1}{2}\bar{\kappa}^{\mu\nu}+\mathring{F}^{\mu\nu}\right)\mathring{F}_{\mu\nu}~{}.$
(C.7)

Under this transformation, the distortion tensor of the Weyl connection is
given by

$\displaystyle\tensor{\kappa}{{}_{\mu\nu}^{\alpha}}\rightarrow\delta_{\mu}^{\alpha}A_{\nu}+\delta_{\nu}^{\alpha}\left(A_{\mu}+B_{\mu}\right)-g_{\mu\nu}A^{\alpha}.$
(C.8)

If $B_{\mu}=-A_{\mu}$, we obtain the connection compatible with the metric but
with non-zero torsion which was the second of the three cases extracted from
Eq. (A.10) with $(c_{1},c_{2},c_{3})=(1,0,1)$. Furthermore, if the projective
transformation is applied to the LC connection, the resulting connection
corresponds to the third example, which had non-vanishing torsion and non-
metricity tensors, $(c_{1},c_{2},c_{3})=(0,1,0)$. For $B_{\mu}=-A_{\mu}$, we
have

$\displaystyle\mathring{F}_{\mu\nu}=-\frac{1}{4}\bar{\kappa}_{\mu\nu}~{},$
(C.9)

hence the $\mathcal{G}$ with torsion and vanishing non-metricity is

$\displaystyle\mathcal{G}=\mathcal{G}^{\textrm{Weyl}}-\frac{1-\Upsilon}{4}\bar{\kappa}^{\mu\nu}\bar{\kappa}_{\mu\nu}~{}.$
(C.10)

We see that $\mathcal{G}$ and $\mathcal{G}^{\textrm{Weyl}}$ only differ by the
term $\bar{\kappa}^{\mu\nu}\bar{\kappa}_{\mu\nu}$. Actually, we can see that
the new $\mathcal{G}$ does not depend on $\Upsilon$

$\displaystyle\mathcal{G}=\mathring{\mathcal{G}}+2\mathring{R}\kappa-8\mathring{R}^{\mu\nu}\tilde{\kappa}_{\mu\nu}+2\mathring{R}^{\alpha\beta\mu\nu}\bar{\kappa}_{\alpha\beta\mu\nu}+\kappa^{2}-4\tilde{\kappa}^{\mu\nu}\tilde{\kappa}_{\mu\nu}+\bar{\kappa}^{\alpha\beta\mu\nu}\bar{\kappa}_{\alpha\beta\mu\nu}-\frac{1}{4}\bar{\kappa}^{\mu\nu}\bar{\kappa}_{\mu\nu}~{}.$
(C.11)

So the calculations are the same as in the case of the Weyl connection if we
simply set $\Upsilon=1$ (see Eq. (2.16)).

If $A_{i}(t)=0$, the $\Upsilon$ term does not contribute to the homogeneous
and tensor perturbations equations as it multiplies $\bar{\kappa}_{\mu\nu}$
(which vanishes in that case; see Eq. (2.21)). Consequently, for a background
Weyl vector that has solely a non-zero temporal component (so it does not
break spatial isotropy), the results derived in this work are the same as
those obtained assuming the $(c_{1},c_{2},c_{3})=(1,0,1)$ connection with non-
zero torsion and metric compatible.

## References

  * [1] A. G. Riess et al. [Supernova Search Team], _Observational evidence from supernovae for an accelerating universe and a cosmological constant_ , _Astron. J._ 116 1009-1038 (1998) [astro-ph/9805201].
  * [2] S. Perlmutter et al. [Supernova Cosmology Project], _Measurements of $\Omega$ and $\Lambda$ from $42$ high redshift supernovae_, _Astrophys. J._ 517 565-586 (1999) [astro-ph/9812133].
  * [3] N. Aghanim et al. [Planck Collaboration], _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 A6 (2020) [Erratum: Astron. Astrophys. 652, C4 (2021)] [1807.06209].
  * [4] S. M. Carroll, _The Cosmological constant_ , _Living Rev. Rel._ 4, 1 (2001) [astro-ph/0004075].
  * [5] J. Martin, _Everything You Always Wanted To Know About The Cosmological Constant Problem (But Were Afraid To Ask)_ , _Comptes Rendus Physique_ 13, 566-665 (2012) [1205.3365].
  * [6] A. Joyce, B. Jain, J. Khoury and M. Trodden, _Beyond the Cosmological Standard Model_ , _Phys. Rept._ 568, 1-98 (2015) [1407.0059].
  * [7] S. Weinberg, _The Cosmological Constant Problem_ , _Rev. Mod. Phys._ 61, 1-23 (1989).
  * [8] S. Weinberg, _The Cosmological constant problems_ , _Contribution to 4th International Symposium on Sources and Detection of Dark Matter in the Universe_ (DM 2000) (2000) [astro-ph/0005265].
  * [9] R. R. Caldwell, R. Dave and P. J. Steinhardt, _Cosmological imprint of an energy component with general equation of state_ , _Phys. Rev. Lett._ 80, 1582-1585 (1998) [astro-ph/9708069].
  * [10] I. Zlatev, Li-Min Wang and P. J. Steinhardt, _Quintessence, cosmic coincidence, and the cosmological constant_ , _Phys. Rev. Lett._ 82, 896-899 (1999) [astro-ph/9807002].
  * [11] E. J. Copeland, M. Sami and S. Tsujikawa, _Dynamics of dark energy_ , _Int. J. Mod. Phys. D_ 15, 1753-1936 (2006) [hep-th/0603057].
  * [12] S. Tsujikawa, _Dark energy: investigation and modeling_ , _ASSL_ Vol. 370, 331-402, Springer, Dordrecht (2011) [1004.1493].
  * [13] S. Tsujikawa, _Quintessence: A Review_ , _Class. Quant. Grav._ 30, 214003 (2013) [1304.1961].
  * [14] P. J. E. Peebles and A. Vilenkin, _Quintessential inflation_ , _Phys. Rev. D_ 59, 063505 (1999) [astro-ph/9810509].
  * [15] K. Dimopoulos and J. W. F. Valle, _Modeling Quintessential Inflation_ , _Astropart. Phys._ 18 287-306 (2002) [astro-ph/0111417].
  * [16] J. de Haro, _On the viability of quintessential inflation models from observational data_ , _Gen. Rel. Grav._ 49, 1, 6 (2017) [1602.07138].
  * [17] K. Dimopoulos, _Jointly modelling Cosmic Inflation and Dark Energy_ , _J. Phys. Conf. Ser._ 2105, 1, 012001 (2021) [2106.14966].
  * [18] A. A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91, 99-102 (1980).
  * [19] A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 347-356 (1981).
  * [20] A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 389-393 (1982).
  * [21] S. Capozziello, _Curvature quintessence_ , _Int. J. Mod. Phys. D_ 11, 483-492 (2002) [gr-qc/0201033].
  * [22] S. M. Carroll, V. Duvvuri, M. Trodden and M. S. Turner, _Is cosmic speed - up due to new gravitational physics?_ , _Phys. Rev. D_ 70, 043528 (2004) [astro-ph/0306438].
  * [23] T. Koivisto and D. F. Mota, _Cosmology and Astrophysical Constraints of Gauss-Bonnet Dark Energy_ , _Phys. Lett. B_ 644, 104-108 (2007) [astro-ph/0606078].
  * [24] S. Nojiri, S. D. Odintsov and V. K. Oikonomou, _Modified Gravity Theories on a Nutshell: Inflation, Bounce and Late-time Evolution_ , _Phys. Rept._ 692, 1-104 (2017) [1705.11098].
  * [25] C. van de Bruck, K. Dimopoulos, C. Longden and C. Owen, _Gauss-Bonnet-coupled Quintessential Inflation_ , 1707.06839.
  * [26] J. de Haro and L. A. Saló, _A Review of Quintessential Inflation_ , _Galaxies_ 9, 4, 73 (2021)[2108.11144].
  * [27] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116, 6, 061102 (2016) [1602.03837].
  * [28] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _GW151226: Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence_ , _Phys. Rev. Lett._ 116, 24, 241103 (2016) [1606.04855].
  * [29] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119, 16, 161101 (2017) [1710.05832].
  * [30] J. Sakstein and B. Jain, _Implications of the Neutron Star Merger GW170817 for Cosmological Scalar-Tensor Theories_ , _Phys. Rev. Lett._ 119, 25, 251303 (2017) [1710.05893].
  * [31] J. M. Ezquiaga and M. Zumalacárregui, _Dark Energy After GW170817: Dead Ends and the Road Ahead_ , _Phys. Rev. Lett._ 119, 25, 251304 (2017) [1710.05901].
  * [32] T. Baker, E. Bellini, P. G. Ferreira, M. Lagos, J. Noller and I. Sawicki, _Strong Constraints on Cosmological Gravity from GW170817 and GRB 170817A_ , _Phys. Rev. Lett._ 119, 25, 251301 (2017) [1710.06394].
  * [33] S. Arai and A. Nishizawa, _Generalized framework for testing gravity with gravitational-wave propagation. II. Constraints on Horndeski theory_ , _Phys. Rev. D_ 97, 10, 104038 (2018) [1711.03776].
  * [34] P. Creminelli and F. Vernizzi, _Dark Energy after GW170817 and GRB170817A_ , _Phys. Rev. Lett._ 119, 25, 251302 (2017) [1710.05877].
  * [35] D. Langlois, R. Saito, D. Yamauchi and K. Noui, _Scalar-tensor theories and modified gravity in the wake of GW170817_ , _Phys. Rev. D_ 97, 6, 061501 (2018) [1711.07403].
  * [36] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _Multi-messenger Observations of a Binary Neutron Star Merger_ , _Astrophys. J. Lett._ 848, 2, L12 (2017) [1710.05833].
  * [37] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _Gravitational Waves and Gamma-rays from a Binary Neutron Star Merger: GW170817 and GRB 170817A_ , _Astrophys. J. Lett._ 848, 2, L13 (2017) [1710.05834].
  * [38] G. W. Horndeski, _Second-order scalar-tensor field equations in a four-dimensional space_ , _Int. J. Theor. Phys._ 10, 363-384 (1974).
  * [39] C. de Rham and S. Melville, _Gravitational Rainbows: LIGO and Dark Energy at its Cutoff_ , _Phys. Rev. Lett._ 121, 22, 221101 (2018) [1806.09417].
  * [40] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _Generalized G-inflation: Inflation with the most general second-order field equations_ , _Prog. Theor. Phys._ 126, 511-529 (2011) [1105.5723].
  * [41] P. G. S. Fernandes, P. Carrilho, T. Clifton and D. J. Mulryne, _The 4D Einstein-Gauss-Bonnet theory of gravity: a review_ , _Class. Quant. Grav._ 39, 6, 063001 (2022) [2202.13908].
  * [42] S. D. Odintsov and V. K. Oikonomou, _Inflationary Phenomenology of Einstein Gauss-Bonnet Gravity Compatible with GW170817_ , _Phys. Lett. B_ 797, 134874 (2019) [1908.07555].
  * [43] S. D. Odintsov, V. K. Oikonomou and F. P. Fronimos, _Rectifying Einstein-Gauss-Bonnet Inflation in View of GW170817_ , _Nucl. Phys. B_ 958, 115135 (2020) [2003.13724].
  * [44] S. D. Odintsov and V. K. Oikonomou, _Swampland implications of GW170817-compatible Einstein-Gauss-Bonnet gravity_ , _Phys. Lett. B_ 805, 135437 (2020) [2004.00479].
  * [45] J. Hjorth et al., _The Distance to NGC 4993: The Host Galaxy of the Gravitational-wave Event GW170817_ , _Astrophys. J. Lett._ 848, 2, L31 (2017) [1710.05856].
  * [46] J. J. Terente Díaz, K. Dimopoulos, M. Karčiauskas and A. Racioppi, _Gauss-Bonnet Dark Energy and the speed of gravitational waves_ , _JCAP_ 10, 031 (2023) [2307.06163].
  * [47] L. N. Granda and D. F. Jimenez, _The speed of gravitational waves and power-law solutions in a scalar-tensor model_ , _Astropart. Phys._ 103, 115-121 (2018) [1802.03781].
  * [48] Y. Gong, E. Papantonopoulos and Z. Yi, _Constraints on scalar-tensor theory of gravity by the recent observational results on gravitational waves_ , _Eur. Phys. J. C_ 78, 9, 738 (2018) [1711.04102].
  * [49] H. Mohseni Sadjadi, _Scalar-Gauss-Bonnet model, the coincidence problem and the gravitational wave speed_ , 2309.07816.
  * [50] H. Mohseni Sadjadi, _Non-minimally coupled quintessence in the Gauss-Bonnet model, symmetry breaking, and cosmic acceleration_ , 2310.03697.
  * [51] T. Helpin and M. S. Volkov, _Varying the Horndeski Lagrangian within the Palatini approach_ , _JCAP_ 01, 044 (2020) [1906.07607].
  * [52] T. Helpin and M. S. Volkov, _A metric-affine version of the Horndeski theory_ , _Int. J. Mod. Phys. A_ 35, 02n03, 2040010 (2020) [1911.12768].
  * [53] M. Kubota, Kin-Ya Oda, K. Shimada and M. Yamaguchi, _Cosmological Perturbations in Palatini Formalism_ , _JCAP_ 03, 006 (2021) [2010.07867].
  * [54] Yu-Qi. Dong and Yu-Xiao Liu, _Polarization modes of gravitational waves in Palatini-Horndeski theory_ , _Phys. Rev. D_ 105, 6, 6 (2022) [2111.07352].
  * [55] Yu-Qi Dong, Yu-Qiang Liu and Yu-Xiao Liu, _Constraining Palatini-Horndeski theory with gravitational waves after GW170817_ , _Eur. Phys. J. C_ 83, 8, 702 (2023) [2211.12056].
  * [56] K. Dimopoulos and S. Sánchez López, _Quintessential inflation in Palatini $f(R)$ gravity_, _Phys. Rev. D_ 103, 4, 043533 (2021) [2012.06831].
  * [57] K. Dimopoulos, A. Karam, S. Sánchez López and E. Tomberg, _Palatini $R^{2}$ quintessential inflation_, _JCAP._ 10, 076 (2022) [2206.14117].
  * [58] M. Kumar Sharma and S. Sur, _Growth of matter perturbations in an interacting dark energy scenario emerging from metric-scalar-torsion couplings_ , 2102.01525.
  * [59] J. B. Jimenez and T. S. Koivisto, _Extended Gauss-Bonnet gravities in Weyl geometry_ , _Class. Quant. Grav._ 31, 135002 (2014) [1402.1846].
  * [60] E. J. Copeland, A. R. Liddle and D. Wands, _Exponential potentials and cosmological scaling solutions_ , _Phys. Rev. D_ 57, 4686-4690 (1998) [gr-qc/9711068].
  * [61] P. G. Ferreira and M. Joyce, _Cosmology with a primordial scaling field_ , _Phys. Rev. D_ 58, 023503 (1998) [astro-ph/9711102].
  * [62] T. Barreiro, E. J. Copeland and N. J. Nunes, _Quintessence arising from exponential potentials_ , _Phys. Rev. D_ 61, 127301 (2000) [astro-ph/9910214].
  * [63] S. M. Carroll, _Lecture Notes On General Relativity_ , gr-qc/9712019.
  * [64] S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud and E. Di Valentino, _Teleparallel gravity: from theory to cosmology_ , _Rept. Prog. Phys._ 86, 2, 026901 (2023) [2106.13793].
  * [65] M. Borunda, B. Janssen and M. Bastero-Gil, _Palatini versus metric formulation in higher curvature gravity_ , _JCAP_ 11, 008 (2008) [0804.4440].
  * [66] T. Clifton, P. G. Ferreira, A. Padilla and C. Skordis, _Modified Gravity and Cosmology_ , _Phys. Rept._ 513, 1-189 (2012) [1106.2476].
  * [67] D. Baumann, _Inflation_ , _Physics of the Large and the Small_ 523-686 (2011) [0907.5424].
  * [68] T. Koivisto and D. F. Mota, _Gauss-Bonnet Quintessence: Background Evolution, Large Scale Structure and Cosmological Constraints_ , _Phys. Rev. D_ 75, 023518 (2006) [hep-th/0609155].
  * [69] D. H. Weinberg, M. J. Mortonson, D. J. Eisenstein, C. Hirata, A. G. Riess and E. Rozo, _Observational Probes of Cosmic Acceleration_ , _Phys. Rept._ 530, 87-255 (2013) [1201.2434].
  * [70] L. Jarv, P. Kuusk and M. Saal, _Remarks on (super-)accelerating cosmological models in general scalar-tensor theory_ , _Proc. Est. Acad. Sci. Phys. Math._ 59, 306-312 (2010) [0903.4357].
  * [71] H. Mohseni Sadjadi, _On cosmic acceleration in four dimensional Einstein-Gauss-Bonnet gravity_ , _Phys. Dark Univ._ 30, 100728 (2020) [2005.10024].
  * [72] B. P. Abbott et al. [LIGO Scientific, Virgo et al. Collaboration], _Properties of the binary neutron star merger GW170817_ , _Phys. Rev. X_ 9, 1, 011001 (2019) [1805.11579].
  * [73] S. Nojiri, S. D. Odintsov and M. Sasaki, _Gauss-Bonnet dark energy_ , _Phys. Rev. D_ 71, 123509 (2005) [hep-th/0504052].
  * [74] K. Dimopoulos, _Can a vector field be responsible for the curvature perturbation in the Universe?_ , _Phys. Rev. D_ 74, 083502 (2006) [hep-ph/0607229].
  * [75] N. Dadhich and J. M. Pons, _On the equivalence of the Einstein-Hilbert and the Einstein-Palatini formulations of general relativity for an arbitrary connection_ , _Gen. Rel. Grav._ 44, 2337-2352 (2012) [1010.0869].
  * [76] A. N. Bernal, B. Janssen, A. Jimenez-Cano, J. A. Orejuela, M. Sanchez and P. Sanchez-Moreno, _On the (non-)uniqueness of the Levi-Civita solution in the Einstein-Hilbert-Palatini formalism_ , _Phys. Lett. B_ 768, 280-287 (2017) [1606.08756].
  * [77] J. A. Schouten, _Ricci-Calculus: An Introduction to Tensor Analysis and Its Geometrical Applications_ , _Springer-Verlag_ , Berlin and Heidelberg (1954).