# Tracker behaviour of quintom dark energy and the Hubble tension

Nandan Roy<EMAIL_ADDRESS>(Corresponding Author) Centre for
Theoretical Physics & Natural Philosophy, Mahidol University, Nakhonsawan
Campus, Phayuha Khiri, Nakhonsawan 60130, Thailand  L. Arturo Ureña-López
<EMAIL_ADDRESS>Departamento de Física, DCI, Campus León, Universidad de
Guanajuato, 37150, León, Guanajuato, México

###### Abstract

We study the dynamics of the quintom dark energy model using state-of-the-art
cosmological observations. The set of equations has been converted into an
autonomous system using suitable transformations of the variables. We have
discussed the fixed points of the model and the general phase-space behavior,
in particular, in finding the existence of the tracker solutions for this
model. The observations suggest that at late times the phantom field should
dominate the dark energy sector with an approximately $15\%$ share to the
quintessence counterpart, and with both fields tracking the background at
early times. A Bayesian model comparison with $Lambda$CDM has also been done
by computing the Bayes factor and a positive preference has been obtained for
the quintom model. Although not fully resolved, the Hubble tension can be
reduced to $2.6\sigma$ when compared with the value of $H_{0}$ reported in
Riess  _et al._ (2022) and to $1.6\sigma$ when compared with that of Scolnic
_et al._ (2023).

##  I Introduction

Although the cosmological constant is the simplest and most consistent
explanation for the accelerated expansion of the universe Riess  _et al._
(1998); Perlmutter  _et al._ (1999); Aghanim  _et al._ (2018), its constant
nature gives rise to some famous problems in cosmology, such as the fine-
tuning problem and the cosmological constant problem Weinberg (1989); Peebles
and Ratra (2003); Padmanabhan (2003). Furthermore, as cosmology enters the era
of ultra-precision measurement, the constraint on some key cosmological
parameters puts the cosmological constant in more challenges. Currently, the
most significant tension in cosmology is the Hubble tension. There is a
discrepancy between the measurements of the Hubble parameter from early
universe data Aghanim  _et al._ (2018); Alam  _et al._ (2017a); Joudaki  _et
al._ (2019) and late-time data Riess  _et al._ (2016, 2019); Wong  _et al._
(2019); Freedman  _et al._ (2019) at the level of $5\sigma$. For a more in-
depth discussion of this topic, see Schöneberg  _et al._ (2022); Verde  _et
al._ (2019). One solution to the problems associated with the cosmological
constant is considering dynamical dark energy models that evolve over time
Amendola and Tsujikawa (2010); Bamba  _et al._ (2012); Copeland  _et al._
(2006); Padmanabhan (2006). The equation of state (EOS) for dark energy is
fixed at $-1$ for the cosmological constant, whereas for dynamical dark energy
models it varies over time. This has also been supported by the fact that the
current estimated value of EoS of dark energy is $w_{DE}\leq-1$ Aghanim  _et
al._ (2018).

Quite a large number of dynamical dark energy models have been proposed as
alternatives to the cosmological constant. Scalar fields with different types
of kinetic terms and potentials are often considered candidates of the
dynamical dark energy Copeland  _et al._ (2006) as, for example, quintessence,
phantom, $K$-essence, etc. Moreover, a wide variety of other models have been
suggested that can reduce or alleviate Hubble tension, such as the early dark
energy model Karwal and Kamionkowski (2016), the running vacuum model Sola
Peracaula  _et al._ (2023); Rezaei and Sola Peracaula (2022); Rezaei  _et al._
(2019, 2021), phantom-crossing models Di Valentino  _et al._ (2020); Roy  _et
al._ (2022), etc.

For quintessence scalar field models in which a canonical scalar field is
minimally coupled to gravity, the allowed value of the EOS of dark energy is
$-1\leq w_{DE}\leq 0$ Ratra and Peebles (1988); Tsujikawa (2013). A value of
EOS $w_{DE}<-1$ is not allowed for quintessence models. The consideration of a
phantom scalar field or $K$-essence models can allow us to achieve an EOS that
is less than $-1$ Nojiri  _et al._ (2005); Caldwell  _et al._ (2003); Caldwell
(2002); Roy and Bhadra (2018) but crossing the critical limit $w_{DE}=-1$
cannot also be achieved in these models. Studies have suggested that $w_{DE}$
crosses over at $-1$ Gerardi  _et al._ (2019); Johri (2004); Perković and
Štefančić (2019); Roy  _et al._ (2022) which means that at present
$w_{DE}\leq-1$ but in the past $w_{DE}>-1$.

In order to achieve a phantom barrier crossing, one possible solution is to
incorporate both quintessence and phantom scalar fields as components of the
universe. These phenomenological models are referred to as quintom models Guo
_et al._ (2005); Cai  _et al._ (2010); Zhao and Zhang (2006); Setare and
Saridakis (2008); Zhao  _et al._ (2005); Setare and Saridakis (2009). These
models have been successful in replicating the phantom divider crossing
behavior of $w_{DE}$ Zhao and Zhang (2006); Panpanich  _et al._ (2019),
suggesting that quintom models may be able to address the Hubble tension.
However, a comprehensive comparison with recent cosmological observations is
necessary to confirm this.

In this work, we have considered a general quintom model without any specific
form of the potential. Instead, we have considered the general
parameterization of the potentials proposed in Roy  _et al._ (2018) and later
used in Roy and Bamba (2019); Linares Cedeño  _et al._ (2021a); Urena-Lopez
and Roy (2020); Roy (2023) to study quintessence and phantom models. We have
constructed an autonomous system considering the appropriate transformation of
variables. As shown in our previous work for the quintessence field Urena-
Lopez and Roy (2020) we have considered the polar transformation and for the
phantom field, we have considered the hyperbolic transformation Linares Cedeño
_et al._ (2021a). Similarly to the previous work, we can identify the
existence of the tracker solutions for the quintom field with four different
scenarios.

Tracker solutions are important in alleviating scalar-field dark energy models
from coincidence and fine-tuning problems Ratra and Peebles (1988); Steinhardt
_et al._ (1999). These solutions have a remarkable property where the energy
density of the scalar field follows the background energy density, behaves
like an attractor solution for a wide range of initial conditions. We have
used recent cosmological observations to constrain the cosmological
parameters. To compare the performance of our model with the $\Lambda$CDM
model, we have used the concept of Bayesian model selection by calculating the
Bayes factor. Additionally, we have provided an update on the current status
of the Hubble tension in this model.

The paper is structured as follows. Section II deals with the mathematical
setup of the quintom model with general parameterization of the potentials.
Section III describes the critical points and the behavior of the tracker
solutions and the estimation of the initial condition. The constraint on the
cosmological parameters from the current observations is discussed in Section
IV. We summarize our findings in Section V.

##  II Mathematical background

A quintom dark energy model consists of two minimally coupled scalar fields:
one scalar field is the canonical scalar field ($\phi$), named quintessence,
and another is the phantom scalar field ($\sigma$). The action of such a
scenario is given by

$\displaystyle S=\int
d^{4}x\sqrt{-g}\bigg{[}\frac{1}{2\kappa^{2}}R+\cal{L}_{\text{DE}}+\cal{L}_{\text{M}}\bigg{]}~{},$
(1)

with the DE Lagrangian $\cal{L}_{\text{DE}}$ given explicitly by

$\begin{split}{\cal{L}_{\text{DE}}}=-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-
V_{\phi}(\phi)\\\
+\frac{1}{2}g^{\mu\nu}\partial_{\mu}\sigma\partial_{\nu}\sigma-
V_{\sigma}(\sigma)\,.\end{split}$ (2)

Here, $\kappa^{2}\equiv 8\pi G$ is the gravitational coupling and
$V_{\phi}(\phi)$, $V_{\sigma}(\sigma)$ are the respective potentials for the
canonical and phantom field, respectively. The term $\cal{L}_{\text{M}}$
corresponds to the standard matter content of the universe. We consider all
components of the universe to be connected by a barotropic equation of state
(EoS) $\rho_{i}=w_{i}p_{i}$. For example, $w_{i}=1/3$ for radiation and
$w_{i}=0$ for matter.

###  II.1 Quintom equations of motion

Once we consider the universe to be flat, homogeneous, and isotropic, the
Friedmann equations are written as follows,

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{\kappa^{2}}{3}\Big{(}\sum_{i}\rho_{i}+\rho_{\phi}+\rho_{\sigma}\Big{)}\,,$
(3a) $\displaystyle\dot{H}$ $\displaystyle=$
$\displaystyle-\frac{\kappa^{2}}{2}\Big{[}\sum_{i}(\rho_{i}+p_{i})+\rho_{\phi}+p_{\phi}+\rho_{\sigma}+p_{\sigma}\Big{]}\,,$
(3b)

where $H=\dot{a}/a$ is the Hubble parameter and $a$ is the scale factor of the
universe. The respective Klein-Gordon equations of the two constituents of the
quintom model are:

$\displaystyle\ddot{\phi}+3H\dot{\phi}+\frac{\partial
V_{\phi}(\phi)}{\partial\phi}$ $\displaystyle=$ $\displaystyle 0\,,$ (4a)
$\displaystyle\ddot{\sigma}+3H\dot{\sigma}-\frac{\partial
V_{\sigma}(\sigma)}{\partial\sigma}$ $\displaystyle=$ $\displaystyle 0\,.$
(4b)

In the above expressions, $p_{\phi}$ and $\rho_{\phi}$ are, respectively, the
pressure and density of the canonical scalar field, and similarly $p_{\sigma}$
and $\rho_{\sigma}$ are the pressure and density of the phantom field. Their
mathematical expressions are given by:

$\displaystyle\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+V_{\phi}(\phi)\,,\quad
p_{\phi}=\frac{1}{2}\dot{\phi}^{2}-V_{\phi}(\phi)\,,$ (5a)
$\displaystyle\rho_{\sigma}=-\frac{1}{2}\dot{\sigma}^{2}+V_{\sigma}(\sigma)\,,\quad
p_{\sigma}=-\frac{1}{2}\dot{\sigma}^{2}-V_{\sigma}(\sigma)\,.$ (5b)

Since in a quintom model, the DE sector is the combination of the canonical
and phantom fields and there is no interaction considered, its total density
and pressure can be written in the following way;

$\rho_{DE}\equiv\rho_{\phi}+\rho_{\sigma}\,,\quad p_{DE}\equiv
p_{\phi}+p_{\sigma}\,,$ (6a) while its total equation of state is given by
$w_{DE}\equiv\frac{p_{DE}}{\rho_{DE}}=\frac{p_{\phi}+p_{\sigma}}{\rho_{\phi}+\rho_{\sigma}}\,.$
(6b)

Also, using the conservation equations, the total density can be written as
$\rho_{tot}\equiv\sum_{i}\rho_{i}+\rho_{\phi}+\rho_{\sigma}$, obtaining:

$\dot{\rho}_{tot}+3H(1+w_{tot})\rho_{tot}=0\,,$ (7)

with the corresponding total EoS,

$w_{tot}=\frac{p_{\phi}+p_{\sigma}+\sum_{i}p_{i}}{\rho_{\phi}+\rho_{\sigma}+\sum_{i}\rho_{i}}=w_{\phi}\Omega_{\phi}+w_{\sigma}\Omega_{\sigma}+\sum_{i}w_{i}\Omega_{i}\,.$
(8)

The individual EoS parameters are defined as

$\displaystyle w_{\phi}=\frac{p_{\phi}}{\rho_{\phi}}\,,\quad
w_{\sigma}=\frac{p_{\sigma}}{\rho_{\sigma}}\,,\quad
w_{i}=\frac{p_{i}}{\rho_{i}}\,,$ (9)

together with the corresponding density parameters,

$\displaystyle\Omega_{\phi}\equiv\frac{\rho_{\phi}}{\rho_{tot}}~{},~{}~{}\Omega_{\sigma}\equiv\frac{\rho_{\sigma}}{\rho_{tot}}~{},~{}~{}\Omega_{i}\equiv\frac{\rho_{i}}{\rho_{tot}}~{}.$
(10)

Notice that the total DE density parameter is given by the combination
$\Omega_{\phi}+\Omega_{\sigma}\equiv\Omega_{DE}$, whereas the combination of
all density parameters results in the Friedmann constraint from Eq. (3a):
$\Omega_{\phi}+\Omega_{\sigma}+\Omega_{m}+\Omega_{r}=1$.

###  II.2 Autonomous system of equations

To write the system (4) as a set of autonomous equations, we define the
following set of dimensionless polar variables,111For references on these
changes of variables of scalar fields see for instance Copeland  _et al._
(1998); Roy and Bhadra (2018); Roy  _et al._ (2018); Roy and Bamba (2019).

$\displaystyle\frac{\kappa\dot{\phi}}{\sqrt{6}H}=\Omega_{\phi}^{1/2}\sin(\theta_{\phi}/2)\,,\quad\frac{\kappa
V_{\phi}^{1/2}}{\sqrt{3}H}=\Omega_{\phi}^{1/2}\cos(\theta_{\phi}/2)\,,$ (11a)
$\displaystyle
y_{1\phi}\equiv-2\sqrt{2}\frac{\partial_{\phi}V_{\phi}^{1/2}}{H}\,,\quad
y_{2\phi}\equiv-4\sqrt{3}\frac{\partial^{2}_{\phi}V_{\phi}^{1/2}}{\kappa
H}\,,$ (11b)

and hyperbolic ones,

$\displaystyle\frac{\kappa\dot{\sigma}}{\sqrt{6}H}=$
$\displaystyle\Omega_{\sigma}^{1/2}\sinh(\theta_{\sigma}/2)\,,\quad\frac{\kappa
V_{\sigma}^{1/2}}{\sqrt{3}H}\,=\Omega_{\sigma}^{1/2}\cosh(\theta_{\sigma}/2)\,,$
(12a) $\displaystyle y_{1\sigma}\equiv$
$\displaystyle-2\sqrt{2}\frac{\partial_{\phi}V_{\sigma}^{1/2}}{H}\,,\quad
y_{2\sigma}\equiv-4\sqrt{3}\frac{\partial^{2}_{\phi}V_{\sigma}^{1/2}}{\kappa
H}\,,$ (12b)

where the variables with suffix $\phi$ represent the variables related to the
canonical field and with $\sigma$ represent the phantom field. As a result,
the field equations (4) are rewritten as the following set of equations,

$\displaystyle\theta_{\phi}^{\prime}$ $\displaystyle=$
$\displaystyle-3\sin\theta_{\phi}+y_{1\phi}\,,$ (13a) $\displaystyle
y_{1\phi}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\gamma_{tot}y_{1\phi}+\Omega_{\phi}^{1/2}\sin(\theta/2)y_{2\phi}\,,$
(13b) $\displaystyle\Omega_{\phi}^{\prime}$ $\displaystyle=$ $\displaystyle
3(\gamma_{tot}-\gamma_{\phi})\Omega_{\phi}\,,$ (13c)
$\displaystyle\theta_{\sigma}^{\prime}$ $\displaystyle=$
$\displaystyle-3\sinh\theta_{\sigma}-y_{1\sigma}\,,$ (13d) $\displaystyle
y_{1\sigma}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{3}{2}\gamma_{tot}y_{1\sigma}+\Omega_{\sigma}^{1/2}\sinh(\theta_{\sigma}/2)y_{2\sigma}\,,$
(13e) $\displaystyle\Omega_{\sigma}^{\prime}$ $\displaystyle=$ $\displaystyle
3(\gamma_{tot}-\gamma_{\sigma})\Omega_{\sigma}\,.$ (13f)

From Eqs. (9), the barotropic EOS of the canonical scalar field can be written
in terms of the polar variable as
$\gamma_{\phi}=1+w_{\phi}=1-\cos\theta_{\phi}$, and then $0<\gamma_{\phi}<2$,
while that of the phantom field is
$\gamma_{\sigma}=1+w_{\sigma}=1-\cosh\theta_{\sigma}$, from which
$\gamma_{\sigma}<0$. The barotropic EOS parameter of the total dark energy,
see Eq. (6b), is given by

$\gamma_{DE}=\frac{\gamma_{\phi}\Omega_{\phi}+\gamma_{\sigma}\Omega_{\sigma}}{\Omega_{\phi}+\Omega_{\sigma}}\,.$
(14)

One can see that the system of equations (13) is not closed until one
specifies the functional form of $y_{2\phi}$ and $y_{2\sigma}$. These two
variables, by definition, carry information about quintessence and phantom
potentials. In general, there are two different ways to consider the
functional forms of the variables $y_{2}$. First, one can consider a
particular form of potential and then calculate the corresponding $y_{2}$, and
secondly, one assumes a particular form of $y_{2}$ in terms of variables $y$
and $y_{1}$, and integrate back to get the actual form of the potential.

In this work, we use the parameterization of $y_{2}$ proposed in Roy  _et al._
(2018), which includes a large class of dark energy potentials (see Table 1
and 2 in Roy  _et al._ (2018)). Thus, we have for the quintom potential,

$\displaystyle y_{2\phi}$ $\displaystyle=$
$\displaystyle\alpha_{0}y_{\phi}+\alpha_{1}y_{1\phi}+\alpha_{2}y^{2}_{1\phi}/y_{\phi}\,,$
(15a) $\displaystyle y_{2\sigma}$ $\displaystyle=$
$\displaystyle\beta_{0}y_{\sigma}+\beta_{1}y_{1\sigma}+\beta_{2}y^{2}_{1\sigma}/y_{\sigma}\,,$
(15b)

where $\alpha$ and $\beta$ are just constant parameters.

As explained in Roy  _et al._ (2018), the parameters $\alpha$ and $\beta$ are
called active parameters of the potentials, as they are the only ones that
appear explicitly in the equations of motion. There are other parameters in
the potentials, but these are instead called passive parameters, and their
values can be connected to the initial values of the dynamical variables
$\theta_{\phi,\sigma}$, $y_{1\phi,\sigma}$ and $\Omega_{\phi,\sigma}$.

Though the above-mentioned parametrization covers a wide range of potentials,
there are still some potentials which cannot be written in the above form. For
example, $V(\phi)=\sum_{i}^{n}\Lambda_{i}\exp(\lambda_{i}\phi)$ Copeland  _et
al._ (1999); Ivashchuk  _et al._ (2003) and
$V(\phi)=V_{0}\tanh^{2}(\frac{\lambda_{1}\phi}{m_{p}})\cosh(\frac{\lambda_{2}\phi}{m_{p}})$
Bag  _et al._ (2018) are not included in the list. However, examining the
parameterizations given in (15) allows us to include a wide range of
potentials in our system, which are frequently used in the relevant
literature.

##  III Critical points with tracking behavior

Here we discuss the critical points of the quintom model that are of physical
interest, particularly with tracking behavior. This study will be useful for
understanding the combined dynamics of the quintom field and our particular
choice of initial conditions.

###  III.1 Critical points and tracker conditions

Critical points are obtained by simultaneously solving the equations in (13)
when all derivatives of the variables on the left side are zero. From
equations (13a) and (13d), one can see from
$\theta_{\phi}^{\prime}=\theta_{\sigma}^{\prime}=0$, that the critical
conditions $y_{1\phi c}=3\sin{\theta_{\phi_{c}}}$ and $y_{1\sigma
c}=-3\sinh{\theta_{\sigma_{c}}}$ are common to all fixed points. Using the
forms of $y_{2\phi}$ and $y_{2\sigma}$ given in equations (15), one can derive
from $y_{1\phi}^{\prime}=y_{1\sigma}^{\prime}=0$ the following two equations;

$\displaystyle\left(\gamma_{tot}+\frac{\alpha_{0}}{9}\Omega_{\phi c}\right.$
$\displaystyle\left.\quad+\frac{2}{3}\alpha_{1}\Omega_{\phi
c}^{1/2}\sin(\theta_{c}/2)+4\alpha_{2}\sin^{2}(\theta_{c}/2)\right)\sin\theta_{c}=0\,,$
(16)

and

$\displaystyle\left(\gamma_{tot}-\frac{\beta_{0}}{9}\Omega_{\sigma c}\right.$
$\displaystyle\left.\quad+\frac{2}{3}\beta_{1}\Omega_{\sigma
c}^{1/2}\sinh\left(\theta_{c}/2\right)-4\beta_{2}\sinh^{2}\left(\theta_{c}/2\right)\right)\sinh\theta_{\sigma
c}=0\,.$ (17)

In general, critical points can be classified into three different classes
depending on the dominant energy densities. The fluid dominated fixed points
($\Omega_{\phi c}=\Omega_{\sigma c}=0$), the quintessence ($\Omega_{\phi
c}=1$), and the phantom field ($\Omega_{\phi c}=1$) dominated fixed points. In
Appendix A.1 we have listed all possible fixed points of the system and
discussed the stability of these fixed points. Here, we are going to discuss
in detail the fluid-dominated fixed points as they correspond to the tracking
behavior of the scalar fields.

It has been shown in Urena-Lopez and Roy (2020); Linares Cedeño  _et al._
(2021a) that Eqs. (III.1) and (III.1), for the quintessence and phantom
models, respectively, can give us a generalized tracking behavior for any
choice of the parameters $\alpha$ and $\beta$ once we consider
$\Omega_{\phi}\simeq\Omega_{\sigma}\simeq 0$ in the early universe. However,
the equivalent simplest tracker condition can be calculated considering
$\alpha_{0}=\alpha_{1}=0$, $\beta_{0}=\beta_{1}=0$ and using them in equations
(III.1) and (III.1). For the quintessence field and the phantom field, the
tracker condition will be, respectively,

$\gamma_{\phi c}=-\gamma_{tot}/(2\alpha_{2})\,,\quad\gamma_{\sigma
c}=-\gamma_{tot}/(2\beta_{2})\,.$ (18)

For the case of the quintessence field, an upper bound can be placed on the
parameter $\alpha_{2}$ by considering the fact $0<\gamma_{\phi}<2$ and this
will translate into $\alpha_{2}<-\gamma_{tot}/4$. For the tracking to start at
the early radiation domination $\gamma_{tot}=4/3$ and $\alpha_{2}<-1/3$. For
the phantom field, the condition would simply be $\beta_{2}>0$.

###  III.2 Fluid-dominated critical points

In Table 1 we have listed all fixed points dominated by background fluids
$m_{1}$ to $m_{4}$. The fluid-dominated fixed points are characterized by
$\Omega_{\phi c}=\Omega_{\sigma c}=0$, which means that the universe is
dominated by the background fluids in their corresponding eras, such as
radiation domination, matter domination, and the contribution to the total
energy density of the quintessence and phantom fields is negligible. During
radiation domination $\gamma_{tot}=4/3$ and matter domination
$\gamma_{tot}=1$. One can see from Table 1 that there are four matter-
dominated fixed points corresponding to four physically interesting scenarios
depending on the tracking behavior of the scalar fields.

Table 1:  List of fluid dominated fixed points. For fixed point $m_{1}$, no tracker behaviour of the scalar field is obtained, for $m_{2}$ (phantom) and and $m_{3}$ (quintessence) only one field shows tracker behaviour. For $m_{4}$ fixed point both the scalar fields shows tracker behaviour. | Fixed
---
Points
$\Omega_{\phi_{c}}$ | $\sin\theta_{\phi_{c}}$ | $\Omega_{\sigma_{c}}$ | $\sinh\theta_{\sigma_{c}}$
$m_{1}$ | 0 | $n\pi$ | 0 | 0
$m_{2}$ | 0 | $n\pi$ | 0 | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{\gamma_{tot}}{4\beta_{2}}$
$m_{3}$ | 0 | $\sin^{2}(\theta_{\phi_{c}}/2)=-\frac{\gamma_{tot}}{4\alpha_{2}}$ | 0 | 0
$m_{4}$ | 0 | $\sin^{2}(\theta_{\phi_{c}}/2)=-\frac{\gamma_{tot}}{4\alpha_{2}}$ | 0 | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{\gamma_{tot}}{4\beta_{2}}$

The fixed point $m_{1}$ corresponds to a scenario where the quintessence and
phantom fields behave like a cosmological constant and there is no tracking
phenomenon from any fluid. In Fig. 1a we have plotted the EoS of the total
dark energy and its components for the solutions that originate from the fixed
point $m_{1}$ for various choices of the parameters $\alpha$ and $\beta$.
Here, we make sure that $\alpha_{2}>-1/3$ and $\beta_{2}<0$ to turn off the
tracker behavior. As expected from the analytical point of view, there is no
tracker behavior and the total dark energy EoS $w_{DE}$ remains
indistinguishable from the cosmological constant, but at a late time it starts
to evolve and currently $w_{DE}>-1$.

(a)  (b)  (c)  (d)

Figure 1:  The following plots show the equations of state (EoS) for the total
dark energy and its component scalar fields (quintessence and phantom) under
different choices of $(\alpha,\beta)$ parameters. In (a), the solution that
originated at the fixed point $m_{1}$ is displayed where the scalar field does
not show the tracking behavior. In (b), the solutions from the fixed point
$m_{2}$ are shown where only the phantom field shows tracking behavior. In
(c), the solutions from the fixed point $m_{3}$ are shown where only the
quintessence field tracks. Finally, in (d), the solutions from the $m_{4}$
fixed point are shown, in which both fields track the background.

At $m_{2}$ the dominating fluid is tracked by the phantom field, but the
quintessence field remains indistinguishable from the cosmological constant in
both the radiation and matter dominated eras. In Fig. 1b, we show that the
solutions originate at the critical point $m_{2}$ for different parameters
$\alpha,\beta$ with the choice of $\alpha_{2}>-1/3$ and $\beta_{2}>0$. The
expected behavior described above in the critical point analysis can be seen
in the evolution of the numerical solutions. For this case, currently
$w_{DE}>-1$

Similarly, the $m_{3}$ quintessence field tracks the fluid counterpart, while
the phantom field remains indistinguishable from the cosmological constant
during radiation and matter domination, but it shows deviation at the late
time. In Fig. 1c, we have plotted the dynamics of the EoS of the quintom field
together with the components for different values of the parameters
$\alpha,\beta$, with $\alpha_{2}<-1/3$ and $\beta_{2}<0$, so that the
quintessence field satisfies the tracker condition while the phantom field
does not.

The most interesting case is the solutions originating from the fixed point
$m_{4}$, in which both the quintessence and the phantom fields track the
background by satisfying the conditions $\alpha_{2}<-1/3$ and $\beta_{2}>0$.
In Fig. 1d, we have plotted the EoS of both the scalar fields and the total
dark energy. It can be seen that both quintessence and phantom fields track
the background, and the corresponding EoS is decided by the tracker condition
they follow. The total dark energy EoS shows interesting behavior; deep in the
radiation-dominated era, it tracks the evolution of the quintessence field
followed by a phantom barrier crossing and the current value of $w_{DE}<-1$.

Figure 2:  Plots of EoS of the quintessence field, the phantom field and the
total DE EoS for different choices of the $f$ parameter when both the fields
are tracking the background. See the text for more details.

It is interesting to note here that the fixed points $m_{1}$ to $m_{3}$ are
reported in the previous section while discussing the fixed points of
quintessence and phantom domination separately and also in Urena-Lopez and Roy
(2020); Linares Cedeño  _et al._ (2021a). The fixed point $m_{4}$ appears when
both the phantom and quintom fields are switched on in the evolution of the
universe. The choice of parameters $\alpha_{2}$ and $\beta_{2}$, for all the
figures mentioned in this section, is also within the posterior range obtained
from the statistical analysis in the next section.

###  III.3 Initial conditions

Here, we discuss the method for finding the initial condition of the system by
solving Eqs. (13) using the tracker condition. An estimation of the proper
initial conditions is necessary for the implementation of the model in the
Boltzmann code called the Cosmic Linear Anisotropy Solving System (CLASS)
Lesgourgues (2011a, b, c); Lesgourgues and Tram (2011).

During the era of matter and radiation domination, both quintessence and the
phantom field contribute negligibly to the energy density, and the total EoS
remains constant at $\gamma_{tot}=4/3$ for radiation and $\gamma_{tot}=1$ for
matter domination. We can use this fact and the tracking condition for the
quintessence and phantom fields (see Eqs. (18)) to estimate the initial
conditions for which a viable evolution can be achieved. We have used the
scheme mentioned above to approximate the solution of the equations in Eq.
(13) at the radiation and matter domination, and match them at the radiation
matter equality to estimate the initial conditions. The initial condition for
the quintessence field would be the following;

$\displaystyle\cos\theta_{i}$ $\displaystyle=1+\frac{2}{3\alpha_{2}}\quad
y_{1i}=3\sin\theta_{i}$ (19a) $\displaystyle\Omega_{\phi i}$
$\displaystyle=a_{i}^{4\left(1+1/2\alpha_{2}\right)}\left(\frac{\Omega_{m0}}{\Omega_{r0}}\right)^{1+1/2\alpha_{2}}\Omega_{\phi_{0}}$
For the phantom field, the initial conditions are,
$\displaystyle\cosh\theta_{i}$ $\displaystyle=1+\frac{2}{3\beta_{2}},\quad
y_{1i}=-3\sinh\theta_{i}$ (19b) $\displaystyle\Omega_{\phi i}$
$\displaystyle=a_{i}^{4\left(1+1/2\beta_{2}\right)}\left(\frac{\Omega_{m0}}{\Omega_{r0}}\right)^{1+1/2\beta_{2}}\Omega_{\sigma_{0}}$

where $\Omega_{\phi 0}=A_{\phi}(1-f)\Omega_{DE0}$ and $\Omega_{\sigma
0}=A_{\sigma}f\Omega_{DE0}$. In addition, $f$ is an extra sampling parameter
which we introduce to have better control over the contributions of
$\Omega_{\phi i}$ and $\Omega_{\sigma i}$ to the total density.

The background equations of motion (13), together with the initial conditions
(19), were incorporated into an amended version of the code CLASS following
the prescription already described in previous work, for example Urena-Lopez
and Roy (2020); Linares Cedeño  _et al._ (2021a); Roy  _et al._ (2018). The
present values of the matter density and radiation densities are represented
by $\Omega_{r0}$ and $\Omega_{m0}$, respectively.

One difference with respect to single-scalar field models is that we need to
adjust the initial conditions of the quintom fields, so that the solution of
the equations of motion delivers the right present value of the total DE
density parameter; i.e. $\Omega_{DE0}$, as the sum of the present value of the
quintessence energy density $\Omega_{\phi 0}$ and the phantom energy density
$\Omega_{\sigma 0}$. For that, we let the internal algorithm of CLASS find the
appropriate value of the coefficients $A_{\phi}$ and $A_{\sigma}$ for a given
choice of the parameter $f$. The initial value of the scale factor is
typically considered $a_{i}\simeq 10^{-14}$ in the CLASS code.

In Figure 2, we have displayed the development of the EOS of both the scalar
field and the total dark energy EOS for different values of the $f$ parameter.
It is evident that the individual EOS of the fields is not affected by
different choices of $f$, however, the total dark energy EOS is impacted and
the epoch of phantom barrier crossing is sensitive to it. This is because the
definition of the $w_{tot}$ (see Eq. (8)) involves the $\Omega_{\phi}$ and
$\Omega_{\sigma}$. A lower value of $f$, which corresponds to $\Omega_{\phi
i}>>\Omega_{\sigma i}$, causes the crossing of the phantom barrier to be
shifted further towards the matter-dominated era.

##  IV Constraints on cosmological parameters

We have used current cosmological observations to constrain the quintom model
presented here by using an amended version of the Boltzmann code class and the
cosmological MCMC sampler MontePython (v3.5) Audren  _et al._ (2013). The
cosmological observations that have been used for this are PantheonPlus Brout
_et al._ (2022) together with SH0ES priorRiess  _et al._ (2022), BAO (BOSS
DR12 Alam  _et al._ (2017b), eBOSS DR14 (Lya) Cuceu  _et al._ (2019)) and the
WiggleZ galaxy survey Kazin  _et al._ (2014). We have also used a compressed
Planck likelihood by following the proposal in Arendse  _et al._ (2020) (see
their Appendix A). The parameter constraining accuracy of compressed Planck
likelihood is the same as the full Planck data, but consumes less time and
computational resources.

We have considered flat priors for the cosmological parameters
$100~{}\omega_{b}:[1.9,2.5],\omega_{cdm}:[0.095,0.145]$. Following the
suggestion in Sabti  _et al._ (2022), we fix the angular scale of the sound
horizon $\theta_{s}$ to the Planck CMB value of $\theta_{s}=1.04110$ Aghanim
_et al._ (2018), and the current value of the Hubble parameter $H_{0}$ as the
derived parameter. Since $\theta_{s}$ is determined by the angular scales of
the acoustic peaks in the temperature and polarization spectra as a purely
geometrical quantity, and its determination is mostly independent of the
underlying physics considered for the CMB era, fixing it does not make any
significant changes to the results.

A flat prior of $[-20,20]$ has been considered for the parameters $\alpha_{2}$
and $\beta_{2}$. Since it has already been shown in Roy  _et al._ (2018); Roy
and Bamba (2019); Ureña-López (2019); Linares Cedeño  _et al._ (2021b) that
the parameters $\alpha$ and $\beta$ remain unconstrained for quintessence and
phantom, we consider $\alpha_{0}=\alpha_{1}=0$ and $\beta_{0}=\beta_{1}=0$.
Parameters $\alpha_{2}$ and $\beta_{2}$ are the most interesting, as they
determine the tracking behavior of the scalar fields. One can also see from
our choice of priors that it includes all four tracking scenarios described in
the previous section.

The constraints obtained from the observations are given in Table 2, with
their mean value and the 68%CL. For parameter estimation, we have used the
publicly available MCMC analyzer code GetDistLewis (2019).

Table 2:  68%CL constraints on the cosmological parameters for the $\Lambda$CDM and the quintom models using the combining data sets of Pantheon Plus Survey Brout  _et al._ (2022) with SH0ES priorRiess  _et al._ (2022), Baryon Acoustic Oscillations from the BOSS DR12 survey Alam  _et al._ (2017b), Lyman-alpha forest data from the eBOSS DR14 survey Cuceu  _et al._ (2019), and the WiggleZ galaxy survey Kazin  _et al._ (2014), along with the compressed Planck likelihood. Parameter |  $\Lambda$CDM | Quintom
---|---|---
$10^{-2}\omega{}_{b}$ | $2.260\pm 0.015$ | $2.250^{+0.015}_{-0.013}$
$\omega_{cdm}$ | $0.11719^{+0.00086}_{-0.00077}$ | $0.11832^{+0.00091}_{-0.0011}$
$H_{0}$ | $68.35^{+0.33}_{-0.37}$ | $69.56^{+0.21}_{-0.96}$
$\Omega_{DE}$ | $0.6992\pm 0.0064$ | $0.7070^{+0.0040}_{-0.0082}$
$\Omega_{\phi}$ | – | $0.0966^{-0.0055}_{-0.14}$
$\Omega_{\sigma}$ | – | $0.61^{+0.12}_{-0.39}$
$w_{DE}$ | -1 | $-1.047^{+0.035}_{-0.0098}$
$w_{\phi}$ | – | $-0.9536^{+0.0061}_{-0.064}$
$w_{\sigma}$ | – | $-1.058^{+0.044}_{-0.0076}$
${\alpha_{2}}$ | – | $-2^{+12}_{-17}$
${\beta_{2}}$ | – | $9.1^{+3.5}_{-5.8}$
$\delta{\chi_{min}^{2}}$ | 0 | $-6.0$
$\ln{B_{Q\Lambda}}$ | - | $1.21$

In Fig. 3, we show the contour plot of the posterior of Hubble parameter
$H_{0}$ and $\Omega_{m}$, with the quintom model in red and the $\Lambda$CDM
model in blue. The green band shows the $1\sigma$ (deep green) and $2\sigma$
(light green) constraints on $H_{0}$ from Riess  _et al._ (2022), whereas the
region covered by the black dashed lines shows the $1\sigma$ and $2\sigma$
constraints obtained from Scolnic  _et al._ (2023). There is a significant
shift of the current value of the Hubble parameter towards the higher value.
We use the following estimator of the tension given in Camarena and Marra
(2018) to quantify the $H_{0}$ tension;

$T_{H0}=\frac{\left|H_{0}-H_{0}^{loc}\right|}{\sqrt{\sigma_{H0}^{2}+\sigma_{\mathrm{loc}}^{2}}},$

where $T_{H0}$ is the amount of tension and $H_{0}$ is the mean and
$\sigma_{H0}^{2}$ is the variance of the posterior $p(H_{0})$.

Once we compare our results with the current value of the Hubble parameter
($H_{0}=73.04\pm 1.04\mathrm{~{}km}\mathrm{~{}s}^{-1}\mathrm{Mpc}^{-1}$)
reported in Riess  _et al._ (2022), from the Hubble Space Telescope and the
SH0ES data, the tension decreases to $T_{H0}\simeq 2.85\sigma$. Compared to
the Hubble value ($H_{0}=73.22\pm
2.06\mathrm{~{}km}\mathrm{~{}s}^{-1}\mathrm{Mpc}^{-1}$) reported in Scolnic
_et al._ (2023), from the standardized TRGB and Type Ia supernova data sets,
the tension reduces to $T_{H0}\simeq 1.67\sigma$. According to Table IV in
Camarena and Marra (2018) for the first case, the tension is moderate in
nature and for the second case the tension is weak.

Figure 3:  Two dimensional posterior plot of the $H_{0}$ versus $\Omega_{m}$
parameters. The horizontal gray band shows the constraint on the $H_{0}$
parameter from the SH0ES collaboration Riess  _et al._ (2019). See the text
for more details.

The triangular plots in Fig. 4 show the 1D and 2D posteriors of the
cosmological parameters
$H_{0},10^{-2}w_{b},w_{cdm},\Omega_{DE},\Omega_{\phi},\Omega_{\sigma}$. A
comparison with the standard $\Lambda$CDM case has been shown by plotting the
constraints on the parameters in blue contours. In Fig. 5, we have shown the
posterior distribution of the total dark energy EoS $w_{DE}$, the quintessence
EoS $w_{\phi}$, the phantom EoS $w_{\sigma}$, and of the parameter $f$.

The value of $w_{DE}$ remains in the phantom region (see Table 2). An
interesting fact to be observed here from Fig. 6 is that $\alpha_{2}$ remains
unconstrained, whereas $\beta_{2}\simeq 9.1^{+9.9}_{-7.8}\,(2\sigma)$
indicates the tracking behavior of the phantom field. The unconstrained
$\alpha_{2}$ can be explained by the fact that the current observations
require the EoS of the dark energy to be $w_{DE}<-1$ and the quintessence
field should not dominate the late-time dynamics.

Though the phantom field dominates late-time cosmology, the contribution of
the quintessence field to the current budget is non-negligible, as it is
approximately $15\%$ of the total dark energy budget. There could be two
physical scenarios related to the unconstrained nature of $\alpha_{2}$. The
quintessence field tracks the background during radiation domination, as
mentioned in Sec. III.2 and Fig. 1d. This case will always be associated with
a phantom barrier crossing. In the other scenario, the quintessence field does
not track the background, the total $w_{DE}$ remains indistinguishable from
the cosmological constant case $w_{DE}=-1$, and at the late time becomes
phantom in nature (Sec. III.2 and Fig. 1b).

As can be seen by comparing Figs. 1c and 1b, when the phantom field tracks the
background but the quintessence field does not, the current value of the EoS
satisfies $w_{DE}>-1$. Since observational constraints suggest the current
value of $w_{DE}<-1$, it is more reasonable for both fields to track the
background, even though observations cannot constrain $\alpha_{2}$. Due to the
lower bound of $\beta_{2}>1.3$, it can be checked that the epoch of the
crossing of the phantom barrier should occur after the radiation-matter
equality. In Fig. 6, we have also shown the posterior distribution of the
parameter $f$, which is related to the distribution of energy between the
phantom and the quintessence at the initial stage. Although the dynamics of
the scalar fields are sensitive to the $f$ parameter, $f$ remains
unconstrained by observations. The possible reason behind this might be that
for any particular choice of the $f$ there are some values of
$\alpha_{2},\beta_{2}$ that help the quintom model to satisfy the
observational constraints.

We have used $\chi^{2}_{min}$ (see Table 2) to compare the $\Lambda$CDM model
with the quintom model. The difference in value between the models,
$\Delta\chi^{2}_{min}=\chi^{2}_{Quintom}-\chi^{2}_{\Lambda CDM}=-6$, suggests
a better fit to the data by the quintom model compared to $\Lambda$CDM. We
also computed the Bayes factor to be certain of the preference for the model
chosen by the data, so that the presence of extra model parameters for the
quintom model does not affect it. We have calculated $\ln
B_{Q\Lambda}=\ln\mathcal{Z}_{Q}-\ln\mathcal{Z}_{\Lambda}$, where $\mathcal{Z}$
is Bayesian evidence, and the suffix $Q$ and $\Lambda$ represent the quintom
and $\Lambda$CDM models, respectively.

In general, Jeffrey’s scale is used to find the preference of one model over
another. The strength of preference of the quintom model over $\Lambda$CDM
will depend on the value of $\ln B_{Q\Lambda}$. The preference is negative if
$\ln B_{Q\Lambda}<1$, and if $\ln B_{Q\Lambda}>1;>2.5;>5.0$ the preference is
positive, moderate, and strong, respectively. For more details on Bayesian
model selection, see Trotta (2007). Since the calculation of the Bayes factor
is computationally challenging, we have used the publicly available code
MCEvidence Heavens  _et al._ (2017). This code can be used directly to
calculate the Bayes factor from the MCMC chains generated by MontePython. In
the comparison of Quintom with $\Lambda$ CDM, we find $\ln B_{Q\Lambda}=1.21$.
This shows that the quintom model would be positively favored by the
observations compared to $\Lambda$CDM.

Figure 4:  Triangular plot displays 1D and 2D posteriors of the cosmological
parameters using Pantheon Scolnic  _et al._ (2018), BAO (DR12 Alam  _et al._
(2017b), DR14 Cuceu  _et al._ (2019)), WiggleZ Kazin  _et al._ (2014)
considering SH0ES and compressed Planck prior. The blue contours and curves
are for $\Lambda$CDM and the red ones are for the quintom model. Figure 5:
Triangle plot showing 1D and 2D posteriors of the EoS of the total dark
energy($w_{DE}$) and its components $w_{\phi}$ and $w_{\sigma}$.  Figure 6:
Triangle plot showing 1D and 2D posteriors of the two tracking parameters
$\alpha_{2}$ for the quintessence and $\beta_{2}$ for the phantom scalar
field. Notice the unconstrained nature of parameters $f$ and $\alpha_{2}$, and
the existence of a preferred value for $\beta_{2}$. See the text for more
details.

##  V Conclusions

In this work, we have studied the dynamics of the quintom scalar field model
and used the current cosmological observations to constrain the cosmological
parameters. Quintom dark energy model consists of two scalar fields; one is
the canonical scalar field called quintessence and the other is the phantom
scalar field. For both of these scalar fields, the choice of potentials
remains arbitrary, since there is no observational preference for any
particular form of the potential Roy  _et al._ (2018); Linares Cedeño  _et
al._ (2021a); Roy and Bamba (2019), and we have then considered the general
parameterization of potentials proposed in Roy  _et al._ (2018). The set of
equations has been converted into an autonomous system using suitable
transformations of the variables. For the quintessence sector, a polar
transformation has been considered, while a hyperbolic one was required to
incorporate the negative kinetic energy of the phantom field.

We have discussed the fixed points of the model and the general phase-space
behavior, in particular, in finding the existence of the tracker solutions for
the quintom model. There could be four physically interesting scenarios: no
tracking behavior of either field, either one field tracks the background
while another does not, or both fields track the background. We have also
found the tracking conditions for quintessence and phantom fields separately,
which has previously been reported in the literature. Of the four cases, the
most interesting is the case in which both fields track the background, with
the quintessence field contributing more to the dark energy density at early
times, whereas the phanton field dominates the universe currently. If both
fields track the background, there should be a phantom barrier crossing of the
total dark energy EoS which is most likely to happen after the matter-
radiation equality era.

The available cosmological observations have been used to constrain the
cosmological parameters. Our choice of priors includes both the tracking and
no-tracking behavior of the quintom and the phantom field. The observations
suggest that at late times the phantom field should dominate the dark energy
sector with an approximately $15\%$ share to the quintessence counterpart, and
with both fields tracking the background at early times.

We have also reported the current status of the Hubble tension in the quintom
model. Although it is not possible to completely resolve the Hubble tension in
this model, it can be reduced to the level of $2.6\sigma$ compared to the
Hubble value reported in Riess  _et al._ (2022), and to $1.6\sigma$ when
compared with Scolnic  _et al._ (2023), making it a moderate and week tension,
respectively, comparable to the strong tension of $\Lambda$ CDM. We have also
studied the selection of Bayesian models by computing the Bayes factor to
compare with $\Lambda$CDM and found a positive preference for the quintom
model.

For completeness, a thorough dynamical system analysis of the model is
included in the Appendix. It is possible that the late-time attractors are
either quintessence- or phantom-field dominated. Since current observations
indicate that the universe is currently phantom dominated, the phantom
attractor points are more pertinent. It would be beneficial to incorporate new
data sets and determine whether the preference for phantom domination
continues, as well as to gain a better understanding of the active parameters
of our model. This research will be conducted in the future and reported
elsewhere.

###### Acknowledgements.

We acknowledge the use of the Chalawan High-Performance Computing cluster,
operated and maintained by the National Astronomical Research Institute of
Thailand (NARIT). The research of LAU-L was partially supported by Programa
para el Desarrollo Profesional Docente; Dirección de Apoyo a la Investigación
y al Posgrado, Universidad de Guanajuato; CONACyT México under Grants No.
A1-S-17899, No. 286897, No. 297771, No. 304001; and the Instituto Avanzado de
Cosmología Collaboration.

##  Appendix A Stability Analysis of the System

Here, we discuss the stability of the general phase space of the quintom
models. To perform the stability analysis, we need to transfer the systems to
a slightly different form by considering $\Omega_{\phi}^{1/2}=r_{\phi}$ and
$\Omega_{\sigma}^{1/2}=r_{\sigma}$. This is done to avoid the singularity in
the Jacobian matrix associated with the fixed points that have
$\Omega_{\phi}=\Omega_{\sigma}=0$. Hence, Eqs. (13) reduce to the following
set of equations,

$\displaystyle\theta_{\phi}^{\prime}$ $\displaystyle=$
$\displaystyle-3\sin\theta_{\phi}+y_{1\phi}\,,$ (20a) $\displaystyle
y_{1\phi}^{\prime}$ $\displaystyle=$ $\displaystyle
3/2\gamma_{tot}y_{1\phi}+r_{\phi}\sin(\theta_{\phi}/2)y_{2\phi}\,,$ (20b)
$\displaystyle r_{\phi}^{\prime}$ $\displaystyle=$ $\displaystyle
3/2(\gamma_{tot}-\gamma_{\phi})r_{\phi}\,,$ (20c)
$\displaystyle\theta_{\sigma}^{\prime}$ $\displaystyle=$
$\displaystyle-3\sinh\theta_{\sigma}-y_{1\sigma}\,,$ (20d) $\displaystyle
y_{1\sigma}^{\prime}$ $\displaystyle=$ $\displaystyle
3/2\gamma_{tot}y_{1\sigma}+r_{\sigma}\sinh(\theta_{\sigma}/2)y_{2\sigma}\,,$
(20e) $\displaystyle r_{\sigma}^{\prime}$ $\displaystyle=$ $\displaystyle
3/2(\gamma_{tot}-\gamma_{\sigma})r_{\sigma}\,.$ (20f)

###  A.1 Fixed Points and Eigenvalues

The fixed points of the systems are listed in Table 3, and were obtained
assuming that the contribution of the radiation component is negligible,
$\Omega_{r}\ll 1$, so the total barotropic EOS is given by
$\gamma_{tot}=1+(\gamma_{\phi}-1)r^{2}_{\phi}+(\gamma_{\sigma}-1)r^{2}_{\sigma}$.
The values of the coefficients $A$ and $B$ are as follows:

$\displaystyle A_{\pm}=$
$\displaystyle\frac{-\alpha_{1}\pm\sqrt{\alpha_{1}^{2}-2\alpha_{0}(1+2\alpha_{2})}}{6(1+2\alpha_{2})}$
(21a) $\displaystyle B_{\pm}=$
$\displaystyle\frac{\beta_{1}\pm\sqrt{\beta_{1}^{2}-2\beta_{0}(1+2\beta_{2})}}{6(1+2\beta_{2})}$
(21b)

The corresponding values of the cosmological parameters at the aforementioned
fixed points are given in Table 4.

Table 3: Complete list of fixed points in the autonomous system (20). Fixed points with the name $m_{i}$ correspond to the matter dominated fixed point, $q_{i}$ corresponds to the quintessence dominated fixed point and $p_{i}$ for phatom domination. | Fixed
---
Points
$r_{\phi_{c}}$ | $\sin\theta_{\phi_{c}}$ | $r_{\sigma_{c}}$ | $\sinh\theta_{\sigma_{c}}$
$m_{1}$ | 0 | $n\pi$ | 0 | 0
$m_{2}$ | 0 | $n\pi$ | 0 | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{\gamma_{tot}}{4\beta_{2}}$
$m_{3}$ | 0 | $\sin^{2}(\theta_{\phi_{c}}/2)=-\frac{\gamma_{tot}}{4\alpha_{2}}$ | 0 | 0
$m_{4}$ | 0 | $\sin^{2}(\theta_{\phi_{c}}/2)=-\frac{\gamma_{tot}}{4\alpha_{2}}$ | 0 | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{\gamma_{tot}}{4\beta_{2}}$
$q_{1}$ | 1 | $n\pi$ | 0 | 0
$q_{2}$ | 1 | $n\pi$ | 0 |  | $\sinh^{2}(\theta_{\sigma c}/2)$
---
$=\frac{1}{2\beta_{2}}\sin^{2}(\theta_{\phi c}/2)$
$q_{3}$ | 1 |  | $\sin(\theta_{\phi_{c}}/2)=A_{+}$
---
0 | 0
$q_{4}$ | 1 |  | $\sin(\theta_{\phi_{c}}/2)=A_{-}$
---
0 | 0
$q_{5}$ | 1 |  | $\sin(\theta_{\phi_{c}}/2)=A_{+}$
---
0 |  | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{{A^{2}_{+}}}{2\beta_{2}}$
---
$q_{6}$ | 1 |  | $\sin(\theta_{\phi_{c}}/2)=A_{-}$
---
0 |  | $\sinh^{2}(\theta_{\sigma_{c}}/2)=\frac{{A^{2}_{-}}}{2\beta_{2}}$
---
$p_{1}$ | 0 | $n\pi$ | 1 | 0
$p_{2}$ | 0 | $n\pi$ | 1 | $\sinh(\theta_{\sigma c}/2)=B_{+}$
$p_{3}$ | 0 | $n\pi$ | 1 | $\sinh(\theta_{\sigma c}/2)=B_{-}$
$p_{4}$ | 0 |  | $\sin^{2}(\theta_{\phi_{c}}/2)=\frac{B^{2}_{+}}{2\alpha_{2}}$
---
1 | $\sinh(\theta_{\sigma c}/2)=B_{+}$
$p_{5}$ | 0 |  | $\sin^{2}(\theta_{\phi_{c}}/2)=\frac{B^{2}_{-}}{2\alpha_{2}}$
---
1 | $\sinh(\theta_{\sigma c}/2)=B_{-}$
Table 4: Value of different cosmological parameters; $\Omega_{\phi},\Omega_{\sigma},\gamma_{\phi},\gamma_{\sigma},\gamma_{tot}$ corresponding to different fixed points given in Table 3. | Fixed
---
Points
$\Omega_{\phi}$ | $\gamma_{\phi}$ | $\Omega_{\sigma}$ | $\gamma_{\sigma}$ | $\gamma_{tot}$
$m_{1}$ | 0 |  | 0; for $n$ even
---
2; for $n$ odd
0 | 0 | 1
$m_{2}$ | 0 |  | 0; for $n$ even
---
2; for $n$ odd
0 | $-\frac{\gamma_{tot}}{2\beta_{2}}$ | 1
$m_{3}$ | 0 | $-\frac{\gamma_{tot}}{2\alpha_{2}}$ | 0 | 0 | 1
$m_{4}$ | 0 | $-\frac{\gamma_{tot}}{2\alpha_{2}}$ | 0 | $-\frac{\gamma_{tot}}{2\beta_{2}}$ | 1
$q_{1}$ | 1 |  | 0; for $n$ even
---
2; for $n$ odd
0 | 0 | $\gamma_{\phi}$
$q_{2}$ | 1 |  | 0; for $n$ even
---
2; for $n$ odd
0 | $-\frac{\gamma_{\phi}}{2\beta_{2}}$ | $\gamma_{\phi}$
$q_{3},q_{4}$ | 1 | $2A^{2}_{\pm}$ | 0 | 0 | $\gamma_{\phi}$
$q_{5},q_{6}$ | 1 | $2A^{2}_{\pm}$ | 0 | $-\frac{A^{2}_{\pm}}{\beta_{2}}$ | $\gamma_{\phi}$
$p_{1}$ | 0 |  | 0; for $n$ even
---
2; for $n$ odd
1 | 0 | $\gamma_{\sigma}$
$p_{2},p_{3}$ | 0 |  | 0; for $n$ even
---
2; for $n$ odd
1 | $-2B^{2}_{\pm}$ | $\gamma_{\sigma}$
$p_{4},p_{5}$ | 0 | $\frac{B^{2}_{\pm}}{2\alpha_{2}}$ | 1 | $-2B^{2}_{\pm}$ | $\gamma_{\sigma}$
Table 5: The table displays eigenvalues for fixed points listed in Table 3. | Fixed Points
---
Eigenvalues
$m_{1}$ | $-3,-3,\frac{3}{2},\frac{3}{2},\frac{3}{2},\frac{3}{2}$
$m_{2}$ |  | $-3,\frac{3}{2},\frac{3}{2},\frac{3}{2}\left(1+\frac{1}{2\beta_{2}}\right),\frac{3\left(-1-3\beta_{2}-\sqrt{1+2\beta_{2}-7\beta_{2}^{2}}\right)}{4\beta_{2}},\frac{3\left(-1-3\beta_{2}+\sqrt{1+2\beta_{2}-7\beta_{2}^{2}}\right)}{4\beta_{2}}$
---
$m_{3}$ | $-3,\frac{3}{2},\frac{3}{2},\frac{3}{2}\left(1-\frac{1}{2\alpha_{2}}\right),\frac{3\left(1+\alpha_{2}-\sqrt{1-10\alpha_{2}+33\alpha_{2}^{2}}\right)}{4\alpha_{2}},\frac{3\left(1+\alpha_{2}+\sqrt{1-10\alpha_{2}+33\alpha_{2}^{2}}\right)}{4\alpha_{2}}$
$m_{4}$ | $0,0,-\frac{3}{2}(2+\frac{1}{\alpha_{2}}),\frac{3}{4}\left(2+\frac{1}{\alpha_{2}}\right),-\frac{3\left(1+3\beta_{2}+\sqrt{1+2\beta_{2}+9\beta_{2}^{2}}\right)}{4\beta_{2}},\frac{3\left(-1-3\beta_{2}+\sqrt{1+2\beta_{2}+9\beta_{2}^{2}}\right)}{4\beta_{2}}$
$q_{1},q_{2}$ | $-3,-3,0,0,\frac{1}{2}\left(-3-\sqrt{9+2\alpha_{0}}\right),\frac{1}{2}\left(-3+\sqrt{9+2\alpha_{0}}\right)$
$q_{3},q_{4}$ | $-3,3A_{\pm}^{2},3A_{\pm}^{2},-3\left(1-2A_{\pm}^{2}\right),Q_{1},Q_{2}$
$q_{5},q_{6}$ | $-3\left(1-2A_{\pm}^{2}\right),\frac{3}{2}\left(2A_{\pm}^{2}+\frac{A_{\pm}}{\beta_{2}}\right),Q_{1},Q_{2},Q_{3},Q_{4}$
$p_{1}$ | $-3,-3,0,0,\frac{1}{2}\left(-3-\sqrt{9-2\beta_{0}}\right),\frac{1}{2}\left(-3+\sqrt{9-2\beta_{0}}\right)$
$p_{2},p_{3}$ | $-3,-3B_{\pm}^{2},-3B_{\pm}^{2},-3\left(1+2B_{\pm}^{2}\right),P_{1},P_{2}$
$p_{4},p_{5}$ | $-3\left(1+2B_{\pm}^{2}\right),\frac{3}{2}\left(-2B_{\pm}^{2}-\frac{B_{\pm}}{\alpha_{2}}\right),P_{1},P_{2},P_{3},P_{4}$

Fixed points can be classified into three different classes: i) matter
dominated, ii) quintessence dominated, and iii) phantom dominated (see Table
4). Fixed points dominated by matter are represented by $m_{i}$, points
dominated by quintessence are represented by $q_{i}$, and points dominated by
phantoms are represented by $p_{i}$.

###  A.2 Eigenvalues

The eigenvalues of these fixed points are listed in Table 5. To find them, we
have considered $n$ to be even for the fixed points where $\sin\theta_{\phi
c}=n\pi$, and the corresponding $\gamma_{\phi}=0$ that is of physical
interest. Odd-$n$ leads to $\gamma_{\phi}=2$ in which the quintessence field
behaves as a stiff fluid. For the complete expressions of
$Q_{1},Q_{2},Q_{3},Q_{4},P_{1},P_{2},P_{3},P_{4}$ see Sec. A.2.

Certain fixed points may not be hyperbolic, meaning that some of their
eigenvalues are zero. This means that the standard linear stability analysis
cannot be used to determine their stability. If one or more of the non-zero
eigenvalues are positive, then the center manifold theorem implies that the
fixed point is likely to be unstable. If all the non-zero eigenvalues are
negative, further investigation is necessary to determine the stability of the
fixed point.

###  A.3 Stability of the matter dominated points

The points $m_{1}$ to $m_{3}$ have a positive eigenvalue, making them unstable
fixed points. The point $m_{4}$ is a nonhyperbolic fixed point, as it has zero
eigenvalues. It can be easily seen that for the fixed points mentioned, all
eigenvalues cannot be simultaneously negative for any combination of
$\alpha_{2}$ and $\beta_{2}$. Consequently, they are also unstable.

###  A.4 Stability of the quintessence dominated points

The fixed points $q_{1}$ and $q_{2}$ may be late-time attractors, as all their
eigenvalues, apart from zero, can be negative for $-9/2\leq\alpha_{0}<0$. On
the other hand, $q_{3}$ and $q_{4}$ are unstable, as they have positive
eigenvalues. It is difficult to make a general statement about the stability
of $q_{5}$ and $q_{6}$, as the forms of some of their eigenvalues are very
complex. The stability of these fixed points is dependent on the selection of
the parameters $\alpha$ and $\beta$ in a complicated manner.

###  A.5 Stability of the phantom dominated points

The fixed point $p_{1}$ may have all negative eigenvalues in the range of
$0<\beta_{0}\leq 9/2$, making it stable in this particular interval. The
eigenvalues of the fixed points $p_{2}$ to $p_{5}$ are complex functions of
the parameters $\alpha$ and $\beta$, making their stability difficult to
determine.

$Q_{1,2}=\frac{1}{4}\left(-6+18A_{\pm}^{2}+2A_{\pm}\alpha_{1}+24A_{\pm}^{2}\alpha_{2}\right)\\\
\pm\frac{1}{4}\left(\sqrt{\left(6-18A_{\pm}^{2}-2A_{\pm}\alpha_{1}-24A_{\pm}^{2}\alpha_{2}\right){}^{2}-8\left(-54A_{\pm}^{2}+72A_{\pm}^{4}-\alpha_{0}+2A_{\pm}^{2}\alpha_{0}-12A_{\pm}\alpha_{1}+18A_{\pm}^{3}\alpha_{1}-108A_{\pm}^{2}\alpha_{2}+144A_{\pm}^{4}\alpha_{2}\right)}\right)$
(22)

$Q_{3,4}=-\frac{3A_{\pm}+3\beta_{2}+6A_{\pm}\beta_{2}-3A_{\pm}^{2}\beta_{2}\pm\sqrt{\left(-3A_{\pm}-3\beta_{2}-6A_{\pm}\beta_{2}+3A_{\pm}^{2}\beta_{2}\right){}^{2}+4\beta_{2}\left(-18A_{\pm}^{2}+9A_{\pm}^{3}-27A_{\pm}\beta_{2}+9A_{\pm}^{2}\beta_{2}\right)}}{2\beta_{2}}$
(23)

$P_{1,2}=\frac{1}{4}\left(-6-18B_{\pm}^{2}+2B_{\pm}\beta_{1}-24B_{\pm}^{2}\beta_{2}\right)\\\
\pm\frac{1}{4}\left(\sqrt{\left(6+18B_{\pm}^{2}-2B_{\pm}\beta_{1}+24B_{\pm}^{2}\beta_{2}\right){}^{2}-8\left(54B_{\pm}^{2}+72B_{\pm}^{4}+\beta_{0}+2B_{\pm}^{2}\beta_{0}-12B_{\pm}\beta_{1}-18B_{\pm}^{3}\beta_{1}+108B_{\pm}^{2}\beta_{2}+144B_{\pm}^{4}\beta_{2}\right)}\right)$
(24)

$P_{3,4}=\frac{3B_{\pm}-3\alpha_{2}+6B_{\pm}\alpha_{2}-3B_{\pm}^{2}\alpha_{2}\pm\sqrt{4\alpha_{2}\left(-18B_{\pm}^{2}+9B_{\pm}^{3}+27B_{\pm}\alpha_{2}-9B_{\pm}^{2}\alpha_{2}\right)+\left(3B_{\pm}-3\alpha_{2}+6B_{\pm}\alpha_{2}-3B_{\pm}^{2}\alpha_{2}\right){}^{2}}}{2\alpha_{2}}$
(25)

###  A.6 General discussion on stability

The stability of fixed points is contingent on the selection of active
parameters $\alpha,\beta$. A qualitative analysis of the entire phase space of
the system was conducted, which revealed that all matter-dominated fixed
points are unstable, suggesting that the universe may have gone through some
of these fixed points. The probable late-time attractors of the universe are
either quintessence-dominated or phantom-dominated. However, the analysis of a
general case is complicated due to the dependency of the eigenvalues on the
values of $\alpha$ and $\beta$. To simplify this, a particular functional form
of the quintessence and phantom potentials could be chosen, which would fix
the values of their active parameters and eigenvalues. This would enable a
detailed study of the dynamics of any particular case of interest.

## References

  * Riess  _et al._ (2022) A. G. Riess _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Scolnic  _et al._ (2023) D. Scolnic, A. G. Riess, J. Wu, S. Li, G. S. Anand, R. Beaton, S. Casertano, R. I. Anderson, S. Dhawan, and X. Ke, Astrophys. J. Lett. 954, L31 (2023), arXiv:2304.06693 [astro-ph.CO] .
  * Riess  _et al._ (1998) A. G. Riess _et al._ (Supernova Search Team), Astron. J. 116, 1009 (1998), arXiv:astro-ph/9805201 [astro-ph] .
  * Perlmutter  _et al._ (1999) S. Perlmutter _et al._ (Supernova Cosmology Project), Astrophys. J. 517, 565 (1999), arXiv:astro-ph/9812133 [astro-ph] .
  * Aghanim  _et al._ (2018) N. Aghanim _et al._ (Planck), (2018), arXiv:1807.06209 [astro-ph.CO] .
  * Weinberg (1989) S. Weinberg, Rev. Mod. Phys. 61, 1 (1989), [,569(1988)].
  * Peebles and Ratra (2003) P. J. E. Peebles and B. Ratra, Rev. Mod. Phys. 75, 559 (2003), [,592(2002)], arXiv:astro-ph/0207347 [astro-ph] .
  * Padmanabhan (2003) T. Padmanabhan, Phys. Rept. 380, 235 (2003), arXiv:hep-th/0212290 [hep-th] .
  * Alam  _et al._ (2017a) S. Alam _et al._ (BOSS), Mon. Not. Roy. Astron. Soc. 470, 2617 (2017a), arXiv:1607.03155 [astro-ph.CO] .
  * Joudaki  _et al._ (2019) S. Joudaki _et al._ , (2019), arXiv:1906.09262 [astro-ph.CO] .
  * Riess  _et al._ (2016) A. G. Riess _et al._ , Astrophys. J. 826, 56 (2016), arXiv:1604.01424 [astro-ph.CO] .
  * Riess  _et al._ (2019) A. G. Riess, S. Casertano, W. Yuan, L. M. Macri, and D. Scolnic, Astrophys. J. 876, 85 (2019), arXiv:1903.07603 [astro-ph.CO] .
  * Wong  _et al._ (2019) K. C. Wong _et al._ , (2019), arXiv:1907.04869 [astro-ph.CO] .
  * Freedman  _et al._ (2019) W. L. Freedman _et al._ ,  (2019), 10.3847/1538-4357/ab2f73, arXiv:1907.05922 [astro-ph.CO] .
  * Schöneberg  _et al._ (2022) N. Schöneberg, G. F. Abellán, A. P. Sánchez, S. J. Witte, V. Poulin, and J. Lesgourgues, Physics Reports 984, 1 (2022).
  * Verde  _et al._ (2019) L. Verde, T. Treu, and A. G. Riess (2019) arXiv:1907.10625 [astro-ph.CO] .
  * Amendola and Tsujikawa (2010) L. Amendola and S. Tsujikawa,  _Dark energy: theory and observations_ (Cambridge University Press, 2010).
  * Bamba  _et al._ (2012) K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
  * Copeland  _et al._ (2006) E. J. Copeland, M. Sami, and S. Tsujikawa, International Journal of Modern Physics D 15, 1753 (2006), arXiv:hep-th/0603057 [hep-th] .
  * Padmanabhan (2006) T. Padmanabhan, in  _Albert Einstein Century International Conference_, American Institute of Physics Conference Series, Vol. 861, edited by J.-M. Alimi and A. Füzfa (2006) pp. 179–196, arXiv:astro-ph/0603114 [astro-ph] .
  * Karwal and Kamionkowski (2016) T. Karwal and M. Kamionkowski, Phys. Rev. D 94, 103523 (2016), arXiv:1608.01309 [astro-ph.CO] .
  * Sola Peracaula  _et al._ (2023) J. Sola Peracaula, A. Gomez-Valent, J. de Cruz Perez, and C. Moreno-Pulido, Universe 9, 262 (2023), arXiv:2304.11157 [astro-ph.CO] .
  * Rezaei and Sola Peracaula (2022) M. Rezaei and J. Sola Peracaula, Eur. Phys. J. C 82, 765 (2022), arXiv:2207.14250 [astro-ph.CO] .
  * Rezaei  _et al._ (2019) M. Rezaei, M. Malekjani, and J. Sola, Phys. Rev. D 100, 023539 (2019), arXiv:1905.00100 [gr-qc] .
  * Rezaei  _et al._ (2021) M. Rezaei, J. Solà Peracaula, and M. Malekjani, Mon. Not. Roy. Astron. Soc. 509, 2593 (2021), arXiv:2108.06255 [astro-ph.CO] .
  * Di Valentino  _et al._ (2020) E. Di Valentino, A. Mukherjee, and A. A. Sen, (2020), arXiv:2005.12587 [astro-ph.CO] .
  * Roy  _et al._ (2022) N. Roy, S. Goswami, and S. Das, Phys. Dark Univ. 36, 101037 (2022), arXiv:2201.09306 [astro-ph.CO] .
  * Ratra and Peebles (1988) B. Ratra and P. J. E. Peebles, Phys. Rev. D37, 3406 (1988).
  * Tsujikawa (2013) S. Tsujikawa, Class. Quant. Grav. 30, 214003 (2013), arXiv:1304.1961 [gr-qc] .
  * Nojiri  _et al._ (2005) S. Nojiri, S. D. Odintsov, and S. Tsujikawa, Phys. Rev. D71, 063004 (2005), arXiv:hep-th/0501025 [hep-th] .
  * Caldwell  _et al._ (2003) R. R. Caldwell, M. Kamionkowski, and N. N. Weinberg, Phys. Rev. Lett. 91, 071301 (2003), arXiv:astro-ph/0302506 [astro-ph] .
  * Caldwell (2002) R. R. Caldwell, Phys. Lett. B545, 23 (2002), arXiv:astro-ph/9908168 [astro-ph] .
  * Roy and Bhadra (2018) N. Roy and N. Bhadra, Journal of Cosmology and Astroparticle Physics 2018, 002 (2018).
  * Gerardi  _et al._ (2019) F. Gerardi, M. Martinelli, and A. Silvestri, Journal of Cosmology and Astroparticle Physics 2019, 042 (2019).
  * Johri (2004) V. B. Johri, Physical Review D 70, 041303 (2004).
  * Perković and Štefančić (2019) D. Perković and H. Štefančić, International Journal of Modern Physics D 28, 1950045 (2019).
  * Guo  _et al._ (2005) Z.-K. Guo, Y.-S. Piao, X. Zhang, and Y.-Z. Zhang, Physics Letters B 608, 177 (2005).
  * Cai  _et al._ (2010) Y.-F. Cai, E. N. Saridakis, M. R. Setare, and J.-Q. Xia, Physics Reports 493, 1 (2010).
  * Zhao and Zhang (2006) W. Zhao and Y. Zhang, Physical Review D 73, 123509 (2006).
  * Setare and Saridakis (2008) M. Setare and E. Saridakis, Physics Letters B 668, 177 (2008).
  * Zhao  _et al._ (2005) G.-B. Zhao, J.-Q. Xia, M. Li, B. Feng, and X. Zhang, Physical Review D 72, 123515 (2005).
  * Setare and Saridakis (2009) M. Setare and E. Saridakis, International Journal of Modern Physics D 18, 549 (2009).
  * Panpanich  _et al._ (2019) S. Panpanich, P. Burikham, S. Ponglertsakul, and L. Tannukij, arXiv preprint arXiv:1908.03324 (2019).
  * Roy  _et al._ (2018) N. Roy, A. X. Gonzalez-Morales, and L. A. Ureña-López, Phys. Rev. D 98, 063530 (2018), arXiv:1803.09204 [gr-qc] .
  * Roy and Bamba (2019) N. Roy and K. Bamba, Physical Review D 99, 123520 (2019).
  * Linares Cedeño  _et al._ (2021a) F. X. Linares Cedeño, N. Roy, and L. A. Ureña López, Phys. Rev. D 104, 123502 (2021a), arXiv:2105.07103 [astro-ph.CO] .
  * Urena-Lopez and Roy (2020) L. A. Urena-Lopez and N. Roy, (2020), arXiv:2007.08873 [astro-ph.CO] .
  * Roy (2023) N. Roy, Gen. Rel. Grav. 55, 115 (2023), arXiv:2302.10509 [astro-ph.CO] .
  * Steinhardt  _et al._ (1999) P. J. Steinhardt, L.-M. Wang, and I. Zlatev, Phys. Rev. D59, 123504 (1999), arXiv:astro-ph/9812313 [astro-ph] .
  * Copeland  _et al._ (1998) E. J. Copeland, A. R. Liddle, and D. Wands, Phys. Rev. D57, 4686 (1998), arXiv:gr-qc/9711068 [gr-qc] .
  * Copeland  _et al._ (1999) E. J. Copeland, A. Mazumdar, and N. Nunes, Physical Review D 60, 083506 (1999).
  * Ivashchuk  _et al._ (2003) V. D. Ivashchuk, V. N. Melnikov, and A. B. Selivanov, Journal of High Energy Physics 2003, 059 (2003).
  * Bag  _et al._ (2018) S. Bag, S. S. Mishra, and V. Sahni, JCAP 1808, 009 (2018), arXiv:1709.09193 [gr-qc] .
  * Lesgourgues (2011a) J. Lesgourgues, arXiv preprint arXiv:1104.2932 (2011a).
  * Lesgourgues (2011b) J. Lesgourgues, (2011b), arXiv:1104.2932 [astro-ph.IM] .
  * Lesgourgues (2011c) J. Lesgourgues, (2011c), arXiv:1104.2934 [astro-ph.CO] .
  * Lesgourgues and Tram (2011) J. Lesgourgues and T. Tram, JCAP 1109, 032 (2011), arXiv:1104.2935 [astro-ph.CO] .
  * Audren  _et al._ (2013) B. Audren, J. Lesgourgues, K. Benabed, and S. Prunet, JCAP 1302, 001 (2013), arXiv:1210.7183 [astro-ph.CO] .
  * Brout  _et al._ (2022) D. Brout _et al._ , Astrophys. J. 938, 110 (2022), arXiv:2202.04077 [astro-ph.CO] .
  * Alam  _et al._ (2017b) S. Alam, M. Ata, S. Bailey, F. Beutler, D. Bizyaev, J. A. Blazek, A. S. Bolton, J. R. Brownstein, A. Burden, C.-H. Chuang, and et al., Monthly Notices of the Royal Astronomical Society 470, 2617–2652 (2017b).
  * Cuceu  _et al._ (2019) A. Cuceu, J. Farr, P. Lemos, and A. Font-Ribera, Journal of Cosmology and Astroparticle Physics 2019, 044–044 (2019).
  * Kazin  _et al._ (2014) E. A. Kazin, J. Koda, C. Blake, N. Padmanabhan, S. Brough, M. Colless, C. Contreras, W. Couch, S. Croom, D. J. Croton, and et al., Monthly Notices of the Royal Astronomical Society 441, 3524–3542 (2014).
  * Arendse  _et al._ (2020) N. Arendse _et al._ , Astron. Astrophys. 639, A57 (2020), arXiv:1909.07986 [astro-ph.CO] .
  * Sabti  _et al._ (2022) N. Sabti, J. B. Muñoz, and D. Blas, Phys. Rev. D 105, 043518 (2022), arXiv:2110.13168 [astro-ph.CO] .
  * Ureña-López (2019) L. A. Ureña-López, JCAP 1906, 009 (2019), arXiv:1904.03318 [astro-ph.CO] .
  * Linares Cedeño  _et al._ (2021b) F. X. Linares Cedeño, A. X. González-Morales, and L. A. Ureña López, JCAP 01, 051 (2021b), arXiv:2006.05037 [astro-ph.CO] .
  * Lewis (2019) A. Lewis,  (2019), arXiv:1910.13970 [astro-ph.IM] .
  * Camarena and Marra (2018) D. Camarena and V. Marra, Phys. Rev. D 98, 023537 (2018), arXiv:1805.09900 [astro-ph.CO] .
  * Trotta (2007) R. Trotta, Mon. Not. Roy. Astron. Soc. 378, 72 (2007), arXiv:astro-ph/0504022 .
  * Heavens  _et al._ (2017) A. Heavens, Y. Fantaye, A. Mootoovaloo, H. Eggers, Z. Hosenie, S. Kroon, and E. Sellentin, (2017), arXiv:1704.03472 [stat.CO] .
  * Scolnic  _et al._ (2018) D. Scolnic _et al._ , Astrophys. J. 859, 101 (2018), arXiv:1710.00845 [astro-ph.CO] .


* Vermeulen and Vansteelandt (2015) [author] Vermeulen, KarelK. and Vansteelandt, StijnS. (2015). Bias-Reduced Doubly Robust Estimation. Journal of the American Statistical Association 110 1024–1036.
  * Wager and Walther (2016) [author] Wager, StefanS. and Walther, GuentherG. (2016). Adaptive Concentration of Regression Trees, with Application to Random Forests. arXiv preprint arXiv:1503.06388.
  * Wang, Agarwal and Dudik (2017) [author] Wang, Yu-XiangY.-X., Agarwal, AlekhA. and Dudik, MiroslavM. (2017). Optimal and adaptive off-policy evaluation in contextual bandits. In Proceedings of the 34th International Conference on Machine Learning 3589–3597.
  * Wang, Foster and Kakade (2020) [author] Wang, RuosongR., Foster, Dean P.D. P. and Kakade, Sham M.S. M. (2020). What are the Statistical Limits of Offline RL with Linear Function Approximation?. arXiv preprint arXiv:2010.11895.
  * Wang and Zubizarreta (2020) [author] Wang, YixinY. and Zubizarreta, Jose RJ. R. (2020). Minimal dispersion approximately balancing weights: asymptotic properties and practical considerations. Biometrika 107 93–105.
  * Wang et al. (2012) [author] Wang, LuL., Rotnitzky, AndreaA., Lin, XihongX., Millikan, Randall ER. E. and Thall, Peter FP. F. (2012). Evaluation of viable dynamic treatment regimes in a sequentially randomized trial of advanced prostate cancer. Journal of the American Statistical Association 107 493–508.
  * Wang et al. (2018) [author] Wang, LanL., Zhou, YuY., Song, RuiR. and Sherwood, BenB. (2018). Quantile-optimal treatment regimes. Journal of the American Statistical Association 113 1243–1254.
  * Wu and Wang (2021) [author] Wu, YunanY. and Wang, LanL. (2021). Resampling-based confidence intervals for model-free robust inference on optimal treatment regimes. Biometrics 77 465–476.
  * Xie, Ma and Wang (2019) [author] Xie, TengyangT., Ma, YifeiY. and Wang, Yu-XiangY.-X. (2019). Towards Optimal Off-Policy Evaluation for Reinforcement Learning with Marginalized Importance Sampling. In Advances in Neural Information Processing Systems 32 9665–9675.
  * Xie et al. (2021) [author] Xie, TengyangT., Cheng, Ching-AnC.-A., Jiang, NanN., Mineiro, PaulP. and Agarwal, AlekhA. (2021). Bellman-consistent Pessimism for Offline Reinforcement Learning. arXiv preprint arXiv:2106.06926.
  * Xu et al. (2018) [author] Xu, ZheZ., Li, ZhixinZ., Guan, QingwenQ., Zhang, DingshuiD., Li, QiangQ., Nan, JunxiaoJ., Liu, ChunyangC., Bian, WeiW. and Ye, JiepingJ. (2018). Large-scale order dispatch in on-demand ride-hailing platforms: A learning and planning approach. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining 905–913.
  * Yang et al. (2020) [author] Yang, MengjiaoM., Nachum, OfirO., Dai, BoB., Li, LihongL. and Schuurmans, DaleD. (2020). Off-policy evaluation via the regularized lagrangian. arXiv preprint arXiv:2007.03438.
  * Yin, Bai and Wang (2021) [author] Yin, MingM., Bai, YuY. and Wang, Yu-XiangY.-X. (2021). Near-Optimal Provable Uniform Convergence in Offline Policy Evaluation for Reinforcement Learning. In International Conference on Artificial Intelligence and Statistics 1567–1575. PMLR.
  * Yin and Wang (2020) [author] Yin, MingM. and Wang, Yu-XiangY.-X. (2020). Asymptotically Efficient Off-Policy Evaluation for Tabular Reinforcement Learning. In Proceedings of the 23nd International Workshop on Artificial Intelligence and Statistics.
  * Young, Hernán and Robins (2014) [author] Young, Jessica GJ. G., Hernán, Miguel AM. A. and Robins, James MJ. M. (2014). Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data. Epidemiologic methods 3 1–19.
  * Yu et al. (2020) [author] Yu, TianheT., Thomas, GarrettG., Yu, LantaoL., Ermon, StefanoS., Zou, James YJ. Y., Levine, SergeyS., Finn, ChelseaC. and Ma, TengyuT. (2020). MOPO: Model-based Offline Policy Optimization. In Advances in Neural Information Processing Systems 33 14129–14142.
  * Zhang and Bareinboim (2016) [author] Zhang, JunzheJ. and Bareinboim, EliasE. (2016). Markov decision processes with unobserved confounders: A causal approach Technical Report, Technical report, Technical Report R-23, Purdue AI Lab.
  * Zhang and Bareinboim (2021) [author] Zhang, J.J. and Bareinboim, E.E. (2021). Non-Parametric Methods for Partial Identification of Causal Effects. Columbia CausalAI Laboratory Technical Report (R-72).
  * Zhang et al. (2012) [author] Zhang, BaqunB., Tsiatis, Anastasios AA. A., Laber, Eric BE. B. and Davidian, MarieM. (2012). A robust method for estimating optimal treatment regimes. Biometrics 68 1010–1018.
  * Zhang et al. (2013) [author] Zhang, BaqunB., Tsiatis, Anastasios A.A. A., Laber, Eric B.E. B. and Davidian, MarieM. (2013). Robust estimation of optimal dynamic treatment regimes for sequential treatment decisions. Biometrika 100 681–694.
  * Zhang et al. (2015) [author] Zhang, YichiY., Laber, Eric BE. B., Tsiatis, AnastasiosA. and Davidian, MarieM. (2015). Using decision lists to construct interpretable and parsimonious treatment regimes. Biometrics 71 895–904.
  * Zhang et al. (2018) [author] Zhang, YichiY., Laber, Eric BE. B., Davidian, MarieM. and Tsiatis, Anastasios AA. A. (2018). Interpretable dynamic treatment regimes. Journal of the American Statistical Association 113 1541–1549.
  * Zhang et al. (2021) [author] Zhang, Michael RM. R., Paine, Tom LeT. L., Nachum, OfirO., Paduraru, CosminC., Tucker, GeorgeG., Wang, ZiyuZ. and Norouzi, MohammadM. (2021). Autoregressive Dynamics Models for Offline Policy Evaluation and Optimization. arXiv preprint arXiv:2104.13877.
  * Zhang et al. (2022) [author] Zhang, RuiqiR., Zhang, XuezhouX., Ni, ChengzhuoC. and Wang, MengdiM. (2022). Off-Policy Fitted Q-Evaluation with Differentiable Function Approximators: Z-Estimation and Inference Theory. arXiv preprint arXiv:2202.04970.
  * Zhao et al. (2012a) [author] Zhao, YingqiY., Zeng, DonglinD., Rush, A JohnA. J. and Kosorok, Michael RM. R. (2012a). Estimating individualized treatment rules using outcome weighted learning. Journal of the American Statistical Association 107 1106–1118.
  * Zhao et al. (2012b) [author] Zhao, YingqiY., Zeng, DonglinD., Rush, A. JohnA. J. and Kosorok, Michael RM. R. (2012b). Estimating Individualized Treatment Rules Using Outcome Weighted Learning. Journal of the American Statistical Association 107 1106–1118.
  * Zheng and van Der Laan (2011) [author] Zheng, WenjingW. and van Der Laan, Mark JM. J. (2011). Cross-Validated Targeted Minimum-Loss-Based Estimation. In Targeted Learning: Causal Inference for Observational and Experimental Data. Springer Series in Statistics 459–474. Springer New York, New York, NY.
  * Zhu, Zeng and Song (2019) [author] Zhu, WenshengW., Zeng, DonglinD. and Song, RuiR. (2019). Proper inference for value function in high-dimensional Q-learning for dynamic treatment regimes. Journal of the American Statistical Association 114 1404–1417.
  * Zhu et al. (2017) [author] Zhu, RuoqingR., Zhao, Ying-QiY.-Q., Chen, GuanhuaG., Ma, ShuanggeS. and Zhao, HongyuH. (2017). Greedy outcome weighted tree learning of optimal personalized treatment rules. Biometrics 73 391–400.