Constraints on the Minimally Extended Varying Speed of Light Model Using Pantheon+ Dataset

1. Introduction

In the Robertson-Walker (RW) metric, the expanding Universe is depicted as progressing from one hypersurface to another, with the scale factor increasing naturally, resulting in the cosmological redshift of various physical quantities, such as mass density, wavelength, and temperature. One can derive this metric from the cosmological principle and Weyl’s postulate. However, the RW model lacks a mechanism to determine cosmological time dilation (TD). The standard model of cosmology (SMC) makes an additional assumption, asserting that the speed of light is constant (c). This assumption arises from the dependence of SMC on general relativity (GR), which assumes that c is invariant. As a result, the cosmological TD between two hypersurfaces at $t=t_1$ and $t=t_2$ vary in proportion to the inverse of the scale factors $a\left(t\right)$ at those specific times. However, if one allows a time-varying speed of light as proposed in this paper, this relationship may not hold anymore. Thus, establishing TD depends on experimental observations. Given the theoretical absence of cosmological TD, the relationship can be considered as a general function $f\left(a\right)$ of the scale factor, allowing the speed of light to be expressed as $c\left(t_1\right)=\frac{f\left(a_2\right)}{f\left(a_1\right)}\frac{a\left(t_1\right)}{a\left(t_2\right)}c\left(t_2\right)$ [1,2,3].

Thus, the speed of light could change over cosmic time or remain constant within an expanding Universe, contingent upon its relationship with cosmological TD. Without explicit laws governing TD, the speed of light in the RW metric could potentially vary with cosmological time, akin to other physical properties such as mass density, temperature, and fundamental constants like the Planck constant [4]. However, to construct a coherent model around this concept, the varying speed of light (VSL) needs to be incorporated into the Einstein field equations (EFEs) and resolved for solutions. Previous studies, notably within the framework of the minimally extended VSL (meVSL) model, have addressed such scenarios [1,2,3]. One can freely select a local value for the speed of light as it merely entails a scaling of length units. As long as this local value remains constant on a given time hypersurface, it satisfies special relativity (SR) to be consistent with local physics laws. Newton’s gravitational constant, G, could potentially vary. To avoid trivial unit rescaling, one must examine the concurrent variation of c, G, and possibly other physical constants [1].

Numerous endeavors have aimed to measure cosmological TD. One approach involves directly observing TD by analyzing the decay time of distant supernova (SN) light curves and spectra [5,6,7,8,9]. Another method entails measuring TD by examining the stretching of peak-to-peak timescales of gamma-ray bursters (GRBs) [10,11,12,13,14,15,16,17,18]. Efforts have also been made to detect the TD effect in the light curves of quasars (QSOs) located at cosmological distances [19,20]. So far, no definitive detection of cosmic TD has been achieved, with conflicting results from different measurements.

The assessment of whether the speed of light alters due to the Universe’s expansion hinges solely on observation. Therefore, it is vital to observe how any changes in the speed of light affect cosmic scales. To achieve this, diverse observational methods have been utilized.

Among these methods, first, there is the cosmic distance duality relation CDDR method. Etherington’s theorem, derived from the geodesic deviation equation, establishes reciprocity between the area distances of galaxies and observers, linked by the redshift factor $(1+z)$ under geometric invariance [21]. This theorem, applicable in spacetimes where photons follow null geodesics, forms the basis for the CDDR. By relating area distances to angular and luminosity distances, the CDDR offers a means to test the validity of the SMC [22,23]. Various tests of the CDDR using astrophysical and cosmological observations have been conducted to constraints on VSL models [24,25,26,27,28,29,30,31,32,33,34]. Our analysis of the meVSL model suggests a potential deviation from the standard CDDR based on current data. However, with different priors for certain cosmological parameters, the current dataset aligns with the SMC, indicating no deviation from the expected CDDR [1,32]. Therefore, acquiring more precise data is essential to thoroughly investigate any deviations from the established CDDR and reaffirm the viability of the meVSL model.

Second, there is the Cosmic Chronometer (CC) method. The CC method involves observing two passively evolving galaxies, typically elliptical galaxies, assumed to have formed at the same cosmic epoch but observed at different redshifts [35]. This approach offers a model-independent means of measuring the Hubble parameter, $H\left(z\right)$, as a function of redshift, derived from spectroscopic surveys with high precision (${\sigma }_z\le 0.001$). The expansion rate of the meVSL model, or the Hubble parameter $H\left(z\right)$, is determined from the differential age evolution of the Universe $\Delta t$ within a given redshift interval ($dz$).

$$\begin{array}{cc}\hfill H\left(z\right)\equiv \frac{\dot{a}}{a}& =-\frac{1}{1+z}\frac{dz}{dt}\approx -\frac{1}{1+z}\frac{\Delta z}{\Delta t}=H{\left(z\right)}^{\left(\mathrm{SMC}\right)}{(1+z)}^{-b/4}\hfill \end{array}$$

$$\begin{array}{c}\hfill =H_{0}E{\left(z\right)}^{\left(\mathrm{SMC}\right)}{(1+z)}^{-b/4},\end{array}$$

(1)

where $E^{\left(\mathrm{SMC}\right)}$ is the normalized Hubble parameter in the SMC model. Various methods exist for measuring $\Delta t$, including predicting its age based on the chemical composition of a stellar population or utilizing spectroscopic observables like the 4000 Åbreak, known to be linearly related to the age of the stellar population [36]. Unlike many cosmological measurements that rely on integrated distances, the CC method determines the expansion rate $H\left(z\right)$ as a function of the redshift–time derivative $dz/dt$, making it a potent tool for testing different cosmological models [37,38,39,40,41,42,43,44,45,46,47,48,49].. This method proves particularly valuable for investigating VSL models [33]. Both minimum ${\chi }^2$-analysis and maximum-likelihood analysis using the most recent CC data have been performed to constrain the parameter b of the meVSL model. However, the precision of the current CC data is insufficient to distinguish between the meVSL model and the SMC [50].

Third, there are cosmological TDs observed from supernovae. The luminosity curve (LC) of a supernova (SN) offers valuable insights into its evolution, aiding in the classification and understanding of its properties. LC analysis helps determine crucial parameters such as peak luminosity, time to peak brightness, and rate of decline, particularly for SNe Ia, essential as standard candles in cosmology. Comparing LCs across distances enables investigation into cosmic expansion and TD, contributing to significant discoveries like accelerated expansion and dark energy. Wilson’s method involves comparing LCs of nearby and distant SNe, revealing TD effects due to light travel time through space [51]. This information is crucial for studying the Universe’s expansion rate and testing cosmological models, involving data collection, mathematical modeling, and comparison of observed TD with theoretical predictions. We derived a TD formula within the meVSL model, analyzing data from 13 high-redshift SNe Ia to determine the exponent b as $b=0.198\pm 0.415$ [9,52]. While less precise than CC, our analysis indicates consistency with both SMC and the meVSL model. Thus, distinguishing between the two based on SNe TD data is challenging.

Fourth, there is the cosmography method. It employs a kinematic description of the Universe’s evolution based on the cosmological principle, emphasizing the dynamics of cosmic expansion. As a model-independent framework, it offers flexibility in managing cosmological parameters, allowing for generalized analysis unconstrained by preconceived models. By focusing on the later stages of cosmic evolution and utilizing Taylor expansions tailored to the observable domain where $z\ll 1$, cosmography imposes constraints on the present-day Universe. We adapt late-time cosmography to incorporate meVSL models [53].

Lastly, in this paper, we aim to discuss the constraint on any evidence of the cosmic variation of the speed of light using Pantheon+ data [54]. The Pantheon compilation comprises a total of 1048 SNe Ia spanning a redshift range of $0.01$ to $2.3$. It includes 365 spectroscopically confirmed SNe Ia from the Pan-STARRS1 (PS1) Medium Deep Survey, combined with a subset of 279 PS1 SNe Ia (with redshifts ranging from $0.03$ to $0.68$) with reliable distance estimates obtained from various sources such as SDSS, SNLS, and HST samples. Cosmological models fitted to minimize the ${\chi }^2$ for flat $\Lambda$CDM and $\omega$CDM models, without accounting for systematic uncertainties on ${\mathsf{\Omega }}_{m0}$, yield values of $0.284\pm 0.012$ and $0.350\pm 0.035$, respectively. The 1-$\sigma$ constraint on $\omega$ for the $\omega$CDM model is $-1.251\pm 0.144$. The Pantheon dataset allows for a precise constraint of approximately 10% on ${\mathsf{\Omega }}_{m0}$ and 12% on $\omega$ for the flat $\omega$CDM model, and about 4% on ${\mathsf{\Omega }}_{m0}$ for the flat $\Lambda$CDM model.

In Section 2, we give a brief overview of the meVSL model, delving into its theoretical foundations and implications for cosmological phenomena. Specifically, we explore the luminosity distance predictions derived from this model, shedding light on its unique characteristics and potential effects for observational data. Following this exposition, in Section 3, we embark on an in-depth analysis aimed at elucidating the constraints imposed on the variation of the speed of light c and the gravitational constant G. Leveraging the comprehensive Pantheon dataset, we scrutinize these parameters across a spectrum of dark energy models, discerning their impacts on cosmological dynamics and structure formation. Finally, in Section , we distill our findings and insights into a comprehensive summary, drawing actionable conclusions and outlining avenues for future research and exploration in the realm of cosmology and fundamental physics.

2. Summary for the meVSL

The conceptualization of the four-dimensional spacetime of a spatially homogeneous and isotropic, expanding universe entails envisioning it as a seamless continuum composed of homogeneous and isotropic spatial hypersurfaces evolving dynamically over time. At the heart of this framework lies the Robertson-Walker (RW) metric, elegantly expressed by

$$\begin{array}{cc}\hfill d{s}^2& =-{\tilde{c}}^{2}d{t}^2+a^2\left(t\right){\gamma }_{ij}d{x}^{i}d{x}^j=-{\tilde{c}}^{2}d{t}^2+a^2\left(t\right)\left(\frac{d{r}^2}{1-k{r}^2}+r^{2}d{\mathsf{\Omega }}^2\right)\hfill \end{array}$$

$$\begin{array}{c}\hfill =-{\tilde{c}}^{2}d{t}^2+a^2\left(t\right)\left[d{\chi }^2+f_k^2\left(\chi \right)d{\mathsf{\Omega }}^2\right]\equiv -{\tilde{c}}^{2}d{t}^2+a{\left(t\right)}^{2}d{l}_{3\mathrm{D}}^2,\end{array}$$

(2)

where $\tilde{c}$ denotes the speed of light and ${\gamma }_{ij}\left(\overrightarrow{x}\right)$ signifies the time-independent spatial metric defining the hypersurface, while $a\left(t\right)$ governs the scale factor dictating the relationship between physical distance and comoving distance. The derivation of redshift involves utilizing the geodesic equation for a light wave, where $d{s}^2=0$ as represented by Equation (2). The consistency of $d{l}_{3\mathrm{D}}$ over time is maintained by exclusively employing comoving coordinates. Building upon this foundation, we arrive at the expression for radial light signals:

$$\begin{array}{cc}\hfill d{l}_{3\mathrm{D}}& =\frac{c\left(t_i\right)d{t}_i}{a\left(t_i\right)}:\frac{c_{1}d{t}_1}{a_1}=\frac{c_{2}d{t}_2}{a_2}\Rightarrow \left\{\begin{array}{cc}{c}_1=c_2=c\hfill & \mathrm{if}\frac{d{t}_1}{a_1}=\frac{d{t}_2}{a_2}\mathrm{SMC}\hfill \\ c_1={\left(\frac{a_1}{a_2}\right)}^{\frac{b}{4}}{c}_2\hfill & \mathrm{if}\frac{d{t}_1}{a_1^{1-\frac{b}{4}}}=\frac{d{t}_2}{a_2^{1-\frac{b}{4}}}\mathrm{meVSL}\hfill \end{array}\right.,\hfill \end{array}$$

(3)

where $d{t}_i=\nu \left(t_i\right)$ denotes the time interval between successive crests of light at $t_i$ (i.e., the inverse of the frequency ${\nu }_i$ at $t_i$), and $f\left(a_i\right)$ represents an arbitrary function of $a\left(t_i\right)$. Thus, the RW metric naturally allows the VSL models if we remove the traditional assumption on the cosmological TD [1,2,3].

Additionally, the introduction of $\chi =D_{\mathrm{c}}$ as the comoving distance as

$$\begin{array}{c}\hfill D_{\mathrm{C}}\left(z\right)\equiv {\int }_0^r\frac{d{r}^{\prime }}{\sqrt{1-k{r}^2}}=\frac{{\tilde{c}}_0}{H_0}{\int }_0^z\frac{d{z}^{\prime }}{E^{\left(\mathrm{GR}\right)}\left(z^{\prime }\right)}\equiv \frac{{\tilde{c}}_0}{H_0}{d}_{\mathrm{C}}\left(z\right),\end{array}$$

(4)

and $f_k\left(\chi \right)=sinh\left(\sqrt{-k}\chi \right)/\sqrt{-k}=D_{\mathrm{M}}$ as the transverse comoving distance

$$\begin{array}{c}\hfill D_{\mathrm{M}}\left(z\right)=D_{\mathrm{M}}^{\left(\mathrm{GR}\right)}\left(z\right)=\left\{\begin{array}{cc}\frac{{\tilde{c}}_0}{H_0}\frac{1}{\sqrt{{\mathsf{\Omega }}_{K0}}}sinh\left(\sqrt{{\mathsf{\Omega }}_{k0}}\frac{H_0}{{\tilde{c}}_0}{D}_{\mathrm{C}}\right)\hfill & {\mathsf{\Omega }}_{k0}>0\hfill \\ D_{\mathrm{C}}\hfill & {\mathsf{\Omega }}_{k0}=0\hfill \\ \frac{{\tilde{c}}_0}{H_0}\frac{1}{\sqrt{|{\mathsf{\Omega }}_{K0}|}}sin\left(\sqrt{|{\mathsf{\Omega }}_{k0}|}\frac{H_0}{{\tilde{c}}_0}{D}_{\mathrm{C}}\right)\hfill & {\mathsf{\Omega }}_{k0}<0\hfill \end{array}\right.,\end{array}$$

(5)

with ${\tilde{c}}_0=2.9979\times {10}^5$ km/s and $H_0=100h$ km/Mpc/s representing the present values of the speed of light and the Hubble parameter, respectively [53].

Within this conceptual framework, the timelike worldlines of constant space delineate the threading, while the spacelike hypersurfaces of constant time define the slicing within the four-dimensional spacetime. Each spacelike threading corresponds to a homogeneous universe, with the slicing being orthogonal to these universes, offering a natural arrangement conducive to the definition of constant physical quantities such as density, temperature, and the speed of light on each spacelike hypersurface. Thus, our choice of coordinates emerges organically, rendering alternative considerations unnecessary. Furthermore, the derivation of the Ricci tensors and Ricci scalar curvature from the provided metric in Eq. (2) [1] further enriches the understanding of the underlying spacetime dynamics and its mathematical representation.

One can define threading as the timelike worldlines of constant space while determining the slicing by the spacelike hypersurface of constant time within the four-dimensional spacetime. Each spacelike threading corresponds to a homogeneous universe, and slicing is orthogonal to these universes. On each spacelike hypersurface, one can establish constant physical quantities such as density, temperature, and the speed of light, making our choice of coordinates inherently natural, thereby obviating the need for alternative considerations. Ricci tensors and Ricci scalar curvature are derived from the provided metric in Eq. (2) [1]. Specifically, the expressions for $R_{00}$ and $R_{ii}$ are given by

$$\begin{array}{cc}\hfill & R_{00}=-\frac{3}{{\tilde{c}}^2}\left(\frac{\ddot{a}}{a}-H^2\frac{dln\tilde{c}}{dlna}\right),R_{ii}=\frac{g_{ii}}{{\tilde{c}}^2}\left(2\frac{{\dot{a}}^2}{a^2}+\frac{\ddot{a}}{a}+2k\frac{{\tilde{c}}^2}{a^2}-H^2\frac{dln\tilde{c}}{dlna}\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & R=\frac{6}{{\tilde{c}}^2}\left(\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}+k\frac{{\tilde{c}}^2}{a^2}-H^2\frac{dln\tilde{c}}{dlna}\right).\hfill \end{array}$$

(7)

In cosmology, one treats matter as a perfect fluid, defined by its total mass density $\rho$ and isotropic pressure P. This density includes both the rest-mass density, measured in the fluid’s rest frame, and the mass content of the internal elastic energy density. Within the framework of GR, the stress-energy tensor describes this perfect fluid, given by

$$\begin{array}{c}\hfill T^{\mu \nu }=\left(\rho +\frac{P}{c^2}\right)U^{\mu }{U}^{\nu }+P{g}^{\mu \nu },\end{array}$$

(8)

where $U^{\mu }$ represents its four-velocity. When the fluid is in motion, a set of fundamental observers is considered comoving with it, characterized by a four-velocity denoted as $U^{\mu }=(c,0,0,0)$. Once one establishes the metric and the stress-energy tensor, the subsequent step involves solving the Einstein Field Equations (EFEs) to elucidate the dynamics of the scale factor in the metric. These equations govern the dynamics of expansion, including the speed and acceleration of the Universe’s expansion as observed between two fundamental observers. Thus, The energy-momentum tensor of the i-component perfect fluid with the equation of state ${\omega }_i$ is given by [1]

$$\begin{array}{c}\hfill T_{\left(i\right)\mu }^{\nu }=\mathrm{diag}\left(-{\rho }_i{\tilde{c}}^2,P_i,P_i,P_i\right),\mathrm{with}{\rho }_i{\tilde{c}}^2={\rho }_{i0}{\tilde{c}}_0^2{a}^{-3(1+{\omega }_i)},\end{array}$$

(9)

where ${\tilde{c}}_0$ is the present value of the speed of light, ${\rho }_{i0}$ is the present value of mass density of the i-component, and we use $a_0=1$. One can derive Friedmann equations including the cosmological constant from Eqs. (6)-(9)

$$\begin{array}{cc}\hfill & \frac{{\dot{a}}^2}{a^2}+k\frac{{\tilde{c}}^2}{a^2}-\frac{\Lambda {\tilde{c}}^2}{3}=\frac{8\pi \tilde{G}}{3}\sum _i{\rho }_i,\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & \frac{{\dot{a}}^2}{a^2}+2\frac{\ddot{a}}{a}+k\frac{{\tilde{c}}^2}{a^2}-\Lambda {\tilde{c}}^2-2H^2\frac{dln\tilde{c}}{dlna}=-8\pi \tilde{G}\sum _i\frac{P_i}{{\tilde{c}}^2}=-8\pi \tilde{G}\sum _i{\omega }_i{\rho }_i,\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \frac{\ddot{a}}{a}=-\frac{4\pi \tilde{G}}{3}\sum _i\left(1+3{\omega }_i\right){\rho }_i+\frac{\Lambda {\tilde{c}}^2}{3}+H^2\frac{dln\tilde{c}}{dlna}.\hfill \end{array}$$

(12)

2.1. Luminosity Distance

The distance modulus, denoted by $\mu =m-M$, represents the discrepancy between the apparent magnitude m (ideally corrected for interstellar absorption effects) and the absolute magnitude M of an astronomical entity. It is linked to the luminosity distance $D_{\mathrm{L}}$ in parsecs through the formula

$$\begin{array}{c}\hfill \mu =5{log}_{10}\left[\frac{D_{\mathrm{L}}}{1\mathrm{Mpc}}\right]+25.\end{array}$$

(13)

This definition proves convenient as the observed brightness of a light source correlates with its distance according to the inverse square law, and brightnesses are typically expressed in magnitudes. Absolute magnitude M denotes the apparent magnitude of an object when viewed from a distance of 10 parsecs. The relationship between magnitudes and flux F is given by

$$\begin{array}{c}\hfill m=-2.5{log}_{10}F\left(D_{\mathrm{L}}\right),M=-2.5{log}_{10}F\left(10\right).\end{array}$$

(14)

The Chevallier-Polarski-Linder (CPL) parametrization, which assumes $\omega$ to be a linear function of the scale factor a, is presented by [55,56]

$$\begin{array}{c}\hfill \omega ={\omega }_0+(1-a){\omega }_a.\end{array}$$

(15)

The expression for E(GR) in Eq. (2) for a flat Universe using the CPL parametrization is

$$\begin{array}{c}\hfill E^{\left(\mathrm{SMC}\right)}=\sqrt{{\mathsf{\Omega }}_{\mathrm{m}0}{a}^{-3}+\left(1-{\mathsf{\Omega }}_{\mathrm{m}0}\right)a^{-3(1+{\omega }_0+{\omega }_a)}{e}^{-3{\omega }_a(1-a)}}.\end{array}$$

(16)

To determine the luminosity distance in the meVSL model, we need to reevaluate its fundamental definition. Here, the observed luminosity $L_0$ detected at the present epoch differs from the absolute luminosity $L_s$ of the source emitted at redshift z. Conservation of flux from the source to the observed point is

$$\begin{array}{c}\hfill \mathcal{F}=\frac{L_s}{4\pi D_{\mathrm{L}}^2\left(z\right)}=\frac{L_0}{4\pi D_{\mathrm{M}}^2\left(z_0\right)}.\end{array}$$

(17)

The absolute luminosity, $L_s\equiv \Delta E_1/\Delta t_1$, represents the ratio of the emitted light energy $\Delta E_1$ to the emission time interval $\Delta t_1$. Similarly, one can denote the observed luminosity as $L_0=\Delta E_0/\Delta t_0$. Consequently, one can rewrite the luminosity distance using Eq. (17) as

$$\begin{array}{c}\hfill D_{\mathrm{L}}^2\left(z\right)=\frac{L_s}{L_0}{D}_{\mathrm{M}}^2\left(z_0\right)=\frac{\Delta E_1}{\Delta E_0}\frac{\Delta t_0}{\Delta t_1}{D}_{\mathrm{M}}^2\left(z_0\right)={\left(1+z\right)}^{2-\frac{b}{4}}{D}_{\mathrm{M}}^2\left(z_0\right),\end{array}$$

(18)

where we employ

$$\begin{array}{cc}\hfill \frac{\Delta E_1}{\Delta E_0}& =\frac{{\tilde{h}}_1{\tilde{\nu }}_1}{{\tilde{h}}_0{\tilde{\nu }}_0}=\frac{{\tilde{\nu }}_1^{\left(\mathrm{GR}\right)}}{{\tilde{\nu }}_0^{\left(\mathrm{GR}\right)}}=(1+z),\frac{\Delta t_0}{\Delta t_1}=\frac{{\tilde{\nu }}_1}{{\tilde{\nu }}_0}=\frac{{\tilde{\nu }}_1^{\left(\mathrm{GR}\right)}{(1+z)}^{-b/4}}{{\tilde{\nu }}_0^{\left(\mathrm{GR}\right)}}={(1+z)}^{1-\frac{b}{4}}.\hfill \end{array}$$

(19)

This relation also holds for the angular diameter distance $D_{\mathrm{A}}$. Consequently, the luminosity distance in the meVSL model is given by

$$\begin{array}{c}\hfill D_{\mathrm{L}}\left(z\right)={\left(1+z\right)}^{1-\frac{b}{8}}{D}_{\mathrm{M}}\left(z\right)={\left(1+z\right)}^{2-\frac{b}{8}}{D}_{\mathrm{A}}\left(z\right).\end{array}$$

(20)

Under this premise, the modification of the absolute magnitude of SNe Ia is expressed as

$$\begin{array}{c}\hfill \mathrm{M}-{\mathrm{M}}_0=-2.5log\left[\frac{L}{L_0}\right]=\frac{5}{4}blog\left[a\right],\end{array}$$

(21)

where the subscript 0 denotes the local value of M. Thus, the distance modulus of meVSL, $\mu \equiv m-M$, is written as:

$$\begin{array}{cc}\hfill \mu \left(z\right)& =m-\mathrm{M}=5{log}_{10}\left[\frac{D_{\mathrm{L}}}{\mathrm{Mpc}}\right]+25,D_{\mathrm{M}}=\frac{{\tilde{c}}_0}{H_0}{d}_{\mathrm{L}}\left(z\right),d_{\mathrm{L}}\left(z\right)\equiv {\int }_0^z\frac{d{z}^{\prime }}{E^{\left(\mathrm{GR}\right)}\left(z^{\prime }\right)}.\hfill \end{array}$$

(22)

The theoretically predicted apparent magnitude $m_{\mathrm{th}}$ is given by

$$\begin{array}{cc}\hfill m_{\mathrm{th}}\left(z\right)& =M+5log\left[\frac{D_{\mathrm{L}}}{1\mathrm{Mpc}}\right]+25+\frac{5}{4}blog\left[a\right]=M+42.3841+5log\left[\frac{d_{\mathrm{L}}}{h{(1+z)}^{b/4}}\right].\hfill \end{array}$$

(23)

2.2. Analysis

The chi-squared (${\chi }^2$) represents a weighted summation of squared deviations, given by

$$\begin{array}{c}\hfill {\chi }^2=\sum _{i,j}\left(m_{i,\mathrm{obs}}-m_{i,\mathrm{th}}\right)C_{ij}^{-1}\left(m_{i,\mathrm{obs}}-m_{j,\mathrm{th}}\right),\end{array}$$

(24)

where $m_{i,\mathrm{obs}}$ signifies the observed apparent magnitude, $m_{i,\mathrm{th}}$ represents the theoretical apparent magnitude of SNe Ia at the redshift $z_i$ as defined in Eq. (22), and $C_{ij}=D_{ij}+C_{\mathrm{sys}}$ denotes the covariance matrix. Here, $D_{ij}={\sigma }_i^2{\delta }_{ij}$ stands for the variance of each observation, and $C_{\mathrm{sys}}$ is a non-diagonal matrix associated with systematic uncertainties. The reduced chi-square statistic defined as chi-square per degree of freedom is used extensively in the goodness of fit testing

$$\begin{array}{c}\hfill {\chi }_{\nu }^2=\frac{{\chi }^2}{\nu },\end{array}$$

(25)

where the degree of freedom $\nu =N-p$, signifies the number of observations N minus the number of fitted parameters p. As a heuristic, when the variance of the measurement error is known a priori, a ${\chi }_{\nu }^2\gg 1$ suggests a substandard model fit. Conversely, a ${\chi }_{\nu }^2>1$ implies that the fit has not adequately captured the data (or that the error variance has been underestimated). Ideally, a ${\chi }_{\nu }^2$ around 1 indicates that the correspondence between observations and estimates is congruent with the error variance.

3. A Bound for the Variation of c

The previous constraint on the temporal variation of the speed of light c was derived from the variation in the radius of a planet [57]. However, this constraint stemmed from the analysis of the time-varying radius of Mercury [58] using a specific model known as the covariant variable speed of light theory proposed by Magueijo [59]. Consequently, this constraint cannot be directly applied within the framework of the meVSL model.

To explore the variation of the speed of light over time, we turn to SNe Ia, which serve as reliable standard candles for probing the cosmic expansion rate in the late Universe. Our investigation focuses on utilizing data from the Pantheon SNe Ia catalog [54]. Specifically, we delve into two primary models: the $\omega$CDM model, which assumes a constant $\omega$, and the CPL dark energy model. Through this analysis, we aim to shed light on the potential temporal evolution of the speed of light and its implications within these cosmological frameworks.

3.1. $\tilde{c}$ for $\omega$CDM

We explore the $\omega$CDM models, specifically those with ${\omega }_a=0$ as specified in Eq. (15). Utilizing a maximum likelihood analysis, we examine various models characterized by varying cosmological parameters. We show the results of these analyses in Table 1, which reveals intriguing insights into the relationships between these parameters.

Within this table, we uncover noteworthy patterns and correlations among cosmological parameters. Notably, when ${\omega }_0$ is held constant, both $M_0$ and h exhibit degeneracy, as do ${\mathsf{\Omega }}_{m0}$ and b. Consequently, fixing ${\mathsf{\Omega }}_{m0}$ results in nearly identical values for $M_0$ and h, with only b values varying. Conversely, when we fix the value of h, only the $M_0$ values change, while ${\mathsf{\Omega }}_{m0}$ and b become irrelevant in this context.

Furthermore, we observe that among the $\omega$CDM models derived from the Pantheon+ data, those exhibiting discernible time variations in the speed of light are specifically characterized by ${\omega }_0=-1$ and $0.30\le {\mathsf{\Omega }}_{m0}$. This observation underscores the intricate interplay between cosmological parameters and their implications for understanding the temporal evolution of fundamental physical constants.

3.1.1. ${\omega }_0=-1$ without Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

We begin our analysis by examining the $\Lambda$CDM model (where $b=0$), the results of which are presented in the first row of Table 1. In this model, the 1-$\sigma$ level values of $M_0$, ${\mathsf{\Omega }}_{m0}$, and h are $(-19.6512,-19.0732)$, $(0.273,0.297)$, and $(0.605,0.791)$, respectively.

Moving on to the meVSL models (where $b\ne 0$), within a 1-$\sigma$ error range, the cosmological parameters for the model with ${\omega }_0=-1$ are as follows: $-19.6605\le M_0\le -19.0851$, $0.188\le {\mathsf{\Omega }}_{m0}\le 0.410$, $0.603\le h\le 0.787$, and $-0.218\le b\le 0.168$. However, compared to the $\Lambda$CDM model, the constraint on ${\mathsf{\Omega }}_{m0}$ in the meVSL model is notably weaker.

In the meVSL framework, the cosmological evolution of the speed of light is described by $\tilde{c}={\tilde{c}}_0{(1+z)}^{-b/4}$. The best-fit value of b is $-0.025$, indicating that the speed of light in the past was greater than its present value. However, the 1-$\sigma$ range of b encompasses both negative and positive values, suggesting that one cannot definitively conclude on the variation of the speed of light in this model.

3.1.2. ${\omega }_0=-1$ with Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

We then hold ${\mathsf{\Omega }}_{m0}$ constant and perform a chi-square test on these models. Notably, the value of h remains unchanged even as we vary ${\mathsf{\Omega }}_{m0}$ across the range of $0.28$ to $0.32$. As ${\mathsf{\Omega }}_{m0}$ increases, both $M_0$ and b decrease. Specifically, the best-fit value of b is positive only for ${\mathsf{\Omega }}_{m0}=0.28$, while for other values within this range, the best-fit values of b are negative. Within the 1-$\sigma$ error range, the values of b exhibit both positive and negative trends for ${\mathsf{\Omega }}_{m0}=0.28$ and $0.29$. However, for $0.30\le {\mathsf{\Omega }}_{m0}\le 0.32$, the 1-$\sigma$ region consistently yields negative values for b, indicating a decrease in the speed of light over cosmic time. These trends are illustrated in Fig. Figure 1. In the left panel of Fig. Figure 1, we depict the cosmic evolution of the best-fit value of $\tilde{c}/{\tilde{c}}_0$ alongside its 1-$\sigma$ errors for the ${\omega }_0=-1$ and ${\mathsf{\Omega }}_{m0}=0.29$ model. Here, the best-fit value of b is $-0.009$, suggesting an increase in the speed of light with increasing z. However, the uncertainty in b allows for both negative and positive values within 1-$\sigma$ error, leading to ambiguity regarding the variation in the speed of light for this model. For ${\mathsf{\Omega }}_{m0}\ge 0.30$, the best-fit value of b is consistently negative within 1-$\sigma$ error, indicating a monotonically decreasing speed of light over cosmic time. In the right panel of Fig. Figure 1, we illustrate the cosmological evolution of $\tilde{c}/{\tilde{c}}_0$ for the model with ${\omega }_0=-1$ and ${\mathsf{\Omega }}_{m0}=0.30$. Here, the best-fit value of b is $-0.027$, corroborating the decrease in the speed of light over time. Furthermore, the ratio of the time variation of the speed of light to its value at the present epoch, expressed as $\dot{\tilde{c}}0/{\tilde{c}}_0=\frac{b}{4}{H}_0$, is constrained within the 1-$\sigma$ ranges of $(-8.76,-0.89)$, $(-11.80,-3.93)$, and $(-14.84,-6.98)$ for $\mathsf{\Omega }\mathrm{m}0=0.30,0.31$, and 0.32, respectively. These constraints represent significant improvements over those provided in [57].

3.1.3. ${\omega }_0=-1$ with Fixing h

We fix the value of h and proceed with a maximum likelihood analysis. In these scenarios, both ${\mathsf{\Omega }}_{\mathrm{m}0}$ ($=0.299\pm 0.111$) and b ($=-0.025\pm 0.193$) remain constant despite varying h from $0.6376$ to $0.74$. As h increases, $M_0$ also increases accordingly. The best-fit value of b remains consistent at $-0.025$ across these models, indicating a decreasing trend of the speed of light over cosmic time. However, within the 1-$\sigma$ error margin, b values exhibit both positive and negative trends. Consequently, it is plausible to conclude that there is no significant variation in the speed of light within these models.

3.1.4. Fixing h

We maintain a constant value for h and conduct a maximum likelihood analysis. In this scenario, as we vary h from $0.6376$ to $0.74$, the best-fit values of $M_0$, ${\omega }_0$, and b all increase. However, the best-fit value of ${\mathsf{\Omega }}_{\mathrm{m}0}$ decreases with increasing h. Additionally, as h increases, so does $M_0$. The best-fit values of b consistently remain negative, indicating a decrease in the speed of light over cosmic time in these models. However, within the 1-$\sigma$ error range, b values can be both positive and negative, suggesting no discernible variation in the speed of light for these models.

3.1.5. Fixing ${\omega }_0$

We conduct a maximum likelihood analysis for $\omega$CDM models without fixing any cosmological parameters except ${\omega }_0$. We vary the value of ${\omega }_0$ from -0.9 to -1.1. As ${\omega }_0$ decreases, ${\mathsf{\Omega }}_{m0}$ also decreases while b increases. However, there are no consistent trends observed in the changes of $M_0$ and h across these models. For models with ${\omega }_0\ge -1.0$, the best-fit values of b are negative. Conversely, we obtain positive best-fit values of b for models with ${\omega }_0\le -1.0$. Nonetheless, within the 1-$\sigma$ error range, the values of b for all models include both positive and negative values, suggesting no clear evidence of variations in the speed of light in these models.

3.2. $\tilde{c}$ for CPL

In this subsection, we repeat the maximum likelihood analysis for the CPL models. We explore various models across different values of cosmological parameters and summarize the results of this analysis shown in Table 2. Initially, we investigate GR (i.e., $b=0$) and subsequently extend the analysis to include the meVSL models (i.e., $b\ne 0$). Unlike the $\omega$CDM models, obtaining viable values for ${\mathsf{\Omega }}_{\mathrm{m}0}$ from this analysis without any prior constraints proves challenging. Consequently, we opt to perform the maximum likelihood analysis exclusively for fixed values of ${\mathsf{\Omega }}_{\mathrm{m}0}$.

3.2.1. $b=0$

We delve into CPL models within the framework of GR. Initially, we set $b=0$ and analyze three distinct scenarios: ${\omega }_0=-1$, ${\omega }_a=0$, and without fixing ${\omega }_0$ and ${\omega }_a$. For the ${\omega }_a=0$ model, we derive constraints of $0.312\le {\mathsf{\Omega }}_{\mathrm{m}0}\le 0.384$, $-1.37\le {\omega }_0\le -1.07$, and $0.610\le h\le 0.798$ at the 68

When ${\omega }_0=-1.0$, the range of ${\omega }_a$ spans $0.94\le {\omega }_a\le 1.78$, while ${\mathsf{\Omega }}_{\mathrm{m}0}$ ranges from $0.109$ to $0.217$ at the 68

Next, we explore meVSL models employing the CPL parameterization of dark energy, indicating $b\ne 0$ models.

3.2.2. ${\omega }_0=-0.95$ with Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

We conduct a maximum likelihood analysis for ${\omega }_0=-0.95$ models, varying ${\mathsf{\Omega }}_{\mathrm{m}0}$ from $0.28$ to $0.32$. Across these models, the best-fit values of h remain consistent at approximately $0.69$, with a 1-$\sigma$ error of $0.07$. As ${\mathsf{\Omega }}_{\mathrm{m}0}$ increases, so does the best-fit value of ${\omega }_a$, while the best-fit values of b decrease.

At a 68% confidence level, all b-values are negative for the given range of ${\mathsf{\Omega }}_{\mathrm{m}0}$, ranging from $-0.344$ to $-0.273$. This indicates that the speed of light is a monotonically decreasing function of cosmic time, with the rate of decrease accelerating as ${\mathsf{\Omega }}_{\mathrm{m}0}$ increases.

At $z=3$, $\tilde{c}$ can exceed ${\tilde{c}}_0$ by approximately ${10}^{13}$% when ${\mathsf{\Omega }}_{\mathrm{m}0}$ ranges from $0.28$ to $0.32$. The 1-$\sigma$ ranges of ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ (${10}^{-12}{\mathrm{yr}}^{-1}$) are ($-5.05,-4.55$) , (-5.38 , -4.88) , (-5.69 , -5.20) , (-5.99 , -5.50) , and (-6.31 , -5.82), respectively.

3.2.3. ${\omega }_0=-1.0$ with Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

We perform a comprehensive analysis for ${\omega }_0=-1.0$ models. Across these models, the best-fit values of h range from $0.689$ to $0.697$ for the given values of ${\mathsf{\Omega }}_{\mathrm{m}0}$. Meanwhile, the best-fit values of ${\omega }_a$ and b span $1.73\le {\omega }_a\le 1.97$ and $-0.293\le b\le -0.215$, respectively.

Similar to the ${\omega }_0=-0.95$ models, the best-fit values of b decrease as ${\mathsf{\Omega }}_{\mathrm{m}0}$ increases, with all b-values being negative at a 68% confidence level. Consequently, the speed of light exhibits a monotonically decreasing trend with cosmic time in these models. Notably, compared to the ${\omega }_0=-0.95$ model, the ${\omega }_0=-1.0$ model exhibits a slightly slower rate of decrease in the speed of light.

At $z=3$, $\tilde{c}$ exceeds ${\tilde{c}}_0$ by approximately $8\times 10$% when ${\mathsf{\Omega }}_{\mathrm{m}0}$ ranges from $0.28$ to $0.32$. Furthermore, the 1-$\sigma$ ranges of ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ (${10}^{-12}{\mathrm{yr}}^{-1}$) are ($-4.05,-3.52$) , (-4.45 , -3.92) , (-4.77 , -4.24) , (-5.09 , -4.56) , and (-5.42 , -4.89), respectively.

3.2.4. ${\omega }_0=-1.05$ with Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

In this subsection, we undertake a maximum likelihood analysis for ${\omega }_0=-1.05$ models. Across these models, the best-fit values of h remain consistently near $0.69$ for all specified values of ${\mathsf{\Omega }}_{\mathrm{m}0}$. The range of best-fit values for ${\omega }_a$ spans from $1.57$ to $1.84$, while for b, it extends from $-0.240$ to $-0.154$. As with previous models, the values of b decrease with increasing ${\mathsf{\Omega }}_{\mathrm{m}0}$, ranging from $-0.169$ to $-0.139$ for ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28$, and from $-0.255$ to $-0.225$ for ${\mathsf{\Omega }}_{\mathrm{m}0}=0.32$ within 1-$\sigma$ error. Hence, the speed of light maintains its monotonic decrease over cosmic time in these models. At $z=3$, $\tilde{c}$ can exceed ${\tilde{c}}_0$ by approximately 6 (9) % for ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28$ ($0.32$). Compared to the ${\omega }_0=-0.95$ and $-1.0$ models, the rate of decrease in the speed of light is slightly smaller in this model. The 1-$\sigma$ ranges of ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ (${10}^{-12}{\mathrm{yr}}^{-1}$) are ($-2.98,-2.45$) , (-3.44 , -2.82) , (-3.70 , -3.17) , (-4.01 , -3.49) , and (-4.50 , -3.97), respectively.

3.2.5. CPL with Fixing ${\mathsf{\Omega }}_{\mathrm{m}0}$

We analyze the Pantheon data without constraining values of ${\omega }_0$ and ${\omega }_a$ within the range $0.28\le {\mathsf{\Omega }}_{\mathrm{m}0}\le 0.32$. Across these models, the best-fit values of h remain approximately constant at $0.71$ for all specified values of ${\mathsf{\Omega }}_{\mathrm{m}0}$. The best-fit values of both ${\omega }_0$ and ${\omega }_a$ range as ($-1.18,0.72$), ($-1.18,0.76$), and ($-1.20,0.71$) for ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28,0.30$, and $0.32$, respectively. Of particular interest is the model with ${\mathsf{\Omega }}_{\mathrm{m}0}=0.30$, where the best-fit value of b is nearly zero, varying within the range $-0.015$ to $0.015$ at a 68 % confidence level. This suggests no significant time variation in the speed of light in this model. However, for ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28$ and $0.32$, the 1-$\sigma$ values of b indicate a clear monotonic decrease (increase) in the speed of light with cosmic time, with ranges of $0.023\le b\le 0.053$ and $-0.039\le b\le -0.009$, respectively.

These trends are illustrated in Fig. Figure 2. In the left panel, the cosmological evolution of $\tilde{c}/{\tilde{c}}_0$ as a function of z for the ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28$ model is depicted. The dashed line represents the best-fit value of b, while the solid lines indicate the 1-$\sigma$ errors. The monotonically decreasing behavior of $\tilde{c}/{\tilde{c}}_0$ with increasing redshift is evident due to the positive values of b. At $z=3$, $\tilde{c}$ decreases by approximately $0.8$ ($1.8$) % within 1-$\sigma$ error. The 1-$\sigma$ range of ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ (${10}^{-13}{\mathrm{yr}}^{-1}$) falls between $4.15$ and $9.57$ at a 68 % confidence level.

On the other hand, the right panel of Fig. Figure 2 illustrates the model with ${\mathsf{\Omega }}_{\mathrm{m}0}=0.32$, where both the best-fit and 1-$\sigma$ error values of b are negative. Consequently, $\tilde{c}/{\tilde{c}}_0$ increases monotonically with redshift. At $z=3$, $\tilde{c}$ increases by about $0.3$ ($1.4$) % at a 68 % confidence level. The 1-$\sigma$ range of ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ (${10}^{-13}{\mathrm{yr}}^{-1}$) spans from $(-7.05,-1.63)$.

3.2.6. CPL with (without) Fixing h

The analysis is conducted with fixed values of h while leaving other parameters unconstrained. For $h=0.6736(0.74)$, the best-fit values and 1-$\sigma$ errors for ${\omega }_0$ and ${\omega }_a$ are $-0.99\pm 0.02(-0.99\pm 0.01)$ and $1.31\pm 0.11(1.45\pm 0.07)$, respectively. Notably, in these models, the values of b are positive. However, the resulting values for ${\mathsf{\Omega }}_{\mathrm{m}0}$ fall outside of the viable range, yielding $-0.028\le {\mathsf{\Omega }}_{\mathrm{m}0}\le 0.094(-0.148\le {\mathsf{\Omega }}_{\mathrm{m}0}\le -0.076)$ for $h=0.6736(0.74)$, which are deemed nonviable. Moreover, without constraining ${\mathsf{\Omega }}_{\mathrm{m}0}$, the matter density contrast is estimated to be $-0.14\le {\mathsf{\Omega }}_{\mathrm{m}0}\le 0.032$, again falling outside of acceptable bounds. These results underscore the importance of appropriately constraining cosmological parameters to ensure the viability of the models.

3.3. $\dot{c}$ and $\dot{G}$

Building upon the findings from previous subSection 3.1 and Section 3.2, we delve into the exploration of viable meVSL models across different dark energy scenarios, extracting valuable constraints on both cosmological and model parameters. Through this analysis, we can derive estimates for the temporal evolution of the speed of light, a fundamental aspect within the meVSL framework. Interestingly, in the meVSL paradigm, not only does the speed of light undergo cosmological evolution, but also the gravitational constant, characterized by $\tilde{G}={\tilde{G}}_0{(1+z)}^{-b}$. Leveraging the constraints obtained on b-values, we can ascertain bounds on the present value of the relative temporal variation of the gravitational constant, ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$.

To contextualize these findings within the broader scope of observational constraints, we juxtapose our results with existing bounds on ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ from diverse sources, as delineated in Table 3. Notably, the analysis of lunar laser ranging (LLR) data emerges as particularly stringent, yielding the lowest bounds on ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$. Conversely, the orbital period rate of pulsars offers the largest bounds, estimated at $2.3\times {10}^{-11}{\mathrm{yr}}^{-1}$. Collectively, these observations paint a picture where ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ is tentatively within the order of ${10}^{-12}{\mathrm{yr}}^{-1}$. Such insights not only enrich our understanding of cosmological dynamics but also pave the way for future investigations into the fundamental nature of physical constants.

The table displays the outcomes regarding the temporal variations of both the speed of light and the gravitational constant within meVSL models, considering various dark energy scenarios. We define $\Delta \tilde{c}(z=3)$ as the percentage deviation between the speed of light’s value at redshift $z=3$ and its current value, expressed as $\Delta \tilde{c}(z=3)\equiv \left(\tilde{c}(z=3)-{\tilde{c}}_0\right)/{\tilde{c}}_0\times 100(\%)$ within a 1-$\sigma$ uncertainty range. Similarly, $\Delta \tilde{G}(z=3)$ indicates the percentage deviation between the gravitational constant’s values at $z=3$ and $z=0$. The ratio of the temporal variation of the speed of light to its current value is denoted by ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$, while ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ represents the current ratio of the gravitational constant’s temporal variation to its value.

In the analysis, positive values of the best-fit parameter b and its 68% confidence level values are only observed for the CPL dark energy model when ${\mathsf{\Omega }}_{\mathrm{m}0}=0.28$. Consequently, both ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ and ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ exhibit positivity in this specific model. Conversely, all other viable models derived from the Pantheon data yield negative b values, resulting in negative values for both ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ and ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ in these scenarios.

The ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ values are on the order of ${10}^{-13}{\mathrm{yr}}^{-1}$, while ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ values are on the order of ${10}^{-12}{\mathrm{yr}}^{-1}$ for $\omega$CDM models and CPL models with varying ${\omega }_0$ and ${\omega }_a$, as delineated in Table 4. However, for CPL models with fixed ${\omega }_0$, the ${\dot{\tilde{c}}}_0/{\tilde{c}}_0$ values are approximately on the order of ${10}^{-12}{\mathrm{yr}}^{-1}$, while ${\dot{\tilde{G}}}_0/{\tilde{G}}_0$ values tend to be on the order of ${10}^{-11}{\mathrm{yr}}^{-1}$.

4. Discussion

The Pantheon+ data provides constraints on cosmological and model parameters with a statistical precision of about 10 %. Leveraging this data, we perform a maximum likelihood analysis to constrain dark energy models within the framework of the meVSL model. Through this analysis, we identify several viable $\omega$CDM and CPL dark energy models and derive constraints on the parameter b, which governs the evolution of physical constants in the universe. The obtained constraints indicate that the relative temporal variations of the speed of light and the gravitational constant lie within the ranges $-64.0\le {\dot{\tilde{c}}}_0/{\tilde{c}}_0\left({10}^{-13}{\mathrm{yr}}^{-1}\right)\le -0.39$ and $-25.6\le {\dot{\tilde{G}}}_0/{\tilde{G}}_0\left({10}^{-12}{\mathrm{yr}}^{-1}\right)\le -0.36$ for most viable models. Consequently, our analysis suggests that according to current Pantheon data, both the speed of light and the gravitational constant were greater in the past and have monotonically increased with redshift, z.

While additional cosmological observations such as CMB and BAO could potentially further constrain cosmological and model parameters, integrating these datasets would necessitate a reanalysis based on the theoretical framework of the meVSL model. However, this endeavor falls outside the scope of the present manuscript and is deferred to future investigations.