# Observational constraints and dynamical analysis of Kaniadakis horizon-
entropy cosmology

A. Hernández-Almada ID 1 , Genly Leon ID 2,3, Juan Magaña ID 4, Miguel A.
García-Aspeitia ID 5, V. Motta ID 6, Emmanuel N. Saridakis ID 7,8,9, Kuralay
Yesmakhanova ID 9,10, Alfredo D. Millano ID 2
1Facultad de Ingeniería, Universidad Autónoma de Querétaro, Centro
Universitario Cerro de las Campanas, 76010, Santiago de Querétaro, México,
2Departamento de Matemáticas, Universidad Católica del Norte, Avda. Angamos
0610, Casilla 1280 Antofagasta, Chile,
3Institute of Systems Science, Durban University of Technology, PO Box 1334,
Durban 4000, South Africa
4Instituto de Astrofísica & Centro de Astro-Ingeniería, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna, 4860, Santiago, Chile,
5 Depto. de Física y Matemáticas, Universidad Iberoamericana Ciudad de México,
Prolongación Paseo
de la Reforma 880, México D. F. 01219, México
6Instituto de Física y Astronomía, Facultad de Ciencias, Universidad de
Valparaíso,
Avda. Gran Bretaña 1111, Valparaíso, Chile,
7Institute for Astronomy, Astrophysics, Space Applications and Remote Sensing,
National Observatory of Athens, Lofos Nymfon, 11852 Athens, Greece,
8CAS Key Laboratory for Researches in Galaxies and Cosmology, Department of
Astronomy,
University of Science and Technology of China, Hefei, Anhui 230026, P.R.
China,
9Ratbay Myrzakulov Eurasian International Centre for Theoretical Physics, Nur-
Sultan 010009, Kazakhstan,
10Ratbay Myrzakulov Eurasian International Centre for Theoretical Physics,
Eurasian National University, Nur-Sultan Astana 010008, Kazakhstan Contact
e-mail<EMAIL_ADDRESS>
###### Abstract

We study the scenario of Kanadiakis horizon entropy cosmology which arises
from the application of the gravity-thermodynamics conjecture using the
Kaniadakis modified entropy. The resulting modified Friedmann equations
contain extra terms that constitute an effective dark energy sector. We use
data from Cosmic chronometers, Supernova Type Ia, HII galaxies, Strong lensing
systems, and Baryon acoustic oscillations observations and we apply a Bayesian
Markov Chain Monte Carlo analysis to construct the likelihood contours for the
model parameters. We find that the Kaniadakis parameter is constrained around
0, namely, around the value where the standard Bekenstein-Hawking is
recovered. Concerning the normalized Hubble parameter, we find
$h=0.708^{+0.012}_{-0.011}$, a result that is independently verified by
applying the $\mathbf{\mathbb{H}}0(z)$ diagnostic and, thus, we conclude that
the scenario at hand can alleviate the $H_{0}$ tension problem. Regarding the
transition redshift, the reconstruction of the cosmographic parameters gives
$z_{T}=0.715^{+0.042}_{-0.041}$. Furthermore, we apply the AICc, BIC and DIC
information criteria and we find that in most datasets the scenario is
statistical equivalent to $\Lambda$CDM one. Moreover, we examine the Big Bang
Nucleosynthesis (BBN) and we show that the scenario satisfies the
corresponding requirements. Additionally, we perform a phase-space analysis,
and we show that the Universe past attractor is the matter-dominated epoch,
while at late times the Universe results in the dark-energy-dominated
solution. Finally, we show that Kanadiakis horizon entropy cosmology accepts
heteroclinic sequences, but it cannot exhibit bounce and turnaround solutions.

###### keywords:

theory, dark energy, cosmological parameters.

††pubyear: 2022††pagerange: Observational constraints and dynamical analysis
of Kaniadakis horizon-entropy cosmology–References

##  1 Introduction

Cosmology is one of the most exciting adventures in the human endeavour: the
understanding of the origin, evolution and future of our Universe by combining
the physics at micro- and macro-scales in a joint framework. In particular,
the Universe acceleration is one of the most important puzzles. Discovered by
the Supernovae team, through type Ia supernovae (SNIa, Riess et al., 1998), it
is also confirmed by the acoustic peaks of the Cosmic Microwave Background
Radiation (CMB, Spergel et al., 2003), and recently tested with large-scale
structure measurements (Nadathur et al., 2020). Tackling the Universe
acceleration is indeed a complicated issue due to the attractive nature of
gravity in the General Relativity (GR) framework. The first approach on its
explanation is through the addition of the cosmological constant (CC, Carroll,
2001). Thus, using the continuity equation and assuming a constant energy
density ($\dot{\rho}=0$), it is possible to conclude that the equation of
state (EoS) is $w=-1$, which is in concordance with what is expected for a
fluid that accelerates the Universe. Another main characteristic is that the
energy density of the CC must be subdominant in order to obtain a late and
non-violent acceleration. From the Quantum Field Theory (QFT) viewpoint, the
CC can be explained by the addition of quantum vacuum fluctuations (QVF)
associated to the space-time. Therefore, the expansion of space-time implies
an increase of QVF, maintaining a constant energy density. However, when we
calculate the energy density of the QVF, the result is in complete
disagreement with the observed one (see Zel’dovich et al., 1968; Weinberg,
1989), this is the so-called fine tuning problem (Addazi et al., 2021). In
addition to the cosmological constant problem, a $4.2\sigma$ tension in the
current Hubble parameter value ($H_{0}$) measured by Supernova $H_{0}$ for the
Equation of State (SH0ES) collaboration (Riess et al., 2021) and the one
obtained by Planck collaboration under the $\Lambda$-Cold Dark Matter
($\Lambda$CDM) scenario (Aghanim et al., 2020) has been recently observed.

The above mentioned problems that afflicts the understanding of the CC in this
framework has driven the community to propose other approaches like scalar
fields, dynamical dark energy, viscous fluids, Chaplygin gas (Hernández-
Almada, A. et al., 2019) or modifications to GR such as unimodular gravity
(García-Aspeitia et al., 2019, 2021), Einstein-Gauss-Bonet (García-Aspeitia &
Hernández-Almada, 2021), Brane Worlds (Garcia-Aspeitia et al., 2018), among
others (see Saridakis et al., 2021; Motta et al., 2021, for a compilation of
the mentioned previous models). Despite the fact that we can have numerous
models and scenarios describing the late time acceleration, at the end of the
day the detailed confrontation with observations, alongside theoretical
consistency, will be the main method for their validation.

One interesting alternative to investigate the dynamics of the Universe is the
gravity-thermodynamics approach (Jacobson, 1995; Padmanabhan, 2005, 2010).
This takes advantage of the first law of thermodynamics and the standard
Bekenstein-Hawking entropy-area relation for black holes, which applied to the
apparent cosmological horizon leads to the Friedmann equations (Frolov &
Kofman, 2003; Cai & Kim, 2005; Akbar & Cai, 2007; Cai & Cao, 2007). Hence, the
application of this conjecture with various alternative entropy relations
leads to modified Friedmann equations whose additional terms can source the
cosmic acceleration (Cai & Cao, 2007; Akbar & Cai, 2006; Paranjape et al.,
2006; Sheykhi et al., 2007; Jamil et al., 2010b; Cai & Ohta, 2010; Wang et
al., 2010; Jamil et al., 2010a; Sheykhi, 2010b, a; Gim et al., 2014; Fan & Lu,
2015; Lymperis & Saridakis, 2018; Sheykhi, 2018; Saridakis, 2020).

One interesting class of extended entropies arises through generalizations for
the Boltzmann-Gibbs statistics. In particular, one modifies the classical
entropy of a system given by $S=-k_{B}\sum_{i}p_{i}\ln{p_{i}}$, where $p_{i}$
is the probability of a system to be within a microstate, through non-
extensive analysis resulting into the Tsallis entropy (Tsallis, 1988; Lyra &
Tsallis, 1998), through quantum-gravitational considerations resulting into
Barrow entropy (Barrow, 2020), or through relativistic extensions resulting
into Kaniadakis entropy (Kaniadakis, 2002, 2005). Hence, the application of
the gravity-thermodynamics conjecture using the above modified horizon entropy
gives rise to modified cosmological scenarios. In Lymperis & Saridakis (2018)
this was performed in the framework of Tsallis entropy, in Saridakis (2020);
Saridakis & Basilakos (2021); Leon et al. (2021); Barrow et al. (2021a)
applied it for the Barrow entropy case, and recently (Lymperis et al., 2021)
use it in the Kaniadakis entropy frame.

In the present work, we investigate Kaniadakis horizon entropy cosmology by
confronting it to observational data from Cosmic Chronometers, Supernovae Type
I (SNIa), HII galaxies, strong lensing systems (SLS), and baryon acoustic
oscillations (BAO) observations. Additionally, we perform a complete dynamical
system analysis in order to etxract information of the local and global
features of the cosmological evolution.

The outline of the paper is as follows: Sec. 2 introduces the framework of the
Kaniadakis cosmology. In Sec. 3 we infer the cosmological parameter under the
Kaniadakis model using five observational datasets mentioned above. The
section 4 presents a stability analysis around the equilibrium points of the
dynamical system under the Kaniadakis cosmology. Finally, we discuss and
summarize our results in Section 5. From now on we use natural units in which
$\hbar=k_{B}=c=1$.

##  2 Kaniadakis horizon entropy cosmology

In this section we present the scenario of Kaniadakis horizon entropy
cosmology (Lymperis et al., 2021), namely the modified Friedmann equations
arising from the application of the gravity-thermodynamics conjecture using
the extended Kaniadakis entropy.

Kaniadakis entropy is an one-parameter generalization of the classical
entropy, given as $S_{K}=-k_{{}_{B}}\sum_{i}n_{i}\,\ln_{{}_{\\{{\scriptstyle
K}\\}}}\\!n_{i}$ (Kaniadakis, 2002, 2005), with $k_{{}_{B}}$ the Boltzmann
constant and where $\ln_{{}_{\\{{\scriptstyle K}\\}}}\\!x=(x^{K}-x^{-K})/2K$.
In such a framework the dimensionless parameter $-1<K<1$ quantifies the
relativistic deviations from standard statistical mechanics, and the latter is
recovered in the limit $K\rightarrow 0$. Kaniadakis entropy can be expressed
as (Abreu et al., 2016, 2018; Abreu & Ananias Neto, 2021)

$S_{K}=-k_{{}_{B}}\sum^{W}_{i=1}\frac{P^{1+K}_{i}-P^{1-K}_{i}}{2K},$ (1)

with $P_{i}$ the probability of a specific microstate and $W$ the total
configuration number. When we apply it in the black-hole framework, we obtain
(Drepanou et al., 2021; Moradpour et al., 2020; Lymperis et al., 2021)

$S_{K}=\frac{1}{K}\sinh{(KS_{BH})},$ (2)

where

$S_{BH}=\frac{1}{4G}A,$ (3)

is the usual Bekenstein-Hawking entropy, with $A$ the horizon area and $G$ is
the gravitational constant. Hence, in the limit $K\rightarrow 0$ we recover
Bekenstein-Hawking entropy.

Let is now apply the gravity-thermodynamics conjecture using Kaniadakis
entropy. We consider a flat homogeneous and isotropic Friedmann-Robertson-
Walker (FRW) metric of the form

$ds^{2}=-dt^{2}+a^{2}(t)\left(dr^{2}+r^{2}d\Omega^{2}\right),$ (4)

where $d\Omega^{2}\equiv d\theta^{2}+\sin^{2}\theta d\varphi^{2}$ is the solid
angle, and $a(t)$ the scale factor. In this setup, the first law of
thermodynamics is interpreted in terms of the heat/energy that flows through
the apparent horizon of the Universe (Jacobson, 1995; Padmanabhan, 2005,
2010), which in the case of flat geometry is just (Bak & Rey, 2000; Frolov &
Kofman, 2003; Cai & Kim, 2005; Cai et al., 2009):

$r_{a}=\frac{1}{H},$ (5)

where $H=\dot{a}/a$ is the Hubble parameter (dots denote derivatives with
respect to $t$). Concerning the horizon temperature, this is given by the
standard relation (Gibbons & Hawking, 1977):

$T=\frac{1}{2\pi r_{a}}.$ (6)

Hence, the first law of thermodynamics is just $-dE=TdS$, where
$-dE=A(\rho_{m}+p_{m})Hr_{a}dt$ is the energy flow through the horizon during
a time interval $dt$ in the case of a Universe filled with a matter perfect
fluid with energy density $\rho_{m}$ and pressure $p_{m}$ (Cai & Kim, 2005).
Differentiating (2) we obtain

$dS_{K}=\frac{8\pi}{4G}\cosh{\left(K\frac{\pi}{GH^{2}}\right)}r_{a}\dot{r}_{a}dt,$
(7)

where we have used that $A=4\pi r_{a}^{2}=4\pi/H^{2}$.

Inserting (3), (6), and (7) into the first law of thermodynamics, alongside
with the relation $\dot{r}_{a}=-\dot{H}/H^{2}$, we acquire (Lymperis et al.,
2021)

$-4\pi G(\rho_{m}+p_{m})=\dot{H}\cosh{\left[K\frac{\pi}{GH^{2}}\right]}.$ (8)

Thus, using the matter conservation equation

$\dot{\rho}_{m}+3H(\rho_{m}+p_{m})=0,$ (9)

the integration of (8) leads to

$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\frac{8\pi
G}{3}\rho_{m}=H^{2}\cosh{\left[K\frac{\pi}{GH^{2}}\right]}-\frac{K\pi}{G}\text{shi}{\left[K\frac{\pi}{GH^{2}}\right]}-\frac{\Lambda}{3},$
(10)

with $\Lambda$ the integration constant and ${\rm
shi}(x)\equiv\int_{0}^{x}\sinh(x^{\prime})/x^{\prime}dx^{\prime}$, a
mathematical odd function of $x$ with no discontinuities.

The above equations (8) and (10) are the modified Friedmann equations in the
scenario of Kaniadakis horizon entropy cosmology (Lymperis et al., 2021). As
expected, in the limit $K\rightarrow 0$ they turn into the standard ones. We
can re-write them as

$\displaystyle H^{2}$ $\displaystyle=\frac{8\pi G}{3}(\rho_{m}+\rho_{DE}),$
(11) $\displaystyle\dot{H}$ $\displaystyle=-4\pi
G(\rho_{m}+\rho_{DE}+p_{m}+p_{DE}),$ (12)

where we have introduced an effective dark-energy sector, with energy density
and pressure respectively of the form

$\displaystyle\rho_{DE}=\frac{3}{8\pi
G}\Bigg{\\{}\frac{\Lambda}{3}+H^{2}\left[1-\cosh\left(\frac{\pi
K}{GH^{2}}\right)\right]+\frac{\pi K}{G}{\rm shi}\left(\frac{\pi
K}{GH^{2}}\right)\Bigg{\\}},$ (13) $\displaystyle p_{DE}=-\frac{1}{8\pi
G}\Bigg{\\{}\Lambda+(3H^{2}+2\dot{H})\left[1-\cosh\left(\frac{\pi
K}{GH^{2}}\right)\right]+\frac{3\pi K}{G}{\rm shi}\left(\frac{\pi
K}{GH^{2}}\right)\Bigg{\\}}.$ (14)

Therefore, the EoS parameter for the dark energy sector is

$\displaystyle w_{DE}=$
$\displaystyle-1-2\dot{H}\left[1\\!-\\!\cosh\left(\frac{\pi
K}{GH^{2}}\right)\\!\right]\times$
$\displaystyle\left\\{\Lambda+3H^{2}\left[1\\!-\\!\cosh\left(\frac{\pi
K}{GH^{2}}\right)\\!+\\!\frac{3\pi K}{G}{\rm shi}\left(\frac{\pi
K}{GH^{2}}\right)\right]\right\\}^{-1}.$ (15)

In the general case, the cosmological equations are the two Friedmann
equations (11) and (12), alongside the matter conservation equation (9). For
convenience we focus on on the dust case, namely we consider $p_{m}=0$. It
proves convenient to express the equations in terms of dimensionless
variables. Introducing the density parameters

$\displaystyle\Omega_{\Lambda}\equiv\frac{\Lambda}{3H^{2}},\quad\Omega_{m}\equiv\frac{8\pi
G\rho_{m}}{3H^{2}},$ (16)

the normalized Hubble function

$\displaystyle E\equiv\frac{H}{H_{0}},$ (17)

with $H_{0}$ the Hubble parameter at the present scale factor $a_{0}$, and
defining the dimensionless parameter $\beta$ as

$\displaystyle\beta\equiv\frac{K\pi}{GH_{0}^{2}},$ (18)

then, the cosmological equations are expressed as

$\displaystyle\Omega_{\Lambda}^{\prime}(N)=3\Omega_{\Lambda}\Omega_{m}\text{sech}\left(\frac{\beta}{E^{2}}\right),$
(19)
$\displaystyle\Omega_{m}^{\prime}(N)=3\Omega_{m}\left[\Omega_{m}\text{sech}\left(\frac{\beta}{E^{2}}\right)-1\right],$
(20) $\displaystyle
E^{\prime}(N)=-\frac{3}{2}E\Omega_{m}\text{sech}\left(\frac{\beta}{E^{2}}\right),$
(21)

where primes denote derivatives with respect to the e-foldings number
$N=\ln(a/a_{0})$ (and thus $f^{\prime}=\dot{f}/H$). Note that using the above
variables, the first Friedmann equation (11) gives rise to the constraint

$\displaystyle\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)+E^{2}\left[-\cosh\left(\frac{\beta}{E^{2}}\right)+\Omega_{\Lambda}+\Omega_{m}\right]=0,$
(22)

which allows us to eliminate $\Omega_{\Lambda}$ in terms of $\Omega_{m}$ and
$E$. Finally, note that for the effective dark energy density parameter, in
the general case we have

$\Omega_{DE}=1-\Omega_{m}=1+\beta\text{shi}\left(\beta\right)-\cosh\left(\beta\right)+\Omega_{\Lambda}.$
(23)

Lastly, it proves convenient to introduce the deceleration parameter $q(z)$,
and the cosmographic jerk parameter $j(z)$, which are defined as

$\displaystyle q:=$ $\displaystyle-1-\frac{E^{\prime}}{E},$ (24)
$\displaystyle j:=$ $\displaystyle q(2q+1)-q^{\prime},$ (25)

where $j=1$ recovers the case of a cosmological constant.

In the scenario of Kaniadakis horizon entropy cosmology one may have the
general integration constant $\Lambda$, which will play the role of an
explicit cosmological constant, or one may set it to zero and thus require for
the extra $K$-dependent terms to drive the Universe acceleration. Since the
corresponding equation structure (which will be used later for the the
Bayesian statistical analysis and the dynamical system approach) is different
in the two cases, we examine them separately in the following subsections.

###  2.1 Case I: $\Lambda\neq 0$

In the general case where $\Lambda\neq 0$, i.e. $\Omega_{\Lambda}\neq 0$, we
use the constraint equation (22) to obtain the reduced dynamical system

$\displaystyle\Omega_{m}^{\prime}(N)=3\Omega_{m}\left[\Omega_{m}\text{sech}\left(\frac{\beta}{E^{2}}\right)-1\right],$
(26) $\displaystyle
E^{\prime}(N)=-\frac{3}{2}E\Omega_{m}\text{sech}\left(\frac{\beta}{E^{2}}\right).$
(27)

This is integrable, with

$\displaystyle\Omega_{m}(E)=\cosh\left(\frac{\beta}{E^{2}}\right)+\frac{-\cosh(\beta)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)+\beta\text{shi}(\beta)+\Omega_{m}^{(0)}}{E^{2}},$
(28) $\displaystyle\Omega_{m}(1)=\Omega_{m}^{(0)},$ (29)

and

$\displaystyle
E^{\prime}(N)=-\frac{3}{2}E-\frac{3\text{sech}\left(\frac{\beta}{E^{2}}\right)\left(-\cosh(\beta)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)+\beta\text{shi}(\beta)+\Omega_{m}^{(0)}\right)}{2E},$
(30)

where $\Omega_{m}(N=0)=\Omega_{m}^{(0)}$ and $E(N=0)=1$.

The equation (30) is easily integrated to give

$\displaystyle{3}N(E)=$
$\displaystyle-\ln\Bigg{(}1-\frac{\cosh(\beta)-E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)+\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}(\beta)}{\Omega_{m}^{(0)}}\Bigg{)},$
(31)

introducing as dynamical variable the redshift $e^{N}=a=(1+z)^{-1}$, where
$z=0$ and $a_{0}=1$ for current time, the previous equation leads to

$\displaystyle(z+1)^{{3}}=1-\frac{\left[\cosh(\beta)-E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)+\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}(\beta)\right]}{\Omega_{m}^{(0)}}.$
(32)

Evaluating (23) at present time gives

$\Omega_{DE}^{(0)}=1-\Omega_{m}^{(0)}=1+\beta\text{shi}\left(\beta\right)-\cosh\left(\beta\right)+\Omega_{\Lambda}^{(0)},$
(33)

and combining it with (32) we have

$\displaystyle\Omega_{m}^{(0)}(z+1)^{3}+\Omega_{\Lambda}^{(0)}=E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right).$
(34)

We mention here that in the general case where $\Lambda\neq 0$ from the above
we obtain the relation between $\beta$, $\Omega_{m}^{(0)}$ and
$\Omega_{\Lambda}^{(0)}$ as

$\Omega_{m}^{(0)}=\cosh\left(\beta\right)-\beta\text{shi}\left(\beta\right)-\Omega_{\Lambda}^{(0)}.$
(35)

We expand (33), (34) up to third order around $\beta=0$, resulting in

$\Omega_{m}^{(0)}+\Omega_{\Lambda}^{(0)}=1-\frac{\beta^{2}}{2}+O\left(\beta^{4}\right),$
(36)

and

$\displaystyle\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}$
$\displaystyle\approx E^{2}-\frac{\beta^{2}}{2E^{2}}+O\left(\beta^{4}\right).$
(37)

Hence, we obtain four roots:

$\displaystyle
E_{1,2}=\mp\frac{\sqrt{\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}-\sqrt{2\beta^{2}+\left[\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}\right]^{2}}}}{\sqrt{2}},$
(38) $\displaystyle
E_{3,4}=\mp\frac{\sqrt{\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}+\sqrt{2\beta^{2}+\left[\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}\right]^{2}}}}{\sqrt{2}}.$
(39)

Solutions $E_{1}$ and $E_{2}$ are complex, and $E_{3}$ is negative. Thus, the
only physical solution is $E_{4}$. In the following, instead of using the
exact implicit formula for $E$ given in (34), we will consider the
approximation $E_{4}$ in (39).

###  2.2 Case II: $\Lambda=0$

In the case where an explicit cosmological constant is absent, namely
$\Lambda=0$, i.e $\Omega_{\Lambda}=0$, the general system (26), (27) reduces
to

$E^{\prime}(N)=-\frac{3}{2}E+\frac{3\beta\text{sech}\left(\frac{\beta}{E^{2}}\right)\text{shi}\left(\frac{\beta}{E^{2}}\right)}{2E}.$
(40)

The last equation is easily integrated to give

$\displaystyle
N(E)=-\frac{1}{3}\ln\left[\frac{E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)}{\cosh(\beta)-\beta\text{shi}(\beta)}\right],$
(41)

which, using the redshift, implies

$\displaystyle(z+1)^{{3}}=\frac{E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)}{\cosh(\beta)-\beta\text{shi}(\beta)}.$
(42)

Hence, using

$\Omega_{m}^{(0)}={\cosh\left(\beta\right)-\beta\text{shi}\left(\beta\right)},$
(43)

as it arises from (35) for $\Omega_{\Lambda}^{(0)}=0$, we obtain

$\displaystyle\Omega_{m}^{(0)}(z+1)^{{3}}={E^{2}\cosh\left(\frac{\beta}{E^{2}}\right)-\beta\text{shi}\left(\frac{\beta}{E^{2}}\right)}.$
(44)

Expanding (44) up to third order around $\beta=0$, we result to

$\displaystyle\Omega_{m}^{(0)}(z+1)^{{3}}$ $\displaystyle\approx
E^{2}-\frac{\beta^{2}}{2E^{2}}+{\mathcal{O}}\left(\beta^{4}\right),$ (45)

and thus at present times gives
$\Omega_{m}^{(0)}=1-\frac{\beta^{2}}{2}+{\mathcal{O}}\left(\beta^{4}\right)$.
Therefore, we obtain four roots:

$\displaystyle
E_{1,2}=\mp\frac{\sqrt{{\Omega_{m}^{(0)}}(z+1)^{3}-\sqrt{2\beta^{2}+{\Omega_{m}^{(0)}}^{2}(z+1)^{6}}}}{\sqrt{2}},$
(46) $\displaystyle
E_{3,4}=\mp\frac{\sqrt{{\Omega_{m}^{(0)}}(z+1)^{3}+\sqrt{2\beta^{2}+{\Omega_{m}^{(0)}}^{2}(z+1)^{6}}}}{\sqrt{2}}.$
(47)

Similarly to the previous case, roots $E_{1}$ and $E_{2}$ are complex while
$E_{3}$ is negative. Consequently, the only physical solution is $E_{4}$ in
(47).

##  3 Observational constraints

In this section we confront the scenario of Kaniadakis horizon-entropy
cosmology with observations. We are interested in extracting the bounds on the
parameter phase-space $\Theta=\\{h,\Omega_{m}^{(0)},\beta\\}$ and
$\\{h,\beta\\}$, particularly on the parameter $\beta$, which is related to
the Kaniadakis basic parameter $K$. For convenience, we focus on the
physically interested case of dust matter, namely we set $w_{m}=0$.

###  3.1 Datasets and methodology

We will employ the most commonly used datasets.

  * •

Observational Hubble Data (OHD). The sample contains $31$ cosmological-
independent measurements of the Hubble parameter in the redshift range
$0.07<z<1.965$ from passive elliptic galaxies, the so-called cosmic
chronometers (Moresco et al., 2016).

  * •

Pantheon Supernova Type Ia sample (SNIa). We use 1048 data points of the
distance modulus, $\mu(z)_{\mathrm{SNIa}}$, of high-redshift SNIa in the
redshift range $0.001<z<2.3$ (Scolnic et al., 2018).

  * •

HII galaxies (HIIG). It contains a total of $181$ data points of the distance
modulus $\mu_{\mathrm{HIIG}}(z)$ estimated from the Balmer line luminosity-
velocity dispersion relation for HII galaxies spanning the redshift region
$0.01<z<2.6$ (González-Morán et al., 2021).

  * •

Strong lensing systems (SLS). We use the sample by Amante et al. (2020) which
contains $143$ strong lensing systems by elliptical galaxies with measurements
of the redshift for the lens and the source, spectroscopic velocity dispersion
and the Einstein radius. These quantities allow us to construct an
observational distance ratio within the region $0.5\leq D^{obs}\leq 1$ .

  * •

Baryon acoustic oscillations (BAO). We consider $6$ correlated data points of
the imprint of baryon acoustic oscillation in the size of the sound horizon in
clustering and power spectrum of galaxies measured by Percival et al. (2010);
Blake et al. (2011); Beutler et al. (2011) and collected by Giostri et al.
(2012).

We would like to mention here that other cosmological observations could be
included in the parameter estimation too, for instance the CMB data. To
perform such analysis in a robust way, a full perturbation approach is needed
in order to obtain the linear Einstein-Boltzmann equations. Nevertheless this
is beyond the scope of the present work. An alternative approach would be to
use the distance priors from Planck 2018 based on slight deviations from
$\Lambda$CDM, such as the $w$CDM model (Chen et al., 2019). However, since
this procedure could lead to biased constraints, in the following we prefer
not to use the CMB dataset.

The inference of the cosmological parameters under Kaniadakis horizon entropy
cosmology for both scenarios ($\Lambda\neq 0$ and $\Lambda=0$) is performed by
a Bayesian Markov Chain Monte Carlo (MCMC) approach and the emcee Python
module (Foreman-Mackey et al., 2013). We set $3000$ chains with $250$ steps
each, and consider uniform priors in the ranges: $h:[0.2,1]$,
$\Omega_{m}^{(0)}:[0,1]$, $\beta:[-\pi,\pi]$. The burn-in phase is stopped up
to obtain convergence according to the auto-correlation time criterion. Then,
we build a Gaussian log-likelihood as the merit-of-function to minimize
through the equation $-2\ln(\mathcal{L})\varpropto\chi^{2}$, where $\chi$ is
the chi-square function given by

$\chi^{2}_{\rm
uncorr}=\sum_{i}^{N_{dat}}\left(\frac{\mathcal{D}-\mathcal{M}}{\sigma_{\mathcal{D}}}\right)^{2}\,,$
(48)

for the samples OHD, HIIG, SLS because the measurements are considered to be
uncorrelated. $N_{dat}$ is the number of points of dataset $\mathcal{D}$,
$\sigma_{\mathcal{D}}$ is the estimated uncertainty for each dataset, and
$\mathcal{M}$ represents the theoretical quantity of that observable based on
$E_{4}$ presented in (39) and (47) for $\Lambda\neq 0$ and $\Lambda=0$ models
respectively. As SNIa and BAO datasets contain correlated points, the figure
of merit is built as

$\chi^{2}_{\rm corr}=\Delta\vec{x}\cdot{\rm C}^{-1}\cdot\Delta\vec{x}^{T}\,,$
(49)

where $\Delta\vec{x}$ is the difference between the observational and
theoretical quantities, and $\rm C^{-1}$ is the covariance matrix. It is worth
to mention that a nuisance parameter is presented in the SNIa data and it is
convenient to marginalize over it to reduce the uncertainties. Thus, the
figure of merit for SNIa data is

$\chi^{2}_{\rm SNIa}=a+\ln\left(\frac{e}{2\pi}\right)-\frac{b^{2}}{e}\,,$ (50)

where $a$, $b$, and $e$ are functions of $\Delta\vec{x}$ and $\rm C^{-1}$. For
more details on these expressions see Motta et al. (2021).

Finally, we perform a joint analysis through the sum of the function-of-merits
of each data sample, namely

$\chi^{2}_{\rm Joint}=\chi^{2}_{\rm OHD}+\chi^{2}_{\rm SLS}+\chi^{2}_{\rm
HIIG}+\chi^{2}_{\rm SNIa}+\chi^{2}_{\rm BAO},$ (51)

where subscripts indicate the dataset under consideration.

###  3.2 Results

Performing the full confrontation of the scenario we construct the
corresponding log-likelihood contours at $68\%$ ($1\sigma$) and $99.7\%$
($3\sigma$) confidence level (CL), and we present them in Fig. 1 alongside the
1D posterior distribution. Moreover, in Table 1 we show the mean values and
the uncertainties at $1\sigma$ confidence level for the parameters $h$,
$\Omega_{m}^{(0)}$ and $\beta$ for both $\Lambda\neq 0$ and $\Lambda=0$ cases.

Figure 1:  Two-dimensional log-likelihood contours at 68% and 99.7% confidence level (CL), alongside the corresponding 1D posterior distribution of the free parameters, in Kaniadakis horizon entropy cosmology, for $\Lambda\neq 0$ (upper panel) and $\Lambda=0$ (lower panel). We use the various datasets described in the text, as well as the joint analysis.  Sample | $\chi^{2}_{\rm min}$ | $h$ | $\Omega_{m}^{(0)}$ | $\beta$ |  $\Delta$AICc |  $\Delta$BIC |  $\Delta$DIC
---|---|---|---|---|---|---|---
Case $\Lambda\neq 0$
OHD | $19.25$ | $0.699^{+0.033}_{-0.034}$ | $0.354^{+0.072}_{-0.061}$ | $-0.004^{+1.259}_{-1.255}$ | $7.6$ | $8.2$ | $-4.5$
BAO | $2.91$ | $0.599^{+0.272}_{-0.270}$ | $0.302^{+0.027}_{-0.023}$ | $-0.016^{+1.596}_{-1.594}$ | $14.1$ | $1.9$ | $0.2$
SLS | $216.52$ | $0.608^{+0.268}_{-0.277}$ | $0.077^{+0.064}_{-0.042}$ | $-0.006^{+2.550}_{-2.526}$ | $5.5$ | $8.3$ | $-5.2$
HIIG | $452.96$ | $0.722^{+0.018}_{-0.018}$ | $0.408^{+0.151}_{-0.137}$ | $0.043^{+2.510}_{-2.562}$ | $19.3$ | $22.3$ | $-19.6$
SNIa | $1042.99$ | $0.598^{+0.273}_{-0.270}$ | $0.359^{+0.126}_{-0.055}$ | $0.009^{+1.108}_{-1.124}$ | $9.0$ | $14.0$ | $-14.6$
Joint | $1743.48$ | $0.708^{+0.012}_{-0.011}$ | $0.283^{+0.016}_{-0.015}$ | $-0.011^{+0.517}_{-0.507}$ | $2.9$ | $8.1$ | $0.5$
Case $\Lambda=0$
OHD | $14.56$ | $0.701^{+0.029}_{-0.030}$ | $0.353^{+0.057}_{-0.050}$ | $1.138^{+0.043}_{-0.051}$ | $0.5$ | $0.0$ | $0.0$
BAO | $2.33$ | $0.602^{+0.272}_{-0.273}$ | $0.297^{+0.023}_{-0.021}$ | $1.186^{+0.017}_{-0.020}$ | $3.6$ | $-0.4$ | $-0.3$
SLS | $212.86$ | $0.596^{+0.276}_{-0.270}$ | $0.057^{+0.031}_{-0.027}$ | $1.373^{+0.020}_{-0.023}$ | $-0.2$ | $-0.3$ | $0.0$
HIIG | $435.64$ | $0.721^{+0.018}_{-0.018}$ | $0.298^{+0.050}_{-0.045}$ | $1.185^{+0.038}_{-0.043}$ | $-0.1$ | $-0.2$ | $-0.1$
SNIa | $1036.48$ | $0.596^{+0.275}_{-0.271}$ | $0.402^{+0.023}_{-0.023}$ | $1.093^{+0.021}_{-0.021}$ | $0.5$ | $0.5$ | $0.5$
Joint | $1753.03$ | $0.715^{+0.012}_{-0.012}$ | $0.326^{+0.015}_{-0.015}$ | $1.161^{+0.013}_{-0.013}$ | $10.4$ | $10.4$ | $10.4$
Table 1: Best-fit values and their $68\%$ CL uncertainties for Kaniadakis
horizon entropy cosmology with $\Lambda\neq 0$ (upper panel) and $\Lambda=0$
(lower panel) employing the data sets: OHD ($31$ data points), BAO (6 data
points), SLS (143 data points), HIIG (181 data points), SNIa (1048 data
points) and the joint analysis of them.

Figure 2: Upper panel, left to right: reconstruction of the $H(z)$, $q(z)$,
and $j(z)$, in Kaniadakis horizon entropy cosmology with $\Lambda\neq 0$.
Lower panel: same as before for the case $\Lambda=0$. We have used the bound
obtained from the joint analysis, and the shaded regions denote the
uncertainties at $1\sigma$. For completeness, the red square-points represent
the results of $\Lambda$CDM cosmology with $h=0.6766$ and
$\Omega_{m}^{(0)}=0.3111$ (Aghanim et al., 2020).

As can be seen, the bounds estimated from each data sample are consistent
among themselves, although the SLS dataset provides lower values for
$\Omega_{m}^{(0)}$. The joint constraints $h=0.708^{+0.012}_{-0.011}$
($h=0.715^{+0.012}_{-0.012}$) for the case $\Lambda\neq 0$ ($\Lambda=0$), are
consistent at $2.67\sigma$ ($3\sigma$) with the one estimated from the CMB
anisotropies (Aghanim et al., 2020) and at $1.74\sigma$ ($1.36\sigma$) with
the one from SH0ES (Riess et al., 2019). Hence, the scenario of Kaniadakis
horizon entropy cosmology can offer an alleviation to the $H_{0}$ tension,
providing a value in between its local measurements, and it indirect
estimation from the early stages of the Universe.

Concerning Kaniadakis parameter, we find that, when $\Lambda\neq 0$, the
combination of the data samples constrains $\beta=-0.011^{+0.515}_{-0.507}$,
namely $\beta$ is constrained around 0 as expected, the value in which
Kaniadakis entropy becomes the standard Bekenstein-Hawking one. However, when
$\Lambda=0$, the joint constraint yields $\beta=1.161^{+0.013}_{-0.013}$,
which is expected, as we mentioned above, because in the absence of an
explicit cosmological constant one needs a significant deviation from standard
cosmology to describe the Universe acceleration. Finally, note that due to
equation (43) that holds in the $\Lambda=0$ case, we acquire a correlation
between $\Omega_{m}^{(0)}$ and the Kaniadakis parameter $\beta$ in the lower
panel of Fig. 1.

Let us make a comment on the predicted entropy today, since this is possible
to be calculated through Eq. (2). According to our model, and imposing for the
horizon area of our Universe its present value, we arrive at the value
$S_{K}\sim 1.44\times 10^{99}$ m2Kg s-2K-1 (for $\Lambda\neq 0$) and
$S_{K}\sim 3.15\times 10^{99}$m2Kg s-2K-1 (for $\Lambda=0$). In comparison,
for the standard Bekenstein-Hawking entropy, we have $S_{BH}\sim 2.83\times
10^{99}$ m2Kg s-2K-1 ($\Lambda\neq 0$) and $S_{BH}\sim 2.79\times 10^{99}$m2Kg
s-2K-1 ($\Lambda=0$). Therefore, the corresponding ratio is
$S_{K}/S_{BH}\simeq 0.5$ for $\Lambda\neq 0$ and $S_{K}/S_{BH}\simeq 1.12$ for
$\Lambda=0$, which implies a small difference between Kaniadakis and
Bekenstein-Hawking entropies.

Due to the competitive qualities of the fits obtained from both scenarios, it
would be interesting to statistically compare them with the concordance
$\Lambda$CDM cosmology. In order to achieve this, we apply the standard
criteria, namely the Akaike information criterion corrected for small samples
(AICc, Akaike, 1974; Sugiura, 1978; Hurvich & Tsai, 1989) and the Bayesian
information criterion (BIC, Schwarz, 1978), since $\Lambda\neq 0$ model
contains one extra free parameter over $\Lambda$CDM. The AICc and BIC are
defined as ${\rm AICc}=\chi^{2}_{\text{min}}+2k+(2k^{2}+2k)/(N-k-1)$ and ${\rm
BIC}=\chi^{2}_{\text{min}}+k\ln(N)$ respectively, where
$\chi^{2}_{\text{min}}$ is the minimum of the $\chi^{2}$ function, $N$ is the
size of the dataset and $k$ is the number of free parameters. Following the
rules described in Hernández-Almada et al. (2022), we find that $\Lambda=0$
model and $\Lambda$CDM are statistically equivalent based on AICc
($\Delta$AICc$<4$), when the sample are treated separately, but show a strong
evidence against ($6<{\rm BIC}<10$) the scenario when the joint analysis is
applied. On the other hand, although AICc suggests that $\Lambda\neq 0$ model
and $\Lambda$CDM are statistically equivalent in the joint analysis, BIC
indicates that there is a strong evidence against the candidate model.
Additionally, for the two models we find that the $\Lambda=0$ case is
preferred by separate datasets, while the $\Lambda\neq 0$ case is
statistically preferred for the combined data analysis.

For completeness, we additionally calculate the Deviance information criterion
(DIC, Spiegelhalter et al., 2002; Kunz et al., 2006; Liddle, 2007). This is
defined as ${\rm DIC}=D(\bar{\theta})+2p_{D}$, where
$D(\bar{\theta})=\chi^{2}(\bar{\theta})$ is the Bayesian deviation,
$p_{D}=\bar{D}(\theta)-D(\bar{\theta})$ is the Bayesian complexity, which
represents the number of effective degrees of freedom, and $\bar{\theta}$ is
the mean value of the parameters. The advantage of DIC is its use of the full
log-likelihood sample instead of only the maximum log-likelihood (or minimum
$\chi^{2}$) as AICc and BIC do. Based on the Jeffreys scale (Jeffreys, 1961),
for $\Delta{\rm DIC}<2$ both models are statistical equivalent. In contrast,
$2<\Delta{\rm DIC}<6$ suggests a moderate tension between models, being the
one with lower value of DIC the best one, and $\Delta{\rm DIC}>10$ implies a
strong tension between the two models. We find that the $\Lambda\neq 0$ case
and $\Lambda$CDM scenario are statistical equivalent for BAO, they have a
moderate tension for OHD and SLS, and a strong tension for HIIG and SNIa. On
the other hand, the $\Lambda\neq 0$ case and $\Lambda$CDM are statistical
equivalent for OHD, BAO, HIIG, and SNIa. In summary, we confirm the results
obtained for the Joint analysis by AICc and BIC for both $\Lambda\neq 0$ and
$\Lambda=0$ models. It is worth to mention that when a posterior distribution
presents a bimodal shape or is asymmetric for a parameter, $p_{D}$ yields
negative values and thus DIC may not be a good criterion. This situation is
mainly presented for $\beta$ in the $\Lambda\neq 0$ case in separate datasets.

As a next step, we use the constraints from the joint analysis to reconstruct
the three cosmographic parameters, namely the Hubble, $H(z)$, the
deceleration, $q(z)$, and jerk, $j(z)$, parameters according to (24), (25).
The cosmic evolution of parameters is shown in Fig. 2. Thus, we report the
current values of $q_{0}=-0.610^{+0.028}_{-0.035}$
($-0.708^{+0.016}_{-0.016}$) for the deceleration parameter, and
$j_{0}=1.041^{+0.051}_{-0.036}$ ($1.137^{+0.014}_{-0.013}$) for the jerk
parameter for the $\Lambda\neq 0$ ($\Lambda=0$) scenario. Furthermore, the
transition redshift between the deceleration and the acceleration stages is
estimated to be $z_{T}=0.715^{+0.042}_{-0.041}$ ($0.652^{+0.032}_{-0.031}$),
which is in agreement with the one obtained by $\Lambda$CDM as shown in Fig.
2. Note that the jerk parameter evolution reveals the dynamical equation of
state of the effective dark energy.

Finally, to investigate in more detail the Hubble tension, we apply a new
diagnostic, called $\mathbf{\mathbb{H}}0(z)$ diagnostic, defined by (Krishnan
et al., 2021)

$\mathbf{\mathbb{H}}0(z)=\frac{H(z)}{E_{\rm\Lambda CDM}(z)}\,,$ (52)

where $H(z)$ is the Hubble function evolution in a given cosmological scenario
alternative to $\Lambda$CDM, and $E_{\rm\Lambda CDM}(z)$ is the dimensionless
Hubble parameter of $\Lambda$CDM paradigm. This diagnostic measures a possible
deviation of $H_{0}$ from its $\Lambda$CDM value. Concerning a flat
$\Lambda$CDM, a non-constant path of $\mathbf{\mathbb{H}}0(z)$ within error
bars suggests a modification of the Planck-$\Lambda$CDM scenario. In Fig. 3 we
depict the obtained results.

As we observe, there is an agreement within $1\sigma$ between
flat-$\Lambda$CDM cosmology and Kaniadakis cosmology for $z\gtrsim 0.7$ in the
$\Lambda\neq 0$ case, and for $0.7\lesssim z\lesssim 1.3$ in the $\Lambda=0$
case. Additionally, it is interesting that the current value
$\mathbf{\mathbb{H}}0(z=0)$ for both models are consistent with the one
obtained by SH0ES (Riess et al., 2019), and that the $\Lambda\neq 0$ model has
a trend to the Planck value in the past (Aghanim et al., 2020). This is
another verification that the scenario of Kaniadakis horizon entropy cosmology
may offer an alleviation to the $H_{0}$ tension. Nevertheless, to further
investigate whether both Kaniadakis models can alleviate the Hubble tension, a
parameter estimation using the linear perturbation equations together with CMB
data should be performed.

We close this section by investigating one important process in every
cosmological scenario: the Big Bang Nucleosynthesis (BBN), since the
production of light elements in the early Universe can be affected in non-
standard cosmologies (Pospelov & Pradler, 2010; Barrow et al., 2021b).
Considering that the freeze-out of the light elements occurs when the weak
interaction rates are lower than $H(z)$, a simple test to guarantee that the
BBN is not spoiled is to require that the deviation $\delta H(z)$ with respect
to the standard Hubble expansion rate at the BBN epoch should be small.
Although in the Friedmann equations mentioned above we have not included a
radiation component, this can be added and we can perform the analysis by
expanding $E_{4}$ around $\beta=0$ in Eq. (39) (resp. Eq. (47)) and neglecting
the fourth order error terms, resulting to

$\displaystyle
E(z)=\left\\{\begin{array}[]{c}\underbrace{\sqrt{\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}}}_{\text{value
from}\;\Lambda{\rm
CDM}}+\underbrace{\frac{\beta^{2}}{4}\left[\Omega_{m}^{(0)}(z+1)^{{3}}+\Omega_{\Lambda}^{(0)}\right]^{-3/2}}_{\text{correction}},\;\Lambda\neq
0\\\ \underbrace{\sqrt{\Omega_{m}^{(0)}}(z+1)^{{3/2}}}_{\text{value
from}\;\Lambda{\rm
CDM}}+\underbrace{\frac{\beta^{2}}{4}\left[\Omega_{m}^{(0)}\right]^{-3/2}(z+1)^{{-1/2}}}_{\text{correction}},\;\Lambda=0\end{array}\right..$
(55)

Hence, we can study both Kaniadakis models ($\Lambda\neq 0$ and $\Lambda=0$)
at $z\sim 10^{10}$ (approximately BBN era). We find that for $\Lambda\neq 0$,
the model is consistent with the Big Bang Nucleosynthesis (BBN) constraints,
since the correction term at $z\sim 10^{10}$ is of the order of $\sim
10^{-49}$, dominating the standard cosmology and not producing significant
effects in the formation of light elements. In the case $\Lambda=0$, the
correction is larger, and calculations at $z\sim 10^{10}$ are of the order
$\sim 10^{-5}$. However, such corrections are still subdominant, allowing the
production of light elements. A further analysis could be performed following
Capozziello et al. (2017); Barrow et al. (2021b); Asimakis et al. (2021).

Figure 3: The $\mathbf{\mathbb{H}}0(z)$ diagnostic for Kaniadakis horizon
entropy cosmology with $\Lambda\neq 0$ (upper panel) and $\Lambda=0$ (lower
panel). We have used the bound obtained from the joint analysis, and the
shaded regions denote the uncertainties at $1\sigma$. For completeness, the
red square-points represent the results of $\Lambda$CDM cosmology with
$h=0.6766$ and $\Omega_{m}^{(0)}=0.3111$ (Aghanim et al., 2020).

##  4 Dynamical system and stability analysis

In this section we perform a full dynamical system analysis in order to
investigate the global dynamics of cosmological scenarios, and obtain
information on the Universe evolution independently of the initial conditions.
In the dynamical system formulation, one starts from local analysis of the
differential equation $\mathbf{x}^{\prime}(\tau)={\bf X}(\mathbf{x})$, where
$\mathbf{x}$ is the state vector, and $\tau$ a convenient time variable, near
an equilibrium point $\mathbf{x}=\bar{\mathbf{x}}$, and progressively extends
the investigated regions of the phase and of the parameter space. Assuming
that the vector field ${\bf X}(\mathbf{x})$ has continuous partial
derivatives, the process of determining the local behavior is based on the
linear approximation of the vector field ${{\bf X}(\mathbf{x})\approx{\bf
DX}(\bar{\mathbf{x}})(\mathbf{x}-\bar{\mathbf{x}})}$ where ${\bf
DX}(\bar{\mathbf{x}})$ is the Jacobian of the vector field at the equilibrium
point $\bar{\mathbf{x}}$, which is referred to as the linearization of the
dynamical equations at the equilibrium point. In this neighborhood we acquire
the system $\mathbf{x}^{\prime}(\tau)={{\bf
DX}(\bar{\mathbf{x}})(\mathbf{x}-\bar{\mathbf{x}})}$. Each of the equilibrium
points can be classified according to the real parts of the eigenvalues of
${\bf DX}(\bar{\mathbf{x}})$ (if none of these are zero). Thus, this approach
provides a general description of the phase space of all possible solutions of
the system, their equilibrium points and stability, as well as the asymptotic
solutions (Wainwright & Ellis, 1997; Ferreira & Joyce, 1997; Copeland et al.,
1998; Perko, 2000; Coley, 2003; Copeland et al., 2006; Chen et al., 2009;
Cotsakis & Kittou, 2013; Giambo & Miritzis, 2010; Papagiannopoulos et al.,
2022). If some real parts of the eigenvalues are zero, the equilibrium point
is nonhyperbolic, and the analysis through linearization fails. Then, we use
numerical tools for the analysis.

In the following subsections we perform the global dynamical system analysis
for the two cases, namely $\Lambda\neq 0$ and $\Lambda=0$.

###  4.1 Case I: $\Lambda\neq 0$

Defining the dimensionless variables $\theta,T$ as

$\displaystyle\theta=\arctan\left(1-\frac{8\pi
G\rho_{DE}}{3H^{2}}\right),\;\theta\in\left[-\frac{\pi}{2},\frac{\pi}{2}\right],\quad
T=\frac{H_{0}}{H+H_{0}},$ (56)

with

$\Omega_{m}:=\frac{8\pi G\rho_{m}}{3H^{2}}=\tan\left(\theta\right),$ (57)

then equation (10) becomes

$\displaystyle\frac{\Lambda}{3H_{0}^{2}}=\frac{(1-T)^{2}\left\\{\cosh\left[\frac{T^{2}\beta}{(1-T)^{2}}\right]-\tan(\theta)\right\\}}{T^{2}}-\beta\text{shi}\left[\frac{T^{2}\beta}{(1-T)^{2}}\right].$
(58)

It proves convenient to introduce the new time derivative as

$\displaystyle f^{\prime}\equiv\frac{df}{d\tau}=\frac{\cosh\left(\frac{\pi
K}{GH^{2}}\right)}{H}\dot{f}.$ (59)

Therefore, we finally extract the dynamical system

$\displaystyle{\theta}^{\prime}(\tau)=3\sin(\theta)\left\\{\sin(\theta)-\cos(\theta)\cosh\left[\frac{T^{2}\beta}{(1-T)^{2}}\right]\right\\},$
(60) $\displaystyle T^{\prime}(\tau)=\frac{3}{2}(1-T)T\tan(\theta).$ (61)

Note that this system diverges at $T=1$ and at $\theta=\pm\pi/2$.

Lastly, the deceleration parameter (24) is written as

$q:=-1-\frac{\dot{H}}{H^{2}}=-1+\frac{3}{2}\tan(\theta)\text{sech}\left[\frac{\beta
T^{2}}{(1-T)^{2}}\right].$ (62)

Note that, for an expanding universe ($H>0$), we have that $T\in[0,1]$, while
$\theta$ is a periodic coordinate with period $\pi$, and thus we can set
$\theta\in\left[-{\pi}/{2},{\pi}/{2}\right]$ (modulo a periodic shift
$c\pi,\;c\in\mathbb{Z}$). Moreover, the physical condition
$0\leq\Omega_{m}\leq 1$ implies that the region of physical interest is
$\theta\in[0,\pi/4]$ (modulo a periodic shift $c\pi,\;c\in\mathbb{Z}$). The
non-physical region $\Omega_{m}>1$ is $\theta\in(\pi/4,\pi/2]$ (modulo a
periodic shift $c\pi,\;c\in\mathbb{Z}$). Hence, we have obtained a global
phase-space formulation. For the representation of the flow of (60) and (61),
we integrate in the variables $T,\theta$ and project in a compact set using
the “cylinder-adapted” coordinates

$\mathbf{S}:\begin{cases}x=\cos(\theta),\\\ y=\sin(\theta),\\\
z=T,\end{cases}$ (63)

with $0\leq T\leq 1,\theta\in[-\pi,\pi]$, with inverse
$\theta=\arctan\left({y}/{x}\right)$, and $T=z$. Thus, the region of physical
interest is $\theta\in[0,\pi/4]$, modulo a periodic shift
$c\pi,\;c\in\mathbb{Z}$.

Label | $\theta$ | $T$ | Existence | Stability
---|---|---|---|---
$dS_{+}$ | $2\pi c_{1}$ | arbitrary | $c_{1}\in\mathbb{Z}$ | stable
$dS_{+}^{(0)}$ | $2\pi c_{1}$ | $0$ | $c_{1}\in\mathbb{Z}$ | stable
$dS_{+}^{(1)}$ | $2\pi c_{1}$ | $1$ | $c_{1}\in\mathbb{Z}$ | stable
$dS_{-}$ | $\pi(1+2c_{1})$ | arbitrary | $c_{1}\in\mathbb{Z}$ | stable
$dS_{-}^{(0)}$ | $\pi(1+2c_{1})$ | $0$ | $c_{1}\in\mathbb{Z}$ | stable
$dS_{-}^{(1)}$ | $\pi(1+2c_{1})$ | $1$ | $c_{1}\in\mathbb{Z}$ | stable
$M_{-}^{(0)}$ | $2\pi c_{1}-\frac{3\pi}{4}$ | $0$ | $c_{1}\in\mathbb{Z}$ | unstable
$M_{+}^{(0)}$ | $2\pi c_{1}+\frac{\pi}{4}$ | $0$ | $c_{1}\in\mathbb{Z}$ | unstable
$M_{-}^{(1)}$ | $2\pi c_{1}-\frac{3\pi}{4}$ | $1$ | $c_{1}\in\mathbb{Z},\beta=0$ | saddle
$M_{+}^{(1)}$ | $2\pi c_{1}+\frac{\pi}{4}$ | $1$ | $c_{1}\in\mathbb{Z},\beta=0$ | saddle
Table 2: The equilibrium points of the dynamical system (60) and (61) of
Kaniadakis horizon entropy cosmology with $\Lambda\neq 0$. We use $dS$ to
denote the de Sitter, dark-energy dominated solutions, and $M$ to denote the
matter-dominated ones.

We proceed by extracting the equilibrium points and characterizing their
stability. There are two equivalent hyperbolic equilibrium points $M_{\pm}$
for which $q=1/2$, i.e. they are associated with matter domination, and two
equilibrium points $dS_{\pm}$ corresponding to dark-energy dominated de-Sitter
solutions for which $q=-1$. The equilibrium points of the dynamical system
(60) and (61) of Kaniadakis horizon entropy cosmology with $\Lambda\neq 0$ are
presented in Table 2.

Figure 4: Phase-space plot of the dynamical system (60) and (61) of Kaniadakis
horizon entropy cosmology with $\Lambda\neq 0$, for the best-fit value of
Kaniadakis parameter obtained by the observational analysis, namely for
$\beta=-0.011$. Upper panel: unwrapped solution space. Lower panel: projection
over the cylinder $\mathbf{S}$ defined in Cartesian coordinates $(x,y,z)$
through (63). At late times the Universe results in a dark-energy dominated,
de Sitter solution, while the past attractor is the matter-dominated epoch.

In Fig. 4 we display an unwrapped solution space of the system (60) and (61)
(upper panel), and the projection over the cylinder $\mathbf{S}$, defined in
Cartesian coordinates $(x,y,z)$ by (63), for the best fit value $\beta=-0.011$
obtained through the observational analysis. For the points that are non-
hyperbolic, their stability is analyzed numerically. The two dashed lines,
indicated by $dS_{-}$ (blue) and $dS_{+}$ (red), are the late-time de Sitter
attractors. The early-time attractors are $M_{\pm}^{(0)}$ for which $q=1/2$,
and they correspond to matter-dominated solutions. Hence, at late times the
Universe results in a dark-energy dominated solution, while the past attractor
of the Universe is the matter-dominated epoch. At the intersection of the
invariant set $T=1$ with the singular lines $\theta=\pm\pi/2$ we obtain the
equilibrium points $L_{\pm}$. Considering that equations (60) and (61) diverge
at $L_{\pm}$, we should introduce suitable variables for the analysis.

For the analysis at $T=1$ it proves convenient to define the variable

$\Phi=\left\\{1+\exp\left[\frac{|\beta|T^{2}}{(1-T)^{2}}\right]\right\\}^{-1},\Phi\in[0,1],$
(64)

as well as the time rescaling

$f^{\prime}\equiv\frac{df}{d\zeta}=(1-\Phi)^{2}\frac{df}{d\eta}=\frac{\tanh^{2}\left(\frac{\pi|K|}{2GH^{2}}\right)+1}{4H}\dot{f}.$
(65)

Hence, using also the variable $\theta$ from (56), we finally obtain the
autonomous system

$\displaystyle\theta^{\prime}(\zeta)=-\frac{3}{2}\sin(\theta)\left[2(\Phi-1)\Phi(\sin(\theta)+\cos(\theta))+\cos(\theta)\right],$
(66)
$\displaystyle\Phi^{\prime}(\zeta)=3(1-\Phi)^{2}\Phi^{2}\tan(\theta)\ln\left[\Phi/(1-\Phi)\right].$
(67)

Figure 5: Phase-space plot of the system (66)-(67) of Kaniadakis horizon
entropy cosmology with $\Lambda\neq 0$, for dust matter. The late attractor
corresponds to $\theta=0$, and thus to a dark-energy domimated solution with
$\Omega_{DE}=1$.

In Fig. 5 we depict the phase-space flow of the system (66)-(67).
Asymptotically, $\theta\rightarrow 0$ and $\Phi$ tends to a constant
$\Phi_{0}$. Therefore, the late attractor corresponds to the dark-energy
dominated solution with $\Omega_{DE}=1$. The current values
$\Phi_{0}=1/(e^{\left|\beta\right|}+1)$,
$\theta_{0}=\text{arctan}\left(\Omega_{m}^{(0)}\right)$ leads to the de Sitter
solution $a(t)=e^{H_{0}(t-t_{U})}$.

####  4.1.1 Heteroclinic sequences

In the phase portrait of a dynamical system, a heteroclinic orbit is a path in
phase space that joins two different equilibrium points. If the equilibrium
points at the start and end of the orbit are the same, the orbit is a
homoclinic orbit (Guckenheimer & Holmes, 1983).

From the above analysis we can see that the invariant sets $T=0$ and $T=1$ are
of interest in the determination of possible heteroclinic sequences. The
direction of the flow can be determined by considering the monotonic function

$M_{1}=\frac{T}{1-T},\quad M_{1}^{\prime}(\tau)=\frac{3\tan(\theta)}{2}M_{1}.$
(68)

In if $\tan(\theta)<0$, the orbits move from $T=1$ to $T=0$, and if
$\tan(\theta)>0$, the orbits move from $T=0$ to $T=1$. In the invariant
manifold $T=0$ ($H\rightarrow\infty$) the dynamics is given by the one-
dimensional flow

${\theta}^{\prime}(\tau)=3\sin(\theta)\left(\sin(\theta)-\cos(\theta)\right).$
(69)

In Fig. 6 we present the one-dimensional dynamical system (69), in which we
can see the heteroclinic sequences $M_{+}^{(0)}\to dS_{+}^{(0)}$ and
$M_{-}^{(0)}\to dS_{-}^{(0)}$.

Figure 6: Phase-space diagram of the one-dimensional dynamical system (69) of
Kaniadakis horizon entropy cosmology with $\Lambda\neq 0$, for dust matter and
any value $\beta$.

Similarly, analyzing the one-dimensional flow in the invariant set $T=1$, that
corresponds to $\Phi=0$, we find that the dynamics on this invariant set is
given by the one-dimensional dynamical system

$\theta^{\prime}(\zeta)=-\frac{3}{2}\sin(\theta)\cos(\theta),$ (70)

Figure 7: Phase-space diagram of the one-dimensional dynamical system (70) of
Kaniadakis horizon entropy cosmology with $\Lambda\neq 0$, for dust matter and
any value $\beta$.

which has a behavior shown in Fig. 7, where the heteroclinic sequences
$L_{+}\to dS_{+}^{(1)}$ and $L_{-}\to dS_{-}^{(1)}$ are presented. Finally, to
find heteroclinic sequences $M_{\pm}^{(0)}\to dS_{\pm}^{(1)}$, the
intersection of the unstable manifold of $M_{\pm}^{(0)}$ with the stable
manifold of $dS_{\pm}^{(1)}$ should be analyzed. Since the former is
$\mathbb{R}^{2}$, then it is required to examine the stable manifold of
$dS_{\pm}^{(1)}$. This is given locally by the graph

$\left\\{(\Phi,\theta)\in\mathbb{R}^{2}:\Phi=h(\theta),h(0)=0,h^{\prime}(0)=0\right\\},|\theta|<\delta,$
(71)

for $\delta>0$ suitably small. By the invariance of the stable manifold we
obtain the quasilinear differential equation for $h$ given by

$\displaystyle\frac{3}{2}\sin(\theta)h^{\prime}(\theta)(2(h(\theta)-1)h(\theta)(\sin(\theta)+\cos(\theta))+\cos(\theta))$
$\displaystyle+3(h(\theta)-1)^{2}h(\theta)^{2}\tan(\theta)\ln\left(\frac{h(\theta)}{1-h(\theta)}\right)=0.$
(72)

Introducing the ansatz
$h(\theta)=a_{1}\theta^{2}+a_{2}\theta^{3}+a_{3}\theta^{4}+\ldots,$ we obtain
$a_{i}=0$ at any order. Therefore, the dynamics at the stable manifold of
$dS_{+}^{(1)}$ is given by equation (70). Then, it is easy to construct
heteroclinic sequences $M_{+}^{(0)}\to dS_{+}^{(1)}$ which pass near the
singularity $L_{+}$ by assuming for instance the initial value
$(\Phi,\theta)=(\varepsilon,\pi/4)$, $\varepsilon\approx 0$ and evolving the
system back and forward in $\zeta$. Similar arguments can be used to construct
heteroclinic sequences $M_{-}^{(0)}\to dS_{-}^{(1)}$, which pass near the
singularity $L_{-}$, with the initial value
$(\Phi,\theta)=(\varepsilon,-3\pi/4)$, $\varepsilon\approx 0$.

Summarizing, for $0\leq\theta\leq\pi/2$ (the physical region is
$0\leq\theta\leq\pi/4$), there exists the heteroclinic sequences
$M_{+}^{(0)}\;(\Omega_{m}\rightarrow 1,H\rightarrow\infty)\to
L_{+}(\Omega_{m}\rightarrow+\infty,H\rightarrow 0)\to dS_{+}^{(1)}(\text{de
Sitter},\Omega_{m}\rightarrow 0,H\rightarrow 0)$ and $M_{+}^{(0)}\to
dS_{+}^{(0)}(\text{de Sitter},\Omega_{m}\rightarrow 0,H\rightarrow\infty)$,
and in the region $-\pi\leq\theta\leq-\pi/2$ (the physical region is
$-\pi\leq\theta\leq-3\pi/4$), there exists the heteroclinic sequences
$M_{-}^{(0)}(\Omega_{m}\rightarrow 1,H\rightarrow\infty)\to
L_{-}(\Omega_{m}\rightarrow-\infty,H\rightarrow 0)\to dS_{-}^{(1)}(\text{de
Sitter},\Omega_{m}\rightarrow 0,H\rightarrow 0)$ and $M_{-}^{(0)}\to
dS_{-}^{(0)}(\text{de Sitter},\Omega_{m}\rightarrow 0,H\rightarrow\infty)$.

####  4.1.2 Bounce and a turnaround

Another interesting cosmological possibility is the possible existence of a
bounce and a turnaround (Saridakis, 2009; Cai et al., 2012; Zhu et al., 2021).
Let us assume that, for the state vector $(a,H,R)$, the field equations can be
written as

$\displaystyle\dot{a}=aH,$ (73)
$\displaystyle\dot{H}=\frac{1}{6}\left(R-12H^{2}\right),$ (74)
$\displaystyle\dot{R}=g(a,H,R),$ (75)

such that the function $g(a,H,R)$ satisfies $g(a,H,R)=-g(a,-H,R)$. Hence, the
system (73), (74) and (75) is invariant under time inversion $t\mapsto-t$ if
also $H\mapsto-H$ and $R\mapsto R$, and by definition $a\geq 0$. Those
solutions can be related to symmetric cyclic solutions with respect to the
origin, chosen to correspond to the possible bounce point
$t_{\text{bounce}}=0$. Therefore, if the bounce exists, the system (73), (74)
and (75) is a reversible system in the sense that it has a reversing symmetry
under time inversion.

Let us consider the simplest case where there is exactly one bouncing and
exactly one turnaround point. Note that both at the bounce and turnaround
points we have $H=0$. In this case, the line connecting these points and
corresponding to $H=0$ defines a plane that separates all points on the
trajectory in this phase space to the ones corresponding to either the
expanding ($H>0$) or contracting ($H<0$) phase. As discussed in Pavlović &
Sossich (2021), it is natural that in cyclic models the value of the Ricci
scalar would approach its maximum around the bounce and since $\dot{H}>0$,
from (74) it follows that this maximum Ricci scalar value is positive, and
moreover that $\ddot{H}=0$ at the bounce.

In summary, at the bounce we have $R=R_{\text{bounce}}>0$, $H=0$, and
$a=a_{\text{min}}$. The bounce is then followed by a phase in which
$\dot{H}>0$, $H>0$, $\dot{a}>0$, and $\dot{R}<0$. Then the Universe enters the
phase characterized by $\dot{H}<0$ and approaches the turnaround point which
is determined by $H=0$, $a=a_{\text{max}}$, and $R=R_{\text{turnaround}}<0$,
where the last condition follows from Eq. (74).

To obtain $g(a,H,R)$ in Eq. (75), we use Eqs. (8), (9), (11), (12), and (14),
where for simplicity we focus on the dust case. Therefore, we obtain

$\displaystyle R=6\dot{H}+12H^{2}=-24\pi Gp_{DE}$
$\displaystyle=3\Bigg{\\{}-8\pi G\rho_{m}\left[\text{sech}\left(\frac{\beta
H_{0}^{2}}{H^{2}}\right)-1\right]$ $\displaystyle+3\beta
H_{0}^{2}\text{shi}\left(\frac{H_{0}^{2}\beta}{H^{2}}\right)-3H^{2}\left[\cosh\left(\frac{\beta
H_{0}^{2}}{H^{2}}\right)-1\right]+\Lambda\Bigg{\\}}.$ (76)

Introducing the dimensionless quantities $E=\frac{H}{H_{0}}$ and
$\mathcal{R}=\frac{R}{12H_{0}^{2}}$, and using $z$ as the independent
variable, we extract the general system for $(a,E,\mathcal{R})$:

$\displaystyle\frac{da}{dz}=-a^{2},$ (77)
$\displaystyle\frac{dE}{dz}=-2\left(\mathcal{R}-E^{2}\right)\frac{a}{E},$ (78)
$\displaystyle\frac{d\mathcal{R}}{dz}=-\frac{9\beta{\Omega_{m}^{(0)}}^{2}\tanh\left(\frac{\beta}{E^{2}}\right)\text{sech}^{2}\left(\frac{\beta}{E^{2}}\right)}{4E^{4}a^{5}}.$
(79)

Finally, in order to examine whether the above requirements are fulfilled in
the present scenario, we use the best fit values $\beta=-0.011$ and
$\Omega_{m}^{(0)}=0.283$ in (77), (78) and (79), and we find

$\displaystyle\frac{da}{dz}=-a^{2},$ (80)
$\displaystyle\frac{dE}{dz}=-2\left(\mathcal{R}-E^{2}\right)\frac{a}{E},$ (81)
$\displaystyle\frac{d\mathcal{R}}{dz}=-\frac{0.0019822\tanh\left(\frac{0.011}{E^{2}}\right)\text{sech}^{2}\left(\frac{0.011}{E^{2}}\right)}{E^{4}a^{5}}.$
(82)

In the case of dust matter and $\Lambda\neq 0$. This system cannot satisfy the
above requirements, and hence the present scenario cannot exhibit bounce and
turnaround solutions.

###  4.2 Case II: $\Lambda=0$

In the case $\Lambda=0$, eq. (10) becomes

$\frac{8\pi G\rho_{m}}{3H^{2}}=-\frac{3\pi
K\text{shi}\left(\frac{K\pi}{GH^{2}}\right)}{3GH^{2}}+\cosh\left(\frac{\pi
K}{GH^{2}}\right).$ (83)

Figure 8: Phase-space diagram of the one-dimensional dynamical system (84) of
Kaniadakis horizon entropy cosmology with $\Lambda=0$, for dust matter and for
the best-fit value of Kaniadakis parameter obtained by the observational
analysis, namely for $\beta=1.161$. The physical region is $0\leq T<1$. The
equilibrium point $T=0$ is unstable, dominated by dark-energy, and de Sitter
equilibrium point $T=T_{c}\approx 0.521$ is stable.

This expression is used as a definition of $\rho_{m}$. If $\beta\neq 0$, and
re-scaling the time derivative $d/d\nu=(1-T)d/d\tau$, we obtain

$\displaystyle T^{\prime}(\nu)=\frac{3}{2}(1-T)^{2}T\cosh\left(\frac{\beta
T^{2}}{(1-T)^{2}}\right)-\frac{3}{2}\beta
T^{3}\text{shi}\left(\frac{T^{2}\beta}{(1-T)^{2}}\right).$ (84)

The equilibrium points of (84) are $T=0$, which is unstable, and the
equilibrium point $T=T_{c}$, where $T_{c}$ is a solution of the transcendental
equation  $(1-T_{c})^{2}\cosh\left({\beta
T_{c}^{2}}/{(1-T_{c})^{2}}\right)-\beta
T_{c}^{2}\text{shi}\left({T_{c}^{2}\beta}/{(1-T_{c})^{2}}\right)=0$ ,
$0<T_{c}<1$, corresponding to de Sitter solution $a(t)\propto
e^{H_{0}t\left(\frac{1}{T_{c}}-1\right)}$, is stable. In Fig. 8 we depict a
phase-space plot of the one-dimensional dynamical system (84) of Kaniadakis
horizon entropy cosmology with $\Lambda=0$, for dust matter, and the join
value $\beta=1.161$. Note that all orbits originate from the invariant subset
$T=0$, classically related to the initial singularity with
$H\rightarrow\infty$. The late-time attractor is $T=T_{c}\approx 0.521$, and
it corresponds to de Sitter solution.

Finally, in order to examine whether the present scenario exhibits a bounce,
we use the best fit values $\beta=1.161$ and $\Omega_{m}^{(0)}=0.326$ in (77),
(78) and (79) and we find

$\displaystyle\frac{da}{dz}=-a^{2},$ (85)
$\displaystyle\frac{dE}{dz}=-2\left(\mathcal{R}-E^{2}\right)\frac{a}{E},$ (86)
$\displaystyle\frac{d\mathcal{R}}{dz}=\frac{0.277619\tanh\left(\frac{1.161}{E^{2}}\right)\text{sech}^{2}\left(\frac{1.161}{E^{2}}\right)}{E^{4}a^{5}},$
(87)

in the case of dust matter and $\Lambda=0$, we deduce that the system cannot
fulfill the bounce requirements, and therefore it cannot exhibit bounce and
turnaround solutions.

##  5 Summary and discussion

The present work was devoted to explore Kanadiakis horizon entropy cosmology
which arises from the application of the gravity-thermodynamics conjecture
using the Kaniadakis modified entropy. The resulting modified Friedmann
equations contain extra terms that constitute an effective dark energy sector.
Moreover, we used data from Observational Hubble Data, Supernova Type Ia, HII
galaxies, Strong Lensing Systems, and Baryon Acoustic Oscillations
observations, and we applied a Bayesian Markov Chain Monte Carlo analysis in
order to construct the likelihood contours for the model parameters.

Regarding the Kaniadakis parameter $\beta$, we found that it is constrained
around 0, namely around the value in which standard Bekenstein-Hawking is
recovered. Furthermore, the present matter density parameter
$\Omega_{m}^{(0)}$ is consistent with the expected value from $\Lambda$CDM
scenario, having a lower value for the $\Lambda\neq 0$ case and a slightly
higher value for the $\Lambda=0$ case.

However, the interesting result comes from the constraint on the normalized
Hubble parameter $h$. In particular, for $\Lambda\neq 0$ we extracted
$h=0.708^{+0.012}_{-0.011}$ while for $\Lambda=0$ we found
$h=0.715^{+0.012}_{-0.012}$. Thus, the obtained value of $H_{0}$ for
$\Lambda\neq 0$ deviates $2.67\sigma$ from the Planck value and $1.74\sigma$
from the SH0ES one, while in the $\Lambda=0$ case the deviation is $3\sigma$
from the Planck value and $1.36\sigma$ from the SH0ES one. Additionally, in
order to verify this result in an independent way, we performed the
$\mathbf{\mathbb{H}}0(z)$ diagnostic. Hence, our analysis reveals that
Kaniadakis horizon entropy cosmology is an interesting candidate to alleviate
the $H_{0}$ tension problem. This is one of the main results of the present
work.

We proceeded by investigating the cosmographic parameters, namely the
deceleration and jerk ones, by using the data in order to reconstruct them in
the redshift region $0<z<2.5$. As we showed, the transition from deceleration
to acceleration happens at $z_{T}=0.715^{+0.042}_{-0.041}$ for the
$\Lambda\neq 0$ case and at $z_{T}=0.652^{+0.032}_{-0.031}$ for the
$\Lambda=0$ case, in agreement within $1\sigma$ with that found in Herrera-
Zamorano et al. (2020) for $\Lambda$CDM cosmology. Furthermore, we applied the
AICc and BIC information criteria and we found that although AICc suggests
that $\Lambda\neq 0$ model and $\Lambda$CDM are statistically equivalent in
the joint analysis, BIC indicates that there is a strong evidence against the
candidate model. Lastly, applying the DIC criterion we found that the
$\Lambda\neq 0$ case and $\Lambda$CDM are statistical equivalent for BAO, they
have a moderate tension for OHD and SLS, and a strong tension for HIIG and
SNIa datasets, while the $\Lambda\neq 0$ case and $\Lambda$CDM are statistical
equivalent for all datasets.

Finally, we performed a detailed dynamical-system analysis, providing a
general description of the phase-space of all possible solutions of the
system, their equilibrium points and stability, as well as the late-time
asymptotic behavior. As we showed, the Universe past attractor is the matter-
dominated epoch, while at late times the Universe results in the dark-energy-
dominated solution, for both $\Lambda=0$, and $\Lambda\neq 0$ cases. Moreover,
we showed that the scenario accepts heteroclinic sequences, but it cannot lead
to bounce and turnaround solutions.

In summary, the scenario of Kaniadakis horizon entropy cosmology exhibits very
interesting phenomenology and is in agreement with observational behavior.
Hence, it can be an interesting candidate for the description of Nature.

## Acknowledgments

We thank the anonymous referee for thoughtful remarks and suggestions. A.H.A.
thanks to the PRODEP project, Mexico for resources and financial support and
thanks also to the support from Luis Aguilar, Alejandro de León, Carlos
Flores, and Jair García of the Laboratorio Nacional de Visualización
Científica Avanzada. G.L. was funded by Agencia Nacional de Investigación y
Desarrollo - ANID for financial support through the program FONDECYT
Iniciación grant no. 11180126 and by Vicerrectoría de Investigación y
Desarrollo Tecnológico at UCN. J.M. acknowledges the support from ANID project
Basal AFB-170002 and ANID REDES 190147\. M.A.G.-A. acknowledges support from
Universidad Iberoamericana that support with the SNI grant, ANID REDES
(190147), Cátedra Marcos Moshinsky and Instituto Avanzado de Cosmología (IAC).
V.M. acknowledges support from Centro de Astrofísica de Valparaíso and ANID
REDES 190147. This work is partially supported by the Ministry of Education
and Science of the Republic of Kazakhstan, Grant AP08856912. A.D. Millano was
supported by Agencia Nacional de Investigación y Desarrollo - ANID-
Subdirección de Capital Humano/Doctorado Nacional/año 2020- folio 21200837 and
by Vicerrectoría de Investigación y Desarrollo Tecnológico at UCN.

## Data Availability

The data underlying this article were cited in Section 3.1.

## References

  * Abreu & Ananias Neto (2021) Abreu E. M. C., Ananias Neto J., 2021, EPL, 133, 49001
  * Abreu et al. (2016) Abreu E. M. C., Ananias Neto J., Barboza E. M., Nunes R. C., 2016, EPL, 114, 55001
  * Abreu et al. (2018) Abreu E. M. C., Neto J. A., Mendes A. C. R., Bonilla A., 2018, EPL, 121, 45002
  * Addazi et al. (2021) Addazi A., et al., 2021 (arXiv:2111.05659)
  * Aghanim et al. (2020) Aghanim N., et al., 2020, Astron. Astrophys., 641, A6
  * Akaike (1974) Akaike H., 1974, IEEE Transactions on Automatic Control, 19, 716
  * Akbar & Cai (2006) Akbar M., Cai R.-G., 2006, Phys. Lett. B, 635, 7
  * Akbar & Cai (2007) Akbar M., Cai R.-G., 2007, Phys. Rev. D, 75, 084003
  * Amante et al. (2020) Amante M. H., Magaña J., Motta V., García-Aspeitia M. A., Verdugo T., 2020, MNRAS, 498, 6013–6033
  * Asimakis et al. (2021) Asimakis P., Basilakos S., Mavromatos N. E., Saridakis E. N., 2021, Big Bang Nucleosynthesis constraints on higher-order modified gravities (arXiv:2112.10863)
  * Bak & Rey (2000) Bak D., Rey S.-J., 2000, Class. Quant. Grav., 17, L83
  * Barrow (2020) Barrow J. D., 2020, Phys. Lett. B, 808, 135643
  * Barrow et al. (2021a) Barrow J. D., Basilakos S., Saridakis E. N., 2021a, Phys. Lett. B, 815, 136134
  * Barrow et al. (2021b) Barrow J. D., Basilakos S., Saridakis E. N., 2021b, Physics Letters B, 815, 136134
  * Beutler et al. (2011) Beutler F., et al., 2011, MNRAS, 416, 3017
  * Blake et al. (2011) Blake C., et al., 2011, MNRAS, 418, 1707
  * Cai & Cao (2007) Cai R.-G., Cao L.-M., 2007, Phys. Rev. D, 75, 064008
  * Cai & Kim (2005) Cai R.-G., Kim S. P., 2005, JHEP, 02, 050
  * Cai & Ohta (2010) Cai R.-G., Ohta N., 2010, Phys. Rev. D, 81, 084061
  * Cai et al. (2009) Cai R.-G., Cao L.-M., Hu Y.-P., 2009, Class. Quant. Grav., 26, 155018
  * Cai et al. (2012) Cai Y.-F., Gao C., Saridakis E. N., 2012, JCAP, 10, 048
  * Capozziello et al. (2017) Capozziello S., Lambiase G., Saridakis E. N., 2017, Eur. Phys. J. C, 77, 576
  * Carroll (2001) Carroll S. M., 2001, Living Rev. Rel., 4, 1
  * Chen et al. (2009) Chen X.-m., Gong Y.-g., Saridakis E. N., 2009, JCAP, 0904, 001
  * Chen et al. (2019) Chen L., Huang Q.-G., Wang K., 2019, J. Cosmology Astropart. Phys., 2019, 028
  * Coley (2003) Coley A. A., 2003, Dynamical systems and cosmology. Vol. 291, Kluwer, Dordrecht, Netherlands, doi:10.1007/978-94-017-0327-7
  * Copeland et al. (1998) Copeland E. J., Liddle A. R., Wands D., 1998, Phys. Rev., D57, 4686
  * Copeland et al. (2006) Copeland E. J., Sami M., Tsujikawa S., 2006, Int. J. Mod. Phys., D15, 1753
  * Cotsakis & Kittou (2013) Cotsakis S., Kittou G., 2013, Phys. Rev., D88, 083514
  * Drepanou et al. (2021) Drepanou N., Lymperis A., Saridakis E. N., Yesmakhanova K., 2021 (arXiv:2109.09181)
  * Fan & Lu (2015) Fan Z.-Y., Lu H., 2015, Phys. Rev. D, 91, 064009
  * Ferreira & Joyce (1997) Ferreira P. G., Joyce M., 1997, Phys. Rev. Lett., 79, 4740
  * Foreman-Mackey et al. (2013) Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, PASP, 125, 306
  * Frolov & Kofman (2003) Frolov A. V., Kofman L., 2003, JCAP, 05, 009
  * García-Aspeitia & Hernández-Almada (2021) García-Aspeitia M. A., Hernández-Almada A., 2021, Phys. Dark Univ., 32, 100799
  * Garcia-Aspeitia et al. (2018) Garcia-Aspeitia M. A., Hernandez-Almada A., Magaña J., Amante M. H., Motta V., Martínez-Robles C., 2018, Phys. Rev. D, 97, 101301
  * García-Aspeitia et al. (2019) García-Aspeitia M. A., Martínez-Robles C., Hernández-Almada A., Magaña J., Motta V., 2019, Phys. Rev., D99, 123525
  * García-Aspeitia et al. (2021) García-Aspeitia M. A., Hernández-Almada A., Magaña J., Motta V., 2021, Phys. Dark Univ., 32, 100840
  * Giambo & Miritzis (2010) Giambo R., Miritzis J., 2010, Class. Quant. Grav., 27, 095003
  * Gibbons & Hawking (1977) Gibbons G. W., Hawking S. W., 1977, Phys. Rev. D, 15, 2738
  * Gim et al. (2014) Gim Y., Kim W., Yi S.-H., 2014, JHEP, 07, 002
  * Giostri et al. (2012) Giostri R., dos Santos M. V., Waga I., Reis R., Calvão M., Lago B. L., 2012, Journal of Cosmology and Astroparticle Physics, 2012, 027
  * González-Morán et al. (2021) González-Morán A. L., et al., 2021, MNRAS, 505, 1441
  * Guckenheimer & Holmes (1983) Guckenheimer J., Holmes P., 1983, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields. (Applied Mathematical Sciences Vol. 42), Springer
  * Hernández-Almada, A. et al. (2019) Hernández-Almada, A. Magaña, Juan García-Aspeitia, Miguel A. Motta, V. 2019, Eur. Phys. J. C, 79, 12
  * Hernández-Almada et al. (2022) Hernández-Almada A., Leon G., Magaña J., García-Aspeitia M. A., Motta V., Saridakis E. N., Yesmakhanova K., 2022, Monthly Notices of the Royal Astronomical Society, 511, 4147–4158
  * Herrera-Zamorano et al. (2020) Herrera-Zamorano L., Hernández-Almada A., García-Aspeitia M. A., 2020, Eur. Phys. J. C, 80, 637
  * Hurvich & Tsai (1989) Hurvich C. M., Tsai C. L., 1989, Biometrika, 76, 297
  * Jacobson (1995) Jacobson T., 1995, Phys. Rev. Lett., 75, 1260
  * Jamil et al. (2010a) Jamil M., Saridakis E. N., Setare M. R., 2010a, JCAP, 11, 032
  * Jamil et al. (2010b) Jamil M., Saridakis E. N., Setare M. R., 2010b, Phys. Rev. D, 81, 023007
  * Jeffreys (1961) Jeffreys H., 1961, Theory of Probability, third edn. Oxford, Oxford, England
  * Kaniadakis (2002) Kaniadakis G., 2002, Phys. Rev. E, 66, 056125
  * Kaniadakis (2005) Kaniadakis G., 2005, Phys. Rev. E, 72, 036108
  * Krishnan et al. (2021) Krishnan C., Ó Colgáin E., Sheikh-Jabbari M., Yang T., 2021, Phys. Rev. D, 103
  * Kunz et al. (2006) Kunz M., Trotta R., Parkinson D. R., 2006, Phys. Rev. D, 74, 023503
  * Leon et al. (2021) Leon G., Magaña J., Hernández-Almada A., García-Aspeitia M. A., Verdugo T., Motta V., 2021, JCAP, 12, 032
  * Liddle (2007) Liddle A., 2007, MNRAS: Letters, 377
  * Lymperis & Saridakis (2018) Lymperis A., Saridakis E. N., 2018, Eur. Phys. J. C, 78, 993
  * Lymperis et al. (2021) Lymperis A., Basilakos S., Saridakis E. N., 2021, Eur. Phys. J. C, 81, 1037
  * Lyra & Tsallis (1998) Lyra M. L., Tsallis C., 1998, Phys. Rev. Lett., 80, 53
  * Moradpour et al. (2020) Moradpour H., Ziaie A. H., Kord Zangeneh M., 2020, Eur. Phys. J. C, 80, 732
  * Moresco et al. (2016) Moresco M., et al., 2016, JCAP, 1605, 014
  * Motta et al. (2021) Motta V., García-Aspeitia M. A., Hernández-Almada A., Magaña J., Verdugo T., 2021, Universe, 7, 163
  * Nadathur et al. (2020) Nadathur S., Percival W. J., Beutler F., Winther H., 2020, Phys. Rev. Lett., 124, 221301
  * Padmanabhan (2005) Padmanabhan T., 2005, Phys. Rept., 406, 49
  * Padmanabhan (2010) Padmanabhan T., 2010, Rept. Prog. Phys., 73, 046901
  * Papagiannopoulos et al. (2022) Papagiannopoulos G., Basilakos S., Saridakis E. N., 2022, Dynamical system analysis of Myrzakulov gravity (arXiv:2202.10871)
  * Paranjape et al. (2006) Paranjape A., Sarkar S., Padmanabhan T., 2006, Phys. Rev. D, 74, 104015
  * Pavlović & Sossich (2021) Pavlović P., Sossich M., 2021, Phys. Rev. D, 103, 023529
  * Percival et al. (2010) Percival W. J., et al., 2010, MNRAS, 401, 2148
  * Perko (2000) Perko L., 2000, Differential Equations and Dynamical Systems, Third Edition. Springer
  * Pospelov & Pradler (2010) Pospelov M., Pradler J., 2010, Annual Review of Nuclear and Particle Science, 60, 539
  * Riess et al. (1998) Riess A. G., Filippenko A. V., Challis P., Clocchiatti A., Diercks A., et al., 1998, The Astronomical Journal, 116, 1009
  * Riess et al. (2019) Riess A. G., Casertano S., Yuan W., Macri L. M., Scolnic D., 2019, The Astrophysical Journal, 876, 85
  * Riess et al. (2021) Riess A. G., Casertano S., Yuan W., Bowers J. B., Macri L., Zinn J. C., Scolnic D., 2021, Astrophys. J. Lett., 908, L6
  * Saridakis (2009) Saridakis E. N., 2009, Nucl. Phys. B, 808, 224
  * Saridakis (2020) Saridakis E. N., 2020, JCAP, 07, 031
  * Saridakis & Basilakos (2021) Saridakis E. N., Basilakos S., 2021, Eur. Phys. J. C, 81, 7
  * Saridakis et al. (2021) Saridakis E. N., et al., 2021 (arXiv:2105.12582)
  * Schwarz (1978) Schwarz G., 1978, Ann. Statist., 6, 461
  * Scolnic et al. (2018) Scolnic D. M., et al., 2018, Astrophys. J., 859, 101
  * Sheykhi (2010a) Sheykhi A., 2010a, Eur. Phys. J. C, 69, 265
  * Sheykhi (2010b) Sheykhi A., 2010b, Phys. Rev. D, 81, 104011
  * Sheykhi (2018) Sheykhi A., 2018, Phys. Lett. B, 785, 118
  * Sheykhi et al. (2007) Sheykhi A., Wang B., Cai R.-G., 2007, Nucl. Phys. B, 779, 1
  * Spergel et al. (2003) Spergel D. N., et al., 2003, Astrophys. J. Suppl., 148, 175
  * Spiegelhalter et al. (2002) Spiegelhalter D. J., Best N. G., Carlin B. P., Van Der Linde A., 2002, Journal of the Royal Statistical Society: Series B (Statistical Methodology), 64, 583
  * Sugiura (1978) Sugiura N., 1978, Communications in Statistics - Theory and Methods, 7, 13
  * Tsallis (1988) Tsallis C., 1988, Journal of Statistical Physics, 52, 479
  * Wainwright & Ellis (1997) Wainwright J., Ellis G. F. R., 1997, Dynamical Systems in Cosmology. Cambridge University Press
  * Wang et al. (2010) Wang M., Jing J., Ding C., Chen S., 2010, Phys. Rev. D, 81, 083006
  * Weinberg (1989) Weinberg S., 1989, Rev. Mod. Phys., 61, 1
  * Zel’dovich et al. (1968) Zel’dovich Y., Krasinski A., Zeldovich Y., 1968, Sov. Phys. Usp., 11, 381
  * Zhu et al. (2021) Zhu M., Ilyas A., Zheng Y., Cai Y.-F., Saridakis E. N., 2021, JCAP, 11, 045