Taxonomy of Dark Energy Models

The accelerated expansion of the Universe is one of the main discoveries of the past decades, indicating the presence of an unknown component: the dark energy. Evidence of its presence is being gathered by a succession of observational experiments with increasing precision in its measurements. However, the most accepted model for explaining the dynamic of our Universe, the so-called Lambda cold dark matter, face several problems related to the nature of such energy component. This has lead to a growing exploration of alternative models attempting to solve those drawbacks. In this review, we briefly summarize the characteristics of a (non-exhaustive) list of dark energy models as well as some of the most used cosmological samples. Next, we discuss how to constrain each model's parameters using observational data. Finally, we summarize the status of dark energy modeling.


I. INTRODUCTION
One of the most important discoveries in modern cosmology since the Universe expansion by George Lemaître [1] and Edwin Hubble [2], and the detection of the cosmic microwave background (CMB) [3,4] is the accelerated expansion at late times. This cosmic acceleration was detected at the end of the 90s, when two independent collaborations where observing high redshift type Ia supernovae (SNIa) to measure the curvature and deceleration parameter of the Universe [5,6]; both groups established a cosmological constant model of the Universe, with Ω m ≈ 0.3, and Ω Λ ≈ 0.7. Later, this result was confirmed by different cosmological probes, for example the CMB [7,8], baryon acoustic oscillations [e.g., 9,10], and gravitational lensing, strong [e.g, 11] and weak [e.g., 12].
In the last decades, substantial effort has been made to understand the nature of DE [e.g. see 13, for a recent review] with the ΛCDM model being the cornerstone of modern cosmology and the simplest one, composed by two dominant components, known as cold dark matter (DM) and cosmological constant, and three subdominant components identified as baryons, photons and neutrinos. The ΛCDM model is not only effective at background level (isotropic and homogeneous Universe), but also robust at linear perturbations, having accurate predictions for the matter power spectrum and the small differences for photons temperature observed in the CMB [14,15]. As mentioned previously, the ΛCDM model is composed by a ∼ 32% of cold dark matter which is essential for the structure formation, being the most suitable explanation for the rotation curves of galaxies [see for instance [16][17][18]. It is assumed that dark matter is a non relativistic particle and the preferred candidates are the particles emerging from the Supersymmetric theory [19,20]. However, the lack of evidence of these particles strengthen the proposition of other candidates as Axions, Ultra light scalar field dark matter, among others [21][22][23]. Despite its remarkable achievements of the λCDM model, it has important afflictions when describing the nature of the CC through the idea of a quantum vacuum fluctuations [24,25]. This idea generates, from the theoretical point of view, an error of ∼ 10 120 orders of magnitude and it is known as the fine tuning problem [26]. In addition, the CC has the coincidence problem [26], i.e. why the Universe accelerate at late epochs and not before of after? On top of that, recent observations from Planck and Supernovae Type Ia differ in their values for H 0 1 , introducing tension between different observations at different redshifts, [29][30][31][32][33]. The community attribute the problem to the ΛCDM model and in particular to the CC, therefore new approaches are proposed to alleviate the discrepancy among observations [see also [34][35][36].
In the context of these tensions and the problems with CC, a plethora of alternative dark energy models have been explored. Our aim is to review a set of those models consisting of fluids that can be described by a variety of formulations. Some of them can be described by different fluids and their particles that compose them, like scalar fields, while others that do not need extra fluids require modifications to general theory of relativity [GTR, [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51] (for recent overviews of GTR, see, e.g., [52][53][54] and references therein). Specifically, the first category contains models avoiding the idea of associating the Universe acceleration with quantum vacuum fluctuations (like in CC) and, thus assume a fluid expression manageable by the quantum field theory. The second one, consist of models introducing extensions to GTR, generating a Universe acceleration without the addition of extra fluids. However, the current overabundance of models is overwhelming, posing difficulties to decide which is the best candidate to compete against the ΛCDM model.
Recent years have also witnessed the increase of observational surveys designed to obtain precise measurements aimed to probe DE nature. With this improvement in the instrument sensitivity came the mechanism to discriminate between models of DE, i.e. tensions between different probes could lead to rule out some of them. In this context, we consider widely used samples (such as SNIa [5], CMB [55], baryon acoustic oscillations) as well as other recent compilations (strong gravitational lens systems [56][57][58], starburst galaxies [59][60][61][62][63][64], and observational Hubble parameters [58,65].
The outline of the review is as follows: Section II summarizes the basic equations, Section III presents the cosmological samples that we use to constraint the different theoretical models, in Section IV we describe the DE models, together with the constraints of the free parameters of each models. Finally, in Section V we give a discussion and conclusion of the models mentioned, discussing the promising models and what is the contribution to the understanding of the Universe acceleration.

II. BASIC EQUATIONS FOR THE BACKGROUND COSMOLOGY
Modern cosmology is based on the general theory of relativity whose master equation is the field equation given by where G µν ≡ R µν − g µν R/2 is known as the Einstein tensor, composed by the Ricci tensor (R µν ), the Ricci scalar (R) and the metric tensor (g µν ). The right side of Eq. (1) shows the Newton gravitational constant (G) and the energy momentum tensor. Hereafter we will use natural units (c = = 1), unless we explicitly mention the opposite. Hereafter, we focus our attention in the background cosmology considering a homogeneous and isotropic Universe. In this vein, we use the standard line element of Friedmann-Lemaitre-Robertson-Walker (FLRW) with flat geometry k = 0, as ds 2 = −dt 2 + a(t) 2 (dr 2 + r 2 dΩ 2 ), (2) where dΩ 2 ≡ dθ 2 + sin 2 θdϕ 2 and a(t) is the scale factor. The energy-momentum tensor is written as always as being p, ρ, u µ and g µν the pressure, density, four-velocity of the fluid and the metric tensor, respectively. The covariant derivative of the energy momentum tensor ∇ µ T µν = 0, generates the continuity equation, given by where the sum is over all the species and dots stand for derivatives with respect to time.
Through the Einstein field equations, we arrive to the Friedmann equation, which is a first order non-linear differential equation composed by the scale factor and the densities of the species. Therefore, the equation takes the form where H is known as the Hubble parameter, which indicates the Universe expansion rate. Hence, it is possible to define the dimensionless Friedmann equation in the form where Ω(z) ≡ 8πGρ(z)/3H 2 0 is known as the density parameter, a function of the redshift (z), and the sum runs for the components of the Universe used in the model 2 , H 0 is the Hubble constant 3 . With the previous equations, it is possible to define the flat condition as notice that the density parameters are evaluated at z = 0 and they are denoted with the label '0'. The comoving distance from the observer to redshift z is given by (units of c are recovered in the following equations) Therefore, we define the luminosity distance, denoted as d L (z), as where c is the speed of light. The angular diameter distance is related to the comoving distance as In addition, the angular diameter distance between two objects at redshifts z 1 and z 2 (z 1 < z 2 ) is given by III. COSMOLOGICAL SAMPLES Observational samples are required not only to investigate the consistency of any theoretical model prediction but also to discern between different models. In this section we review the most common samples used as "standard candles" to probe those models. Given a cosmological model, the astronomical distance to an object is generally measured using the redshift to the object and the distance-redshift relationship provided by the model (see equation 9). However, when the goal is to infer the cosmological model, an independent distance measurement is needed. The methods to obtain these independent measurements are based on empirical relationships and, thus, the astronomical objects following those relations are called standard candles.
Knowing the intrinsic flux of a standard candle 4 , we can obtain a relation between the apparent magnitude (m, observed) and the absolute magnitude (M , acquired from an empirical relation) In the following sections we include a general description of some of the most used cosmological samples with emphasis in the quantification of the goodness-of-fit (or figure-of-merit) by defining a function (χ 2 ) associated to the sample errors or covariances, that is applied to investigate the cosmological models presented in Section IV.

A. Type Ia Supernovae
Type Ia supernova is believed to originate by a white dwarf accreating matter from a companion star. When the white dwarf exceeds the Chandrasekhar mass limit (∼ 1.4M , where M is 1 solar mass), there is a collapse and subsequent explosion. As all the SNIa have roughly the same luminosity 5 , they can be used as standard candles.
Samples of SNIa (e.g. [31,48,[66][67][68]) provide distance modulus measurements (see Eq. 12). As the measurements in this kind of samples are correlated, it is convenient to build the chi square function as where e = ∆1 T · Cov −1 P · ∆1, and ∆μ is the vector of residuals between the theoretical distance modulus and the observed one, ∆1 = (1, 1, . . . , 1) T , Cov P is the covariance matrix formed by adding the systematic and statistic uncertainties, i.e. Cov P = Cov P,sys + Cov P,stat . The super-index T on the above expressions denotes the transpose of the vectors.
The theoretical distance modulus is estimated by where M is a nuisance parameter which has been marginalized in (13).

B. Baryon Acoustic Oscillations
Another way to establish a constraint of model parameters is through the standard rules known as Baryon Acoustic Oscillations (BAO). These are primordial signatures in the matter power spectrum produced by the interaction between baryons and photons in a hot plasma in the pre-recombination epoch.
The theoretical BAO angular scale (θ th ) is estimated as where d L and D A are written in Eqs. (9) and (10), respectively. The parameter r drag indicates the sound horizon at baryon drag epoch. The comoving sound horizon, r s (z), is defined as where the sound speed c s (z) = 1/ 3 1 +R b / (1 + z) , withR b = 31500 Ω b h 2 (T CM B /2.7K) −4 , and T CM B is the CMB temperature. The redshift z drag at the baryon drag epoch is well fitted with the formula proposed by Eisenstein and Hu [69] where where Ω m0 and Ω b0 are the matter component (dark matter plus baryons) and baryon component at z = 0 respectively. For this work, we set the r drag = 147.21 ± 0.23 obtained by Planck collaboration [15]. Notice that, as BAO data points are estimated using r drag , which depends on the cosmological model, they could be considered as biased.
The most recent compilation of transversal BAO measurements θ BAO (z) is presented in [70]. A total of 15 measurements [71][72][73][74][75] were obtained using the data realeases (DR), DR7, DR10, DR11, DR12, DR12Q (quasars), of Sloan Digital Sky Survey (SDSS) [76]. As transversal angular BAO points are considered uncorrelated, the chi square function is built as where θ i BAO ± σ θ i BAO is the BAO angular scale, N is the number of data and its uncertainty at 68% measured at z i . It is worth to mention that there is a sample of 6 correlated data points, with their associated covariance matrix, collected in [77] and measured by [78][79][80]. In this case, the chi square function is where X is the difference between the theoretical and observational quantities of d A (z drag )/D V (z i ) measured at the redshift z i , and Cov −1 is the inverse covariance matrix (see [77] for details), the dilation scale (D V ) is defined as [81] where d A (z) = (1 + z)D A (z) is the comoving angular-diameter distance.

C. Cosmic Microwave Background Radiation
In the early Universe baryons and photons are coupled, leading to coherent oscillations that are observed in the power spectrum 6 of the CMB (e.g. WMAP [8,82], Planck [83,84]). This is a powerful probe due to its ability to estimate the cosmological parameters with high precision [85]. The information of the CMB acoustic peaks can be condensed in three quantities, their distance posteriors: the acoustic scale (l A ), the shift parameter (R), and the decoupling redshift (z * ). Several authors have proved that these quantities are almost independent of the DE model considered and, thus they can be used to test the parameters of alternative cosmologies [86][87][88]. The acoustic scale is defined as where r s is the sound horizon (Eq. (17)) at the redshift of decoupling z * given by Hu and Sugiyama [89], where The shift parameter is defined as [90] where Ω m0 include the baryon and DM components. Thus, the χ 2 for the CMB data is constructed as where Cov −1 CMB is the inverse covariance matrix of the distance posteriors and the superscripts th and obs refer to the theoretical and observational estimations respectively. To infer the parameters of the alternative cosmologies we employ the distance posteriors of WMAP [82] and Planck [91].

D. Observational Hubble Parameter
The Hubble parameter is estimated mostly by using the differential age (DA, [65]) methodology and from BAO measurements. The former method consists of measuring the age between pairs of passive evolving galaxies (dubbed cosmic chronometers) with similar metallicity and separated by a small redshift interval 7 with redshift z 2.0. Thus, the expansion rate is written as where dz is measured with high accuracy 8 . The OHD from the DA method are considered cosmological independent measurements. On the other hand, the OHD from BAO surveys are non-homogeneous since they depend on the cosmological model selected. By taking a unique value for r drg in these data, an OHD homogeneous sample can be obtained [see 68, for further details]. The observational Hubble parameter data (OHD) represents the most direct way to constrain the parameter space to mimic the observational expansion rate, the chi square function can be expressed as where H th (z i ) is the theoretical estimate using (6) or a generalization, H obs (z i ) ± σ i obs is the observational Hubble parameter (from DA, (non)-homogeneous BAO points, or the joint of them) with its uncertainty at the redshift z i , and N is the number of points used. , where δθ E and δσ are the uncertainties of the Einstein radius and the observed line-of-sight (1D) velocity dispersion, respectively. The theoretical counterpart is estimated by the ratio A corrective parameter f is often introduced in Eq. (32) to take into account possible systematic differences among systems (e.g. elliptical instead of spherical profile for the lens halo, line-of-sight stellar velocity dispersion as opposed to the dark matter halo velocity dispersion, steeper mass distribution profile, see for example [99,100]).

F. Ionized Gas in Starburst Galaxies
Authors [59-63, and references therein] argued that the correlation between the measured luminosity L and the inferred velocity dispersion σ of the ionized gas (e.g. Hβ, Hα, [OIII] emission lines) in extreme starburst galaxies (i.e. containing a population of O and/or B stars) may be used as a cosmological tracer to constrain cosmological model parameters. Compilations provide apparent magnitude, emission line luminosity and velocity dispersion (e.g. [63,64]) and the chi square function is estimated as where In the above expressions, µ i obs ± σ i obs is the observed distance modulus with its uncertainty at redshift z i . The theoretical estimate at the redshift z is obtained by using (9) and where µ 0 is a nuisance parameter which has been marginalized.

G. Joint Analysis
Testing the consistency of a given cosmological model requires a range of observational samples with complementary sensitivity to the cosmological parameters. Constraint on those parameters is usually achieved by combining several samples also known as joint analysis.
For instance, a Bayesian Markov Chain Monte Carlo (MCMC) analysis is able to constrain the phase-space parameter Θ of a cosmological model given a number of cosmological samples. In general, the procedure consists of using emcee Python package [101] for two phases: the burn-in and the MCMC. The first is performed with a certain number of steps to achieve the convergence of the chains according to the Gelman-Rubin criterion [102]. The second phase is performed with an appropriate number of steps for sampling the confidence regions. Additionally, for each model, the priors (flat or Gaussian) of the parameters are chosen according to values provided in the literature. For the joint analysis, the figure-of-merit to be optimized is given by where the χ data represents the name of the different samples. In general, the joint analysis is calculated using the combination of at least three data samples, but ideally it should contain all of them.

IV. TAXONOMY OF DARK ENERGY MODELS
This section is dedicated to describe the different constrictions through the samples mentioned previously for the different DE models studied in literature. We divide our study in DE models linked to a fluid with the capability of accelerating the Universe and models in which the Einstein field equations of General Theory of Relativity are modified.

A. Accelerating Universe Fluids
In this subsection we summarize all those models that involve a fluid enabling a late acceleration without modifications to GTR.

The ΛCDM Model
The ΛCDM is the consensus model dominated by a cold dark matter and a cosmological constant component with subdominant species of baryons and relativistic particles (photons and neutrinos), being not only the most favoured by diverse observations but also the simplest. The dimensionless Friedmann function is given by Eq. (6), which reads The flatness condition is satisfied and written in the form Therefore, the CC density parameter can be written in terms of matter, while radiation takes the form where N ef f = 3.04 is the standard number of relativistic species [104] and h = H 0 /100 km s −1 Mpc −1 , where H 0 = 67.66 ± 0.42km s −1 Mpc −1 with Planck [55], while H 0 = 73.2 ± 1.3km s −1 Mpc −1 with Riess [105], presenting a tension between the observations and known as H 0 tension. Regarding the matter density parameter, the value is constrained as Ω m0 = 0.3111 ± 0.0056, using Planck satellite [55], which is a combination of baryonic and dark matter. Despite the model achievements, ΛCDM is afflicted with several problems, like the nature of CC [24,25], the σ 8 tension [29,30] and the H 0 tension [31][32][33].

The ωCDM Model
This model is the simplest extension of the CC. The dark energy has a constant equation of state (EoS) but it deviates from w = −1, and should satisfy ω < −1/3 to obtain an accelerated Universe. The equation E(z) can be written as: The ωCDM constrains are obtained assuming flat priors on the parameters. Table I presents the mean values for the ωCDM parameters using using independently OHD (31 data from cosmic chronometers), CMB (Planck) and SNIa (Pantheon) and the joint of them. Fig. 1 shows the mean value curve of the H(z) function (top panel) for the wCDM model using these data. The bottom panel shows the constraint contours at 1σ, 2σ, and 3σ confidence levels. Notice that, although SNIa data is not able to constrain the h parameter, the three different samples provide consistent constraints on the Ω m0 − ω space. Indeed, the joint analysis provides stringent constraints which are consistent with those of the ΛCDM model.

Dark Energy Parameterizations
The natural alternatives to the ωCDM is to consider DE varies with redshift through a parameterization w(z). These functions are proposed phenomenologically to mimic the behaviour of the CC at late times. In the following we present some of these models for a Universe containing dark and baryonic matter, radiation, and dark energy.
• The Chevallier-Polarski-Linder parametrization [CPL, 106,107].-An approach to study dynamical DE models is through a parametrization of its EoS. The dimensionless Hubble parameter E(z) for this Universe is given by We compute Ω r0 in the same form as is given by Eq. (44).
The density parameter for DE is written as Ω de = 1 − Ω m0 − Ω r0 , and the function f de (z) depends on w(z) as where ρ de (z) is the energy density of DE at redshift z, and ρ de (0) is its present value. One of the most popular parameterization is proposed by [106,108], and reads as where ω 0 is the EoS at redshift z = 0 and ω 1 = dw/dz| z=0 . Although this function is widely used it has a divergence problem when z = −1. The function f de (z) for the CPL parametrization is The h, Ω m0 , ω 0 , ω 1 parameters are constrained using the OHD from cosmic chronometers [109]. Figure 2 shows the reconstruction of H(z) using the best fit of the MCMC analysis: h = 0.73 +0.10 −0.08 , Ω m0 = 0.29 +0.09 −0.08 , ω 0 = −1.51 +0. 80 −0.91 , and ω 1 = −0.20 +2.38 −2.53 . The confidence contours of the parameters at 1σ, 2σ, and 3σ are also shown.
• The Jassal-Bagla-Padmanabhan (JBP) parametrization.-Jassal et al. [110] proposed the following ansatz to parametrize the dark energy EoS where ω 0 is the EoS at redshift z = 0 and • The Barbosa-Alcaniz (BA) parametrization.-Barboza and Alcaniz [111] considered a EoS given by: This ansatz behaves linearly at low redshifts as w 0 + w 1 , and w → w 0 + w 1 z when z → ∞. In addition, w(z) is well-behaved for all epochs of the Universe. For instance, the DE dynamics in the future, at z = −1, can be investigated without dealing with a divergence. Solving the integral in Eq. (47) and using Eq. (52) results in: • Feng-Shen-Li-Li [FSLL,112] parametrizations.-The authors suggested two dark energy EoS given by: Both functions have the advantage of being divergence-free throughout the entire cosmic evolution, even at z = −1. At low redshifts, w(z) behaves as w 0 + w 1 z and w 0 + w 1 z 2 for FSLLI and FSLLII, respectively. In addition, when z → ∞, the EoS has the same value (w 0 ) as the present epoch for FSLLI and w 0 +w 1 for FSLLII.
• Sendra-Lazkoz [SL,113] introduced new polynomial parameterizations to reduce the parameter correlation, so they can be better constrained by the observations at low redshifts. One of these parameterizations is given by: where the constants are defined as c 1 = (16w 0 − 9w 0.5 + 7)/4, and c 2 = −3w 0 + (9w 0.5 − 3)/4, and w 0.5 is the value of the EoS at z = 0.5. This w(z) function is well-behaved at higher redshifts as (−1 − 8w 0 + 9w 0.5 )/2. The substitution of Eq. (57) into Eq. (47) results: Notice that, although DE parameterizations are common and they could solve the coincidence problem, there is not a unique way to choose the form of the function. Furthermore, in many cases there are not strong arguments to justify the functional form by an association with a first-principles theory of quantum fields or gravity. A different approach, which is model-independent, consist of, for example investigating the cosmographic parameters that characterize the kinematics of the cosmic expansion [e.g., [114][115][116][117][118]. Some authors have used the Hubble parameter, the deceleration parameter (q(a) = −äa/ȧ 2 ), or even higher order derivatives of the scale factor a, such as Jerk and Snap [e.g., 119,120]. By estimating these cosmographic parameters using cosmological data, it is possible to associate its features to a given DE model [see 119,[121][122][123][124][125].
The cosmological constrains for the aforementioned models are obtained assuming flat priors on the DE parameters and a Gaussian prior on h.

Chaplygin-Like Fluid
One point of view for studying the DE and DM problems is through the unified dark fluids approach, which is known as Chaplygin gas [see, for instance, [129][130][131][132][133]. An example of this is the well-known generalized Chaplygin gas [134,135] described by the EoS p = −Aρ −α where A and α are constants (the case α = 1 is the original model proposed by S. Chaplygin [136]). This fluid behaves as DM at early epoch and DE at late times and may have its origen from the Nambu-Goto d-brane action. Although this interesting formulation reproduces the accelerated expansion of the Universe, it presents flaws to describe the CMB anisotropies [137]. In this context, an alternative to Chaplygin gas was proposed by Hova and Yang [138], dubbed generalized Chaplygin gas-like, with EoS being sinc(x) ≡ sin(x)/x and ρ df the dark fluid density, which plays the role of the mixture of DE and DM densities. In this case, µ is a dimensionless parameter constrained as µ 0.688 to be consistent with the stellar age bound 9 and ρ df 0 is the present energy density of this fluid, constrained in terms of the density parameter as Ω df 0 ∼ 0.96 in Hova and Yang [138]. It behaves as a CC in the late stage of the universe and as DM at the matter domination epoch. The evolution of the EoS of the dark fluid is given by where ξ(z) ≡ arctan[(z + 1) −3 tan λ] and λ ≡ µπ/2. To explore the universe dynamics in this context, we consider a general FLRW metric including baryonic and radiation components, hence we write the Friedmann and acceleration equations as From Eq. (6) we have [138,139] where Ω df 0 ≡ 8πGρ df 0 /3H 2 0 is the density parameter associated with the Chaplygin gas-like fluid, Ω i0 and ω i are the density parameters and the EoS for baryonic matter and radiation (according to Eq. 44), Ω k ≡ −k/H 2 0 is the curvature density parameter and H 0 = h × 100 km s −1 Mpc −1 . In addition, from (7) we have the constraint Figure 4 shows the best fit curve (top panel) for the Chaplygin-like gas when the curvature term is neglected using the OHD, SNIa and OHD+SNIa (Joint) samples. Additionally, 2D contours at 1σ, 2σ, and 3σ confidence level (CL) and 1D posterior distribution of the free parameters are displayed for each sample. Table III presents the best fit values and their uncertainties at 1σ for the model free parameters.

Viscous Model
The accelerated expansion of the Universe may be also described by considering dissipative effects in the Universe components, mainly in the matter component. The bulk viscosity coefficient, which satisfies the cosmological principle, is introduced in the energy-momentum tensor, Eq. 3, as an effective pressure as p →p = p + Π where Π = −3ξH based on the Eckart formalism. Under this argument, several models for ξ have been addressed such as: Probably this model, where ρ m is the energy density of dust matter and ξ 0 , s are constants, is the simplest one that successfully reproduce the late accelerated stage of the Universe. Some studies that consider a single fluid in the Universe are presented in [140] (see for example [141] for a case in a causal theory). Additionally, there are other works that include several components such as radiation and DE [67].
• ξ = ξ(z). In spite of the success of the previous model at late epochs of the Universe, it has problems in early epochs because ξ diverges. This motivates the use of alternative viscosity models such as those proposed by [142], in particular polynomial forms of the redshift.
• ξ = A cosh(bE −n ) and ξ = A tanh(bE −n ). Alternatively, more complex models are investigated in [142] by proposing the viscosity as a hyperbolic function of the dimensionless Hubble parameter E. When a single dust matter fluid is contained in the Universe, the dimensionless Hubble parameter is obtained using Eq. (6). Then, we obtain the following system to be solved where λ(z) = ξ(z)κ 2 /3H 2 0 . For λ(z) = λ 0 + λ 1 (1 + z) n , we obtain, where Fig . 5 shows the best fit curve of H(z) (top panel) using non-homogeneous OHD+SNIa data. Two-dimensional contours at at 1σ, 2σ, 3σ CL of the free model parameters are also displayed at the bottom panel for OHD, SNIa, OHD+SNIa data. Additionally, we include the best fit curves and 2D contours when λ(z) = 1/3 tanh(bE −n ) and   A generalized form of the previous model (Eq. (65)) considers one additional fluid to the dust matter component. Authors in [143] analyze the case that includes a DE component with EoS w = −1 to describe the late time stage of the Universe. Notice that the radiation component at this time can be considered negligible. In this case, the Hubble parameter is given by where E(0) = Ω(0) = 1, Ω(z) = Ω m0 (1 + z) 3 + Ω de0 , and 2 F 1 is the hypergeometric function. Based on the results obtained in [142], Eq. (67) is obtained assuming n = −2. It is straightforward that the case for a constant viscosity coefficient is obtained when λ 1 = 0. Figure 6 displays the best fit curves (top panel) of Eq. (67)

Interacting Viscous Models
A generalized case of the viscous models presented in the previous section is to consider a flat FLRW Universe which contains a non-perfect fluid as dust matter (dm) component that interacts with a perfect fluid as the DE where ρ r , ρ dm , and ρ de are the relativistic species, dust matter and dark energy densities, respectively. Notice that the DE component behaves as CC when γ de = 0. In particular, the typical ansatz for the viscosity coefficient is considered and given by where ρ dm0 is the dm density at present epoch and ξ is a free parameter with units of [ξ] =[eV]. It is convenient to use the dimensionless parameter of ξ defined as ξ 0 = √ 3ξ/κρ dm0 . Additionally, the interacting term Q is consider to be [144] where β is a free parameter. It is straightforward that an Universe with only viscosity effects is obtained when β = 0. Figure 7 shows the best fit curve (top panel) to OHD data for interacting viscous model (ξ 0 = 0, β = 0), viscous model (β = 0), interacting model (ξ 0 = 0) and ΛCDM, respectively. 2D contours at 1σ, 2σ, 3σ CL and 1D posterior distributions of the free model parameters are presented at the bottom panel. Authors in [67] found that the energy density dynamics of the mentioned models are similar to the evolution of ΛCDM. Table VII reports the best fit values and their uncertainties at 1σ for the free parameters of IVM, IM, VM and LCDM models.

Phenomenological Emergent Dark Energy Model
The phenomenological emergent dark energy model (PEDE) was proposed by [37] and assume that the DE is negligible at early times, emerging at late times. These kind of models are known as emergent and contribute to  We consider a FLRW metric which contains matter (m, dark matter plus baryons), radiation (r), and PEDE. The dynamics of this Universe is described by the Friedmann equation (5) and the continuity equation for each component By solving Eqs. (75a), (75b), (75c) we can rewrite the Eq. (6) in terms of the density parameters and redshift, as where Ω DE (z) = Ω DE0 f (z), where Ω DE follows Eq. (7). Notice that [37] propose a phenomenological functional form for f (z), described by Eq. (47) and hence Ω DE (z) as 10 where Ω DE → 0 at z → ∞ and Ω DE → 1.4 at z → −1. Notice that Therefore, the dimensionless Friedmann equation results as

Generalized Emergent Dark Energy
Recently, [38] proposed a generalisation for the PEDE model, also known as Generalized Emergent Dark Energy Model (GEDE) model, by introducing where z t is a transition redshift, Ω DE (z t ) = Ω m0 (1 + z t ) 3 , ∆ is an appropriate dimensionless non-negative free parameter with the characteristic that if ∆ = 0 the ΛCDM model is recovered, and when ∆ = 1 and z t = 0 the previously PEDE model is obtained. As z t can be related to Ω m0 and ∆, then z t is not a free parameter. Notice that the DE density parameter is given by 1+tanh(∆log 10 (1+zt)) .
The GEDE Friedmann equation is written as The results obtained from the MCMC analysis using the same data as PEDE model are shown in Fig. 9, presenting the best fit curve for H(z) confronting with the OHD data and the constraints for Ω (0) m and ∆ which is the free parameter for GEDE.

B. Modifications to General Theory of Relativity
In this subsection, we present models that modify the GTR in order to obtain a late Universe acceleration.

Constant Brane Tension
Brane world models are inspired by the seminal papers of [146,147] in which they assume a four dimensional manifold called the brane immersed in a five dimensional Anti-d'Sitter space time called the bulk. The mentioned configuration is a via to understand the hierarchy problem but could be also extended to describe the cosmology. The main parameter of the theory is called the brane tension, which differentiate between the high and low energy physics involved and becomes a free parameter that needs to be constrained by different cosmological samples, for this model in particular, the brane tension is constant, hence, we call it a Constant Brane Tension (CBT) model.
First of all, we introduce the Einstein's field equation projected onto the brane where T µν is defined in Eq.(3) of the matter trapped in the brane, G µν is the classical Einstein's tensor described by (1) and the rest of terms in the right and left sides of this equation are explicitly given by: Here G N is the Newton's gravitational constant, λ is the previously mentioned brane tension, κ (4) and κ (5) are the fourand five-dimensional coupling constants of gravity, respectively. The tensor Π µν represents the quadratic corrections on the brane generated by the energy-momentum tensor, F µν gives the contributions of the energy-momentum tensor in the bulk, which is projected onto the brane through the unit normal vector n A . The tensor ξ µν provides the contribution of the five-dimensional Weyl's tensor projected onto the brane manifold [148] 11 . It is worth to note that non-local corrections are negligible in cosmological cases [40], under the assumption of a AdS (5) bulk.
To derive the Friedmann equations under the modified field equations, we consider an homogeneous and isotropic Universe in which a line element is given by Eq. (2). We consider radiation and dark matter components as perfect fluids in the brane. We assume that the bulk has no matter component. Using Eqs. (83), we obtain the modified Friedmann equation: Notice that ρ i is the energy density for the radiation, dark matter and DE. It is worth to notice that the low energy regime, i.e. the canonical Friedmann equation, is recovered when ρ i /2λ → 0. Crossed terms were not used in the Friedmann equation, i.e. there is not interaction between different species. In addition, if we consider, for instance, that the bulk black hole mass vanishes, the bulk geometry reduces to AdS 5 and ρ = 0 [40,41]. Thus, the Friedmann equation can be written as: The above equation can be expressed in terms of the density parameters through the dimensionless Friedmann equation +M Ω 2 0m (1 + z) 6 + Ω 2 0r (1 + z) 8 + Ω 2 0de (1 + z) 6(1+ω de ) , where being ρ crit the Universe critical density. Notice that when M → 0, the canonical Friedmann equation with w de is recovered. If w de = ω Λ = −1, i.e. the DE is the CC, we obtain the traditional ΛCDM dynamics. At early times the brane dynamics dominate over other terms in the Universe, but is negligible at late time. Indeed, given a value for the brane tension, we can infer the limits of high and low energies in terms of the redshift: z + 1 i (λ/ρ 0i ) 1/3(1+ωi) and z + 1 i (λ/ρ 0i ) 1/3(1+ωi) respectively. For example, in matter domination epoch, the previous expressions can be rewritten as: z (λ/ρ 0m ) 1/3 − 1 and z (λ/ρ 0m ) 1/3 − 1, for high and low energy limits, respectively. Table X shows the best-fit for the different free parameters together with the estimated χ 2 . Fig. 10 shows the severe tension between the different cosmological samples (OHD, SNIa, SLS, HIIG, BAO), hence concluding that the model is not viable to replace the ΛCDM model unless an additional DE component is included, which defeat the purpose of choosing this model.

Variable Brane Tension
A natural extension to the previous model is the one called variable brane tension (BVT). The framework is the same as the on in Section IV B 1, but now an extra degree of freedom is assumed, a brane tension emerge as a function of the redshift. Naturally, the model resolve the problems associated with the presence of the brane tension at early epochs but also generates a CC with five dimensional origins. We briefly discuss the theoretical framework of a BVT model which was previously studied in [43]. We start from the BVT field equation as where where G µν is described in (1), ξ µν is a non-local Weyl tensor decomposed in its irreducibility form which also contains U is the non-local energy density, P µν is the non-local anisotropic stress tensor, u α is the four-velocity and µν ≡ g µν + u µ u ν . In addition T µν is the standard energy-momentum tensor and Π µν contains a quadratic form of the energy-momentum tensor. Notice that the corrective terms that comes from brane world is contingent to the brane tension defined by λ, which in this model is not a constant. Therefore, the low energy limit is considered when λ → ∞ recovering the traditional field equation of GR, while in the other limit λ → 0 extra terms play a preponderant role. Finally, notice that in this case we do not consider extra fields onto the bulk, neglecting the terms that come from F µν and only considering those fields living in the brane. Therefore, if we introduce the previously line element in Eq. (89) together with the perfect fluid energy-momentum tensor (Eq. (3)) we have the following Friedmann equation [43] E(z) 2 = Ω 0m (z + 1) 3 + Ω 0r (z + 1) 4 + M λ(z) [Ω 2 0m (z + 1) 6 + Ω 2 0r (z + 1) 8 ], here we have already considered matter and radiation components, where their evolution come from the conservation of the energy-momentum tensor (∇ µ T µν = 0) and their separability from the quadratic part of the field equation. The brane tension evolves homogeneous and isotropically because it only on the temporal function and can be chosen using other physical assumptions. In addition, the brane tension is not directly coupled with the continuity equation of the fluids and it is defined asλ(z) ≡ λ 0 λ(z), beingλ(z) a dimensionless function that can be selected appropriately. Moreover, M ≡ 3H 2 0 /16πGλ 0 and, under the flatness condition, we have the constriction Regarding the choice of theλ(z) function, we pick a polynomial form asλ(z) = (z + 1) n , where n is a free parameter and n ∈ R. Garcia-Aspeitia et al. [43] discuss the inspiration for this function, arguing that it could be a generalization of the Eötvös law, similar functions can be found in tracker behavior for scalar fields.
Using a joint analysis that contains OHD, CMB, BAO and SNIa observations (see Fig. 11 and Table XI), [43] found out n = 6.19. ± 0.12. Notice that the result is consistent with predictions because if we use Eq. (92), neglecting (z + 1) 8 , the term M will behave as a CC at late times but with extra dimensions origin.

Unimodular Gravity
Unimodular Gravity (UG) is a remarkable proposition to tackle the problem of the CC by limiting the metric in the following way √ −g = ξ, where ξ is a constant, restricting the field equations at only nine linear independent equations and the field equation is trace-free [44]. The possibility to integrate the line element of FLRW give us the opportunity to obtain clues about the nature of CC, tracing its presence at epochs of reionization [46]. UG can be described by the following field equation where all the tensors are the standards of GR and G is the Newton's gravitational constant.
In order to study the background cosmology, we consider an isotropic, homogeneous FLRW metric (2), the perfect fluid energy momentum tensor is written as show Eq. (3). Hence, we have [44,45] where the dots stands for time derivative. In addition, a general conservation for UG theory is now written in the form Without independently assuming the energy momentum conservation (∇ µ T µν = 0), the Eq.
where the sum is over all the species in the Universe and j ≡ ... a /aH 3 is the Jerk Parameter (JP) [118,149], well known in cosmography and proposed by [46] for the study of UG.
On the other hand, the integral-transcendent-Friedmann equation can be computed with the help of Eq. (5) and (97), obtaining the Friedmann equation as where the non-canonical extra term in Eq. (98), i.e. the UG correction to the Friedmann and acceleration equations, is defined in the form where the sum runs over the different species in the Universe and a ini is some constant initial value.  According to [46], it is plausible to consider an ansatz for the JP in terms of the redshift with the following characteristics where w is the EoS for any fluid. If we choose Ω 0i → Ω 0r as the radiation density parameter 12 and w → w r = 1/3 as the EoS of radiation, then the functional form reproduce the ΛCDM jerk parameter in all eras [46]. From the previous equation and Eq. (98) it is possible to deduce where Ω 0exs ≡ w r Ω 0r . Notice that the source of the Universe acceleration is the constant term in the previous equation (101), where we can naturally relate Ω 0Λ → Ω 0exs (z ini + 1) 4 . Since our choice in Eq. (100) depends on the EoS and the energy density parameter of the radiation, the constant inherits those terms. Fig. 12 shows the constraints for z ini and h considering the SNIa, BAO, OHD, CMB and a joint data, with the best-fit values in Table XII [see 46, 103, for details]. The results suggest a z ini = 11.47, which is in the reionization era. The interpretation of the result is that UG is an emergent DE theory, with DE arising during the reionization period.

Einstein-Gauss-Bonet
Einstein-Gauss-Bonet (EGB) is a recent proposition that modifies the geometrical part of the field equations [47], maintaining the continuity equation in its original form. In [48] they constrained the free parameter through diverse observations finding results compatibles with the cosmological standard model. However, the authors also found that specific values of the free parameter could generate an eternal acceleration, even in epochs in which is not expected (reionization, nucleosynthesis, etc). Other mathematical flaws of the EGB model have just been found [150][151][152][153][154].
The action of the EGB gravity can be written in the form [47] where Λ is an effective cosmological constant, R is the Ricci scalar, L m is the matter Lagrangian, α is an appropriate free parameter, G = 6R µν [µν R ρσ ρσ] is the Gauss-Bonnet contribution to the Einstein-Hilbert action and d + 1 is considered in the limit when lim d→3 d + 1 as presented in [47]. Minimizing the action, the field equation can be written as Notice that when α = 0, the standard Einstein field equation with a CC is recovered.
In order to study the background cosmology, we assume a flat FLRW given by Eq.
In order to constrain theᾱ parameter, we divide the problem in two branches through Eq. (105). Therefore, if we only consider the branch where we have a real value of E(z), then we have [48] where is the standard cosmological model. Eq. (106) is constrained to the condition E(0) = 1, having the following relation [48] Notice that whenᾱ → 0, in (106) the standard Friedmann equation is recovered. Fig. 13 show the MCMC analysis implemented using OHD, SNIa, SLS, HIIG and BAO samples. The upper panel shows the reconstruction of H(z) and the bottom panel presents the CL contours for the free parameter of the theorȳ α. In addition, best-fits for the free parameters if the model are presented in Table XIII in conjunction with the χ 2 parameter (see details in [48]).

Cardassian Models
Cardassian 13 models are a phenomenological form to describe the late time acceleration of the Universe that could be theoretically sustained under the assumption of extra dimensions. The models assume an extra function in the Friedmann equations whose form is related to those of the polytropic fluids. Therefore, just the presence of matter and radiation in this functional form automatically produce a late time Universe acceleration without the need for a CC. The main disadvantage is the apparition of additional terms that must be constrained with observations and the difficulties to describe the model from theoretical arguments. It is important to mention that the Cardassian models are divided into the Original Cardassian (OC) and the Modified Polytropic Cardassian (MPC).
The theoretical details of OC model, introduced by Ref. [49], is given by the following Friedmann equation  where n is a free parameter. Despite OC model rise from phenomenological assumption in order to obtain an accelerated Universe, it is important to notice that the mathematical structure can be strictly obtained from brane world models with n-branes immersed in a five dimensional bulk.
In addition, the MPC model [50], takes the form where here l and n are free parameters. Table XIV gives the minimum chi-square and mean values for the OC and MPC parameters using the samples OHD (DA), SNIa (compressed JLA), and the joint analysis of these data sets. Fig. 14 shows the best fit of H(z) for the original (left panel) and modified polytropic Cardassian models (right), respectively. The figures also show the confidence contours (bottom panels) for both Cardassian models using SNIa (cJLA, compressed JLA), OHD uniformed with the sound horizon estimation from Planck data, and the Joint analysis of both data (denoted J3). It is worth to note that both models can reproduce the dynamics of the Hubble measurements.

V. DISCUSSION AND CONCLUSIONS
In this review we presented a brief and non-exhaustive review on DE models that however summarizes some important scientific results, settling the ground for the ideas exposed in this work.
Our motivation to explore alternatives to ΛCDM lies on the problems that afflicts the cosmological constant and the recently observed tension in the σ 8 and H 0 parameters estimated from different samples. Notice how the CC is directly related to the problem of the Universe acceleration and, therefore a competitive model should be implemented. However, regarding the H 0 tension, it is not clear whether it is related to dark energy or not [see for example 156], although some authors assume that this is the case [33].
For each dark energy presented in this review, we have summarized its main characteristics, its theoretical structure and in its capabilities to fit the data provided by modern observations. Parameterizations are always the standard form to tackle the problem of Universe acceleration but there is not a unique way of choosing the form of the function. Furthermore, in many cases there are not solid arguments to justify the chosen form and it is anfractuous to associate the parameterization with some quantum field or with a model that modifies the GR. Models like Chaplygin and Viscous encompass the dark energy and dark matter contribution in just one theoretical framework through a diffusion function in the continuity equations. The theoretical background is robust, with its equations deduced from a quantum field theory in which a scalar field is involved. In addition, the diverse data samples tend to prefer the mentioned models over others with extra complexities.
We also described models with a late apparition of the dark energy (CC always exist in the Universe evolution). The PEDE, GEDE and UG models allow us to estimate the birth of DE at reionization epoch. These kind of models are valuable because their mathematical expression is able to resolve the tension between Planck and SNIa data. Observational constraints also favour them over other models. Furthermore, extra dimensional models contain a theoretical complexity that undermines its recognition, for instance, the RS model has severe faults when it is contrasted with observations. However, the addition of extra degrees of freedom to obtain a variable brane tension reduce the disagreement with observations and shed light on the nature of CC, i.e. it comes from the existence of extra dimensions. The disadvantage is that the CC problem is now dragged to the extra dimensions scenario, introducing difficulties to calculate the expected value in our four-dimensional Universe. EGB model surge as a curiosity in recent literature but, after a more rigorous mathematical exam, there are several notorious flaws in its theoretical arguments. In addition, from a cosmological point of view, it is based on an spurious early acceleration that could cause severe problems with the well-known characteristics of the Universe in epochs such as the nucleosynthesis. Constraints on EGB obtained with the SLS data sample and applied to dynamical systems point out that the early acceleration never ends. Finally, Cardassian models are phenomenological models that could be justified by the assumption of extra dimensions. Although this complicate the equations, they have the advantage of avoiding tensions when contrasted with observations. Notice that, as expected, Cardassian models tend to reproduce the CC, not showing a dynamical DE.
As a final remark, we would like to reflect on the H 0 tension considering the different models discussed in this review. In Table XV we present the compilation of the dimensionless Hubble parameter (h) for all the models with their respective selected samples and priors. We have used a Gaussian prior in most of the cases using Riess or Planck data [31,55], while only in four models we consider a flat prior. For the purpose of the following discussion, we will compare only those models with flat priors which are CPL, UG and both Cardassian models. Only those models using joint samples that include CMB are able to constraint H 0 and, thus the best-fit value is in agreement with Planck's result. However this is not the final word about the H 0 tension, [156] has recently suggested this might be related to a misunderstanding of the distance ladder measurements (i.e. a need for a better agreement between the SNIa absolute magnitude and the Cepheid-based distance ladder) instead of an 'exotic late-time physics'.
Finally, we summarize our results in the following way. The Brane model with constant tension deduced from the RS paradigm and the EGB have several failures which may call into question the viability of both models. Variable brane tension, Cardassian and viscous models can be considered as promising having problems with the complexity of their theoretical background. Parameterizations, Chaplygin, UG, PEDE, and GEDE are excellent competitors against the standard paradigm, with a strong possibility of resolving the H 0 tension and contributing to elucidate the nature of the dark energy component, the Universe acceleration, and its possible consequences in the fate of our Universe.