© The Authors 2026

1. Introduction

The standard cosmological model, Λ cold dark matter (ΛCDM), has achieved remarkable success in explaining a wide range of observations, from the cosmic microwave background (CMB) anisotropies (Aghanim et al. 2020) to the large-scale structure of the Universe (Peebles 1993). However, with the advent of precision cosmology, increasing tensions – most notably the Hubble tension (H₀) – have emerged between early-Universe predictions and late-Universe measurements (Riess et al. 2022; Di Valentino et al. 2021; Abdalla et al. 2022). Before attributing these discrepancies to exotic dark energy or modified gravity, it is imperative to scrutinize the fundamental observational principles that underpin our distance measurements (Perivolaropoulos & Skara 2022). Among these, the cosmic distance-duality relation (CDDR) stands as a cornerstone of observational geometry (Etherington 1933).

First derived by Etherington in 1933, the CDDR, also known as the reciprocity theorem, relates the luminosity distance D_L and the angular diameter distance D_A of a source at redshift z via the simple identity,

$$\begin{array}{c}\hfill \frac{D_L}{D_A}{(1+z)}^{-2}=1.\end{array}$$(1)

The validity of this relation rests on three crucial theoretical pillars: (i) the space-time is described by a Riemannian metric theory of gravity (Ellis 1971, 2007; ii) photons propagate along unique null geodesics (Uzan et al. 2004); and (iii) the photon number is conserved during propagation (Bassett & Kunz 2004). Consequently, any observational violation of the CDDR would signal fundamental revisions to physics (Barua et al. 2026). Possible mechanisms include nonmetric theories of gravity (Hees et al. 2014; Vagnozzi et al. 2023), the existence of gray dust in the intergalactic medium (Corasaniti 2006; Holanda et al. 2013), or the coupling of photons with axion-like particles into which they might oscillate (Tiwari 2017; Keil et al. 2026).

To test this relation, a phenomenological parameter η is commonly introduced to quantify potential deviations, parameterized as D_L = D_A(1 + z)^2 + η (Holanda et al. 2010; Li et al. 2011). Standard cosmology predicts η = 0 (Ellis 2007). Observational constraints on η require dual measurements of D_L and D_A at the same redshift (Holanda et al. 2012). While type Ia supernovae (SNe Ia) are the premier standard candles for D_L (Scolnic et al. 2022), obtaining model-independent D_A measurements remains challenging. Galaxy clusters (GCs) provide a powerful solution (De Filippis et al. 2005). By combining X-ray surface brightness observations with the Sunyaev-Zel’dovich effect (SZE) – the inverse Compton scattering of CMB photons by hot intracluster gas (Sunyaev & Zeldovich 1972) – the angular diameter distance to clusters can be determined directly, independent of the cosmic expansion history (Cavaliere & Fusco-Femiano 1976; Bonamente et al. 2006).

Historically, tests of the CDDR have often relied on assuming a specific cosmological model (typically flat ΛCDM) to align datasets or to reconstruct distance functions (Uzan et al. 2004; Holanda et al. 2010). Within these model-dependent frameworks, deviations from the CDDR are frequently interpreted as signatures of cosmic opacity, a topic extensively explored in various cosmological contexts (e.g., Hu et al. 2017). Furthermore, such tests have been extended to alternative cosmological scenarios, including the R_h = ct universe (Hu & Wang 2018). However, these model-dependent approaches inevitably introduce circularity and bias if the assumed background cosmology deviates from reality (Liang et al. 2013). To circumvent this, model-independent reconstruction techniques, such as Gaussian processes (Seikel et al. 2012; Zhang 2014; Keil et al. 2026), polynomial fits (Montiel et al. 2021), and Padé approximants (Barua et al. 2026), have been widely adopted. Nevertheless, these methods can suffer from overfitting or boundary artifacts, particularly when data are sparse at high redshifts (Kanodia et al. 2026).

A more transparent approach is the “matched pair” technique pioneered by Holanda et al. (2010, 2012), where D_L and D_A measurements are strictly paired based on redshift proximity, thereby eliminating the dependence on the fiducial model entirely. To further improve the robustness of this technique, Zhou et al. (2021) introduced the dimensionless distance error consistency method, which allows for an adaptive matching tolerance based on data quality. This rigorous approach is arguably superior to simple redshift cuts and has been successfully applied in recent work (Hu 2023; Hu et al. 2023). However, a common limitation in these previous analyses – including those utilizing the advanced matching criterion – is the neglect of the full systematic covariance matrix of SNe Ia during the analysis. In this work, we combined this advanced matching strategy with a rigorous treatment of the full Pantheon+ covariance matrix to ensure a stringent test.

A critical, yet often overlooked systematic in CDDR tests is the intrinsic evolution of the probes themselves (Tutusaus et al. 2017). While SNe Ia are standardized to a high degree, recent analyses suggest that their absolute magnitude (M_B) may evolve with redshift due to progenitor age or metallicity drift (Kang et al. 2020; Nicolas et al. 2021; Scolnic et al. 2022). Neglecting this evolution can mimic a violation of the CDDR, leading to spurious detections of η ≠ 0 or biased constraints on opacity parameters (Tutusaus et al. 2017; Kanodia et al. 2026). Therefore, a rigorous test must break the degeneracy between the fundamental space-time geometry (the CDDR parameter) and the astrophysical systematics (the source evolution) (Alfano & Luongo 2026).

In this work, we present a stringent, model-independent test of the CDDR using a carefully selected sample of 38 GCs, providing D_A measurements via the X-ray/Sunyaev-Zel’dovich (SZ) technique (Bonamente et al. 2006). We paired these clusters with the latest SNe Ia from the Pantheon+ compilation (Scolnic et al. 2022), which offers significantly reduced statistical and systematic uncertainties compared to previous samples. Unlike previous studies that often fix the supernova absolute magnitude, we adopted a joint analysis framework using Markov chain Monte Carlo (MCMC) methods. We simultaneously constrained the CDDR parameter η and the redshift-dependent evolution of the SNe Ia absolute magnitude, parameterized as M_B(z) = M₀ + ε ⋅ z. By utilizing the covariance information of the Pantheon+ sample through the extracted sub-covariance matrix for the matched SNe Ia, this study aims to provide robust, degeneracy-free constraints on the validity of the Etherington reciprocity theorem.

In addition to the baseline cluster–SNe Ia matched-pair analysis, we further incorporated recent DESI baryon acoustic oscillation (BAO) measurements in the form of D_M/r_d (Adame et al. 2025a,b). Since D_A = D_M/(1 + z), these BAO measurements constrain the CDDR parameter η through the same generalized reciprocity relation as the cluster D_A measurements and may equivalently be regarded as four additional angular-diameter-distance data points once r_d is specified. This extension therefore preserves the model-independent nature of the test while providing additional geometric leverage for the joint analysis.

The structure of this paper is outlined as follows. Section 2 details the observational datasets utilized in this study, along with the specific methodology employed for data pairing. In Sect. 3, we present the statistical framework and report the numerical constraints obtained from our analysis. Finally, Sect. 4 summarizes our findings and discusses their implications.

2. Data and methodology

In this section, we provide a detailed description of the observational datasets employed in our analysis, including the GC sample derived from X-ray and SZE measurements and our type Ia supernova (SNe Ia) compilation. We also introduce the DESI 2024 BAO measurements used as an extension of the baseline dataset. Furthermore, we outline the matching procedure used to construct pairs of angular diameter and luminosity distances at identical redshifts, which forms the basis of our model-independent test.

2.1. Galaxy cluster sample

For the angular diameter distance (D_A) measurements, we utilized the sample of 38 GCs compiled and analyzed by Bonamente et al. (2006). This sample spans a redshift range of 0.14 ≤ z ≤ 0.89 and represents one of the most robust datasets for SZE-based distance determinations.

The distances to these clusters were derived by combining X-ray surface brightness observations with SZE measurements. The X-ray data were obtained from the Chandra X-ray Observatory, which provides high-resolution imaging and spectroscopy essential for characterizing the intracluster medium (ICM). The SZE data, which measure the spectral distortion of the CMB caused by inverse Compton scattering off hot ICM electrons, were collected using the interferometric arrays at the Owens Valley Radio Observatory (OVRO) and the Berkeley-Illinois-Maryland Association (BIMA).

The fundamental principle of this distance measurement technique relies on the different dependencies of the SZE intensity and X-ray surface brightness on the electron density (n_e). The SZE temperature decrement is proportional to the line-of-sight integral of the electron pressure (ΔT_SZE ∝ ∫n_eT_e dl), whereas the X-ray surface brightness is proportional to the line-of-sight integral of the electron density squared (S_X ∝ ∫n_e²Λ_ee dl). By assuming a geometric model for the cluster (e.g., spherical symmetry) and modeling the radial profiles of the gas density and temperature, one can solve for the angular diameter distance D_A directly, independent of the cosmic distance ladder or the expansion history of the Universe.

In the analysis by Bonamente et al. (2006), a hydrostatic equilibrium model was employed to account for radial variations in plasma density, temperature, and metal abundance, thereby reducing systematic uncertainties associated with simpler isothermal β-models. The final sample consists of 38 clusters with high-quality joint X-ray and SZE detections, providing a reliable dataset for testing the CDDR. It is worth noting that the derivation of D_A from X-ray and SZE observations relies on standard local physics, including photon propagation and the thermodynamic properties of the ICM, which are typically formulated within a metric theory of gravity such as general relativity. However, this procedure does not rely on any assumptions about the global expansion history, and therefore remains independent of any specific background cosmological model.

2.2. Type Ia supernova sample

For the luminosity distance (D_L) measurements, we employed the latest Pantheon+ compilation, as presented by Brout et al. (2022) and Scolnic et al. (2022). Pantheon+ represents a significant expansion and improvement over the original Pantheon sample, comprising 1701 light curves of 1550 distinct spectroscopically confirmed SNe Ia. The sample covers a broad redshift range of 0.001 < z < 2.26, integrating data from multiple surveys including CfA1-4, CSP, PS1, SDSS, SNLS, and high-redshift samples from HST.

A key feature of Pantheon+ is its substantial increase in the number of low-redshift supernovae and its comprehensive treatment of systematic uncertainties. The compilation provides a detailed covariance matrix (C_{stat + sys}) that accounts for uncertainties arising from calibration, photometric zero points, intrinsic scatter, selection bias, and peculiar velocities. This robust quantification of errors is crucial for our joint analysis, as it allows us to propagate systematic correlations through the relevant sub-covariance matrix extracted from the full Pantheon+ covariance when constraining the CDDR parameter.

In the context of our matched-pair analysis, the Pantheon+ dataset is particularly advantageous because it includes the SH0ES (Supernovae and H₀ for the equation of state of dark energy) subsample SNe Ia in galaxies that also host Cepheid variables. While the SH0ES calibration is typically used to fix the absolute magnitude M_B for H₀ measurements, in this work, we treated the absolute magnitude and its potential redshift evolution (ε) as free nuisance parameters. We utilized the corrected apparent magnitudes (m_B^corr) provided by the Pantheon+ analysis, which were standardized for light-curve shape (x₁), color (c), and host-galaxy mass step corrections.

2.3. Data pairing

To perform a model-independent test of the CDDR, we paired each GC (providing D_A) with a type Ia supernova (providing D_L) at (nearly) the same redshift. A key challenge in such analyses is the redshift mismatch between independent datasets, which can introduce nonnegligible systematic biases if not properly controlled.

Rather than applying a simple redshift threshold (e.g., Δz < 0.005), we adopted a distance-based matching criterion, which more directly reflects the impact of redshift differences on cosmological distance measurements. This approach ensures that the pairing is performed in a physically meaningful way, minimizing the propagation of mismatch errors into the final constraints.

To implement this, we assumed a fiducial flat ΛCDM cosmology with Ω_m = 0.3 solely for the purpose of defining the matching tolerance. This assumption did not affect the model-independence of the final CDDR test, as it was only used to determine the relative distance differences between objects, rather than derive the distances themselves. We computed the line-of-sight comoving distance, D_C(z), for both the GCs and the supernovae, and required that a valid pair satisfies

$$\begin{array}{c}\hfill \frac{|D_C(z_{\mathrm{cluster}})-D_C(z_{\mathrm{SNe}})|}{D_C(z_{\mathrm{cluster}})}\le 5\%,\end{array}$$(2)

which defines our baseline matching tolerance.

This dimensionless distance-based criterion has been shown to be more robust than simple redshift cuts, as it naturally accounts for the redshift dependence of cosmological distances. It has been successfully applied in previous studies, for example in constraining the mass density profiles of strong gravitational lensing systems (Hu 2023) and in model-independent calibrations of SNe Ia absolute magnitudes (Hu et al. 2023).

To assess the robustness of our results against the choice of matching criterion, we additionally considered a stricter selection with a tolerance of ΔD_C/D_C ≤ 3% and compared the resulting constraints with those obtained from the baseline sample. From the pool of potential matches satisfying Eq. (2), we constructed the final sample using a greedy selection algorithm. All candidate pairs were ranked according to their absolute distance difference |D_C(z_cluster)−D_C(z_SNe)|, and the pair with the smallest difference was selected iteratively. Once a pair was selected, both the corresponding GC and supernova were removed from the pool to ensure that each object was used at most once (sampling without replacement). This procedure minimizes the overall mismatch across the dataset and yields a well-defined set of independent pairs.

Using this method, we obtain a final matched sample of 38 GC–SNe Ia pairs, as illustrated in Fig. 1. In parallel with the pair selection, we extracted the corresponding 38 × 38 sub-covariance matrix from the full Pantheon+ covariance matrix (C_{stat + sys}). This step preserves both statistical and systematic correlations among the selected SNe Ia, which is essential for a consistent likelihood analysis.

Fig. 1. Redshift mismatch (Δz = zcluster − zSN) between the paired GCs and SNe Ia used in the baseline matched sample.

The final matched sample consists of 38 GC–SNe Ia pairs. The detailed properties of these pairs, including cluster redshifts, angular diameter distances, and the matched SNe Ia corrected magnitudes, are listed in Table A.1 in Appendix A.

2.4. DESI BAO supplement

To address the request for additional angular-diameter-distance information, we supplemented the baseline cluster-only analysis with four DESI 2024 BAO measurements of D_M/r_d at z = 0.510, 0.706, 0.930, and 1.317 (Adame et al. 2025a,b). These values were selected from the DESI Gaussian BAO summary products, and the corresponding 4 × 4 sub-covariance matrix was extracted from the full DESI covariance matrix. In our implementation this sub-covariance is diagonal for the selected set of D_M/r_d points, but we nevertheless retained the formal matrix treatment for completeness. The corresponding BAO measurements are listed in Table 1.

The BAO extension was used as a supplement to the matched-pair analysis rather than as a replacement for it. We considered two versions of this extension: one in which the sound horizon r_d is treated as a free parameter, and one in which a Gaussian Planck prior is imposed on r_d. This allowed us to separate the constraining power of the DESI points themselves from that of an external early-Universe calibration.

We emphasize that the four DESI BAO measurements of D_M/r_d constrain the CDDR parameter η through the same generalized reciprocity relation D_L = D_A(1 + z)²(1 + ηz) as the cluster D_A sample. Since D_A = D_M/(1 + z), the BAO points may equivalently be regarded as four additional angular-diameter-distance measurements, once r_d is specified by a prior or treated as a free parameter, and added directly to the matched-pair sample. We explicitly verified this equivalence (see Sect. 3.1 below): a fully symmetric formulation in which the four DESI points are converted to D_A and concatenated with the 38 cluster D_A measurements in a unified matched-pair likelihood yields parameter constraints statistically indistinguishable from those obtained with our D_M/r_d-space likelihood. The two formulations differ only by a Jacobian transformation. The DESI BAO points therefore play the same physical role as the clusters in constraining η; it is only their limited number (four points vs. 38 clusters) that makes their statistical leverage on η alone relatively modest. The main numerical effect of including DESI is to tighten constraints in the (M₀, ε, r_d) subspace, which then indirectly sharpens the marginalized posterior on η through the parameter correlations.

3. Analysis and results

3.1. Statistical method

We adopted a joint likelihood analysis to simultaneously constrain the CDDR parameter η, the SNe Ia absolute magnitude M_B, and its evolutionary parameter ε. The distance duality relation was parameterized as

$$\begin{array}{c}\hfill D_L(z)=D_A(z){(1+z)}^2(1+\eta z).\end{array}$$(3)

Ideally, if the reciprocity theorem holds, η = 0.

For the SNe Ia, the observed distance modulus including redshift evolution is written as

$$\begin{array}{c}\hfill {\mu }_{\mathrm{SN}}^{\mathrm{obs}}(z_i)=m_B^{\mathrm{corr}}-(M_0+\epsilon z_i),\end{array}$$(4)

where M₀ is the absolute magnitude at z = 0, and ε quantifies a possible linear evolution. For the GCs, the observed angular diameter distance D_A was converted to a distance modulus prediction using the modified CDDR parameter η:

$$\begin{array}{c}\hfill {\mu }_{\mathrm{cluster}}^{\mathrm{th}}(z_i,\eta )=5{log}_{10}\left[D_A^{\mathrm{obs}}(z_i){(1+z_i)}^2(1+\eta z_i)\right]+25.\end{array}$$(5)

The GC angular diameter distance measurements possess asymmetric statistical uncertainties (σ_DA, + and σ_DA, −). To rigorously incorporate this non-Gaussian feature, we utilized a split-normal (or piecewise Gaussian) likelihood formulation. First, we propagated these asymmetric errors from the linear distance space (D_A) to the distance modulus space (μ) using error propagation: σ_μ, ± ≈ (5/ln10)(σ_DA, ±/D_A).

In our joint analysis, the total covariance matrix C_tot combines the supernova sub-covariance matrix (extracted from the full Pantheon+ covariance) and the cluster diagonal variances. The asymmetric uncertainties in clusters arise from the nonlinear propagation of systematic errors in X-ray temperature and SZE measurements (Bonamente et al. 2006). To strictly preserve this information, the variance term for the i-th cluster, ${\sigma }_{i,\mathrm{cluster}}^2$, was selected dynamically based on the residual direction:

$$\begin{array}{c}\hfill {\sigma }_{i,\mathrm{cluster}}=\{\begin{array}{cc}{\sigma }_{\mu ,i,+}\hfill & \text{if}{\mu }_{\mathrm{SN}}^{\mathrm{obs}}(z_i)>{\mu }_{\mathrm{cluster}}^{\mathrm{th}}(z_i).\hfill \\ {\sigma }_{\mu ,i,-}\hfill & \text{otherwise}.\hfill \end{array}\end{array}$$(6)

This formulation ensures that the likelihood properly accounts for the asymmetric error distribution of the cluster measurements.

We define the residuals vector as

$$\begin{array}{c}\hfill \mathbf{\Delta }={\mathit{\mu }}_{\mathrm{SN}}^{\mathrm{obs}}-{\mathit{\mu }}_{\mathrm{cluster}}^{\mathrm{th}}.\end{array}$$(7)

We assumed that the observational uncertainties of the GCs and the SNe Ia are independent, since they arise from different instruments and physical processes (X-ray and SZ observations versus optical light curves). Therefore, no cross-covariance term between cluster and SN measurements is included.

The final log-likelihood function is constructed as

$$\begin{array}{c}\hfill ln\mathcal{L}=-\frac{1}{2}[{\mathbf{\Delta }}^T{\mathbf{C}}_{\mathrm{tot}}^{-1}\mathbf{\Delta }+ln(det{\mathbf{C}}_{\mathrm{tot}})+Nln(2\pi )].\end{array}$$(8)

We performed the parameter estimation using the MCMC method with the emcee Python package. We assumed flat priors for all parameters: M₀ ∈ [ − 20.5, −18.0], ε ∈ [ − 1.0, 1.0], and η ∈ [ − 2.0, 2.0].

For the DESI-extended analysis, the parameter set becomes {M₀, ε, η, r_d}. We verified that the 38 SNe Ia matched to the GCs are distinct from the four SNe Ia used for the BAO anchor, ensuring no object-level overlap.

Although both subsets originate from the Pantheon+ compilation, a residual cross-covariance may arise from shared calibration systematics. However, given the small number of BAO-anchoring SNe Ia and the relatively large uncertainties of the cluster measurements, this contribution is subdominant and does not affect the parameter inference at the current level of precision. Therefore, we neglected this cross-covariance and treated χ_cl² and χ_BAO² as additive.

In our implementation, the theoretical prediction for the BAO observable was constructed from the SN-inferred luminosity distance combined with the deformed CDDR relation. The SNe Ia distance modulus at the DESI redshifts is written as

$$\begin{array}{c}\hfill {\mu }_{\mathrm{SN}}(z)=m_B-(M_0+\epsilon z),\end{array}$$(9)

from which the luminosity distance is obtained as

$$\begin{array}{c}\hfill D_L(z)={10}^{({\mu }_{\mathrm{SN}}-25)/5}.\end{array}$$(10)

Using the generalized CDDR relation, the transverse comoving distance is

$$\begin{array}{c}\hfill D_M(z)=\frac{D_L(z)}{(1+z)(1+\eta z)}·\end{array}$$(11)

The theoretical BAO observable is therefore

$$\begin{array}{c}\hfill {\left(\frac{D_M}{r_d}\right)}_{\mathrm{th}}=\frac{1}{r_d}\frac{D_L(z)}{(1+z)(1+\eta z)}·\end{array}$$(12)

The residual vector is defined as

$$\begin{array}{c}\hfill {\mathrm{\Delta }}_{\mathrm{BAO}}={\left(\frac{D_M}{r_d}\right)}_{\mathrm{obs}}-{\left(\frac{D_M}{r_d}\right)}_{\mathrm{th}}.\end{array}$$(13)

The total covariance matrix is constructed as

$$\begin{array}{c}\hfill C_{\mathrm{BAO}}=C_{\mathrm{DESI}}+J{C}_{\mathrm{SN}}{J}^T,\end{array}$$(14)

where C_DESI = diag(σ_i²) and the Jacobian matrix J propagates the SN covariance into the BAO observable space,

$$\begin{array}{c}\hfill J_{\mathit{ii}}=\frac{ln10}{5}{\left(\frac{D_M}{r_d}\right)}_i.\end{array}$$(15)

The BAO log-likelihood is written as

$$\begin{array}{c}\hfill ln{\mathcal{L}}_{\mathrm{BAO}}=-\frac{1}{2}[{\mathrm{\Delta }}_{\mathrm{BAO}}^T{C}_{\mathrm{BAO}}^{-1}{\mathrm{\Delta }}_{\mathrm{BAO}}+lndet{C}_{\mathrm{BAO}}+Nln(2\pi )].\end{array}$$(16)

We considered both the free-r_d case and the case with a Planck prior:

$$\begin{array}{c}\hfill {\chi }_{\mathrm{prior}}^2={\left(\frac{r_d-147.09}{0.26}\right)}^2.\end{array}$$(17)

We note that the BAO log-likelihood above, which is written in D_M/r_d space, is equivalent – up to a Jacobian transformation – to a fully symmetric formulation in which the BAO points are first converted into angular-diameter distances via

$$\begin{array}{c}\hfill D_A^{\mathrm{BAO}}(z)={\left(\frac{D_M}{r_d}\right)}_{\mathrm{obs}}\frac{r_d}{1+z},\end{array}$$(18)

and then concatenated with the cluster D_A measurements in a single, unified matched-pair likelihood of the form of Eq. (8). We verified that both formulations are mathematically equivalent and yield statistically indistinguishable posteriors on the model parameters. Therefore, our choice of the D_M/r_d-space likelihood is purely a matter of convenience, as the DESI products are natively released in D_M/r_d space with their covariance defined there.

3.2. Baseline matched-pair constraints

The constraints for the baseline 5% matched sample, at the 68% confidence level (CL), are

$$\begin{array}{cc}\hfill M_0& =-19.{460}_{-0.124}^{+0.126},\hfill \\ \hfill \epsilon & =-0.{184}_{-0.574}^{+0.724},\hfill \\ \hfill \eta & =0.{050}_{-0.307}^{+0.348}.\hfill \end{array}$$(19)

The best-fit log-likelihood is lnℒ_best = −44.56, with the Akaike information criterion (AIC) = 95.12 and the Bayesian information criterion (BIC) = 100.03.

These results were derived from the 38 matched GC–SNe Ia pairs spanning the redshift range 0.14 ≲ z ≲ 0.89. These baseline constraints already show that neither a CDDR violation nor a significant redshift evolution of the standardized SNe Ia absolute magnitude is required by the data. The posterior asymmetry visible in the 1D distributions arises naturally from the split-normal treatment of the cluster uncertainties and is therefore a feature of the data model rather than an MCMC artifact. The corresponding constraints are summarized in Table 2 and illustrated in Fig. 3. The quality of the fit is further illustrated by the Hubble diagram shown in Fig. 2.

Fig. 2. Results of our joint analysis. Upper panel: Hubble diagram for the 38 matched pairs. The red circles represent the distance moduli of GCs (derived from DA assuming η = 0), and the blue squares represent the matched Pantheon+ SNe Ia. The dashed line indicates the prediction of the standard flat ΛCDM model (Ωm = 0.3) for reference. Lower panel: Distance-modulus residuals (Δμ = μSNe − μCluster) as a function of redshift. The data points scatter around zero with no significant redshift dependence, supporting the validity of the CDDR and the absence of strong SNe Ia evolution. The error bars represent the asymmetric uncertainties of the cluster distance moduli. The Pantheon+ SNe Ia uncertainties are not shown individually for clarity, but are fully accounted for in the covariance matrix used in the likelihood analysis.

Fig. 3. 1D and 2D marginalized posterior distributions for the parameters M0, ε, and η. The contours represent the 68% (1σ) and 95% (2σ) confidence levels, respectively. Diagonal panels: Marginalized probability density functions, with the vertical dashed lines indicating the median values. A slight degeneracy is observed between the evolution parameter ε and the CDDR-violation parameter η, but the results remain consistent with the standard CDDR (η = 0) within 1σ, indicating no evidence for CDDR violation.

3.3. Robustness against the matching tolerance

To test the robustness of our results against the choice of the matching threshold, we repeated the full cluster-only matched-pair analysis with a stricter tolerance of ΔD_C/D_C = 3%. In this case, 37 clusters remain in the sample, and the resulting constraints are

$$\begin{array}{cc}\hfill M_0& =-19.{411}_{-0.124}^{+0.134},\hfill \\ \hfill \epsilon & =-0.{154}_{-0.595}^{+0.721},\hfill \\ \hfill \eta & =-0.{044}_{-0.289}^{+0.343}.\hfill \end{array}$$(20)

The best-fit likelihood is lnℒ_best = −43.01, with AIC = 92.03 and BIC = 96.86. This test directly probes the sensitivity of the matched-pair method to the redshift-matching criterion. These results are statistically indistinguishable from the baseline 5% case. A comparison between the 5% and 3% matching cases is presented in Table 2. Since all 38 clusters are retained for any tolerance ≥5%, the corresponding matched dataset and inferred constraints are identical to the baseline case. This shows that the conclusions are not sensitive to the exact choice of matching tolerance once the threshold is at or above 5%. This stability is also clearly illustrated in Figs. 3 and 4.

Fig. 4. Posterior constraints for the stricter ΔDC/DC ≤ 3% matched sample, showing consistency with the baseline 5% results.

3.4. Supplementary constraints from DESI BAO

We next augmented the matched-pair analysis with four DESI 2024 BAO measurements. When the sound horizon is left completely free, we obtain

$$\begin{array}{cc}\hfill M_0& =-19.{504}_{-0.101}^{+0.094},\hfill \\ \hfill \epsilon & =-0.{221}_{-0.515}^{+0.659},\hfill \\ \hfill \eta & =0.{123}_{-0.286}^{+0.333},\hfill \\ \hfill r_d& =140.{56}_{-8.11}^{+8.12}\mathrm{Mpc}.\hfill \end{array}$$(21)

This result shows that the four selected BAO points provide some additional geometric information but do not by themselves tightly constrain r_d. This is expected given the limited number of BAO points and the fact that the observable D_M/r_d primarily constrains the overall distance scale rather than the CDDR deformation directly. Most importantly, the inferred values of η and ε remain fully consistent with the baseline matched-pair analysis. If instead we impose a Gaussian Planck prior on r_d, the posterior becomes

$$\begin{array}{cc}\hfill M_0& =-19.{481}_{-0.093}^{+0.086},\hfill \\ \hfill \epsilon & =-0.{239}_{-0.519}^{+0.645},\hfill \\ \hfill \eta & =0.{079}_{-0.267}^{+0.309},\hfill \\ \hfill r_d& =147.078\pm 0.259\mathrm{Mpc}.\hfill \end{array}$$(22)

The corresponding best-fit values are lnℒ_best = −51.48, AIC = 110.95, and BIC = 117.90 for the free-r_d case, and lnℒ_best = −51.83, AIC = 111.66, and BIC = 118.61 for the Planck-prior case. The combined constraints from the baseline, tolerance, and DESI-extended analyses are summarized in Table 3. We also repeated the DESI+Planck analysis for the stricter 3% matched sample, obtaining $M_0=-19.{472}_{-0.093}^{+0.086}$, $\epsilon =-0.{249}_{-0.506}^{+0.639}$, $\eta =0.{069}_{-0.265}^{+0.292}$, and r_d = 147.084 ± 0.259, again confirming the stability of the conclusions. The consistency of the results before and after including DESI BAO data is illustrated in Fig. 5.

Fig. 5. Posterior constraints after including DESI BAO data with the sound horizon rd treated as a free parameter (prior-free case), demonstrating that the inferred parameters remain consistent with the cluster-only analysis.

3.5. Comment on spline reconstruction

To assess alternative, continuous distance representations, we explored a cubic-spline reconstruction using the full Pantheon+ sample. Our motivation was to test whether the matched-pair selection itself could be bypassed in favor of a continuous reconstruction. In practice, however, we found that the spline-based results depend sensitively on the adopted parameterization. One implementation led to an apparently tight but strongly shifted posterior for ε, while an alternative parameterization produced pathological drifts in the absolute scale. This indicates that, for the present dataset, the spline nodes compete with M₀ and ε in a way that introduces additional nonphysical degeneracies.

Because our scientific goal is to test the CDDR while preserving a transparent interpretation of the nuisance parameters, we did not adopt the spline reconstruction as part of the baseline inference. Instead, we regard it as an exploratory exercise whose main lesson is methodological: with the current sample size and error budget, the matched-pair approach is more stable and easier to interpret physically.

4. Conclusions and discussions

The original motivation for this work was to test the CDDR while simultaneously marginalizing over a possible redshift evolution of the standardized SNe Ia absolute magnitude. This remains one of the key physical aspects of the analysis: the parameters η and ε are not independent in practice, and neglecting one can bias the inference on the other.

A key aspect of our analysis is the intrinsic degeneracy between the CDDR-violation parameter η and the SNe Ia evolution parameter ε. This degeneracy is physically well motivated: a violation of photon conservation (e.g., negative η) can be observationally mimicked by an apparent brightening of SNe Ia at high redshifts (positive ε). Previous studies often fixed ε = 0, leading to apparently tighter constraints on η. However, such assumptions risk introducing biases if source evolution is present. By jointly fitting η and ε, our constraints are necessarily broader but significantly more robust against astrophysical systematics.

Our analysis can be summarized in three aspects. First, we explicitly clarify that the cluster-based distances are independent of the assumed background cosmology, but still depend on astrophysical assumptions about cluster structure and hydrostatic equilibrium. Second, we demonstrate that the matched-pair result is robust against the choice of matching tolerance: tightening the threshold from 5% to 3% changes the sample from 38 to 37 pairs but leaves the inferred posterior for η statistically unchanged. Third, we supplement the cluster-only analysis with four DESI BAO measurements which, once r_d is specified, constrain η through the same generalized reciprocity relation as the cluster D_A data and may equivalently be treated as four additional D_A measurements. We show that the central conclusion is unchanged both when r_d is kept free and when it is anchored by a Planck prior.

Our main conclusions can be summarized as follows:

1. We find no statistically significant evidence for a violation of the CDDR. In all matched-pair and DESI-augmented analyses, the parameter η remains consistent with zero at the 1σ level.

2. We find no significant evidence for redshift evolution of the standardized SNe Ia absolute magnitude. The parameter ε remains consistent with zero in all matched-pair analyses.

3. Our conclusions are robust against the matching tolerance. The stricter 3% threshold produces constraints fully consistent with the baseline 5% case.

4. The four DESI BAO points are, in our framework, equivalent to four additional angular-diameter-distance measurements (once r_d is specified) and therefore constrain η in the same way as the cluster D_A sample. Their inclusion modestly improves the geometric leverage of the analysis but does not qualitatively alter the inferred values of η and ε. These conclusions are stable irrespective of whether the DESI data are written in D_M/r_d-space or in μ-space, with the two formulations being equivalent up to a Jacobian transformation.

5. Exploratory spline reconstructions are highly parameterization-dependent in the current data regime and introduce additional nonphysical degeneracies. For this reason, the matched-pair method remains the most transparent and stable baseline method for the present study.

An additional by-product of the analysis is the calibration of M₀ in a framework that does not rely on the local distance ladder. In the baseline matched-pair analysis, we obtain $M_0=-19.{460}_{-0.124}^{+0.126}$, and the extended analyses give consistent values. This provides an independent consistency check on the standard SNe Ia calibration, although current cluster uncertainties remain too large for a precision determination competitive with the local distance ladder. We emphasize that all tested extensions, including stricter matching criteria and the inclusion of DESI BAO data, yield statistically consistent results, reinforcing the robustness of the matched-pair approach.

Future progress will result primarily from larger and cleaner cluster samples. Surveys such as eROSITA are expected to expand the number of well-characterized clusters dramatically (Merloni et al. 2012), while LSST and other time-domain surveys will vastly increase the number of SNe Ia available for matched or near-matched analyses. In that regime it may become worthwhile to revisit reconstruction-based methods. For the present data quality, however, the matched-pair framework offers the best compromise between transparency, robustness, and physical interpretability.

Data availability

The data and code underlying this article are publicly available on Zenodo at https://doi.org/10.5281/zenodo.20736413 and on GitHub at https://github.com/HUJIAN0000/CDDR-MatchedPairs

Acknowledgments

We acknowledge the use of the Pantheon+ supernova data and the GC compilation. This work is supported by Yunnan Youth Basic Research Projects 202501AT070439, the National Natural Science Foundation of China (No.12473029), Dali Expert Workstation of Rainer Spurzem, Yunnan Academician Workstation of Wang Jingxiu (202005AF150025), China Manned Space Project (No. CMS-CSST-2021-A08), Guanghe Foundation (No. ghfund202407013470), Jiangsu Funding Program for Excellent Postdoctoral Talent (20220ZB59), and China Postdoctoral Science Foundation (2022M721561).

References

Appendix A: Matched galaxy cluster and SNe Ia sample

In this appendix, we present the catalog of the 38 matched pairs used in our baseline analysis. For each GC, we list its redshift (z_cl), angular diameter distance (D_A) with asymmetric uncertainties, and the derived distance modulus (μ_cl). The matched SNe Ia properties, including redshift (z_SN), SNe Ia name (CID), and corrected apparent magnitude (m_B), are also provided.

All Tables

All Figures

	Fig. 1. Redshift mismatch (Δz = zcluster − zSN) between the paired GCs and SNe Ia used in the baseline matched sample.
In the text

Fig. 2. Results of our joint analysis. Upper panel: Hubble diagram for the 38 matched pairs. The red circles represent the distance moduli of GCs (derived from DA assuming η = 0), and the blue squares represent the matched Pantheon+ SNe Ia. The dashed line indicates the prediction of the standard flat ΛCDM model (Ωm = 0.3) for reference. Lower panel: Distance-modulus residuals (Δμ = μSNe − μCluster) as a function of redshift. The data points scatter around zero with no significant redshift dependence, supporting the validity of the CDDR and the absence of strong SNe Ia evolution. The error bars represent the asymmetric uncertainties of the cluster distance moduli. The Pantheon+ SNe Ia uncertainties are not shown individually for clarity, but are fully accounted for in the covariance matrix used in the likelihood analysis.

In the text

Fig. 3. 1D and 2D marginalized posterior distributions for the parameters M0, ε, and η. The contours represent the 68% (1σ) and 95% (2σ) confidence levels, respectively. Diagonal panels: Marginalized probability density functions, with the vertical dashed lines indicating the median values. A slight degeneracy is observed between the evolution parameter ε and the CDDR-violation parameter η, but the results remain consistent with the standard CDDR (η = 0) within 1σ, indicating no evidence for CDDR violation.

In the text

	Fig. 4. Posterior constraints for the stricter ΔDC/DC ≤ 3% matched sample, showing consistency with the baseline 5% results.
In the text

	Fig. 5. Posterior constraints after including DESI BAO data with the sound horizon rd treated as a free parameter (prior-free case), demonstrating that the inferred parameters remain consistent with the cluster-only analysis.
In the text