Testing the Cosmic Distance Duality Relation Using Strong Gravitational Lensing Time Delays and Type Ia Supernovae

SGL occurs when light rays from a distant source are deflected by an intervening massive object, creating multiple images of the source. The travel time of light from the source to the observer depends on both the path length and the gravitational potential it traverses. For a single lens plane, we can quantify the extra time a light ray takes compared to an unlensed path using the excess time delay function (Wong et al., 2020):

where $\boldsymbol{\theta}$ represents the image position in the lens plane, $\boldsymbol{\beta}$ is the source position, $D_{\Delta t}$ is the time-delay distance, and $\psi(\boldsymbol{\theta})$ denotes the lens potential. The time-delay distance is a composite distance measure defined as (Refsdal, 1964; Schneider et al., 1992; Suyu et al., 2010)

where $z_{d}$ is the lens redshift, and $D_{d}$, $D_{s}$, and $D_{ds}$ are the angular diameter distances to the lens, to the source, and between the lens and source, respectively.

In systems with multiple images, we can measure the time delay $\Delta t_{ij}$ between two images $i$ and $j$ (Wong et al., 2020):

This equation encapsulates the difference in excess time delays between the two images.

By carefully monitoring the flux variations of multiple images and constructing accurate lens models, the time delays and $D_{\Delta t}$ can be measured. This approach provides a unique method for probing cosmological parameters, especially $H_{0}$, independent of other techniques like the cosmic distance ladder.

In our work, seven systems were considered initially. However, one system with a redshift of 2.375 was excluded as it exceeded the maximum redshift (z = 2.26) of the SNe Ia sample, making direct calibration impossible. Our analysis utilizes a carefully selected sample of six well-studied SGL systems with precisely measured time delays between the lensed images released by the H0LiCOW collaboration: B1608+656 (Suyu et al., 2010; Jee et al., 2019), RXJ1131-1231 (Chen et al., 2019; Suyu et al., 2013), HE 0435-1223 (Wong et al., 2017; Chen et al., 2019), SDSS 1206+4332 (Birrer et al., 2019), WFI2033-4723 (Rusu et al., 2020), and PG 1115+080 (Chen et al., 2019).

These SGL systems have been meticulously monitored over several years, resulting in high-precision time delay measurements (Wong et al., 2020). A key strength of this dataset is that each system has been individually modeled, taking into account the specific mass distribution of the lens galaxy (Wong et al., 2020). This tailored approach significantly reduces potential biases in the measurements and provides robust estimates of both the time delays and their associated uncertainties.

Integrating the SGL data with the SNe Ia dataset for the DDR test presents a challenge due to the complex nature of the SGL time delay likelihood functions. Each SGL system has its own distinct likelihood, making the consideration of SNe Ia uncertainties difficult. To address this, we developed a modified version of the likelihood analysis based on the publicly available H0LiCOW code¹¹1https://github.com/shsuyu/H0LiCOW-public. Our adaptations to this code allowed for the appropriate integration of SNe Ia error contributions into the SGL time delay likelihoods. The theoretical formulation of these modified likelihood functions is detailed in Appendix B of this paper, providing a transparent and reproducible methodology for our analysis.

In this paper, we employ the SNe Ia Pantheon+ compilation, which contains 1701 light curves of 1550 unique SNe Ia in the redshift range $0.001<z<2.26$ (Brout et al., 2022), to test the DDR. Compared to the original Pantheon compilation (Scolnic et al., 2018), the sample size of the Pantheon+ compilation has greatly increased, and the treatments of systematic uncertainties in redshifts, peculiar velocities, photometric calibration, and intrinsic scatter model of SNe Ia have been improved.

It should be noted that not all the SNe Ia in the Pantheon+ compilation are included in the Pantheon+ compilation. Since the sensitivity of peculiar velocities is very large at low redshift ($z<0.008$), which may lead to biased results, we adopt the processing treatments used by Brout et al. (2022), namely, removing the data points in the redshift range $z<0.01$.

For each SN Ia, the observed distance modulus is given by

where $m_{B}$ is the observed magnitude in the rest-frame B-band. The theoretical distance modulus $\mu_{\mathrm{th}}$ is defined as

where ${D_{L}(z)}$ is the luminosity distance.

To reconstruct the distance modulus $m_{B}$ of SNe Ia, we employ a nonparametric reconstruction technique based on artificial neural networks (ANN). The ANN method, implemented using the REFANN (Wang et al., 2020) Python code, allows us to reconstruct a function from data without assuming a specific parameterization. This approach has been widely applied in various cosmological studies (Qi et al., 2023; Dialektopoulos et al., 2022; Benisty et al., 2023).

REFANN offers several key advantages that are crucial for our analysis of the cosmic distance duality relation. As a data-driven approach, REFANN can reconstruct the distance–redshift relation without assuming any specific functional form, allowing us to capture potential anomalies or subtle features that might indicate violations of the distance duality relation. It provides robust error estimation at interpolated points, which is essential for our analysis of potential DDR violations.

While other nonparametric reconstruction methods, such as Gaussian process (GP) regression, have been widely used in cosmological studies (Seikel et al., 2012a, b; Benisty et al., 2023; Qi et al., 2019a, 2023), REFANN offers unique benefits for our specific research goals. Unlike GP, which requires the specification of a kernel function and can be computationally intensive for large datasets, REFANN adapts more flexibly to the underlying data structure without such constraints. Furthermore, REFANN has demonstrated superior performance in handling non-Gaussian errors and capturing complex, nonlinear relationships in cosmological data, which are crucial aspects of our DDR analysis.

REFANN’s flexibility in handling complex data relationships and its proven effectiveness in cosmological studies make it particularly well-suited for our research goals. These characteristics enable us to perform a more comprehensive and reliable test of the cosmic distance duality relation compared to traditional smoothing methods and other nonparametric approaches.

The optimal ANN model used for reconstructing $m_{B}$ consists of one hidden layer with 4096 neurons in total. By training the ANN on the Pantheon+ dataset, we obtain a data-driven reconstruction of the distance modulus as a function of redshift. This reconstruction enables us to directly extract the luminosity distance $D_{L}(z)$ at the redshift of each strong gravitational lensing sample, facilitating a direct comparison with the angular diameter distances $D_{A}(z)$ derived from the lensing observations.

The ANN reconstruction of $m_{B}$ provides a flexible and model-independent approach to estimate the luminosity distances of SNe Ia. By using the expressive power of neural networks, the ANN method can capture complex nonlinear relationships between the distance modulus and redshift, allowing for a more accurate representation of the underlying cosmological distance–redshift relation.

Furthermore, the training and optimization processes of ANN can introduce additional uncertainties and sensitivities. The choice of hyperparameters, such as the network architecture and training strategy, can impact the flexibility and generalization capabilities of the ANN model. These factors may contribute to the observed variations in the reconstructed distance moduli and the associated confidence regions.

Despite these challenges, the ANN reconstruction method provides a valuable tool for testing the DDR using SNe Ia and strong gravitational lensing data. By reconstructing the distance modulus $m_{B}$ in a nonparametric method, we can obtain a model-independent estimate of the luminosity distances, which can be directly compared with the angular diameter distances derived from lensing observations. This approach allows us to probe the validity of the DDR and explore potential deviations from the standard cosmological model without relying on specific parametrizations or assumptions about the underlying cosmology. The reconstruction of the distance modulus $m_{B}$ is shown in Figure 1.

To investigate potential violations of the DDR, we introduce a parameterization that allows for deviations from the standard relation,

where $\eta(z)$ quantifies the deviation from the standard DDR. We consider three different parameterizations for $\eta(z)$ (Yang et al., 2019):

where $\eta_{0}$ is a constant that quantifies the magnitude of the potential DDR violation. In the standard cosmological model, $\eta_{0}=0$, and any statistically significant deviation from zero would indicate a violation of the DDR.

With the reconstructed $m_{B}$ by ANN method, we extract the $m_{B}$ values at the redshifts corresponding to the SGL systems and convert them to luminosity distances ($D_{L}$), associating the $\eta(z)$ parameter and the absolute magnitude of SNe Ia ($M_{B}$). This allows us to directly compare the $D_{L}$ values derived from SNe Ia with the angular diameter distances ($D_{A}$) obtained from the SGL measurements.

Using the $D_{L}$ obtained from supernova observations, we derive the $D_{A}$ by applying DDR with a potential violation parameter $\eta(z)$:

For each strong lensing system, we use these converted $D_{A}$ values to calculate the angular diameter distances to the lens ($D_{d}=D_{A}(z_{d})$) and to the source ($D_{s}=D_{A}(z_{s})$). We then compute the ratio of the lens-source distance to the source distance ($D_{ds}/D_{s}$) using the equation

Finally, we calculate a theoretical time-delay distance ($D_{\Delta t}$) using

This derived $D_{\Delta t}$ is then compared with the observed time-delay distance from strong lensing measurements. The detailed mathematical formulations for converting luminosity distances to angular diameter distances and time-delay distances, as well as the associated error propagation calculations, are provided in Appendix A. These equations form the foundation of our analysis, allowing us to test the cosmic distance duality relation using the combination of strong gravitational lensing time delays and SNe Ia data.

When $M_{B}$ is allowed to vary as a free parameter, we find interesting correlations between $\eta_{0}$ and $M_{B}$ for all three parameterizations. Figure 1 shows the two-dimensional contour plots and one-dimensional posterior distributions for $\eta_{0}$ and $M_{B}$ for each parameterization.

For the $P_{1}$ model, where $\eta(z)=1+\eta_{0}z$, our analysis yields the following constraints: $\eta_{0}=0.20^{+0.11}_{-0.21}$ and $M_{B}=-18.86^{+0.21}_{-0.56}$ mag. The constraint on $\eta_{0}$ is consistent with zero at approximately 1$\sigma$ level, suggesting no strong evidence for a violation of the DDR. The constraint on $M_{B}$ is notably different from the commonly accepted value of approximately $-19.3$ mag. This shift could be attributed to the correlation between $M_{B}$ and $\eta_{0}$, as allowing for DDR violations affects the luminosity distance calculations and consequently the inferred absolute magnitude of SNe Ia.

For the $P_{2}$ model, where $\eta(z)=1+\eta_{0}\frac{z}{1+z}$, we obtain the following constraints: $\eta_{0}=0.19^{+0.17}_{-0.23}$ and $M_{B}=-19.00^{+0.22}_{-0.56}$ mag. Similar to the $P_{1}$ model, the constraint on $\eta_{0}$ suggests no strong evidence for a violation of the DDR. The constraint on $M_{B}$ is closer to the commonly accepted value of approximately $-19.3$ mag compared to the $P_{1}$ model, but still shows a slight shift toward dimmer absolute magnitudes.

The similarity in the central values of $\eta_{0}$ between the $P_{1}$ and $P_{2}$ models ($0.2$ and $0.19$, respectively) suggests a degree of robustness in the detected trend toward a small positive DDR violation. However, the large error bars for both $\eta_{0}$ and $M_{B}$ across the two models indicate that this trend is not statistically significant. These results emphasize the importance of considering multiple parameterizations when testing for DDR violations, as different models can provide complementary insights into the nature of any potential deviations from the standard DDR.

For the $P_{3}$ model, where $\eta(z)=1+\eta_{0}\log(1+z)$, our analysis yields the following constraints: $\eta_{0}=0.17^{+0.18}_{-0.23}$ and $M_{B}=-19.08^{+0.31}_{-0.49}$ mag. These results provide further insights into potential DDR violations and complement our findings from the $P_{1}$ and $P_{2}$ models. The constraint on $\eta_{0}$ remains consistent with zero within approximately 1$\sigma$, continuing the trend observed in the previous models.

Notably, all three parameterizations are consistent with $\eta_{0}=0$ within 1$\sigma$ uncertainties, providing support for the standard cosmological model and the validity of the DDR. However, the relatively large uncertainties in $\eta_{0}$ leave room for potential small deviations from the standard DDR. These results underscore the importance of considering multiple parameterizations in DDR violation tests. All three models point toward similar conclusions, demonstrating the robustness of our DDR violation detection capabilities.

The observed correlations between $\eta_{0}$ and $M_{B}$ have important implications for cosmological studies. They suggest that any potential violation of the DDR could affect our estimates of the intrinsic luminosity of SNe Ia, which in turn could impact measurements of cosmological parameters, including the Hubble constant. This interaction between DDR violations and standard candle calibrations underscores the need for careful consideration of these effects in precision cosmology.

An interesting feature of our results when $M_{B}$ is treated as a free parameter is the notably nonsmooth nature of the contours in the $\eta_{0}$–$M_{B}$ plane, as shown in Figure 2. This irregularity is particularly striking and warrants further discussion. To ensure that the nonsmooth contours in the $\eta_{0}$–$M_{B}$ plane are not due to insufficient MCMC sampling, we performed extensive convergence tests. We tested the MCMC sampling with 10,000, 20,000, and 40,000 points. The persistence of the nonsmooth features across these tests, with consistent results regardless of the number of sampling points, confirms that they are inherent to the likelihood surface and not artifacts of the sampling process. This robustness check strengthens our confidence in the complex relationship observed between $\eta_{0}$ and $M_{B}$, suggesting it is a genuine feature of our model rather than a sampling artifact.

The multimodal structure in the $\eta_{0}$–$M_{B}$ plane may be fundamentally connected to the absolute distance calibrations in our analysis. The time-delay distance $D_{\Delta t}$ is a combination of three distances (the angular diameter distances between observer-lens, lens-source, and observer-source). The model prediction $\mu_{D_{\Delta t}}$ in our likelihood analysis is calibrated using SNe Ia distances, which critically depends on the absolute magnitude $M_{B}$. Without proper $M_{B}$ calibration, we cannot obtain absolute distances, and consequently cannot derive the model prediction for time-delay distance.

Meanwhile, the observational time-delay distance $D_{\Delta t}$ from strong lensing measurements also contains absolute distance information. When constructing the likelihood using the residuals between these two absolute distances, the multimodal structure could emerge if there exists some tension between these distance measurements. Specifically, the high-density regions in the posterior distribution might represent different possible reconciliations between the distance scales from SNe Ia and strong lensing observations.

The multimodal structure could indicate either a potential systematic difference in the distance scales derived from these two distinct observational probes, or a genuine deviation from the standard distance duality relation, as parameterized by $\eta_{0}$. We cannot definitively distinguish between these possibilities with current data alone.

The complexity in the $\eta_{0}$–$M_{B}$ posterior is therefore not merely a statistical feature but could be revealing important physical information about the consistency between different distance measurements in cosmology. This interpretation suggests that the observed structure warrants further investigation with additional independent distance indicators to better understand its physical origin.

To further investigate the robustness of our constraints on DDR violations, we performed a second analysis with $M_{B}$ fixed to $-19.253\pm 0.027$ mag from the Cepheid-calibrated distance ladder method (Riess et al., 2022). Figure 3 shows the posterior distributions of $\eta_{0}$ for each of the three parameterizations under this fixed $M_{B}$ scenario. The results for each parameterization are $\eta_{0}=-0.02\pm 0.13$ for $P_{1}$, $\eta_{0}=-0.03^{+0.23}_{-0.54}$ for $P_{2}$, and $\eta_{0}=-0.19^{+0.24}_{-0.37}$ for $P_{3}$.

While $P_{1}$ and $P_{2}$ yield results very close to zero, $P_{3}$ shows a more pronounced negative value but is still consistent with zero within 1 $\sigma$ uncertainties. This suggests that the logarithmic form of $P_{3}$ might be more sensitive to potential DDR violations when $M_{B}$ is fixed. Notably, all three parameterizations yield results consistent with $\eta_{0}=0$ within 1$\sigma$ uncertainties, providing strong support for the validity of the standard DDR. This consistency across different models strengthens our confidence in the robustness of the DDR.

Interestingly, the central values of $\eta_{0}$ have changed from positive (when $M_{B}$ was free) to negative for all parameterizations. This sign flip highlights the strong degeneracy between $\eta_{0}$ and $M_{B}$, demonstrating how fixing $M_{B}$ can significantly impact our inferences about potential DDR violations. While fixing $M_{B}$ provides a more direct probe of DDR violations, it may also introduce biases if the fixed value is not precisely correct. The complementary information from both approaches provides a more comprehensive understanding of potential DDR violations and their cosmological implications.

A notable feature in our analysis is the increase in $\eta_{0}$ uncertainties when a prior is applied to $M_{B}$ (Table 1). This behavior is closely tied to the structure of the $\eta_{0}$–$M_{B}$ parameter space, as revealed by our contour plots. These plots show two distinct high-density regions. When $M_{B}$ is allowed to vary freely, the best-fit value of $\eta_{0}$ falls within the upper, more compact high-density region, resulting in smaller uncertainties. This region corresponds to $M_{B}$ values of approximately $-19.0$ mag or slightly larger. However, fixing $M_{B}$ to $-19.253$ mag shifts the best-fit $\eta_{0}$ to the lower, more extended high-density region, leading to larger uncertainties. This shift also explains the change in $\eta_{0}$’s central value from positive to negative or near-zero when $M_{B}$ is fixed. These results highlight the complex interplay between $\eta_{0}$ and $M_{B}$, underscoring the importance of carefully considering priors in analyses of the distance duality relation.