# Improving data-driven model-independent reconstructions and updated
constraints on dark energy models from Horndeski cosmology

Mauricio Reyes  and Celia Escamilla-Rivera  0000-0002-8929-250X

###### Abstract

In light of the statistical performance of cosmological observations, in this
work we present an improvement on the Gaussian reconstruction of the Hubble
parameter data $H(z)$ from Cosmic Chronometers, Supernovae Type Ia and
Clustering Galaxies in a model-independent way in order to use them to study
new constraints in the Horndeski theory of gravity. First, we have found that
the prior used to calibrate the Pantheon supernovae data significantly affects
the reconstructions, leading to a 13$\sigma$ tension on the $H_{0}$ value.
Second, according to the $\chi^{2}$-statistics, the reconstruction carried out
by the Pantheon data calibrated using the $H_{0}$ value measured by The
Carnegie-Chicago Hubble Program is the reconstruction which fits best the
observations of Cosmic Chronometers and Clustering of Galaxies datasets.
Finally, we use our reconstructions of $H(z)$ to impose model-independent
constraints in dark energy scenarios as Quintessence and K-essence from
general cosmological viable Horndeski models, landscape in where we found that
a Horndeski model of the K-essence type can reproduce the reconstructions of
the late expansion of the universe within 2$\sigma$.

##  1 Introduction

Currently, local measurements of the Hubble-Lemaître parameter from supernovae
type Ia [1, 2] and lensing time delays [3] are in disagreement with the value
inferred from a $\Lambda$CDM fit to the Cosmic Microwave Background (CMB) [4].
This discrepancy has not been easily explained by any obvious systematic
effect [5, 6] in either measurement. According to this, an increasing
discussion on the cosmology community [7, 8] is focusing on the possibility
that this Hubble-Lemaître tension may be indicating new physics beyond the
concordance model [9].

However, theoretical proposals for the Hubble-Lemaître tension are not easily
within reach, since one of the biggest challenge remains in the calculations
of the angular scale of the acoustic peaks in the CMB power spectrum, which
fix the ratio of the sound horizon at decoupling to the distance to the CMB
surface of last scatter. Furthermore, some solutions involve early-time
modifications on the theory of gravity itself, which change the sound horizon,
or even late-time changes to the expansion rate. On the side of late-time
proposals, we can mention: a phantom-like dark energy (DE) component [10], a
vacuum phase transition [11], or interacting DE [12]. Moreover, these models
are tightly constrained [13, 14] by late-time observables, especially those
from baryon acoustic oscillations (BAO) [15]. Model-independent
parameterizations of the late-time expansion history are similarly constrained
[16], while others proposals aim to put theoretical background to
parameterizations modified gravity at the same level as other
parameterizations into the pipeline and analysis of observational data and
forecasts [17]. On the early-time resolution, it is possible to reduce the
sound horizon with radiation energy density, which can be constrained by BAO
and by the CMB power spectrum [14]. As an extension, it is also possible to
address the Hubble tension through a direct modification of gravity [14, 8]
that is active around late and early times along the Hubble flow, e.g. in
Galileon models [18, 19, 20, 21] and also in extensions as Teleparallel
Gravity [22, 23, 24]. Furthermore, in [25] was shown a method to reconstruct
Horndeski cosmological models using the Effective Field Theory approach [26]
and interpolation methods with a piece-wise fifth order spline.

The discussed Hubble tension is one of the several issues111Despite great
efforts, dark matter remains undetected and the nature of the Cosmological
Constant via dark energy continues to have several issued associated with it.
why it has recently become of interest the study of modified theories of
General Relativity (GR). It has been of particular interest the study of
Hordesnki theory [27, 28, 29], this is the most general scalar-tensor theory
in four dimensions that is built only from the metric tensor $g_{\mu\nu}$ and
a scalar field $\phi$ that is able of avoiding Ostrogradsky instabilities
[30]. The theory has four free functions $G_{2},G_{3},G_{4},G_{5}$ that depend
on the scalar field $\phi$ and its kinetic term
$2X=-\partial_{\mu}\phi\partial^{\mu}\phi$. The correct choice of the $G_{i}$
functions allows us to recover various theories including models such as
Quintessence, Brans–Dicke, k-essence and others [31]. At the present time
observations [30] have reduced the theory to $G_{5}=0$ and $G_{4}$ to be a
function that just depends on the scalar field $G_{4}(\phi)$. Given the
generality of this theory, if the tension $H_{0}$ can be solved by resorting
to an extension to GR, Horndeski theory will be able to describe that
extension, at least effectively with a proper choice of $G_{i}$. In this work
it will be of particular interest to find a combination of $G_{i}$’s that
could be a candidate to solve the $H_{0}$ tension [32].

Furthermore, when consider data from a specific survey we have to deal with a
problem of model-dependency according to the nature of the observations, e.g.
observations from CMB are biased by a $\Lambda$CDM model. This problem has
motivated the study of other data treatment methods that extract important
information directly from the observational data, as for example promising
alternatives related to non-parametric reconstructions of the cosmic expansion
[33, 34, 35, 36]. These kinds of approaches attempt to reconstruct the
cosmological evolution directly from observations without establishing an
association with a theoretical physical model. In this line of thought, we can
divide these approaches as model-dependent and model-independent methods.
Model-dependent methods are useful when the relationship between the variables
of the landscape under study is known, and their objective is to constrain the
parameters of the chosen cosmological model. Moreover, when there is not any
hint about the explicit form of this relationship we need to propose it at
hand, which can leads to biased results. On the other hand, model-independent
methods provide a generic trend of the parameter of interest when the
relationship between the parameters is unknown or there is little feedback
about it. These kinds of methods have become very recurrent given their
characteristics for enhancing scatter and other diagnostic analyses whose can
display the underlying structure in the data. Some example in that regard is
the Gaussian Process (GP) [37], which is a Bayesian statistical approach that
allows to obtain a function that describes a specific data set without an
assuming a particular cosmological model. But of course, the price we pay to
have a model-independent method is the need of a prior for the function that
we seek to reconstruct and a likelihood for the observational data. Once this
information is given we can calculate the posterior distribution for the
functions that can describe the observations.

According to the ideas discussed above, we will study the late-time expansion
of the universe in a model-independent way through the reconstruction of the
late observational data. To perform this analysis, we will make use of
observational data that involve an improved GP reconstruction version of
$H(z)$. Once we have at hand the reconstructed late expansion of the universe,
we proceed to use our new reconstructions to found new cosmological
constraints in Hordenski models, in particular we will study a Quintessence
and a K-essence model and develop a possible solution of $H_{0}$ tension at
late-times. In this line of research, interesting analyses have been done in
[32] where it was found that early modified gravity techniques can reconcile
the $H_{0}$ value by increasing the expansion rate in the epoch of matter-
radiation equality. Also extensions up to cubic covariant Galileon model has
been study in [38]. Furthermore, surviving Horndeski effective field theory of
dark energy has been treated at all redshift with the Gravitational Waves
constraints [39].

We outlined this paper as follows: in Sec. 2 we explain our new GP methodology
on the observational datasets to perform model-independent reconstructions of
$H(z)$. Also, we include the description of GP kernels used to perform such
reconstructions. In Sec. 3 we consider a model-independent methodology to
constraint dark energy density at late-times showing that it is possible to
report a 1$\sigma$ statistical uncertainty with a numerical error associated
with the calculation of this density. In Sec. 4 we employ these model-
independent reconstructions to derive new constraints on specific relations of
the Galileon $G_{2}$ from Horndeski theory of gravity. We found that under
certain initial conditions our system can be reduced to a one-parameter to be
fitted with the observations. Finally, we present our results and conclusions
in Sec. 5.

##  2 Improving GP reconstructions for the late-time data

Gaussian Processes (GP) employed in cosmology provide a method to extract
information from observations without assuming a particular model for the
description of the data, which makes them a powerful tool, since they can help
us to understand the nature of the late time expansion of the universe [37,
40]. From a statistical point of view, GP is a stochastic process formed by a
collection of random variables $\boldsymbol{x}$, such that every finite linear
combination of these variables has a multivariate Gaussian distribution

$\boldsymbol{x}\sim\mathcal{N}(\mu,\Sigma),$ (2.1)

with $\mu$ being the vector that contains the mean value of the random
variables and $\Sigma$ the covariance matrix between the random variables. A
Gaussian process for a function $f(z)$ is described by a covariance function
(kernel) $k(\widetilde{z},z)$ and by its mean value. The covariance function
expresses the fact that the function evaluated at one point $z$ is not
independent of the function evaluated at other point $\widetilde{z}$, but
rather these are correlated and the relationship is given by the covariance
function. There is a wide range of possible covariance functions222 As long as
a reasonable covariance function is chosen for the problem at hand. In our
case we need the covariance function to be continuous and differentiable at
least one time. Also, we expect that the correlation between the function
evaluated at two points decreases with the distance between them., moreover it
has been shown that the choice of a kernel does not significantly influence
the derived value of the cosmological parameters to be determined by the
reconstruction [35, 36, 34], i.e. even when the mean value of these parameters
varies when we change the covariance function, these values coincide within
1$\sigma$. This motivates us to consider the following two different kernels
for our analysis:

$\displaystyle k(z,z^{\prime})$ $\displaystyle=$
$\displaystyle\sigma_{f}^{2}\exp\left(\frac{-(z-z^{\prime})^{2}}{2l^{2}}\right),\quad\rightarrow\quad\text{Squared
Exponential}$ (2.2) $\displaystyle k(z,z^{\prime})$ $\displaystyle=$
$\displaystyle\sigma_{f}^{2}\left(\frac{l_{f}}{(z-z^{\prime})^{2}+l_{f}^{2}}\right),\quad\rightarrow\quad\text{Cauchy
function}$ (2.3)

where $\sigma_{f}$ determines the width of the reconstructed function and $l$
determines how the correlation between two points change, a large $l$ results
in a smooth reconstructed function, while a small $l$ gives an oscillating
function [37].

GP have been widely used in the literature for different purposes, ranging
from the study of growth rate of structure $f\sigma_{8}(z)$ [41], to obtain
constraints on Teleparallel gravity models [35], and even some authors have
found a way to use the reconstructions to test the $\Lambda$CDM model, in [37]
they use the reconstructions of $H(z)$ to look for deviations from the
equation of state of the cosmological constant $\Lambda$, in [42] they found a
way to test the concordance model using a diagnostic function $O_{m}$ [42,
43]. To proceed our reconstructions we require the following samples:

  1. 1.

Pantheon SNeIa sample. Consists of 1048 SNeIa in 40 bins compressed.
Currently, this sample is the largest spectroscopically supernovae Ia data
set. These astrophysical objects can give determinations of the distance
modulus $\mu$, its theoretical prediction is related to the luminosity
distance $d_{L}$

$\mu(z)=5\log{\left[\frac{d_{L}(z)}{1\text{Mpc}}\right]}+25,$ (2.4)

where the luminosity distance is given in Mpc. In this expression we need to
add the nuisance parameter $M$, as an unknown offset sum of the supernovae
absolute magnitude (with other possible systematics), which can be degenerate
with the $H_{0}$ value. We assume a spatial flatness and $d_{L}$ can be
related to the comoving distance $D$ using
$d_{L}(z)=\frac{c}{H_{0}}(1+z)D(z),$ where $c$ is the speed of light, and
obtain

$D(z)=\frac{H_{0}}{c}(1+z)^{-1}10^{\frac{\mu(z)}{5}-5}.$ (2.5)

The normalised Hubble function $H(z)/H_{0}$ is derived from the inverse of the
derivative of $D(z)$ with respect to $z$
$D(z)=\int^{z}_{0}H_{0}d\tilde{z}/H(\tilde{z}),$ where $H_{0}$ is the Hubble
value consider as a prior to normalise $D(z)$. For our sample, we calibrated
this data by using four different priors for $H_{0}$:

     * •

Cefeids [1] : $H_{0}^{\text{R}}=73.24\pm 1.74\text{km/s/Mpc}$.

     * •

Planck 2018 [4]: $H_{0}^{\text{Planck}}=67.4\pm 0.5\text{km/s/Mpc}$.

     * •

H0LiCOW [44]: $H_{0}^{\text{HW}}=73.3^{+1.7}_{-1.8}\text{km/s/Mpc}$.

     * •

TRGB [45]: $H_{0}^{\text{TRGB}}=69.8\pm 1.9\text{km/s/Mpc}$.

  2. 2.

Cosmic Chronometers sample (CC). This data provide model independent
measurements for H(z). These measurements are obtained by comparing the
stellar evolution of galaxies with low star formation and similar metallicity
[46]. Currently there are 31 model-independent measurements [47] of $H(z)$
obtained through this method.

  3. 3.

Clustering Galaxies sample (CL). Usually, this sample analyse the baryon
acoustic oscillations (BAO) with a two-point correlation function, the
analysis allows them to get measurements [48] of the product $H(z)r_{d}$,
where $r_{d}$ is the sound horizon at the drag epoch. If we assume a fiducial
cosmology they can estimate $r_{d}$, and then break the degeneracy between
$r_{d}$ and $H(z)$. Currently, there are 20 measurements of $H(z)$ obtained
with this method [47]. However, we must be careful working with this data,
since is biased due to an underlying $\Lambda$CDM cosmology.

To perform the reconstructions of $H(z)$, it is necessary to transform the
measurements. We start by taking the data of the Pantheon sample, once is
calibrated it gives us information about the distance module $\mu(z)$, which
is related to the luminous distance through the following equation

$d_{L}=10^{(\mu(z)-25)/5},$ (2.6)

we substitute the data of $\mu(z)$ in (2.6), and we obtain

$d_{L}=\overline{d}_{L}+\sigma_{d_{L}},$ (2.7)

where $\overline{d}_{L}$ is the mean value of $d_{L}$ and $\sigma_{d_{L}}$ is
the uncertainty that results from the error propagation. In order to simplify
the analysis, it is useful to define the following quantity

$d_{p}\equiv\frac{d_{L}}{(1+z)c},$ (2.8)

where $c$ is the speed of light. When we substitute (2.7) in (2.8) allows us
to transform the data from $\mu(z)$ into information about the function
$d_{p}(z)$. Moreover, if we assume that the universe is spatially flat and we
make use of (2.8) we can show that

$d_{p}^{\prime}=\frac{1}{H(z)},$ (2.9)

which allows us to transform the measurements of CC and CL to information
about $d_{p}^{\prime}(z)$. To do so, we substitute the measurements of $H(z)$
in (2.9) to obtain

$d^{\prime}_{p}(z)=\frac{1}{\overline{H}(z)}\pm\frac{\sigma_{H}(z)}{\overline{H}^{2}(z)},$
(2.10)

where $\overline{H}(z)$ is the mean value of the observations of $H(z)$ with
redshift $z$ and $\sigma_{H}(z)$ is its corresponding uncertainty at
1$\sigma$.

Once the process of transforming the measurements is performed we can notice
that Pantheon data gives information about the function $d_{p}(z)$ and the CC
and CL measurements give information about $d^{\prime}p(z)$. We can combine
these data sets and apply a GP. With this process we can reconstruct the
function $d_{p}$ as a function of $z$ in the redshift range of the observable,
e.g. for Pantheon compilation $0<z<2.36$, using this function we can get
$H(z)$ from

$H(z)=\frac{1}{d^{\prime}_{p}(z)}=\frac{1}{\overline{d^{\prime}}_{p}(z)}\pm\frac{\sigma_{d^{\prime}_{p}}(z)}{\overline{d^{\prime}}^{2}_{p}(z)},$
(2.11)

where $\overline{d^{\prime}_{p}}$ is the mean value of the derivative of the
reconstruction and $\sigma_{d^{\prime}_{p}}$ is the standard deviation of
$d^{\prime}_{p}$. In Figure 1 we show the results for each one of the
reconstructions using the kernels (2.2)-(2.3) for three sets of data: (1)
Pantheon+$H_{0}$ prior, (2) Pantheon+$H_{0}$ prior + CC, and (3)
Pantheon+$H_{0}$ prior + CC + CL.

Furthermore, since the reconstructions have been calibrated with different
priors for $H_{0}$, we can compute the $\chi^{2}$-statistics for each case.
Therefore, we proceed to define the following quantities:

$\displaystyle\chi^{2}_{\text{Pantheon}}=\sum_{i}^{N_{\text{Pantheon}}}\frac{\left(\mu(z_{i})_{\text{obs}}-\mu(z_{i})_{\text{recons}}\right)^{2}}{\sigma(z_{i})_{\text{obs,Pantheon}}^{2}},$
(2.12)
$\displaystyle\chi^{2}_{\text{CC}}=\sum_{i}^{N_{\text{CC}}}\frac{\left(H(z_{i})_{\text{obs}}-H(z_{i})_{\text{recons}}\right)^{2}}{\sigma(z_{i})_{\text{obs,CL}}^{2}},$
(2.13)
$\displaystyle\chi^{2}_{\text{CL}}=\sum_{i}^{N_{\text{CL}}}\frac{\left(H(z_{i})_{\text{obs}}-H(z_{i})_{\text{recons}}\right)^{2}}{\sigma(z_{i})_{\text{obs,CL}}^{2}},$
(2.14)
$\displaystyle\chi^{2}_{H_{0}}=\frac{\left(H_{0,\text{obs}}-H_{0,\text{recons}}\right)^{2}}{\sigma(z_{i})_{\text{obs},H_{0}}^{2}},$
(2.15)

where $N$ is the number of $D$ data point, $D_{\text{obs}}(z_{i})$ is the data
point observed with redshift $z_{i}$ , $\sigma_{i}$ is the uncertainty
associated with each measurement, $D_{\text{recons}}(z_{i})$ is the value
predicted by the reconstruction, and we define the following quantity

$\chi^{2}_{\text{Total}}=\chi^{2}_{\text{Pantheon}}+\chi^{2}_{\text{CC}}+\chi^{2}_{\text{CL}}+\chi^{2}_{H_{0}}.$
(2.16)

Figure 1: Reconstructions of $H(z)$ with the squared exponential kernel (2.2)
and its corresponding uncertainties up to 1$\sigma$ with the four $H_{0}$
prior. From top to bottom: Using $H^{\text{R}}_{0}$, using
$H_{0}^{\text{HW}}$, using $H_{0}^{\text{Planck}}$, and using
$H_{0}^{\text{TRGB}}$. We denote: (1) Pantheon+$H_{0}$ prior (green olive
color), (2) Pantheon+$H_{0}$ prior+ CC (red color), and (3) Pantheon+$H_{0}$
prior + CC + CL (blue color). We include also the observational data given by
the green olive color and red color points which describe CL and CC
measurements, respectively.

Figure 2:  Reconstructions of $H(z)$ with the Cauchy function kernel (2.3) and
its corresponding uncertainties up to 1$\sigma$ with the four $H_{0}$ prior.
From top to bottom: Using $H^{\text{R}}_{0}$, using $H_{0}^{\text{HW}}$, using
$H_{0}^{\text{Planck}}$, and using $H_{0}^{\text{TRGB}}$. We denote: (1)
Pantheon+$H_{0}$ prior (green olive color), (2) Pantheon+$H_{0}$ prior+ CC
(red color), and (3) Pantheon+$H_{0}$ prior + CC + CL (blue color). We include
also the observational data given by the green olive color and red color
points which describe CL and CC measurements, respectively.

In Table 1 we show the results for $\chi^{2}_{\text{Total}}$ (2.16) for each
reconstruction case. Notice that the reconstruction that better describes the
data sets is the one calibrated with $H_{0}^{\text{TRGB}}$, therefore this
reconstruction can guide us in the search for a model that could solve the
$H_{0}$ tension with the observations of the late-time universe.

$\chi^{2}_{\text{Total}}$(2.16) for the different reconstructions.
---
Data sets |  Squared exponential kernel (2.2) |  Cauchy kernel (2.3)
Pantheon+ $H_{0}^{\text{R}}$ \+ CC+ CL |  1151.88 |  1150.53
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC+ CL |  1065.38 |  1063.97
Pantheon+ $H_{0}^{\text{Planck}}$ \+ CC+ CL |  1138.33 |  1137.24
Pantheon + $H_{0}^{\text{TRGB}}$ \+ CC+ CL |  952.90 |  952.04
Table 1: $\chi^{2}$-statistics for the full data sets reconstructed Pantheon+
$H_{0}$ \+ CC+ CL (first column) using the $H_{0}$ priors for the squared
exponential kernel (second column) and the Cauchy kernel (third column), using
(2.16).

###  2.1 GP calibrated reconstruction methods and comparison tests

To compare our model-independent reconstruction with other proposals in this
direction, several methods have been presented in the literature to obtain
reconstructions of the Hubble parameter [36, 35, 34, 49]. Moreover, some of
them present some disadvantages with respect to our methodology:

  1. 1.

SNIa/QSO cosmokinetic reconstruction method [34]. In this proposal a GP
reconstruction is based to cover a large redshift range by using quasars [50,
51] calibrated measurements with Gamma Ray Bursts [52, 52, 53]. In comparison
to our proposal in the SNIa/QSO cosmokinetic reconstruction method is not
indicated the value of $H_{0}$ used to calibrate Pantheon data, therefore, a
straightforward bias emerges as a consequence of this prior.

  2. 2.

Covariance functions reconstruction method [35]. In this research the main
goal was to obtain constraints on Teleparallel Gravity models [23] through GP
using the CC data set, BAO data [54, 55, 56], and the Pantheon compressed data
set [57] including 15 SNeIa measured by CANDELS and CLASH Multy-Cycle Treasury
program (hereafter, PMCT). This data set compresses all the information of
SNeIa in only six data points for $E(z)=H(z)/H_{0}$. However, in this method
it was used only five points, since the last one is non-Gaussian distributed.
Since the PMCT data is given in terms of $E(z)$, the method requires the
calibration of this sample before being reconstructed. To calibrate the PCMT
data, first they perform a GP on the CC data, this allows them to get a value
for $H_{0}$. Using this value of $H_{0}$ they calibrate the PCMT data. Second,
it is possible to perform a GP combining CC+PMCT data to obtain a new $H_{0}$
value and afterwards they use this value to re-calibrate the PMCT sample.
Finally, the process is repeat it until the obtained value of $H_{0}$
converges. It is important to notice that this method found that $\Lambda$CDM
is within the 1$\sigma$ region for all the reconstructions, regardless of the
data sets used in their analysis.

  3. 3.

Reconstruction with weighted polynomial regression method [36]. In this
proposal was presented a method to obtain constraints on the late-cosmic
expansion. Based on CC+PMCT data samples and under the consideration to
circumvent the use of $H(z)$ measurements obtained from BAO data. It is
possible to calibrate the PMCT data and use the same procedure as in the GP
covariance functions reconstruction method. The results of this method suggest
that the reconstructions have a preference for low values of $H_{0}$.

The methods discussed in (2) and (3) have two main disadvantages with respect
to our proposal:

  * •

The number of data points used in methods (2) and (3) to perform the
reconstructions is below 50 points, which is much lower than the number of
points used in our analysis, which correspond to 31 points for the CC data
set, 20 points for CL, 1048 for Pantheon and one prior for calibration. The
low number of points used in their reconstruction results in a 1$\sigma$
uncertainty for their reconstructions significantly bigger than ours, for
example for every reconstruction performed with the full data set our
uncertainty in $H_{0}$ is reduced by at least a factor of 2 with respect to
these methods.

  * •

The iterative process carried out to calibrate the PMCT data causes a
preference for small values of $H_{0}$, to see this let us define the
following [58] quantity

$T\equiv\frac{|H_{0i}-H_{0j}|}{\sqrt{\sigma^{2}_{i}+\sigma^{2}_{j}}},$ (2.17)

where $H_{0i}$ and $H_{0j}$ are the mean values of the Hubble parameter and
$\sigma_{i}$, $\sigma_{j}$ are the corresponding 1$\sigma$ uncertainties. By
considering $H_{0}^{\text{Planck}}$, $H_{0}^{\text{R}}$ and
$H_{0}^{\text{Prior}}=69.39\pm 5.18$333We are using the $H_{0}$ for the Cauchy
kernel obtained form Cosmic Chronometers obtained in [36]., we calculate the
tension between them:

    * –

Tension between $H_{0}^{\text{Planck}}$ and $H_{0}^{\text{Prior}}$: $T=0.382$.

    * –

Tension between $H_{0}^{\text{R}}$ and $H_{0}^{\text{Prior}}$: $T=0.76$.

We can see that the prior gives a priority to $H_{0}^{\text{Planck}}$ (low
values of $H_{0}$). By repeating the process several times this preference
increases and seems to be confirmed by their final results reported. We avoid
this problem in our proposal by calibrating Pantheon data for every choice of
$H_{0}$ that we consider. This is extremely important for our work since we
are studying the Hubble tension effects.

In Tables 2-3 we report our results for $H_{0}$, along with a comparison with
three methods described. The results are present as follows: In the first
column of Table 2 we show the data sets used to perform our reconstructions,
in the second column we show the value of $H_{0}$ obtained from each of the
reconstructions using the squared exponential kernel, and in columns 3, 4 and
5, we present a comparison with methods mentioned. Table 3 follows has the
same structure as Table 2, the only difference is that the Table 3 presents
the results using the Cauchy kernel, finally the symbols appearing in the
tables denotes the following:

  * $\bigstar$: denotes that our result agrees with the Reference reported at the top of the Table.

  * $\bigotimes$: denotes that our result does not match the Reference at the top of the Table.

  * $-$: denotes that our result was not reported in the Reference at the top of the Table.

  * $\bullet$: denotes that our result has not been previously reported in the literature.

$H_{0}$ results for the squared exponential covariance function (2.2)
---
Data sets |  Our method |  Method (1) [34] |  Method (2) [35] |  Method (3) [36]
CC |  $67.32\pm 4.74$ |  – |  $\bigstar$ |  $\bigstar$
CC + CL |  $67.08\pm 3.01$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{R}$ |  $72.94\pm 0.40$ |  – |  – |  –
Pantheon+ $H_{0}^{R}$ \+ CC |  $72.83\pm 0.41$ |  – |  $\bigstar$ |  –
Pantheon+ $H_{0}^{R}$ \+ CC + CL |  $73.48\pm 0.36$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{HW}}$ |  $73.00\pm 0.41$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC |  $72.97\pm 0.42$ |  – |  $70.85\pm 1.19$ $\bigotimes$ |  –
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC + CL |  $73.63\pm 0.38$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{Planck}}$ |  $67.11\pm 0.37$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{Planck}}$ \+ CC |  $67.10\pm 0.38$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{Planck}}$ \+ CC + CL |  $66.77\pm 0.33$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ |  $69.53\pm 0.41$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ \+ CC |  $69.51\pm 0.42$ |  – |  $\bigstar$ |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ \+ CC + CL |  $69.53\pm 0.37$ $\bullet$ |  – |  – |  –
Table 2: Values of $H_{0}$ [km/s/Mpc] with their corresponding 1$\sigma$
uncertainties inferred via the GP reconstructions. Notice that our result for
the Pantheon+ $H_{0}^{\text{HW}}$\+ CC data set is at $1.68\sigma$ tension
with respect method (3) [35].  $H_{0}$ results for the Cauchy covariance
function (2.3)
---
Data sets |  Our method |  Method (1) [34] |  Method (2) [35] |  Method (3) [36]
CC |  $69.43\pm 5.23$ |  – |  $\bigstar$ |  $\bigstar$
CC + CL |  $67.41\pm 3.17$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{R}$ |  $72.92\pm 0.42$ |  – |  – |  –
Pantheon+ $H_{0}^{R}$ \+ CC |  $72.86\pm 0.42$ |  – |  $\bigstar$ |  –
Pantheon+ $H_{0}^{R}$ \+ CC + CL |  $73.51\pm 0.37$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{HW}}$ |  $72.98\pm 0.43$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC |  $73.00\pm 0.44$ |  – |  $70.89\pm 1.20$ $\bigotimes$ |  –
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC + CL |  $73.67\pm 0.39$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{Planck}}$ |  $67.10\pm 0.39$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{Planck}}$ \+ CC |  $67.13\pm 0.39$ $\bullet$ |  – |  – |  –
Pantheon + $H_{0}^{\text{Planck}}$ \+ CC + CL |  $66.81\pm 0.35$ $\bullet$ |  - |  – |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ |  $69.52\pm 0.43$ $\bullet$ |  – |  – |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ \+ CC |  $69.55\pm 0.44$ |  – |  $\bigstar$ |  –
Pantheon+ $H_{0}^{\text{TRGB}}$ \+ CC + CL |  $69.58\pm 0.38$ $\bullet$ |  – |  – |  –
Table 3: Values of $H_{0}$ [km/s/Mpc] with their corresponding 1$\sigma$
uncertainties inferred via the GP reconstructions. Notice that our result for
the Pantheon+ $H_{0}^{\text{HW}}$\+ CC data set is at 1.65$\sigma$ tension
with respect method (3) [35].

##  3 Model-independent constraints on density parameters

Since the Hubble constant and the dark energy density at present time are
correlated [4], it is expected that the Hubble tension also manifests itself
as a tension between the value of the dark energy density inferred from the
fit of $\Lambda$CDM to the observations of the CMB and the value obtained with
the reconstructions. To obtain the dark energy density from the
reconstructions we will assume the following: a spatially flat universe with
matter content modeled as a dust-like fluid, in addition we will assume that
dark energy can be modeled by an effective fluid with a parametrised equation
of state $\omega(z)$. Under these assumptions the standard Friedmann equations
have the following form

$\displaystyle
H^{2}=H_{0}^{2}[\Omega_{m}^{0}(1+z)^{3}+\Omega_{\text{eff}}^{0}f(z)],$ (3.1)
$\displaystyle\frac{\ddot{a}}{a}=-\frac{H^{2}}{2}\left(\Omega_{m}+\Omega_{\text{eff}}(1+3\omega)\right),$
(3.2) $\displaystyle\dot{\rho}+3H(1+P)=0,$ (3.3)

where the index “0” denotes quantities evaluated at present cosmic time.
Following the standard definition, $\Omega_{\text{eff}}$ denotes the dark
energy density parameter and it changes as function of redshift as
$\Omega_{\text{eff}}(z)=\Omega_{\text{eff}}^{0}f(z)$, with
$f(z)=\exp\left(3\int_{0}^{z}\frac{1+\omega}{1+z^{\prime}}dz^{\prime}\right)$.
In particular for $\Lambda$CDM $f(z)=1$ for all z. According to the constraint
$\Omega_{m}+\Omega_{\text{eff}}=1$, we can rewrite (3.1) as

$H^{2}=H^{2}_{0}[(1-\Omega_{\text{eff}}^{0})(1+z)^{3}+\Omega_{\text{eff}}^{0}f(z)].$
(3.4)

If we evaluate (3.4) at $z=0$, we have that the function $f(z)$ must obey
$f(z=0)=1$. With this ansatz we can rewrite (3.1) as follows

$H^{2}-H^{2}_{0}(1+z)^{3}=H_{0}^{2}\Omega_{\text{eff}}^{0}[f(z)-(1+z)^{3}].$
(3.5)

Notice that if $z=0$, both sides of the above equation vanish and therefore we
cannot get any information from $\Omega_{\text{eff}}^{0}$. However, we expect
that $H(z)$ is a smooth (differentiable) function, then $f(z)$ does not have
abrupt changes, i.e. $f(z\approx 0)\approx 1$. Moreover, in the limit
$z\rightarrow 0$ we found that

$H^{2}-H^{2}_{0}(1+z)^{3}\approx
H_{0}^{2}\Omega_{\text{eff}}^{0}[1-(1+z)^{3}],$ (3.6)

therefore

$\Omega_{\text{eff}}^{0}\approx\frac{H^{2}-H^{2}_{0}(1+z)^{3}}{H_{0}^{2}[1-(1+z)^{3}]}.$
(3.7)

We proceed to define $\Delta z$ as

$\Delta z=\frac{z_{f}-z_{i}}{N},$ (3.8)

where $N$ is the number of points placed in the reconstruction, and the
subscripts $f$, $i$ denote final and initial values of $z$, respectively.
Taking $z_{i}=0$, we see that the smallest $z$ that we have access with the
reconstruction is $\Delta z$, therefore

$\Omega_{\text{eff}}^{0}=\frac{H^{2}-H^{2}_{0}(1+\Delta
z)^{3}}{H_{0}^{2}[1-(1+\Delta z)^{3}]}.$ (3.9)

We can obtain the error propagation of the latter equation by writing $H_{0}$
and $H(\Delta z)$ as

$\displaystyle H_{0}$ $\displaystyle=$
$\displaystyle\overline{H}_{0}\pm\sigma_{H_{0}},$ (3.10) $\displaystyle
H(\Delta z)$ $\displaystyle=$
$\displaystyle\overline{H}_{1}\pm\sigma_{H_{1}},$ (3.11)

where the first term (denoted by a bar on top it) represents the mean value
,and the second is the 1$\sigma$ uncertainty of $H_{0}$ and $H(\Delta z)$,
respectively. With these definitions we can perform the error propagation to
obtain

$\Omega_{\text{effe}}^{0}\approx\frac{\sqrt{2}}{(1+\Delta
z)^{3}-1}\left(\frac{\overline{H}_{1}^{2}}{\overline{H}_{0}^{4}}\sigma_{H_{1}}^{2}+\frac{\overline{H}_{1}^{4}}{\overline{H}_{0}^{6}}\sigma_{H_{0}}^{2}\right)^{1/2},$
(3.12)

by definition $\Omega_{\text{effe}}^{0}$ is the $1\sigma$ error for
$\Omega_{\text{eff}}^{0}$. We can rewrite this error in terms of the number of
points $N$ in the reconstruction

$\Omega_{\text{
effe}}^{0}\approx\frac{\sqrt{2}N^{3}}{(N+z_{f})^{3}-N^{3}}\left(\frac{\overline{H}_{1}^{2}}{\overline{H}_{0}^{4}}\sigma_{H_{1}}^{2}+\frac{\overline{H}_{1}^{4}}{\overline{H}_{0}^{6}}\sigma_{H_{0}}^{2}\right)^{1/2},$
(3.13)

from where it is clear that the error depends on the number $N$ of points
placed in the reconstruction, which leads us to the following issue: we cannot
increase the number of data points in the reconstruction arbitrarily since
this would increase the statistical uncertainty in the calculation of
$\Omega_{\text{eff}}^{0}$, but either we can make $N$ too small or (3.7) would
not remain valid.

Figure 3: In blue it is shown how $\Omega_{\text{eff}}^{0}$ varies as a
function of the number of points $N$ placed in the reconstruction, while in
red it is shown how $\Delta z$ varies as function of $N$. The intersection
indicates which is the $N$ that minimizes uncertainty and at the same time
gives us a $\Delta z$ small enough to have a decent approximation. We get that
$N_{\text{optimal}}\equiv N_{0}\approx 43$, which corresponds to $\Delta
z=2.37/N_{0}.$

In order to solve this issue we define the following quantity:

$E\equiv\frac{H_{\text{approx}}(\Delta z)}{H_{\text{recons}}(\Delta z)},$
(3.14)

where $H_{\text{approx}}$ is the $H(z)$ given by (3.6) and $H_{\text{recons}}$
is the $H(z)$ given by the reconstruction. If $E=1$, (3.7) is completely
valid. If $E=1$, (3.7) is completely valid, if $E<1$ we will underestimate the
value of $\Omega_{\text{effe}}^{0}$ and if $E>1$ we will overestimate the
value $\Omega_{\text{effe}}^{0}$, since the Eq. (3.7) is no longer valid.
Therefore, $E$ gives us a measure of how good the approximation given by (3.7)
is. We present the results for $\Omega_{\text{eff}}$ in the following form:

$\Omega_{\text{eff}}^{0}=\overline{\Omega}_{\text{eff}}^{0}\pm\Omega_{\text{effs}}^{0}\pm\Omega_{\text{effn}}^{0},$
(3.15)

where $\overline{\Omega}_{\text{eff}}^{0}$ is the mean value of
$\Omega_{\text{eff}}^{0}$, $\Omega_{\text{effs}}^{0}$ is the uncertainty at
1$\sigma$ and $\Omega_{\text{effn}}^{0}$ is the error due to the approximation
given by (3.7). This error is obtained through (3.14). In Table 4 we show the
results of calculating (3.15) with the reconstructions. It is interesting to
note that $\Omega_{\text{effn}}^{0}<<\Omega_{\text{effs}}^{0}$ regardless of
the reconstruction.

Constrains on $\Omega_{\text{eff}}^{0}$ via Gaussian Processes
---
Data sets used |  Squared exponential kernel (2.2) |  Cauchy kernel (2.3)
CC |  $0.651\pm 0.194\pm 8.29\%$ |  $0.641\pm 0.191\pm 8.30\%$
CC+ CL |  $0.654\pm 0.156\pm 6.90\%$ |  $0.660\pm 0.161\pm 7.45\%$
Pantheon+ $H_{0}^{\text{R}}$ \+ CC+ CL |  $0.708\pm 0.052\pm 3.40\%$ |  $0.705\pm 0.054\pm 3.34\%$
Pantheon+ $H_{0}^{\text{HW}}$ \+ CC+ CL |  $0.712\pm 0.054\pm 3.46\%$ |  $0.709\pm 0.055\pm 3.41\%$
Pantheon+ $H_{0}^{\text{Planck}}$ \+ CC+ CL |  $0.663\pm 0.053\pm 2.15\%$ |  $0.644\pm 0.054\pm 2.33\%$
Pantheon + $H_{0}^{\text{TRGB}}$ \+ CC+ CL |  $0.664\pm 0.056\pm 2.67\%$ |  $0.668\pm 0.059\pm 2.73\%$
Table 4: Obtained values for the density parameter $\Omega_{\text{eff}}^{0}$
with their corresponding 1$\sigma$ uncertainties. The first column shows the
data sets used to perform each reconstruction. The second column shows the
obtained value for $\Omega_{\text{eff}}^{0}$ using the squared exponential
kernel. The third column shows the obtained value for
$\Omega_{\text{eff}}^{0}$ using the Cauchy kernel. The first value in the
second and third columns represents the mean value of
$\Omega_{\text{eff}}^{0}$, the second value is the 1$\sigma$ statistical
uncertainty and the third value is the numerical error associated with the
calculation of $\Omega_{\text{eff}}^{0}$.

Using the constraint equation $\Omega_{m}^{0}+\Omega_{\Lambda}^{0}=1$, and the
reconstruction that minimizes the $\chi^{2}$ statistic in Table 1, it can be
shown that

$\displaystyle\Omega_{m}^{0}=0.336\pm 0.0563\pm
2.67\%,\quad\rightarrow\quad\text{Squared Exponential},$ (3.16)
$\displaystyle\Omega_{m}^{0}=0.332\pm 0.0594\pm
2.73\%,\quad\rightarrow\quad\text{Cauchy function},$ (3.17)

both values coincide with the $\Omega_{m}^{0}$ reported by the Planck 2018
collaboration, so there is not a tension effect on this parameter. Although
this is due to the size of the 1$\sigma$ uncertainties in
$\Omega_{\text{eff}}^{0}$.

##  4 Updated constraints on dark energy models from Horndeski Cosmology

In the last few years, Horndeski theory of gravity [30] has been studied for
scalar-field dark energy [59, 60, 61, 62, 63] where it has been show that its
solutions can be ghost-free [64] and undergo tachyonic instability roles [65].
Also, Galileon cosmology from Horndeski theory has been done in Refs. [66, 19,
67, 68]. Furthermore, extensions up to quartic order have been considered in
Ref. [69, 70].

The Horndeski theory is the most general covariant scalar-tensor theory with a
single scalar field whose Lagrange equation has derivatives of the metric and
the scalar field up to second order. The action of this theory is given by

$S[g_{ab},\phi]=\int d^{4}x\sqrt{-g}\left[\sum^{5}_{i=2}\frac{1}{8\pi
G}\mathcal{L}[g_{ab},\phi]+\mathcal{L}_{m}[g_{ab}]\right],$ (4.1)

where the four Lagrangians $\mathcal{L}_{i}$, are given by:

$\displaystyle\mathcal{L}_{2}$ $\displaystyle=$ $\displaystyle G_{2}(\phi,X),$
(4.2) $\displaystyle\mathcal{L}_{3}$ $\displaystyle=$ $\displaystyle-
G_{3}(\phi,X)\square\phi,$ (4.3) $\displaystyle\mathcal{L}_{4}$
$\displaystyle=$ $\displaystyle
G_{4}(\phi,X)R+G_{4X}\left[(\square\phi)^{2}-\phi_{;ab}\phi^{;ab}\right],$
(4.4) $\displaystyle\mathcal{L}_{5}$ $\displaystyle=$ $\displaystyle
G_{5}(\phi,X)G_{ab}\phi^{;ab}-\frac{G_{5X}(\phi,X)}{6}\left[(\square\phi)^{3}+2\phi_{;\mu}^{\nu}\phi_{;\nu}^{\alpha}\phi_{;\alpha}^{\mu}-3\phi_{;ab}\phi^{;ab}\square\phi\right],$
(4.5)

with $G_{2}$, $G_{3}$, $G_{4}$ and $G_{5}$ arbitrary functions of the scalar
field $\phi$ and its kinetic term $2X=-\partial_{\mu}\phi\partial^{\mu}\phi$.
In the above equations $R$ the Ricci scalar, $G_{\mu\nu}$ the Einstein tensor.
We use the subscripts $X$ and $\phi$ to denote partial derivatives, with
respect to $X$ and $\phi$. The $G_{i}$ functions are completely arbitrary,
however, in the next section will be explained how we have constrained these
functions.

The detection of the gravitational waves GW170817 [71], the Gamma Ray Burst
GRB 170817A and the a transient Swope Supernova Survey SSS17a [71, 72] places
a strong constraint on the speed at which gravitational waves propagate [73,
74, 75]

$-3\times 10^{-15}<c_{GW}-1<7\times 10^{-16}.$ (4.6)

This constraint motivates to only consider only a subclass of Horndeski
models, that is, those that meet the condition $c_{GW}=1$, these models [30]
are given by $G_{4X}=0$ and $G_{5}=0$. The Lagrangian for this models has the
form

$\mathcal{L}=G_{2}(\phi,X)-G_{3}(\phi,X)\square\phi+G_{4}(\phi)R.$ (4.7)

We analyze the Lagrangian (4.7) using a spatially flat, homogeneous and
isotropic Friedmann-Lemaitre-Robertson-Walker metric,
$ds^{2}=dt^{2}-a(t)^{2}(dx^{2}+dy^{2}+dz^{2})$, in addition to describe the
matter content of the Universe, we use the energy-momentum tensor $T_{\mu\nu}$
for a perfect fluid $T_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ and since
in this work we are focusing on cosmic late-time aspects, the density of
radiation can be neglected, then the Friedmann equations are given by

$\displaystyle H^{2}=\frac{8\pi G}{3}\rho_{m}+\varepsilon,$ (4.8)
$\displaystyle\dot{\varepsilon}+3H(\varepsilon+P)=0,$ (4.9)

where

$\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\frac{1}{3}\left[-G_{2}+2X(G_{2X}-G_{3\phi})\right]+2H\dot{\phi}\left[-G_{4\phi}+XG_{3X}\right]+H^{2}\left[1-2G_{4}\right],$
(4.10) $\displaystyle P$ $\displaystyle=$
$\displaystyle\frac{2}{3}\dot{H}[1-2G_{4}]-H^{2}(1-2G_{4})+\frac{2}{3}(H\dot{\phi}+\ddot{\phi})\left[G_{4\phi}-XG_{3X}\right]$
(4.11)
$\displaystyle+\frac{1}{3}\left[G_{2}-2X(G_{3\phi}-2G_{4\phi\phi})\right]+\frac{2H\dot{\phi}}{3}\left[G_{4\phi}+XG_{3X}\right].$

From this point forward we will analyse two different cases for the $G_{i}$’s
using the potential

$V=C\exp(n\phi),$ (4.12)

with $n$ and $C$ free parameters.

###  4.1 Quintessence

In this first example we will analyse the model

$\displaystyle G_{2}(\phi,X)=X-V(\phi),\quad G_{3}=0,\quad G_{4}=1/2,$ (4.13)

which is a Horndeski model equivalent to a Quintessence model [27]. We can
obtain the field equations for this model if we substitute (4.13) in
(4.8)-(4.11) to obtain

$\displaystyle H^{2}=\frac{8\pi G}{3}\rho_{m}+\varepsilon,$ (4.14)
$\displaystyle\dot{\varepsilon}+3H(\varepsilon+P)=0,$ (4.15)

where

$\displaystyle\varepsilon=\frac{1}{3}[-G_{2}+2XG_{2X}],$ (4.16) $\displaystyle
P=\frac{1}{3}[G_{2}].$ (4.17)

In this model, the effective equation of state for the scalar field is given
by

$\omega_{\phi}=\frac{\frac{1}{2}\dot{\phi}^{2}-V(\phi)}{\frac{1}{2}\dot{\phi}^{2}+V(\phi)},$
(4.18)

which in the limit where $\dot{\phi}\approx 0$, $\omega_{\phi}\approx-1.$
Since we need that our model reduces to $\Lambda$CDM at least for $z=0$, it is
necessary that $\dot{\phi}(z=0)=0$, which represents an initial condition for
the scalar field $\phi$. We can obtain the value of $C$ from (4.12) by
evaluating (4.14) at $z=0$

$V(\phi(z=0))=3H_{0}^{2}[1-\Omega_{m}^{0}]\equiv C.$ (4.19)

Substituting (4.16) and (4.17) in the continuity equation (4.15), it is
possible to show that the scalar field must satisfy the Klein-Gordon (KG)
equation

$\ddot{\phi}+3H\dot{\phi}+V_{\phi}=0.$ (4.20)

It is useful to change the KG from the cosmic time to redshift:

$\phi^{\prime\prime}[H(1+z)]^{2}+\phi^{\prime}[H^{\prime}(1+z)+H](1+z)H-3H^{2}(1+z)\phi^{\prime}+V_{\phi}=0.$
(4.21)

This equation requires two initial conditions: the first one was obtained
previously, i.e $\dot{\phi}(z=0)=0$ and the second changes depending on the
potential that we choose, specifically for potentials of the type (4.12) must
satisfy $\phi(z=0)=0$. To obtain the value of $n$ that best fits the
reconstructions, we minimize the following quantity

$\chi^{2}=\sum_{i}^{N_{\text{recons}}}\frac{\left(H_{\text{recons}}(z_{i})-H_{\text{Q}}(z_{i};n)\right)^{2}}{\sigma_{\text{recons}}(z_{i})^{2}},$
(4.22)

where $N$ represents the number of points in the reconstruction,
$H_{\text{recons}}(z_{i})$ is the Hubble parameter of the reconstruction
evaluated in the redshift $z_{i}$, $\sigma_{\text{recons}}(z_{i})$ is the
1$\sigma$ uncertainty of the reconstruction and $H_{\text{Q}}(z_{i};n)$ is the
Hubble parameter for Quintessence. Notice that the value of $C$ can be
obtained by considering the condition (4.19) evaluate at $z=0$, therefore our
free parameter will be $n$. For the squared exponential kernel we obtain
$n=1.0\times 10^{-9}$, while for the Cauchy kernel we find $n=1.1\times
10^{-9}$.

In Fig. 4 we show the evolution of $H(z)$ by solving (4.21) for the potential
(4.12) with the mean value of $H_{0}$ obtained from the table 2 and the mean
value of $\Omega^{0}_{m}$ obtained from the table 3, that correspond to the
reconstruction Pantheon + $H^{\text{TRGB}}_{0}$ \+ CC + CL. This result is
compared with the Hubble parameter obtained from the reconstruction named
Pantheon + $H^{\text{TRGB}}_{0}$ \+ CC + CL.

Figure 4: Evolution of $H(z)$ obtained by solving (4.21) with the potential
$V=C\exp(n\phi)$, compared with the $H(z)$ obtained from the reconstruction
with the squared exponential (left) and Cauchy (right) kernels at 2$\sigma$
level. The values of $n$ are obtained minimizing (4.22).

Notice that the Quintessence model is within the 2$\sigma$ confidence
intervals for the reconstructions on $0<z<2$, however for $z>2$, this model
cannot describe the reconstructions on this confidence interval.

###  4.2 K-essence

In this section we take the model

$G_{2}=X-V,\quad G_{3}=-c_{1}V,\quad G_{4}=1/2,$ (4.23)

where $c_{1}$ is a positive constant with units of $s^{2}$, in particular, we
choose $c_{1}=1s^{2}$. If $G_{3}$ is an increasing function on the redshift,
then values of $c_{1}$ close to 1$s^{2}$ tend to decrease the expansion rate
of the universe with respect the expansion rate given by $\Lambda$CDM.
Meanwhile, for large values of $c_{1}$ our model goes to $\Lambda$CDM. Notice
that this model can be reduced to a K-essence model by performing an
integration by parts of the corresponding action (4.1), where the term $G_{3}$
vanishes444See Equation (2.2) from [76].:

$G_{2}=X-V-2Xc_{1}V_{\phi},\quad G_{4}=1/2,.$ (4.24)

The Friedmann equations for this model are given by

$\displaystyle H^{2}=\frac{8\pi G}{3}\rho_{m}+\varepsilon,$ (4.25)
$\displaystyle\dot{\varepsilon}+3H(\varepsilon+P)=0,$ (4.26)

with

$\displaystyle\varepsilon=\frac{1}{3}[V+X-2XV_{\phi}c_{1}],$ (4.27)
$\displaystyle P=\frac{1}{3}[X-V-2XV_{\phi}c_{1}].$ (4.28)

Therefore, the equation of state for the scalar field has the form

$\displaystyle\omega_{\phi}=\frac{X-V-2XV_{\phi}c_{1}}{V+X-2XV_{\phi}c_{1}}.$
(4.29)

If $\dot{\phi}(z=0)=0$, (4.29) reduces to $\omega_{\phi}=-1$. We can obtain
the initial condition for $\phi$ if we evaluate (4.25) at $z=0$. In
particular, for the potential $V=C\exp(n\phi)$, it is necessary that
$\phi(z=0)=0$, with $C=3H_{0}^{2}[1-\Omega_{m}^{0}]$. To obtain the evolution
of $\phi(z)$, it is useful to rewrite (4.26) as

$-\frac{1}{3}\frac{\mathrm{d}}{\mathrm{d}z}\left[V+H^{2}(1+z)^{2}\phi^{\prime
2}(\frac{1}{2}-V_{\phi}c_{1})\right]+(1+z)\left[1-2V_{\phi}c_{1}\right]H^{2}\phi^{\prime
2}=0.$ (4.30)

To obtain the value of $n$ that best fits the reconstructions, we minimize the
following quantity

$\chi^{2}=\sum_{i}^{N_{\text{recons}}}\frac{\left[(H_{\text{recons}}(z_{i})-H_{\text{EQ}}(z_{i};n)\right]^{2}}{\sigma_{\text{recons}}(z_{i})^{2}},$
(4.31)

where $N$ represents the number of points in the reconstruction,
$H_{\text{recons}}(z_{i})$ is the Hubble parameter of the reconstruction
evaluated at redshift $z_{i}$, $\sigma_{\text{recons}}(z_{i})$ is the
1$\sigma$ uncertainty of the reconstruction and $H_{\text{EQ}}(z_{i};n)$ is
the Hubble parameter for K-essence. For the squared exponential kernel we get
$n=79.31$ while for the Cauchy kernel $n=77.41$.

In Fig. 5 we show the evolution of the Hubble parameter obtained by solving
(4.30) for the potential (4.12) with the mean value of $H_{0}$ and the mean
value of $\Omega^{0}_{m}$ obtained from the reconstruction carried out with
the Pantheon + $H^{\text{TRGB}}_{0}$ \+ CC + CL data. This result is compared
with the $H(z)$ obtained from the reconstruction carried out with the data of
Pantheon + $H^{\text{TRGB}}_{0}$ \+ CC + CL.

Reduced $\chi^{2}$ at 1$\sigma$
---
Model |  Squared exponential kernel (2.2) |  Cauchy kernel (2.3)
Quintessence |  $9.73$ |  $8.47$
K-essence |  $7.78$ |  $6.70$
Table 5:  Reduced $\chi^{2}$-statistics between the Pantheon +
$H^{\text{TRGB}}_{0}$ \+ CC + CL reconstruction and the Horndeski model using
the squared exponential and Cauchy kernel, respectively. Reduced $\chi^{2}$ at
2$\sigma$
---
Model |  Squared exponential kernel (2.2) |  Cauchy kernel (2.3)
Quintessence |  $2.43$ |  $2.11$
K-essence |  $1.94$ |  $1.67$
Table 6: Reduced $\chi^{2}$-statistics between the Horndeski model and the
Pantheon + $H^{\text{TRGB}}_{0}$ \+ CC + CL reconstruction taking the
uncertainties at the two sigma interval using the squared exponential and
Cauchy kernel, respectively.

Table 5 shows a comparison between the value of the reduced
$\chi^{2}$-statistic for each model. It is observed that both models have a
high value of the reduced $\chi^{2}$-statistics, this may be due to a possible
overfitting of the reconstructions to the observations in the region $0<z<1$,
this results are reported in Table 6 as a comparison of the reduced
$\chi^{2}$-statistic assuming the value of the uncertainties of the
reconstruction at 2$\sigma$ level. It is observed that in this case we have a
reasonable value for the reduced $\chi^{2}$-statistic for the K-essence model.

Figure 5:  Evolution of $H(z)$ obtained from (4.30) with the potential
$V=C\exp(n\phi)$, compared with the $H(z)$ obtained from the reconstruction
with the squared exponential (left) and Cauchy (right) kernels at 2$\sigma$
level. The values of $n$ are obtained minimizing the equation (4.31).

##  5 Conclusions

In this paper we improved reconstructions of the late time expansion of the
universe using Gaussian Processes. We show the possibility of using our
reconstructions to study the Hubble tension and to obtain constrains on the
Galileon $G_{2}$ from Horndeski theory of gravity. Our work presents an
improvement over previous methods reported in the literature, and we describe
that the use of compressed data sets to perform Gaussian Processes is not the
best option, since the methods that use this type of data sets such as the
Covariance functions reconstruction method and the Reconstruction with
weighted polynomial regression method have uncertainties in $H_{0}$ twice
larger than ours.

We showed that the obtained values for $H_{0}$ from each reconstruction reach
differences up to 13$\sigma$ depending on the prior used to calibrate the
SNeIa data given by Pantheon. Using the $\chi^{2}$-statistics, we found that
the reconstruction that best fit the observations is the combination Pantheon
+ $H^{\text{TRGB}}_{0}$ \+ CC + CL, with Pantheon dataset calibrated using
$H_{0}$ measured by The Carnegie-Chicago Hubble Program [45]. Our results
reports that the $H_{0}$ value that best describes the observations of the
late universe is $H_{0}=69.53\pm 0.37$, and this value links a tension of
3.4$\sigma$ with respect the value of $\Lambda$CDM inferred from Planck 2018
data. However, it is still possible to improve the quality of the
reconstructions presented in this work by recalibrating the Clustering of
Galaxies sample.

Furthermore, assuming a parametric form for the dark energy equation of state
we obtain new constraints on the density parameter for this exotic fluid at
the present time. By defining (3.14) we derive a way to estimate the numerical
error associated with the $\Omega_{\text{eff}}^{0}$ calculation. We also show
that the density parameters obtained from the reconstructions are not in
tension with those obtained by the Planck 2018 collaboration assuming a flat
$\Lambda$CDM, although this is mainly due to the size of our 1$\sigma$
uncertainties.

Finally, we showed that Horndeski models as Quintessence type under the
potential (4.12) are not able to fit the reconstructions of the Hubble flow
within the 2$\sigma$ confidence interval, therefore there is a possibility
that these models are not good candidates to describe the cosmic late time
observations. The K-essence model discussed in Sec.4.2 with an exponential-
type potential fits reconstructions at 2$\sigma$ confidence interval and
therefore represents a good candidate to describe the late time expansion.
Furthermore, in contrast with the results presented in [77], where an
exponential potential is analysed555Also so-called Modulus model. with a fixed
positive $n$ prior value, in our analysis we obtained best fit positive values
for $n$ in both Quintessence and K-essence cases.

It is expected that the K-essence model modify the history of the large-scale
structure formation, and therefore it will be interesting to analyse the first
order linear perturbations of this model since these will tell us if it is
also a candidate to solve the $S_{8}$ tension [78]. These analyses are an
ongoing project, and they will be reported elsewhere.

## Acknowledgments

CE-R acknowledges the Royal Astronomical Society as FRAS 10147. CE-R and M.R
are supported by DGAPA-PAPIIT-UNAM Project IA100220. This article is also
based upon work from COST action CA18108, supported by COST (European
Cooperation in Science and Technology). The simulations were performed in
Centro de cómputo Tochtli-ICN-UNAM. The authors thank the anonymous reviewer
whose comments/suggestions helped improve and clarify this manuscript.

## References

  * [1] Adam G. Riess et al.  A 2.4% Determination of the Local Value of the Hubble Constant.  Astrophys. J., 826(1):56, 2016.
  * [2] Adam G. Riess et al.  Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant.  Astrophys. J., 861(2):126, 2018.
  * [3] S. Birrer et al.  H0LiCOW - IX. Cosmographic analysis of the doubly imaged quasar SDSS 1206+4332 and a new measurement of the Hubble constant.  Mon. Not. Roy. Astron. Soc., 484:4726, 2019.
  * [4] N. Aghanim et al.  Planck 2018 results. VI. Cosmological parameters.  Astron. Astrophys., 641:A6, 2020.
  * [5] N. Aghanim et al.  Planck intermediate results. LI. Features in the cosmic microwave background temperature power spectrum and shifts in cosmological parameters.  Astron. Astrophys., 607:A95, 2017.
  * [6] Kevin Aylor, MacKenzie Joy, Lloyd Knox, Marius Millea, Srinivasan Raghunathan, and W. L. Kimmy Wu.  Sounds Discordant: Classical Distance Ladder & $\Lambda$CDM -based Determinations of the Cosmological Sound Horizon.  Astrophys. J., 874(1):4, 2019.
  * [7] Eleonora Di Valentino, Olga Mena, Supriya Pan, Luca Visinelli, Weiqiang Yang, Alessandro Melchiorri, David F Mota, Adam G Riess, and Joseph Silk.  In the Realm of the Hubble tension $-$ a Review of Solutions.  arXiv preprint arXiv:2103.01183, 2021.
  * [8] Eleonora Di Valentino, Luis A Anchordoqui, Yacine Ali-Haimoud, Luca Amendola, Nikki Arendse, Marika Asgari, Mario Ballardini, Elia Battistelli, Micol Benetti, Simon Birrer, et al.  Cosmology intertwined II: the Hubble constant tension.  arXiv preprint arXiv:2008.11284, 2020.
  * [9] Wendy L. Freedman.  Cosmology at a Crossroads.  Nature Astron., 1:0121, 2017.
  * [10] Eleonora Di Valentino, Alessandro Melchiorri, Eric V. Linder, and Joseph Silk.  Constraining Dark Energy Dynamics in Extended Parameter Space.  Phys. Rev. D, 96(2):023523, 2017.
  * [11] Abdolali Banihashemi, Nima Khosravi, and Amir H. Shirazi.  Ginzburg-Landau Theory of Dark Energy: A Framework to Study Both Temporal and Spatial Cosmological Tensions Simultaneously.  Phys. Rev. D, 99(8):083509, 2019.
  * [12] Suresh Kumar and Rafael C. Nunes.  Probing the interaction between dark matter and dark energy in the presence of massive neutrinos.  Phys. Rev. D, 94(12):123511, 2016.
  * [13] G. E. Addison, D. J. Watts, C. L. Bennett, M. Halpern, G. Hinshaw, and J. L. Weiland.  Elucidating $\Lambda$CDM: Impact of Baryon Acoustic Oscillation Measurements on the Hubble Constant Discrepancy.  Astrophys. J., 853(2):119, 2018.
  * [14] Eleonora Di Valentino, Alessandro Melchiorri, and Olga Mena.  Can interacting dark energy solve the $H_{0}$ tension?  Phys. Rev. D, 96(4):043503, 2017.
  * [15] Shadab Alam et al.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample.  Mon. Not. Roy. Astron. Soc., 470(3):2617–2652, 2017.
  * [16] Gong-Bo Zhao et al.  Dynamical dark energy in light of the latest observations.  Nature Astron., 1(9):627–632, 2017.
  * [17] Luisa G. Jaime, Mariana Jaber, and Celia Escamilla-Rivera.  New parametrized equation of state for dark energy surveys.  Phys. Rev. D, 98(8):083530, 2018.
  * [18] Janina Renk, Miguel Zumalacárregui, Francesco Montanari, and Alexandre Barreira.  Galileon gravity in light of ISW, CMB, BAO and H0 data.  Journal of Cosmology and Astroparticle Physics, 2017(10):020–020, Oct 2017.
  * [19] Simone Peirone, Noemi Frusciante, Bin Hu, Marco Raveri, and Alessandra Silvestri.  Do current cosmological observations rule out all covariant galileons?  Physical Review D, 97(6), Mar 2018.
  * [20] Simone Peirone, Giampaolo Benevento, Noemi Frusciante, and Shinji Tsujikawa.  Cosmological data favor Galileon ghost condensate over $\Lambda$CDM.  Phys. Rev. D, 100(6):063540, 2019.
  * [21] Noemi Frusciante, Simone Peirone, Luís Atayde, and Antonio De Felice.  Phenomenology of the generalized cubic covariant galileon model and cosmological bounds.  Physical Review D, 101(6), Mar 2020.
  * [22] Rafael C. Nunes.  Structure formation in $f(T)$ gravity and a solution for $H_{0}$ tension.  JCAP, 05:052, 2018.
  * [23] Celia Escamilla-Rivera and Jackson Levi Said.  Cosmological viable models in $f(T,B)$ theory as solutions to the $H_{0}$ tension.  Class. Quant. Grav., 37(16):165002, 2020.
  * [24] Geovanny A Rave-Franco, Celia Escamilla-Rivera, and Jackson Levi Said.  Dynamical complexity of the teleparallel gravity cosmology.  arXiv preprint arXiv:2101.06347, 2021.
  * [25] Marco Raveri.  Reconstructing Gravity on Cosmological Scales.  Phys. Rev. D, 101(8):083524, 2020.
  * [26] Giulia Gubitosi, Federico Piazza, and Filippo Vernizzi.  The Effective Field Theory of Dark Energy.  JCAP, 02:032, 2013.
  * [27] Emilio Bellini, Ignacy Sawicki, and Miguel Zumalacárregui.  hi_class background evolution, initial conditions and approximation schemes.  Journal of Cosmology and Astroparticle Physics, 2020(02):008, 2020\.
  * [28] Guillermo Ballesteros, Alessio Notari, and Fabrizio Rompineve.  The H0 tension: $\Delta$ GN vs. $\Delta$ Neff.  Journal of Cosmology and Astroparticle Physics, 2020(11):024, 2020\.
  * [29] William T Emond, Antoine Lehébel, and Paul M Saffin.  Black holes in self-tuning cubic Horndeski cosmology.  Physical Review D, 101(8):084008, 2020.
  * [30] Tsutomu Kobayashi.  Horndeski theory and beyond: a review.  Reports on Progress in Physics, 82(8):086901, 2019.
  * [31] Ryotaro Kase and Shinji Tsujikawa.  Dark energy in Horndeski theories after GW170817: A review.  International Journal of Modern Physics D, 28(05):1942005, 2019\.
  * [32] Miguel Zumalacarregui.  Gravity in the Era of Equality: Towards solutions to the Hubble problem without fine-tuned initial conditions.  Phys. Rev. D, 102(2):023523, 2020.
  * [33] Ariadna Montiel, Ruth Lazkoz, Irene Sendra, Celia Escamilla-Rivera, and Vincenzo Salzano.  Nonparametric reconstruction of the cosmic expansion with local regression smoothing and simulation extrapolation.  Phys. Rev. D, 89(4):043007, 2014.
  * [34] Ahmad Mehrabi and Spyros Basilakos.  Does $\Lambda$CDM really be in tension with the Hubble diagram data?  The European Physical Journal C, 80(7):1–8, 2020.
  * [35] Rebecca Briffa, Salvatore Capozziello, Jackson Levi Said, Jurgen Mifsud, and Emmanuel N. Saridakis.  Constraining teleparallel gravity through Gaussian processes.  Class. Quant. Grav., 38(5):055007, 2020.
  * [36] Adrià Gómez-Valent and Luca Amendola.  H0 from cosmic chronometers and Type Ia supernovae, with Gaussian Processes and the novel Weighted Polynomial Regression method.  Journal of Cosmology and Astroparticle Physics, 2018(04):051, 2018\.
  * [37] Marina Seikel, Chris Clarkson, and Mathew Smith.  Reconstruction of dark energy and expansion dynamics using Gaussian processes.  Journal of Cosmology and Astroparticle Physics, 2012(06):036, 2012\.
  * [38] Noemi Frusciante, Simone Peirone, Luis Atayde, and Antonio De Felice.  Phenomenology of the generalized cubic covariant Galileon model and cosmological bounds.  Phys. Rev. D, 101(6):064001, 2020.
  * [39] Noemi Frusciante, Simone Peirone, Santiago Casas, and Nelson A. Lima.  Cosmology of surviving Horndeski theory: The road ahead.  Physical Review D, 99(6), Mar 2019.
  * [40] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * [41] Ming-Jian Zhang and Hong Li.  Gaussian processes reconstruction of dark energy from observational data.  The European Physical Journal C, 78(6):1–12, 2018.
  * [42] Marina Seikel, Sahba Yahya, Roy Maartens, and Chris Clarkson.  Using H(z) data as a probe of the concordance model.  Physical Review D, 86(8):083001, 2012.
  * [43] Celia Escamilla-Rivera and Júlio C Fabris.  Nonparametric reconstruction of the Om diagnostic to test $\Lambda$CDM.  Galaxies, 4(4):76, 2016.
  * [44] Kenneth C Wong, Sherry H Suyu, Geoff CF Chen, Cristian E Rusu, Martin Millon, Dominique Sluse, Vivien Bonvin, Christopher D Fassnacht, Stefan Taubenberger, Matthew W Auger, et al.  H0LiCOW–XIII. A 2.4 percent measurement of H0 from lensed quasars: 5.3 $\sigma$ tension between early-and late-Universe probes.  Monthly Notices of the Royal Astronomical Society, 498(1):1420–1439, 2020.
  * [45] Wendy L Freedman, Barry F Madore, Dylan Hatt, Taylor J Hoyt, In Sung Jang, Rachael L Beaton, Christopher R Burns, Myung Gyoon Lee, Andrew J Monson, Jillian R Neeley, et al.  The Carnegie-Chicago Hubble Program. VIII. An independent determination of the Hubble constant based on the tip of the red giant branch.  The Astrophysical Journal, 882(1):34, 2019.
  * [46] Raul Jimenez and Abraham Loeb.  Constraining cosmological parameters based on relative galaxy ages.  The Astrophysical Journal, 573(1):37, 2002.
  * [47] Juan Magana, Mario H Amante, Miguel A Garcia-Aspeitia, and V Motta.  The Cardassian expansion revisited: constraints from updated Hubble parameter measurements and type Ia supernova data.  Monthly Notices of the Royal Astronomical Society, 476(1):1036–1049, 2018.
  * [48] Yuting Wang, Gong-Bo Zhao, Chia-Hsun Chuang, Ashley J Ross, Will J Percival, Héctor Gil-Marín, Antonio J Cuesta, Francisco-Shu Kitaura, Sergio Rodriguez-Torres, Joel R Brownstein, et al.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: tomographic BAO analysis of DR12 combined sample in configuration space.  Monthly Notices of the Royal Astronomical Society, 469(3):3762–3774, 2017.
  * [49] Celia Escamilla-Rivera and Julio C. Fabris.  Nonparametric reconstruction of the Om diagnostic to test $\Lambda$CDM.  Galaxies, 4(4):76, 2016.
  * [50] Guido Risaliti and Elisabeta Lusso.  Cosmological constraints from the hubble diagram of quasars at high redshifts.  Nature Astronomy, 3(3):272–277, 2019.
  * [51] Guido Risaliti and Elisabeta Lusso.  A hubble diagram for quasars.  The Astrophysical Journal, 815(1):33, 2015.
  * [52] Marek Demianski, Ester Piedipalumbo, Disha Sawant, and Lorenzo Amati.  Cosmology with gamma-ray bursts-II. Cosmography challenges and cosmological scenarios for the accelerated Universe.  Astronomy & Astrophysics, 598:A113, 2017.
  * [53] Lorenzo Amati and MASSIMO DELLA VALLE.  Measuring cosmological parameters with gamma ray bursts.  International Journal of Modern Physics D, 22(14):1330028, 2013\.
  * [54] Shadab Alam, Metin Ata, Stephen Bailey, Florian Beutler, Dmitry Bizyaev, Jonathan A Blazek, Adam S Bolton, Joel R Brownstein, Angela Burden, Chia-Hsun Chuang, et al.  The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample.  Monthly Notices of the Royal Astronomical Society, 470(3):2617–2652, 2017.
  * [55] Hélion Du Mas Des Bourboux, Jean-Marc Le Goff, Michael Blomqvist, Julien Guy, James Rich, Christophe Yèche, Julian E Bautista, Étienne Burtin, Kyle S Dawson, Daniel J Eisenstein, et al.  Baryon acoustic oscillations from the complete SDSS-III Ly$\alpha$-quasar cross-correlation function at z= 2.4.  Astronomy & Astrophysics, 608:A130, 2017.
  * [56] Gong-Bo Zhao, Yuting Wang, Shun Saito, Héctor Gil-Marín, Will J Percival, Dandan Wang, Chia-Hsun Chuang, Rossana Ruggeri, Eva-Maria Mueller, Fangzhou Zhu, et al.  The clustering of the SDSS-IV extended Baryon Oscillation Spectroscopic Survey DR14 quasar sample: a tomographic measurement of cosmic structure growth and expansion rate based on optimal redshift weights.  Monthly Notices of the Royal Astronomical Society, 482(3):3497–3513, 2019.
  * [57] Adam G Riess, Steven A Rodney, Daniel M Scolnic, Daniel L Shafer, Louis-Gregory Strolger, Henry C Ferguson, Marc Postman, Or Graur, Dan Maoz, Saurabh W Jha, et al.  Type Ia supernova distances at redshift$>1.5$ from the Hubble Space Telescope multi-cycle treasury programs: the early expansion rate.  The Astrophysical Journal, 853(2):126, 2018.
  * [58] David Camarena and Valerio Marra.  Impact of the cosmic variance on H0 on cosmological analyses.  Physical Review D, 98(2):023537, 2018.
  * [59] Miguel Zumalacárregui.  Gravity in the era of equality: towards solutions to the hubble problem without fine-tuned initial conditions.  Physical Review D, 102(2):023523, 2020.
  * [60] Emilio Bellini, Antonio J Cuesta, Raul Jimenez, and Licia Verde.  Constraints on deviations from $\Lambda$CDM within Horndeski gravity.  Journal of Cosmology and Astroparticle Physics, 2016(02):053, 2016\.
  * [61] Carlos García-García, Emilio Bellini, Pedro G Ferreira, Dina Traykova, and Miguel Zumalacárregui.  Theoretical priors in scalar-tensor cosmologies: Thawing quintessence.  Physical Review D, 101(6):063508, 2020.
  * [62] Bilguun Bayarsaikhan, Seoktae Koh, Enkhbat Tsedenbaljir, and Gansukh Tumurtushaa.  Constraints on dark energy models from the horndeski theory.  Journal of Cosmology and Astroparticle Physics, 2020(11):057–057, Nov 2020.
  * [63] Noemi Frusciante, Ryotaro Kase, Nelson J. Nunes, and Shinji Tsujikawa.  Most general cubic-order Horndeski Lagrangian allowing for scaling solutions and the application to dark energy.  Phys. Rev. D, 98(12):123517, 2018.
  * [64] Thomas Helpin and Mikhail S. Volkov.  Varying the horndeski lagrangian within the palatini approach.  Journal of Cosmology and Astroparticle Physics, 2020(01):044–044, Jan 2020.
  * [65] Noemi Frusciante, Georgios Papadomanolakis, Simone Peirone, and Alessandra Silvestri.  The role of the tachyonic instability in Horndeski gravity.  JCAP, 02:029, 2019.
  * [66] Nathan Chow and Justin Khoury.  Galileon cosmology.  Physical Review D, 80(2), Jul 2009.
  * [67] Antonio De Felice and Shinji Tsujikawa.  Generalized Galileon cosmology.  Phys. Rev. D, 84:124029, 2011.
  * [68] Joe Kennedy, Lucas Lombriser, and Andy Taylor.  Reconstructing Horndeski theories from phenomenological modified gravity and dark energy models on cosmological scales.  Phys. Rev. D, 98(4):044051, 2018.
  * [69] C. Deffayet, G. Esposito-Farèse, and A. Vikman.  Covariant galileon.  Physical Review D, 79(8), Apr 2009.
  * [70] Alexandre Barreira, Baojiu Li, Carlton M Baugh, and Silvia Pascoli.  Spherical collapse in galileon gravity: fifth force solutions, halo mass function and halo bias.  Journal of Cosmology and Astroparticle Physics, 2013(11):056–056, Nov 2013.
  * [71] B et al Abbott.  Gw170817: Observation of gravitational waves from a binary neutron star inspiral.  Physical Review Letters, 119(16), Oct 2017.
  * [72] D. A. Coulter, R. J. Foley, C. D. Kilpatrick, M. R. Drout, A. L. Piro, B. J. Shappee, M. R. Siebert, J. D. Simon, N. Ulloa, D. Kasen, and et al.  Swope Supernova Survey 2017a (SSS17a), the optical counterpart to a gravitational wave source.  Science, 358(6370):1556–1558, Oct 2017.
  * [73] T. et al Baker.  Strong constraints on cosmological gravity from gw170817 and grb 170817a.  Physical Review Letters, 119(25), Dec 2017.
  * [74] Paolo Creminelli and Filippo Vernizzi.  Dark Energy after GW170817 and GRB170817A.  Physical Review Letters, 119(25), Dec 2017.
  * [75] Jose María Ezquiaga and Miguel Zumalacárregui.  Dark Energy After GW170817: Dead Ends and the Road Ahead.  Physical Review Letters, 119(25), Dec 2017.
  * [76] Cedric Deffayet, Oriol Pujolas, Ignacy Sawicki, and Alexander Vikman.  Imperfect Dark Energy from Kinetic Gravity Braiding.  JCAP, 10:026, 2010.
  * [77] Carlos García-García, Emilio Bellini, Pedro G. Ferreira, Dina Traykova, and Miguel Zumalacárregui.  Theoretical priors in scalar-tensor cosmologies: Thawing quintessence.  Physical Review D, 101(6), Mar 2020.
  * [78] Eleonora Di Valentino, Luis A Anchordoqui, Yacine Ali-Haimoud, Luca Amendola, Nikki Arendse, Marika Asgari, Mario Ballardini, Elia Battistelli, Micol Benetti, Simon Birrer, et al.  Cosmology Intertwined III: $f\sigma{{}_{8}}$ and $S_{8}$.  arXiv preprint arXiv:2008.11285, 2020.