The Hubble Tension, The $M$ Crisis of Late Time $H(z)$ Deformation Models and the Reconstruction of Quintessence Lagrangians

We present a detailed and pedagogical analysis of recent cosmological data, including CMB, BAO, SnIa and the recent local measurement of $H_0$. We thus obtain constraints on the parameters of standard dark energy parameterizations, including $\Lambda CDM$, and $H(z)$ deformation models such as $wCDM$ (constant equation of state $w$ of dark energy), and the CPL model (corresponding to the evolving dark energy equation-of-state parameter $w(z) = w_0 + w_a \frac{z}{1+z}$). The fitted parameters include the dark matter density $\Omega_{0m}$, the SnIa absolute magnitude $M$, the Hubble constant $H_0$ and the dark energy parameters (e.g., $w$ for $wCDM$). All models considered lead to a best-fit value of $M$ that is inconsistent with the locally determined value obtained by Cepheid calibrators ($M$ tension). We then use the best-fit dark energy parameters to reconstruct the quintessence Lagrangian that would be able to reproduce these best-fit parameterizations. Due to the derived late phantom behavior of the best-fit dark energy equation-of-state parameter $w(z)$, the reconstructed quintessence models have a negative kinetic term and are therefore plagued with instabilities.


Introduction
The success of the current standard cosmological model ΛCDM has been challenged by the mismatch of the value of the Hubble constant obtained by different cosmological probes.
The latest constraint as obtained by the Planck Collaboration [1] is: while the latest local determination as obtained by the SH0ES collaboration [2] is: This 9% mismatch corresponding to a tension of more than 4σ constitutes the well-known Hubble crisis.
This mismatch is equivalent to the mismatch of the Pantheon SnIa absolute magnitudes, which when calibrated using the CMB sound horizon and propagated via BAO measurements to low z (inverse distance ladder, z 1100) have a value [3] M P18 = −19.401 ± 0.027 mag , while when calibrated using the Cepheid stars have a value [4] M R20 = −19.244 ± 0.037 mag . (4) These measurements are in tension at a level of about 3.5 σ, and this is illustrated in Figure 1 [5] assuming Planck/ΛCDM luminosity distance. The data are inconsistent with the M R20 of Equation (4), but they become consistent if there is a transition in the absolute magnitude with amplitude ∆M ≈ −0.2.
Note that in this model there is no transition in the Hubble function H(z).
The SnIa absolute magnitudes M obtained from the distance modulus equation with a luminosity distance obtained from the Planck/ΛCDM Hubble expansion (z > 0.01) are in tension with the SnIa absolute magnitudes obtained from the Cepheid calibrators at z < 0.01 ( Figure 1). However, since the measurements are made at different redshift ranges, the discrepancy can be reconciled by assuming a transition of the absolute magnitude M at z 0.01 by ∆M 0.02, as shown in Figure 1. This transition corresponds to brighter SnIa at high z (early times) compared to low z (late times) since the SnIa absolute Luminosity L is connected with the absolute magnitude as H(z) deformations obtained by dynamical dark energy have been used as possible approaches to the Hubble tension. In this class of models, H(z) becomes deformed from its ΛCDM form (6) in such a way that the CMB anisotropy spectrum remains practically invariant while the Hubble parameter H 0 is shifted to the H R20 0 value (e.g., [6]). This class of models has been shown to require the presence of phantom dark energy at least at late times and has three important problems: • It worsens the growth tension of the ΛCDM model as it indicates larger values of the parameters σ 8 and Ω 0m than indicated by dynamical cosmological probes [6,7] • It provides a worse fit than ΛCDM to low z geometric probes such as SnIa and BAO [6]. • As in the case of ΛCDM , it favors a lower value of the SnIa absolute magnitude M than the local Cepheid calibrators.
One of the goals of the present analysis is to investigate in some more detail the last problem of the H(z) deformation models. In addition to ΛCDM , we consider a generic class of H(z) deformation models with dynamical dark energy and identify their best-fit parameter values, including the SnIa absolute magnitude M and the Hubble parameter H 0 . The cosmological data used include CMB shift parameters, BAO, SnIa Pantheon data and the local determination of H 0 data point (Equation (2)). We thus compare these models with ΛCDM with respect to the quality of fit and the best-fit parameter values with emphasis on the best-fit values of H 0 and M. We thus find the extent to which these models suffer from the M and H 0 tensions and to what extent they imply the presence of phantom dark energy at late times.
Another goal of the present analysis is to identify possible quintessence Lagrangians that are able to reproduce the identified best-fit parametrizations in the context of a physical field theory model. We also test such quintessence models for instabilities and unphysical properties.
The structure of this analysis is the following: In the next section, we provide a pedagogical review of the cosmological data used in our analysis and the statistical techniques utilized. All of our codes are based on Mathematica and are publicly available. In Section 3, we describe the data analysis performed. In Section 4, we present the dark energy models considered, including field theoretic quintessence and dark energy parameterizations. We also present the likelihood contours of the parameters of the models considered and the quality of their fit to the data. In Section 5, we reconstruct a scalar field quintessence Lagrangian that can potentially reproduce the best-fit forms of the parameterizations considered in Section 4. Finally, in Section 6, we conclude, summarize our main results and discuss possible extensions of the present analysis. In the appendices, we provide a pedagogical review of basic cosmological concepts and describe our notation. We also present some derivations of equations used in our analysis.

Cosmological Data-Parameters
Cosmological models are described by parameters whose number for most models ranges between 4 and 20 [8]. The most common parameters used and the ones involved in the present analysis are: • h: The dimensionless Hubble parameter defined as: H 0 = 100 h · kms −1 Mpc −1 .
• Ω m,0 : Present value of the matter density parameter.
• Ω b : The baryon density parameter. • w(z): The equation-of-state parameter. It is also common to use two or more parameters (w 0 , w a , ...) to define it. For example, in the CPL model [9,10] w(z) = w 0 + w a /(1 + z).
From the w(z) parametrization, it is straightforward to obtain the dark energy density parameter Ω DE (z) (see Appendix A). • M: In the present analysis, we also consider the SnIa absolute magnitude M. This parameter can be constrained using either a combination of cosmological data (SnIa, BAO and CMB) at z > 0.01 or Cepheid calibrators at z < 0.01. The root of the Hubble crisis lies in the mismatch of the values of M obtained by the above two distinct approaches, as discussed below.
Our goal is to impose constraints on the values of the parameters of these models using observational data and identify the implications of these values for cosmology in general and for quintessence models in particular. The types of cosmological data considered are Type Ia supernovae, the CMB shift parameters and BAO measurements, which are discussed in what follows. We could also use other data such as cosmic chronometers, i.e., measurements of the Hubble parameter at different redshifts, but they have big error bars, and we do not think they have much constraining power. Thus, we use the data that we think are the most important.

Supernovae as Distance Indicators
Recent data coming especially from distant supernovae indicate that the expansion of the Universe is accelerating. These data come from distance modulus measurements of a certain type of supernovae: Type Ia or SnIa. They can be used as distance indicators (standard candles).
A supernova is a very energetic explosion, which releases vast amounts of energy, as electromagnetic radiation, in a relatively short period of time. Supernovae are extremely luminous and due to the burst of radiation they emit, they can light up their whole galaxy for weeks. A scale to compare the energy that one supernova can unleash is approximately the energy that our Sun will radiate in its whole lifetime (10 44 J). A supernova can only emerge by two mechanisms: • The collapse of the core of a massive star. Such a star has a core mass higher than the Chandrasekhar limit, which is 1.4 solar masses (1.4 M ). • The abrupt re-ignition of nuclear fusion in a compact star (white dwarfs, neutron stars and black holes). In order to have re-ignition, additional energy is required to raise the temperature in the stars core. The star can obtain this energy either by a merger or by accretion.
Supernovae are classified according to their light curves and their absorption line of different chemical elements in their spectra [11,12]. If the spectrum of a supernova includes a spectral line of hydrogen, it is classified as Type II. Otherwise, it is classified as Type I. Now, if the spectrum of a Type I supernova contains a single ionized silicon at 615 nm, it is called a Type Ia, while if it contains a line of non-ionized helium at 587.6 nm, it is called a Type Ib. Otherwise, if it lacks both of these lines, it is called a Type Ic. There is also a way to subdivide the type II supernovae on the basis of their light curves, but it will not be of our interest. Lastly, only Type Ia supernovae emerge from the second mechanism, while every other type emerges from the first.
Type Ia supernovae emerge as a result of the second mechanism, and more specifically, they emerge in binary star systems when one of the companion stars has a mass lower than the Chandrasekhar limit and thus, ends up as a white dwarf. Once the other companion star reaches its red giant phase, the white dwarf starts accreting matter from it due to its gravitational field. When it reaches the Chandrasekhar limit, the degeneracy pressure holding it in equilibrium can no longer keep up with the ever-increasing gravitational pressure and the star shrinks and thus raises its temperature. Eventually, the temperature reaches a point where carbon fusion can take place, which leads to a violent explosion that can be detected by a light curve and has the form shown in Figure 2b. The initial fusion begins a runaway thermonuclear process known as carbon detonation. During this process, large quantities of the radioactive isotope nickel-56 ( 56 Ni) are produced. As it is radioactive, it undergoes positron decay to the stable isotope iron-56 ( 56 Fe) [15]. This decay chain can be simplified as follows: 56 28 Ni → 56 27 The radiation produced during this process has the form of short wavelength gamma rays. This radiation does not contribute immediately to the light curve but only after it has increased its wavelength through interactions with the supernova ejecta. However, not all the gamma rays produced can interact with the ejecta in order to lower their energy, so a percentage of them diffuses through the ejecta. These gamma rays do not contribute to the light curve at a phenomenon called gamma-ray leakage. A qualitative description of the diffusion process can be summarized in four phases, each contributing to the light curve, as shown in Figure 3a.
(a) (b) Figure 3. (a) Typical light curve of Type Ia supernovae along with the instantaneous power produced by the radioactive decay chain of 56 Ni divided into four phases. Adapted from: [15] (b) Typical Type Ia supernova light curves for different amounts of nickel-56 mass produced. Adapted from: [15].
Phase I: At early times, the outer layers of the ejecta are hot and densely packed and have high opacity to radiation of all wavelengths. However, this instantaneous luminosity observed is only a small fraction of the energy radiated from the decay happening in the center of the ejecta and thus contributes to the light curve as a small initial brightness.
Phase II: As the ejecta expand and disperse, its opacity to longer wavelength radiation falls until it becomes completely translucent, and radiation in the UV, optical and infrared range can escape. In this phase, the light curve rises until it reaches its peak.
Phase III: As the ejecta become fully translucent, the trapped radiation, which has been produced earlier, can escape. This corresponds to a rise of luminosity above the instantaneous power from the radioactive decay until this excess amount of radiation has all escaped. This phase corresponds to the part of the light curve right after the peak.
Phase IV: The point where the opacity of the ejecta becomes small enough that most of the short wavelength radiation can escape and produce gamma ray leaks, which do not contribute to the light curve. This can be seen in the part of the light curve where the observed luminosity falls again under the instantaneous power from the radioactive decay.
Furthermore, at late times, in the so-called nebular phase, a significant distribution to the light curve is made by the positrons, which in total can carry about 3.5% of the total decay energy [16].

•
The brightness of the light curve at its maximum is proportional to the mass of the synthesized 56 Ni [15][16][17]. This feature can be seen in Figure 3b. • The light curve width depends on the optical opacity of the ejecta, the total ejected mass and the kinetic energy of the explosion [16], which is calculated as the difference between the energy produced by nuclear fusion and the gravitational binding energy of the white dwarf progenitor [15]. The latter contributes only weakly to the light curve width [16].

Light Curve Features in Modified Gravity
Many studies suggest that the absolute luminosity at peak L peak increases as G decreases. In [18], the assumption: is made, which means that a slow decrease in G results in dimmer distant supernovae than predicted for a standard scenario. In [19], the assumption that there is some relation between the peak luminosity and the Chandrasekhar mass is made and the hypothesis: is tested, with γ > 0 of order unity. These studies have been used to find possible resolutions to the Hubble crisis [20] along with the f σ 8 tension [21]. However, there is a recent study by Wright and Li [15], in which a semi-analytical model is used to calculate the light curves of Type Ia supernovae, which suggests that the SNIa peak luminosity decreases with Chandrasekhar mass. More specifically, the relations between M Ni and G and M Ch and G are being tested as to how they affect the peak luminosity. It seems that the dominant effect by varying G in the peak luminosity is through the Chandrasekhar mass M Ch as the light curves produced with different values of nickel-56 mass M Ni closely match the G = G 0 light curves due to the relationship between M Ni and κ. It is shown that M Ch increases as G decreases, which, as seen in Figure 4, corresponds to a wider light curve but also to a lower peak luminosity compared to the G = G 0 light curve.

Observations
The absolute luminosity of Type Ia supernovae is almost constant at the peak of brightness, so the distance to it can be determined by measuring its apparent luminosity. Thus, they are a kind of standard candle (by which luminosity can be measured using observations). However, in reality, things are not so easy. The intrinsic spread in absolute magnitudes is too large to produce any meaningful cosmological constraints.
On the other hand, at the end of the 1990s, a high-quality sample of local supernovae (z << 1) helped towards the correlation of the absolute magnitude M and the width of the light curve. Therefore, if someone measures the apparent magnitude of a Type Ia supernova and the width of its light curve, they can predict its absolute magnitude. Thus, the universal SNIa absolute magnitude is nothing but the corrected magnitude in terms of the light curve width, as seen in Figure 5. Therefore, a more appropriate term when referring to them is not standard but standardizable candles.
However, there is another correction that we need to apply, since the redshift increases when we observe different parts of the power spectrum (broader, brighter SNe appear slightly bluer than they should be [22]). This correction is named K-correction, and we always assume that it has already been included in the estimation of the apparent magnitude that we use.
Type Ia supernovae are the preferred distance indicators as: (a) (b) Figure 5. (a) Light curves of some local Type Ia supernovae as measured. (b) The light curves of the same supernovae corrected in terms of their width. Adapted from: [14].

•
They are the most common type of supernova in the Universe.
• They are extremely luminous, as at their peak luminosity they can reach an absolute magnitude of about: M ≈ −19, which is about the absolute magnitude of a bright galaxy. • They have a relatively small dispersion of peak absolute magnitude. • Their explosion mechanism is fairly uniform and well understood and, according to known physics, has no cosmic evolution. • There are a lot of local SNeIa that we can use to test their physics and calibrate the absolute magnitude for the distant ones.

Baryonic Acoustic Oscillation Measurements (BAO)
The high redshift Universe (z > 1100) was pretty much homogeneous except from some really small perturbations (all four species were perturbed approximately by the same fractional amount) and consisted mainly of four density components: dark matter, baryons, photons and neutrinos. Photons and baryons were tightly coupled due to Compton scattering. Neutrinos do not interact and move too fast, so gravity cannot stop them. On the other hand, dark matter becomes attracted and falls into the perturbations' overdensities, due to gravity. The perturbations of the photon-baryon fluid (they are coupled) have both overdensity and overpressure. This overpressure creates an expanding sound wave moving with the speed of sound c s of that time. Thus, the perturbation of photons/baryons is carried outwards.

Recombination
As the universe cools down, there is a point (approximately at redshift z ≈ 1320) that protons and electrons can combine to form hydrogen: e − + p + ←→ H + γ. At this point, photons do not scatter as efficiently and start to decouple, while the sound speed drops and thus the sound wave begins to slow down.

Photon Decoupling
The same process continues until the photons decouple completely. Then photons' perturbation begins to smooth out and the sound speed of the baryon perturbation is reduced so much that the pressure wave stops. Thus, we are left with the original dark matter perturbation surrounded by the baryon perturbation in a shell, as shown in Figure 6.

Beyond Photon Decoupling
As can be seen in Figure 7e,f, the two perturbations left attract each other and they start to mix, eventually coming together. Although, the acoustic peak perturbation is lower than the dark matter one since dark matter mass is dominated by baryons.

Structure Formation
Galaxies form in matter (both dark matter and baryons) overdensities. Most galaxies appear at the original perturbation. However, a 1% enhancement of galaxies appears at the acoustic scale and can be seen in the galaxy correlation function.
The time that baryons are released from the drag of the photons is known as the drag epoch. We can obtain the redshift z d of this epoch, which is shortly after recombination, using the fitting formula introduced in [24]: where: Photon-baryon fluid acoustic waves propagate at the sound speed: where: It is also easy to prove the equation: that we will use in the data analysis later. Thus: The sound horizon depends on: • The epoch of recombination, which affects the drag epoch z d . • The expansion of the Universe, H(z).

•
The baryon-to-photon ratio, which affects c s . Adapted from: [25], while the original animation can be found here: [26].

BAO Measurements
The BAO scale can be found as a peak in the galaxy correlation function or equivalently as damped oscillations in the large-scale structure power spectrum, as seen in Figure 8. In spectroscopic surveys, we observe the angular and redshift distributions of galaxies as a power spectrum P(k ⊥ , k ) in the redshift space, where k ⊥ and k are the wavenumbers perpendicular and parallel to the direction of light, respectively. Two possible measurements that we can perform are in the line-of-sight dimension or in the transverse direction measuring the ratios: respectively, with d A being the angular diameter distance and r s (z ) the sound horizon at the decoupling epoch. Through these observables, we can define the distances: The first helps us measure the expansion rate H(z), while the second is just the comoving angular diameter distance. Here, we can see from Equation (A79) that the comoving angular diameter distance D M corresponds to the metric distance d m , which we defined in Equation ( A72) We can also use the spherically averaged spectrum to obtain a combined distance scale ratio: 1/3 (22) which relates with the distance: However, as we saw earlier, r s can change due to deviations of the cosmological parameters. This means that BAO measurements do not really constrain the values mentioned above, but they constrain the values: Figure 8. (a) The large-scale redshift-space correlation function of the SDSS LRG sample. Adapted from: [27]. (b) Damped oscillations in the large-scale structure power spectrum. Adapted from: [28].
where r f id s is the sound horizon in the context of the fiducial cosmology assumed in the construction of the large-scale structure correlation function.

CMB Measurement
The CMB, shown in Figure 9a, can be treated as a BAO measurement at z = z = 1090, as shown in Figure 9b, measuring the angular scale of the sound horizon at a high redshift [29]. Thus, similar to the BAO measurements, we can define the characteristic angle of the location of the peaks: where r s (z) is the comoving sound horizon at redshift z and z is the redshift to the photondecoupling surface. The multipole l (for the full analysis, see [12]) corresponding to the angle θ A is: Now, the comoving angular diameter distance at that epoch D M (z ) can be expressed as: where R is the CMB shift parameter.
Shift Parameter R From References [6,31,32], we know that the CMB power spectrum, which can be seen in Figure 10, will be almost identical if these parameters are fixed:

•
The physical density parameters for matter ω m , baryons ω b , radiation ω r and curvature ω k . • The primordial fluctuation spectrum. • The flat-universe comoving angular diameter distance to the recombination surface D M (z ).
Furthermore, the product: √ ω m · D M (z ) does not depend on H 0 since: since: H 0 ≡ 100 · h · kms −1 Mpc −1 . This combination is the well-known shift parameter: In general, the shift parameter is defined as [31,32]: where: S(x) is defined in Equation (A3) and y is: In a flat universe: Thus: However, it is also common in the bibliography [12,33] to define the shift parameter as: in natural units (NU), or as: in SI units.
In this analysis, the definition used for the shift parameter is the latter one.

Photon Decoupling Epoch
In order to study the CMB, we should be able to find the value of the redshift to the photon-decoupling surface z . In spite of the complexity of such a procedure, there is a fitting formula derived by Hu and Sugiyama [33,35]: where: This fitting formula is designed to be valid at a percent level, assuming that the parameters Ω b , Ω m are inside the ranges: 0.0025 Ω b 0.25 (41) 0.025 Ω m 0.64 (42)

Bayes Theorem
Bayes Theorem states that: where: • p(θ|x): The posterior probability distribution p(θ|x) for the parameters θ and data x.
It is the probability that the parameters will obtain certain values after completing the experiment and making some assumptions [36]. • p(x|θ): It is called likelihood, and we also refer to it as L(x; θ) • p(θ): The prior probability distribution p(θ) for the parameters θ. It expresses what we know about the parameters before performing the experiment, including the results of previous experiments or theory. For example, we know that the age of the Universe must be positive.
In the absence of any previous information, it is common to adopt the principle of indifference and assume that all values of the parameters are equally likely and take p(θ) = constant. As a bound, someone can either use some finite bounds or use infinite bounds and work with an unnormalized prior. This prior is called a flat prior. • p(x): The evidence.

Maximum Likelihood
Consider data consisting of M independent measurements X i,obs with known standard deviation σ i . Consider also a theoretical model prediction X th (p) to be tested by the data, with p representing the parameters used in the model. The χ 2 function of the model as a function of its parameters is defined as [11]: In reality, χ 2 quantifies the discrepancy between the predictions of the theoretical model and the observations at a particular value of the parameter p. Thus, a small value of χ 2 indicates a good fit. The values p 0 that minimize the χ 2 (p) function are called best-fit parameters.

Non-Independent Measurements Covariance
If X and Y are two random variables, their covariance, denoted as cov(X, Y), is defined by [37]: where X and Y are the mean values of X and Y, respectively. In general, for discrete random variables: with h(x i , y j ) referring to the joint distribution of X and Y.
For continuous random variables: where, again, f (x, y) is the joint distribution of X and Y. It is evident that:

Correlation
The correlation of X and Y is denoted as ρ(X, Y) and is defined by: while for discrete variables, we have: Covariance Matrix Having that in mind, we can construct the covariance matrix for the measurement vector X i as: As we can see, the covariance matrix is always symmetric and square.
When we do not have independent data, we cannot define the χ 2 function as before since we should incorporate the covariance matrix in its definition. Thus, χ 2 (p) takes the form:

Chi-by-Eye
There is a useful "chi-by-eye" rule, which states that a fit is good and believable when the minimum χ 2 is roughly equal to the number of data minus the number of parameters, and it is increasingly true for a large number of data [38].

Likelihood Function
We define the corresponding likelihood function L as: The likelihood function maximizes as χ 2 (p) minimizes. The parameters that result in a higher value of the likelihood function are more likely to be the true parameters.
If the model has two parameters, we can visualize the likelihood function as a surface sitting above the parameter space while, in general, for a model with N parameters, the likelihood function takes the shape of a hypersurface spanned by the parameter vector p [39].

Fisher Matrix
Sometimes we want to estimate the errors of the parameters from the likelihood function. At first, we assume a flat prior so we can identify the posterior through the likelihood. Then, we can expand the log-likelihood lnL as a Taylor series around its maximum [38]: Thus, we expand around the maximum log-likelihood ln L(p 0 ): for every p j . Furthermore, we stop the expansion to the quadratic term, which implies that we say that the likelihood surface is a multivariate Gaussian. Thus:

Hessian Matrix
The Hessian Matrix is defined as: It encloses information about the errors of the parameters and their covariance. If the Hessian is not diagonal, the parameter estimates are correlated, which means that they have a similar effect on the data, and it can be difficult to discern them using the data. However, it is not certain that the parameters themselves will be correlated.

Conditional Error
If all parameters are fixed except, e.g., p i , we can estimate its error as: This is called a conditional error. It is the minimum error bar attainable on p i , when the other parameters are known [36], but it is uninteresting and almost never used.

Fisher Matrix
In 1935, Fisher answered the question of how accurately someone can measure model parameters from a given dataset. The Fisher information matrix is defined as: or as: In practice, we choose a fiducial model and compute the Fisher matrix using [40]. Substituting: in the Equation (56), for a one parameter case, we obtain: Here, by identifying that 2∆ ln L is really ∆χ 2 , we see that the 1σ displacement for the parameter p i , when 2∆ ln L = 1, is: which is analogous to the conditional error. In general, where F −1 refers to the inverse of the Fisher information matrix, and σ p i is the expected marginal error [36,38]. The Fisher matrix approach always gives you an optimistic estimate of the errors and this is obvious from the inequalities above, which are forms of the Cramér-Rao inequality [36,38]. The inequality becomes an equality only if the likelihood is Gaussian.

Marginalization
Most of the time, even though we use many parameters in a model, we may not have interest in all of them or we may have included some nuisance parameters, i.e., parameters whose values we know with limited accuracy.
Thus, we want to report the confidence level of the cosmological parameters that are interesting to us regardless of the value of the uninteresting/nuisance ones. This happens through a process called marginalization, where we marginalize over them [38]. For example, in cosmological models with dark energy, a nuisance parameter is the Hubble constant H 0 , or in cosmological models with curvature, we want to know the value of Ω k regardless of the values of Ω m or Ω Λ . To perform the marginalization process, we need, at first, to estimate a prior distribution π(ν). A reasonable choice is to use a Gaussian probability density function with a mean value of ν 0 (the most likely value) and variance σ ν .
In this way, we can build a posterior likelihood function that will depend only on the interesting parameters p and not on the parameter ν: While if we consider the reasonable choice of π(ν) (Gaussian), we obtain: Thus, in the same way as before, by minimizing χ 2 (p), we can find the best-fit parameters [11]. Figure 11. Probability distribution for the case M = 2 (two parameters a 0 and a 1 ). Furthermore, there are three different confidence regions all at the same confidence level. The first region is defined by the vertical lines and represents a 68% confidence interval for the variable a 0 without regard to the value of a 1 , while the second one is defined by the horizontal lines and represents a 68% confidence interval for the variable a 1 without regard to the value of a 0 . The third one is the ellipse, which shows a 68% confidence interval for a 0 and a 1 , jointly. Adapted from: [41].

Confidence Limits
It is common practice not to present every detail of the probability distribution of errors in parameter estimation but to summarize it in the form of confidence limits. The full probability distribution is a function defined in the M-dimensional space of parameters p (M = number of parameters).
A confidence region is a region in M-dimensional space that contains a certain percentage of the total probability distribution. The ideal is to find a small region that contains a large percentage of the total probability distribution.
When we perform an analysis, we are free to pick both the confidence level and the shape of the confidence region. The only requirement is that the region we pick contains the stated percentage of probability. The most commonly used percentages are: 68.27%, 95.45%, 99.73%, which correspond to standard deviations 1σ, 2σ, 3σ from the most likely value, while the convention when we want to choose a shape for a confidence region is: for one dimension, we use a line segment centered on the measured value, and for two or higher dimensions, it is most common to use ellipses or ellipsoids.
The whole point of the confidence level is to inspire confidence; for example, when we say that we have a confidence region with a confidence limit of 99%, there is a 99% chance that the true parameter falls within this region around the measured value. An example is shown in Figure 11. Confidence regions derived using a constant ∆χ 2 region as a boundary. Here, the shape of the confidence region is chosen to be an ellipse. It is important to note that the intervals AA', BB', CC' are the ones that contain the percentage of the normally distributed data that correspond to the respective ∆χ 2 . Adapted from: [41].

Constant χ 2 Boundaries as Confidence Limits
In order to obtain the minimum value of χ 2 for the observed data set, we use the value p = p 0 for the parameters. If the vector of the parameters' values p is perturbed away from p 0 , χ 2 increases. The region in which χ 2 does not increase more than a fixed amount ∆χ 2 defines some M-dimensional confidence region around p 0 . If ∆χ 2 is a large number, the confidence region will be large, while smaller values for ∆χ 2 correspond to smaller regions.
Although we are free to use whatever value we want for ∆χ 2 , there are some special values that correspond to the most commonly used confidence limits: 68.27%, 95.45%, 99.73%. These values do not only depend on the preferred confidence limits but also on the number of parameters M (or equivalently, the degrees of freedom of the model) as shown in Table 1. The scale of the likelihood contours as the value of ∆χ 2 changes is shown in Figure 12.

Errors
When we use a model, it is really important to know not only the favored values obtained by the data but also their respective errors. In order to simplify the procedure of finding the errors and have a good approximation of them, we can use a linear model of the form: with errors in both axes: σ x i , σ y i At first, we need to calculate the χ 2 function for this model. To do, that we need to know the weighted sum of variances, which in this case is the variance of the linear combination: y i − a − bx i of the random variables x i , y i [41].
Thus: Figure 13. Standard errors for the parameters a and b found as the projections of the ellipse, which is the confidence region boundary with ∆χ 2 = 1, on the a and b axes. Adapted from: [41].
The goal is to minimize χ 2 (a, b) with respect to a and b, so we need to solve: in order to find the best-fit values for the parameters a and b. As we saw earlier, a confidence region boundary where χ 2 takes a value greater than its minimum: χ 2 + ∆χ 2 with ∆χ 2 = 1, defines a 1σ region. Thus, by taking the projections onto the a and b axes, we obtain the standard errors for the parameters a and b, respectively, as shown in Figure 13. In the linear case, we can obtain these projections by Taylor expanding χ 2 as follows: which due to Equations (71) and (72), which are true at the minimum of χ 2 , and the fact that: Our goal is to solve the equation ∆χ 2 = 1 numerically for all the parameters of our model. As a result of this, in combination with Equation (62), we find the components of the Fisher matrix, and from the diagonal components of its inverse, we obtain the marginal errors for the parameters.

χ 2 for CMB Data
In this analysis, we use the central values for R, vl A and Ω b h 2 from Planck 2018 for a flat universe. The data can be found in [33] and are expressed in terms of a data vector and a covariance matrix. The data vector is: and the corresponding covariance matrix is: The only remaining task is to find the inverse matrix of C v and define the χ 2 CMB . It is useful to define a new vector: Thus, the χ 2 function takes the form:

6dFGS and WiggleZ
Both surveys give a measurement for the ratio: From 6dFGS we obtain: d z (z = 0.106) = 0.336 ± 0.015, while from WiggleZ, we obtain three correlated measurements: with the inverse of the covariance matrix given as: Having four measurements, we produce the data vectors: and the total inverse covariance matrix: The vector used to define the χ 2 function is: and the χ 2 function for these measurements is:

SDSS
From this survey, we obtain a measurement of: where: r d = r s (z d ) (86) and its value in the context of the fiducial cosmology is: Thus, we obtain the data vectors: The χ 2 function: χ 2 SDSS is obtained by using Equation (44) for the parameter D V r d , which is essentially the ratio: 1 d z . Thus: In this analysis, d A is used instead of D M so the new data vector for the first measurement is: The corresponding inverse covariance matrix for u 1 and u 2 is: The vector used to define the χ 2 function is: and the χ 2 function for these measurements is: Finally, the total χ 2 function for the BAO data is the sum of the individual χ 2 functions for the data from the four surveys. Thus:

χ 2 for SNIa Data
For the SNIa data, we use the Pantheon dataset. A representation of a subset of the Pantheon data called the binned Pantheon dataset is shown in Figure 14. The data vector used for the definition of χ 2 is: where the apparent magnitude m th is obtained using Equation (A87) or equivalently as: where: is the Hubble-free luminosity distance [42]. However, the parameters M and H 0 are degenerate. Therefore, it is common either to marginalize them as nuisance parameters [5,43,44] to make the analysis of a cosmological model using only SNIa data. Then, the m th is defined in terms of M as: Thus, the appropriate χ 2 function is: where the total covariance matrix C total ij is obtained as: with D ij being the diagonal matrix: and C systematic is a non-diagonal matrix associated with the systematic uncertainties that emerge from the bias corrections method [5,42] and can also be found in the same sources with the dataset.

Total χ 2 function
As in the BAO analysis Equation (95), the total χ 2 function is obtained as the sum of all the individual χ 2 functions used in the analysis. Thus: which is the function that we want to minimize in order to find the best-fit parameters for a given cosmological model that we want to test. Useful information about data analysis can be found in Appendix C.

Dark Energy Models
The discovery of cosmic acceleration [45,46] has lead to a wide range of approaches for its description and its modeling. The main approaches include either the introduction of a fluid in the context of general relativity with a parametrized equation-of-state w(z) in the energy momentum tensor or the introduction of new scalar fields in the context of GR or modified gravity. In the present analysis, we assume a flat universe. This assumption is valid since the Planck collaboration [1] derived a value for the curvature density using the joint constraint from the CMB and BAO measurements: Ω k = 0.001 ± 0.002, which is consistent with a flat universe.

Spatially-Flat ΛCDM Model
The simplest form of dark energy is the cosmological constant Λ [47]. It corresponds to a time-independent energy density: obtained from an ideal fluid with equation of state w = −1 and non-relativistic matter, i.e., cold dark matter (CDM). This corresponds to the ΛCDM model. Assuming a flat universe in the presence of matter and a cosmological constant, Equation (A35) gives: or in terms of the dimensionless density parameters: with: This is the minimal standard model of cosmology, since the predicted Hubble parameter H(z) has only two free parameters that can be constrained through observations. Even though ΛCDM is simple and is not yet excluded by the observations, it faces some challenges both theoretically and observationally.

Theoretical Challenges
The two most notable theoretical challenges are the fine tuning problem [48,49] and the coincidence problem [50,51].

Fine-Tuning Problem
Since we now observe the cosmic acceleration, we require that the cosmological constant Λ is of the same order of magnitude as the square of the present Hubble parameter H 0 .
Using the unit conversion for km to Mpc and frequency s −1 to the corresponding energy in eV, we obtain: This value corresponds to an energy density: which by substituting h ≈ 0.7 and G ∝ m −2 pl gives: We can assume that this energy density comes from the vacuum energy of empty space, which is mathematically equivalent to the cosmological constant [8,47]. We can obtain this vacuum energy by summing over all the zero-point energies of some field with mass m, momentum k and frequency ω up to a cut-off scale k max >> m [12].
Since the integral is dominated by the large k modes (k >> m): A reasonable cut-off scale is the Planck scale (m pl ) up to which general relativity is believed to hold. Thus: We can see that: which means that the predicted value is 121 orders of magnitude larger than the observed one.
Another cut-off scale that can be used is the one when supersymmetry (if it exists) breaks at around 1 TeV [11]. Then: which is a value 56 orders of magnitude larger than the observed one. This discrepancy between the much larger initial values for the vacuum energy is also called the "smallness" problem [48] and is also referred to as the cosmological constant problem/puzzle [11,47].

Coincidence Problem
The lack of explanation for why dark energy has the same order of magnitude as nonrelativistic matter density at the present epoch, which is obvious in Figure A1, is called the coincidence problem [11]. This is quite bizarre since ρ m scales as a −3 , while ρ Λ scales as a 0 and: which indicates that Λ was negligible in the past and will dominate the future. If the cosmological constant is considered to be an initial condition in the early Universe, it seems really unlikely that Λ should have a value comparable to matter at the present cosmological epoch while galaxies and other large-scale structures have formed [11]. A really interesting fact is that if the cosmological constant had a couple orders of magnitude higher energy density, large-scale structure would not have formed, while if it was some orders of magnitude lower, it would not have been detectable. A possible solution to this problem is the anthropic principle.

Anthropic Principle
The general idea of this principle is that physical theories should take into consideration the existence of life on Earth [12].
Carter was the one that used the expression "Anthropic Principle" and proposed two variants for it [52]. The first is the weak anthropic principle, which states the spacetime position of life in the universe is privileged to the extent of being compatible with our existence as observers, while the second one, the strong anthropic principle, states that the universe and the fundamental physical constants must be such as to admit the creation of observers within it at some point.

Observational Challenges
In general, these observational challenges occur when different observations favor different values for the same parameter, and the most notable ones are the H 0 and growth tension. A more detailed analysis for these challenges along with other tensions and curiosities concerning the ΛCDM model can be found in: [53].

H 0 Tension
As discussed in the Introduction, using a ΛCDM background cosmology and data from the CMB and BAO measurements, a value of H 0 = 67.4 ± 0.5 kms −1 Mpc −1 is obtained [1] for the Hubble constant. On the other hand, using local measurements coming from SnIa data leads to a value of: H 0 = 74.03 ± 1.42 kms −1 Mpc −1 [54]. In general, this discrepancy can range from 4.4σ to over 5σ depending on which local data are used [4,55]. More detailed presentations of the subject along with reviews of the solutions that have been proposed can be found in [56,57].

Growth Tension
The amplitude of the primordial power spectrum, measured through the parameter σ 8 , which is the linear amplitude of matter fluctuations on scales of 8 h −1 Mpc, and the matter density parameter Ω 0m are the two parameters that affect the growth rate and magnitude of linear cosmological perturbations [42]. Dynamical probes for the cosmological perturbations' growth rate, which measure it directly, such as weak lensing [58][59][60] and redshift space distortion [21,61,62], indicate that the observed growth rate is weaker than the theoretical prediction obtained in the context of a ΛCDM background using the observed background expansion rate measured through Type Ia supernovae, BAO and CMB data, which are geometric probes. Dynamical probes favor smaller values for both Ω 0m and σ 8 than the geometric ones, and this discrepancy varies from around 2-3 σ [63-66] to even more than 5σ [67] depending on the model parametrization and the dataset used. However, if the constraints from the CMB on the ΛCDM background are not taken into account, leaving only the constraints from the SNIa data, the tension decreases to a level below 2σ [68]. The tension also decreases when marginalized confidence contours are used [69] and can be resolved completely within the Minimal theory of massive gravity scenario [70]. A brief analysis of the growth tension along with some proposed solutions can be found in [71].

Fitting the ΛCDM parameters: Maximum Likelihood BAO and CMB Data
In the Mathematica code used, we construct the χ 2 function for the model: χ 2 (Ω 0m , h) where we use the value: Ω b h 2 = 0.02236 from [1], with which we perform two separate analyses: one solely with BAO data and one with both the BAO and the CMB data. At first, we use only the BAO data: with χ 2 BAO defined in Equation (95). Using both the BAO and CMB data, we obtain: with χ 2 CMB defined in Equation (78).

SNIa Data
In this analysis, we use the Pantheon dataset and try to constrain the parameter Ω 0,m and the degenerate combination M. Thus, we construct the χ 2 function as: Pantheon defined in Equation (101).

Combined Data
At last, we combine all the data, and the new χ 2 function is defined as: since BAO and CMB measurements explicitly constrain the Hubble constant. Implementing the maximum likelihood method [38,41,72], i.e., minimizing the χ 2 function over all free parameters, we obtain the best-fit values for each data combination, which can be seen in Table 2.

Results
The constraints obtained using the Pantheon dataset are in agreement with corresponding previous studies [5,42,73].
From Figure 15, we can see that the BAO measurements cannot significantly constrain the parameters h and Ω m , while the addition of the CMB measurements dramatically improves the constraints. As expected, the preferred value of h is consistent with (1) and inconsistent with (2), which demonstrates the H 0 tension. Since the best-fit value of M remains practically unaffected by the addition of BAO and CMB data, we conclude that the M-tension remains present. Figure 16 shows the constraints obtained from the Pantheon dataset. When the measurements from the geometric probes (BAO and CMB) are added, the constraints on Ω m are dramatically improved, while the constraints on M are practically unaffected.  Figure 16. The 1σ − 3σ confidence contours in the parametric space (M − Ω 0m ).The red contour corresponds to the full Pantheon dataset while the blue contour corresponds to the Pantheon+CMB+BAO data combination. As expected, when the CMB+BAO data are included the constraints on Ω 0m improve dramatically.

Spatially-Flat wCDM Model
This generic parametrization of the dark energy equation-of-state parameter is of the form which is assumed to be an arbitrary constant. When w 0 = −1, the WCDM model reduces to the ΛCDM model. Assuming a flat universe in presence of matter and a spatially-homogeneous fluid with w < −1/3, Equation (A35) gives: or in terms of the dimensionless density parameters: with: In order to impose constraints on the model parameters, we use the combined data from the Pantheon dataset and BAO and CMB measurements. However, instead of M, we consider M and h separately, and thus, the m th used in the Pantheon dataset analysis is the one defined in Equation (97). Thus, the constructed χ 2 function is defined as: Implementing the maximum likelihood method, we obtain the best-fit values, which can be seen in Table 3. In order to produce the contour plots for the wCDM model, we vary all four parameters in the χ 2 function and then we show a two-dimensional projection of the ellipsoid produced in the four-dimensional parameter space.

Results
When we compare the values favored by the data for the wCDM model with the respective values for the ΛCDM model, we can see that it favors smaller values for the equation-ofstate w 0 and the matter density Ω 0,m along with higher values for the dimensionless Hubble parameter h, while giving, at the same time, a slightly better fit that can be seen from the lower value of χ 2 . In addition, the favored range of the SnIa absolute magnitude M remains inconsistent with the value range indicated by local Cepheid calibrators (Equation (4)), i.e., the M-tension remains along with the Hubble tension, as shown in Figure 17.

Chevallier-Polarski-Linder (CPL) Parametrization
A commonly used ansatz, proposed by Chevalier-Polarski [9] and Linder [10], allows for dynamical dark energy and is based on an expansion of w(a) around the present value of the scale factor a = 1: This can be derived as the Taylor series expansion of w(a) around a = 1 including only linear terms: where w denotes the derivative in terms of the scale factor a.
In terms of the redshift z, the equation-of-state w(z) takes the form: Its present value w(a 0 ) is: and the value of its slope: Again, implementing the maximum likelihood method, we obtain the best-fit values, which can be seen in Table 3. In order to produce the contour plots for the CPL model, we vary all five parameters in the χ 2 function, and then we show a two-dimensional projection of the ellipsoid produced in the five-dimensional parameter space.     Figure 19. The 1σ − 3σ confidence contours in the parametric spaces: (w a − Ω 0m ), (w a − M), (w a − h) and (w a − w 0 ) for the CPL model.

Results
When we compare the values favored by the data for the CPL model with the respective values for the ΛCDM and wCDM models, we can see that CPL favors even smaller values for the equation of state w 0 at the present time. Furthermore, it gives a value for the matter density Ω 0,m around that predicted from the wCDM analysis with for the dimensionless Hubble parameter h, indicating also the presence of the Hubble tension. The M tension also remains since the preferred value of M is significantly lower than the corresponding Cepheid value, as seen in Figure 18. Furthermore, the positive values preferred for w a , as can be seen in Figure 19, indicate a crossing of the phantom divide line at w = −1. Finally, despite the additional parameters involved, the quality of fit is similar to that of ΛCDM and wCDM since χ 2 is similar in all three parametrizations, as shown in Table 3.

Adding the Local H 0 Determination
Taking into account the locally determined the value of the Hubble constant (2)) (hereafter abbreviated as "the Riess data point") corresponds to adding a term in χ 2 as follows to observe how the best-fit parameters are affected. Thus, the new χ 2 function is: The new contours produced with the use of the new χ 2 function for every model can be seen in Figures 20 and 21, while we obtain the best-fit values, which can be seen in Table 4

Results ΛCDM
In Figure 20 and in Table 4, we can see that adding the Riess point leads to slightly smaller values for both M and matter density Ω 0,m . The quality of fit also worsens significantly since χ 2 increases by about 22 units compared to the case when the Riess point is not included. wCDM In Figure 21a and in Table 4 we can see that the addition of the Riess point indicates that the best-fit absolute remains in tension with the value indicated by local Cepheid calibrators M R20 . In addition, the new data combination prefers a somewhat higher value for the Hubble constant and a lower value for the matter density. However, the Hubble and M tensions remain despite the introduction of the Riess point. Lastly, wCDM in this case favors a lower value for the equation of state w, thus moving further away from the cosmological constant value w = −1. The quality of fit also worsens significantly since χ 2 increases by about 15 units compared to the case when the Riess point is not included. This worsening is not as bad as in the case of ΛCDM but indicates that the Riess point is clearly not consistent with the rest of the cosmological data even in the context of wCDM. Figure 22a, it is clear that the addition of the Riess point has similar effects in the CPL as in wCDM model; somewhat higher values for the Hubble constant, lower values for the matter density, slightly higher SNIa absolute magnitude, and lower value for the equation of state at the present epoch w 0 . Lastly, the parameter w a remains almost unaffected as the change in its best-fit value is really small. Clearly, however, the Hubble and M tensions remain for this H(z) deformation model.

CPL From
The quality of fit is very similar as in the case of wCDM despite the extra parameter. It is also significantly worse compared to the case when the Riess point is not included. This worsening is not as bad as in the case of ΛCDM but indicates that the Riess point is clearly not consistent with the rest of the cosmological data also in the context of CPL.  The decrease of the preferred value of Ω 0m in all cases when the Riess point is included is anticipated due to the property of the CMB anisotropy spectrum to favor a fixed value of Ω 0m h 2 . Thus, we anticipate that the introduction of the Riess point, which tends to raise h 2 , would also tend to lower the favored value of Ω 0m in order to maintain the product Ω 0m h 2 as approximately constant and consistent with the observed CMB anisotropy spectrum.

Scalar Field Dark Energy Models (Quintessence)
In the context of quintessence [74], we consider a self-interacting canonical scalar field φ minimally coupled to gravity playing the role of dark energy [75]. This scalar field is described by the Lagrangian density: where V(φ) is the potential energy density of the field φ.
The stress-energy tensor can be obtained as [12]: which reduces to: At this point, it is important to calculate the terms: where we used: dg = gg ab dg ab The detailed proof can be found at Appendix D. Furthermore: Thus, the stress-energy tensor takes the form: which reduces to: or by substituting L, we obtain: Assuming that the scalar field is close to spatially uniform on cosmological scales, we can neglect its spatial derivatives ∂ i φ compared to its time derivativesφ. Thus: which means that T µν is diagonal.
We can obtain the energy density and the pressure of the field as: Thus, the equation-of-state parameter is: In general, forφ >> V(φ), i.e., for a kinetic energy dominated field: while forφ << V(φ), i.e., for a potential energy dominated field: while we can effectively obtain the cosmological constant w = −1 forφ = 0. Quintessence can play the role of dark energy if: However, this is not enough since dark energy domination today requires w = −1 not only now but for an extended period of time (roughly between z ≈ 1 and now). Thus, it is a requirement that the conditionφ < V(φ) holds for a while. This can happen if the time derivative of this condition is also fulfilled: where V (φ) = dV(φ)/dφ. In summary, a scalar field can play the role of dark energy if: These are the slow-roll conditions [76]. The time evolution of the scalar field is determined by the Klein-Gordon equation:φ which can be obtained by the variation of the action: The detailed proof can be found in Appendix E. Assuming a flat universe in the presence of matter and the quintessence Equation (A35) gives: or in terms of the dimensionless density parameters: where we used Equation (A25) for the dynamical energy density ρ φ . Furthermore, This model has Ω 0m as a parameter, along with the number of parameters that are used to parametrize the equation of state w(z).
There are two ways to approach the effect that different dark energy models have in the cosmological expansion: either calculate the equation of state for some specified theory and then its effects on the cosmological expansion or start from the observations of the cosmological expansion and then reconstruct the scalar field physics responsible for the effects observed.
The latter approach is made difficult due to some issues [75]: • Noisiness of measurements of the expansion. • Translation from the measured quantity to ρ DE and w through one or two derivatives. • Range of the scale factor or equivalently redshift coverage: z = 1 a − 1 However, since we already have the constraints for the model parameters, this is the right approach, and others have used it before [77][78][79][80][81][82][83][84][85].

Reconstruction Equations
At first we have the parametrization of the equation of state w: From this equation, we see that in reality we can assume that wCDM and ΛCDM are special cases of the CPL model with w a = 0 and w 0 = const and w 0 = −1, respectively. Thus, we can reconstruct their fields, too.
Furthermore, we can obtain the energy density using Equation (A25): Then, the Hubble parameter H(z) is: From Equations (146) and (147), we can obtain the potential in terms of the redshift: and the kinetic term: In this way, we obtain the field φ in terms of the redshift as: Furthermore, we can reconstruct the potential in terms of the field V(φ) by constructing points with the values of V(z) and φ(z) calculated for many different values of the redshift and then plotting them.
Thus, the only remaining objective is to substitute the best-fit values for the model parameters obtained through maximum likelihood estimation.

Results
We can see from Figure 23a that the best-fit value for the wCDM model gives a phantom field, well explained here [86], i.e., a field, with w < −1 which should have a negative kinetic term.
Furthermore, we can see that the data favor a value for the equation of state at the present epoch w 0 that makes the CPL model a phantom field and a value for w a that makes it raise that value with redshift. Thus, at some point (at around z ≈ 0.4) it crosses the phantom divide line (w = −1), i.e., the line that separates the physics that obeys the null energy condition (ρ + p ≥ 0) from the physics that violates it [75]. This feature is responsible for the weird behavior in the plots of the fields φ and V(φ) in Figure 24.
It is expected that both models will break down trying to explain this phantom regime since, as seen in Equation (150), w is allowed to take values higher than −1. In general, though, ghost solutions are problematic because they lead to instabilities [87,88].

Discussion and Conclusions
We have presented constraints of the parameters of some generic dark energy parameterizations, including ΛCDM , wCDM and CPL. In the model parameters, we included the SnIa absolute magnitude M. In order to constrain these parameters, we used the SnIa Pantheon dataset, CMB shift parameters and BAO data. We also included the local determination of H 0 as a data point and demonstrated that the Hubble tension persists for these parameterizations even after the local H 0 point is included in the data. This tension manifests itself in three ways: • The best-fit value of H 0 in the context of all these models is not consistent with the local determination of H 0 shown in Equation (2). • The best-fit value of the SnIa absolute magnitude M is not consistent with the value of M determined by the local Cepheid calibrators. • The quality of fit of all these parameterizations becomes significantly worse when the local determination of the H 0 point is included.
We conclude that these dark energy parameterizations are unable to resolve the Hubble tension. Another argument adding to the inability of quintessence models to ease the Hubble tension can be found in [89]. It is stated that, being a w > −1 model (excluding the cosmological constant), it can only make the Hubble tension worse. In general, late time models trying to resolve the cosmological tensions do not succeed, as can also be seen in [7], where it is shown that late time approaches worsen the growth tension and in [6,90] where they worsen the fit in comparison with the ΛCDM, or even the Hubble tension when BAO data are used. However, since the Hubble tension also manifests itself as an M-tension [91,92], a transition of M can be used for these models to fully resolve the Hubble tension, as shown in Figure 1 and discussed in detail in [20,93].
In addition to finding the best-fit form of H(z) corresponding to these parametrizations, we have reconstructed the quintessence Lagrangian that would reproduce the observed form of H(z). Due to the phantom nature of the best-fit dark energy equation-of-state parameter, we found that the reconstructed Lagrangian has a negative kinetic term and is therefore plagued with instabilities. A possible extension of this work is the consideration of more parameterizations of the equation-of-state parameter w(z) and the comparison of the quality of fit provided by each one, using the most recent cosmological data. A few interesting parameterizations include the linear [94]: w(z) = w 0 + w a z, the Alcaniz-Barbosa [95]: w(z) = w 0 + w a z(1+z) 1+z 2 , the sqrt [83]: w(z) = w 0 + w a z √ 1+z 2 , the Sine [96]: w(z) = w 0 + w a sin z and the Jassal-Bagna-Padmanabhan [97]: w(z) = w 0 + w a z (1+z) 2 . Another interesting extension of the present analysis is the consideration of modified gravity scalar tensor Lagrangians, which can accommodate more naturally the phantom behavior of the Hubble expansion without instabilities [98].

Data Availability Statement:
The SNIa data used in this paper are from the Pantheon dataset, also used in the analysis [5] and can be found in the Github repository [99] and on the website [100]. The Pantheon dataset consists of 1048 data points coming from six different probes and covering the redshift range 0.01 < z < 2.3. This publicly available dataset contains the name of each SNIa, its redshift both in the CMB and the heliocentric frame and the observed apparent magnitude m obs along with its corresponding error σ m obs . The apparent magnitude m obs is reported after the application of a K-correction and a correction in terms of the light curve width along with some other corrections due to biases from simulations of the SNIa [42]. The BAO observational data that we will use in this analysis are the data from 6dFGS and WiggleZ found in [101], from SDSS (Data Release 7 and 11) found in [102,103] and from Ly-α measurements found in [104]. For the CMB observational data, we use the central values for R,l A and Ω b h 2 from Planck 2018 for a flat universe. The data can be found in [33] and are expressed in terms of a data vector and a covariance matrix. Lastly, the numerical data files for the reproduction of the figures can be found in the M_crisis Github repository under the MIT license.

Acknowledgments:
We thank Lavrentios Kazantzidis and George Alestas for useful discussions and their help with the code.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results

Abbreviations
The following abbreviations are used in this manuscript:

. FRW Metric
The cosmic metric is well approximated by the Friedmann-Robertson-Walker metric [105,106]: The curvature parameter k can be −1, 0, 1. • a(t) is the cosmic scale factor. • The cosmic time t is the proper time measured by a free-falling observer. • The coordinates r, θ, φ are comoving coordinates.
Another form of the FRW metric is: where: The time is again the proper time measured by a free-falling observer, while χ is a new radial coordinate (in which r, θ, φ is a comoving coordinate).
The cosmic dynamics is determined by the Friedmann equations Friedmann Equation: ȧ a

Friedmann Acceleration Equation:
a where ρ and p correspond to the sum of all contributions to the energy density and pressure of the contents of the universe. The Friedmann equation in terms of the Hubble parameter takes the form: The critical density is defined for a given value of the Hubble parameter as The quantities evaluated at present are denoted by a subscript '0'. For example, the age of the Universe evaluated today is expressed as t = t 0 and is the critical density today. The dimensionless density parameters are defined as follows: where the index I corresponds to a given density component (e.g., matter, radiation or dark energy). Thus, the curvature density parameter is Using the stress-energy tensor of a perfect fluid as seen by a comoving observer: and the covariant conservation equation for it: the continuity equation is then obtained aṡ ρ + 3ȧ a (ρ + p) = 0 (A14)

Constant Equation of State
Introducing the constant equation of state: w = p/ρ, which can be used to parametrize most cosmological fluids, it is easy to show that ρ scales as: In the case of cold dark matter (w = 0), radiation (w = 1/3) and vacuum energy or cosmological constant Λ (w = −1), we obtain: Figure A1. Evolution of the energy densities in the Universe. Adapted from: [105].

Dynamical Equation of State
Introducing a time dependent equation of state: in Equation (A14) we obtain: which by dividing both sides of the equation by ρ(t) takes the form: Thus: dρ which we can integrate as follows: Here, we can make the substitution: which corresponds to the new differential: Thus: Finally we obtain: There are many ways to choose two components to study, but an interesting choice is a universe where matter and radiation densities are comparable.
We denote the value of the scale factor when matter and radiation were equally important as: which was shortly before the cosmic microwave background was released. We can write the total energy density at that time in the form: To help with the calculations, we introduce conformal time: Using the above the Friedmann Equations, take the form: with primes denoting the derivatives with respect to conformal time.
Solving for a, we obtain: where: The general form of the Friedmann equations for multiple components is: The equations of state are: p i = w i ρ i , where w i is the dimensionless equation-of-state parameter of the i-th component. Now, considering a universe with matter, radiation, curvature and a cosmological constant (Λ), we obtain: while the density parameters at the present time satisfy the relation: Appendix A.2. Redshift z Figure A2. Diagram showing two null (photon) geodesics. In our case where the photons are emitted radially, we have θ E = θ R and φ E = φ R . Adapted from [106].
Consider a photon emitted (radially) from a distant galaxy at time t = t E and received by us on earth at t = t R . Along the photon path, we have ds = dθ = dφ = 0, so from the FRW metric for a flat universe Equation (A2), we obtain: Since we have an incoming photon: The photon is emitted at a comoving distance χ E and is observed here on Earth at a comoving distance χ R = 0: Now, consider another photon emitted at time t = t E + δt E and received at t = t R + δt R . For that photon: Combining Equations (A41) and (A42), we obtain: Adding t E +δt E t E cdt a(t) in both sides of the equation: which gives: If δt R is really small, we can consider a(t) to be constant. Therefore, the above equation gives: Thus, for the redshift, we obtain: which by using: c = lν takes the form: Finally: In cosmology we define the redshift parameter z as the fractional shift in wavelength of a photon emitted by a distant galaxy at time t and observed on Earth today at t 0 . Thus, by replacing t R with t 0 and t E simply with t, we obtain: which gives: using that: a(t 0 ) = a 0 = 1.
To better understand the concept of redshift and scale factor, consider the following example. When we observe a galaxy with a redshift z = 2, we are observing it as it was at time t = t E when the Universe had a scale factor a(t E ) = 1/3, and thus, it was 1/3 of its present size. According to Hubble's Law, named after Edwin Hubble, who discovered the expansion of the Universe in 1929, galaxies appear to recede from us with a recession speed proportional to their distance from us [106]: The proportionality constant H 0 is called the Hubble constant and measures the current expansion of the Universe. which, by evaluating the integral, gives:

Metric Distance
The distance multiplying the solid angle dΩ 2 = dθ 2 + sin 2 θdφ 2 is the metric distance: which means that for a flat universe, the metric distance is simply equal to the comoving distance χ: distance is not observable but is useful in defining observable distances. Figure A3. Figure that helps with the definition of angular diameter distance. Adapted from: [105] Comoving/Coordinate Distance The comoving distance between two objects in the Universe is the distance between them that remains constant with the epoch, when they are both moving only with the Hubble flow, i.e., they move solely due to the expansion of the Universe.

Appendix A.4.2. Observable Distances
Physical Distance It is the actual proper distance r p between two objects in the Universe that can be measured by a physical ruler. It is related to the comoving distance d co (z) through:

Angular Diameter Distance
Let us assume that we observe an object at a comoving distance χ and that the photons that we observe today were emitted at time t 1 . Assuming that the object has a known physical size D and someone on Earth measures its angular size to be δθ, we can define the angular diameter distance as: which for small δθ 1 is the Euclidean formula for its distance. Through the FRW metric, we can find the relation between the physical transverse size of the object and its angular size on the sky: Hence: The angular diameter distance is really useful because of objects of known physical size D, which we call "standard rulers" as, for example, the fluctuations in the CMB.

Luminosity Distance
Let us consider a luminous cosmological object at a fixed comoving distance χ, which emits at an absolute luminosity L. In a static Euclidean space, the relation between the absolute luminosity and the observed flux F would be: since the power radiated by the luminous object is distributed in the spherical surface with radius χ. In an FRW spacetime, however, the result must be modified. This happens because: • At the time t 0 that we observe the light from the object, the proper area of a sphere drawn around a supernova and passing through the Earth is 4πd 2 m .
• The rate at which we detect photons from the object is reduced compared to the rate that they are emitted, by the redshift factor: The energy of the photons is also being redshifted, so the energy that we observe them to have is reduced compared to the one they had when they were emitted by the same redshift factor: 1 + z. Thus, we obtain the formula: We define the luminosity distance to be: Hence, we find that: The luminosity distance proves to be very useful because of objects called standard candles, which are objects of "known" absolute luminosity, or more appropriately of an absolute luminosity that we can estimate independently of their distance and apparent luminosity. Such objects are variable stars called cepheids or a special type of supernovae (SN) called Type Ia supernovae.

Distance Modulus
The apparent magnitude m of an astronomical object is defined as: where F is the apparent flux of the object, and F re f is a reference flux.
The absolute magnitude M is defined as the apparent magnitude the object would have if it was 10 pc away from the observer: The distance modulus µ is defined as the difference between them and after some calculations can be expressed as: where d L is the luminosity distance.
If the distance is expressed in Mpc, then we can write µ in the form: This kind of distance is very useful when we use data from supernovae.

Appendix B. Theoretical Background
As we discussed earlier, in 1929, Hubble showed that the Universe is expanding, as every other galaxy appeared receding from us with a recession speed analogous to its distance. The proportionality constant is H 0 = 100 h kms −1 Mpc −1 . We have already talked about observable distances, and thus, we can identify that the distance that was measured and used in Hubble's Law is actually the physical distance r p , so: However, some later results taken by probing supernovae in high redshifts showed that this linear relation does not hold anymore, and it seems that the expansion is speeding up. Thus, the latest data imply that we live not only in an expanding universe, but in a universe with an accelerating expansion.
Appendix B.1. Dark Energy In a universe described by general relativity and that is matter-dominated (as it was thought to be by cosmologists in the last century), one would expect the expansion to be slowing down due to the influence of gravity. Therefore, the second derivative of the expansion was named the deceleration parameter: Figure A6. The remnants of Supernova 1604 or Kepler's Supernova. It was a Type Ia supernova that occurred in our galaxy in the constellation of Ophiuchus. It was observable by the naked eye and was named by Johannes Kepler who described it in: De Stella Nova. Adapted from: [110]. q(z) ≡ −ä aH 2 (A90) We can see that matter can only lead to decelerating expansion using the acceleration equation and its equation of state:ä As does radiation:ä since a(t), ρ m and ρ r are positive. In the context of general relativity, the only way we can obtain an accelerating expansion is by assuming that there is an additional component in our Universe called "dark energy". In order to have an equation that describes the universe as it is, we need dark energy to have [111]: • A positive energy density (ρ X > 0), assuming that the universe is flat. • A negative pressure (p X < 0), which can cancel out gravity and potentially lead to accelerating expansion.
Adding the dark energy term in the Friedmann acceleration equation, yields: with: p X = wρ X (A96)

Incomplete Gamma Function
The incomplete gamma function is defined using the gamma function as: x 0 e −t t a−1 dt (A103) for a > 0. It is also common to use Q(a, x), the compliment of P(a, x), which is also called an incomplete gamma function: x n (A109)

Error Function
The error function is a special case of the incomplete gamma function and is defined as: The error function has the following properties: The probability density function p(x) for a Gaussian distribution with mean µ and standard deviation σ is: The variance of this distribution is σ 2 .
for ∆χ 2 = 9.72 and ν = 4: A useful way to obtain these values is through Mathematica by using the code in Figure  A9. det(e A ) = e a 0 · e a 1 · · · e a n = e a 0 +a 0 +...a n = e tr(A) (A122) Assuming a new matrix B with: In this case we have: B = g ab ⇒ B −1 = g ab and det B = g. Thus: dg g = g ab dg ab ⇒ dg = gg ab dg ab (A126)

Appendix E. Derivation of the Klein-Gordon Equation
The general scalar field action in Riemann spacetime is: In a spacetime with signature (+, −, −, −), the Lagrangian density of φ is: For any region Ω, we consider variation of the field: which vanishes on the surface Γ(Ω), bounding the region Ω.