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

Abstract

In the context of the minimally extended varying speed of light (meVSL) model, both the absolute magnitude and the luminosity distance of type Ia supernovae (SNe Ia) deviate from those predicted by general relativity (GR). Using data from the Pantheon+ survey, we assess the plausibility of various dark energy models within the framework of meVSL. Both the constant equation of state (EoS) of the dark energy model ($\omega$CDM) and the Chevallier–Polarski–Linder (CPL) parameterization model ($\omega ={\omega }_0+{\omega }_a(1-a)$) indicate potential variations in the cosmic speed of light at the $1-\sigma$ confidence level. For ${\mathrm{\Omega }}_{\mathrm{m}0}=0.30,0.31$, and 0.32 with $({\omega }_0,{\omega }_a)=(-1,0)$, the $1-\sigma$ range of ${\dot{c}}_0/c_0\left({10}^{-13}{\mathrm{yr}}^{-1}\right)$ is (−8.76, −0.89), (−11.8, 3.93), and (−14.8, −6.98), respectively. Meanwhile, the $1-\sigma$ range of ${\dot{c}}_0/c_0\left({10}^{-12}{\mathrm{yr}}^{-1}\right)$ for CPL dark energy models with $-1.05\le {\omega }_0\le -0.95$ and $0.28\le {\mathrm{\Omega }}_{\mathrm{m}0}\le 0.32$ is (−6.31, −2.98). The value of c at $z=3$ can exceed that of the present by 0.2∼3% for $\omega$CDM models and 5∼13% for CPL models. Additionally, for viable models except for the CPL model with ${\mathrm{\Omega }}_{\mathrm{m}0}=0.28$, we find $-25.6\le {\dot{G}}_0/G_0\left({10}^{-12}{\mathrm{yr}}^{-1}\right)\le -0.36$. For this particular model, we obtain an increasing rate of the gravitational constant within the range $1.65\le {\dot{G}}_0/G_0\left({10}^{-12}{\mathrm{yr}}^{-1}\right)\le 3.79$. We obtain some models that do not require dark matter energy density through statistical interpretation. However, this is merely an effect of the degeneracy between model parameters and energy density and does not imply that dark matter is unnecessary.

1. Introduction

The principles of physics should remain unchanged regardless of the chosen units or measurement instruments. This is exemplified by dimensionless quantities such as the fine structure constant $\alpha$, which retains its value across different unit systems as demonstrated in the Standard Model of particle physics. Conversely, dimensional constants like ℏ, c, G, e, and k are human conventions whose numerical values vary depending on the units employed. Consequently, only dimensionless constants can be considered truly fundamental. Therefore, investigating potential variations over time in dimensionless fundamental constants is a valid scientific pursuit, while variations in dimensional constants like c or G depend on the chosen units and may result in discrepancies among observers. This perspective is supported by various studies [1,2,3,4], albeit within the framework of a static universe or one existing at the present epoch [5,6,7].

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. The redshift, defined as a function of the time-evolving cosmic scale factor $a\left(t\right)$, yields positive values for z in our expanding universe. Estimating the redshift of a galaxy involves analyzing the emission lines emitted by glowing gas within the galaxy. For example, the H$\alpha$ line, a red Balmer line of neutral hydrogen, has a rest wavelength of 6562 Å. If the observed wavelength of this line presently measures 8100Å, it indicates that the galaxy is positioned at $z=0.234$ (i.e., $a=0.81$). Consequently, in an expanding universe, the value of dimensional quantities, such as wavelength, varies depending on the observation time (i.e., cosmic time). Additionally, it has been observed that the temperature of the CMB decreases with the age of the universe, scaling inversely with the scale factor as $T=T_0{a}^{-1}$ [8,9,10,11]. 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$ varies 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 a function of the scale factor, allowing the speed of light to be expressed as $c\left(t_1\right)={(a\left(t_1\right)/a\left(t_2\right))}^{b/4}c\left(t_2\right)$ in the meVSL model, where b is a constant [5,6,7]. Also, to preserve the Einstein Field Equations (EFEs), the Einstein constant, denoted as $\kappa =8\pi G/c^4$, must remain constant even if both c and G undergo cosmological evolution. In the meVSL model, where $c=c_0{a}^{b/4}$, this implies $G=G_0{a}^b$, where $c_0$ and $G_0$ represent the current values of the speed of light and the gravitational constant, respectively [5,6,7].

In addition to models like meVSL, there are models, such as co-varying physical couplings (CPCs) [12,13,14], where physical constants vary with cosmic time. The CPC retains the EFEs, with G, c, and $\mathrm{\Lambda }$ treated as functions of spacetime. The interaction between the Bianchi identity and the requirement of stress–energy tensor conservation complicates the potential variations of these constants, which are constrained to co-vary according to the General Constraint (GC). Unlike meVSL, this model includes the dynamics of physical constants through the adoption of GC.

Special relativity (SR)’s universal Lorentz covariance, grounded in Minkowski spacetime, adequately upholds its principles [15]. In contrast, within GR, an inertial frame (IF) refers to one that is freely falling. While Lorentz invariant (LI) spacetime intervals can be established between events, the definition of a global time in GR is impeded by the absence of a universal IF. However, a global time can be delineated for the universe, satisfying the cosmological principle (CP), enabling a foliation of spacetime into non-intersecting spacelike 3D surfaces. This description pertains to the universe modeled by the RW metric [8,9,10,11]. The LI varying speed of light (VSL) model is plausible if the speed of light, denoted as c, remains locally constant at each given epoch but varies in cosmic time [5,6,7]. In other words, in an expanding universe, if the speed of light is expressed as a function of the scale factor, $c\left[a\right]$, then although its value changes akin to wavelengths at different epochs, such as $a_1$ and $a_2$, it maintains a constant local value at each epoch, ensuring LI and thereby preserving the validity of quantum mechanics and electromagnetism in accordance with SR during every epoch. 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 [5,6,7,16]. 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 [5,6,7]. 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 the 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 [5,6,7].

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 [17,18,19,20,21]. Another method entails measuring TD by examining the stretching of peak-to-peak timescales of gamma-ray bursters (GRBs) [22,23,24,25,26,27,28,29,30]. Efforts have also been made to detect the TD effect in the light curves of quasars (QSOs) located at cosmological distances [31,32]. 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 [5,7].

First, there is the cosmic distance duality relation CDDR method among these methods. 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 [33]. 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 [34,35]. Various tests of the CDDR using astrophysical and cosmological observations have been conducted to constrain VSL models [13,36,37,38,39,40,41,42,43,44,45]. 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 [5,44]. 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 [46]. 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$) [5,6,7]

$$\begin{array}{ll}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}\\ & =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 [47]. 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 [48,49,50,51,52,53,54,55,56,57,58,59,60]. This method proves particularly valuable for investigating VSL models [45]. 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 [61].

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 [62]. This information is crucial for studying the universe’s expansion rate and testing cosmological models, involving data collection, mathematical modeling, and the comparison of observed TD with theoretical predictions. We derive 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$ [21,63]. 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 [64].

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 [65]. 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 $\mathrm{\Lambda }$CDM and $\omega$CDM models, without accounting for systematic uncertainties on ${\mathrm{\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 ${\mathrm{\Omega }}_{m0}$ and 12% on $\omega$ for the flat $\omega$CDM model, and about 4% on ${\mathrm{\Omega }}_{m0}$ for the flat $\mathrm{\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. In Section 3, we elucidate the statistical methods employed in this paper. Additionally, a portion of the Mathematica file used for this purpose is included in Appendix A. In Section 4, we provide a brief explanation of the temporal variations of the speed of light and Newton’s constant used in this paper, along with a concise interpretation for the $\mathrm{\Lambda }$CDM model for comparison. The investigation into models capable of implementing variations in the speed of light across various $\omega$CDM models is conducted in Section 5. Research on potential models within the CPL framework of the meVSL model is performed in Section 6. In Section 7, we offer various possible constraints on the temporal variations of the speed of light and Newton’s constant within the meVSL model, comparing them with existing observational constraints. Finally, in Section 8, 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 cosmic time [8,9,10,11]. At the heart of this framework lies the RW metric, elegantly expressed by

$$\begin{array}{ll}d{s}^2& =-c{\left(t\right)}^{2}d{t}^2+a^2\left(t\right){\gamma }_{ij}d{x}^{i}d{x}^j=-c{\left(t\right)}^{2}d{t}^2+a^2\left(t\right)\left(\frac{d{r}^2}{1-k{r}^2}+r^{2}d{\mathrm{\Omega }}^2\right)\\ & =-c{\left(t\right)}^{2}d{t}^2+a^2\left(t\right)\left[d{\chi }^2+f_k^2\left(\chi \right)d{\mathrm{\Omega }}^2\right]\equiv -c{\left(t\right)}^{2}d{t}^2+a{\left(t\right)}^{2}d{l}_{3\mathrm{D}}^2,\end{array}$$

(2)

where $c\left(t\right)$ 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 the physical distance and comoving distance. We adopt that the speed of light is a function of cosmic time, deviating from the conventional RW metric. We already show that there is no contradiction in this assumption as long as we adopt the CP and Weyl’s postulate [6,7]. The derivation of the 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}}& \hfill =\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=1/\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 b characterizes the deviation of c from the constant value. Thus, the RW metric naturally allows the VSL models if we remove the traditional assumption on the cosmological TD [5,6,7].

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}^{\prime 2}}}=\frac{c_0}{H_0}{\int }_0^z\frac{d{z}^{\prime }}{E^{\left(\mathrm{SMC}\right)}\left(z^{\prime }\right)}\equiv \frac{c_0}{H_0}{d}_{\mathrm{C}}\left(z\right),\end{array}$$

(4)

where $d_{\mathrm{L}}$ is the so-called Hubble free luminosity distance 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{SMC}\right)}\left(z\right)=\left\{\begin{array}{cc}\hfill \frac{c_0}{H_0}\frac{1}{\sqrt{{\mathrm{\Omega }}_{k0}}}sinh\left(\sqrt{{\mathrm{\Omega }}_{k0}}\frac{H_0}{{\tilde{c}}_0}{D}_{\mathrm{C}}\right)\hfill & \hfill {\mathrm{\Omega }}_{k0}>0\hfill \\ \hfill D_{\mathrm{C}}\hfill & \hfill {\mathrm{\Omega }}_{k0}=0\hfill \\ \hfill \frac{c_0}{H_0}\frac{1}{\sqrt{|{\mathrm{\Omega }}_{k0}|}}sin\left(\sqrt{|{\mathrm{\Omega }}_{k0}|}\frac{H_0}{{\tilde{c}}_0}{D}_{\mathrm{C}}\right)\hfill & \hfill {\mathrm{\Omega }}_{k0}<0\hfill \end{array}\right.,\end{array}$$

(5)

with $c_0=2.9979\times {10}^5$ km/s and $H_0=100$ h km/Mpc/s representing the present values of the speed of light and the Hubble parameter, respectively [64,66]. In a recent local measurement by the SH0ES collaboration, utilizing Cepheid-calibrated SNeIa, the Hubble constant was reported to be approximately 73 km/s/Mpc [67], contrasting with the prediction of around 67 km/s/Mpc by the standard $\mathrm{\Lambda }$CDM model based on observations of the CMB [68].

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 at a given epoch, with the slicing being orthogonal to these hypersurfaces, 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 [8,9,10,11]. 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 Equation (2) further enriches the understanding of the underlying spacetime dynamics and its mathematical representation [5]. Specifically, the expressions for $R_{00}$ and $R_{ii}$ are given by

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

(6)

$$\begin{array}{cc}\hfill & \hfill R=\frac{6}{c^2}\left(\frac{\ddot{a}}{a}+\frac{{\dot{a}}^2}{a^2}+k\frac{c^2}{a^2}-H^2\frac{dlnc}{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 is its rest frame mass density for a perfect fluid. 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 (i.e., comoving) observers is considered comoving with it, characterized by a four-velocity denoted as $U^{\mu }=(c,0,0,0)$ [8,10]. 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 [5]

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

(9)

where $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 general DE from Equations (6)–(9)

$$\begin{array}{l}\frac{{\dot{a}}^2}{a^2}+k\frac{c^2}{a^2}=\sum _i\frac{8\pi G}{3}{\rho }_i=\frac{8\pi G}{3}\left[{\rho }_{\mathrm{r}}\left(a\right)+{\rho }_{\mathrm{m}}\left(a\right)+{\rho }_{\mathrm{DE}}\left(a\right)\right],\\ \frac{{\dot{a}}^2}{a^2}+2\frac{\ddot{a}}{a}+k\frac{c^2}{a^2}-2H^2\frac{dlnc}{dlna}=-8\pi G\sum _i\frac{P_i}{c^2}=-8\pi G\sum _i{\omega }_i{\rho }_i\end{array}$$

(10)

$$\begin{array}{l}=-8\pi G\left[\frac{1}{3}{\rho }_{\mathrm{r}}\left(a\right)+{\omega }_{\mathrm{DE}}{\rho }_{\mathrm{DE}}\left(a\right)\right],\\ \frac{\ddot{a}}{a}=-\frac{4\pi G}{3}\sum _i\left(1+3{\omega }_i\right){\rho }_i+H^2\frac{dlnc}{dlna}\end{array}$$

(11)

$$=-\frac{4\pi G}{3}\left(2{\rho }_{\mathrm{r}}\left(a\right)+{\rho }_{\mathrm{m}}\left(a\right)+\left(1+{\omega }_{\mathrm{DE}}\right){\rho }_{\mathrm{DE}}\left(a\right)\right)+\frac{b}{4}{H}^2,$$

(12)

$${\rho }_{\mathrm{r}}\left(a\right)={\rho }_{\mathrm{r}0}{a}^{-4-\frac{b}{2}},{\rho }_{\mathrm{m}}\left(a\right)={\rho }_{\mathrm{m}0}{a}^{-3-\frac{b}{2}},{\rho }_{\mathrm{DE}}\left(a\right)={\rho }_{\mathrm{DE}0}{a}^{-3(1+{\omega }_0+{\omega }_a)-\frac{b}{2}}{e}^{-3{\omega }_a(1-a)},$$

(13)

$${\omega }_{\mathrm{DE}}={\omega }_0+{\omega }_a\left(1-a\right)=\left\{\begin{array}{cc}{\omega }_0,{\omega }_a\ne 0\hfill & \mathrm{CPL}\hfill \\ {\omega }_a=0\hfill & \omega \hfill \end{array}\right.,$$

(14)

where ${\rho }_{\mathrm{r}0}$, ${\rho }_{\mathrm{m}0}$, and ${\rho }_{\mathrm{D}\mathrm{E}0}$ represent the mass density of radiation (photon and neutrino), matter (baryon and DM), and dark energy, respectively, at the present epoch [5,6]. The Chevallier–Polarski–Linder (CPL) parametrization, which assumes ${\omega }_{\mathrm{DE}}$ to be a linear function of the scale factor a, is presented in Equation (14) [69,70]. We also define the $\omega$ model when ${\omega }_a=0$.

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}$$

(15)

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 $\mathcal{F}$ is given by

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

(16)

The expression for E^(SMC) in Equation (1) for a flat universe (using the CPL parametrization) is obtained from Equation (10)

$$\begin{array}{c}\hfill \frac{H^{\left(\mathrm{SMC}\right)}}{H_0}\equiv E^{\left(\mathrm{SMC}\right)}\simeq \sqrt{{\mathrm{\Omega }}_{\mathrm{m}0}{a}^{-3}+\left(1-{\mathrm{\Omega }}_{\mathrm{m}0}\right)a^{-3(1+{\omega }_0+{\omega }_a)}{e}^{-3{\omega }_a(1-a)}},\end{array}$$

(17)

where we ignore the radiation and the curvature contribution because we analyze the late-time universe. 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}$$

(18)

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 Equation (18) as [5,6]

$$\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}$$

(19)

where we employ

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

(20)

where we use the cosmic evolution relation for the Planck constant $h_i=h_0{a}_i^{-b/4}$. In the meVSL model, an expanding universe must adhere to adiabatic expansion, leading to the cosmological evolution of the Planck constant [16]. The first law of thermodynamics, which ensures energy conservation, requires that the entropy of the universe remains unchanged. Consequently, the Planck constant should evolve as $h\left(a_i\right)=h_0{a}_i^{-b/4}$ in this model [5,6,16]. This relation also holds for the angular diameter distance $D_{\mathrm{A}}$. Consequently, the luminosity distance in the meVSL model is given by [5]

$$\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}$$

(21)

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

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

(22)

where the subscript 0 denotes the local value of M. The last equality of Equation (22) can be obtained from the following [5]. SNe Ia are nuclear explosions of white dwarfs (WDs) in binary systems, where the WD accretes matter from a companion until it approaches the Chandrasekhar limit. This limit represents the maximum mass a WD can have before electron degeneracy pressure fails to counteract the gravitational collapse. For WDs, this limit is typically around $1.4$ solar masses. If a WD exceeds this mass, it can collapse into a neutron star or black hole, while those below remain stable.

The Chandrasekhar mass limit M_Ch is determined by the equation of state for an ideal Fermi gas, with a constant ${\omega }_0^3\approx 2.018$ (a constant related to the solution for the so-called Land–Emden equation), the average molecular weight per electron ${\mu }_e$, and the mass of the hydrogen atom $m_{\mathrm{H}}$ [5]

$$\begin{array}{c}\hfill {\mathrm{M}}_{\mathrm{Ch}}=\frac{{\omega }_0^3\sqrt{3\pi }}{2}{\left(\frac{\hslash c}{G}\right)}^{\frac{3}{2}}\frac{1}{{\left({\mu }_e{m}_{\mathrm{H}}\right)}^2}\equiv {\mathrm{M}}_{\mathrm{Ch}0}{a}^{-\frac{b}{2}}.\end{array}$$

(23)

The peak luminosity of SNe Ia is proportional to the mass of synthesized nickel, which is a fraction of the Chandrasekhar mass. Consequently, the absolute magnitude of SNe Ia, which measures luminosity, is related to the Chandrasekhar mass and the total amount of nickel synthesized $L\propto {\mathrm{M}}_{\mathrm{Ch}}\propto a^{-\frac{b}{2}}$. The absolute magnitude M is given by $M\propto -2.5log\left[L\right]\propto \frac{5}{4}blog\left[a\right]$. Thus, the distance modulus of meVSL $\mu \equiv m-M$ is written as:

$$\begin{array}{cc}\hfill \mu \left(z\right)& \hfill =m-M=5{log}_{10}\left[\frac{D_{\mathrm{L}}}{\mathrm{Mpc}}\right]+25,D_{\mathrm{L}}=\frac{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{SMC}\right)}\left(z^{\prime }\right)},\hfill \end{array}$$

(24)

where $E^{\left(\mathrm{SMC}\right)}$ is in Equation (17). The theoretically predicted apparent magnitude $m_{\mathrm{th}}$ can be obtained from Equations (22) and (24)

$$\begin{array}{ll}{m}_{\mathrm{th}}\left(z\right)& =M_0+5log\left[\frac{D_{\mathrm{L}}}{1\mathrm{Mpc}}\right]+25+\frac{5}{4}blog\left[a\right]\\ & =M_0+42.3841-5log\left[\mathrm{h}\right]-\frac{15}{8}blog\left[(1+z)\right]+5log\left[d_{\mathrm{L}}\right]\\ & \equiv \mathcal{M}+5log\left[d_{\mathrm{L}}\right].\end{array}$$

(25)

3. Statistical Analysis

We examine constraints on the evidence for cosmic variation in the speed of light using the Pantheon+ dataset [65]. The Pantheon compilation includes 1048 SNe Ia, covering a redshift range $0.01\le z\le 2.3$. This dataset incorporates 365 spectroscopically confirmed SNe Ia from the Pan-STARRS1 (PS1) Medium Deep Survey, along with a subset of 279 PS1 SNe Ia (with redshifts ranging from $0.03$ to $0.68$) with reliable distance estimates derived from various sources, including SDSS, SNLS, and HST samples. To determine cosmological parameters using $H\left(z\right)$, higher-redshift SNe Ia are employed, and the degenerate parameters $\mathcal{M}$ are typically marginalized as nuisance parameters [65,71]. For example, Equation (25) is used to construct and minimize ${\overline{\chi }}^2\left({\mathrm{\Omega }}_{\mathrm{m}0}\right)\equiv \int d\mathcal{M}{\chi }^2(\mathcal{M},{\mathrm{\Omega }}_{\mathrm{m}0})$, where the degenerate combination is provided in the same equation [72,73]. However, marginalizing the parameter $\mathcal{M}$ can result in the loss of valuable physical information regarding potential spatial variations of $H_0$ and/or temporal variations of the absolute magnitude $\mathcal{M}$. For instance, an $\mathcal{M}$ value that evolves with redshift, resulting in low $\mathcal{M}$ values at low z, could indicate either higher local values of $H_0$ due to a local matter underdensity or lower values of the absolute magnitude $\mathcal{M}$ in recent cosmological times, possibly caused by a time variation of Newton’s constant. Specifically, since we are studying the potential variation of the speed of light with cosmic time, we do not marginalize other parameters in this manuscript when we use the minimal ${\chi }^2$ method [74,75,76,77]. 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}$$

(26)

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 Equation (24), 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. $m_{\mathrm{th}}$ is a function of $M, \mathrm{h}, {\mathrm{\Omega }}_{\mathrm{m}0}, {\omega }_0,$ and ${\omega }_a$. 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}$$

(27)

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. Since all models have nearly identical reduced chi-square values, this implies that all models are statistically equally viable. We have included a portion of the Mathematica code used in this paper for reference in the Appendix A.

4. Bounds for the Variation of $\mathit{c}$ for Different Models

The previous constraint on the temporal variation of the speed of light c was derived from the variation in the radius of a planet [78]. However, this constraint stemmed from the analysis of the time-varying radius of Mercury [79] using a specific model known as the covariant variable speed of light theory proposed by Magueijo [80]. 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 [65]. 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.

First, in the first row of

Table 1. Best-fit values and their corresponding $1-\sigma$ uncertainties for both $\mathrm{\Lambda }$CDM (${\omega }_a=0$ and $b=0$) and $\omega$CDM models (with ${\omega }_a=0$) are presented. Models highlighted in green indicate potential meVSL models for $\omega$CDM.

Models	Submodels	${\mathit{M}}_0$	${\mathit{\omega }}_0$	${\mathbf{\Omega }}_{\mathit{m}0}$	h	b	$\mathit{\nu }$	${\mathit{\chi }}_{\mathit{\nu }}^2$
$\mathrm{\Lambda }$CDM	fixing $M_0$	$-19.3500$	$-1$	$0.285\pm 0.012$	$0.702\pm 0.002$	0	1046	$0.988$
		$-19.4500$	$-1$	$0.285\pm 0.012$	$0.670\pm 0.002$	0	1046	$0.989$
		$-19.5500$	$-1$	$0.285\pm 0.012$	$0.640\pm 0.002$	0	1046	$0.989$
		$-19.3622\pm 0.2890$	$-1$	$0.285\pm 0.012$	$0.698\pm 0.093$	0	1045	$0.989$
${\omega }_0=-1$	fixing ${\mathrm{\Omega }}_{\mathrm{m}0}$	$-19.3556\pm 0.2899$	$-1$	$0.28$	$0.700\pm 0.093$	$0.009\pm 0.022$	1045	$0.989$
		$-19.3558\pm 0.2899$	$-1$	$0.29$	$0.700\pm 0.093$	$-0.009\pm 0.022$	1045	$0.989$
		$-19.3561\pm 0.2899$	$-1$	$0.30$	$0.700\pm 0.093$	$-0.027\pm 0.022$	1045	$0.989$
		$-19.3563\pm 0.2899$	$-1$	$0.31$	$0.700\pm 0.093$	$-0.044\pm 0.022$	1045	$0.989$
	Section 5.1	$-19.3566\pm 0.2898$	$-1$	$0.32$	$0.700\pm 0.093$	$-0.061\pm 0.022$	1045	$0.989$
	fixing h	$-19.4395\pm 0.0072$	$-1$	$0.299\pm 0.111$	$0.6736$	$-0.025\pm 0.193$	1045	$0.989$
	Section 5.3	$-19.2353\pm 0.0072$	$-1$	$0.299\pm 0.111$	$0.74$	$-0.025\pm 0.193$	1045	$0.989$
$\omega$CDM	fixing h	$-19.4525\pm 0.0071$	$-1.23\pm 0.05$	$0.380\pm 0.085$	$0.6736$	$-0.057\pm 0.159$	1044	$0.988$
	Section 5.4	$-19.2481\pm 0.0071$	$-1.22\pm 0.05$	$0.378\pm 0.087$	$0.74$	$-0.055\pm 0.163$	1044	$0.988$
	fixing ${\omega }_0$	$-19.3655\pm 0.2893$	$-0.9$	$0.311\pm 0.106$	$0.695\pm 0.093$	$-0.106\pm 0.163$	1044	$0.993$
		$-19.3659\pm 0.2847$	$-0.95$	$0.301\pm 0.092$	$0.696\pm 0.091$	$-0.058\pm 0.152$	1044	$0.992$
		$-19.3728\pm 0.2877$	$-1.0$	$0.299\pm 0.111$	$0.695\pm 0.092$	$-0.025\pm 0.193$	1044	$0.990$
		$-19.3676\pm 0.2908$	$-1.05$	$0.290\pm 0.110$	$0.697\pm 0.093$	$0.023\pm 0.305$	1044	$0.989$
	Section 5.5	$-19.3642\pm 0.2895$	$-1.1$	$0.288\pm 0.134$	$0.699\pm 0.093$	$0.054\pm 0.262$	1044	$0.988$
	No fixing	$-19.5513\pm 0.1284$	$-1.22\pm 0.04$	$0.335\pm 0.088$	$0.644\pm 0.038$	$0.022\pm 0.176$	1043	$0.989$

, we considered the $\mathrm{\Lambda }$CDM models (${\omega }_0=-1,{\omega }_a=0$, and $b=0$). When varying the local absolute magnitude $M_0$ from −19.35 to −19.55, we found that the best-fit value of ${\mathrm{\Omega }}_{\mathrm{m}0}$ remains unchanged at 0.285, while the best-fit values of h decrease from $0.702$ to $0.640$. For this model, the $1-\sigma$ values for $M_0$, ${\mathrm{\Omega }}_{\mathrm{m}0}$, and h are $(-19.6512,-19.0732)$, $(0.273,0.297)$, and $(0.605,0.791)$, respectively. This suggests the interesting possibility that the h value derived from SNe Ia data can be similar to that of Planck data, potentially alleviating the Hubble tension. If we also allow the value of $M_0$ to vary, its best-fit value becomes $-19.36$, and the best-fit value of h is $0.698$. Similar results were discussed in the reference [81]. $\mathrm{\Lambda }$CDM refers to the case where we exclude VSL models. We present this model solely for comparison purposes. In our manuscript, we investigate the meVSL model for $\omega$CDM and CPL models.

The main interest of this manuscript is whether the speed of light can vary in models with values of cosmological parameters similar to those of the SMC. Therefore, we will limit our discussions to such models. Among the CPL models, there are cases where the dark matter density is zero. There are also interesting papers that have interpreted these results [82,83]. The so-called CCC (covarying coupling constants) + TL (tired light) cosmology obtains ${\mathrm{\Omega }}_{\mathrm{m}0}=0$ using the baryonic acoustic oscillation data [82]. One can find the origin of dark matter and dark energy as a time-dependent conformal scale factor in the Scale Invariant Vacuum (SIV) paradigm [83]. However, we will treat this outcome as merely a result of statistical analysis and will not consider any further interpretation or implication. Although this could be an interesting topic for another paper, as mentioned earlier, in this manuscript, we will focus primarily on models where the speed of light can vary with cosmological parameters consistent with the SMC. Also, as is evident from Equations (17) and (25), we acknowledge the degeneracy relationship between ${\mathrm{\Omega }}_{\mathrm{m}0}$ and the ${\omega }_a$ term in the exponential exponent. Thus, we accept that the value of ${\mathrm{\Omega }}_{\mathrm{m}0}$ can vary due to this term, and we refrain from further physical interpretations or alternative explanations.

5. $\mathit{c}$ for $\mathit{\omega }$CDM

We explore the $\omega$CDM models ${\omega }_a=0$ in Equation (14). 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 significant patterns and correlations among cosmological parameters. Notably, when ${\omega }_0$ is held constant, both $M_0$ and $\mathrm{h}$ exhibit degeneracy as do ${\mathrm{\Omega }}_{\mathrm{m}0}$ and b. Consequently, fixing ${\mathrm{\Omega }}_{\mathrm{m}0}$ results in nearly identical values for $M_0$ and h, with only b values varying. Conversely, when we fix the value of $\mathrm{h}$, only the $M_0$ values change, while ${\mathrm{\Omega }}_{\mathrm{m}0}$ and b become irrelevant in this context.

Furthermore, we observe that among the $\omega$CDM models derived from the Pantheon+ data, those displaying noticeable time variations in the speed of light are specifically characterized by ${\omega }_0=-1$ and ${\mathrm{\Omega }}_{\mathrm{m}0}\ge 0.30$. This observation highlights the complex interplay between cosmological parameters and their implications for understanding the temporal evolution of fundamental physical constants.

5.1. ${\omega }_0=-1$ with Fixing ${\mathrm{\Omega }}_{m0}$

We keep ${\mathrm{\Omega }}_{\mathrm{m}0}$ constant and perform a ${\chi }^2$ test on these models. Notably, the value of h remains unchanged even as we vary ${\mathrm{\Omega }}_{\mathrm{m}0}$ within the range of $0.28$ to $0.32$. As ${\mathrm{\Omega }}_{\mathrm{m}0}$ increases, both $M_0$ and b decrease. In particular, the best-fit value of b is positive only for ${\mathrm{\Omega }}_{\mathrm{m}0}=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 show both positive and negative trends for ${\mathrm{\Omega }}_{\mathrm{m}0}=0.28$ and $0.29$. However, for $0.30\le {\mathrm{\Omega }}_{\mathrm{m}0}\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 Figure 1.

In panel a of Figure 1, we depict the cosmic evolution of the best-fit value of $c/c_0$ along with its $1-\sigma$ errors for the model with ${\omega }_0=-1$ and ${\mathrm{\Omega }}_{\mathrm{m}0}=0.29$. 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 the $1-\sigma$ error range, leading to ambiguity regarding the variation in the speed of light for this model.

For ${\mathrm{\Omega }}_{\mathrm{m}0}\ge 0.30$, the best-fit value of b is consistently negative within the $1-\sigma$ error range, indicating a monotonically decreasing speed of light over cosmic time. In panel b of Figure 1, we show the cosmological evolution of $c/c_0$ for the model with ${\omega }_0=-1$ and ${\mathrm{\Omega }}_{\mathrm{m}0}=0.30$. In this case, the best-fit value of b is $-0.027$, further supporting the decrease in the speed of light over time.

Additionally, the ratio of the time variation of the speed of light to its present value, expressed as ${\dot{c}}_0/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 ${\mathrm{\Omega }}_{\mathrm{m}0}=0.30$, 0.31, and 0.32, respectively. These constraints represent significant improvements over those reported in [78].

5.2. ${\omega }_0=-1$ without Fixing ${\mathrm{\Omega }}_{m0}$

Within the $1-\sigma$ error range, the cosmological parameters for the model with ${\omega }_0=-1$ are $-19.6605\le M_0\le -19.0851$, $0.188\le {\mathrm{\Omega }}_{\mathrm{m}0}\le 0.410$, $0.603\le h\le 0.787$, and $-0.218\le b\le 0.168$. Compared to the $\mathrm{\Lambda }$CDM model, the constraint on ${\mathrm{\Omega }}_{\mathrm{m}0}$ in the meVSL model is significantly weaker. It is because, as can be understood from Equation (25), both h and b contribute to $m_{\mathrm{th}}$, allowing the effect of changes in h to yield similar results from changes in b.

In the meVSL framework, the cosmological evolution of the speed of light is described by $c=c_0{(1+z)}^{-b/4}$. The best-fit value of b is $-0.025$, suggesting that the speed of light was higher in the past than it is today. However, the $1-\sigma$ range of b includes both negative and positive values, making it inconclusive to definitively determine the variation of the speed of light in this model.

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

We keep the value of h fixed and perform a maximum likelihood analysis. In these scenarios, both ${\mathrm{\Omega }}_{\mathrm{m}0}(=0.299\pm 0.111)$ and $b(=-0.025\pm 0.193)$ remain stable despite variations in h ranging from $0.6376$ to $0.74$. As h increases, $M_0$ also experiences a corresponding increase. The best-fit value of b remains consistent at $-0.025$ across these models, indicating a consistent decreasing trend in the speed of light over cosmic time. We also understand from Equation (25) that the variation in $M_0$ is insensitive to changes in h. However, within the $1-\sigma$ error margin, the values of b 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.

5.4. Fixing h

We perform a maximum likelihood analysis while varying other parameters under a fixed value of h. Under this condition, while we systematically vary h from $0.6736$ to $0.74$, the best-fit values of parameters $M_0$, ${\omega }_0$, and b all show an increasing trend. However, concurrently, the best-fit value of ${\mathrm{\Omega }}_{\mathrm{m}0}$ exhibits a decrease as h increases. Furthermore, as h increases, $M_0$ also shows an upward trend. Notably, the best-fit values of b consistently remain negative across these models, showing a consistent decrease in the speed of light over cosmic time. Nonetheless, within the $1-\sigma$ error range, the values of b demonstrate a variability, with some values being positive and others negative, suggesting no clear discernible pattern in the variation of the speed of light for these models.

5.5. Fixing ${\omega }_0$

We conduct a maximum likelihood analysis for $\omega$CDM models, with no fixed cosmological parameters except ${\omega }_0$. The range of ${\omega }_0$ is varied from $-0.9$ to $-1.1$. As ${\omega }_0$ decreases, ${\mathrm{\Omega }}_{\mathrm{m}0}$ decreases while b increases. However, consistent trends are not 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, while positive best-fit values of b are obtained for models with ${\omega }_0\le -1.0$. Nevertheless, 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.

5.6. Without Fixing

Finally, we conduct a maximum likelihood analysis without fixing any parameters. In this case, the best-fit value of $M_0$ decreases with increasing ${\mathrm{\Omega }}_{\mathrm{m}0}$, reaching $0.335$, while the values of ${\omega }_0$ and h decrease to $-1.22$ and $0.644$, respectively. These are shown in the last row of Table 1.

6. $\mathit{c}$ for CPL

In this section, we perform a maximum likelihood analysis for the CPL models, exploring various scenarios across different cosmological parameter values. The results of this analysis are summarized in

Table 2. Presented below are the best-fit values and their corresponding $1-\sigma$ errors for cosmological parameters in CPL models. Only models highlighted in green (or cyan) indicate the potential for exhibiting time variations in the speed of light. The b values highlighted in green indicate negative values, and those in cyan indicate positive values. These represent a clear monotonic decrease or increase in the speed of light with cosmic time, respectively.

Models	${\mathit{M}}_0$	${\mathit{\omega }}_0$	${\mathit{\omega }}_{\mathit{a}}$	${\mathbf{\Omega }}_{\mathit{m}0}$	h	b	$\mathit{\nu }$	${\mathit{\chi }}_{\mathit{\nu }}^2$
$b=0$	$-19.3617\pm 0.2914$	$-1$	$1.36\pm 0.42$	$0.163\pm 0.054$	$0.704\pm 0.094$	0	1044	$0.987$
	$-19.3582\pm 0.2915$	$-1.22\pm 0.15$	0	$0.348\pm 0.036$	$0.704\pm 0.094$	0	1044	$0.988$
Section 6.1	$-19.3325\pm 0.0908$	$-0.76\pm 0.02$	$1.30\pm 0.25$	$-0.062\pm 0.068$	$0.712\pm 0.030$	0	1043	$0.988$
$b\ne 0$	$-19.2231\pm 0.1985$	$-1$	$1.44\pm 0.08$	$-0.102\pm 0.065$	$0.748\pm 0.068$	$0.464\pm 0.156$	1043	$0.988$
Section 6.2	$-19.5513\pm 0.1284$	$-1.22\pm 0.04$	0	$0.335\pm 0.088$	$0.644\pm 0.038$	$0.022\pm 0.176$	1043	$0.989$
${\omega }_0=-0.95$	$-19.4115\pm 0.2313$		$1.85\pm 0.22$	$0.28$	$0.688\pm 0.073$	$-0.273\pm 0.014$	1044	$0.987$
	$-19.4065\pm 0.2324$		$1.90\pm 0.23$	$0.29$	$0.690\pm 0.074$	$-0.291\pm 0.014$	1044	$0.987$
	$-19.4110\pm 0.2328$	$-0.95$	$1.95\pm 0.23$	$0.30$	$0.689\pm 0.074$	$-0.309\pm 0.014$	1044	$0.987$
	$-19.4114\pm 0.2335$		$2.01\pm 0.24$	$0.31$	$0.689\pm 0.074$	$-0.326\pm 0.014$	1044	$0.987$
Section 6.3	$-19.4119\pm 0.2342$		$2.07\pm 0.24$	$0.32$	$0.688\pm 0.074$	$-0.344\pm 0.014$	1044	$0.987$
${\omega }_0=-1.0$	$-19.4108\pm 0.2322$		$1.73\pm 0.23$	$0.28$	$0.689\pm 0.074$	$-0.215\pm 0.015$	1044	$0.988$
	$-19.3848\pm 0.2349$		$1.79\pm 0.24$	$0.29$	$0.697\pm 0.075$	$-0.235\pm 0.015$	1044	$0.988$
	$-19.4036\pm 0.2344$	$-1.0$	$1.85\pm 0.24$	$0.30$	$0.691\pm 0.075$	$-0.255\pm 0.015$	1044	$0.988$
	$-19.4100\pm 0.2346$		$1.91\pm 0.24$	$0.31$	$0.689\pm 0.074$	$-0.274\pm 0.015$	1044	$0.988$
Section 6.4	$-19.4107\pm 0.2353$		$1.97\pm 0.25$	$0.32$	$0.689\pm 0.075$	$-0.293\pm 0.015$	1044	$0.988$
${\omega }_0=-1.05$	$-19.4038\pm 0.2339$		$1.57\pm 0.24$	$0.28$	$0.691\pm 0.074$	$-0.154\pm 0.015$	1044	$0.988$
	$-19.4076\pm 0.2344$		$1.62\pm 0.25$	$0.29$	$0.689\pm 0.074$	$-0.175\pm 0.015$	1044	$0.988$
	$-19.4079\pm 0.2351$	$-1.05$	$1.68\pm 0.25$	$0.30$	$0.689\pm 0.075$	$-0.195\pm 0.015$	1044	$0.988$
	$-19.4096\pm 0.2356$		$1.73\pm 0.25$	$0.31$	$0.689\pm 0.075$	$-0.213\pm 0.015$	1044	$0.988$
Section 6.5	$-19.4082\pm 0.2367$		$1.84\pm 0.26$	$0.32$	$0.690\pm 0.075$	$-0.240\pm 0.015$	1044	$0.988$
fixing ${\mathrm{\Omega }}_{m0}$	$-19.3499\pm 0.1176$	$-1.18\pm 0.04$	$0.72\pm 0.23$	$0.28$	$0.707\pm 0.038$	$0.038\pm 0.015$	1043	$0.988$
	$-19.3465\pm 0.1199$	$-1.18\pm 0.04$	$0.76\pm 0.25$	$0.30$	$0.708\pm 0.039$	$0.001\pm 0.015$	1043	$0.989$
Section 6.6	$-19.3478\pm 0.1225$	$-1.20\pm 0.04$	$0.71\pm 0.26$	$0.32$	$0.708\pm 0.040$	$-0.024\pm 0.015$	1043	$0.989$
fixing h	$-19.4544\pm 0.0072$	$-0.99\pm 0.02$	$1.31\pm 0.11$	$0.033\pm 0.061$	$0.6736$	$0.207\pm 0.120$	1043	$0.988$
Section 6.7	$-19.2475\pm 0.0078$	$-0.99\pm 0.01$	$1.45\pm 0.07$	$-0.112\pm 0.036$	$0.74$	$0.460\pm 0.094$	1043	$0.988$
No fixing	$-19.3500\pm 0.0857$	$-0.76\pm 0.02$	$1.30\pm 0.10$	$-0.054\pm 0.086$	$0.706\pm 0.028$	$-0.023\pm 0.108$	1042	$0.989$

. Initially, we examine the case of GR, represented by $b=0$, and then extend the analysis to include the meVSL models, where b is allowed to vary. Unlike the $\omega$CDM models, obtaining viable values for ${\mathrm{\Omega }}_{\mathrm{m}0}$ from this analysis without any prior constraints proves challenging. Therefore, we choose to exclusively perform the maximum likelihood analysis for fixed values of ${\mathrm{\Omega }}_{\mathrm{m}0}$.

6.1. $b=0$

We investigate CPL models within the framework of GR, initially setting $b=0$. We analyze three scenarios: ${\omega }_0=-1$, ${\omega }_a=0$, and without fixing ${\omega }_0$ and ${\omega }_a$. When ${\omega }_0=-1.0$, the range of ${\omega }_a$ spans $0.94\le {\omega }_a\le 1.78$, while ${\mathrm{\Omega }}_{\mathrm{m}0}$ ranges from $0.109$ to $0.217$ at the 68% confidence level. For the ${\omega }_a=0$ model, constraints at the 68% confidence level are $0.312\le {\mathrm{\Omega }}_{\mathrm{m}0}\le 0.384$, $-1.37\le {\omega }_0\le -1.07$, and $0.610\le \mathrm{h}\le 0.798$.

However, allowing ${\omega }_a$ to vary yields excessively small values of ${\mathrm{\Omega }}_{\mathrm{m}0}$, making them impractical as viable models. Permitting both ${\omega }_0$ and ${\omega }_a$ to vary, the ranges of cosmological parameters at the 68% confidence level are $-0.78\le {\omega }_0\le -0.74$, $1.05\le {\omega }_a\le 1.55$, and $-0.13\le {\mathrm{\Omega }}_{\mathrm{m}0}\le 0.006$. Negative values of ${\mathrm{\Omega }}_{\mathrm{m}0}$ arise from statistical analysis, but physically, negative mass is not meaningful in the SMC. Therefore, we exclude consideration of this model.

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

6.2. Fixing ${\omega }_0=-1$ or ${\omega }_a=0$

First, fixing ${\omega }_0=-1$ and allowing other variables to vary, we obtain the values of $M_0$, $\omega$, ${\mathrm{\Omega }}_{\mathrm{m}0}$, h, and b within the $1-\sigma$ confidence level as ($-19.4216,-19.0246$), ($-1.36,-1.52$), ($-0.167,-0.037$), ($0.68,0.816$), and ($0.308,0.62$), respectively. Since this model also yields negative ${\mathrm{\Omega }}_{\mathrm{m}0}$, we exclude consideration of this model.

Next, fixing ${\omega }_a=0$ and allowing other variables to vary, we conduct a maximum likelihood analysis. In this case, we obtain the values of $M_0$, $\omega$, ${\mathrm{\Omega }}_{\mathrm{m}0}$, h, and b within the $1-\sigma$ confidence level as ($-19.6797,-19.4229$), ($-1.26,-1.18$), ($0.247,0.423$), ($0.606,0.682$), and ($-0.154,0.198$), respectively. This model exhibits a relatively small value of $M_0$ compared to other models.

6.3. ${\omega }_0=-0.95$ with Fixing ${\mathrm{\Omega }}_{m0}$

We perform a maximum likelihood analysis for models with ${\omega }_0=-0.95$, while varying ${\mathrm{\Omega }}_{\mathrm{m}0}$ from $0.28$ to $0.32$. Across these models, the best-fit values of h consistently hover around $0.69$, with a $1-\sigma$ error margin of $0.07$. With increasing ${\mathrm{\Omega }}_{\mathrm{m}0}$, the best-fit value of ${\omega }_a$ also rises, whereas the best-fit values of b decline.

At a 68% confidence level, all values of b fall into the negative range for the specified ${\mathrm{\Omega }}_{\mathrm{m}0}$ range, spanning from $-0.344$ to $-0.273$. This suggests that the speed of light decreases monotonically over cosmic time, with its rate of decrease accelerating as ${\mathrm{\Omega }}_{\mathrm{m}0}$ increases.

At $z=3$, the speed of light can exceed $c_0$ by approximately ${10}^{13}$% when ${\mathrm{\Omega }}_{\mathrm{m}0}$ ranges from $0.28$ to $0.32$. The $1-\sigma$ ranges of ${\dot{c}}_0/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.

6.4. ${\omega }_0=-1.0$ with Fixing ${\mathrm{\Omega }}_{m0}$

We conduct an extensive analysis focusing on ${\omega }_0=-1.0$ models. Across these models, the best-fit values of h range from $0.689$ to $0.697$ for the specified values of ${\mathrm{\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 with increasing ${\mathrm{\Omega }}_{\mathrm{m}0}$, with all b values falling into the negative range at a 68% confidence level. This indicates a monotonically decreasing trend in the speed of light over 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$, the speed of light exceeds $c_0$ by approximately $8\times 10$% when ${\mathrm{\Omega }}_{\mathrm{m}0}$ ranges from $0.28$ to $0.32$. Furthermore, the $1-\sigma$ ranges of ${\dot{c}}_0/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.

6.5. ${\omega }_0=-1.05$ with Fixing ${\mathrm{\Omega }}_{m0}$

In this subsection, we perform a maximum likelihood analysis for models with ${\omega }_0=-1.05$. Across these models, the best-fit values of h consistently hover around $0.69$ for all specified values of ${\mathrm{\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$. Similar to previous models, the values of b decrease with increasing ${\mathrm{\Omega }}_{\mathrm{m}0}$, ranging from $-0.169$ to $-0.139$ for ${\mathrm{\Omega }}_{\mathrm{m}0}=0.28$ and from $-0.255$ to $-0.225$ for ${\mathrm{\Omega }}_{\mathrm{m}0}=0.32$ within a $1-\sigma$ error. Thus, the speed of light continues its monotonic decrease over cosmic time in these models. At $z=3$, c can exceed $c_0$ by approximately $6\%$ ($9\%$) for ${\mathrm{\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{c}}_0/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.

6.6. CPL with Fixing ${\mathrm{\Omega }}_{m0}$

We analyze the Pantheon data without constraining values of ${\omega }_0$ and ${\omega }_a$ within the range $0.28\le {\mathrm{\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 ${\mathrm{\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 ${\mathrm{\Omega }}_{\mathrm{m}0}=0.28,0.30$, and $0.32$, respectively. Of particular interest is the model with ${\mathrm{\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 ${\mathrm{\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 Figure 2. In the left panel, the cosmological evolution of $c/c_0$ as a function of z for the ${\mathrm{\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 $c/c_0$ with increasing redshift is evident due to the positive values of b. At $z=3$, c decreases by approximately $0.8$ ($1.8$)% within $1-\sigma$ error. The $1-\sigma$ range of ${\dot{c}}_0/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 Figure 2 illustrates the model with ${\mathrm{\Omega }}_{\mathrm{m}0}=0.32$, where both the best-fit and $1-\sigma$ error values of b are negative. Consequently, $c/c_0$ increases monotonically with redshift. At $z=3$, c increases by about $0.3$ ($1.4$)% at a 68% confidence level. The $1-\sigma$ range of ${\dot{c}}_0/c_0$ (${10}^{-13}{\mathrm{yr}}^{-1}$) spans from $(-7.05,-1.63)$.

6.7. CPL with (without) Fixing h

The analysis is conducted with fixed values of h while leaving other parameters unconstrained. For h $=0.6736\left(0.74\right)$, 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 ${\mathrm{\Omega }}_{\mathrm{m}0}$ fall outside of the viable range, yielding $-0.028\le {\mathrm{\Omega }}_{\mathrm{m}0}\le 0.094(-0.148\le {\mathrm{\Omega }}_{\mathrm{m}0}\le -0.076)$ for $h=0.6736\left(0.74\right)$, which are deemed nonviable. Moreover, without constraining ${\mathrm{\Omega }}_{\mathrm{m}0}$, the matter density contrast is estimated to be $-0.14\le {\mathrm{\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.

7. $\dot{\mathit{c}}$ and $\dot{\mathit{G}}$

Expanding on the insights gained from the previous Section 5 and Section 6, we embark on an exploration of viable meVSL models across diverse dark energy scenarios, aiming to extract valuable constraints on both cosmological and model parameters. This analysis allows us to derive estimates for the temporal evolution of the speed of light, a fundamental aspect within the meVSL framework. Notably, in the meVSL paradigm, the speed of light undergoes cosmological evolution, mirroring the behavior of the gravitational constant, which is characterized by $G=G_0{(1+z)}^{-b}$. By 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{G}}_0/G_0$.

To place these findings within the broader context of observational constraints, we compare our results with existing bounds on ${\dot{G}}_0/G_0$ from various sources as summarized in

Table 3. Table provides the latest $1-\sigma$ observational constraints on the present rate of change of the gravitational constant, ${\dot{G}}_0/G_0$. Here, “WD” refers to white dwarf observations, “BBN” signifies Big Bang nucleosynthesis, “LLR” denotes lunar laser ranging data, and “GWs” represents gravitational waves.

obs		${\dot{\mathit{G}}}_0/{\mathit{G}}_0\left({10}^{-12}{\mathbf{yr}}^{-1}\right)$	Ref.
pulsars		23	[84]
WD	cooling	−1.8	[85]
	pulsation	−130	[86]
BBN		−0.3∼0.4	[87]
		−3.6∼4.5	[88]
LLR		−0.5∼0.9	[89]
		−0.005∼0.147	[90]
SNe Ia		−30∼73	[91]
		3	[92]
GWs	LIGO	70	[93]
	LISA	$0.7$	[94]

. Notably, the analysis of lunar laser ranging (LLR) data stands out for its stringency, yielding the most stringent bounds on ${\dot{G}}_0/G_0$. In contrast, the orbital period rate of pulsars offers the widest bounds, estimated at $2.3\times {10}^{-11}{\mathrm{yr}}^{-1}$. Taken together, these observations suggest that ${\dot{G}}_0/G_0$ is tentatively within the order of ${10}^{-12}{\mathrm{yr}}^{-1}$. Such insights not only deepen our understanding of cosmological dynamics but also pave the way for future investigations into the fundamental nature of physical constants.

Table 4. Time variations of the speed of light and that of the gravitational constant at a 68% confidence level for viable models. We denote $A\equiv {10}^{-13}{\mathrm{yr}}^{-1}$.

${\mathit{\omega }}_0$	${\mathit{\omega }}_{\mathit{a}}$	${\mathbf{\Omega }}_{\mathbf{m}0}$	$\mathbf{\Delta }\mathit{c}(\mathit{z}=3)(\%)$	${\dot{\mathit{c}}}_0/{\mathit{c}}_0\left[\mathit{A}\right]$	$\mathbf{\Delta }\mathit{G}(\mathit{z}=3)(\%)$	${\dot{\mathit{G}}}_0/{\mathit{G}}_0\left[10\mathit{A}\right]$	${\mathit{M}}_0$	b
$-1$	0	$0.30$	0.2∼1.7	$-8.76$∼$-0.89$	$0.7\sim 7.0$	$-3.51$∼$-0.36$	$-19.3561\pm 0.2899$	$-0.027\pm 0.022$
		$0.31$	0.8∼2.3	$-1.18$∼$-0.39$	3.1∼9.6	$-4.72$∼$-1.57$	$-19.3563\pm 0.2899$	$-0.044\pm 0.022$
		$0.32$	1.4∼2.9	$-1.48$∼$-0.70$	5.6∼12.2	$-5.94$∼$-2.79$	$-19.3566\pm 0.2898$	$-0.061\pm 0.022$
$-0.95$	$1.85\pm 0.22$	$0.28$	9.4∼10.5	$-51.3$∼$-46.3$	43.2∼48.9	$-20.5$∼$-18.5$	$-19.4115\pm 0.2313$	$-0.273\pm 0.014$
	$1.95\pm 0.23$	$0.30$	$10.8$∼$11.8$	$-57.8$∼$-52.8$	$50.5$∼$56.5$	$-23.1$∼$-21.1$	$-19.4110\pm 0.2328$	$-0.309\pm 0.014$
	$2.07\pm 0.24$	$0.32$	$12.1$∼$13.2$	$-64.0$∼$-59.0$	$58.0$∼$64.3$	$-25.6$∼$-23.6$	$-19.4119\pm 0.2342$	$-0.344\pm 0.014$
$-1.0$	$1.73$±$0.23$	$0.28$	$7.2$∼$8.3$	$-41.1$∼$-35.8$	$32.0$∼$37.6$	$-16.5$∼$-14.3$	$-19.4108\pm 0.2322$	$-0.215\pm 0.015$
	$1.85\pm 0.24$	$0.30$	$8.7$∼$9.8$	$-48.3$∼$-42.9$	$39.5$∼$45.4$	$-19.3$∼$-17.2$	$-19.4036\pm 0.2344$	$-0.255\pm 0.015$
	$1.97\pm 0.25$	$0.32$	$10.1$∼$11.3$	$-55.1$∼$-49.7$	$47.0$∼$53.2$	$-22.0$∼$-19.9$	$-19.4107\pm 0.2353$	$-0.293\pm 0.015$
$-1.05$	$1.57\pm 0.24$	$0.28$	$4.9$∼$6.0$	$-30.2$∼$-24.9$	$21.3$∼$26.4$	$-12.1$∼$-9.9$	$-19.4038\pm 0.2339$	$-0.154\pm 0.015$
	$1.68\pm 0.25$	$0.30$	$6.4$∼$7.5$	$-37.6$∼$-32.2$	$28.3$∼$33.8$	$-15.0$∼$-12.9$	$-19.4079\pm 0.2351$	$-0.195\pm 0.015$
	$1.84\pm 0.26$	$0.32$	$8.1$∼$9.2$	$-45.6$∼$-40.2$	$36.6$∼$42.4$	$-18.2$∼$-16.1$	$-19.4082\pm 0.2367$	$-0.240\pm 0.015$
$-1.18\pm 0.04$	$0.72\pm 0.23$	$0.28$	$0.8$∼$1.8$	$4.11$∼$9.48$	$3.1$∼$7.1$	$1.65$∼$3.79$	$-19.3499\pm 0.1176$	$0.038\pm 0.015$
$-1.20\pm 0.04$	$0.71\pm 0.26$	$0.32$	$0.3$∼$1.4$	$-6.98$∼$-1.61$	$1.3$∼$5.6$	$-2.79$∼$-0.64$	$-19.3478\pm 0.1225$	$-0.024\pm 0.015$

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 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 c(z=3)\equiv \left(c(z=3)-c_0\right)/c_0\times 100(\%)$ within a $1-\sigma$ uncertainty range. Similarly, $\Delta 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{c}}_0/c_0$, while ${\dot{G}}_0/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 ${\mathrm{\Omega }}_{\mathrm{m}0}=0.28$. Consequently, both ${\dot{c}}_0/c_0$ and ${\dot{G}}_0/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{c}}_0/c_0$ and ${\dot{G}}_0/G_0$ in these scenarios.

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

Despite having three different constraints on ${\dot{G}}_0/G_0$ derived from SNe Ia data, including this work, there are notable discrepancies among the results [91,92]. These discrepancies stem from the varying methodologies and datasets used. For instance, one approach employs the standard candle method to define the redshift–distance relation, establishing a broad upper limit on $|\dot{G}/G|$ and exploring various parameterizations and impacts [91]. In contrast, another method focuses on the intrinsic properties of Type Ia supernovae (SNIa) and their dependency on the Chandrasekhar mass $M_{\mathrm{Ch}}\propto G^{-3/2}$ to track the variation of G across different redshifts [92]. Expanding on the insights gained from the previous Section 3 and Section 5, we embark on an exploration of viable meVSL models across diverse dark energy scenarios, aiming to extract valuable constraints on both cosmological and model parameters. This analysis allows us to derive estimates for the temporal evolution of the speed of light, a fundamental aspect within the meVSL framework. Notably, in the meVSL paradigm, the speed of light undergoes cosmological evolution, mirroring the behavior of the gravitational constant, which is characterized by $G=G_0{(1+z)}^{-b}$. By 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{G}}_0/G_0$.

8. Discussion

The Pantheon+ data provide constraints on cosmological and model parameters with a statistical precision of about 10%. Leveraging this dataset, we perform a maximum likelihood analysis to constrain dark energy models within the framework of the modified varying speed of light (meVSL) model. This analysis allows us to 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 constraints obtained from our analysis indicate that the relative temporal variations of the speed of light and the gravitational constant lie within the ranges $-64.0\le {\dot{c}}_0/c_0\left({10}^{-13}{\mathrm{yr}}^{-1}\right)\le -0.39$ and $-25.6\le {\dot{G}}_0/G_0\left({10}^{-12}{\mathrm{yr}}^{-1}\right)\le -0.36$, respectively, for most viable models. These findings suggest that, according to the current Pantheon data, both the speed of light and the gravitational constant were greater in the past and have decreased monotonically with redshift z.

Among the CPL models we consider, some do not require a dark matter energy density. We understand this phenomenon as a result of the degeneracy between ${\omega }_a$ and ${\mathrm{\Omega }}_{\mathrm{m}0}$. We focus our viable model consideration on meVSL models where the cosmological parameters closely match those of the Standard Cosmological Model (SCM). However, models with ${\mathrm{\Omega }}_{\mathrm{m}0}$ approaching zero could become an interesting topic of study, requiring a deeper understanding and further research. For the purposes of this paper, we defer to the references cited within the main text to substantiate these points.

While additional cosmological observations, such as those from the Cosmic Microwave Background (CMB) and Baryon Acoustic Oscillations (BAO), could provide further constraints on cosmological and model parameters, integrating these datasets would necessitate a reanalysis within the theoretical framework of the meVSL model. This comprehensive task falls outside the scope of the present manuscript and is deferred to future investigations. Future studies incorporating CMB and BAO data could significantly refine the constraints on the parameters governing the evolution of physical constants, thereby enhancing our understanding of the meVSL model and its implications for cosmology.