1 Introduction↩︎

Observations from the Wilkinson Microwave Anisotropy Probe (WMAP) [1], [2], the Supernova Cosmology Project [3], [4], the Sloan Digital Sky Survey (SDSS) [5]–[7], the Planck mission [8], and X-ray studies [9] consistently indicate that the Universe is currently undergoing accelerated expansion. To explain this phenomenon, a hypothetical component called Dark Energy (DE), with a negative pressure, has been introduced. The Cosmological Constant \(\Lambda_{CC}\) is the simplest candidate, yet it faces conceptual challenges such as the cosmological constant and coincidence problems, motivating the development of alternative dynamical models [10]–[12].

From a theoretical perspective, \(\Lambda_{CC}\) can be incorporated as a constant term in Einstein’s field equations. Quantum Field Theory calculations using Planck- or electroweak-scale cut-offs predict vacuum energy densities vastly larger than observed, by factors of \(10^{123}\) and \(10^{55}\), respectively. The lack of a natural symmetry to suppress \(\Lambda_{CC}\) gives rise to the cosmological constant problem, while the coincidence problem questions why matter and DE densities are comparable today [10]–[13].

Within the standard cosmological model, DE accounts for approximately two-thirds of the current energy density \(\rho_{\text{tot}}\) [14], with the remaining fraction composed mainly of Dark Matter (DM) and baryons. Despite precise measurements, the microscopic origin of DE remains elusive.

Dynamical DE models, where the equation of state parameter \(\omega_D\) evolves over time, provide a flexible framework consistent with observational data. Examples include scalar field scenarios such as quintessence [15]–[17], k-essence [18]–[20], tachyon [21]–[23], phantom [24]–[26], dilaton [27], [28], and quintom models [29]–[32]. Interacting DE models, including those based on the Chaplygin gas [33]–[35], the Agegraphic Dark Energy (ADE) and the New ADE (NADE) [36], [37], have also been explored.

Holographic approaches provide a complementary perspective, based on the principle that the entropy of a system scales with its boundary area rather than volume [38]–[40]. Holographic Dark Energy (HDE), first proposed by Li [41], postulates a DE density \[\rho_D= 3\alpha M_p^2 L^{-2},\] where \(\alpha\) is a dimensionless constant and \(M_p = (8\pi G_N)^{-1/2}\) is the reduced Planck mass. Cohen et al. [42] initially suggested that the vacuum energy should be bounded to prevent black hole formation, but the naive choice \(\rho_\Lambda \propto H^2\) fails to drive acceleration. Using the future event horizon as the infrared cutoff, HDE models successfully reproduce late-time cosmic acceleration [41]. Extensions include the Holographic Ricci Dark Energy, with \(L \propto R^{-1/2}\) [43], and the Granda–Oliveros formulation, where \(\rho_\Lambda\) depends on both \(H\) and \(\dot{H}\) [44]–[46]. These models have been extensively confronted with supernova, CMB, and BAO observations [47]–[55] and are reviewed in [56]–[76].

We now want to study two different Dark Energy energy density models which involve higher time derivatives of the Hubble parameter \(H\).
The first one is given by: \[\begin{align} \rho_D = 3 \left[ \alpha \left(\frac{\dddot{H}}{H^2}\right) +\beta \left(\frac{\ddot{H}}{H}\right) + \gamma \dot{H}+\delta H^2\right] \label{cz1} \end{align}\tag{1}\] where the quantities \(\alpha\), \(\beta\), \(\gamma\) and \(\delta\) are dimensionless parameters. For mathematical simplicity, we set the reduced Planck mass to unity, i.e. \(M_p = 1\). We note that the inverse of the Hubble parameter squared \(H^{-2}\) and the inverse of the Hubble parameter \(H^{-1}\) are included in the first and second terms to guarantee consistency of physical dimensions across all terms.

The cosmological behavior and main features of this DE model depend crucially on the four parameters of the model. The energy density given in Eq. (1 ) can be regarded as a generalization of several previously proposed DE models. For instance, by setting \(\alpha = 0\), we recover the energy density introduced in Chen & Jing [77] and studied in other subsequent works. Instead, in the limiting case of \(\alpha = \beta=0\), the energy density given in Eq. (1 ) reduces to the energy density of DE with the Granda-Oliveros cut-off [78]. Furthermore, for the particular choice \(\alpha=\beta=0\), \(\gamma=1\), and \(\delta=2\), the model reproduces the dark energy (DE) density with infrared (IR) cut-off determined by the average radius associated with the Ricci scalar, valid in a spatially flat Universe (\(k=0\)). Since the present model introduces an additional free parameter, it provides a more general framework than the Ricci Dark Energy (RDE) scenario. Similar DE models have been investigated in detail in [79]–[81].
The general expression of the second model we consider is given by: \[\begin{align} \rho_D = 3 \left[ \alpha \left(\frac{\dddot{H}}{H^2}\right) - \zeta\left( \frac{\ddot{H}\dot{H}}{H^3} \right) +\beta \left(\frac{\ddot{H}}{H}\right) + \gamma \dot{H}+\delta H^2\right] \label{cz2} \end{align}\tag{2}\] where \(\alpha\), \(\beta\), \(\gamma\), \(\delta\) and \(\zeta\) are five constant parameters. In the limiting case of \(\zeta=0\), we recover the first model introduced.
Moreover, in the limiting case of \(\alpha=\zeta=0\), we recover the model studied in Chen & Jing [77], while for \(\alpha=\beta=\zeta=0\) we obtain the HDE energy density model with Granda-Oliveros cut-off. Furthermore, for \(\alpha=\beta=\zeta=0\), \(\gamma=1\) and \(\delta=2\) we obtain the HDE model with cut-off given by the average radius of the Ricci scalar for a spatially flat Universe.
In the following sections, we derive key cosmological quantities for a power-law form of the scale factor in these two models, including the dark energy density, dark energy pressure, the equation of state parameter, the deceleration parameter, and the evolution of the dark energy fractional density. Additionally, we compute the age of the Universe for selected values of the relevant parameters.
The paper is structured as follow. In Section II, we study the first model introduced. In Section III, we study the second model we consider. In Section IV, we evaluate the present age of the Universe for both models, exploring different sets of model parameters. Finally, in Section V, we write the conclusions of this work.

2 Holographic Model in a Non-Flat Universe↩︎

In this Section, we outline the main features of the first Dark Energy (DE) model under investigation and derive several fundamental cosmological quantities.
We consider a Universe described by a spatially homogeneous and isotropic Friedmann–Lemaitre–Robertson–Walker (FLRW) spacetime, whose line element can be written as: \[\begin{align} ds^2 = -dt^2 + a^2(t)\left[\frac{dr^2}{1 - kr^2} + r^2 (d\theta^2 + \sin^2\theta\, d\varphi^2)\right], \label{6} \end{align}\tag{3}\] where \(t\) denotes the cosmic time, \(a(t)\) is the scale factor describing the expansion history of the Universe, \(r\) is the comoving radial coordinate, and \(\theta\) and \(\varphi\) are the usual angular coordinates in spherical symmetry, with \(0 \leq \theta \leq \pi\) and \(0 \leq \varphi < 2\pi\). The constant \(k\) represents the spatial curvature parameter and can take the values \(-1\), \(0\), or \(+1\), corresponding to open, flat, and closed geometries, respectively.

The evolution of a homogeneous and isotropic Universe within the framework of General Relativity is determined by the Friedmann equations, which, in the presence of both dark energy (DE) and dark matter (DM), can be written as: \[\begin{align} H^2 + \frac{k}{a^2} &=& \frac{1}{3M^2_p}\left( \rho_D + \rho_m \right), \tag{4} \\ \dot{H} + 2H^2 + \frac{k}{a^2} &=& \frac{8\pi G}{6}\, p_D, \tag{5} \end{align}\] where \(H = \dot{a}/a\) is the Hubble parameter, \(\rho_D\) and \(p_D\) denote the energy density and pressure of dark energy, respectively, and \(\rho_m\) is the energy density of pressureless dark matter (i.e., \(p_m = 0\)).
The fractional energy densities of the matter component, the dark energy component, and the spatial curvature are defined as: \[\begin{align} \Omega_m &=& \frac{\rho_m}{\rho_{\text{cr}}} = \frac{\rho_m}{3M^2_p H^2}, \tag{6} \\ \Omega_D &=& \frac{\rho_D}{\rho_{\text{cr}}} = \frac{\rho_D}{3M^2_p H^2}, \tag{7} \\ \Omega_k &=& \frac{k}{a^2 H^2}, \tag{8} \end{align}\] where \(\rho_{\text{cr}}\) is the critical energy density required for a spatially flat Universe, and is given by: \[\begin{align} \rho_{\text{cr}} = 3M^2_p H^2. \end{align}\]

Using the results of Eqs. (6 ) and (7 ), the Friedmann equation given in Eq. (4 ) can be rewritten in the following dimensionless form: \[\begin{align} \Omega_m + \Omega_D + \Omega_k = 1. \label{11} \end{align}\tag{9}\] In order to ensure the validity of the Bianchi identity, or equivalently the local conservation of energy-momentum, the total energy density \(\rho_{\text{tot}} = \rho_D + \rho_m\) must satisfy the continuity equation: \[\begin{align} \dot{\rho}_{\text{tot}} + 3H\left( \rho_{\text{tot}} + p_{\text{tot}} \right) = 0, \label{12old} \end{align}\tag{10}\] where \(\rho_{\text{tot}}\) and \(p_{\text{tot}}\) denote the total energy density and total pressure of the cosmic fluid, respectively, and are defined by: \[\begin{align} \rho_{\text{tot}} &=& \rho_m + \rho_D, \\ p_{\text{tot}} &=& p_D, \end{align}\] since DM is assumed to be pressureless.
The continuity equation written in Eq. (10 ) can alternatively be expressed in terms of the total equation of state (EoS) parameter \(\omega_{\text{tot}} = p_{\text{tot}} / \rho_{\text{tot}}\) as: \[\begin{align} \dot{\rho}_{\text{tot}} + 3H\left(1 + \omega_{\text{tot}} \right) \rho_{\text{tot}} = 0. \label{12} \end{align}\tag{11}\]

Since the energy densities of DM and DE are assumed to be conserved separately, we have that Eq. (11 ) can be decomposed into two independent continuity equations. In the non-interacting case, these read: \[\begin{align} \dot{\rho}_D + 3H\left( 1 + \omega_D \right)\rho_D &=& 0, \tag{12} \\ \dot{\rho}_m + 3H\rho_m &=& 0. \tag{13} \end{align}\] Using the general definition of the EoS parameter \(w\) for DE: \[\begin{align} \omega_D = \frac{p_D}{\rho_D}, \label{mammud} \end{align}\tag{14}\] we can equivalently express Eq. (12 ) in the following form: \[\begin{align} \dot{\rho}_D + 3H\left( p_D + \rho_D \right) = 0. \label{12de} \end{align}\tag{15}\]

Furthermore, Eqs. (12 ), (13 ) and (15 ) can be rewritten in terms of the new variable \(x = \ln a\), using the relation \(\frac{d}{dt} = H \frac{d}{dx}\). Defining derivatives with respect to \(x\) by a prime, we obtain: \[\begin{align} \rho'_D + 3\left( 1 + \omega_D \right)\rho_D &=& 0, \tag{16} \\ \rho'_D + 3\left( p_D + \rho_D \right) &=& 0, \tag{17} \\ \rho'_m + 3\rho_m &=& 0. \tag{18} \end{align}\]

Since dark energy (DE) accounts for roughly two-thirds of the present-day total energy density of the Universe, while its contribution was essentially negligible in the early stages after the Big Bang, it is reasonable to postulate that the DE density evolves with the expansion of the Universe. In this context, it is natural to explore DE models in which the energy density depends on the Hubble parameter \(H\) and its derivatives with respect to the cosmic time \(t\), given that \(H\) encapsulates the expansion rate of the Universe.

3 First Model↩︎

We not study the first model of energy density considered in this paper.
In this work, we consider a power-law form of the scale factor \(a(t)\) given by the following relation: \[a(t) = b_0 t^n, \label{scale}\tag{19}\] where \(b_0\) and \(n\) represent two positive constants. We will make some considerations about the value of \(b_0\) later in the paper.
We now want to calculate some quantities using the expression of the scale factor given in Eq. (19 ).
Using Eq. (19 ) in the expression of \(\rho_D\) given in Eq. (1 ), we obtain the following expression for \(\rho_D\) as function of the time: \[\begin{align} \rho_{D,1}(t) = \frac{3}{t^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right).\label{cz3} \end{align}\tag{20}\] We now want to calculate the expression of the EoS paramater for DE.
From the continuity equation given in Eq. (12 ), we find the general expression: \[\begin{align} \omega_{D,1} = -1 - \frac{\dot{\rho}_{D,1}}{3H\rho_{D,1}} \end{align}\]

From Eq. (20 ), we obtain: \[\begin{align} \dot{\rho}_{D,1}(t) = -\frac{6}{t^3} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right). \end{align}\] Moreover, we have that: \[\begin{align} H(t) = \frac{n}{t} \end{align}\] Therefore, we obtain the following relation for \(\omega_D\): \[\begin{align} \omega_{D,1} = -1+\frac{2}{3n} \label{eos} \end{align}\tag{21}\] This expression implies that \(\omega_D\) is always greater than \(-1\), indicating a quintessence-like behavior. In particular, for small values of \(n\), the EoS parameter deviates significantly from the cosmological constant limit, reflecting a less dominant dark energy component. For \(n=1\), the equation of state (EoS) parameter of dark energy is \[\begin{align} \omega_{D,1} = -\frac{1}{3}. \end{align}\] This value corresponds to the critical limit between decelerated and accelerated expansion of the Universe:

• If \(\omega_{D,1} > -1/3\) (i.e. \(< n <1\)), the expansion is decelerating.

• If \(\omega_{D,1} < -1/3\), (i.e. (\(n>1\))) the expansion is accelerating.

Thus, for \(n=1\), the model predicts a marginal expansion at the boundary between acceleration and deceleration.

As \(n\) increases, \(\omega_D\) approaches \(-1\), effectively mimicking a cosmological constant. Therefore, the parameter \(n\) controls the deviation of dark energy from a pure \(\Lambda\)-like behavior, with higher values of \(n\) corresponding to an accelerated expansion that closely resembles a \(\Lambda\)CDM scenario. Overall, the model predicts a non-phantom, accelerating universe driven by quintessence-like dark energy.
In Fig. (1) we plot the behavior of \(\omega_{D,1}\) obtained in Eq. (21 ) in the range of values of \(n \in \left[0.1-10\right]\).

We now derive the expression of the pressure of DE \(p_{D,1}\).
From the continuity equation for DE, we obtain the following general expression for \(p_{D,1}\): \[\begin{align} p_{D,1}= -\rho_{D,1} - \frac{\dot{\rho}_{D,1}}{3H} \end{align}\] Using the expressions of \(\rho_{D,1}\), \(\dot{\rho}_{D,1}\) and \(H\) we derived before, we can write: \[\begin{align} p_{D,1}(t) &=& \left( \frac{2}{3n} -1 \right)\rho_{D,1}(t)\nonumber \\ &=& \left( \frac{2}{n} -3 \right)\left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)t^{-2} \end{align}\]

We now want to obtain the expression of the Hubble parameter squared as function of the time. For the energy density of DM, from the continuity equation obtain in Eq. (13 ), we obtain: \[\begin{align} \rho_m= \rho_{m0}a^{-3} \end{align}\] Using the expression of \(a(t)\) we have chosen in this paper, we can write: \[\begin{align} \rho_m&=& \rho_{m0}b_0^{-3}t^{-3n}\nonumber \\ &=&\frac{\rho_{m0}}{b_0^{3}t^{3n}} \end{align}\] which leads to \[\begin{align} \Omega_m&=& \Omega_{m0}b_0^{-3}t^{-3n}\nonumber \\ &=&\frac{\Omega_{m0}}{b_0^{3}t^{3n}} \end{align}\] Generally, for the curvature energy density, we have: \[\begin{align} \rho_k= ka^{-2} \end{align}\] Using the expression of the scale factor we consider in this paper, we can write: \[\begin{align} \rho_k&=& k\cdot b_0^{-3}t^{-3n}\nonumber \\ &=&\frac{k}{b_0^{3}t^{3n}} \end{align}\] which leads to: \[\begin{align} \Omega_k&=& \Omega_{k0}\cdot b_0^{-3}t^{-3n}\nonumber \\ &=&\frac{\Omega_{k0}}{b_0^{3}t^{3n}} \end{align}\] We then obtain the following expression for \(H^2(t)\): \[\begin{align} H_1^2(t)= H_0^2\left\{ \frac{\Omega_{m0}}{b_0^{3}t^{3n}}+\frac{\Omega_{k0}}{b_0^{2}t^{2n}}+\frac{3}{H_0^2t^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)\right\} \end{align}\]

We now want to write the equations we derived as functions of the redshift \(z\).
The energy density of DE can be written as: \[\begin{align} \rho_{D,1}(z) \;=\; 3\, \left(-\frac{6\alpha}{n}+2\beta-n\gamma+\delta n^2\right)b_0^{2\over n}\,(1+z)^{2\over n} \end{align}\]

The pressure of DE \(p_{D,1}\) is given by: \[\begin{align} p_{D,1}(z) &=& \left(\frac{2}{3n}-1\right) \rho_D(z) \nonumber\\ &=& \left(\frac{2}{n}-3\right) \left(-\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right) b_0^{2\over n} (1+z)^{2\over n}. \end{align}\]

The Hubble parameter squared as function of \(z\) is given by: \[\begin{align} H_1^2(z)= H_0^2\left\{ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2+ \frac{3}{H_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{2\over n}\right\}\label{unical1} \end{align}\tag{22}\] Finally, we can obtain the deceleration parameter \(q\) from the relation: \[\begin{align} q_1(z)= -1+\left(\frac{1+z}{2}\right)\frac{1}{h_1^2(z)}\frac{dh_1^2(z)}{dz} \end{align}\] where \(h_1^2=H_1^2/H_0^2\) From the result of Eq. (22 ), we can write: \[\begin{align} \frac{dh_1^2(z)}{dz} = 3\Omega_{m0}(1+z)^2+2\Omega_{k0}(1+z)+ \frac{6}{nH_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{\frac{2}{n}-1} \end{align}\] Therefore, the final expression of \(q_1(z)\) is given by: \[\begin{align} q_1(z)&=&-1+\frac{1}{2}\cdot \left[ 3\Omega_{m0}(1+z)^3+2\Omega_{k0}(1+z)^2+ \frac{6}{nH_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{\frac{2}{n}} \right]\times \nonumber \\ &&\left[ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2+ \frac{3}{H_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{\frac{2}{n}} \right] \end{align}\]

We can now make some considerations about \(b_0\).
Considering the powe law form of the scale factor we defined in Eq. (19 ), we obtain that the Hubble parameter can be written as follows: \[\begin{align} H(t) = \frac{n}{t}. \end{align}\] If we adopt the conventional normalization \(a_0=1\), then the present time \(t_0\) satisfies \[\begin{align} 1 = a_0 = b_0 t_0^n \;\;\Rightarrow\;\; t_0 = b_0^{-1/n}. \end{align}\] Therefore, \[\begin{align} H_0 = H(t_0) = \frac{n}{t_0} = n\,b_0^{1/n}. \label{unical2} \end{align}\tag{23}\]

From Eq. (23 ), we obtain \[\begin{align} b_0=\left(\frac{H_0}{n}\right)^n \rightarrow b_0^{2\over n} = \left(\frac{H_0}{n}\right)^{2} \end{align}\]

Using the result we obtained for \(b_0^{2\over n}\), the energy density of DE assumes the following form: \[\begin{align} \rho_{D,1}(z) &=& 3\left(-\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2\right) \left(\frac{H_0}{n}\right)^{2}(1+z)^{2\over n}\nonumber \\ &=&3\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right) H_0^{2}(1+z)^{2\over n} \end{align}\]

The Hubble parameter squared can be written as: \[\begin{align} H_1^2(z)= H_0^2\left\{ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2+3 \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right) \,(1+z)^{2\over n}\right\} \end{align}\] while the pressure of DE can be rewritten as: \[\begin{align} p_{D,1}(z) &=& \left( \frac{2}{3n}-1\right) \rho_D(z) \nonumber\\ &=& \left( \frac{2}{n}-3\right) \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)H_0^2(1+z)^{2\over n}. \end{align}\]

Finally, the deceleration parameter is given by: \[\begin{align} q_1(z)&=&-1+\frac{1}{2}\cdot \left[ 3\Omega_{m0}(1+z)^3+2\Omega_{k0}(1+z)^2+ \frac{6}{n} \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\,(1+z)^{\frac{2}{n}} \right]\times \nonumber \\ &&\left[ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2+ 3 \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\,(1+z)^{\frac{2}{n}} \right]^{-1} \end{align}\]

The present day value of \(q_1\) is given by: \[\begin{align} q_{1,0}&=&-1+\frac{1}{2}\cdot \left[ 3\Omega_{m0}+2\Omega_{k0}+ \frac{6}{n}\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right) \right]\times \nonumber \\ &&\left[ \Omega_{m0}+\Omega_{k0}+ 3 \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right) \right]^{-1} \end{align}\]

We now want to find the evolutionary form of the fractional energy density of DE for this model.
Using the expression of \(\rho_{D,1}(z)\), we can write: \[\begin{align} \Omega_{D,1}(z) &=& \frac{\rho_{D,1}(z)}{3H_0^2}\nonumber\\ &=&\left[-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{2\over n}. \end{align}\] Therefore, the evolutionary form of \(\Omega_{D,1}\) is given by: \[\begin{align} \Omega'_{D,1}(z) &=&\frac{2}{n}\left[-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{\frac{2}{n}-1}. \end{align}\] Moreover, we find that the present day value of \(\Omega_D\) is given by: \[\begin{align} \Omega_{D,1,0} &=&\left[-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right]. \end{align}\] We can obtain some hints about the values of the parameters of the model using the value of \(\Omega_{D,1,0}\).
We now that at present time \(\Omega_{D,1,0}\approx 0.685\), therefore we should have: \[\begin{align} -\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \approx 0.685 \end{align}\] Then, once the value on \(n\) is fixed, we obtain a constraint on the values of \(\alpha\), \(\beta\), \(\gamma\) and \(\delta\).
For example, if we choose \(n=2\), we obtain the following set of values: \(\alpha = 1\), \(\beta = 2\), \(\gamma = 0.5\) and \(\delta = 0.685\). Another possible set is: \(\alpha = 0.5\), \(\beta = 1\), \(\gamma = 0.2\) and \(\delta = 0.66\).
If we consider the case with \(n=3\), we obtain the following set of values: \(\alpha = 1\), \(\beta = 2\), \(\gamma = 0.5\) e \(\delta = 0.63\).
If we consider \(n=4\), we obtain the following set of values \((\alpha,\beta,\gamma,\delta)=(1,1,1,0.90375)\).

In Table 1, we have the values of \(q_{1,0}\) for the different values of \(n\) we considered and for the different values of the parameters we derived. Moreover, we considered \(\Omega_{m0}=0.315\) and \(\Omega_{k0}=0.01\).

3.1 Interacting Case↩︎

We now consider the presence of interaction between Dark Sectors.
We now extend our analysis by considering the possibility of an interaction between the dark sectors. This idea refers to scenarios where dark matter (DM) and dark energy (DE) are not entirely independent but may exchange energy or momentum. Such a coupling is often motivated by attempts to address the so-called coincidence problem, namely why the energy densities of DM and DE are of the same order today despite their different evolutionary histories. Allowing for an interaction modifies the standard cosmological dynamics and can leave distinctive observational imprints, such as changes in the expansion history, deviations in structure formation, or shifts in the cosmic microwave background (CMB) anisotropies. These models have therefore been widely studied as possible alternatives or extensions to the concordance \(\Lambda\)CDM framework.

In the presence of such a coupling, the conservation equations for DE and DM are modified as follows: \[\begin{align} \dot{\rho}_{D} + 3H \rho_{D} (1+\omega_{D}) &=& -Q, \tag{24}\\ \dot{\rho}_{m} + 3H \rho_{m} &=& Q, \tag{25} \end{align}\] where \(Q\) specifies the rate of energy transfer between the two sectors. In general, \(Q\) may be a function of several cosmological quantities, including the Hubble parameter \(H\), the deceleration parameter \(q\), and the energy densities \(\rho_{m}\) and \(\rho_{D}\), i.e. \(Q=Q(\rho_{m},\rho_{D},H,q)\). A variety of choices for this function have been considered in the literature. In our study, we adopt the phenomenological form \[Q = 3 d^{2} H \rho_{m}, \label{Q}\tag{26}\] where \(d^{2}\) is a dimensionless constant quantifying the strength of the interaction, often called the transfer rate or coupling parameter [82]–[84].

Observational analyses combining different cosmological probes — such as the Gold SNe Ia sample, CMB data from WMAP, and BAO measurements from SDSS — suggest that \(d^{2}\) should be positive and small. This outcome is in agreement with the requirements imposed by the cosmic coincidence problem as well as with thermodynamical considerations [85]. Additional constraints from CMB anisotropy studies and galaxy cluster observations further indicate the range \(0 < d^{2} < 0.025\) [86]. More generally, the parameter is usually considered within \([0,1]\), with the special case \(d^{2}=0\) reducing to the standard non-interacting FRW cosmology. It is worth stressing that many other functional forms of \(Q\) have been proposed in the literature, each leading to different phenomenological consequences.

The expression of the energy density of DE is the same as in the non-interacting case.
For the energy density of DM, from the continuty equation for the interacting case, we obtain the following expression of \(\rho_m\) as function of the scale factor: \[\begin{align} \rho_{m,I}= \rho_{m0}a^{-3(1-d^2)} \end{align}\] Using the expression of \(a(t)\) we have chosen in this paper, we can write: \[\begin{align} \rho_{m,I}(t)&=& \rho_{m0}b_0^{-3(1-d^2)}t^{-3n(1-d^2)}\nonumber \\ &=&\frac{\rho_{m0}}{b_0^{3(1-d^2)}t^{3n(1-d^2)}} \end{align}\] from which we obtain: \[\begin{align} \Omega_{m,I}(t)&=& \Omega_{m0}b_0^{-3(1-d^2)}t^{-3n(1-d^2)}\nonumber \\ &=&\frac{\Omega_{m0}}{b_0^{3(1-d^2)}t^{3n(1-d^2)}} \end{align}\]

The expression of \(\rho_{m,I}\) for the interacting case as function of the redshift is given by \[\begin{align} \rho_{m,I}(z)&=& \rho_{m0}(1+z)^{3(1-d^2)} \end{align}\] which leads to \[\begin{align} \Omega_{m,I}(z)&=& \Omega_{m0}(1+z)^{3(1-d^2)} \end{align}\]

In this case, the Hubble parameter squared is given by: \[\begin{align} H_{1,I}^2(z)&=& H_0^2\left\{ \Omega_{m0}(1+z)^{3(1-d^2)}+\Omega_{k0}(1+z)^2+ \frac{3}{H_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{2\over n}\right\}\nonumber\\ &=& H_0^2\left\{ \Omega_{m0}(1+z)^{3(1-d^2)}+\Omega_{k0}(1+z)^2+ 3\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\,(1+z)^{2\over n}\right\} \label{unical3} \end{align}\tag{27}\] From the continuity equation given in Eq. (24 ), we obtain the following expression for the EoS parameter for DE: \[\begin{align} \omega_{D,1,I} = -1 - \frac{\dot{\rho}_{D,1}}{3H\rho_{D,1}}- \frac{Q}{3H\rho_{D,1}} \end{align}\] The term \(\frac{Q}{3H\rho_{D,1,I}}\) is given by: \[\begin{align} \frac{Q}{3H\rho_{D,1}} &=& d^2\left(\frac{\rho_{m,I}}{\rho_{D,1}}\right)\nonumber \\ &=&\frac{d^2\Omega_{m0}(1+z)^{3(1-d^2)}}{\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)(1+z)^{2\over n}} \end{align}\] where we used the fact that \(\rho_{m0} = 3H_0^2\Omega_{m0}\).
The final expression of the EoS parameter us then given by: \[\begin{align} \omega_{D,1,I}(z) = -1 + \frac{2}{3n}- \frac{d^2\Omega_{m0}(1+z)^{3(1-d^2)}}{\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)(1+z)^{2\over n}}\label{ruth} \end{align}\tag{28}\] In the limiting case of \(d^2=0\), we obtain the expression derived for the non-interacting case.
In Fig. (2) we plot the behavior of the EoS parameter of DE \(\omega_D\) obtained in Eq. (28 ) in the range of values of \(n \in \left[0.1-10\right]\). We considered \(d^2=0.02\). Moreover, we have chosen \(\Omega_{m0}=0.315\) and we have that \(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta\) is always 0.685 for all the combinations of values we have considered.

The present day value of the EoS is given by: \[\begin{align} \omega_{D,1,I,0} = -1 + \frac{2}{3n}- \frac{d^2\Omega_{m0}}{-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta } \end{align}\] Considering the values we are taking into account, we have that \[\begin{align} \frac{d^2\Omega_{m0}}{-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta }\approx 0.00920 \end{align}\] Therefore, the presence of interaction does not affect the final value of the EoS considerably.
We now want to calculate the expression of the pressure of DE for this case.
From the continuity equation for DE for the non interacting case, we obtain: \[\begin{align} p_{D,1,I}= -\rho_{D,1} - \frac{\dot{\rho}_{D,1}}{3H}-\frac{Q}{3H} \end{align}\] The term \(\frac{Q}{3H}\) is given by: \[\begin{align} \frac{Q}{3H} = d^2\rho_{m,I} = 3d^2H_0^2\Omega_{m0}(1+z)^{3(1-d^2)} \end{align}\]

Therefore, we obtain: \[\begin{align} p_{D,1,I}(z) = \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)H_0^2(1+z)^{2\over n}+3d^2H_0^2\Omega_{m0}(1+z)^{3(1-d^2)}. \end{align}\] In the limiting case of \(d^2=0\), we recover the same result of the non-interacting case.
We now want to obtain the final expression of the deceleration parameter \(q_{1,I}\).
We still use the general expression: \[\begin{align} q_{1,I}(z)= -1+\left(\frac{1+z}{2}\right)\frac{1}{h_{1,I}^2(z)}\frac{dh_{1,I}^2(z)}{dz} \end{align}\] Using the expression of the Hubble parameter squared obtained in Eq. (27 ), we can write: \[\begin{align} \frac{dh_{1,I}^2(z)}{dz} &=& 3(1-d^2) \Omega_{m0}(1+z)^{3(1-d^2)-1} +2\Omega_{k0}(1+z)\nonumber \\ &&+ \frac{6}{nH_0^2} \left( -\frac{6\alpha}{n} + 2\beta - n\gamma + \delta n^2 \right)b_0^{2\over n}\,(1+z)^{\frac{2}{n}-1} \end{align}\] Therefore, the final expression of \(q_{1,I}\) is given by: \[\begin{align} q_{1,I}(z) &=&-1+\frac{1}{2}\cdot \left[ 3(1-d^2)\Omega_{m0}(1+z)^{3(1-d^2)}+2\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \frac{6}{n}\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\,(1+z)^{\frac{2}{n}} \right]\times \nonumber \\ &&\left[ \Omega_{m0}(1+z)^{3(1-d^2)}+\Omega_{k0}(1+z)^2+ 3 \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\,(1+z)^{\frac{2}{n}} \right]^{-1} \end{align}\]

The present day value of \(q_{1,I}\) can be written as: \[\begin{align} q_{1,I,0} &=&-1+\frac{1}{2}\cdot \left[ 3(1-d^2)\Omega_{m0}+2\Omega_{k0}+ \frac{6}{n}\left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\, \right]\times \nonumber \\ &&\left[ \Omega_{m0}+\Omega_{k0}+ 3 \left(-\frac{6\alpha}{n^3} + \frac{2\beta}{n^2} - \frac{\gamma}{n}+ \delta \right)\, \right]^{-1} \end{align}\]

In Table 2, we have the values of \(q_{1,I,0}\) for the different values of \(n\) we considered and for the different values of the parameters we derived. Moreover, we considered \(\Omega_{m0}=0.315\) and \(\Omega_{k0}=0.01\).

4 Second Model↩︎

We now consider the second model of this paper.
Using the scale factor defined in Eq. (19 ), we obtain the following expression for the energy density of DE: \[\begin{align} \rho_{D,2}(z) &=&3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] H_0^{2}(1+z)^{2\over n} \end{align}\] In this case, the Hubble parameter squared can be written as: \[\begin{align} H_2^2(z)= H_0^2\left\{ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2+ 3 \left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right]\,(1+z)^{2\over n}\right\} \end{align}\] The expression of the EoS parameter of DE \(\omega_D\) is the same of the first model: \[\begin{align} \omega_{D,2} = -1+\frac{2}{3n} \end{align}\] Following the same procedure of the first model, we obtain that the pressure of DE is given by: \[\begin{align} p_{D,2}(z) &=& \left(\frac{2}{3n}-1\right) \rho_{D,2}(z) \nonumber\\ &=& \left( \frac{2}{n}-3\right) \left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] H_0^{2}(1+z)^{2\over n}. \end{align}\]

Finally, the deceleration parameter can be written as: \[\begin{align} q_2(z)&=& -1 + \frac{1}{2}\cdot \left\{ 3\Omega_{m0}(1+z)^3+2\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \frac{6}{n}\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \,(1+z)^{\frac{2}{n}} \right\}\times \nonumber \\ &&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left\{ \Omega_{m0}(1+z)^3+\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ 3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \,(1+z)^{\frac{2}{n}} \right\}^{-1} \end{align}\] Therefore, we have that the present-day value of \(q_2\) is given by: \[\begin{align} q_{2,0}&=& -1 + \frac{1}{2}\cdot \left\{ 3\Omega_{m0}+2\Omega_{k0}+ \frac{6}{n}\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \right\}\times \nonumber \\ &&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left\{ \Omega_{m0}+\Omega_{k0}+ 3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \right\}^{-1} \end{align}\] In Table 3, we write the values of \(q_{2,0}\) for the values of \(n\) we studied and for different combinations of values of the parameters we considered. We have also considered \(\Omega_{m0} =0.685\) and \(\Omega_{mk0} =0.001\)

We now want to find the evolutionary form of the fractional energy density of DE.
Using the expression of \(\rho_{D,2}\), we can write \(\Omega_{D,2}(z)\) as follow: \[\begin{align} \Omega_{D,2}(z) &=& \frac{\rho_{D,2}(z)}{3H_0^2}\nonumber\\ &=&\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{2\over n}. \end{align}\] Therefore, the evolutionary form of \(\Omega'_D\) is given by: \[\begin{align} \Omega'_{D,2}(z) &=&\frac{2}{n}\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{\frac{2}{n}-1}. \end{align}\] We also find that the present day value of \(\Omega_{D,2}\) is given by: \[\begin{align} \Omega_{D,2,0} &=&\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right]. \end{align}\] We can obtain some hints about the values of the parameters of the model using the value of \(\Omega_{D0}(z)\).
We now that at present time \(\Omega_{D,0}\approx 0.685\), therefore we should have: \[\begin{align} -\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \approx 0.685 \end{align}\] Therefore, once the value of \(n\) is fixed, we obtain a constraint on the values of \(\alpha\), \(\beta\), \(\gamma\), \(\delta\) and \(\zeta\).
For example, if we choose \(n=2\), we obtain the following set of possible values: \(\alpha = 2\), \(\zeta = 1\), \(\beta = 0.5\), \(\gamma = 0.2\) and \(\delta = 1.535\).
If we choose \(n=3\), we obtain this possible set of values: \(\alpha=1\), \(\beta=1\), \(\zeta=1\), \(\gamma=2\), \(\delta=1.13\).
If we choose \(n=4\), we obtain that \(\alpha=1\), \(\beta=0.13\), \(\gamma=1\), \(\delta=1\) e \(\zeta =0.1\).

4.1 Interacting Case↩︎

We now consider the interacting case.
The expression of \(\rho_D\) is the same obtained in the non-interacting case.
The Hubble parameter squared for this case can be written as: \[\begin{align} H_{2,I}^2(z)&=& H_0^2\left\{ \Omega_{m0}(1+z)^{3(1-d^2)}+\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,+ 3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \,(1+z)^{2\over n}\right\} \end{align}\]

From the continuity equation given in Eq. (24 ), we obtain the following expression for the EoS parameter for DE: \[\begin{align} \omega_{D,2,I} = -1 - \frac{\dot{\rho}_{D,2}}{3H\rho_{D,2}}- \frac{Q}{3H\rho_{D,2}} \end{align}\] The term \(\frac{Q}{3H\rho_{D,2}}\) is given by: \[\begin{align} \frac{Q}{3H\rho_{D,2}} &=& d^2\left(\frac{\rho_{m,I}}{\rho_{D,2}}\right)\nonumber \\ &=&\frac{d^2\Omega_{m0}(1+z)^{3(1-d^2)}}{\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{2\over n}} \end{align}\] where we used the fact that \(\rho_{m0} = 3H_0^2\Omega_{m0}\).
The final expression of the EoS parameter is then given by: \[\begin{align} \omega_{D,2,I} = -1 +\frac{2}{3n}- \frac{d^2\Omega_{m0}(1+z)^{3(1-d^2)}}{\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{2\over n}} \end{align}\] In the limiting case of \(d^2=0\), we obtain the expression derived for the non-interacting case.
The behavior of \(\omega_{D,2,I}\) is the same as \(\omega_{D,1,I}\) since the values of the parameters involved in the interacting term assume similar values.
The present-day value of the EoS is given by: \[\begin{align} \omega_{D,2,I,0} = -1 + \frac{2}{3n}- \frac{d^2\Omega_{m0}}{-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta } \end{align}\] As in the other model, the interaction term does not affect the value of the EoS considerably.

We now want to calculate the expression of the pressure of DE for this case.
From the continuity equation for DE for the non-interacting case, we obtain: \[\begin{align} p_{D,2,I}= -\rho_{D,2} - \frac{\dot{\rho}_{D,2}}{3H}-\frac{Q}{3H} \end{align}\] The term \(\frac{Q}{3H}\) is given by: \[\begin{align} \frac{Q}{3H} = d^2\rho_{m,I} = 3d^2H_0^2\Omega_{m0}(1+z)^{3(1-d^2)} \end{align}\]

Therefore, we obtain: \[\begin{align} p_{D,2,I}(z) = \left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right]H_0^2(1+z)^{2\over n}+3d^2H_0^2\Omega_{m0}(1+z)^{3(1-d^2)}. \end{align}\] In the limiting case of \(d^2=0\), we recover the same result of the non-interacting case.
For the deceleration parameter \(q\), we use the general definition introduced before: \[\begin{align} q_{2,I} = -1+\left(\frac{1+z}{2}\right)\frac{1}{h_{2,I}^2(z)}\frac{dh_{2,I}^2(z)}{dz} \end{align}\] Using the expression of \(H_{2,I}^2(z)\), we obtain: \[\begin{align} \frac{dh_{2,I}^2(z)}{dz} &=& 3(1-d^2) \Omega_{m0}(1+z)^{3(1-d^2)-1} +2\Omega_{k0}(1+z)\nonumber \\ &&+ \frac{6}{n} \left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{\frac{2}{n}-1} \end{align}\] where \(h_{2,I}^2= H_{2,I}^2/H_0^2\).
Therefore, the final expression of \(q_{2,I}\) is given by: \[\begin{align} q_{2,I}&=& -1 + \frac{1}{2}\cdot \left\{ 3(1-d^2)\Omega_{m0}(1+z)^{3(1-d^2)}+2\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \frac{6}{n}\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right](1+z)^{\frac{2}{n}} \right\}\times \nonumber \\ &&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left\{ \Omega_{m0}(1+z)^{3(1-d^2)}+\Omega_{k0}(1+z)^2\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ 3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right]\,(1+z)^{\frac{2}{n}} \right\}^{-1} \end{align}\] The present-day value of \(q_{2,I}\) is given by: \[\begin{align} q_{2,I,0}&=& -1 + \frac{1}{2}\cdot \left\{ 3(1-d^2)\Omega_{m0}+2\Omega_{k0}\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ \frac{6}{n}\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \right\}\times \nonumber \\ &&\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\left\{ \Omega_{m0}+\Omega_{k0}\right.\nonumber \\ &&\left.\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,+ 3\left[-\frac{6\alpha}{n^3} + \frac{2(\zeta+\beta)}{n^2} - \frac{\gamma}{n}+ \delta \right] \right\}^{-1} \end{align}\] In Table 4, we write the values of \(q_{2,I,0}\) for the different values of \(n\) we choose and the combinations of the values of the parameters we considered. Moreover, we considered \(\Omega_{m0} =0.685\), \(\Omega_{mk0} =0.001\) and \(d^2=0.02\).

5 Age of the Present Universe↩︎

The determination of the Universe’s age is one of the cornerstone results of modern cosmology. Within the standard \(\Lambda\)CDM paradigm, this age can be estimated by integrating the Friedmann equations backward in time from the present epoch to the initial singularity, using the observed Hubble expansion rate. The calculated value is highly sensitive to the cosmological parameters, particularly the current Hubble constant \(H_0\), the matter density parameter \(\Omega_m\), and the dark energy contribution \(\Omega_D\). Current measurements, for instance those from the Planck mission, suggest an age of \(t_0 \simeq 13.8 \, \mathrm{Gyr}\), with uncertainties at the percent level. Nevertheless, the ongoing \(H_0\) tension between early- and late-time observations implies that slightly different values of \(H_0\) could shift the inferred cosmic age by several hundred million years. Moreover, alternative cosmological frameworks—such as interacting dark energy models, modifications of general relativity, or scenarios with non-zero spatial curvature—can also affect the theoretical estimate of \(t_0\), providing an additional observational handle to discriminate between competing models of the Universe.
The age of universe can be determined thanks to: \[\begin{align} t_0 -t = -\int_{t_0}^t dt = \int_0^z\frac{dz'}{(1+z')H(z')} \label{ageuni} \end{align}\tag{29}\] We now calculate the present age of the Universe for the models we studied in this paper.
Eq. (29 ) does not have an analytical solution for these models, therefore we solve it numerically considering a value of \(z=100\).
For the first model, in the case with \(n=2\) and \(\alpha = 1\), \(\beta = 2\), \(\gamma = 0.5\) and \(\delta = 0.685\) and for the case with \(n=2\) and \(\alpha = 0.5\), \(\beta = 1\), \(\gamma = 0.2\) and \(\delta = 0.66\), we obtain \(t_0-t\approx 9.992\) Gyr for the non interacting case while \(t_0-t\approx 10.163\) Gyr for the interacting case.
Instead, for \(n=3\) and \(\alpha = 1\), \(\beta = 2\), \(\gamma = 0.5\) e \(\delta = 0.63\), we obtain \(t_0-t \approx 10.452\) Gyr for the non interacting case, while we obtain \(t_0-t \approx 10.644\) Gyr for the interacting case.
For \(n=4\) and \((\alpha,\beta,\gamma,\delta)=(1,1,1,0.90375)\), we obtain \(t_0-t \approx 10.666\) Gyr for the non interacting case and \(t_0-t \approx 10.868\) Gyr for the interacting case.
We now consider the second DE energy density model we studied.
For \(n=2\) and \(\alpha = 2\), \(\zeta = 1\), \(\beta = 0.5\), \(\gamma = 0.2\) and \(\delta = 1.535\), we obtain \(t_0-t \approx 9.992\) Gyr for the non-interacting case and \(t_0-t \approx 10.163\) Gyr for the interacting case.
When \(n=3\) and \(\alpha=1\), \(\beta=1\), \(\zeta=1\), \(\gamma=2\), \(\delta=1.13\), we obtain \(t_0-t \approx 10.452\) Gyr for the non-interacting case and \(t_0-t \approx 10.644\) Gyr for the interacting case.
For \(n=4\) and \(\alpha=1\), \(\beta=0.13\), \(\gamma=1\), \(\delta=1\) e \(\zeta =0.1\), we obtain \(t_0-t \approx 10.666\) Gyr for the non interacting case and \(t_0-t \approx 10.868\) Gyr for the interacting case.

6 Conclusions↩︎

In this paper, we studied two dark energy (DE) models involving higher-order derivatives of the Hubble parameter \(H\). These models can be considered as generalizations of previously studied DE scenarios.

By assuming a power-law form for the scale factor \(a(t)\), we derived expressions for the DE energy density, DE pressure, equation-of-state (EoS) parameter, deceleration parameter, and the evolutionary behavior of the fractional DE energy density for both non-interacting and interacting dark sector scenarios. For the interacting case, we adopted a coupling term proportional to the dark matter (DM) energy density.

We determined several combinations of parameter values for the models corresponding to power-law indices \(n=2\), \(n=3\), and \(n=4\). Using these parameter sets, we calculated the present-day age of the Universe. The resulting values differ from the observationally inferred age of \(\approx 13.8\) Gyr.

Nevertheless, we observe that the calculated age increases with \(n\), and that interacting scenarios consistently yield higher values compared to the non-interacting cases.

Future works can be devoted to studying the search for best fit values for the involved parameters. Furthermore, these models can be studied considering other scale factors in order to better understand their behaviour.

References↩︎