%Version 3.1 December 2024
% See section 11 of the User Manual for version history
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                                 %%
%% Please do not use \input{...} to include other tex files.       %%
%% Submit your LaTeX manuscript as one .tex document.              %%
%%                                                                 %%
%% All additional figures and files should be attached             %%
%% separately and not embedded in the \TeX\ document itself.       %%
%%                                                                 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%\documentclass[referee,sn-basic]{sn-jnl}% referee option is meant for double line spacing

%%=======================================================%%
%% to print line numbers in the margin use lineno option %%
%%=======================================================%%

%%\documentclass[lineno,pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%=========================================================================================%%
%% the documentclass is set to pdflatex as default. You can delete it if not appropriate.  %%
%%=========================================================================================%%

%%\documentclass[sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style

%%Note: the following reference styles support Namedate and Numbered referencing. By default the style follows the most common style. To switch between the options you can add or remove “Numbered” in the optional parenthesis. 
%%The option is available for: sn-basic.bst, sn-chicago.bst%  
 
%%\documentclass[pdflatex,sn-nature]{sn-jnl}% Style for submissions to Nature Portfolio journals
%%\documentclass[pdflatex,sn-basic]{sn-jnl}% Basic Springer Nature Reference Style/Chemistry Reference Style
\documentclass[pdflatex,sn-mathphys-num]{sn-jnl}% Math and Physical Sciences Numbered Reference Style
%%\documentclass[pdflatex,sn-mathphys-ay]{sn-jnl}% Math and Physical Sciences Author Year Reference Style
%%\documentclass[pdflatex,sn-aps]{sn-jnl}% American Physical Society (APS) Reference Style
%%\documentclass[pdflatex,sn-vancouver-num]{sn-jnl}% Vancouver Numbered Reference Style
%%\documentclass[pdflatex,sn-vancouver-ay]{sn-jnl}% Vancouver Author Year Reference Style
%%\documentclass[pdflatex,sn-apa]{sn-jnl}% APA Reference Style
%%\documentclass[pdflatex,sn-chicago]{sn-jnl}% Chicago-based Humanities Reference Style

%%%% Standard Packages
%%<additional latex packages if required can be included here>

\usepackage{graphicx}%
\usepackage{multirow}%
\usepackage{amsmath,amssymb,amsfonts}%
\usepackage{amsthm}%
\usepackage{mathrsfs}%
\usepackage[title]{appendix}%
\usepackage{xcolor}%
\usepackage{textcomp}%
\usepackage{manyfoot}%
\usepackage{booktabs}%
\usepackage{algorithm}%
\usepackage{algorithmicx}%
\usepackage{algpseudocode}%
\usepackage{listings}%
%%%%

%%%%%=============================================================================%%%%
%%%%  Remarks: This template is provided to aid authors with the preparation
%%%%  of original research articles intended for submission to journals published 
%%%%  by Springer Nature. The guidance has been prepared in partnership with 
%%%%  production teams to conform to Springer Nature technical requirements. 
%%%%  Editorial and presentation requirements differ among journal portfolios and 
%%%%  research disciplines. You may find sections in this template are irrelevant 
%%%%  to your work and are empowered to omit any such section if allowed by the 
%%%%  journal you intend to submit to. The submission guidelines and policies 
%%%%  of the journal take precedence. A detailed User Manual is available in the 
%%%%  template package for technical guidance.
%%%%%=============================================================================%%%%

%% as per the requirement new theorem styles can be included as shown below
\theoremstyle{thmstyleone}%
\newtheorem{theorem}{Theorem}%  meant for continuous numbers
%%\newtheorem{theorem}{Theorem}[section]% meant for sectionwise numbers
%% optional argument [theorem] produces theorem numbering sequence instead of independent numbers for Proposition
\newtheorem{proposition}[theorem]{Proposition}% 
%%\newtheorem{proposition}{Proposition}% to get separate numbers for theorem and proposition etc.

\theoremstyle{thmstyletwo}%
\newtheorem{example}{Example}%
\newtheorem{remark}{Remark}%

\theoremstyle{thmstylethree}%
\newtheorem{definition}{Definition}%

\raggedbottom
%%\unnumbered% uncomment this for unnumbered level heads

\begin{document}

\title[Article Title]{Numerical and Statistical Insights into $f(R,T)$ Cosmology: GRP, RK4, and MLE Approaches}

%%=============================================================%%
%% GivenName	-> \fnm{Joergen W.}
%% Particle	-> \spfx{van der} -> surname prefix
%% FamilyName	-> \sur{Ploeg}
%% Suffix	-> \sfx{IV}
%% \author*[1,2]{\fnm{Joergen W.} \spfx{van der} \sur{Ploeg} 
%%  \sfx{IV}}\email{iauthor@gmail.com}
%%=============================================================%%

\author[1]{\fnm{Navya} \sur{Jain}}\email{jainnavya008@gmail.com}

\author*[2]{\fnm{R.K.} \sur{Mishra}}\email{ravkmishra@yahoo.co.in}
%\equalcont{These authors contributed equally to this work.}

\affil[1,2]{\orgdiv{Department of Mathematics}, \orgname{Sant Longowal Institute of Engineering \& Technology}, \orgaddress{\street{Longowal}, \city{City}, \postcode{148106}, \state{Punjab}, \country{India}}}

%%==================================%%
%% Sample for unstructured abstract %%
%%==================================%%

\abstract{This study explores the stability of an $f(R,T)$ gravity-based cosmological model using perturbation techniques and numerical methods. Stability conditions are examined by analyzing the growth rate of small perturbations in density and pressure. Furthermore, the Runge-Kutta fourth-order (RK4) method is employed to numerically track the evolution of these perturbations and validate the theoretical predictions. To ensure consistency with observational data, the model is tested against Hubble parameter measurements and the Pantheon Type Ia supernova dataset using the Maximum Likelihood Estimation (MLE) method. The findings provide deeper insights into the stability of modified gravity models and underscore the significance of combining analytical, numerical, and observational approaches in cosmological studies.}

%%================================%%
%% Sample for structured abstract %%
%%================================%%

% \abstract{\textbf{Purpose:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Methods:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Results:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.
% 
% \textbf{Conclusion:} The abstract serves both as a general introduction to the topic and as a brief, non-technical summary of the main results and their implications. The abstract must not include subheadings (unless expressly permitted in the journal's Instructions to Authors), equations or citations. As a guide the abstract should not exceed 200 words. Most journals do not set a hard limit however authors are advised to check the author instructions for the journal they are submitting to.}

\keywords{Stability, Perturbation theory, Growth rate parameter, $f(R,T)$ gravity, Runge-Kutta method}

%%\pacs[JEL Classification]{D8, H51}

%%\pacs[MSC Classification]{35A01, 65L10, 65L12, 65L20, 65L70}

\maketitle

\section{Introduction}\label{sec1}

Cosmology aims to understand the large-scale structure, evolution, and fundamental forces governing the universe \cite{1,2,3}. Over the years, General Relativity (GR) has served as the cornerstone of modern cosmology, successfully explaining gravitational interactions and the universe's expansion \cite{4}. The standard $\Lambda$CDM model, based on GR, provides a comprehensive framework for describing cosmic dynamics, incorporating cold dark matter (CDM) and a cosmological constant ($\Lambda$) to account for the observed accelerated expansion \cite{5,6}. However, despite its success, the model faces significant theoretical challenges, such as the fine-tuning and coincidence problems, motivating the search for alternative gravitational theories.\\
To address these issues, modified theories of gravity have been proposed, among which $f(R,T)$ gravity has gained considerable attention \cite{7}. This theory extends GR by incorporating dependencies on both the Ricci scalar $R$ and the trace of the energy-momentum tensor $T$, allowing for possible interactions between matter and geometry. Such modifications provide a broader theoretical framework that can potentially explain cosmic acceleration without invoking dark energy explicitly \cite{8,mishra,9,10,11,12,13a,14a,10b}. However, any viable cosmological model must not only fit observational data but also exhibit stability against small perturbations. Stability analysis ensures that minor deviations in cosmic parameters do not lead to unphysical behavior over time, making it a crucial aspect of theoretical cosmology \cite{lifshitz}.\\
Perturbation techniques play a key role in analyzing the stability of cosmological models by examining small fluctuations in density, and pressure\cite{104,105}. In this study, we employ a unique approach by incorporating the growth rate parameters, which provide deeper insights into the evolution of perturbations. Unlike previous works, which primarily focus on linear perturbations \cite{nariai1969gravitational}, our analysis offers a refined perspective on the stability conditions in the context of $f(R,T)$ gravity. Analytical techniques provide a fundamental understanding, serving as a powerful tool for validating theoretical predictions. However, to complement this we employ a well-established fourth order Runge-Kutta (RK4) method, enabling precise numerical tracking of perturbation evolution over time. We believe that this approach not only provide deeper insights into the nonlinear behaviour of perturbation but also provides an independent validation of stability conditions derived from analytical framework \cite{siagian2024comparative}. While extensively employed across numerous disciplines, the utilization of this approach within the current cosmological framework sheds new light on the evolution of perturbations and their influence on the stability of \( f(R,T) \) gravity models.\\
In addition to theoretical development, it is crucial that any viable cosmological model be confronted with empirical data to assess its physical relevance. In this work, we validate our theoretical framework using observational datasets such as the Hubble parameter measurements \cite{moresco2022unveiling} and the Pantheon compilation of Type Ia supernovae \cite{brout2022pantheon+}. For this purpose, we apply the Maximum Likelihood Estimation (MLE) technique, a robust statistical method that identifies optimal model parameters describing the universe's expansion. This data-driven approach helps establish a vital connection between theoretical constructs and observational reality, offering deeper insights into cosmic evolution.\\
To build our analysis, we assume a spatially homogeneous and isotropic cosmological model governed by the Friedmann-Lema\^itre-Robertson-Walker (FLRW) line element:
\begin{equation}\label{eq1}
	ds^{2} = c^2 dt^{2} - a^2(t)(dx^2 + dy^2 + dz^2)
\end{equation}
where \( a(t) \) denotes the scale factor, \( t \) is the cosmic time, and \( (x, y, z) \) represent the co-moving spatial coordinates. The dynamical behavior of the universe?s expansion is described through two key parameters: the Hubble parameter \( H(t) \) and the deceleration parameter \( q(t) \), mathematically defined as
\begin{equation}  
	H(t) = \frac{\dot{a}}{a}, \quad \text{and} \quad q(t) = -\frac{a\ddot{a}}{\dot{a}^2}
\end{equation}
This study follows a systematic approach to analyzing the stability of an $f(R,T)$ cosmological model. First, we derive the field equations governing the model and establish the background evolution of the universe. Next, we introduce perturbations in density and pressure to examine their impact on stability. This approach helps determine whether small perturbations grow uncontrollably or remain bounded over time. To further strengthen our analysis, we employ the RK4 numerical method to verify the stability of our solutions. Finally, we compare our theoretical findings with observational data using the MLE method, ensuring the consistency of our model with empirical evidence. Through this comprehensive investigation, we aim to provide deeper insights into the stability of modified gravity models and their relevance to modern cosmology.
\section{Mathematical Framework}
The \( f(R,T) \) gravity framework generalizes GR by altering its foundational aspects to better tackle certain cosmological and astrophysical challenges. While GR successfully describes gravity as the geometric curvature of spacetime caused by mass-energy, it faces difficulties in explaining phenomena like the universe's accelerated expansion and reconciling with quantum theories. To overcome these shortcomings, Harko et al. introduced the \( f(R,T) \) model in 2011 \cite{7}, which modifies the Einstein-Hilbert action by incorporating arbitrary functions depending on the Ricci scalar \( R \) and the trace of the energy-momentum tensor \( T \). This explicit dependence on \( T \) enables the theory to capture matter-geometry couplings more thoroughly, offering novel insights into how matter influences spacetime dynamics. Within this context, the action is expressed as
\begin{equation}\label{eq2} S = \int \left( \frac{f(R,T)}{16\pi} + \mathcal{L}_m \right) \sqrt{-g}d^4x
\end{equation}
where $g$ is the determinant of the metric tensor $g_{jk}$, $\mathcal{L}_m$ denotes the matter Lagrangian density, $f(R,T)$ is an arbitrary function of $R$ and $T$, and the term $\int\mathcal{L}_m \sqrt{-g}d^{4}x$ replaces the Ricci scalar $R$ of GR, introducing additional degrees of freedom. The stress-energy tensor $T_{jk}$, $T_{jk} = -\frac{2}{\sqrt{-g}} \frac{\delta (\sqrt{-g} \mathcal{L}m)}{\delta g^{jk}},$ incorporates the effects of matter and energy on the gravitational field. Assuming $\mathcal{L}_m$ depends only on $g_{jk}$ and not its derivatives, the stress-energy tensor simplifies as $T{jk} = g_{jk} - 2 \frac{\partial \mathcal{L}_m}{\partial g^{jk}}$. Varying the action $S$ with respect to $g_{jk}$ yields the field equations
\begin{multline}\label{eq4} 
	(8\pi - f_T(R,T)) T_{jk} + \frac{1}{2} f(R,T) g_{jk} = f_R(R,T) R_{jk} + f_T(R,T) \varTheta_{jk} \\ + (g_{jk} \Box - \nabla_j \nabla_k) f_R(R,T)
\end{multline}
where $\nabla_k$ is the covariant derivative, and $\Box$ is the D'Alembertian operator with
$\Box = \frac{1}{c^2} \frac{\partial^2}{\partial t^2} - \bigtriangleup, \quad \text{with} \quad \bigtriangleup = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} + \frac{\partial^2}{\partial z^2}$. The additional term $\varTheta_{jk}$, defined as 
$\varTheta_{jk} = -2 \left[ T_{jk} + g^{\mu\nu} \frac{\partial^2 \mathcal{L}m}{\partial g^{jk} \partial g^{\mu\nu}} \right] + g{jk} \mathcal{L}m,$ depends on the choice of the matter Lagrangian $\mathcal{L}_m$. Adopting a perfect fluid model where $\mathcal{L}_m=-p$, $\varTheta_{jk}$ becomes
$\varTheta{jk} = -2(T_{jk} + \frac{1}{2} p g_{jk})$. Substituting $\varTheta_{jk}$ into Equation \ref{eq4}, the field equations take the form
\begin{multline}\label{eq5} G_{jk} f_R(R,T) = 8\pi T_{jk} + (p + T_{jk}) f_T(R,T) + \frac{1}{2}(f(R,T) - R f_R(R,T)) \\
	- (g_{jk} \Box - \nabla_j \nabla_k) f_R(R,T)
\end{multline}
To facilitate analytical tractability, we adopt the linear form \( f(R,T) = R + 2f(T) \), as introduced by Harko et al., which effectively captures the interaction between the Ricci scalar \( R \) and the trace \( T \) of the energy-momentum tensor. This particular formulation enhances our understanding of gravitational behavior in scenarios involving cosmic acceleration. \\
Employing the FLRW line element together with Einstein's equations, we construct a model to capture the evolution of a universe that is uniform in space and isotropic in all directions. The energy-momentum tensor \( T_{jk} \), representing matter-energy distribution, is modeled as a perfect fluid characterized by density \( \rho \) and pressure \( p \). Under this assumption, it is written as \( T_{jk} = -p g_{jk} + (p + \rho) u_j u_k \), with \( u_j \) being the four-velocity vector satisfying the condition \( g_{jk} u^j u^k = 1 \). Assuming a comoving frame, where the fluid is stationary relative to the coordinates, we have \( u^j = (0, 0, 0, 1) \). For the particular functional form \( f(R,T) = R + 2\gamma T \), the corresponding field equations in \( f(R,T) \) theory, as per equation~\ref{eq5}, reduce significantly when the spacetime geometry is governed by FLRW symmetry and the source is a perfect fluid. The tensor components in this framework are: \( T^1_1 = T^2_2 = T^3_3 = -p \) and \( T^4_4 = \rho \). These lead directly to the modified gravitational field equations under this theoretical setup.
\begin{equation}\label{eq6}
	2\frac{\ddot{a}}{a}+\left( \frac{\dot{a}}{a}\right) ^2=\gamma\rho-(1+3\gamma)p\quad{and}\quad 	3\left( \frac{\dot{a}}{a}\right) ^{2}=(1+3\gamma)\rho-\gamma p
\end{equation}
\section{Stability Analysis}
In cosmology, stability analysis through perturbation techniques is an indispensable tool for assessing the robustness and physical viability of theoretical models, particularly in modified gravity frameworks such as $f(R,T)$ gravity. These techniques help evaluate how small deviations in the metric and matter fields evolve over time, determining the conditions under which a cosmological model remains stable over cosmic timescales \cite{lifshitz,nariai1969gravitational} The introduction of modifications to GR through $f(R,T)$ gravity leads to complex dynamical behaviors, necessitating a thorough perturbative analysis to understand their impact on the evolution of the universe.\\
By linearizing the field equations around a homogeneous and isotropic FLRW background, we derive perturbation equations that enable the identification of potential instabilities. If small perturbations grow unbounded, the model may exhibit unphysical behavior, such as singularities or runaway solutions, making stability analysis a crucial step in validating theoretical predictions. Beyond theoretical interest, stability analysis imposes constraints on viable cosmological models, ensuring consistency with observational data from cosmic microwave background (CMB) measurements \cite{1}, large-scale structure surveys, and other astrophysical observations\cite{2}.\\
In this section, we apply perturbation techniques to analyze the stability of the universe in the $f(R,T)$ framework, focusing on density and pressure fluctuations. Our approach utilizes two methods: the growth rate analysis method, which analytically examines the behavior of perturbations through characteristic equations, and the Runge-Kutta fourth-order method, which numerically solves the perturbation equations to track the stability evolution over time. Since the FLRW metric provides a well-defined background representing the large-scale structure of the universe, it serves as a fundamental setting to study deviations from homogeneity and isotropy.
From equations \ref{eq1} and \ref{eq6}, we derive the continuity equation incorporating $\rho, p,~\text{and}~\frac{\dot{a}}{a}$ as follows:
\begin{equation}\label{eq17}
	(1+3\gamma)\dot{\rho}-\gamma \dot{p}=-3(1+2\gamma)\frac{\dot{a}}{a}(\rho+p).
\end{equation} 
\subsection{Growth Rate Method for Stability Analysis}
The Growth Rate method involves solving this equation to determine the growth rate parameters, which indicate whether perturbations grow, decay, or oscillate over time. This is crucial in determining how the system responds to disturbances and whether it remains stable. As we know, the continuity equation encapsulates the principle of energy conservation in an expanding universe, ensuring that energy density and pressure evolve in accordance with conservation laws. Now differentiating equation \ref{eq17}, we obtain:
\begin{equation}\label{eq18}
	(1+3\gamma)\ddot{\rho}-\gamma \ddot{p}=-3(1+2\gamma)\frac{\dot{a}}{a}(\dot{\rho}+\dot{p})+\frac{3}{2}(1+2\gamma)^{2}(\rho+p)^{2}.
\end{equation}
To analyze the stability of perturbations, we introduce small fluctuations around equilibrium values
\begin{equation}\label{eq19}
	\rho=\rho_{0}+\delta\rho
\end{equation}
\begin{equation}\label{eq}
	p=p_{0}+\delta p
\end{equation}
Here, $\rho_{0}$ and $p_{0}$ represent the equilibrium energy density and pressure, while $\delta\rho$ and $\delta p$ denote small perturbations. Substituting these perturbed quantities into equation \ref{eq18}, we derive a second-order differential equation governing the evolution of perturbations. \\
Assuming the perturbations take an exponential form, i.e., $\delta\rho = \delta\rho_{0}e^{\lambda t}, \quad \delta p = \delta p_{0}e^{\lambda t},$  where $\lambda$ is a growth rate parameter (GRP) and substituting these into the continuity equations, we linearize the system to obtain a set of equations that can be expressed in matrix form as $A\delta = 0$, where $A$ is the perturbation matrix and $\delta$ is the vector of perturbations. Using equations \ref{eq18}-\ref{eq}, we obtain:
\begin{align}
	\begin{bmatrix} 
		(1+3\gamma)\lambda^{2}+X & 0\\ 
		0 & -\gamma\lambda^{2}+X 
	\end{bmatrix}
	\begin{bmatrix}
		\delta\rho \\ \delta p 
	\end{bmatrix} = 
	\begin{bmatrix}
		0 \\ 0
	\end{bmatrix}
\end{align}
where,  
\begin{equation}
	X=3(1+2\gamma)\frac{\dot{a}}{a}\lambda-3(1+2\gamma)^2(\rho_{0}+p_{0}).
\end{equation}
To determine the GRP $\lambda$, we solve the characteristic equation by setting the determinant of the coefficient matrix to zero. This leads to two separate equations:
\begin{equation}\label{eql1}
	(1+3\gamma)\lambda^{2}+3(1+2\gamma)\frac{\dot{a}}{a}\lambda-3(1+2\gamma)^2(\rho_{0}+p_{0})=0,
\end{equation}
\begin{equation}\label{eql2}
	-\gamma\lambda^{2}+3(1+2\gamma)\frac{\dot{a}}{a}\lambda-3(1+2\gamma)^2(\rho_{0}+p_{0})=0.
\end{equation}
By analyzing the signs and magnitudes of the $\lambda$, we can assess the stability of the universe. The conclusions are as follows:
\begin{itemize}
	\item If $\lambda > 0$, at least one mode exhibits exponentially growing perturbations, indicating an unstable system.
	\item If $\lambda = 0$, the system is marginally stable; perturbations neither grow nor decay exponentially, but they may evolve linearly or remain constant.
	\item If $\lambda < 0$, perturbations decay exponentially over time, implying a stable universe.
\end{itemize}
This analysis categorizes the stability of the universe under small perturbations and provides insights into the timescales over which stability or instability manifests, which is critical for predicting the universe's long-term evolution. Since the deceleration parameter (DP) is the most significant parameter in understanding the dynamic evolution of the universe, it plays a crucial role in distinguishing between different phases of cosmic expansion one marked by deceleration (\( q > 0 \)) and the other by acceleration (\( q < 0 \)). Observational data from Type Ia supernovae provide strong evidence for a pivotal cosmic transition \cite{13,14,15}. Berman \cite{berman1988cosmological} studied cosmological models with a constant DP, while Akarsu and Dereli \cite{akarsu2012cosmological} introduced a special linear varying law for DP. More recently, Mishra et al. \cite{mishra2024certain,mishra2024fractional,r} proposed a fractional form of DP. Motivated by the characteristic signature-flipping nature of DP and supported by various studies, we adopt the following form:  
\begin{equation}\label{eqq}
	q(t)=-1+\frac{2\beta(H-1)}{H}
\end{equation}
where \( H \) is the Hubble parameter and \( \beta \) is an arbitrary constant. \\
This choice of \( q(t) \) allows us to obtain a well-defined expression for the scale factor \( a(t) \). Since our continuity equation inherently involves \( a(t) \) and its derivatives, its explicit functional form is crucial for solving the system analytically. Without assuming a specific dependence for \( a(t) \), the system remains underdetermined, making it difficult to extract meaningful physical insights. By selecting a well-motivated form of \( q(t) \), we can systematically derive a consistent evolution of \( a(t) \) and, consequently, the cosmic parameters of interest. Starting from equation~\ref{eqq}, we can derive the expression for the Hubble parameter \( H(t) \) as:
\begin{equation}\label{eqh}
	H(t) = \frac{\dot{a}}{a} = \left( \frac{e^t}{a} \right)^{2\beta}
\end{equation}
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{DP}
	\caption{Temporal evolution of the deceleration parameter (DP), Hubble parameter, and scale factor \( a(t) \) with respect to cosmic time \( t \).}
	\label{Fig1}
\end{figure}
Here, the integration constant is taken to be zero for simplicity. Integrating equation~\ref{eqh}, we obtain the following form of the scale factor:
\begin{equation}\label{eq20}
	a(t) = \left(e^{2\beta t} - 1\right)^{\frac{1}{2\beta}}
\end{equation}
Figure~\ref{Fig1} presents the graphical behavior of the deceleration parameter, Hubble parameter, and scale factor as functions of cosmic time \( t \). As time advances toward large values, the deceleration parameter \( q \) approaches \( -1 \), which is consistent with predictions made by the flat \(\Lambda\)CDM cosmological model. This standard model of cosmology, built on the Big Bang framework, suggests that a substantial fraction of the universe is composed of dark matter.\\
The graph of the Hubble parameter \( H \) demonstrates a sharp decline after the inflationary epoch and eventually trends toward very small values, a pattern observed for various values of the parameter \( \beta \). Furthermore, \( a(t) \) grows unbounded with time, indicating that the universe continues to expand indefinitely. This accelerated expansion is typically attributed to the presence of DE, a mysterious component that dominates the late-time dynamics of the cosmos. By substituting the expression for \( \frac{\dot{a}}{a} \) into the field equations~\ref{eql1}--\ref{eql2} and solving for the eigenvalue \( \lambda \), we derive the following set of results:
\begin{equation}
	\lambda_{1}^{\pm}=\frac{3(1+2\gamma)e^{2\beta t}}{2(1+3\gamma)(e^{2\beta t}-1)}\bigg[-1\pm\sqrt{1+\frac{4}{3}(1+3\gamma)(\rho_{0}+p_{0})\bigg(\frac{e^{2\beta t}-1}{e^{2\beta t}}\bigg)^2}\bigg]
\end{equation}
\begin{equation}
	\lambda_{2}^{\pm}=\frac{3(1+2\gamma)e^{2\beta t}}{2\gamma(e^{2\beta t}-1)}\bigg[1\mp\sqrt{1-\frac{4}{3}\gamma(\rho_{0}+p_{0})\bigg(\frac{e^{2\beta t}-1}{e^{2\beta t}}\bigg)^2}\bigg]
\end{equation}  
Here, \( \lambda_{1}^{\pm} \) and \( \lambda_{2}^{\pm} \) represent the corresponding GRPs.  
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{ME1+}
	\caption{Visualization of the growth rate parameter \( \lambda_{1}^{+} \) with respect to cosmic time \( 0<t<20 \) for a constant \( \gamma \) in the range \( -0.3<\gamma<0.1 \).}
	\label{FigE1+}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{ME1-}
	\caption{Visualization of the growth rate parameter \( \lambda_{1}^{-} \) with respect to cosmic time \( 0<t<20 \) for a constant \( \gamma \) in the range \( -0.3<\gamma<0.1 \).}
	\label{FigE1-}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{ME2+}
	\caption{Visualization of the growth rate parameter \( \lambda_{2}^{+} \) with respect to cosmic time \( 0<t<20 \) for a constant \( \gamma \) in the range \( -0.3<\gamma<0.1 \).}
	\label{FigE2+}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{ME2-}
	\caption{Visualization of the growth rate parameter \( \lambda_{2}^{-} \) with respect to cosmic time \( 0<t<20 \) for a constant \( \gamma \) in the range \( -0.3<\gamma<0.1 \).}
	\label{FigE2-}
\end{figure}
Since \( \rho_{0} \) and \( p_{0} \) are equilibrium values, we consider \( p_{0}=-1.1\rho_{0} \) for a cosmological constant-dominated universe. Figures \ref{FigE1+}-\ref{FigE2-} demonstrate that all the GRPs are negative, indicating stable modes where perturbations decay exponentially over time, bringing the system back toward equilibrium. This suggests that small inhomogeneities do not lead to unbounded growth or collapse, thereby maintaining the system's overall stability. This stability is crucial for understanding the universe's expansion. A stable cosmological model implies that the expansion driven by the cosmological constant or dark energy will continue smoothly without leading to catastrophic instabilities. Thus, within the \( f(R,T) \) framework, our model remains stable and demonstrates resilience under the influence of exponential perturbations. \\ 
In the next subsection, we will employ a numerical technique, namely the Runge-Kutta method, to further analyze the stability of the constructed cosmological model.
\subsection{Runge-Kutta Method for Stability Analysis}
While growth rate analysis provides an analytical understanding of stability, numerical methods such as the Runge-Kutta fourth-order (RK4) method offer a computational approach to track the evolution of perturbations. The RK4 method is particularly useful for solving nonlinear differential equations arising in cosmological perturbations, where analytical solutions may be challenging to obtain. Applying the RK4 method, we numerically integrate the perturbation equations derived from equation \ref{eq18}, tracking the behavior of $\delta\rho$ and $\delta p$ over time. This approach allows us to validate the stability conditions inferred from growth rate analysis by observing whether perturbations remain bounded or diverge. The numerical results provide insights into the dynamical evolution of the system, complementing the theoretical predictions obtained through growth rate analysis.\\
The combined use of growth rate and numerical methods offers a comprehensive stability assessment, ensuring that the perturbation behavior aligns with physically viable cosmological scenarios. In the next section, we present the results of our stability analysis, comparing the outcomes of both methods and discussing their implications for the stability of the universe in the $f(R,T)$ gravity framework.
The RK4 method approximates the solution of an ODE of the form $\frac{dy}{dt} = f(t, y)$ by iteratively updating the solution with intermediate slope evaluations, thereby minimizing truncation errors. Given an initial value $y_n$ at time $t_n$, the next value $y_{n+1}$ is computed as:
\begin{equation}
	y_{n+1} = y_n + \frac{h}{6} (k_1 + 2k_2 + 2k_3 + k_4)
\end{equation}
where the intermediate slopes are defined as:
\begin{align}
	k_1 &= f(t_n, y_n), \nonumber\\
	k_2 &= f\left(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_1 \right), \nonumber\\
	k_3 &= f\left(t_n + \frac{h}{2}, y_n + \frac{h}{2} k_2 \right), \nonumber\\
	k_4 &= f(t_n + h, y_n + h k_3). \nonumber
\end{align}
Here, $h$ represents the step size, and the weighted averaging of the four intermediate slopes ensures both numerical stability and precision. In the context of cosmological perturbation analysis, the RK4 method is employed to solve the coupled differential equations governing the evolution of density and pressure perturbations. By integrating these equations numerically, we can rigorously examine the behavior of perturbations over cosmic time and determine their stability under varying initial conditions and cosmological parameters. \\
This numerical approach complements the growth rate analysis presented in the preceding subsection, offering a more comprehensive investigation of stability. While growth rate analysis provides insights into the asymptotic behavior of perturbations, the RK4 method enables a direct visualization of perturbation evolution, capturing transient effects and nonlinearities that may not be evident in a purely linearized framework. Unlike the growth rate method, which inherently encapsulates perturbation effects without requiring an explicit equation of state, the RK4 approach necessitates an explicit assumption, such as $p = w\rho$.
Using the relation $\frac{\dot{a}}{a} = \left(\frac{e^t}{a}\right)^{2\beta}$ and introducing a small perturbation in density, $\rho = \rho_0 + \delta\rho$, the continuity equation transforms into
\begin{equation}
	(1+3\gamma-\gamma w)\dot{\delta\rho} = -3(1+2\gamma)(1+w)\left(\frac{e^t}{a}\right)^{2\beta}(\rho_0 + \delta\rho)
\end{equation}
For numerical implementation, we assume initial conditions $\delta\rho(0)=1$ and $a(0)=1$, ensuring a well-defined starting point for the integration process. Using the above equation the analytical solution for $\delta\rho$ is obtained as 
\begin{equation}
	\delta\rho=exp\bigg(\frac{-3(1+2\gamma)(1+w)}{(1+3\gamma-\gamma w)2\beta}ln(e^{2\beta t}+a_{0}^{2\beta}-1)\bigg)-\rho_0
\end{equation}
The numerical results obtained via the RK4 method are compared with analytical solutions, as illustrated in Figures \ref{Fig3} and \ref{Fig2}. 
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{RKN1}
	\caption{Density perturbation with respect to cosmic time $(t)$}
	\label{Fig3}
\end{figure}
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{RKN2}
	\caption{Scale factor evolution with respect to cosmic time $(t)$}
	\label{Fig2}
\end{figure}
Table \ref{table1} presents a detailed comparison between the numerical and analytical values of the scale factor $a(t)$ and the density perturbation $\delta\rho(t)$.
\begin{table}[h]
	\setlength{\tabcolsep}{3.5pt}
	\caption{Comparison of Numerical and Analytical values of $a(t)$ and $\delta\rho(t)$}
	\begin{tabular}{|c|c|c|c|c|c|c|c|c|} 
		\hline
		Time & Numerical $a$ & Analytical $a$ & Numerical $\rho$ & Analytical $\rho$ & Error $a$ & Error $\rho$ & {\begin{tabular}[c]{@{}c@{}}
				Relative \\ Error $a$ \end{tabular}} & {\begin{tabular}[c]{@{}c@{}}
				Relative \\ Error $\rho$ \end{tabular}}
		\cr\hline
		0.00 & 1.0000 & 1.0000 & 1.0000 & 1.0000 & 0.0000 & 0.0000 & 0.0000 & 0.0000 \cr\hline
		1.05 & 2.8577 & 2.8577 & 0.57821 & 0.57821 & 1.4938e-10 & 8.4432e-12 & 5.2275e-11 & 1.4603e-11 \cr\hline
		2.11 & 8.2482 & 8.2482 & 0.33258 & 0.33258 & 5.7123e-10 & 1.597e-11 & 6.9255e-11 & 4.8019e-11 \cr\hline
		3.16 & 23.571 & 23.571 & 0.1923 & 0.1923 & 1.7597e-09 & 1.6745e-11 & 7.4656e-11 & 8.7077e-11 \cr\hline
		4.21 & 67.357 & 67.357 & 0.11119 & 0.11119 & 5.1462e-09 & 1.4237e-11 & 7.6402e-11 & 1.2804e-10 \cr\hline
		5.26 & 192.48 & 192.48 & 0.06429 & 0.06429 & 1.4815e-08 & 1.0905e-11 & 7.6966e-11 & 1.6962e-10 \cr\hline
		6.32 & 555.57 & 555.57 & 0.03698 & 0.03698 & 4.2862e-08 & 7.8324e-12 & 7.7149e-11 & 2.118e-10 \cr\hline
		7.37 & 1587.6 & 1587.6 & 0.021382 & 0.021382 & 1.2258e-07 & 5.4235e-12 & 7.7207e-11 & 2.5365e-10 \cr\hline
		8.42 & 4536.9 & 4536.9 & 0.012363 & 0.012363 & 3.5038e-07 & 3.6535e-12 & 7.7228e-11 & 2.9551e-10 \cr\hline
		9.47 & 12965 & 12965 & 0.0071484 & 0.0071484 & 1.0013e-06 & 2.4118e-12 & 7.7233e-11 & 3.3739e-10 \cr\hline
		10.53 & 37421 & 37421 & 0.0041117 & 0.0041117 & 2.8902e-06 & 1.5611e-12 & 7.7233e-11 & 3.7966e-10 \cr\hline
		11.58 & 1.0694e+05 & 1.0694e+05 & 0.0023774 & 0.0023774 & 8.2595e-06 & 1.0022e-12 & 7.7236e-11 & 4.2154e-10 \cr\hline
		12.63 & 3.0559e+05 & 3.0559e+05 & 0.0013746 & 0.0013746 & 2.3602e-05 & 6.3702e-13 & 7.7234e-11 & 4.6341e-10 \cr\hline
		13.68 & 8.7327e+05 & 8.7327e+05 & 0.00079482 & 0.00079482 & 6.7448e-05 & 4.0161e-13 & 7.7236e-11 & 5.0529e-10 \cr\hline
		14.74 & 2.5206e+06 & 2.5206e+06 & 0.00045718 & 0.00045718 & 0.00019468 & 2.5033e-13 & 7.7236e-11 & 5.4756e-10 \cr\hline
		15.79 & 7.2029e+06 & 7.2029e+06 & 0.00026434 & 0.00026434 & 0.00055631 & 1.5581e-13 & 7.7234e-11 & 5.8944e-10 \cr\hline
		16.84 & 2.0583e+07 & 2.0583e+07 & 0.00015284 & 0.00015284 & 0.0015898 & 9.6493e-14 & 7.7238e-11 & 6.3132e-10 \cr\hline
		17.89 & 5.882e+07 & 5.882e+07 & 8.8376e-05 & 8.8376e-05 & 0.0045431 & 5.9494e-14 & 7.7238e-11 & 6.7319e-10 \cr\hline
		18.95 & 1.6978e+08 & 1.6978e+08 & 5.0833e-05 & 5.0833e-05 & 0.013113 & 3.637e-14 & 7.7237e-11 & 7.1547e-10 \cr\hline
		20 & 4.8517e+08 & 4.8517e+08 & 2.9392e-05 & 2.9392e-05 & 0.037473 & 2.226e-14 & 7.7237e-11 & 7.5734e-10  \cr\hline
	\end{tabular}
	\label{table1}
\end{table}
From figure \ref{Fig3}, it is evident that $\delta\rho \to 0$ as time progresses. This indicates that energy density perturbations decay over time, implying that small fluctuations in density do not grow uncontrollably but rather dissipate. Such behavior suggests that the perturbations are stable and that the background cosmological model is resistant to perturbative deviations. Additionally, Figure \ref{Fig2} shows that $a(t) \to \infty$, which is consistent with a continuously expanding universe. The combination of $\delta\rho \to 0$ and $a \to \infty$ further suggests that the universe evolves towards a more homogeneous state, where initial density fluctuations smooth out over time. \\
To assess the accuracy of the numerical approach, we evaluate the error between numerical and analytical solutions using the relative error metric. The relative error is defined as
\begin{equation}
	RE_{a}=\frac{E_{a}}{a_{ana}}, \qquad{where}\qquad E_{a}= |a_{num} - a_{ana}|
\end{equation}
\begin{equation}
	RE_{\delta\rho}=\frac{E_{\delta\rho}}{\delta\rho_{ana}}, \qquad{where}\qquad	E_{\delta\rho} = |\delta\rho_{num} - \delta\rho_{ana}|
\end{equation}
Here $RE_{a}$ and $RE_{\delta\rho}$ represent the relative errors in the scale factor and density perturbation, respectively, while $E_{a}$ and $E_{\delta\rho}$ denote the absolute errors in the scale factor and density perturbation. The terms $a_{ana}$ and $\delta\rho_{ana}$ are the respective analytical values, whereas $a_{num}$ and $\delta\rho_{num}$ are the numerical values. Figure \ref{Fige} depicts the evolution of relative and absolute errors with respect to cosmic time $t$. These error metrics quantify the deviation of the numerical solutions from the exact analytical expressions. The relative errors, though increasing as shown in figure, remain extremely small, reinforcing the accuracy of the numerical solution. The computed errors remain within acceptable bounds, demonstrating that the RK4 method accurately approximates the perturbation dynamics. The small errors validate the robustness of the numerical approach in capturing the stability characteristics of the system.
\begin{figure}[h]
	\centering
	\includegraphics[width=.9\linewidth]{RAerror}
	\caption{Error evolution with respect to cosmic time $(t)$}
	\label{Fige}
\end{figure}
Thus, both numerical and analytical methods independently confirm the stability of perturbations, reinforcing the robustness of the underlying cosmological model. The RK4 method, through its ability to track the dynamical evolution of perturbations, provides a crucial numerical validation of the stability conclusions drawn from the GRP analysis.
\section{MLE-Based Analysis of Observational Datasets}\label{sec5}
In this section, we employ the Maximum Likelihood Estimation technique to analyze two key cosmological datasets: the Hubble parameter measurements compiled by Moresco et al.~\cite{moresco2022unveiling} and the Pantheon Type Ia supernova (SNe Ia) dataset from~\cite{brout2022pantheon+}. MLE is a powerful statistical tool used to estimate model parameters by maximizing the probability of the observed data under a given theoretical framework. This method is particularly advantageous over conventional $\chi^2$ minimization, as it better accommodates datasets and  non-linear models with non-Gaussian or complex error structures by directly computing the likelihood function.
\subsection{Model Fitting to Hubble Parameter Data}
Redshift $z$ serves as a critical observational quantity in cosmology, capturing the wavelength shift in light due to cosmic expansion. Within the context of the FLRW metric, the redshift is related to the cosmic scale factor $a(t)$ by:
\begin{equation} \label{eqa}
	1 + z = \frac{a_0}{a(t)}
\end{equation}
Here, $a_0$ is the scale factor's present day value. As $a(t)$ evolves with $t$, the $t-z$ relation can be written as $\frac{dt}{dz} = -\frac{1}{(1+z)H(z)}$, expressing how time changes with redshift. Utilizing Equations~\ref{eq20} and~\ref{eqa}, we obtain $t-z$ relationship as:
\begin{equation}
	t = \frac{1}{2\beta} \ln \left(1 + \frac{e^{2\beta t_{0}} - 1}{(1+z)^{2\beta}} \right),
\end{equation}
where $t_0$ denotes the current age of the cosmos.
The observational dataset, spanning redshift values from $0.07$ to $1.965$, is utilized to examine the redshift-dependent variation of the Hubble parameter, $H(z)$. The expression characterizing $H(z)$ within our framework is formulated as:
\begin{equation} \label{eq24}
	H(z) = 1 + (H_0 - 1)(1 + z)^{2\beta},
\end{equation}
where $H_0$ is the current Hubble constant, and $\beta$ quantifies the rate of cosmic acceleration specific to our scenario. Assuming the data uncertainties follow a Gaussian distribution, the associated negative log-likelihood function is constructed as:
\begin{equation} \label{eq25}
	\ln \mathcal{L}(H_0, \beta) = -\frac{1}{2} \sum_{i=1}^{32} \frac{\left(H_{\text{obs}, i} - H_{\text{model}, i}\right)^2}{\sigma_i^2},
\end{equation}
where $H_{\text{obs}, i}$ corresponds to the observed value at redshift $z_i$, $H_{\text{model}, i}$ is the computed result from our model, and $\sigma_i$ represents the observational error. By minimizing this objective function, we determine the optimal values of $H_0$ and $\beta$ that best reproduce the empirical data, with each deviation weighted by its respective uncertainty.\\
We present the parameter uncertainties at the $95\%$ confidence level, which are calculated from the covariance matrix derived using the second derivative (Hessian) of the log-likelihood function. Parameter fitting was performed in MATLAB, employing the \texttt{fminunc} optimizer that uses a quasi-Newton technique. With an initial guess provided, the algorithm gradually adjusted the values through successive steps until convergence was achieved. The final optimized parameters are:
\[
H_0 = 61.9837\pm5.2344, \quad \beta = 0.51262\pm0.0853, \quad \mathcal{L} = 8.32828.
\]
To evaluate the comparative performance, we also examine the canonical $\Lambda$CDM model, for which the Hubble function is given as:
\begin{equation} \label{eq26}
	H(z) = H_0 \sqrt{\Omega_m (1 + z)^3 + \Omega_\Lambda},
\end{equation}
where $H_0 = 67.4$~km/s/Mpc, $\Omega_m = 0.315$, and $\Omega_\Lambda = 0.685$, aligned with the findings reported by the Planck 2020 collaboration~\cite{aghanim2020planck}. The optimization details for our model, including metrics over successive iterations, are presented in Table~\ref{table2} to illustrate the convergence process.
\begin{table}[h!t]%
	\caption{Optimality analysis}
	\centering
	\begin{tabular}{|c|c|c|c|c|}
		\hline
		Iteration & Func-count & $f(x)$ & Step-size & First-order optimality  \cr
		\hline
		$0$  & $3$ & $39.5722$ & - & $408$ \cr\hline
		$1$  & $9$ & $12.9252$ & $0.000533987$ & $89.2$ \cr\hline
		$2$  & $12$ & $10.558$ & $1$ & $30.2$ \cr\hline
		$3$  & $15$ & $10.2927$ & $1$ & $4.11$ \cr\hline
		$4$  & $18$ & $10.2878$ & $1$ & $0.705$ \cr\hline
		$5$  & $24$ & $10.2856$ & $1$ & $0.713$ \cr\hline
		$6$  & $27$ & $10.2567$ & $1$ & $5.96$ \cr\hline
		$7$  & $30$ & $10.1998$ & $1$ & $12.5$ \cr\hline
		$8$  & $33$ & $10.044$ & $1$ & $23.3$ \cr\hline
		$9$  & $36$ & $9.71174$ & $1$ & $36$ \cr\hline
		$10$  & $39$ & $9.12469$ & $1$ & $43$ \cr\hline
		$11$  & $42$ & $8.56266$ & $1$ & $30$ \cr\hline
		$12$  & $45$ & $8.36014$ & $1$ & $9.91$ \cr\hline
		$13$  & $48$ & $8.32928$ & $1$ & $1.06$ \cr\hline
		$14$  & $51$ & $8.32828$ & $1$ & $0.0194$ \cr\hline
		$15$  & $54$ & $8.32828$ & $1$ & $0.00415$ \cr\hline
		$16$  & $57$ & $8.32828$ & $1$ & $0.000132$ \cr\hline
	\end{tabular}
	\label{table2}
\end{table}
Figure~\ref{Fig:HubbleData} displays a graphical assessment contrasting the observed Hubble dataset with both the developed theoretical framework and the conventional $\Lambda$CDM scenario. The agreement between empirical measurements and our model outcomes is clearly illustrated, as reinforced by the minimal value of the negative log-likelihood function. Additionally, the figure reflects the sustained accuracy of the proposed framework throughout the entire span of redshift, emphasizing its predictive strength and consistency.
\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{hdata1}
	\caption{Graphical representation of the Hubble observations, including error bars, alongside the optimal fit derived through MLE for the suggested model, and the predictions from the $\Lambda$CDM framework.}
	\label{Fig:HubbleData}
\end{figure}
\subsection{Analysis of Pantheon Type Ia Supernova Dataset}
This analysis incorporates the Pantheon+ SHOES Type Ia Supernova dataset, which includes a total of 1701 measurements of the distance modulus $\mu_{\text{obs}}$ over a redshift interval of $0.00122 \leq z \leq 2.26$~\cite{brout2022pantheon+}. This analysis seeks to estimate the parameters of the theoretical framework by aligning predicted distance moduli with those derived from observations. The model-based distance modulus, denoted by $\mu_{\text{model}}$, takes the form:
\begin{equation} \label{eq27}
	\mu_{\text{model}} = m - M = 5 \log_{10} \left( D_{l}(z) \right) + \mu_0,
\end{equation}
Here, $m$ and $M$ refer to the apparent and intrinsic brightness, respectively, while $\mu_0$ represents a calibration-related parameter not directly tied to cosmology. The corresponding dimensionless luminosity distance, symbolized as $D_l(z)$, is determined using:
\begin{equation} \label{eq28}
	D_l(z) = c(1 + z) \int_0^z \frac{dz^*}{H(z^*)},
\end{equation}
where $H(z^*) = \frac{H(z)}{H_0}$ represents the normalized Hubble parameter, and $z^*$ is the integration variable ranging from 0 to $z$.\\
To estimate optimal values for $H_0$ and $\beta$, we apply the Maximum Likelihood Estimation technique. The likelihood function, assuming Gaussian uncertainties in the observed data, is formulated as:
\begin{equation} \label{eq29}
	\ln \mathcal{L}(H_0, \beta) = -\frac{1}{2} \sum_{i=1}^{1701} \frac{\left( \mu_{\text{obs}, i} - \mu_{\text{model}, i} \right)^2}{\sigma_i^2},
\end{equation}
In this expression, $\mu_{\text{obs}, i}$ corresponds to the experimentally measured distance modulus, $\mu_{\text{model}, i}$ is the model-based estimation, and $\sigma_i$ denotes the associated uncertainty for each observation.
The optimization process involves iteratively tuning $H_0$ and $\beta$ to minimize the negative log-likelihood. Table~\ref{table3} summarizes the optimization procedure, documenting key statistics such as the number of function evaluations, values of the objective function, and the norm of each optimization step.
\begin{table}[h]
	\caption{Optimality analysis}
	\centering
	\begin{tabular}{|c|c|c|c|c|}
		\hline
		Iteration & Func-count & $f(x)$ & First-order optimality & Norm of Step \\
		\hline
		0  & 3  & \( 2.054312 \times 10^{3} \)  & \( 4.888 \times 10^{3} \) & - \\\hline
		1  & 12 & \( 9.285880 \times 10^{2} \)  & \( 1.525 \times 10^{3} \) & \( 3.500 \times 10^{-1} \) \\\hline
		2  & 15 & \( 7.928825 \times 10^{2} \)  & \( 1.733 \times 10^{1} \) & \( 3.850 \times 10^{-1} \) \\\hline
		3  & 18 & \( 7.928210 \times 10^{2} \)  & \( 4.679 \) & \( 4.601 \times 10^{-3} \) \\\hline
		4  & 21 & \( 7.928195 \times 10^{2} \)  & \( 2.123 \times 10^{-2} \) & \( 2.112 \times 10^{-2} \) \\\hline
		5  & 24 & \( 7.928188 \times 10^{2} \)  & \( 1.061 \) & \( 3.201 \times 10^{-2} \) \\\hline
		6  & 27 & \( 7.928145 \times 10^{2} \)  & \( 4.082 \) & \( 2.286 \times 10^{-1} \) \\\hline
		7  & 30 & \( 7.927982 \times 10^{2} \)  & \( 1.254 \times 10^{1} \) & \( 1.240 \) \\\hline
		8  & 33 & \( 7.927851 \times 10^{2} \)  & \( 1.416 \times 10^{1} \) & \( 1.247 \) \\\hline
		9  & 36 & \( 7.927726 \times 10^{2} \)  & \( 8.884 \) & \( 1.344 \) \\\hline
		10 & 39 & \( 7.927690 \times 10^{2} \)  & \( 2.230 \) & \( 2.093 \times 10^{-1} \) \\\hline
		11 & 42 & \( 7.927687 \times 10^{2} \)  & \( 2.205 \times 10^{-2} \) & \( 1.076 \times 10^{-1} \) \\\hline
		12 & 45 & \( 7.927687 \times 10^{2} \)  & \( 2.227 \times 10^{-3} \) & \( 1.239 \times 10^{-2} \) \\\hline
		13 & 48 & \( 7.927687 \times 10^{2} \)  & \( 1.555 \times 10^{-3} \) & \( 3.629 \times 10^{-2} \) \\\hline
		14 & 51 & \( 7.927687 \times 10^{2} \)  & \( 4.570 \times 10^{-5} \) & \( 9.152 \times 10^{-3} \) \\\hline
		15 & 54 & \( 7.927687 \times 10^{2} \)  & \( 8.892 \times 10^{-5} \) & \( 3.991 \times 10^{-5} \) \\\hline
		16 & 62 & \( 7.927687 \times 10^{2} \)  & \( 2.850 \times 10^{-5} \) & \( 1.775 \times 10^{-5} \) \\\hline
		17 & 65 & \( 7.927687 \times 10^{2} \)  & \( 6.747 \times 10^{-5} \) & \( 2.136 \times 10^{-5} \) \\\hline
		18 & 70 & \( 7.927687 \times 10^{2} \)  & \( 7.643 \times 10^{-6} \) & \( 1.172 \times 10^{-6} \) \\\hline
		19 & 81 & \( 7.927687 \times 10^{2} \)  & \( 8.902 \times 10^{-5} \) & \( 7.312 \times 10^{-8} \) \\\hline
		20 & 97 & \( 7.927687 \times 10^{2} \)  & \( 5.910 \times 10^{-5} \) & \( 9.896 \times 10^{-10} \) \\
		\hline
	\end{tabular}
	\label{table3}
\end{table}
The final optimized parameters with $95\%$ confidence level were found to be
\[
H_0 = 63.0994\pm2.0, \quad \beta = 0.60202\pm0.1, \quad \mathcal{L} = 792.7687
\] 
The relatively low value of the negative log-likelihood indicates a strong correspondence between the model predictions and the observational datasets, thereby affirming the reliability of the developed cosmological framework. A detailed comparison among the observational data, the proposed model (as obtained through MLE optimization), and the standard $\Lambda$CDM model is depicted in Figure~\ref{Fig}.
\begin{figure}[h!]
	\centering
	\includegraphics[width=\linewidth]{hdata2}
	\caption{Graphical representation of Pantheon+ SHOES Type Ia Supernova observations, including error bars, alongside the optimal fit derived through MLE for the model, and predictions from the $\Lambda$CDM framework.}
	\label{Fig}
\end{figure}
As seen in the figure, the theoretical outcomes of our model closely follow the observed data, providing further confirmation of its capability to accurately capture the universe's expansion behavior. The efficiency of the fit, supported by the low negative log-likelihood value, strengthens the model's credibility. In summary, the application of the MLE technique has proven to be an effective approach for cosmological parameter estimation. The proposed model, fine-tuned via MLE, demonstrates a superior fit to the observational datasets when compared with the standard $\Lambda$CDM framework, underscoring its robustness and effectiveness.
\section{Conclusion}
In this study, we have systematically analyzed the stability of an $f(R,T)$ gravity-based cosmological model using both perturbation techniques and numerical methods. Our approach involved a twofold validation process: theoretical stability analysis via the GRP and numerical verification through the RK4 method. The GRP analysis clearly demonstrates that perturbations decay over time, confirming that small deviations do not grow uncontrollably but instead dissipate, ensuring the stability of the background cosmology. The RK4 method further supports this conclusion, as the results show that energy density perturbations asymptotically tend toward zero while the scale factor continues to increase indefinitely. The consistently negative values of GRP reinforce that the system remains dynamically stable, preventing any unbounded growth or collapse of small inhomogeneities.\\
In addition to the stability analysis, the proposed model was assessed against observational data using the MLE technique, incorporating both Hubble parameter measurements and the Pantheon compilation of Type Ia supernovae (SNe Ia). MLE was employed due to its effectiveness in handling intricate error structures, making it a more reliable alternative to the conventional $\chi^2$ minimization method. The results exhibit a strong consistency between the theoretical predictions and observational data, thereby enhancing the model?s credibility. The observed minor deviations in parameter estimates across different datasets further illustrate the model?s versatility in fitting various cosmological observations.\\
Furthermore, the optimization analysis presented in the tables demonstrates a structured convergence pattern across iterations. Initially, the function values exhibit significant variations, but as iterations progress, they stabilize, indicating an approach toward an optimal solution. The function order optimality values consistently decrease, confirming improved accuracy in subsequent iterations. Additionally, the norm of the step size reduces systematically, further reinforcing convergence stability. The obtained results validate the effectiveness of the employed numerical optimization method, ensuring robustness and accuracy in solving the present problem. The smooth decline in function values and optimality metrics highlights the reliability of the proposed approach in achieving precise and stable outcomes.\\
These findings collectively affirm that our chosen deceleration parameter function provides a stable and observationally consistent framework within $f(R,T)$ gravity for describing cosmic evolution. The complementary agreement between analytical stability analysis, numerical validation, and observational data fitting underscores the reliability of our approach. Future research could enhance this model by incorporating additional observational constraints and extending the analysis to nonlinear perturbations, particularly in higher-dimensional settings, to gain deeper insights into the long-term dynamics of the universe.

\bmhead{Acknowledgements}

The authors gratefully acknowledge the Department of Mathematics, SLIET, Longowal, for offering the necessary infrastructure and support to carry out this research work.

\section*{Declarations}
\bmhead{Funding}
The authors have no relevant financial or non-financial interests to disclose.
\bmhead{Data availability}
The manuscript titled "Numerical and Statistical Insights into $f(R,T)$ Cosmology: GRP, RK4, and MLE Approaches" integrates theoretical derivations with observational data. Relevant observational datasets are properly referenced throughout the manuscript where applicable. All investigations utilize recognized numerical techniques and standard theoretical frameworks. The referenced literature provides access to the datasets, and readers are invited to reproduce the outcomes, provided proper attribution is maintained.

\bmhead{Conflict of interest} The authors declare that no competing interests are associated with this work.



%%=============================================%%
%% For submissions to Nature Portfolio Journals %%
%% please use the heading ``Extended Data''.   %%
%%=============================================%%

%%=============================================================%%
%% Sample for another appendix section			       %%
%%=============================================================%%

%% \section{Example of another appendix section}\label{secA2}%
%% Appendices may be used for helpful, supporting or essential material that would otherwise 
%% clutter, break up or be distracting to the text. Appendices can consist of sections, figures, 
%% tables and equations etc.


%%===========================================================================================%%
%% If you are submitting to one of the Nature Portfolio journals, using the eJP submission   %%
%% system, please include the references within the manuscript file itself. You may do this  %%
%% by copying the reference list from your .bbl file, paste it into the main manuscript .tex %%
%% file, and delete the associated \verb+\bibliography+ commands.                            %%
%%===========================================================================================%%

\bibliography{sn-bibliography}% common bib file
%% if required, the content of .bbl file can be included here once bbl is generated
%%\input sn-article.bbl

\end{document}
