\documentclass[fleqn,10pt, twocolumn]{wlscirep}

\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{subcaption}
\usepackage{amsmath}
\usepackage{amsbsy}
\usepackage{rotating}
\usepackage{booktabs}
\usepackage{bm}
\hypersetup{colorlinks,citecolor=blue}
\hypersetup{colorlinks=true,linkcolor=magenta,filecolor=magenta,    urlcolor=blue}
\renewcommand{\thefootnote}{\arabic{footnote}}

\title{DESI DR2 Favors Metastable in Emergent Dark Energy with No Evidence for Dark Matter–Dark Energy Interactions}

\author[1,2]{Himanshu Chaudhary}
\author[3,4,5]{Salvatore Capozziello}
\author[6]{Dhruba Jyoti Gogoi}
\author[7]{G. Mustafa}
\author[8]{Chengxun Yuan}

\affil[1]{Department of Physics, Babeș-Bolyai University, Kogălniceanu Street, Cluj-Napoca, 400084, Romania}
\affil[2]{Research Center of Astrophysics and Cosmology, Khazar University, Baku, AZ1096, 41 Mehseti Street, Azerbaijan}

\affil[3]{Dipartimento di Fisica ``E. Pancini'', Università di Napoli ``Federico II'', Complesso Universitario di Monte Sant'Angelo, Edificio G, Via Cinthia, I-80126, Napoli, Italy}
\affil[4]{Istituto Nazionale di Fisica Nucleare (INFN), Sezione di Napoli, Via Cinthia 9, I-80126 Napoli, Italy}
\affil[5]{Scuola Superiore Meridionale, Via Mezzocannone 4, I-80134, Napoli, Italy}

\affil[6]{Department of Physics, Madhabdev University, Narayanpur, Lakhimpur 784164, Assam, India}

\affil[7]{Department of Physics, Zhejiang Normal University, Jinhua 321004, People's Republic of China}

\affil[8]{School of Physics, Harbin Institute of Technology, Harbin 150001, People's Republic of China}

\affil[*]{himanshu.chaudhary@ubbcluj.ro}
\affil[*]{himanshuch1729@gmail.com}
\affil[$\dagger$]{capozziello@na.infn.it}
\affil[$\ddagger$]{moloydhruba@yahoo.in}
\affil[$\S$]{gmustafa3828@gmail.com}
\affil[$\P$]{yuancx@hit.edu.cn}

\affil[+]{All authors contributed equally to this work}

%\keywords{Keyword1, Keyword2, Keyword3}

\begin{abstract}
The nature of dark energy remains one of the most fundamental open questions in modern cosmology. Among the proposed extensions of the standard cosmological model, interacting dark energy scenarios provide a physically motivated framework for exploring possible energy exchange within the dark sector and the dynamical evolution of cosmic acceleration. In this work, we present a comprehensive observational study of the Interacting Metastable Emergent Dark Energy (IMEDE) model using the Dark Energy Spectroscopic Instrument Data Release 2 (DESI DR2) baryon acoustic oscillation measurements, together with Planck PR4 CamSpec cosmic microwave background observations and three independent Type Ia supernova compilations: Pantheon$+$, DES-Dovekie, and Union3. We find that the IMEDE model provides a consistently better description of the DESI DR2 observations than the standard $\Lambda$CDM cosmology, with weak to moderate statistical significance based on the best-fit likelihood analysis and strong to decisive Bayesian evidence across the different dataset combinations. This statistical preference originates from a coherent reduction of the residuals across several independent BAO distance measurements rather than from any individual outlying data point. The reconstructed dark energy equation of state shows a transition from a phantom phase to a quintessence phase, crossing the phantom divide at $z\simeq0.3$--$0.5$, characteristic of a Quintom-B dark energy scenario. Although the IMEDE model does not resolve the $H_0$ or $S_8$ tensions and provides no evidence for an interaction between dark matter and dark energy, it is well supported by the current cosmological data sets and motivates further tests with forthcoming DESI DR3 and next-generation cosmological surveys.
\end{abstract}
\begin{document}

\flushbottom
\maketitle
% * <john.hammersley@gmail.com> 2015-02-09T12:07:31.197Z:
%
%  Click the title above to edit the author information and abstract
%
%\thispagestyle{empty}

%\noindent Please note: Abbreviations should be introduced at the first mention in the main text – no abbreviations lists. Suggested structure of main text (not enforced) is provided below.

\section*{Introduction}
The $\Lambda$ Cold Dark Matter ($\Lambda$CDM) model serves as the standard cosmological paradigm, providing a comprehensive description of large-scale structure, cosmic microwave background anisotropies, and the background expansion history of the Universe. This model assumes a spatially flat geometry dominated by CDM and a cosmological constant ($\Lambda$), which drives the accelerated expansion observed at late times \cite{Planck:2018vyg,SupernovaSearchTeam:1998fmf,Copeland:2006wr,Sahni:1999gb}. Despite its empirical success, the $\Lambda$CDM framework presents significant theoretical difficulties and increasing observational inconsistencies \cite{Park:2017xbl,Hu:2023jqc}.

Theoretically, the model is limited by the cosmological constant problem, characterized by the extensive magnitude difference between the observed $\Lambda$ value and the vacuum energy density predicted by quantum field theory, which introduces a strict fine-tuning requirement \cite{Tian:2019enx,SolaPeracaula:2024nsz,Diaz-Pachon:2024nsq}. Another limitation is the coincidence problem, which concerns the comparable magnitudes of dark energy and dark matter densities strictly in the present epoch \cite{Perivolaropoulos:2021jda,Velten:2014nra}. Observationally, the model shows discrepancies between early-universe predictions and late-time local measurements, specifically the Hubble ($H_0$) tension and the structure growth ($S_8$) tension \cite{di2021realm,verde2019tensions,di2025cosmoverse,Pantos:2026koc,wright2025kids,DES:2026zjp}. These issues have motivated the development of alternative theoretical formulations, such as dynamical dark energy models \cite{Yang:2021eud,Cai:2026swf,Dinda:2024kjf,Chevallier:2000qy,Jassal:2004ej}, interacting dark sectors \cite{Giare:2024smz,Li:2025owk,Clemson:2011an,Li:2014eha,Xia:2016vnp}, and modified theories of gravity \cite{Clifton:2011jh,Dolgov:2003px,KumarSharma:2022qdf,DeFelice:2010aj,Tsujikawa:2007xu,Sharma:2019yix}.

The DESI DR2 provides new insights into the nature of dark energy \cite{abdul2025desi}. When DESI DR2 is combined with CMB observations, the $w_0w_a$CDM model favors a dynamical dark energy scenario over the $\Lambda$CDM model, characterized by the Quintom-B behavior ($w_0>-1$, $w_a<0$, and $w_0+w_a<-1$) \cite{Cai:2025mas}. The inclusion of the Pantheon$+$, Union3, and DESY5 Type Ia supernova datasets further strengthens this preference, increasing the significance to 2.8$\sigma$, 3.8$\sigma$, and 4.2$\sigma$, respectively.  Later, the revised DESY5 dataset was recalibrated and released as DES-Dovekie \cite{DES:2025sig}. Using this updated supernova compilation, the preference for dynamical dark energy remains significant at $3.2\sigma$, although it is reduced compared to the $4.2\sigma$ significance obtained with the original DESY5 dataset.

This preference necessitates the re-evaluation of phenomenological frameworks, including emergent dark energy models. The Phenomenologically Emergent Dark Energy (PEDE) model \cite{li2019simple} was later extended to the Generalized Emergent Dark Energy (GEDE) parameterization \cite{Li:2020ybr} and subsequently to the Modified Emergent Dark Energy model \cite{benaoum2022modified}. In these models, the dark energy density is negligible in the early Universe, preserving the standard radiation- and matter-dominated epochs, and increases only at late times to drive the observed cosmic acceleration. This emergent behavior is introduced to alleviate the coincidence problem while providing a testable framework for investigating the $H_0$ and $S_8$ tensions.

Recently, \cite{Li:2026hwq} introduced the Metastable Emergent Dark Energy (MEDE) model to describe late-time dark energy emergence and accommodate a smooth phantom-crossing ($w=-1$). Constrained by Planck CMB, DESI DR2 BAO, and PantheonPlus data, the authors found a transition redshift of $z_t=0.425^{+0.084}_{-0.120}$ and a variation amplitude of $\Delta =0.87^{+0.29}_{-0.35}$. The MEDE framework demonstrates a statistically significant deviation from a static cosmological constant; it is notably preferred over $\Lambda$CDM ($\Delta \text{DIC} = -9.29$) and performs comparably to the CPL parameterization, offering a highly competitive, dynamical alternative that preserves established early-universe physics.

Another well-motivated extension of the standard cosmological model involves allowing a non-gravitational interaction between dark matter (DM) and dark energy (DE). Interacting dark energy (IDE) models have attracted considerable attention because they preserve the observational successes of the $\Lambda$CDM model while offering a possible mechanism to alleviate some of the existing cosmological tensions \cite{Kumar:2016zpg,Murgia:2016ccp,DiValentino:2017iww,Kumar:2021eev,Pan:2023mie,Benisty:2024lmj}. In these models, the energy exchange between the dark sector components is governed by an interaction term $Q$, usually parameterized through a dimensionless coupling parameter $\xi$. Such interactions can modify both the background expansion history and the growth of cosmic structures, providing a richer phenomenology than non-interacting dark energy models. This motivates extending the MEDE framework by incorporating a dark sector interaction, leading to the Interacting Metastable Emergent Dark Energy (IMEDE) model investigated in this work.

In the present work, we investigate the IMEDE framework through a Markov Chain Monte Carlo (MCMC) analysis using a combination of cosmological observations. Our dataset consists of Baryon Acoustic Oscillation (BAO) measurements from DESI DR2, Cosmic Microwave Background (CMB) observations from \textit{Planck} PR4 together with the \textit{Planck} PR4 and ACT DR6 lensing likelihoods, and three Type Ia supernova compilations: Pantheon$+$, DES-Dovekie, and Union3. This combination enables robust constraints on the late-time expansion history while reducing parameter degeneracies.

The primary objective of this work is to determine whether the IMEDE model can alleviate the $H_0$ and $S_8$ tensions while remaining consistent with current observational constraints. We further investigate whether introducing an interaction between dark matter and dark energy is supported by current cosmological observations and whether the IMEDE parameterization accommodates the dynamical phantom-crossing (Quintom-B) behavior indicated by the DESI DR2 background expansion data. Finally, through the evaluation of the interaction parameter $Q$ together with the metastable dark energy parameters ($z_t$ and $\Delta$), we quantify the statistical preference for the IMEDE model relative to the $\Lambda$CDM paradigm using Bayesian evidence.

The paper is structured as follows: Section~\ref{sec_2} presents the theoretical framework and mathematical formulation of the IMEDE model. Section~\ref{sec_3} describes the observational datasets and the MCMC methodology used in this work. Section~\ref{sec_4} presents the cosmological parameter constraints, investigates the implications for the $H_0$ and $S_8$ tensions, examines the statistical preference for an interaction between dark matter and dark energy, and discusses the evolution of the dark energy equation of state together with the overall cosmological implications of the IMEDE model. Finally, the main conclusions are summarized in Section~\ref{sec_5}.

\section{Theoretical background}\label{sec_2}
A possible interaction between dark matter and dark energy provides a well-motivated extension of the standard cosmological model. In the absence of a fundamental theory describing the dark sector, a phenomenological fluid description is commonly adopted, in which dark matter (DM) and dark energy (DE) are treated as two interacting fluids exchanging energy and momentum. Within this framework, the energy--momentum tensors of the two components are not conserved separately, but satisfy
\begin{equation}
    \nabla_\mu T^{\mu\nu}_{(k)}=\mathcal{Q}^{\nu}_{(k)},
\end{equation}
where the index $k$ denotes either DM or DE, and $\mathcal{Q}^{\nu}_{(k)}$ characterizes the energy--momentum exchange between the two sectors. The conservation of the total energy--momentum tensor, imposed by the Bianchi identities, requires
$\mathcal{Q}^{\nu}_{\rm de}=-\mathcal{Q}^{\nu}_{\rm dm}$.
In this work, we consider only energy exchange,
$\mathcal{Q}^{0}_{\rm de}=-\mathcal{Q}^{0}_{\rm dm}=\mathcal{Q}$,
while the momentum transfer is assumed to vanish, i.e.,
$\mathcal{Q}^{i}_{(k)}=0$.

Assuming dark matter to be cold and pressureless, the background continuity equations become
\begin{equation}
    \dot{\rho}_{\mathrm{dm}} + 3\mathcal{H}\rho_{\mathrm{dm}}
    = a^2 \mathcal{Q}^0_{\mathrm{dm}}
    = a\mathcal{Q},
\end{equation}
\begin{equation}
    \dot{\rho}_{\mathrm{de}} + 3\mathcal{H}(1+w)\rho_{\mathrm{de}}
    = a^2 \mathcal{Q}^0_{\mathrm{de}}
    = -a\mathcal{Q},
\end{equation}
where an overdot denotes differentiation with respect to conformal time, $\mathcal{H}=aH$ is the conformal Hubble parameter, and $w\equiv w(a)$ is the dark energy equation of state. A positive interaction term ($\mathcal{Q}>0$) corresponds to energy transfer from dark energy to dark matter, whereas $\mathcal{Q}<0$ describes the opposite direction of energy flow \cite{Olivares:2005tb}.

Following Refs.~\cite{Hoerning:2023hks,Valiviita:2008iv,Gavela:2010tm,DaCosta:2014vsr,deCesare:2021wmk}, we adopt the interaction
\begin{equation}
    \mathcal{Q}=H Q \rho_{\mathrm{de}},
\end{equation}
where $Q$ denotes the dimensionless coupling constant. The corresponding analytical solutions of the background continuity equations are
\begin{equation}
    \rho_{\mathrm{dm}}
    =\rho_{\mathrm{dm},0}a^{-3}
    +\rho_{\mathrm{de},0}a^{-3}
    Q\int_{1}^{a}
    \frac{a'^2e^{\xi(a')}}{a'^Q}\,
    {\rm d}a',
\end{equation}
\begin{equation}
    \rho_{\mathrm{de}}
    =\rho_{\mathrm{de},0}
    \frac{e^{\xi(a)}}{a^{Q}},
\end{equation}
where
\[
\xi(a)=\int_a^1
\frac{3[1+w(a')]}{a'}
\,{\rm d}a'.
\]

The linear perturbation equations for dark matter and dark energy are solved in the synchronous gauge and follow the formulation presented in Refs.~\cite{Hoerning:2023hks,Valiviita:2008iv,Gavela:2010tm,DaCosta:2014vsr,deCesare:2021wmk,Ma:1995ey},
\begin{align}\label{eq_7}
\dot{\delta}_{\mathrm{dm}}
=-\theta_{\mathrm{dm}}
-\frac{\dot h}{2}
+\mathcal{H}Q\frac{\rho_{\mathrm{de}}}{\rho_{\mathrm{dm}}}
(\delta_{\mathrm{de}}-\delta_{\mathrm{dm}})
&\\
+ Q\frac{\rho_{\mathrm{de}}}{\rho_{\mathrm{dm}}}
    \left(\frac{k v_T}{3}+\frac{\dot h}{6}\right),
\end{align}

\begin{equation}
    \dot{\theta}_{\mathrm{dm}}
    =-\mathcal{H}\theta_{\mathrm{dm}}
    \left(1+Q\frac{\rho_{\mathrm{de}}}{\rho_{\mathrm{dm}}}\right),
\end{equation}

\begin{align}
    \dot{\delta}_{\mathrm{de}}
    =&-(1+w)\left(\theta_{\mathrm{de}}+\frac{\dot h}{2}\right)
    -\frac{3\mathcal{H}\theta_{\mathrm{de}}}{k^2}\dot w
    -Q\left(\frac{k v_T}{3}+\frac{\dot h}{6}\right)\notag\\
    &-3\mathcal{H}(c_s^2-w)
    \left[
    \delta_{\mathrm{de}}
    +\{3(1+w)+Q\}
    \frac{\mathcal{H}\theta_{\mathrm{de}}}{k^2}
    \right],
\end{align}

\begin{equation}\label{eq_11}
    \dot{\theta}_{\mathrm{de}}
    =-\mathcal{H}(1-3c_s^2)\theta_{\mathrm{de}}
    +\frac{Q\mathcal{H}\theta_{\mathrm{de}}}{1+w}(1+c_s^2)
    +\frac{k^2c_s^2}{1+w}\delta_{\mathrm{de}}.
\end{equation}

Here, $\delta_{\mathrm{dm}}$ and $\delta_{\mathrm{de}}$ denote the density contrasts, while $\theta_{\mathrm{dm}}$ and $\theta_{\mathrm{de}}$ are the corresponding velocity divergences of dark matter and dark energy. The parameter $c_s^2$ is the physical sound speed of dark energy, the adiabatic sound speed is taken to be equal to the dark energy equation of state, and $v_T$ represents the centre-of-mass velocity of the total fluid, defined as $(1+w_T)v_T=\sum_a(1+w_a)\Omega_av_a$. The additional contributions in the density perturbation equations arise from the perturbation of the Hubble expansion rate, $\delta H$ \cite{Hoerning:2023hks}.
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Interacting Metastable Emergent Dark Energy}
In this work, we investigate the IMEDE model, in which the MEDE scenario is extended by allowing an interaction between cold dark matter and dark energy. Following \cite{Li:2026hwq}, the dark energy equation of state is given by
\begin{equation}
w(z)=-1-\Delta\tanh\!\left[\log_{10}\!\left(\frac{1+z}{1+z_t}\right)\right].
\end{equation}
Here, $z_t$ denotes the transition redshift at which the equation of state changes most rapidly, while $\Delta$ controls the amplitude of the transition from a phantom phase ($w<-1$) at earlier times to a quintessence phase ($w>-1$) at late times. Consequently, the IMEDE model naturally shows a Quintom-B evolution. In the limit $\Delta\rightarrow0$, the equation of state reduces to $w(z)=-1$, recovering the standard $\Lambda$CDM cosmology. The interaction between dark matter and dark energy is described by the coupling parameter $Q$, allowing us to investigate the combined effects of emergent dark energy dynamics and dark-sector interactions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure*}
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth]{fig_1.pdf}% width
\end{subfigure}
\hfil
\begin{subfigure}{.49\textwidth}
\includegraphics[width=\linewidth]{fig_2.pdf}% width
\end{subfigure}
\caption{This figure shows the CMB temperature power spectrum, $C_\ell^{TT}$ (left panel), together with the corresponding fractional difference, $\Delta C_\ell^{TT}/C_{\ell,\Lambda{\rm CDM}}^{TT}$ (lower-left panel). The right panel shows the linear matter power spectrum, $P(k)$, and its corresponding fractional difference, $\Delta P/P_{\Lambda{\rm CDM}}$ (lower-right panel). The black dashed curve shows the $\Lambda$CDM model, whereas the colored curves correspond to different values of the model parameter ($z_t$, $\Delta$, and $Q$), with all other cosmological parameters fixed.}\label{fig_1}
\end{figure*}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Dataset and Methodology}\label{sec_3}
To constrain the parameters of the IMEDE model, we perform an MCMC analysis \cite{lewis2002cosmological,lewis2013efficient,Neal:2005uqf} using the publicly available Bayesian cosmological parameter inference package {\tt Cobaya}\footnote{\url{https://github.com/CobayaSampler/cobaya}} \cite{torrado2021cobaya}. The convergence of the Markov chains is generated using the Gelman--Rubin statistic, requiring all chains to satisfy $R-1<0.01$ before the sampling is considered converged \cite{gelman1992inference}. The resulting posterior distributions are analyzed with the {\tt GetDist}\footnote{\url{https://github.com/cmbant/getdist}} package \cite{lewis2025getdist}, while the theoretical predictions are computed with the Boltzmann solver Code for Anisotropies in the Microwave Background ({\tt CAMB})\footnote{\url{https://github.com/cmbant/CAMB.git}} \cite{lewis2000efficient,howlett2012cmb}.

The baseline cosmological parameter vector is chosen as $\bm{\theta}_{\Lambda\mathrm{CDM}}
= \{\Omega_{\rm cdm},\,\Omega_{\rm b},\,100\,\theta_{\rm MC},\,\ln(10^{10}A_{\rm s}),\,n_{\rm s},\,\tau\}.$ For the IMEDE model, this parameter space is extended by introducing the transition redshift $z_t$, the transition width $\Delta$, and the interaction coupling parameter $Q$, such that $\bm{\theta}_{\rm IMEDE}
= \{\bm{\theta}_{\Lambda\mathrm{CDM}},\,z_t,\,\Delta,\,Q\}.$

Special care is taken in choosing the prior for the interaction parameter, $Q$. \cite{Gavela:2010tm} have shown that, in interacting dark energy models, the interaction parameter must have the opposite sign to $(1+w)$ to avoid early-time gravitational instabilities. Therefore, models with $w>-1$ require $\xi<0$, whereas models with $w<-1$ require $\xi>0$. Since the IMEDE model predicts a dynamical equation of state that smoothly crosses the phantom divide, we restrict our analysis to the stable prior range, $-3<\xi<-0.001$, corresponding to energy transfer from dark matter to dark energy. We also fix the dark energy sound speed to $c_{s,\mathrm{de}}^{2}=1$, a standard choice that avoids unphysical instabilities on small scales and ensures a stable evolution of cosmological perturbations. Exploring the alternative stability branch is beyond the scope of this work.

To constrain the parameters of the IMEDE model, we consider several combinations of early-and late-time cosmological observations. In particular, we use CamSpec Cosmic Microwave Background likelihood based on the Planck NPIPE PR4 release, Baryon Acoustic Oscillations, and Type Ia supernova data from Pantheon${+}$, DES-Dovekie, and Union3. In the following, we briefly describe each dataset and its implementation in our analysis.

\begin{itemize}
\item \textbf{Cosmic Microwave Background:} We first consider a combination of Cosmic Microwave Background (CMB) datasets. Our analysis includes the Planck temperature (\texttt{TT}), polarization (\texttt{EE}), and temperature--polarization cross-correlation (\texttt{TE}) power spectra, using the low-$\ell$ \texttt{Commander} and \texttt{simall} likelihoods ($\ell<30$), together with the high-$\ell$ \texttt{CamSpec} likelihood ($\ell\geq30$) from the Planck PR4 (NPIPE) release~\footnote{\url{https://github.com/carronj/planck_PR4_lensing.git}} \cite{Efstathiou2021Detailed,rosenberg2022cmb}. In addition, we include the Planck PR4 CMB lensing likelihood \cite{carron2022cmb,carron2022planck} and the ACT DR6 lensing likelihood~\footnote{\url{https://github.com/ACTCollaboration/act_dr6_lenslike.git}} \cite{madhavacheril2024atacama,qu2024atacama}.

\item \textbf{Type Ia Supernova :} Second, we consider three Type~Ia supernova (SNe~Ia) compilations. We first use the Pantheon$+$ sample\footnote{\url{https://github.com/PantheonPlusSH0ES/DataRelease.git}}\cite{scolnic2022pantheon}, comprising 1,701 light curves from 1,550 SNe~Ia over the redshift range $0.001 \leq z \leq 2.26$. To reduce systematic uncertainties arising from peculiar velocities, we exclude supernovae with $z<0.01$, following the standard Pantheon$+$ analysis \cite{brout2022pantheon}. We then include the DES-Dovekie compilation\footnote{\url{https://github.com/des-science/DES-SN5YR.git}} \cite{DES:2025sig}, which consists of 1,820 photometrically calibrated SNe~Ia spanning the redshift interval $0.02 \leq z \leq 1.14$. Finally, we use the Union3 sample\footnote{\url{https://github.com/rubind/union3_release}} \cite{rubin2025union}, containing 2,087 SNe~Ia covering $0.05 < z < 2.26$. Together, these SNe~Ia compilations provide complementary constraints on the late-time expansion history of the Universe.

\item \textbf{Baryon Acoustic Oscillation :} Finally, we use the Baryon Acoustic Oscillation (BAO) measurements from the Dark Energy Spectroscopic Instrument (DESI) Data Release~2 (DR2)\footnote{\url{https://github.com/CobayaSampler/bao_data.git}} \cite{abdul2025desi}. The dataset includes measurements from multiple tracers, namely the Bright Galaxy Sample (BGS), Luminous Red Galaxies (LRG), Emission Line Galaxies (ELG), Quasars (QSO), and the Lyman-$\alpha$ forest, spanning the redshift range $0.3 \leq z \leq 2.33$. The DESI DR2 BAO likelihood is implemented using the compressed distance measurements $D_M/r_d$, $D_H/r_d$, and $D_V/r_d$, where $D_H(z)=c/H(z)$ is the Hubble distance, $D_M(z)=c\int_0^z dz'/H(z')$ is the transverse comoving distance, $D_V(z)=\left[zD_M^2(z)D_H(z)\right]^{1/3}$ is the volume-averaged distance, and $r_d$ denotes the sound horizon at the baryon drag epoch.
\end{itemize}

In addition to constraining the cosmological parameters, we carry out a Bayesian model comparison between the IDMDE and $\Lambda$CDM models. The Bayesian evidence, $\mathcal{Z}$, is evaluated using the {\tt MCEvidence}\footnote{\url{https://github.com/yabebalFantaye/MCEvidence.git}} package \cite{heavens2017marginal,heavens2017no}, accessed through the Cobaya interface implemented in the {\tt wgcosmo} repository\footnote{\url{https://github.com/williamgiare/wgcosmo.git}}. The relative performance of the two models is quantified through the logarithmic Bayes factor, $\ln B
= \ln \mathcal{Z}_{\rm IDMDE} -\ln \mathcal{Z}_{\Lambda{\rm CDM}}.$ The strength of the Bayesian evidence is interpreted according to the revised Jeffreys scale \cite{kass1995bayes,trotta2008bayes}: $\ln B<1$ corresponds to inconclusive evidence, $1\leq\ln B<2.5$ to weak evidence, $2.5\leq\ln B<5$ to moderate evidence, $5\leq\ln B<10$ to strong evidence, and $\ln B\geq10$ to decisive evidence in favor of the $\Lambda$CDM model. Conversely, negative values of $\ln B$ indicate a preference for the $\Lambda$CDM model over IDMDE. The priors chosen for the $\Lambda$CDM and MEDE models are summarized in Table~\ref{tab_1}.

Before constraining the IMEDE model with observational data, we first examine its impact on the Cosmic Microwave Background (CMB) temperature power spectrum and the linear matter power spectrum. These two observables provide complementary information about the evolution of the Universe and help understand how the IMEDE model parameters affect the background expansion and the growth of matter perturbations. Fig~\ref{fig_1} shows the CMB temperature power spectrum (left panels) and the linear matter power spectrum (right panels) for several representative combinations of the IMEDE model parameters $(z_t,\Delta,Q)$. The upper panels display the absolute spectra, while the lower panels show the corresponding fractional differences with respect to the $\Lambda$CDM model.

In the left panels, the CMB temperature power spectrum is only mildly affected by variations in the IMEDE model parameters. As the values of $z_t$, $\Delta$, and $Q$ increase, the low-$\ell$ temperature anisotropies are noticeably enhanced, while the amplitudes of the acoustic peaks are gradually suppressed. In contrast, the positions of the acoustic peaks remain nearly unchanged, indicating that the IMEDE model mainly modifies the amplitude of the temperature fluctuations rather than the acoustic scale. This behavior is also evident in the lower-left panel, where the fractional difference with respect to the $\Lambda$CDM model becomes larger for parameter combinations with larger values of $z_t$, $\Delta$, and $Q$.

The right panels show the corresponding linear matter power spectrum. Unlike the CMB spectrum, the IMEDE model enhances the matter power spectrum over the entire range of scales considered. The enhancement is relatively small on large scales but becomes progressively stronger toward smaller scales (larger $k$), with the largest deviations obtained for larger values of $z_t$, $\Delta$, and $Q$. The fractional differences shown in the lower-right panel clearly demonstrate that the departure from the $\Lambda$CDM prediction increases with both the scale and the magnitude of the IMEDE model parameters.

These results show that the IMEDE model has noticeable effects on both the CMB temperature anisotropies and the linear matter power spectrum. While the CMB temperature spectrum is mainly affected through changes in the amplitude of the acoustic peaks, the matter power spectrum shows a significant enhancement on small scales. Therefore, both observables provide powerful and complementary probes for constraining the IMEDE model parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{
\renewcommand{\arraystretch}{1}
\begin{table}[t] 
    \centering
    \begin{tabular}{|lll|}
    \hline
    parametrization & parameter & prior\\  
    \hline 
    $\rm \Lambda$\textbf{CDM} & $\Omega_\mathrm{cdm}h^2$ & $\mathcal{U}[0.001, 0.99]$ \\   
    & $\Omega_\mathrm{b}h^{2}$ & $\mathcal{U}[0.005, 0.1]$ \\
    & $100\theta_\mathrm{MC}$ & $\mathcal{U}[0.5, 10]$ \\
    & $\ln(10^{10} A_\mathrm{s})$ & $\mathcal{U}[1.61, 3.91]$ \\
    & $n_\mathrm{s}$ & $\mathcal{U}[0.8, 1.2]$ \\
    & $\tau$ & $\mathcal{U}[0.01, 0.8]$ \\
    \hline 
    \textbf{IMEDE} & $z_t$ & $\mathcal{U}[0, 1]$ \\
     & $\Delta$ & $\mathcal{U}[0, 3]$ \\
     & $Q$ & $\mathcal{U}[-3, -0.001]$ \\
    \hline
    \end{tabular}
    \caption{ Parameters and priors used in the analysis. Here $\mathcal{U}[{\rm min, max}]$ denotes a uniform prior over the specified range.}\label{tab_1}
\end{table}
}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{figure}
\centering
\includegraphics[scale=0.32]{fig_3.pdf}
\caption{The figure shows the corner plot of the IMEDE model at 68\% ($1\sigma$) and 95\% ($2\sigma$) confidence levels using DESI DR2 combined with CMB and SNe~Ia datasets (Pantheon$+$, DES-Dovekie, and Union3), shown as superimposed contours for the different dataset combinations.}\label{fig_2}
\end{figure}

\begin{figure*}
\centering
\includegraphics[scale=0.30]{fig_4.pdf}
\caption{The figure shows the BAO distance measurements, $D_M/(r_d\sqrt{z})$, $D_V/(r_d\sqrt{z})$, and $zD_H/(r_d\sqrt{z})$, together with the best-fit predictions of the IMEDE (solid lines) and $\Lambda$CDM (dashed lines) models obtained from DESI DR2 combined with CMB and SNe~Ia datasets (Pantheon$+$, DES-Dovekie, and Union3). The colored points with error bars represent the DESI DR2 BAO observations. The lower panel shows the normalized residuals, $(\mathrm{data}-\mathrm{model})/\sigma$, for both models, while the gray shaded regions correspond to the $\pm1\sigma$ and $\pm2\sigma$ confidence intervals.}\label{fig_3}
\end{figure*}

\begin{figure*}
\centering
\includegraphics[scale=0.42]{fig_5.pdf}
\caption{The figure shows the evolution of $H(z)/(1+z)$ (first row), $q(z)$ (second row), $w(z)$ (third row), and $f_{\mathrm{DE}}(z)$ (fourth row) for the IMEDE and $\Lambda$CDM models as functions of $z$. The solid lines represent the mean predictions, while the light and dark shaded regions correspond to the $1\sigma$ and $2\sigma$ confidence intervals, respectively.}\label{fig_4}
\end{figure*}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{sidewaystable}
\centering
\vspace*{\fill}
\resizebox{\textheight}{!}{%
\begin{tabular}{lcccccccccccccccc}
\hline
\textbf{Parameters} 
& \textbf{CMB + DESI DR2} 
& \textbf{CMB + DESI DR2 + Pantheon$+$} 
& \textbf{CMB + DESI DR2 + DES-Dovekie} 
& \textbf{CMB + DESI DR2 + Union3} \\
\hline

$\Omega_{\mathrm{\text{b}}}h^2$ 
& $0.02229 \pm 0.00013$ 
& $0.02230 \pm 0.00013$
& $0.02230 \pm 0.00013$
& $0.02230 \pm 0.00013$ \\

$\Omega_{\mathrm{\text{cdm}}}h^2$ 
& $0.1156_{-0.0017}^{+0.0033}$ 
& $0.1142_{-0.0019}^{+0.0045}$ 
& $0.1143_{-0.0021}^{+0.0043}$ 
& $0.1152_{-0.0019}^{+0.0035}$  \\

$100\,\theta_{\mathrm{MC}}$ 
& $1.04107 \pm0.00030$ 
& $1.04116_{-0.00036}^{+0.00029}$ 
& $1.04115_{-0.00036}^{+0.00028}$ 
& $1.04110_{-0.00032}^{+0.00028}$  \\

$\tau_{\mathrm{reio}}$ 
& $0.0520 \pm 0.0068$ 
& $0.0544 \pm 0.0071$
& $0.0536 \pm 0.0069$
& $0.0520 \pm 0.0068$ \\

$\ln\bigl(10^{10}\,A_{\mathrm{s}}\bigr)$ 
& $3.034 \pm 0.013$ 
& $3.039 \pm 0.013$
& $3.037 \pm 0.013$
& $3.034 \pm 0.013$ \\

$n_{\mathrm{s}}$ 
& $0.9675 \pm 0.0039$ 
& $0.9678 \pm 0.0039$
& $0.9678 \pm 0.0039$
& $0.9676 \pm 0.0040$ \\

$z_t$ 
& $0.431_{-0.082}^{+0.12}$ 
& $0.316_{-0.12}^{+0.092}$ 
& $0.317_{-0.098}^{+0.079}$ 
& $0.398_{-0.090}^{+0.079}$  \\

$\Delta$ 
& $>1.83$ 
& $1.06_{-0.43}^{+0.37}$ 
& $1.28_{-0.47}^{+0.39}$ 
& $1.84 \pm 0.51$ \\

$Q$ 
& $>-0.112$ 
& $>-0.0980$
& $>-0.111$
& $>-0.116$ \\

\hline
$H_0 (\mathrm{km\,s^{-1}\,Mpc^{-1}})$ 
& $65.2_{-1.7}^{+1.2}$ 
& $67.46 \pm 0.59$
& $67.28 \pm 0.55$ 
& $66.00 \pm 0.81$ \\

$\Omega_{\mathrm{m}}$ 
& $0.326_{-0.015}^{+0.018}$ 
& $0.301_{-0.0071}^{+0.011}$ 
& $0.3032_{-0.0071}^{+0.0096}$ 
& $0.317_{-0.0092}^{+0.010}$  \\

$\sigma_{8}$ 
& $0.829_{-0.036}^{+0.020}$ 
& $0.849_{-0.038}^{+0.016}$ 
& $0.851_{-0.038}^{+0.017}$ 
& $0.839_{-0.036}^{+0.018}$ \\

$S_{8}$ 
& $0.864_{-0.026}^{+0.016}$
& $0.850_{-0.025}^{+0.012}$
& $0.855_{-0.027}^{+0.014}$
& $0.862_{-0.028}^{+0.015}$ \\

$r_d\text{(Mpc)}$ 
& $147.62 \pm 0.23$ 
& $147.65 \pm 0.23$
& $147.65 \pm 0.23$
& $147.62 \pm 0.23$ \\

\hline

$\Delta{\chi^{2}_{\text{MAP}}}$
& -9.30 
& -9.00
& -11.81 
& -14.58
\\

$N\sigma$ 
& 2.23  
& 2.18
& 2.65 
& 3.06
\\

$\ln B_{i,j}$ 
& $7.89$
& $7.12$
& $9.01$
& $10.33$\\

\hline
\end{tabular}
}
\caption{This table shows the constraints on the IMEDE parameters, presenting the mean values and their uncertainties at the 68\% ($1\sigma$) confidence level, obtained using DESI DR2 combined with CMB and SNe~Ia datasets (Pantheon$+$, DES-Dovekie, and Union3).}\label{tab_2}
\end{sidewaystable}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{sidewaystable}
\centering
\vspace*{\fill}
\resizebox{\textheight}{!}{%
\begin{tabular}{lc|ccc|ccc|ccc|ccc}
\hline
\textbf{Distance} &
\textbf{Redshift}
& \multicolumn{3}{c|}{\textbf{CMB + DESI DR2}}
& \multicolumn{3}{c|}{\textbf{CMB + DESI DR2 + Pantheon$+$}}
& \multicolumn{3}{c|}{\textbf{CMB + DESI DR2 + DES-Dovekie}}
& \multicolumn{3}{c|}{\textbf{CMB + DESI DR2 + Union3}}\\
\cline{3-14}
&
&
\textbf{DESI DR2} & \textbf{IMEDE ($\#\sigma$)} & \textbf{$\Lambda$CDM ($\#\sigma$)}
&
\textbf{DESI DR2} & \textbf{IMEDE ($\#\sigma$)} &\textbf{$\Lambda$CDM ($\#\sigma$)}
& \textbf{DESI DR2} & \textbf{IMEDE ($\#\sigma$)} & \textbf{$\Lambda$CDM ($\#\sigma$)}
& \textbf{DESI DR2} & \textbf{IMEDE ($\#\sigma$)} & \textbf{$\Lambda$CDM ($\#\sigma$)}\\
\hline
$D_V(z)/(r_d \sqrt{z})$ & $0.29$ & $14.74 \pm 0.12$ & $+1.30\,\sigma$ & $+1.34\,\sigma$ & $14.74\pm0.12$ & $+1.81\,\sigma$ & $+1.44\,\sigma$ & $14.74\pm0.12$ & $+1.84\,\sigma$ & $+1.34\,\sigma$ & $14.74\pm0.12$ & $+1.77\,\sigma$ & $+1.55\,\sigma$  \\
& $1.48$ & $21.41 \pm 0.32$ & $+0.52\,\sigma$ & $+0.56\,\sigma$ & $21.41 \pm 0.32$ & $+0.48\,\sigma$ & $+0.59\,\sigma$ & $21.41\pm0.32$ & $+0.44\,\sigma$ & $+0.56\,\sigma$ & $21.41\pm0.32$ & $+0.54\,\sigma$ & $+0.61\,\sigma$  \\

\hline
$D_M(z)/(r_d \sqrt{z})$ & $0.51$ & $19.01 \pm 0.22$ & $+1.64\,\sigma$ & $+1.38\,\sigma$ & $19.01 \pm 0.22$ & $+1.73\,\sigma$ & $+1.45\,\sigma$ & $19.01 \pm 0.22$ & $+1.79\,\sigma$ & $+1.38\,\sigma$ & $19.01 \pm 0.22$ & $+1.85\,\sigma$ & $+1.52\,\sigma$  \\
& $0.70$ & $20.73 \pm 0.20$ & $+0.35\,\sigma$ & $-0.37\,\sigma$ & $20.73 \pm 0.20$ & $+0.04\,\sigma$ & $-0.29\,\sigma$ & $20.73 \pm 0.20$ & $+0.14\,\sigma$ & $-0.37\,\sigma$ & $20.73 \pm 0.20$ & $+0.38\,\sigma$ & $-0.21\,\sigma$  \\
& $0.93$ & $22.36 \pm 0.15$ & $-0.002\,\sigma$ & $-1.15\,\sigma$ & $22.36 \pm 0.15$ & $-0.67\,\sigma$ & $-1.05\,\sigma$ & $22.36 \pm 0.15$ & $-0.54\,\sigma$ & $-1.15\,\sigma$ & $22.36 \pm 0.15$ & $-0.10\,\sigma$ & $-0.95\,\sigma$   \\
& $1.32$ & $24.02 \pm 0.26$ & $-0.37\,\sigma$ & $-0.89\,\sigma$ & $24.02 \pm 0.26$ & $-0.73\,\sigma$ & $-0.83\,\sigma$ & $24.02 \pm 0.26$ & $-0.69\,\sigma$ & $-0.89\,\sigma$ & $24.02 \pm 0.26$ & $-0.44\,\sigma$ & $-0.78\,\sigma$ \\
& $2.33$ & $26.19 \pm 0.34$ & $+1.89\,\sigma$ & $+1.84\,\sigma$ & $26.19 \pm 0.34$ & $+1.80\,\sigma$ & $+1.88\,\sigma$ & $26.19 \pm 0.34$ & $+1.78\,\sigma$ & $+1.84\,\sigma$ & $26.19 \pm 0.34$ & $+1.90\,\sigma$ & $+1.90\,\sigma$ \\

\hline
$z D_H(z)/(r_d\sqrt{z})$ & $0.51$ & $15.61 \pm 0.29$ & $-0.46\,\sigma$ & $-1.72\,\sigma$ & $15.61 \pm 0.29$ & $-1.43\,\sigma$ & $-1.69\,\sigma$ & $15.61 \pm 0.29$ & $-1.30\,\sigma$ & $-1.73\,\sigma$ & $15.61 \pm 0.29$ & $-0.86\,\sigma$ & $-1.66\,\sigma$\\
& $0.70$ & $16.27 \pm 0.27$ & $-1.10\,\sigma$ & $-2.11\,\sigma$ & $16.27 \pm 0.27$ & $-1.99\,\sigma$ & $-2.08\,\sigma$ & $16.27 \pm 0.27$ & $-1.91\,\sigma$ & $-2.11\,\sigma$ & $16.27 \pm 0.27$ & $-1.48\,\sigma$ & $-2.06\,\sigma$  \\
& $0.93$ & $17.01 \pm 0.18$ & $+0.90\,\sigma$ & $+0.40\,\sigma$ & $17.01 \pm 0.18$ & $+0.23\,\sigma$ & $+0.44\,\sigma$ & $17.01 \pm 0.18$ & $+0.21\,\sigma$ & $+0.40\,\sigma$ & $17.01 \pm 0.18$ & $+0.61\,\sigma$ & $+0.45\,\sigma$  \\
& $1.32$ & $16.28 \pm 0.25$ & $-0.05\,\sigma$ & $+0.47\,\sigma$ & $16.28 \pm 0.25$ & $+0.11\,\sigma$ & $+0.49\,\sigma$ & $16.28 \pm 0.25$ & $-0.0001\,\sigma$ & $+0.48\,\sigma$ & $16.28 \pm 0.25$ & $+0.02\,\sigma$ & $+0.49\,\sigma$  \\
& $2.33$ & $13.17 \pm 0.15$ & $-1.09\,\sigma$ & $-0.01\,\sigma$ & $13.17 \pm 0.15$ & $-0.52\,\sigma$ & $-0.01\,\sigma$ & $13.17 \pm 0.15$ & $-0.65\,\sigma$ & $-0.01\,\sigma$ & $13.17 \pm 0.15$ & $-0.87\,\sigma$ & $-0.03\,\sigma$  \\
\hline
\end{tabular}
}
\caption{The DESI results (and their 1$\sigma$ errors) are presented for three distinct types of distances investigated by BAO measurements. For each data point, we indicate the consistency between the best-fit predictions of the IDE and $\Lambda$CDM models and the observed data, expressed in units of observational uncertainties ($\# \sigma$).}\label{tab_3}
\end{sidewaystable}

\section{Results}\label{sec_4}
Fig~\ref{fig_2} shows the marginalized constraints on the cosmological parameters of the IMEDE model obtained from the combined analysis of DESI DR2, CMB, and Type Ia supernova data, including the Pantheon$+$, DES-Dovekie, and Union3 catalogs. The off-diagonal panels present the two-dimensional (2D) posterior contours, showing the correlations between pairs of model parameters. The darker and lighter shaded regions correspond to the 68\% ($1\sigma$) and 95\% ($2\sigma$) confidence levels, respectively. The diagonal panels shows the 1D marginalized posterior distributions for each parameter. The corresponding numerical constraints, given as the mean values with their 68\% ($1\sigma$) uncertainties, are listed in Table~\ref{tab_2} for the different combinations of DESI DR2, CMB, and Type Ia supernova datasets.

We first examine the status of the $H_0$ tension in the IMEDE model. It is important to emphasize that this model belongs to the class of interacting dynamical dark energy models, where the interaction occurs only at late times, after the recombination epoch. Consequently, the early-Universe physics remains unchanged, and hence the sound horizon at the drag epoch, $r_d$, also remains unchanged \cite{Bernal:2016gxb,Aylor:2018drw,Knox:2019rjx,Efstathiou:2021ocp,Jiang:2024xnu}. This expectation is confirmed by our parameter constraints, which yield $r_d = 147.62 \pm 0.23$, $147.65 \pm 0.23$, $147.65 \pm 0.23$, and $147.62 \pm 0.23$ Mpc for the CMB + DESI DR2, CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 dataset combinations, respectively.

This point is particularly important because any successful solution of the $H_0$ tension must also account for the value of the sound horizon at the drag epoch, $r_d$. These two quantities are anti-correlated and can be regarded as two sides of the same coin. In general, if a cosmological model predicts a higher value of $H_0$, closer to the local measurement reported by the SH0ES collaboration, $(73.04 \pm 1.04),\mathrm{km,s^{-1},Mpc^{-1}}$~\cite{riess2022comprehensive}, it typically requires a corresponding reduction in $r_d$ of approximately $7\%$~\cite{knox2020hubble}.

Indeed, this behavior can also be observed in the IMEDE model. Despite the different values of $H_0$ inferred from the various dataset combinations, the corresponding values of $r_d$ remain essentially unchanged. This can be clearly seen in Table~\ref{tab_2}, where the constraints on $H_0$ (Row~10) and $r_d$ (Row~14) show that the sound horizon at the drag epoch is nearly identical for all dataset combinations. Therefore, any dynamical dark energy model that aims to alleviate the $H_0$ tension while leaving $r_d$ essentially unchanged is unlikely to provide a complete solution to the $H_0$ tension~\cite{vagnozzi2023seven,jedamzik2021reducing}.

Yet, several attempts have been made to resolve the $H_0$ tension by modifying the expansion history of the Universe around $z \sim 10^4$ to $z \sim 10^3$~\cite{poulin2019early}. However, such modifications typically require a higher physical matter density, as is the case in Early Dark Energy (EDE) models. Although EDE provides a good fit to the BAO + CMB observations and can reduce the $H_0$ tension, it does not provide a satisfactory fit to the SNe~Ia datasets and also worsens the $S_8$ tension~\cite{chaussidon2025early}.

We now discuss the main reason why the IMEDE model is unable to address the $H_0$ tension. As mentioned earlier, the IMEDE model belongs to the class of interacting dynamical dark energy models and modifies only the late-time expansion history of the Universe. In particular, the IMEDE model prefers an equation of state with $w(z)>-1$ at late times, a result that will be discussed in more detail in the following sections. As shown in ~\cite{Vagnozzi:2018jhn,Vagnozzi:2019ezj,Alestas:2020mvb,Banerjee:2020xcn,Lee:2022cyh,Colgain:2025nzf}, changing the dark energy equation of state from $w(z)=-1$ ($\Lambda$CDM) to $w(z)>-1$ (quintessence) generally leads to a lower inferred value of $H_0$, whereas a phantom-like equation of state, $w(z)<-1$, typically favors a higher value of $H_0$. For this reason, the IMEDE model is unable to resolve the $H_0$ tension. 

This behavior is clearly reflected in our results, which give $H_0 = 65.2^{+1.2}_{-1.7}$, $67.46 \pm 0.59$, $67.28 \pm 0.55$, and $66.00 \pm 0.81~\mathrm{km\,s^{-1}\,Mpc^{-1}}$ for the CMB + DESI DR2, CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 dataset combinations, respectively. Quantitatively, the corresponding tensions with the SH0ES measurement are approximately $4.39\sigma$, $4.67\sigma$, $4.90\sigma$, and $6.11\sigma$, respectively. None of these values is consistent with the local SH0ES measurement, indicating that the IMEDE model does not alleviate the $H_0$ tension.

We also investigate the parameter $S_8 \equiv \sigma_8 \left(\Omega_m/0.3\right)^{0.5}$, which quantifies the amplitude of matter density fluctuations. The $S_8$ tension is another important challenge in cosmology. Recent weak-lensing surveys, including the Dark Energy Survey (DES)~Y6 ($S_8 = 0.794^{+0.009}_{-0.012}$)~\cite{DES:2026fyc}, the Kilo-Degree Survey (KiDS)-Legacy (Shear) ($S_8 = 0.815^{+0.016}_{-0.021}$)~\cite{wright2025kids}, and several others (see Table~1 of Ref.~\cite{pantos2026status}), consistently favor lower values of $S_8$ than the combined Planck + ACT + SPT constraint, $S_8 = 0.836^{+0.012}_{-0.013}$~\cite{DES:2026fyc}.

In the IMEDE model, we obtain $S_8 = 0.892^{+0.032}_{-0.059}$, $0.834^{+0.010}_{-0.019}$, $0.8312^{+0.0090}_{-0.016}$, and $0.839^{+0.012}_{-0.023}$ for the CMB + DESI DR2, CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 dataset combinations, respectively. Comparing these results with the KiDS-Legacy (Shear) weak-lensing measurement, the corresponding tensions are $1.75\sigma$, $1.34\sigma$, $1.45\sigma$, and $1.66\sigma$, respectively. Relative to the DES~Y6 weak-lensing constraint, the corresponding tensions are $2.98\sigma$, $2.63\sigma$, $2.65\sigma$, and $2.84\sigma$, respectively. Although the inclusion of Type Ia supernova data reduces the $S_8$ tension compared to the CMB + DESI DR2 dataset alone, the IMEDE model still predicts values of $S_8$ that remain systematically higher than those preferred by current weak-lensing observations, particularly the DES~Y6 results.

We now examine the constraints on the interaction parameter, $Q$, which determines both the strength and the direction of the energy exchange between dark matter and dark energy. In our model, negative values of $Q$ ($\mathcal{Q}<0$) correspond to an energy transfer from dark matter to dark energy ($\mathrm{DM}\rightarrow\mathrm{DE}$), whereas positive values ($\mathcal{Q}>0$) describe the opposite energy flow ($\mathrm{DE}\rightarrow\mathrm{DM}$). The standard noninteracting $\Lambda$CDM cosmology is recovered for $Q=0$. Using the CMB+DESI DR2 dataset, we obtain only a lower bound on the interaction parameter, $Q>-0.112$ (95\% C.L.), indicating no statistically significant preference for a nonzero interaction. The inclusion of the Pantheon$+$, DES-Dovekie, and Union3 Type Ia supernova compilations yields similarly tight lower bounds, $Q>-0.0980$, $Q>-0.111$, and $Q>-0.116$, respectively. The consistent constraints obtained from the different Type Ia supernova compilations indicate that the current cosmological observations do not provide evidence for an interaction between dark matter and dark energy.

We now examine the constraints on the transition redshift, $z_t$, obtained from the different dataset combinations. We find $z_t=0.431^{+0.12}_{-0.082}$, $0.316^{+0.092}_{-0.12}$, $0.317^{+0.079}_{-0.098}$, and $0.398^{+0.079}_{-0.090}$ for the CMB + DESI DR2, CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 dataset combinations, respectively. Since $w(z_t)=-1$, these results indicate that the phantom-divide crossing occurs at relatively low redshift, $z_t\simeq0.3$-$0.4$, during the late-time accelerated expansion of the Universe.

We also constrain the transition parameter $\Delta$, obtaining $\Delta>1.83$, $1.06^{+0.37}_{-0.43}$, $1.28^{+0.39}_{-0.47}$, and $1.84\pm0.51$ for each dataset combination, respectively. The positive values of $\Delta$ indicate an evolution of the dark energy equation of state from a phantom-like phase ($w<-1$) at $z>z_t$ to a quintessence-like phase ($w>-1$) at $z<z_t$. While the CMB + DESI DR2 dataset provides only a lower bound on $\Delta$, the addition of the Type Ia supernova datasets yields tighter constraints on this parameter.

To assess the statistical preference for the IMEDE model relative to $\Lambda$CDM, we compare their best-fit $\chi^2$ values obtained using \textsc{iminuit}~\cite{dembinski2020scikit}, defining $\Delta\chi^2_{\rm MAP} = \chi^2_{\rm IMEDE} - \chi^2_{\Lambda{\rm CDM}}.$
The IMEDE model extends $\Lambda$CDM through three additional parameters, $Q$, $z_t$, and $\Delta$. Under the assumptions of Gaussian likelihoods and a valid null hypothesis ($\Lambda$CDM), Wilks' theorem \cite{wilks1938large} implies that $\Delta\chi^2_{\rm MAP}$ follows a $\chi^2$ distribution with one degree of freedom. The corresponding Gaussian significance, $N\sigma$, is then used to quantify the preference for the interacting dark sector model. We classify the preference as inconclusive ($N\sigma<1$), weak ($1\le N\sigma<2.5$), moderate ($2.5\le N\sigma<5$), strong ($5\le N\sigma<10$), or decisive ($N\sigma\ge10$). A value of $N\sigma=0$ indicates that there is no statistical preference for the IMEDE model over $\Lambda$CDM, as the improvement in the best-fit $\chi^2$ is either absent or too small to be statistically significant.

The CMB + DESI DR2 and CMB + DESI DR2 + Pantheon$+$ dataset combinations yield $\Delta\chi^2_{\rm MAP}=-9.30$ and $-9.00$, corresponding to Gaussian significances of $N\sigma=2.23$ and $2.18$, respectively. According to our adopted classification, these results provide weak statistical evidence in favor of the IMEDE model over $\Lambda$CDM. For the CMB + DESI DR2 + DES-Dovekie and CMB + DESI DR2 + Union3 dataset combinations, we obtain $\Delta\chi^2_{\rm MAP}=-11.81$ and $-14.58$, corresponding to $N\sigma=2.65$ and $3.06$, respectively. These results indicate a moderate statistical preference for the IMEDE model over the standard $\Lambda$CDM cosmology.

The logarithmic Bayes factors are $\ln \mathcal{B}_{{\rm IMEDE},\Lambda{\rm CDM}}=7.89$, $7.12$, $9.01$, and $10.33$ for the CMB + DESI DR2, CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 dataset combinations, respectively. According to the revised Jeffreys' scale, the first three dataset combinations provide strong evidence in favor of the IMEDE model over $\Lambda$CDM, while the CMB + DESI DR2 + Union3 combination yields decisive evidence in favor of IMEDE model.

To further investigate the origin of the Bayesian preference for IDMDE, we compare the best-fit predictions of the IDMDE and $\Lambda$CDM models with the DESI DR2 BAO distance measurements summarized in Table~\ref{tab_3}, while the corresponding residuals are shown in the upper panels of Fig.~\ref{fig_3}. A comparison across all four dataset combinations shows a stable pattern of residuals, indicating that the DESI BAO measurements largely drive the distance constraints, while the inclusion of Pantheon$+$, DES-Dovekie, or Union3 produces only minor shifts in the best-fit predictions.

Overall, both cosmological models provide a good description of the DESI DR2 BAO measurements. However, compared with $\Lambda$CDM, the IMEDE model systematically reduces the residuals for several BAO measurements, particularly those involving the radial distance $zD_H(z)/(r_d\sqrt{z})$ and the transverse comoving distance $D_M(z)/(r_d\sqrt{z})$. This improvement is observed consistently for all four dataset combinations and is not driven by a single BAO measurement.

The largest improvement is found for the radial distance measurement $zD_H(z)/(r_d\sqrt{z})$ at $z=0.70$. For the CMB + DESI DR2 dataset, the residual decreases from approximately $2.11\sigma$ in $\Lambda$CDM to $1.10\sigma$ in the IMEDE model. Similar improvements are found for the CMB + DESI DR2 + Pantheon$+$, CMB + DESI DR2 + DES-Dovekie, and CMB + DESI DR2 + Union3 combinations, where the residuals decrease from about $2.08$--$2.11\sigma$ in $\Lambda$CDM to approximately $1.48$--$1.99\sigma$ in the IMEDE model. Although this measurement remains the most discrepant BAO point, the tension is consistently reduced in the IMEDE scenario.

A similarly significant improvement is obtained for the radial distance measurement at $z=0.51$. For the CMB + DESI DR2 combination, the residual decreases from $1.72\sigma$ in $\Lambda$CDM to only $0.46\sigma$ in the IMEDE model. Even after including the different Type Ia supernova compilations, the IMEDE residuals remain systematically smaller, ranging between $0.86\sigma$ and $1.43\sigma$, whereas the corresponding $\Lambda$CDM residuals remain close to $1.7\sigma$. Likewise, the transverse comoving distance measurements at $z=0.93$ and $z=1.32$ are reproduced substantially better by the IMEDE model, with residuals generally below $0.7\sigma$, compared with values approaching or exceeding $1\sigma$ for $\Lambda$CDM.

Several other BAO measurements are also reproduced slightly better by the IMEDE model. For example, the low-redshift volume-averaged distance $D_V(z)/(r_d\sqrt{z})$ at $z=0.29$ shows a smaller residual for the CMB + DESI DR2 dataset, while the higher-redshift measurements of $D_V(z)/(r_d\sqrt{z})$ at $z=1.48$ and $D_M(z)/(r_d\sqrt{z})$ at $z=2.33$ are described at nearly the same level by both cosmological models. Overall, most of the remaining BAO measurements are reproduced within approximately the $1\sigma$ level by the IMEDE model. These results indicate that the statistical preference for the IMEDE model is not produced by a single outlying BAO measurement. Instead, it originates from a coherent reduction of the residuals across several independent DESI DR2 distance measurements spanning the full redshift range.

In Fig~\ref{fig_4}, we show the reconstructed evolution of $H(z)/(1+z)$ (first row), $q(z)$ (second row), $w(z)$ (third row), and $f_{\mathrm{DE}}(z)$ (fourth row) as functions of redshift for the IMEDE and $\Lambda$CDM models. In each panel, the solid curves denote the mean reconstructed evolution, whereas the darker and lighter shaded regions correspond to the $1\sigma$ and $2\sigma$ confidence intervals, respectively. The first row compares the evolution of $H(z)/(1+z)$ predicted by the IMEDE and $\Lambda$CDM models. For the CMB + DESI DR2 dataset, the IMEDE model predicts values that are slightly lower than those of the $\Lambda$CDM model over the intermediate redshift range. After including the Pantheon$+$, DES-Dovekie, and Union3 supernova samples, the predictions of the two models become similar, with only small residual differences that are most clearly visible in the zoom-in panels. The second row shows the evolution of the deceleration parameter, $q(z)$, for the IMEDE and $\Lambda$CDM models. Both models show the transition from a decelerating ($q>0$) to an accelerating ($q<0$) Universe. Compared to $\Lambda$CDM, the IMEDE model predicts a slightly less negative deceleration parameter at low redshift, indicating a marginally weaker present-day acceleration.

The third row presents the reconstructed evolution of the dark energy equation of state, $w(z)$, for the IMEDE model, while the $\Lambda$CDM model corresponds to the constant value $w=-1$. For all four dataset combinations, the reconstructed equation of state exhibits a transition from the phantom regime ($w<-1$) at low redshift to the quintessence regime ($w>-1$) at higher redshift, crossing the phantom divide, $w=-1$, at $z\sim0.3$--$0.5$. This behavior is characteristic of a Quintom-B dark energy scenario \cite{Cai:2025mas}. The fourth row shows the corresponding evolution of the normalized dark energy density, $f_{\mathrm{DE}}(z)$. While the $\Lambda$CDM model predicts a constant value, $f_{\mathrm{DE}}(z)=1$, the IMEDE model satisfies $f_{\mathrm{DE}}(0)=1$ at the present epoch for all dataset combinations. As the redshift increases, $f_{\mathrm{DE}}(z)$ initially shows a mild enhancement above unity before decreasing below the $\Lambda$CDM prediction, reflecting the dynamical evolution of dark energy in the IMEDE model.

\section{Conclusions}\label{sec_5}
In this work, we have presented a comprehensive observational study of the Interacting Metastable Emergent Dark Energy model using the DESI DR2 baryon acoustic oscillation measurements combined with Planck PR4 CMB observations and three different Type Ia supernova catalog. Our primary objective was to investigate whether current cosmological observations support an interacting metastable emergent dark energy scenario and to examine the preference of the interaction between dark matter and dark energy.

Our statistical analysis consistently favors the IMEDE model over the standard $\Lambda$CDM cosmology for all dataset combinations considered. Based on the Gaussian significance, the CMB + DESI DR2 and CMB + DESI DR2 + Pantheon$+$ combinations provide weak statistical evidence in favor of the IMEDE model, while the CMB + DESI DR2 + DES-Dovekie and CMB + DESI DR2 + Union3 combinations yield a moderate statistical preference. Bayesian model comparison further strengthens this result, providing strong evidence for the IMEDE model for the first three dataset combinations and decisive evidence for the CMB + DESI DR2 + Union3 combination.

The statistical preference for the IMEDE model is not driven by a single anomalous BAO measurement. Instead, it originates from a coherent reduction of the residuals across several independent DESI DR2 distance measurements spanning the full observed redshift range. In particular, the largest improvements occur for the radial BAO measurements at $z=0.51$ and $z=0.70$, together with noticeable improvements in several transverse distance measurements. These results indicate that the IMEDE model provides a more consistent description of the DESI DR2 distance measurements than the standard $\Lambda$CDM cosmology.

The reconstructed expansion history also reveals several distinctive features of the IMEDE scenario. The dark energy equation of state undergoes a transition from a phantom phase at earlier times to a quintessence phase at recent epochs, crossing the phantom divide at $z\simeq0.3$--$0.5$. This Quintom-B behavior naturally emerges from the Interacting Metastable Emergent Dark Energy framework and is accompanied by a smooth evolution of the dark energy density without introducing any significant modification to the early-Universe physics.

Despite these encouraging results, the IMEDE model does not provide a complete solution to either the $H_0$ or the $S_8$ tensions. Since the IMEDE modifies only the late-time evolution of the Universe, the sound horizon at the drag epoch remains essentially unchanged. Consequently, although the inferred values of the Hubble constant vary among the different dataset combinations, they remain significantly lower than the SH0ES determination. Likewise, the predicted values of $S_8$ remain systematically higher than those preferred by current weak-lensing surveys, particularly DES~Y6.

Overall, our results demonstrate that the IMEDE model provides a statistically well-motivated extension of the standard cosmological model and offers a good description of the DESI DR2 observations. The emergence of a Quintom-B dark energy evolution together with the coherent improvement across multiple independent BAO measurements suggests that the observed preference is unlikely to arise from a single statistical fluctuation. If confirmed by forthcoming DESI DR3 observations and future surveys such as Euclid, the Vera C. Rubin Observatory, and the Nancy Grace Roman Space Telescope, these results could provide important evidence for a dynamical and interacting dark sector beyond the standard $\Lambda$CDM paradigm.

\section*{Acknowledgements}
SC acknowledges the Istituto Nazionale di Fisica Nucleare (INFN) Sez. di Napoli,  Iniziative Specifiche QGSKY and MoonLight-2  and the Istituto Nazionale di Alta Matematica (INdAM), gruppo GNFM, for the support. This paper is based upon work from COST Action CA21136 -- Addressing observational tensions in cosmology with systematics and fundamental physics (CosmoVerse), supported by COST (European Cooperation in Science and Technology). This research was funded by the National Natural Science Foundation of China (NSFC) under Grant No. U2541210 and 12505238.

\section*{Data and Code Availability Statements}
All data used in this study are publicly available from the corresponding survey collaborations and repositories cited in the main text. The computational codes and software employed in the analysis are publicly available, and their corresponding GitHub repositories or official source links are provided in the relevant citations throughout the manuscript. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\bibliography{sample}

\end{document}


