© The Authors 2026

1. Introduction

General relativity and the assumption of homogeneity and isotropy of the Universe at large scales are two of the foundations of modern cosmology. The latter leads to the Friedmann–Lemaître–Robertson–Walker (FLRW) metric, under which Einstein’s field equations reduce to the Friedmann equations. The first, combined with the FLRW metric, provides a formula for the luminosity distance–redshift relation, d_L(z), in terms of three parameters: the Hubble constant (H₀), the matter density (Ω_M), and the cosmological constant (Λ) with density parameter Ω_Λ.

Using Hubble diagrams of Type Ia supernovae (SNe Ia), Riess et al. (1998) and Perlmutter et al. (1999) found Ω_Λ >  0 and a deceleration parameter q₀ < 0, indicating that the Universe’s expansion is currently accelerating. In this context, Λ is interpreted as a constant dark energy density that drives cosmic acceleration. Given that models with Λ, when fitted to various observational data, yield Ω_Λ >  0 (Weinberg et al. 2013), the existence of dark energy is widely accepted, consolidating the Λ cold dark matter (ΛCDM) model as the standard model.

Despite its ability to fit observations, the ΛCDM model has problems of fine-tuning, cosmic coincidence, and tensions in cosmological parameters measured with independent experiments (Perivolaropoulos & Skara 2022; Efstathiou 2025). These issues motivate alternative gravitational theories that explain cosmic acceleration without Λ (Koyama 2016; Odintsov et al. 2025). Moreover, cosmic acceleration may be partly an apparent effect associated with the assumption of homogeneity and isotropy and the corresponding use of the FLRW metric (Räsänen 2006; Enqvist 2008; Wiltshire 2009).

Cosmography provides a model-independent framework to study cosmic acceleration, expressing d_L(z) in terms of spatial curvature and kinematic parameters, such as H₀, q₀, and the jerk parameter (j₀). However, this approach relies on the assumption of a homogeneous and isotropic Universe (Hu & Wang 2022).

Among current SN Ia samples, such as Pantheon+ (Brout et al. 2022) and the full five-year sample from the Dark Energy Survey Supernova (DES-SN) program (DES-SN5YR; Sánchez et al. 2024), a significant fraction of SNe originates from deep-field surveys. Analysing deep fields separately allows one to infer direction-dependent parameters without requiring global homogeneity and isotropy. In particular, comparing deep fields across different directions provides a test of isotropy.

Since angular variations in deep fields are expected to be negligible, and assuming that the large-scale density within each field is independent of the comoving radial coordinate, the metric can be approximated as FLRW on a field-by-field basis. Under this assumption, we can measure q₀ for each deep field using the FLRW metric, without assuming global homogeneity and isotropy. The main limitation of this approach is the small number of low-z SNe Ia (z < 0.1) in deep fields, which are crucial for breaking parameter degeneracies (Linder 2006).

An approach to reducing the impact of the lack of low-z SNe on the estimation of q₀ is to use a d_L(z) relation with as few parameters as possible. In the flat ΛCDM model, d_L(z) depends only on H₀ and Ω_M. The latter, which is used to measure q₀, remains relatively stable in the absence of low-z SNe (Brout et al. 2022). To date, the only other d_L(z) relation that accurately fits SN observations with just two parameters is that proposed by Hu et al. (2024a). This relation corresponds to a third-order Padé approximation with H₀ and q₀ as free parameters, while j₀ is fixed to one, as predicted by the flat ΛCDM model (Bochner et al. 2015).

In this work, we introduce an empirical magnitude–redshift relation (m(z)) and the corresponding d_L(z) relation, which provides a straightforward quantitative tool to fit Hubble diagrams of SNe Ia up to at least z ≃ 1.1, using only two parameters and without assuming any theoretical model or spatial curvature. This relation is relatively stable in the absence of low-z SNe, making it suitable for studying SN deep fields without the need to add a low-z sample.

The paper is organised as follows. In Sect. 2, we describe the SN samples and define the deep-field regions. In Sect. 3, we present the m(z) relations used to fit Hubble diagrams. In Sect. 4, we present the results, including the new empirical relation, its comparison with five m(z) relations, and its application to the deep-field regions and other samples. In Sect. 5, we discuss the derived j₀ values, the search for cosmic anisotropies, and future analyses. We summarise our conclusions in Sect. 6.

2. Data samples

Pantheon+¹ contains data from 1550 SNe Ia, drawn from 18 surveys, with 0.001 ≤ z ≤ 2.261. DES-SN5YR² contains data from 1635 photometrically classified SNe Ia from the DES-SN programme with 0.060 ≤ z ≤ 1.121, and 194 SNe Ia from four surveys with 0.025 < z < 0.093, referred to as the low-z sample.

The Pantheon+ and DES-SN5YR samples provide (i) RA and Dec coordinates; (ii) standardised SN magnitudes (m_B^corr, hereafter m), derived from SALT2 (Pantheon+) and SALT3 (DES-SN5YR) light-curve fits and corrected for stretch, colour, host-galaxy mass, and selection bias; (iii) redshifts corrected for the cosmic microwave background (CMB) dipole and peculiar velocities (z_HD, hereafter z); and (iv) the covariance matrix (C) to account for statistical and systematic errors.

Following Brout et al. (2022), we selected SNe with z >  0.01 from the Pantheon+ compilation to minimise the impact of peculiar velocities on the derived cosmological parameters. We refer to this subsample as Pantheon+. Because DES-SN5YR introduced several improvements compared with Pantheon+, its m values are not on the same scale. Specifically, there is a constant offset of 0.04 mag between the selection bias corrections in DES-SN5YR and Pantheon+, which does not affect the cosmological results of each separate analysis (Vincenzi et al. 2025). Hence, we analysed both samples separately.

Fig. 1 shows the sky distribution of SNe in Pantheon+ and DES-SN5YR. The SNe with z ≥ 0.4 are concentrated in small regions of the sky, which reflects the location of the SN deep fields. We defined 14 circular regions for Pantheon+ and four for DES-SN5YR, referred to as XY, which contain virtually all SNe with z ≥ 0.4. In this notation, X denotes the closest Galactic pole and Y indicates the rank of proximity to that pole. Table 1 lists the regions, their centres in (RA, Dec), angular radii (θ), the number of SNe within each region, and the corresponding z ranges.

Fig. 1. Sky distribution of Pantheon+ and DES-SN5YR SNe. Dashed lines indicate the Galactic equator. Plus and cross symbols mark the north and south Galactic poles, respectively. Circular and triangular regions are also shown.

The S1, S2, S4, and S5 regions in Pantheon+ closely match those in DES-SN5YR. In these regions, DES-SN5YR contains 170–490 more SNe than Pantheon+, while 47–100% of the SNe in the Pantheon+ S1, S2, S4, and S5 regions originate from the DES-SN programme and are already included in DES-SN5YR. Therefore, we used only the DES-SN5YR data for these regions.

The remaining ten Pantheon+ regions contain 13–77 SNe each. To increase the sample while minimising the sky area, we defined the triangular N147, N235, and S367 regions. The N147 region includes N1, N4, N7, and the SNe located within the triangle formed by their centres; the same applies for N235 and S367. For S367, we increased the declination of the centres of S3 and S7 by 1.2° to include more SNe without significantly increasing the area. Fig. 1 shows the triangular regions, while Table 1 lists their barycentres, the number of SNe, and the z ranges. We refer to the S1, S2, S4, S5, N147, N235, S367, and N6 regions as the deep-field regions.

3. Methods

3.1. Magnitude–redshift relations

The m(z) relation for a standard candle with absolute magnitude M is given by

$$\begin{array}{c}\hfill m(z)=\mathcal{M}+5log{\mathcal{D}}_L(z),\end{array}$$(1)

where ℳ = M + 5log(c H₀⁻¹ Mpc⁻¹)+25 and 𝒟_L(z) =d_L(z)H₀/c. To derive d_L(z) using either a gravitational theory or a cosmographic expansion, we need the metric of the expanding Universe. The simplest choice is the FLRW metric, characterised by the scale factor a and the curvature parameter k. For this metric,

$$\begin{array}{c}\hfill d_L(z)=(1+z)d_M(z),\end{array}$$(2)

where d_M(z) is the transverse comoving distance. Defining H(z) = ȧ/a and Ω_k = −kc²/H₀², d_M(z) is given by

$$\begin{array}{c}\hfill d_M(z)=\frac{c}{H_0}\{\begin{array}{cc}{\mathrm{\Omega }}_k^{-1/2}sinh({\mathrm{\Omega }}_k^{1/2}{\displaystyle {\int }_0^z}\frac{H_0}{H(z\prime )}\mathrm{d}z\prime )\hfill & {\mathrm{\Omega }}_k>0\hfill \\ {\displaystyle {\int }_0^z}\frac{H_0}{H(z\prime )}\mathrm{d}z\prime \hfill & {\mathrm{\Omega }}_k=0\hfill \\ |{\mathrm{\Omega }}_k{|}^{-1/2}sin(|{\mathrm{\Omega }}_k{|}^{1/2}{\displaystyle {\int }_0^z}\frac{H_0}{H(z\prime )}\mathrm{d}z\prime )\hfill & {\mathrm{\Omega }}_k<0\hfill \end{array}.\end{array}$$(3)

3.1.1. The ΛCDM, flat ΛCDM, and flat wCDM models

In the phenomenological generalisation of the ΛCDM model,

$$\begin{array}{c}\hfill H(z)=H_0\sqrt{{\mathrm{\Omega }}_M{(1+z)}^3+{\mathrm{\Omega }}_k{(1+z)}^2+{\mathrm{\Omega }}_{\mathrm{\Lambda }}{(1+z)}^{3(1+w)}},\end{array}$$(4)

where Ω_M + Ω_Λ + Ω_k = 1 and w is the dark energy equation-of-state parameter. We used the m(z) relations for ΛCDM (w = −1), flat ΛCDM (Ω_k = 0, w = −1), and the phenomenological flat wCDM (Ω_k = 0) models.

3.1.2. Padé cosmography and flat Padé model

The Padé approximation is a method to approximate a function F(z) as the ratio of two polynomials of orders A and B. The rational function is known as the (A,B) Padé approximant, and its coefficients are expressed in terms of those of the Taylor series of F(z) truncated at order A + B (Baker & Graves-Morris 1996).

Hu & Wang (2022), using Pantheon (Scolnic et al. 2018), report that the (2,1) Padé approximant for d_L(z) performs better than cosmographic expansions based on Taylor series. More recently, Hu et al. (2024a) analysed Pantheon+ and found that this approximant outperforms other Padé forms. We therefore adopt it as the best cosmographic approach.

The (2,1) Padé approximant for d_L(z) used by Hu et al. (2024a) was derived by Capozziello et al. (2020) for Ω_k = 0. To compute the expression valid for any Ω_k, we used the third-order Taylor series for d_L(z) (Eq. 2 of Cattoën & Visser 2007) and the procedure of Baker & Graves-Morris (1996), obtaining

$$\begin{array}{c}\hfill d_L(z)=\frac{\mathit{cz}}{H_0}\frac{6(1-q_0)+(5-8q_0-3q_0^2+2\widehat{j_0})z}{6(1-q_0)+2(1-q_0-3q_0^2+\widehat{j_0})z}.\end{array}$$(5)

Here, $\widehat{j_0}\equiv j_0-{\mathrm{\Omega }}_k$, and

$$\begin{array}{c}\hfill q_0={\frac{\mathrm{d}lnH(z)}{\mathrm{d}z}|}_{z=0}-1,j_0=\frac{1}{H_0}{\frac{{\mathrm{d}}^{2}H(z)}{\mathrm{d}{z}^2}|}_{z=0}+q_0^2.\end{array}$$(6)

For Ω_k = 0, Eq. (5) is equivalent to that given in Hu et al. (2024a). We refer to Eq. (5) as the Padé cosmography.

Hu et al. (2024a) also found that the (2,1) Padé approximant for d_L(z) with j₀ = 1 fixed performs better than the Padé cosmography, the flat ΛCDM, ΛCDM, and flat wCDM models. Since j₀ was fixed to unity because this is the value for the flat ΛCDM model, the d_L(z) equation proposed by Hu et al. (2024a) is not cosmographic, but rather model-based. We refer to Eq. (5) assuming $\widehat{j_0}=1$ as the flat Padé model.

3.1.3. Empirical approach

To fit the Hubble diagram, we introduce the ansatz

$$\begin{array}{c}\hfill d_L(z)=\frac{\mathit{cz}}{H_0}(1+z){10}^{f(z)/5},\end{array}$$(7)

where the factor 1 + z is inspired by Eq. (2) and f(z) is a free function that vanishes as z → 0, ensuring d_L(z) =cz/H₀ at low z. The corresponding m(z) relation is

$$\begin{array}{c}\hfill m(z)=\mathcal{M}+5log(z(1+z))+f(z),\end{array}$$(8)

where the dependence of f(z) on z can be empirically determined from the correlation between m − 5log(z(1 + z)) and z.

3.2. Parameter estimation and model selection

Let v denote the vector with the n_p free parameters of m(z). The posterior probability of m(z) is P ∝ p(v)ℒ, where p(v) is the prior and ℒ = e^−χ2/2. Following Conley et al. (2011),

$$\begin{array}{c}\hfill {\chi }^2=\mathrm{\Delta }{\mathbf{m}}^T{\mathbf{C}}^{-1}\mathrm{\Delta }\mathbf{m},\end{array}$$(9)

where Δm is the vector of residuals, whose i-th component is

$$\begin{array}{c}\hfill \mathrm{\Delta }{m}_i=m_i-m(z_i,\mathbf{v}).\end{array}$$(10)

We assumed non-informative priors; in this case, we obtained the best-fit parameter vector v_best that maximises P by minimising χ².

To calculate uncertainties, we used EMCEE (Foreman-Mackey et al. 2013), which samples P using a Markov chain Monte Carlo (MCMC) process. First, we defined the priors as flat distributions: ℳ ∈ (23, 25), Ω_M ∈ (0, 1), Ω_Λ ∈ (0, 1), w ∈ ( − 2, 0), q₀ ∈ ( − 2, 1), and $\widehat{j_0}\in (-4,7)$. These priors are wide enough to consider them as uninformative. Next, we initialised 10n_p walkers in a small Gaussian ball around v_best and ran 10⁵ steps (10⁶ for the flat wCDM model). Then, given an integrated autocorrelation time (τ) provided by EMCEE, we discarded the initial 3τ steps as burn-in and thin by τ/2. With this configuration, we obtain for all parameters a Gelman & Rubin (1992) statistic of $\widehat{R}\approx 1.0$, indicating convergence of the MCMC chains, and an effective sample size n_eff >  2000 (n_eff >  400 for the flat wCDM model). Finally, for each parameter, we adopted the 68.27% confidence interval of its marginalised distribution as the 1σ uncertainty.

To identify which model among a set of candidates best describes the Hubble diagram, we used the corrected Akaike information criterion (AICc; Sugiura 1978) and the Bayesian information criterion (BIC; Schwarz 1978). For each model, we computed AICc = χ_min² + 2n_pN/(N − n_p − 1) and BIC = χ_min² + n_plnN, where χ_min² = χ²(v_best) and N is the number of SNe. Based on these information criteria (ICs), the preferred model is the one with the lowest IC value (IC_min), and the strength of evidence for each model is given by Δ_IC = IC − IC_min. Models with Δ_AICc < 2 have substantial support, while those with 4 < Δ_AICc < 7 have considerably less support (Burnham & Anderson 2002). Similarly, Δ_BIC < 2 indicates little evidence against the model, and 6 < Δ_BIC < 10 indicates strong evidence against the model (Kass & Raftery 1995). Therefore, models with both Δ_AICc < 2 and Δ_BIC < 2 perform comparably, yielding similar fits to the Hubble diagram. Models with larger Δ_AICc are progressively less supported, and those with larger Δ_BIC are increasingly disfavoured.

4. Results

4.1. Empirical relation

Figure 2 shows m − 5log(z(1 + z)) versus z for Pantheon+ and DES-SN5YR. Their weighted Pearson correlation coefficients (r_P; Earp et al. 2019) are close to −0.7, with p-values below 10⁻²²⁰, indicating a moderate-to-strong negative linear correlation. Based on this, we express m − 5log(z(1 + z)) as a linear function of z, with an intercept ℳ and a slope b:

$$\begin{array}{c}\hfill m(z)=\mathcal{M}+bz+5log(z(1+z)).\end{array}$$(11)

Fig. 2. Relation between m − 5log(z(1 + z)) and z for Pantheon+ and DES-SN5YR. Binned data (blue squares) with error bars are shown for visualisation purposes only, both here and throughout the paper.

Thus, f(z) =bz in Eq. (8), which, together with Eq. (7) yields

$$\begin{array}{c}\hfill d_L(z)=\frac{\mathit{cz}}{H_0}(1+z){10}^{bz/5}.\end{array}$$(12)

We refer to Eqs. (11) and (12) as the empirical relation.

The empirical relation exhibits several characteristics. First, it is relatively stable in the absence of low-z SNe (see Sect. 4.5). Second, its parameters and their associated covariance matrix can be computed analytically. Fig. 3 shows the confidence contours for ℳ and b computed analytically and using EMCEE (using the flat prior b ∈ ( − 2, 0)), which are virtually identical. Third, under the FLRW metric, b ∝ q₀ + 1. Indeed, from Eqs. (12) and (2),

$$\begin{array}{c}\hfill d_M(z)=\frac{\mathit{cz}}{H_0}{10}^{bz/5},\end{array}$$(13)

Fig. 3. Confidence contours at the 68.27% and 95.45% levels of the empirical relation parameters, obtained analytically and using EMCEE.

and using Eq. (16) of Li et al. (2020),

$$\begin{array}{c}\hfill H(z)=H_0\frac{\sqrt{{10}^{-2bz/5}+{\mathrm{\Omega }}_k{z}^2}}{1+bln(10)z/5},\end{array}$$(14)

which, together with Eq. (6), yields

$$\begin{array}{c}\hfill q_0=-\frac{2}{5}ln(10)b-1.\end{array}$$(15)

4.2. Constraints on cosmological parameters

Table 2 lists the parameter constraints from Pantheon+ and DES-SN5YR for the six m(z) relations used in this work. The parameters for the flat ΛCDM, ΛCDM, and flat wCDM models are consistent with those reported by the Pantheon+ team (Ω_M = 0.334 ± 0.018 for the flat ΛCDM model, Ω_M = 0.306 ± 0.057 and Ω_Λ = 0.625 ± 0.084 for the ΛCDM model, and ${\mathrm{\Omega }}_M=0.{309}_{-0.069}^{+0.063}$ and w = −0.90 ± 0.14 for the flat wCDM model; Brout et al. 2022) and the DES Collaboration (Ω_M = 0.352 ± 0.017 for the flat ΛCDM model, ${\mathrm{\Omega }}_M=0.{291}_{-0.065}^{+0.063}$ and Ω_Λ = 0.55 ± 0.10 for the ΛCDM model, and ${\mathrm{\Omega }}_M=0.{264}_{-0.096}^{+0.074}$ and $w=-0.{80}_{-0.16}^{+0.14}$ for the flat wCDM model; Camilleri et al. 2024).

For each m(z) relation, the parameters derived from DES-SN5YR and Pantheon+ are consistent within 1σ, except for ℳ, which differs by more than 3.3σ. This discrepancy arises from the constant offset of 0.04 mag between selection bias corrections in DES-SN5YR and Pantheon+ (Vincenzi et al. 2025). After accounting for this offset, the ℳ values derived from DES-SN5YR and Pantheon+ become consistent within 1σ.

4.3. Hubble diagrams and model comparison

Fig. 4 shows the Hubble diagrams for DES-SN5YR and Pantheon+, the best fits for the m(z) relations, and the residuals relative to the empirical relation. For DES-SN5YR, the empirical relation provides a fit comparable to the other m(z) relations. The differences between the m(z) relations and the empirical relation have rms values below 0.009 mag and maximum absolute deviations Δm_max below 0.027 mag. For Pantheon+, the empirical relation begins to diverge from the other m(z) relations at z ≳ 1.4. For z < 1.121 (the highest z in DES-SN5YR), the rms and Δm_max remain below 0.007 and 0.012 mag, respectively.

Fig. 4. Hubble diagrams of DES-SN5YR and Pantheon+, along with the best-fit m(z) relations. The lower panels show the residuals relative to the empirical relation.

Given that the empirical relation accurately reproduces the Hubble diagrams of both samples up to z = 1.121, we assume that it is valid at least up to this redshift. Beyond this limit, we cannot reliably assess its validity because Pantheon+ contains only 19 SNe at z >  1.121. In the following, we restrict our analyses of Pantheon+, including its deep-field regions, to z < 1.121, although analyses without this cut lead to the same conclusions. Table 2 lists the parameter constraints for Pantheon+ (z < 1.121).

Table 3 lists the Δ_AICc and Δ_BIC values. Based on these values, the empirical relation, the flat ΛCDM model, and the flat Padé model provide comparable fits to the Hubble diagrams of DES-SN5YR and Pantheon+ (z < 1.121), with the flat ΛCDM model being only marginally disfavoured for DES-SN5YR (Δ_AICc = 2.2 and Δ_BIC = 2.1). In contrast, the ΛCDM model, the flat wCDM model, and the Padé cosmography are less supported by the AICc and disfavoured by the BIC, reflecting the penalty associated with their additional free parameter.

4.4. Deceleration parameter and Hubble constant

Table 4 lists the q₀ values. For the flat ΛCDM, ΛCDM, and flat wCDM models, we used q₀ = Ω_M/2 + (1 + 3w)Ω_Λ/2. Although developed solely to fit Hubble diagrams, the q₀ values of the empirical relation are consistent with those of the other m(z) relations.

We estimated H₀ using the empirical relation together with Pantheon+ (z < 1.121) and the Cepheid distance moduli (μ_Ceph) from the SH0ES sample (Riess et al. 2022). To do this, we replaced Eq. (10) with Δm_i = m_i − M − μ_Ceph, i for those SNe hosted in galaxies with μ_Ceph. We obtain M = −19.244 ± 0.030, b = −0.604 ± 0.026, and H₀ = 73.4 ± 1.0 km s⁻¹ Mpc⁻¹. This H₀ value is consistent with those obtained by Brout et al. (2022) from Pantheon+ and SH0ES, which range from 73.3 ± 1.1 to 73.6 ± 1.1 km s⁻¹ Mpc⁻¹.

4.5. Stability in the absence of low-z SNe

To evaluate the stability of the best-fit parameters of the six m(z) relations in the absence of low-z SNe, we computed ℳ and q₀ using SNe with z >  z_min, where z_min = 0.01, 0.03, 0.06, 0.10, and 0.15. Fig. 5 shows the resulting (ℳ, q₀) points. We find that the empirical relation, the flat ΛCDM model, and the flat Padé model are the most stable in the absence of low-z SNe. Because the best-fit parameters of the empirical relation are the intercept and slope of the approximately linear correlation between m − 5log(z(1 + z)) and z (see Fig. 2 and Eq. (11)), their stability depends primarily on the z baseline of the data. Therefore, they are expected to vary only mildly when low-z SNe are removed.

Fig. 5. Distribution of q0 versus ℳ for the m(z) relations fitted to DES-SN5YR and Pantheon+ (z < 1.121) with different zmin. Solid curves are the 68.27% confidence contours for zmin = 0.01, while the contour of the empirical relation is shown as a shaded region in each panel for comparison. Although the flat wCDM constraints in the (w, ΩM) plane are highly non-Gaussian, their projection onto q0 leads to nearly symmetric uncertainties, as expected since q0 = QH(z = 0), where the parameter QH(z) is almost free of the w–ΩM degeneracy (Camilleri et al. 2024).

When fitted to DES-SN5YR, the three most stable relations show a gap between cases with z_min ≥ 0.06 and those with z_min ≤ 0.03. For z_min ≥ 0.06, the selection removes at least 85% of the SNe in the low-z sample of DES-SN5YR. This likely reflects differences in the low-z sample relative to the rest of DES-SN5YR, as also discussed by Huang et al. (2025). This effect does not affect the analysis of the S1, S2, S4, and S5 deep-field regions, as they do not contain any SNe from the low-z sample.

4.6. Deep-field regions

Fig. 6 shows the Hubble diagrams for the deep-field regions. To constrain the parameters of the S1, S2, S4, and S5 regions, we selected the m and z values for the SNe in these regions, along with the corresponding entries of the covariance matrix. We then replaced Eq. (10) with Δm_i = m_i − m(z_i, v_j) for the SNe within the j-th region with parameters v_j. We applied the same procedure to the S367, N147, N235, and N6 regions. Table 5 lists the parameter constraints for the empirical relation, the flat ΛCDM model, and the flat Padé model. Fig. 6 shows the best fits.

Fig. 6. Same as Fig. 4, but for the deep-field regions.

Table 6 lists the Δ_AICc and Δ_BIC values for the empirical relation, the flat ΛCDM model, and the flat Padé model. These values indicate that the three m(z) relations provide comparable fits to the Hubble diagrams in the deep-field regions. Since the empirical relation is purely phenomenological and its parameters, uncertainties, and confidence contours can be computed analytically (Sect. 4.1), we used it as a convenient parametrisation to characterise the deep-field regions.

Fig. 7 shows the confidence contours and the marginalised distributions for ℳ and b for the deep-field regions. The b values are consistent within 1.6σ. When accounting for the 0.04 mag offset between the selection bias corrections in DES-SN5YR and Pantheon+, the ℳ values are consistent within 1.6σ, in agreement with isotropy. We obtain similar results when comparing the parameters of the flat ΛCDM and flat Padé models, for which the differences in ℳ, Ω_M, and q₀do not exceed 1.7σ, 1.5σ, and 1.7σ, respectively.

Fig. 7. Confidence contours (68.27%) and marginalised distributions of empirical relation parameters fitted to the deep-field regions.

To test whether fitting the S1, S2, S4, and S5 regions with different parameters for each region (eight parameters in total) provides a significantly better fit than an isotropic fit (two parameters for all regions), we performed a likelihood ratio test (LRT). The test statistic λ = χ_min, iso² − χ_min², where χ_min² is the value for the model with eight parameters (listed in Table 6) and χ_min, iso² is the corresponding value for the isotropic fit, follows a χ² distribution with six degrees of freedom (Wilks 1938). Fitting the empirical relation, we obtain χ_min, iso² = 1451.3, λ = 8.4, and a p-value of 0.21. For the S367, N235, N147, and N6 regions, we find χ_min, iso² = 383.5, λ = 6.7, and a p-value of 0.35. In both cases, the isotropic fit cannot be rejected. We obtain similar results when using the flat ΛCDM and flat Padé models, for which the LRT p-values lie between 0.20 and 0.36. Therefore, we find no evidence for anisotropy in the Hubble diagrams of the deep-field regions.

Table 7 lists the q₀ values for the deep-field regions. The values measured using the empirical relation are consistent within 1.6σ and are lower than zero by more than 4.8σ (99.9998% confidence). These results indicate that the Universe is undergoing accelerated expansion in a statistically consistent manner across all eight deep-field regions, under the assumption that the large-scale densities in these regions are independent of the comoving radial coordinate.

4.7. Application to other samples

4.7.1. DES-Dovekie sample

As this work was nearing completion, Popovic et al. (2025) presented a re-analysis of DES-SN5YR. The resulting recalibrated dataset, DES-Dovekie³, includes improved photometric cross-calibration and a fixed host-galaxy colour law. Since DES-Dovekie supersedes DES-SN5YR (Popovic et al. 2025), its analysis provides the best constraints and overall results from DES-SN combined with a low-z sample.

We fit the empirical relation (b = −0.603 ± 0.020), Padé cosmography ($q_0=-0.{444}_{-0.114}^{+0.072}$, $\widehat{j_0}=0.{805}_{-0.466}^{+0.881}$), flat ΛCDM (Ω_M = 0.329 ± 0.015), ΛCDM (${\mathrm{\Omega }}_M=0.{286}_{-0.062}^{+0.056}$, ${\mathrm{\Omega }}_{\mathrm{\Lambda }}=0.{600}_{-0.102}^{+0.088}$ or ${\mathrm{\Omega }}_k=0.{11}_{-0.14}^{+0.16}$), flat wCDM (${\mathrm{\Omega }}_M=0.{268}_{-0.086}^{+0.063}$, w = −0.85 ± 0.15), and the flat Padé model (q₀ = −0.470 ± 0.022). The parameters are consistent with those derived from DES-SN5YR and Pantheon+ (z < 1.121), as well as with those reported by Popovic et al. (2025) for the flat ΛCDM (Ω_M = 0.330 ± 0.015), ΛCDM (Ω_M = 0.279 ± 0.057, Ω_k = 0.14 ± 0.15), and flat wCDM models (${\mathrm{\Omega }}_M=0.{263}_{-0.078}^{+0.064}$, $w=-0.{838}_{-0.142}^{+0.130}$). The Δ_AICc and Δ_BIC values are similar to those obtained for DES-SN5YR.

The magnitude differences between a flat ΛCDM model with Ω_M from Planck CMB measurements (Ω_M^Planck = 0.3166; Planck Collaboration VI 2020) and the empirical relation have rms = 0.004 mag and a maximum absolute deviation Δm_max = 0.011 mag. In comparison, the empirical relation calibrated with DES-SN5YR and Pantheon+ (z < 1.121) yields (rms,  Δm_max) = (0.012,  0.025) for DES-SN5YR and (0.008,  0.014) for Pantheon+. Hence, the empirical relation shows close agreement with the flat ΛCDM model using Ω_M^Planck, particularly when calibrated using DES-Dovekie and Pantheon+ (z < 1.121).

For the S1, S2, S4, and S5 regions in DES-Dovekie, the b values are −0.555 ± 0.058, −0.612 ± 0.040, −0.474 ± 0.066, and −0.585 ± 0.036, respectively. These values and those from S367, N147, N235, and N6 are mutually consistent within 1.8 σ. Using an LRT, we obtain λ = 6.6 and a p-value of 0.36, indicating that the isotropic fit cannot be rejected. The corresponding q₀ values for S1, S2, S4, and S5 regions are −0.489 ± 0.054, −0.436 ± 0.037, −0.564 ± 0.061, and −0.461 ± 0.033, respectively.

We also constrained parameters using only DES SNe within DES-Dovekie. For this, we selected SNe with z >  0.093, thereby excluding the low-z sample entirely. We obtain b = −0.589 ± 0.026 for the empirical relation, Ω_M = 0.307 ± 0.020 for the flat ΛCDM model, and q₀ = −0.488 ± 0.029 for the flat Padé model. The corresponding ΔIC values are 0.5, 0.0, and 2.1, indicating that the empirical relation and the flat ΛCDM model provide comparable fits to the Hubble diagram.

4.7.2. Amalgame sample

The Amalgame sample (Popovic et al. 2024)⁴ consists of 1792 photometrically classified SNe Ia from the Sloan Digital Sky Survey (Sako et al. 2011) and Pan-STARRS1 (Scolnic et al. 2018), spanning 0.066 < z < 0.680. We obtain b = −0.613 ± 0.041 for the empirical relation, ${\mathrm{\Omega }}_m=0.{334}_{-0.027}^{+0.031}$ for the flat ΛCDM model, and $q_0=-0.{469}_{-0.041}^{+0.045}$ for the flat Padé model. These values are consistent with values obtained from Pantheon+ (z < 1.121), DES-SN5YR, and DES-Dovekie. Fig. 8 shows the Hubble diagram for Amalgame. The ΔIC values for the empirical relation, the flat ΛCDM model, and the flat Padé model are 0.4, 0.0, and 0.6, respectively, showing that all three provide comparable fits to the Hubble diagram.

Fig. 8. Same as Fig. 4, but for the Amalgame SN sample.

4.7.3. Baryon acoustic oscillation distances and cosmic chronometers

Although the empirical relation is designed to fit Hubble diagrams of standard candles, we applied it to available d_M(z) and H(z) datasets assuming the FLRW metric, which allowed us to transform the empirical relation into expressions for d_M(z) (Eq. 13) and H(z) (Eq. 14).

Abdul Karim et al. (2025) reported 13 distance measurements from baryon acoustic oscillations (BAO) in the second data release (DR2) of the Dark Energy Spectroscopic Instrument (DESI), relative to the sound horizon (r_d), at effective redshifts z_eff ranging from 0.295 to 2.33⁵. These include six d_M(z)/r_d values, six d_H(z)/r_d values (with d_H(z) =c/H(z)), and one d_V(z)/r_d value, where d_V(z)≡[zd_M(z)²d_H(z)]^1/3. Fitting Eqs. (13) and (14) to the seven BAO distances with z_eff < 1.121, and assuming Ω_k = 0, we obtain b = −0.584 ± 0.016. For comparison, the b values obtained from Pantheon+ (z < 1.121), DES-SN5YR, DES-Dovekie, and Amalgame, when all SN samples are restricted to z ≥ 0.295, are −0.697 ± 0.050, −0.622 ± 0.035, −0.618 ± 0.033, and −0.688 ± 0.070, respectively. Compared to the BAO constraint, the latter three values are consistent within 1.4σ, while the Pantheon+ (z < 1.121) result is marginally consistent, showing a 2.2σ difference.

Alfano et al. (2026) collected 34 H(z) values measured with cosmic chronometers, spanning 0.07 ≤ z ≤ 1.965. Fitting Eq. (14) to the 27 H(z) values with z < 1.121, and assuming Ω_k = 0, we obtain b = −0.635 ± 0.095. This value is consistent with those from Pantheon+ (z < 1.121), DES-SN5YR, DES-Dovekie, and Amalgame when all SN samples are restricted to z ≥ 0.07, yielding b = −0.617 ± 0.033, −0.607 ± 0.027, −0.591 ± 0.024, and −0.613 ± 0.041, respectively.

5. Discussion

5.1. Jerk parameter

The value of j₀ provides an observational test for the flat ΛCDM model (Bochner et al. 2015). Indeed, for the non-flat wCDM model (Eq. 4), Eq. (6) yields j₀ = 1 − Ω_k + 3(1 + w)(2q₀ + Ω_k − 1)/2, which reduces to j₀ = 1 for the flat ΛCDM model. Therefore, if j₀ deviates from unity, it could indicate that the Universe is not flat, that dark energy is not a cosmological constant, or that Eq. (4) is incorrect. Given that $\widehat{j_0}=1$ for the flat ΛCDM model, the test remains valid if $\widehat{j_0}$ is used instead of j₀.

Table 8 lists the $\widehat{j_0}$ values. For the empirical relation, we use

$$\begin{array}{c}\hfill \widehat{j_0}=\frac{9}{4}{(\frac{5}{9}+q_0)}^2+\frac{5}{9},\end{array}$$(16)

which is derived from Eqs. (6) and (14). For DES-SN5YR and Pantheon+ (z < 1.121), the $\widehat{j_0}$ values are consistent within 1.2σ. In particular, the $\widehat{j_0}$ values for the ΛCDM model, the flat wCDM model, and the Padé cosmography are consistent with unity within 1.2σ. In contrast, the $\widehat{j_0}$ values for the empirical relation are approximately 0.4 below unity, with an overwhelming significance of at least 30σ. Assuming w = −1, these values translate to ${\mathrm{\Omega }}_k=0.{208}_{-0.007}^{+0.005}$ and $0.{202}_{-0.006}^{+0.005}$ for Pantheon+ (z < 1.121) and DES-SN5YR, respectively; or, assuming Ω_k = 0, to w = −0.853 ± 0.001 and −0.854 ± 0.001, respectively. However, the empirical relation is intended to fit the Hubble diagram, not its derivatives. Therefore, although the empirical relation provides q₀ values consistent with those from other m(z) relations (see Table 4), it may misestimate $\widehat{j_0}$.

5.2. Cosmic anisotropies through hemisphere comparison

Pantheon+ has been widely used to search for anisotropies in cosmological parameters using the hemisphere comparison (HC) method (Schwarz & Weinhorst 2007) under different cosmological models (e.g. Mc Conville & Colgáin 2023; Clocchiatti et al. 2024; Hu et al. 2024a,b). The HC method splits the sky into two opposite hemispheres (‘up’ and ‘down’), constrains the cosmological parameter x in each, computes a statistic such as the signal-to-noise ratio (S/N) of Δx = |x_up − x_down|, and repeats the procedure for different directions across the sky.

We applied the HC method to Pantheon+ (z < 1.121) using the empirical relation. We selected a set of directions on the celestial sphere using HEALPY (Górski et al. 2005; Zonca et al. 2019) with a grid resolution parameter N_side = 128. This provided 196.608 equal-area pixels on the sky, corresponding to a set of 98.304 independent directions. For each direction, we computed S/N(Δx) =Δx/(σ_up² + σ_down² − 2σ_up, down)^1/2, where σ_up and σ_down are the uncertainties of x_up and x_down, respectively, and σ_up, down is their covariance. Fig. 9 shows the sky maps of S/N(Δℳ) and S/N(Δb). Their maximum values are 3.7 at (RA, Dec) =(239.2° ,  − 70.17° ) and 3.3 at (RA, Dec) =(212.3° ,  − 39.07° ), respectively.

Fig. 9. Sky maps of S/N(Δℳ) and S/N(Δb). Dots represent the SNe in Pantheon+ (z < 1.121) and circles mark the maximum S/N values.

To evaluate the probabilities of obtaining S/N(Δℳ)≥3.7 and S/N(Δb)≥3.3 by chance when scanning many directions on the sky (the look-elsewhere effect; Bayer & Seljak 2020), we performed 10⁴ simulations. In each simulation, we randomised the RA and Dec coordinates of the SNe and produced the S/N(Δℳ) and S/N(Δb) maps using the HC method and the same directions as above. The probabilities of obtaining S/N(Δℳ)≥3.7 and S/N(Δb)≥3.3 by chance are 6.0% and 11.2%, respectively, corresponding to 1.9σ and 1.6σ significances. Therefore, we find no evidence for anisotropies in ℳ and b.

If the directions of maximum S/N(Δℳ) and S/N(Δb) were real anisotropies, detecting them at a 3σ level with the HC method would require adding 2500 and 4400 SNe to Pantheon+ (z < 1.121), respectively. An alternative approach would be to perform deep-field SN surveys towards the directions of the candidate anisotropies. This is particularly important given the absence of high-z SNe near those directions (see Fig. 9).

5.3. JWST, LSST, CSST, and Roman

Given the small number of SNe at z >  1.121, the validity of the empirical relation beyond this redshift remains to be verified. At present, the James Webb Space Telescope (JWST) is observing SNe Ia at z >  1.7 (Casey et al. 2023; DeCoursey et al. 2025). Siebert et al. (2025) report $\mu =46.{08}_{-0.18}^{+0.19}$ for SN 2025ogs at z = 2.05, while Pierel et al. (2025) report $\mu =46.{10}_{-0.18}^{+0.17}$ for SN 2023aeax at z = 2.15 and $\mu =47.{14}_{-0.24}^{+0.21}$ for SN 2023adsy at z = 2.903. At these redshifts, adopting the parameters inferred from Pantheon+ (Table 2) and H₀ = 70 km s⁻¹ Mpc⁻¹, the empirical relation predicts μ = 45.92, 46.04, and 46.71, whereas the ΛCDM model predicts μ = 46.0, 46.13, and 46.94, respectively. The μ values for SNe 2025ogs and 2023aeax are consistent with both relations, whereas that for SN 2023adsy is more consistent with the ΛCDM model than with the empirical relation. However, due to its red colour, it is uncertain whether SN 2023adsy is representative of high-z SNe Ia or an outlier (Pierel et al. 2024, 2025).

Beyond JWST, future deep-field surveys with the China Space Station Telescope (CSST; Li et al. 2023) and the Nancy Grace Roman Space Telescope (Hounsell et al. 2023) will observe SNe Ia up to z = 1.3 and z = 1.7, respectively. These SNe, combined with those observed by JWST and the 19 SNe at 1.121 < z ≤ 2.261 in Pantheon+, will be essential for determining the maximum z at which the empirical relation remains valid.

The deep rolling surveys of the Vera C. Rubin Observatory Legacy Survey of Space and Time (LSST; Gris et al. 2024) will provide photometry of thousands of SNe Ia with z ≲ 1.1. Thus, the Hubble diagram for each of the five LSST Deep Drilling Fields can be analysed with the empirical relation. Of these fields, the newly defined Euclid Deep Field South (RA = 61.241°, Dec = −48.423°), located at an angular separation of 21° from the S5 field, will enable study of cosmic acceleration in a direction distinct from those analysed here.

6. Conclusions

We present a new empirical m(z) relation based on data from Pantheon+ and DES-SN5YR. This relation fits their Hubble diagrams with only two parameters, ℳ and b, both of which can be computed analytically. Notably, under the FLRW metric, b ∝ q₀ + 1.

For the DES-SN5YR and Pantheon+ Hubble diagrams, the empirical relation provides fits up to z = 1.121 that closely match those obtained with the Padé cosmography and the ΛCDM, flat ΛCDM, flat wCDM, and flat Padé models. The validity of the empirical relation beyond z = 1.121 cannot be reliably assessed given the small number of SNe at z >  1.121. Based on the Δ_IC values, the empirical relation performs better than the ΛCDM model, the flat wCDM model, and the Padé cosmography, and comparably to the flat ΛCDM and flat Padé models.

For DES-SN5YR and Pantheon+ (z < 1.121), we obtain b = −0.628 ± 0.021 and b = −0.606 ± 0.026, respectively. We also applied the empirical relation to DES-Dovekie and Amalgame, obtaining b = −0.603 ± 0.020 and b = −0.613 ± 0.041, respectively. For Pantheon+ (z < 1.121) combined with SH0ES, we obtain H₀ = 73.4 ± 1.0 km s⁻¹ Mpc⁻¹.

The empirical relation, as well as the flat ΛCDM and flat Padé models, remains stable in the absence of low-z SNe. Fitting the empirical relation to the Hubble diagrams of the eight deep-field regions, we find that their fitted ℳ and b parameters are consistent within 1.6σ, in agreement with isotropy. Applying an LRT, we find no evidence of anisotropy. This result is further supported by using the empirical relation in the HC method applied to Pantheon+ (z < 1.121), which finds no statistically significant evidence for anisotropies in ℳ and b.

Furthermore, the b values from the eight deep-field regions translate to q₀ values ranging from −0.6 to −0.4, which are consistent within 1.6σ and lower than zero at 4.8σ. These results strongly support an accelerating universe, under the assumption that the large-scale density in each deep-field region is independent of the comoving radial coordinate.

Assuming an FLRW metric, the empirical relation can also be applied to BAO distances from DESI DR2, and to H(z) values from cosmic chronometers. The derived b values, assuming Ω_k = 0, are largely consistent with those from DES-SN5YR, Pantheon+ (z < 1.121), DES-Dovekie, and Amalgame.

The empirical relation, m(z) =ℳ + bz + 5log(z(1 + z)), or equivalently $d_L(z)=\frac{\mathit{cz}}{H_0}(1+z){10}^{bz/5}$, provides a straightforward quantitative tool for fitting Hubble diagrams of SNe Ia without requiring a low-z sample. It enables tests of basic properties of the Universe in a simplified manner, reducing the computational cost typically associated with physically motivated cosmological models. Our analysis demonstrates that the empirical relation is valid at least up to z ≈ 1.1. Ongoing and upcoming deep SN surveys will allow us to test its validity at higher redshifts.

Acknowledgments

We acknowledge the Director of the Millennium Institute of Astrophysics (MAS), Francisco Förster, for financial assistance through ANID grant ICN12_009, and acknowledge complementary funding from ANID/FONDECYT grant 1251692. OR is supported by the Rubin-Chile Fund under grant DIA2650. We thank Felipe Asenjo and Ariel Órdenes for helpful discussions.

References

All Tables

All Figures

	Fig. 1. Sky distribution of Pantheon+ and DES-SN5YR SNe. Dashed lines indicate the Galactic equator. Plus and cross symbols mark the north and south Galactic poles, respectively. Circular and triangular regions are also shown.
In the text

	Fig. 2. Relation between m − 5log(z(1 + z)) and z for Pantheon+ and DES-SN5YR. Binned data (blue squares) with error bars are shown for visualisation purposes only, both here and throughout the paper.
In the text

	Fig. 3. Confidence contours at the 68.27% and 95.45% levels of the empirical relation parameters, obtained analytically and using EMCEE.
In the text

	Fig. 4. Hubble diagrams of DES-SN5YR and Pantheon+, along with the best-fit m(z) relations. The lower panels show the residuals relative to the empirical relation.
In the text

Fig. 5. Distribution of q0 versus ℳ for the m(z) relations fitted to DES-SN5YR and Pantheon+ (z < 1.121) with different zmin. Solid curves are the 68.27% confidence contours for zmin = 0.01, while the contour of the empirical relation is shown as a shaded region in each panel for comparison. Although the flat wCDM constraints in the (w, ΩM) plane are highly non-Gaussian, their projection onto q0 leads to nearly symmetric uncertainties, as expected since q0 = QH(z = 0), where the parameter QH(z) is almost free of the w–ΩM degeneracy (Camilleri et al. 2024).

In the text

	Fig. 6. Same as Fig. 4, but for the deep-field regions.
In the text

	Fig. 7. Confidence contours (68.27%) and marginalised distributions of empirical relation parameters fitted to the deep-field regions.
In the text

	Fig. 8. Same as Fig. 4, but for the Amalgame SN sample.
In the text

	Fig. 9. Sky maps of S/N(Δℳ) and S/N(Δb). Dots represent the SNe in Pantheon+ (z < 1.121) and circles mark the maximum S/N values.
In the text