© The Authors 2025

1. Introduction

The determination of the Hubble constant H₀ has been one of the most relevant goals in cosmology in recent decades. It defines the universe’s current expansion rate and, therefore, sets the scale and age of the universe today. Several methods to determine H₀ have already been applied (Di Valentino et al. 2021), but there is still a > 4σ tension (Verde et al. 2019) between cosmic microwave background measurements by the Planck Collaboration VI (2020) and distance ladder measurements from the SH0ES (Supernova H₀ for the Equation of State) program (Riess et al. 2022, 2024; Li et al. 2024; Anand et al. 2024) as well as other late universe measurements. The most recent results from distance ladder calibrations with the Tip of the Red Giant Branch (TRGB) from Freedman (2021) (see also Freedman et al. 2019, 2020) are in agreement within 1.3σ with the Planck Collaboration VI (2020) and at a 1.6σ level with Riess et al. (2022). Another usage of the TRGB from Soltis et al. (2021) is in better agreement with Riess et al. (2022).

The Hubble tension hints at either new unknown physics beyond our current model of the universe or problems within the various H₀ measurements or distance calibrations. Therefore, we aim to determine H₀ without distance calibrations. Three such methods are using radiative transfer modeling of type II supernovae (SNe II) based on the tailored-expanding-photosphere method (e.g. Vogl et al. 2020; Dessart & Hillier 2009, 2006; Dessart et al. 2008), gravitational wave sources (e.g. Palmese et al. 2024; Mukherjee et al. 2021; Gayathri et al. 2021), or masers (Pesce et al. 2020). A further direct method is time-delay measurements of distant gravitationally lensed transient objects as first proposed by Refsdal (1964). Such objects can be quasars or SNe, which can show up as multiple images if they are strongly gravitationally lensed. The light rays of these images have shorter or longer light paths as they pass through different parts of the gravitational potential, leading to a time difference between the image appearances. As transient objects show features like peaks and kinks in light and color curves, we can infer the time delay between images of such objects. Together with models of the gravitational lensing system and reconstructions of the mass perturbations along the line of sight, we can calculate H₀.

This technique was already applied by the H0LiCOW (H₀ Lenses in COSMOGRAIL’s Wellspring) program (Suyu et al. 2017; Birrer et al. 2019; Sluse et al. 2019; Chen et al. 2019; Wong et al. 2020; Rusu et al. 2020) collaborating with the COSmological MOnitoring of GRAvItational Lenses (COSMOGRAIL; Eigenbrod et al. 2005; Courbin et al. 2018; Bonvin et al. 2018) and the Strong lensing at High Angular Resolution Program (SHARP; Chen et al. 2019). Combining the time-delay measurements of six quasars, the Hubble constant is determined as $H_0=73.3_{-1.8}^{+1.7}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$ assuming a flat ΛCDM cosmology and physically motivated lens mass distributions (Wong et al. 2020). This value is in agreement with late universe measurements. The STRong-lensing Insights into the Dark Energy Survey (STRIDES) collaboration (Treu et al. 2018) analyzed a seventh lensed quasar system (Shajib et al. 2020) which agrees well with Wong et al. (2020). Systematic effects are further tested, and more lensed quasars are studied by the Time-Delay COSMOgraphy (TDCOSMO) organization (Millon et al. 2020a,b; Gilman et al. 2020; Birrer et al. 2020).

Originally, Refsdal (1964) suggested using lensed supernovae (LSNe) as the transient source for time-delay cosmography. The first spatially resolved LSN called “SN Refsdal”, which is of type II and lensed by a galaxy cluster, was detected 50 years later (Kelly 2015; Kelly et al. 2016a,b). Since then, a few more LSNe have been observed (Goobar et al. 2017, 2022; Chen et al. 2023). The most recently detected LSNe are “SN H0pe” (Frye et al. 2023, 2024), which is confirmed to be of type Ia, and “SN Encore” (Pierel et al. 2024), which appeared in a lensing system that previously hosted another lensed SN, most likely of type Ia (Rodney et al. 2021). Both of these SNe are lensed by a cluster, which is so far more easily detected from deep imaging of galaxy clusters than typically shallower imaging of galaxy-scale lenses. We focus on galaxy-scale lenses in this work since they are expected to be more abundant than cluster-scale lenses at the same imaging depth in future surveys. With SN Refsdal and SN H0pe, the first time-delay cosmography measurements were performed (Kelly et al. 2023; Grillo et al. 2024; Liu & Oguri 2024; Pascale et al. 2024; Pierel et al. 2024). Lensed supernovae are not only interesting in terms of cosmography but also give insights into their explosion mechanisms, as the modeling of the gravitational potential can predict the temporal and spatial location of future images, which gives the opportunity to study them from very early on.

The performance of LSNe Ia for time-delay cosmography was studied in several papers (e.g., Dobler & Keeton 2006; Goldstein et al. 2018; Huber et al. 2019, 2022; Suyu et al. 2020; Huber & Suyu 2024; Arendse et al. 2024). SNe II were not given that much attention, even though they are the most abundant type of SNe by volume (Li et al. 2011) and we expect to detect hundreds of LSNe II with the upcoming Rubin Observatory Legacy Survey of Space and Time (LSST; LSST Science Collaboration 2009) (Wojtak et al. 2019; Goldstein et al. 2018; Goldstein & Nugent 2017). This is the reason we focus on this type in our work, in particular SNe IIP. A major concern for the correct time-delay determination are microlensing effects stemming from small and compact masses within the lens galaxy, such as stars. In paper HOLISMOKES V (Bayer et al. 2021), we found that in the early plateau phase of SNe IIP during at least the first 30 days after explosion, the specific intensity profiles are achromatic as in the case of SNe Ia as well (Goldstein et al. 2018; Huber et al. 2019, 2021). Therefore, we can neglect microlensing effects on the color curves during these phases, as the additional magnification of microlensing cancels out in the color curve.

While we focused on getting time delays from spectra in HOLISMOKES V, this work now focuses on the determination of a time delay based on the presence of a small kink in the color curves during the plateau phase (de Jaeger et al. 2018). As photometric follow-up is much less expensive than a spectroscopic one, this method should be more easily applicable for studying a large sample of objects. An arbitrary example of SN IIP color curves of two lensing images is shown in Fig. 1, illustrating the time delay Δt and a shift in color s, caused by the different extinction for the two images A and B. We use SN data from the Carnegie Supernova Project (CSP) (Hamuy et al. 2006) to construct mock color curves and determine time delays and extinction parameters (Cardelli et al. 1989) with the help of Monte Carlo Markov Chain (MCMC) sampling (Dunkley et al. 2005).

Fig. 1. Exemplary noiseless LSN IIP color curves of two lensing images A and B showing a time shift Δt and a shift in color s due to extinction.

Following the findings in Suyu et al. (2020), we would need cosmology grade time-delay measurements with uncertainties ≲5% of about 20 spatially resolved LSNe Ia with time delays longer than 20 days and including microlensing effects to determine H₀ with an uncertainty of ∼1%. We further impose the assumption of ∼3% lens mass modeling uncertainties as well as ∼3% lens environment uncertainties. This is realistic when spatially resolved kinematics of the foreground lens are available (Yıldırım et al. 2020, 2023). We expect to find more LSNe II than LSNe Ia (Wojtak et al. 2019; Goldstein et al. 2018), mostly of type IIn, during the era of LSST, and we find in HOLISMOKES V that LSNe II are less susceptible to microlensing resulting in less biased time delays, even systems with time-delay precision ≲10% will help to achieve the desired precision on the Hubble constant.

We structure our work as follows: In Sect. 2 we go into detail on dust extinction and the creation of our mock SN IIP color curves. We present our MCMC sampling in Sect. 3 and the results in Sect. 4 applying it to the rest-frame SN mock color curves. After that, we look into the impact of redshift in Sect. 5. We finish this work with a discussion and conclusion in Sect. 6.

2. Type II supernova color curves

de Jaeger et al. (2018) finds that SN II color curves show a characteristic change in their slope within the plateau phase, which lasts from (rest-frame) days of ∼10 to ∼100 after explosion. Before the kink around day 35 (rest-frame) after explosion, the color becomes redder more rapidly than afterwards. This “kink” (as shown in Fig. 1) enables us to do cosmography as well as determine extinction parameters.

In this section, we detail the creation of the mock data used in this work. In Sect. 2.1 we introduce the extinction laws by Cardelli et al. (1989) referred to as Cardelli laws. Sect. 2.2 shows the CSP light curves that we use as archetypical SN IIP light curves. We fit these light curves with orthogonal polynomial functions in Sect. 2.3 to create the mock data in Sect. 2.4.

2.1. Dust extinction and extinction parameters

The multiple images of a lensed supernova are observed in different regions of the lensing galaxy. For that reason, it is necessary to emulate the impact of dust on the colors of the lensed images. For this experiment, we use the Cardelli extinction laws, which provide a flexible model of the extinction in galaxies (Cardelli et al. 1988, 1989). The Cardelli laws consider dust absorption and reemission at longer wavelengths as well as light scattering away from the line of sight. The backscattering of light from surrounding stars into the line of sight is omitted as our source is far away. We further ignore extinction occurring in the Milky Way, as we assume it to be the same for the different lensed images owing to the small angular separation. We remark that the Cardelli laws for mock data creation can be exchanged by any other flexible extinction model in galaxies, such as the more recent extinction laws by Gordon et al. (2023), without impacting the retrieval of the time delay.

We denote the intrinsic magnitude as m_λ, 0 and the observed magnitude as m_λ at the wavelength λ. The wavelength-dependent extinction coefficient is defined as:

$$\begin{array}{c}\hfill A_{\lambda }=m_{\lambda }-m_{\lambda ,0}.\end{array}$$(1)

Extinction laws are typically defined as A_λ/A_V, for which A_V is a normalization defined in the optical at λ ≈ 5500 Å. They can also be expressed via the total-to-selective-extinction ratio (Cardelli et al. 1989):

$$\begin{array}{c}\hfill R_V=\frac{A_V}{E(B-V)}=\frac{A_V}{(A_B-A_V)},\end{array}$$(2)

depending on the color excess:

$$\begin{array}{c}\hfill E(B-V)\equiv A_B-A_V.\end{array}$$(3)

The Cardelli laws (Cardelli et al. 1988, 1989) are generally defined depending only on R_V:

$$\begin{array}{c}\hfill \frac{A_{\lambda }}{A_V}=a(x)+\frac{b(x)}{R_V}·\end{array}$$(4)

The parametric functions a(x) and b(x), where x = λ⁻¹, are determined analytically from spline fits to measured data in the three wavelength regimes, infrared (IR), optical, and ultraviolet (UV) (Cardelli & Clayton 1988; Clayton & Cardelli 1988). The extinction impacts the color curve via a magnitude shift s (see Fig. 1), which is correlated for the case of multiple color curves.

2.2. Observed light curves

Until today, only two spatially resolved lensed SNe II were observed. The lensed core-collapse SN observed in the galaxy cluster Abell 370 has a single epoch observed, which gives us no opportunity to investigate full light curves (Chen et al. 2023). “SN Refsdal” is not of type IIP, but is a 1987A-like SN II (Kelly et al. 2016a). Therefore, we have to create mock color curves based on measured light curves of unlensed SNe IIP. We based our mock color curves on observed light curves from CSP. The data used in this work from CSP (Anderson et al. 2024) are rest-frame u, g, and r band light curves (see Krisciunas et al. 2017 for the definition of the passbands) for SN2005J, which are shown in Fig. 2. We treated these filters as universally applicable to other facilities, which might be used in the future for follow-up observations in the LSST era. The corresponding rest-frame color curves u − g and u − r, used later in the study, show a strong kink mainly because of the larger brightness decrease in the u band compared to the g and r bands. The impact of redshift will be discussed later in Sect. 6.

Fig. 2. u,  g, and r light curves of SN2005J observed by the CSP (Anderson et al. 2024) including examples of orthogonal polynomial fits (solid lines) to the u, g and r band light curves.

2.3. Light curve model

We first derived a parametrized model of the light curves, which we can use to generate mock light curves and color curves of LSNe IIP. We used orthogonal polynomials as the fitting function C(η,t) that is applied to the observed light curves, where η is the vector of fitting parameters. In addition to the orthogonal polynomial fitting procedure, we tested spline fitting in Appendix A to test the sensitivity of our result to the choice of model.

For the fitting with orthogonal polynomial functions, we considered a linear combination of Legendre polynomials P_n as the functional form C(η_model,t) on an interval [t_min, t_max] up to the order of n = 6:

$$\begin{array}{c}\hfill C({\mathit{\eta }}_{\mathrm{model}},t)=c_1{P}_0+c_2{P}_1+c_3{P}_2+\dots +c_{n+1}{P}_n,\end{array}$$(5)

where P_n was defined as:

$$\begin{array}{c}\hfill P_n=\frac{1}{2^{n}n!}\frac{{\mathrm{d}}^n}{\mathrm{d}{u}^n}{(u^2-1)}^n,\end{array}$$(6)

with substitution u = (2t − t_min − t_max)/(t_max − t_min) to ensure that the Legendre polynomials are orthogonal on an interval t = [t_min, t_max] (with u = [ − 1, 1]) and P₀ = 1. The orthogonality interval [t_min, t_max] depends on the time range covered by the color curve of the lensed image observed for the longest period. In our case, we can take either the first or the second image time range, as both are of equal length. This will most likely be the second image for a real measurement, as its light curve can be followed and gathered at an earlier time than the first image. We chose the highest order, n = 6, as it is flexible enough.

The parameter vector is defined as η_mock = (c₁,c₂,c₃,c₄,c₅,c₆,c₇)_{light curve}. For the mock data creation, we fited the measured CSP data several times with C(η_mock,t) to ensure reaching the global minimum and selected the fit with the lowest χ². An example of a light curve fit with up to 6th order orthogonal polynomials to the light curves u, g, and r can be seen in Fig. 2.

2.4. Mock data creation

To generate mock light curves, we assumed systems detected by LSST (Oguri & Marshall 2010) but followed up at a high cadence at a different facility. The parametric SN model, as introduced in Sect. 2.3 was evaluated with an artificial two-day cadence, where each simulated data point can vary around the perfect two-day cadence to simulate realistic observing conditions, adding a value randomly picked from a uniform distribution [−0.1,0.1].

Uncertainties were added following the procedure of Huber et al. (2022) and LSST Science Collaboration (2009), assuming a follow-up facility with a 2 m main mirror. Specifically, at each observational time t_i and for each observational band, a random number g_i was drawn from a Gaussian with zero mean and standard deviation q_i. These numbers were then added to the value determined from the fitted function. The value of q_i was set by the signal-to-noise ratio (S/N) of the observed magnitude, which depends on C(η_mock,t_i) and the limiting magnitude, which we adopted as one magnitude deeper than the LSST-like 5σ depth (Huber et al. 2019). If C(η_mock,t_i) was more than two magnitudes fainter than the 5σ depth, the data point was not generated. Following this approach, we created the mock light curves for the first strongly lensed image A. The mock light curves of the second image, B, were created similarly. However, adding a time and magnitude shift from the filter-dependent extinction and evaluating the model at time steps coinciding with the ones from image A ensures measurements are taken simultaneously for both A and B. For each image, we subtracted the mock light curves in different filters to create the mock color curves d_A and d_B. The color curve uncertainties σ_A and σ_B were calculated by adding the light curve uncertainties in quadrature.

For the main study, we investigated three different scenarios in the mock data creation with SN2005J, which are shown in Table 1:

1. High S/N assuming the uncertainties given by unlensed low redshift observations (Anderson et al. 2024).

2. Low S/N assuming the V band light curve peak of the first image to be at 22 mag and for the second image to be at 23 mag.

3. Low S/N assuming the V band light curve peak of the first image to be at 22 mag and for the second image to be at 23 mag with additionally 35% of the data missing to simulate bad weather conditions and other incidents where observations are not possible.

We used case 1 for method verification, as its high S/N is unrealistic. To mock up case 2 and case 3, we assumed a typical peak absolute brightness at −17 mag in the rest-frame V band (Anderson et al. 2014) for SNe II. In rare cases, the peak absolute magnitude can reach −18 mag. Considering a flat ΛCDM cosmology with H₀ = 72 km s⁻¹ Mpc⁻¹ (Bonvin et al. 2017), Ω_m = 0.26 and Ω_Λ = 0.74 (Oguri & Marshall 2010) and a typical source redshift of lensed SNe II between z = 0.3 and z = 0.6 (Oguri & Marshall 2010), we get the apparent magnitude range of 23.9–25.7 mag. The lens further magnifies this with a typical magnification of μ ≈ 3 for the first image and μ ≈ 1 for the second image (Oguri & Marshall 2010). The first image is, therefore, magnified to around 22.7 magnitudes for the low redshift end, and the second image is not magnified. To keep our cases 2 and 3 realistic but also feasible in terms of observational exposure time, we chose the V band light curve peak of the first image to be at 22 mag and for the second image at 23 mag and adjusted the u, g, and r accordingly, which is still realistic for bright SNe II at the lower redshift boundary. This could also be achievable by observing with a follow-up facility with a main mirror larger than 2 m or deeper 5σ depth than just one magnitude deeper than LSST 5σ depth. To quantify the data loss for scenario 3, we analyzed two years of observations of the quasar RXJ1131-1231 at the 1.2 meter Leonhard Euler Telescope with a two-day cadence from COSMOGRAIL. We find that a loss of ∼35% of observations because of bad weather, downtime of the observatory, or the moon phase is a realistic approximation (Millon et al. 2020b). Reddening due to dust was included, as discussed in Sect. 3.

3. Time-delay and extinction retrieval procedure

The light curve model in Eq. (5) will also be used in the MCMC sampling process in this section, using lower orders than the full 6 orders. For the time-delay and extinction parameter retrieval, we use a custom Metropolis-Hastings sampling algorithm (Metropolis et al. 1953) based on the approach of Dunkley et al. (2005) to ensure full control of the sampling. In this section, we introduce our algorithm.

3.1. Likelihood function

The likelihood function assigns a probability $p(\mathit{\eta },\mathrm{\Delta }t,A_{V,\mathrm{B}}-A_{V,\mathrm{A}})$ to the sampling parameters $\mathit{\eta }={(c_1,c_2,c_3,c_4,c_5,c_6,c_7)}_{\mathrm{color}\mathrm{curve}},\mathrm{\Delta }t,A_{V,\mathrm{B}}-A_{V,\mathrm{A}}$, and for Gaussian uncertainties the likelihood can be expressed as L = A_normexp(−χ²/2) with A_norm being a normalization factor. During the sampling of our MCMC, as stated in Sect. 3.2 we were just interested in minimizing the χ² function and do not regard A_norm. To define χ² in the most general form, we defined the one-dimensional data vectors d_A and d_B containing the data points of the considered color curves of image A and B. We further considered the model data vectors ${\mathit{d}}_{\mathrm{A}}^{\mathrm{model}}$ and ${\mathit{d}}_{\mathrm{B}}^{\mathrm{model}}$. The entries of these one-dimensional vectors were defined combining C(η,t_j) (Eq. 5) into vector form:

$$\begin{array}{c}\hfill d_{\mathrm{A},j}^{\mathrm{model}}=C(\mathit{\eta },t_j),\end{array}$$(7)

and

$$\begin{array}{c}\hfill d_{\mathrm{B},j}^{\mathrm{model}}=C(\mathit{\eta },t_j-\mathrm{\Delta }t)+s,\end{array}$$(8)

where Δt is the time shift, which we sampled together with the parameter vector η and extinction parameters that are related to s. We wrote a general shift s in color between images A and B into this definition, which can be expressed by extinction parameters. Considering generic filters x and y of a color curve x − y, we define the magnitudes m of each image A and B and each filter and get:

$$\begin{array}{c}\hfill s=(m_{x,\mathrm{B}}-m_{y,\mathrm{B}})-(m_{x,\mathrm{A}}-m_{y,\mathrm{A}}).\end{array}$$(9)

For each image (A or B) and filter (x or y), we can write the observed magnitude as:

$$\begin{array}{c}\hfill \begin{array}{c}{m}_{x,\mathrm{A}}=m_{0,x}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{A}}+A_{x,\mathrm{A}},\hfill \\ m_{y,\mathrm{A}}=m_{0,y}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{A}}+A_{y,\mathrm{A}},\hfill \\ m_{x,\mathrm{B}}=m_{0,x}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{B}}+A_{x,\mathrm{B}},\hfill \\ m_{y,\mathrm{B}}=m_{0,y}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{B}}+A_{y,\mathrm{B}},\hfill \end{array}\end{array}$$(10)

where m₀ is the intrinsic magnitude, Δm_lens is the change from lensing magnification, and A is the extinction coefficient. Even though both images have different lensing magnifications, they have the same Δm_lens across filters (i.e. Δm_lens, A = Δm_{lens, x, A} = Δm_{lens, y, A} and Δm_lens, B = Δm_{lens, x, B} = Δm_{lens, y, B}) as lensing is achromatic, neglecting chromatic microlensing effects, which we expect to be small given that the size of photosphere will not change substantially during the plateau phase (Dessart & Hillier 2011), which was also confirmed in our earlier microlensing studies in HOLISMOKES V (Bayer et al. 2021) until around day 30 after explosion. After this epoch, chromatic microlensing might be stronger but still smaller than for SNe Ia. We neglected it in this work, even for later epochs, because of the stable size of the photosphere, but further investigations of microlensing for later epochs of SNe IIP could be of interest in future studies. The additional magnification cancels out, and the equations from Eq. (10) can be inserted into Eq. (9):

$$\begin{array}{cc}\hfill s=& (m_{0,x}-m_{0,y}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{B}}+\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{B}}+A_{x,\mathrm{B}}-A_{y,\mathrm{B}})\hfill \\ \hfill & -(m_{0,x}-m_{0,y}-\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{A}}+\mathrm{\Delta }{m}_{\mathrm{lens},\mathrm{A}}+A_{x,\mathrm{A}}-A_{y,\mathrm{A}})\hfill \\ \hfill =& A_{x,\mathrm{B}}-A_{y,\mathrm{B}}-A_{x,\mathrm{A}}+A_{y,\mathrm{A}}\hfill \\ \hfill =& [(a({\lambda }_x)+\frac{b({\lambda }_x)}{R_{V,\mathrm{B}}})-(a({\lambda }_y)+\frac{b({\lambda }_y)}{R_{V,\mathrm{B}}})]·A_{V,\mathrm{B}}\hfill \\ \hfill & -[(a({\lambda }_x)+\frac{b({\lambda }_x)}{R_{V,\mathrm{A}}})-(a({\lambda }_y)+\frac{b({\lambda }_y)}{R_{V,\mathrm{A}}})]·A_{V,\mathrm{A}}\hfill \\ \hfill =& [(a({\lambda }_x)+\frac{b({\lambda }_x)}{R_V})-(a({\lambda }_y)+\frac{b({\lambda }_y)}{R_V})]·(A_{V,\mathrm{B}}-A_{V,\mathrm{A}}).\hfill \end{array}$$(11)

For the last two steps, we inserted the Cardelli extinction laws from Eq. (4) and fixed R_V, A = R_V, B = R_V = 3.1 which is the standard value for the diffuse interstellar medium (Cardelli et al. 1989). The effective wavelengths used to calculate the parameters a(λ) and b(λ) for the CSP filters are: ${\lambda }_u=3639.3\mathrm{\AA }$, ${\lambda }_g=4765.1\mathrm{\AA }$, and ${\lambda }_r=6223.3\mathrm{\AA }$ (Krisciunas et al. 2017).

The final definition of the model data vector for image B was expressed as follows, where A_V, B − A_V, A is the differential extinction parameter:

$$\begin{array}{c}\hfill \begin{array}{c}{d}_{\mathrm{B},i}^{\mathrm{model}}=\mathit{C}(\mathit{\eta },t_i-\mathrm{\Delta }t)+[(a({\lambda }_x)+\frac{b({\lambda }_x)}{R_V})\hfill \\ -(a({\lambda }_y)+\frac{b({\lambda }_y)}{R_V})]·(A_{V,\mathrm{B}}-A_{V,\mathrm{A}}).\hfill \end{array}\end{array}$$(12)

The χ² was defined in the most general form as:

$$\begin{array}{cc}\hfill {\chi }^2& ={\chi }_{\mathrm{A}}^2+{\chi }_{\mathrm{B}}^2\hfill \\ \hfill & ={({\mathit{d}}_{\mathrm{A}}-{\mathit{d}}_{\mathrm{A}}^{\mathrm{model}})}^T{C}_{\mathrm{D},\mathrm{A}}^{-1}({\mathit{d}}_{\mathrm{A}}-{\mathit{d}}_{\mathrm{A}}^{\mathrm{model}})\hfill \\ \hfill & +{({\mathit{d}}_{\mathrm{B}}-{\mathit{d}}_{\mathrm{B}}^{\mathrm{model}})}^T{C}_{\mathrm{D},\mathrm{B}}^{-1}({\mathit{d}}_{\mathrm{B}}-{\mathit{d}}_{\mathrm{B}}^{\mathrm{model}}),\hfill \end{array}$$(13)

where C_D is the covariance matrix. For the simple case of only one color curve, the covariance matrix had only non-zero entries on the main diagonal, where the uncertainties of each color curve data point were calculated by adding the uncertainties of the two respective light-curve data points in quadrature. This also held for using multiple color curves containing only different filters, e.g., for the two color curves u − g and r − i, as it is a good approximation to assume these filters do not correlate. For the case of using color curves containing the same filter, e.g. the two color curves u − g and u − r, we also have entries on the diagonals within the off-diagonal blocks of the matrix, which are not zero. The inverse of the covariance matrix is calculated using the Cholesky decomposition of the numpy.linalg python module (Harris et al. 2020).

3.2. Sampling technique and prior

We started to sample the posterior probability distribution with initial parameters ${\mathit{\eta }}_0,\mathrm{\Delta }{t}_0$ and (A_V, B − A_V, A)₀ selected with a χ² non-linear fit of the mock color curves based on a Levenberg-Marquart algorithm as implemented in the optimize.curve_fit function of the python module scipy (Virtanen et al. 2020). The next set of parameters η_i of step i was selected by the walker setting the step size. For the initial chain, we defined the trial probability distribution for taking the next step as a multivariate Gaussian with the identity matrix as the covariance matrix. For the subsequent chain, we used the covariance matrix calculated from the previous chain and scaled it by a constant factor to obtain the trial probability distribution with the desired acceptance rate. This ensured that we had an ideal step size for the walker and an acceptance rate around ∼25% after several chains. We applied a broad uniform prior [−100,100] to the time delay. All other parameters had improper uniform priors. To sample the parameters $\}\mathit{\eta },\mathrm{\Delta }t,A_{V,\mathrm{B}}-A_{V,\mathrm{A}}\{$, we followed the Metropolis-Hastings algorithm (Dunkley et al. 2005). We performed a check for convergence as described in Sect. 3.3, and we kept running the chain if convergence was not achieved.

3.3. Convergence testing

During the convergence test, the power spectrum of the samples of each parameter was fitted with the parametric function

$$\begin{array}{c}\hfill P(k)=P_0\frac{{(N{k}^{\ast }/2\pi j)}^{\alpha }}{1+{(N{k}^{\ast }/2\pi j)}^{\alpha }}\end{array}$$(14)

from Dunkley et al. (2005). It considered the chain length N, the wave number k^* at the turn-over of the power spectrum, the power P₀ at wave number k = 0, the Fourier modes j = k(N/(2π)), and an exponent α. A chain converged if the following two criteria were met:

1. j^* > 20, with j^* = k^*(N/(2π)) as the turn-over Fourier mode. This ensured that the minimum wave number k_min was in the regime of white noise P(k)∼k⁰ indicating that the region of high probability was fully sampled.

2. The convergence ratio r for normalized chains should obey r = (P₀/N) < 0.01.

If one of the criteria was not met, we explored an additional 20 000 steps before testing again and repeated this process until convergence was reached. The first 20% of the converged chain was removed as the burn-in phase before convergence testing.

3.4. Polynomial order selection with the Bayesian Information Criterion

We investigated several noise realizations for each mock case to obtain statistically significant results. For case 1 of Table 1, we used 100 noise realizations, and for cases 2 and 3, we used 500 noise realizations, as fluctuations increase with lower S/N. Each noise realization was fitted with different orders of the orthogonal polynomials as presented in Section 2.3, ranging from the linear order 1 to the most complex case of order 6. To select which order to include in the final results, we applied the Bayesian Information Criterion (BIC) to each noise realization and fitted order. The BIC is a simplified criterion for Bayesian model comparison and is defined as:

$$\begin{array}{c}\hfill \mathrm{BIC}=ln(N_{\mathrm{data}})N_{\mathrm{par}}-2ln(\widehat{L}),\end{array}$$(15)

with the number of data points N_data, number of free parameters N_par and maximum likelihood $\widehat{L}$ of the noise realization (Schwarz 1978). Choosing the model with the lowest BIC ensures that not only is the best fitting model considered, but also the complexity of each model is taken into account. For the final combination of all investigated noise realizations, we considered the polynomial order of fitting with the lowest BIC value for each realization.

4. Time-delay and extinction retrieval results

In this section, we present the precision and accuracy of our time shift and differential extinction parameter determination method with lensed SN IIP color curves. The input values in the rest-frame mock data creation for all cases in this section were Δt = 16.3 days and A_V, B − A_V, A = 0.2 mag, which is a realistic value based on the study of the extinction properties of mostly elliptical lensing galaxies (Elíasdóttir et al. 2006). For each of the scenarios 1, 2, and 3 introduced in Sect. 2.4 and Table 1 we investigated the single color curve u − g as well as the combinations of two color curves u − g and u − r. In Appendices A, B, C, and D, we investigated additional scenarios to test the robustness of our method, including other time-shift values.

Before presenting the results, we demonstrate the selection of the orders via the BIC for one case 1 and one case 3 noise realization. All figures about the BIC selection are shown in Appendix E. The high S/N case 1 mock data example with fits of orders 6 down to 2 are shown in Fig. E.1, where we highlight the selected order with a red frame. Further, we list the corresponding BIC values in Table E.1. It is visible that the 5th order fits as well as the 6th order and is favored in the BIC calculation. The lower orders visibly fit worse. For case 3 we show the same orders fit to the mock data in Fig. E.2 as well as the BIC values in Table E.1. For this case, it is not clearly visible which order to favor, but considering the calculated BIC value, the 4th order is selected for this specific noise realization. For each case, this procedure is applied to all noise realizations, resulting in a combination of fits with different orders for the full data set.

The corner plots for mock data case 1 for one color curve u − g and two color curves u − g and u − r are shown in Fig. 3. Cases 2 and 3 are shown in Figs. 4 and 5, respectively. In Fig. 3, both corner plots share the same scale. The corner plots in Figs. 4 and 5 have different scaling, also compared to Fig. 3. The time shifts and differential dust extinctions with uncertainties for all cases are additionally listed in Table 2.

Fig. 3. Corner plots of the MCMC chain results and bar plots of selected orders for 100 noise realizations for case 1. The left panels show corner plots of the time shift and differential dust extinction retrieved from 100 noise realizations for case 1. We plot the contours for the 1σ (68%, solid red line) and 2σ (95%, dashed red line) credible intervals for the differential dust extinction AV, B − AV, A and the time delay Δt. The input values in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag. On the right, we show histograms of the selected orthogonal polynomial orders of fitting using the BIC for the corresponding cases on the left. Panel a: Case 1 with one color curve u − g. Panel b: Orders of fitting selected for case 1 with one color curve u − g. Panel c: Case 1 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 1 with two color curves u − g and u − r. We impose the same order for both color curves.

Fig. 4. Same as Fig. 3, but for 500 noise realizations for case 2. Panel a: Case 2 with one color curve u − g. Panel b: Orders of fitting selected for case 2 with one color curve u − g. Panel c: Case 2 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 2 with two color curves u − g and u − r. We impose the same order for both color curves.

Fig. 5. Same as Fig. 3, but for 500 noise realizations for case 3. Panel a: Case 3 with one color curve u − g. Panel b: Orders of fitting selected for case 3 with one color curve u − g. Panel c: Case 3 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 3 with two color curves u − g and u − r. We impose the same order for both color curves.

Considering case 1 for one and two color curves, we can retrieve the time shift and the differential extinction perfectly with relative uncertainties down to the 1% level for typical delays and extinctions, which proves the functionality of the method. Furthermore, combining multiple color curves helps with constraining both the time delay and the differential dust extinction. For high S/N data in case 1, only high orders of the orthogonal polynomial fitting are selected via the BIC, as shown in Figs. 3b and 3d.

When adjusting the S/N down for case 2, the result gets less precise, as expected. Considering only the u − g color curve, we cannot retrieve the time delay precisely enough to be of cosmology grade as defined in Sect. 1, unless the delays are longer than ∼35 days, to reach an uncertainty ≲10% per lens system. This gets even less precise for the most realistic case 3, including 35% data loss, where uncertainties reach ∼5–6 days. A_V, B − A_V, A cannot be constrained to better than ∼0.7 mag with a single color curve in low S/N cases. By including the second color curve u − r, the constraint on the differential dust extinction is improved substantially to < 0.1 mag in precision, which helps break the degeneracy between Δt and A_V, B − A_V, A and decreases the time delay uncertainty down to ∼1 day. For the lower S/N cases with one color curve, lower orders of fitting are selected via the BIC because the higher uncertainties of the mock data leave more freedom in the functional form of the fit. This gets shifted to higher orders again when a second color curve is taken into account, as it sets more constraints on the fit.

5. Inference for non-zero redshift

Up to now, we investigated rest-frame photometry in the CSP filters u, g, and r. In reality, we would observe most lensed SNe II at redshifts around z = 0.6 (Oguri & Marshall 2010). Assuming a redshift of z = 0.6, the rest-frame u, g, and r bands would shift roughly into the observed r, i, and z bands, respectively. In this section, we investigate how the shape of the kink is affected with the help of K-corrections and the rest-frame UV color curves of the Ultra Violet Optical Telescope (UVOT) on the Neil Gehrels Swift Observatory (SWIFT) (Gehrels et al. 2004; Roming et al. 2005) SNe, which would be approximately observed in the u and g band for z = 0.6.

5.1. K-corrections

For applying K-corrections, we need spectra, but for the CSP and SWIFT data, we have no spectra available to cover the rest-frame UV. We, therefore, use publicly available model spectra of SN IIP with temporal evolution created by Dessart et al. (2013)¹. We applied the K-correction to the photometry according to Hogg et al. (2002) by adding K to the rest-frame magnitude in the rest-frame filter to retrieve the observed magnitude in the observed filter given the function:

$$\begin{array}{c}\hfill \begin{array}{c}K=-2.5log[\frac{1}{1+z}\frac{{\displaystyle \int {\lambda }_{\mathrm{obs}}}f({\lambda }_{\mathrm{obs}})R({\lambda }_{\mathrm{obs}})\mathrm{d}{\lambda }_{\mathrm{obs}}}{{\displaystyle \int {\lambda }_{\mathrm{obs}}}{g}^R({\lambda }_{\mathrm{obs}})R({\lambda }_{\mathrm{obs}})\mathrm{d}{\lambda }_{\mathrm{obs}}}\times \hfill \\ \frac{{\displaystyle \int {\lambda }_{\mathrm{rest}}}{g}^Q({\lambda }_{\mathrm{rest}})Q({\lambda }_{\mathrm{rest}})\mathrm{d}{\lambda }_{\mathrm{rest}}}{{\displaystyle \int {\lambda }_{\mathrm{rest}}}f((1+z){\lambda }_{\mathrm{rest}})Q({\lambda }_{\mathrm{rest}})\mathrm{d}{\lambda }_{\mathrm{rest}}}],\hfill \end{array}\end{array}$$(16)

with the observed wavelength λ_obs, observed flux f(λ_obs), transmission function of the filter in the observed frame R(λ_obs), rest-frame wavelength λ_rest, rest-frame flux f(λ_rest), and transmission function of the filter in the rest-frame Q(λ_rest). g^R(λ_obs) and g^Q(λ_rest) are the spectral flux densities for a magnitude of zero in the AB system in the observed frame and the rest-frame, respectively. We inserted the LSST u and g bands for the rest-frame transmission functions and the r and i bands for the observed-frame transmission function, as u and g are roughly shifted into these bands for z = 0.6. The resulting color curves are shown in Fig. 6.

Fig. 6. Rest-frame u − g color curve and corresponding redshifted r − i color curve created from model spectra. The kink of the rest-frame color curve is slightly weakened by the K-correction.

The strength of the kink was defined by the change in the slopes of the color curve before and after the kink. The stronger the change, the stronger the kink. Applying the K-correction, the kink of the color curve, which occurs in this case around day 95, is slightly weakened. As this change is not too drastic, we expect that redshifted SNe IIP color curves still exhibit a sufficiently strong kink for our method to work, especially when combining multiple color curves. To include redshifted color curves using the observed u and g bands, we will next investigate the evolution of rest-frame UV color curves.

5.2. Rest-frame UV color curves

To estimate how rest-frame UV color curves are redshifted into the range of LSST bands, we investigated the UV color curves with the UVOT filters uvw1 and b of the two SWIFT SNe SN2020jfo and SN2017eaw. We had to limit the SWIFT light curves to the UVOT rest-frame uvw1 and b bands, as only these two bands have good coverage of the epochs around the kink observed in the color curves. The data for SN2020jfo was taken from Sollerman et al. (2021) and for SN2017eaw from Szalai et al. (2019). The rest-frame uvw1 − b color curves of both SNe in Fig. 7 show a particularly strong kink. We found that this results from the long wavelength tail of the uvw1 filter reaching into the visible spectral range.

Fig. 7. uvw1 − b color curves of SWIFT SNe SN2017eaw (blue) (Szalai et al. 2019) and SN2020jfo (green) (Sollerman et al. 2021).

We investigated how cutting this uvw1 filter tail away at ∼3347.5 Å to create a filter uvw1^* affects the kink in color curves, by creating artificial color curves uvw1 − b and uvw1^* − b with the spectra by Dessart et al. (2013). The resulting color curves are shown in Fig. 8.

Fig. 8. Color curves uvw1 − b (blue) and uvw1* − b (orange) created with the Dessart et al. (2013) model spectra using the SWIFT filters uvw1, b, and the modified uvw1*, where the red tail has been cut. The color curve uvw1* − b still shows a sufficiently strong kink for our time delay method to work.

uvw1^* − b shows a weaker kink than uvw1 − b, but still stronger than the u − g rest-frame color curve from Fig. 6. We, therefore, conclude that bluer colors tend to show even stronger kinks for SNe IIP, which helps us constrain the time delay for redshifted SNe IIP.

6. Discussion and conclusion

We developed a Metropolis-Hastings MCMC-based method of fitting lensed SN IIP color curves with orthogonal polynomial functions to retrieve the time delay Δt and the differential dust extinction A_V, B − A_V, A between two lensing images A and B.

In the high-S/N limit, the method works flawlessly using either one or two color curves, and the time delay Δt and the differential dust extinction A_V, B − A_V, A can be retrieved precisely without bias. We note that A_V of each individual image is only determined precisely for the case of one color curve, assuming the knowledge of the macro-magnification. Orthogonal polynomial functions are flexible enough to fit data that may follow different functional forms.

Considering more realistic cases with low S/N and data loss because of observational constraints, the inference of Δt and especially with A_V, B − A_V, A is more uncertain. For one color curve, the precision of the time delay degrades to ≳3 days. Adding the information of a third light curve band helps constrain both parameters. Fitting two color curves calculated from light curves in 3 bands gives more constraints on A_V, B − A_V, A helping to break the degeneracy between the time delay and the differential dust extinction. This also has the advantage that absolute magnification cancels out between the two color curves, assuming we choose the correct R_V and extinction law during sampling.

We are aware that our assumption of knowing perfectly R_V and the correct extinction law is strong and simplifies the determination of Δt and A_V, B − A_V, A. We show in Appendix D that for combining two realistic mock color curves with low S/N and data loss using a different extinction law, Δt can still be determined precisely and without bias, which is the main objective of this work. The accuracy of A_V, B − A_V, A depends more on R_V than on the specific extinction law in the wavelength range of the LSST filters.

Despite all the measures we take to create realistic mock color curves, we anticipate some additional uncertainties in real data. Three additional nuisance effects may be expected for real data: the decrease in S/N because of the additional shot noise from the host galaxy of the SN and the lens, photometric errors due to the blending of the SN images, and the impact of the seeing on the above two effects. The influence of these effects will have to be studied in the future when more observed lensed SNe IIP are available. Another simplification of the existing framework is the use of a non-redshifted reddening. This should not impact our results heavily within the obtained uncertainties, as the observed bands probe a region of the extinction laws that is monotonically decreasing, and the lens will be likely at low redshifts around z = 0.2 − 0.4 (Oguri & Marshall 2010).

Considering all rest-frame cases presented in Sect. 4 as well as Appendix A, B, and C, we can determine several trends:

• The more color curves we take into account during sampling the better we can determine A_V, B − A_V, A, which also results in a better constraint on the time delay, as different color curves can help to break the degeneracy in Δt and A_V, B − A_V, A.

• The data needed for this method to work for rest-frame u, g, and r band should ideally be dedicated follow-up multiband observation with a cadence of not more than two rest-frame days and a depth that is 1 mag deeper than LSST to have a high enough S/N for typical lensed-SN redshifts and magnifications (Huber et al. 2022). In addition to dedicated follow-up observations, the upcoming space-based survey telescope Roman might also have a high enough cadence of ∼5 observed frame days, high resolution, and sensitivity. One should note it is covering the NIR starting roughly at the LSST r band wavelength range, limiting the amount of useful SNe IIP coverage to high-redshift ones as the color curves of redder rest-frame filters show almost no kink, especially considering the range of expected uncertainties (Pierel et al. 2021). Furthermore, we should remember that Roman is a survey that is more likely to cover a lensed SN by chance rather than collecting dedicated data through a follow-up target of opportunity-based observation.

• The K-correction investigation in Sect. 5 of rest-frame UV bands redshifted into the LSST bands u, g, r, and i, shows strong enough kinks to be usable for our MCMC retrieval technique. Bluer colors tend to show stronger kinks for SNe IIP. Hence, the redshift of a lensed SN IIP helps us to get better time-delay measurements with ground-based photometry.

• A higher S/N ratio, a larger number of considered color curves, and lower data loss result in a higher selected order of fitting for the orthogonal polynomials as the data impose higher constraints.

• Cluster-lensed SNe IIP are less susceptible to differential dust extinction, as dust occurs more in star formation regions, and the lensing images might not be in the line of sight of the lens. This makes cluster-lensed SNe IIP desirable for our time-delay retrieval method using color curves as they could provide more precise results compared to galaxy-scale lensed SNe IIP.

The method is precise enough for obtaining cosmology grade time delays for gravitationally lensed SNe IIP with realistic uncertainties around ∼1.0 day. For 20 SNe IIP with time delays longer than 20 days, we will thus be able to constrain H₀ with an uncertainty of ∼1% during the era of LSST (Suyu et al. 2020). Comparing the results of this work to our previous paper HOLISMOKES V (Bayer et al. 2021), which uses a combination of absorption lines within spectra of LSNe IIP, we achieve a similar accuracy and precision of ∼1.0 day.

Acknowledgments

We thank F. Dux, M. Millon, and F. Courbin for discussions on realistic data loss in observations, and J. Anderson and T. de Jaeger for providing us with the observational data of SN2005J and other SNe of the CSP for developing and testing our code, and E. Fitzpatrick for discussions on dust extinction laws. We also thank R. Cañameras for his feedback on our work and the anonymous referee for constructive comments. We acknowledge the use of public data from the Swift data archive. We thank the Max Planck Society for support through the Max Planck Research Group and Fellowship for SHS. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (LENSNOVA: grant agreement No. 771776; COSMICLENS: grant agreement No. 787886). This research is supported in part by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311.

References

Appendix A: Fitting with a different model

For our main study, we used the same functional form (orthogonal polynomials) to generate and fit the mock data. In this appendix, we explore if the method might produce model-induced biased results if the mock data is created with a spline model and sampled with orthogonal polynomials. Further, we do not expect real data to follow the behavior of orthogonal polynomials, nor a spline, and we want to ensure that the model used during sampling is flexible enough to fit the data based on a different model.

For the mock data creation, we created a spline function consisting of two third-order polynomials connected at one knot t_knot defined as:

$$\begin{array}{c}\hfill C(\mathit{\eta },t)=\{\begin{array}{c}{C}_1({\mathit{\eta }}_1,t)=c_1{t}^3+c_2{t}^2\hfill \\ +c_{3}t+c_4\text{for}t\le t_{\mathrm{knot}}\hfill \\ C_2(\mathit{\eta },t)=c_4+c_5{t}^3\hfill \\ +c_6{t}^2+c_{7}t\hfill \\ +(c_1-c_5)t_{\mathrm{knot}}^3\hfill \\ +(c_2-c_6)t_{\mathrm{knot}}^2\hfill \\ +(c_3-c_7)t_{\mathrm{knot}}\text{for}t>t_{\mathrm{knot}},\hfill \end{array}\end{array}$$(A.1)

with the additional condition of the same slope at the knot:

$$\begin{array}{c}\hfill \frac{\partial C_1({\mathit{\eta }}_1,t_{\mathrm{knot}})}{\partial t}=\frac{\partial C_2(\mathit{\eta },t_{\mathrm{knot}})}{\partial t}.\end{array}$$(A.2)

The collective parameter vectors were defined as $\mathit{\eta }=(c_1,c_2,c_3,c_4,c_5,c_6,c_7,t_{\mathrm{knot}})$ and ${\mathit{\eta }}_1=(c_1,c_2,c_3,c_4)$ for C₁(η₁, t). With the slope condition, we can eliminate the parameter c₇ and express it via the other parameters:

$$\begin{array}{c}\hfill c_7=3t_{\mathrm{knot}}^2(c_1-c_5)+2t_{\mathrm{knot}}(c_2-c_6)+c_3.\end{array}$$(A.3)

We inserted Eq. A.3 into Eq. A.1:

$$\begin{array}{c}\hfill C(\mathit{\eta },t)=\{\begin{array}{c}{C}_1({\mathit{\eta }}_1,t)=c_1{t}^3+c_2{t}^2\hfill \\ +c_{3}t+c_4\text{for}t\le t_{\mathrm{knot}}\hfill \\ C_2(\mathit{\eta },t)=c_5{t}^3+c_6{t}^2\hfill \\ +[3t_{\mathrm{knot}}^2(c_1-c_5)\hfill \\ +2t_{\mathrm{knot}}(c_2-c_6)\hfill \\ +c_3]t+(c_1-c_5)t_{\mathrm{knot}}^3\hfill \\ +(c_2-c_6)t_{\mathrm{knot}}^2\hfill \\ +(c_3-[3t_{\mathrm{knot}}^2(c_1-c_5)\hfill \\ +2t_{\mathrm{knot}}(c_2-c_6)\hfill \\ +c_3])t_{\mathrm{knot}}+c_4\text{for}t>t_{\mathrm{knot}},\hfill \end{array}\end{array}$$(A.4)

and redefined the parameter vector $\mathit{\eta }=(c_1,c_2,c_3,c_4,c_5,c_6,t_{\mathrm{knot}})$. Besides the change in generating mock data, we followed the same approach as described in Sect. 3.

We first investigated high S/N with only one color curve u − g and 100 noise realizations, as this simple case will show most easily if a bias occurs. The resulting corner plot is shown in Fig. A.1 for which only the 6th order was selected by the BIC. It has a slight bias in the time shift and the differential extinction, although the bias in Δt is only 0.03 days, which is at the subpercent level for most lensed SN systems of interest in our program.

Fig. A.1. Corner plot for mock data created with spline fitting, including 100 noise realizations with high S/N. The initial input values shown in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag.

We further checked the case of low S/N with one color curve u − g and 500 noise realizations. Because of the low S/N, the result shown in Fig. A.2 is unbiased within the uncertainties. The selected orders are shown in the bar plot in A.3

Fig. A.2. Corner plot for mock data created with spline fitting including 500 noise realizations with low S/N. The initial input values shown in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag.

This proves that for the realistic cases with low S/N, we can retrieve unbiased results even with different models used in mock data creation and sampling.

Fig. A.3. Histogram of the selected orthogonal polynomial orders of fitting using the BIC for the low S/N mock data created with spline fitting.

Appendix B: Completely degenerate case

Another test of the functionality of our code was to create a mock color curve case that was perfectly linear in both images and tried to retrieve a time delay with our method. We expected that the preferred order selected with the BIC of the orthogonal polynomials would be one, and the distribution in the time shift and differential extinction would be completely degenerate. As it can be seen in Fig. B.1, we can retrieve an exact degeneracy up to the borders of the prior, and during the BIC selection, only order one was selected. Even if A_V, B − A_V, A has no prior boundaries as it has an unbound uniform prior, it is limited by the boundaries of the prior for Δt. This test shows that for featureless (linear) color curves, our color curve time-delay determination breaks down.

Fig. B.1. Corner plot combining 100 noise realizations of perfect linear mock color curves u − g. Because of the linearity of the mock data, a complete degeneracy is retrieved.

Appendix C: Different shifts

Our main study focused only on one time shift Δt = 16.3 days. Here, we also investigated two cases with Δt_small = 0.5 days and Δt_large = 50 days to check that any low or high time shift can be retrieved. We did not explore time shifts longer than 50 days, as these can also be detected by eye and shifted manually into the range of Δt< 50 days, which ensures that our code can be applied to any time shift. To create the two plots in Fig. C.1 and Fig. C.2, we followed the same procedure as explained in Sect. 3 with 100 noise realizations and high S/N data, just switching out the time shift in the mock data creation to Δt_small = 0.5 days and Δt_large = 50 days. For both new time shifts, a similar distribution of the fitting orders as for case 1 with one color curve shown in Fig. 3b was selected.

Fig. C.1. Corner plot of 100 noise realizations of high S/N mock data with Δtsmall = 0.5 days and AV, B − AV, A = 0.2 mag

Fig. C.2. Corner plot of 100 noise realizations of high S/N mock data with Δtlarge = 50.0 days and AV, B − AV, A = 0.2 mag

For both cases, we can retrieve the time delay and differential dust extinction, as well as for case 1 considered in Sect. 4 using only one color curve u − g, which confirms that the code can be applied to any time shift within 0.5 days to 50 days.

Besides all the testing presented here in the appendix, we also used the nested sampling Monte Carlo algorithm MLFriends (Buchner 2014, 2019) from the UltraNest² package (Buchner 2021) with our technique and obtained the same results, confirming that our sampling code works correctly. Using our own code gave us even faster results because we did not calculate the Bayesian evidence.

Appendix D: Different extinction laws for mock data and sampling

In the main study, we knew and adopted the correct extinction law to sample the color curves. In reality, this will not be the case, and we might assume R_V wrongly or use an extinction law not fitting the real data. Therefore, we tested three additional cases with 100 noise realizations for the realistic case of low S/N with 35% data loss for two color curves u − g and u − r. All three cases were still based on the same mock data as explained in Sect. 2.4.

The first test case, T1, set R_V= 2.0 for the sampling. The second case, T2, used R_V= 5.0 for sampling. For both T1 and T2, we still used Cardelli extinction laws in the sampling process described in Sect. 3. The third test case, T3, applied more recent extinction laws from Gordon et al. (2023), referred to as Gordon laws, in the sampling, but kept R_V= 3.1. The Gordon laws have a similar approach as the Cardelli laws in Eq. 4, but are based on more recent measurements for a(x) and b(x) and are normalized differently:

$$\begin{array}{c}\hfill \frac{A_{\lambda }}{A_V}=a(x)+\frac{b(x)}{R_V}-\frac{b(x)}{3.1}.\end{array}$$(D.1)

The resulting corner plots for T1, T2, and T3 are shown in Fig. D.1, D.2, and D.3 respectively. The orthogonal polynomial orders selected via the BIC for all cases follow a similar distribution, as shown in Fig. 5d. For T1 and T2, the time delays have a bias of 0.1-0.2 days, which for delays > 20 days are still at the < 1% level. That is within the tolerance to achieve a 1% H₀ measurement. For short delays, the R_V value has a relatively larger impact on the accuracy of the delays. Therefore, R_V should be incorporated in the measurement of the delays. In both test cases of T1 and T2, A_V, B − A_V, A can not be sampled correctly anymore, which is to be expected considering an incorrect R_V value.

Fig. D.1. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves setting RV = 2.0 during MCMC sampling.

Fig. D.2. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves setting RV = 5.0 during MCMC sampling.

For T3, both Δt and A_V, B − A_V, A are retrieved unbiased and precisely because the Cardelli extinction laws and the Gordon laws follow a similar behavior in the considered wavelength range and differ only a bit in the a(x) and b(x) values. If the data covered the ultraviolet or infrared wavelength range of the extinction laws, the results would deviate more. The same would occur if extinction laws that do not follow the Cardelli laws as closely as the Gordon laws would be applied.

Fig. D.3. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves using extinction laws by Gordon et al. (2023) with RV = 3.1 during MCMC sampling.

Appendix E: Figures and table for the BIC selection

Fig. E.1. Orthogonal polynomial fits to one noise realization of case 1 to demonstrate the BIC selection. Panel a: 6th order, Panel b: ${\displaystyle \underline{5^{\mathrm{th}}\text{order}}}$, Panel c: 4th order, Panel d: 3rd order, Panel e: 2nd order; The selected order is underlined and highlighted with a red frame of the plot.

Fig. E.2. Orthogonal polynomial fits to one noise realization of case 3 to demonstrate the BIC selection. For the 3rd order, the retrieved time delay and differential dust extinction deviate strongly from the input, which is visible from the different shapes of the fits in the range of the data points of image A and image B. Panel a: 6th order, Panel b: 5th order, Panel c: ${\displaystyle \underline{4^{\mathrm{th}}\text{order}}}$, Panel d: 3rd order, Panel e: 2nd order; The selected order is underlined and highlighted with a red frame of the plot.

All Tables

All Figures

	Fig. 1. Exemplary noiseless LSN IIP color curves of two lensing images A and B showing a time shift Δt and a shift in color s due to extinction.
In the text

	Fig. 2. u,  g, and r light curves of SN2005J observed by the CSP (Anderson et al. 2024) including examples of orthogonal polynomial fits (solid lines) to the u, g and r band light curves.
In the text

Fig. 3. Corner plots of the MCMC chain results and bar plots of selected orders for 100 noise realizations for case 1. The left panels show corner plots of the time shift and differential dust extinction retrieved from 100 noise realizations for case 1. We plot the contours for the 1σ (68%, solid red line) and 2σ (95%, dashed red line) credible intervals for the differential dust extinction AV, B − AV, A and the time delay Δt. The input values in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag. On the right, we show histograms of the selected orthogonal polynomial orders of fitting using the BIC for the corresponding cases on the left. Panel a: Case 1 with one color curve u − g. Panel b: Orders of fitting selected for case 1 with one color curve u − g. Panel c: Case 1 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 1 with two color curves u − g and u − r. We impose the same order for both color curves.

In the text

	Fig. 4. Same as Fig. 3, but for 500 noise realizations for case 2. Panel a: Case 2 with one color curve u − g. Panel b: Orders of fitting selected for case 2 with one color curve u − g. Panel c: Case 2 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 2 with two color curves u − g and u − r. We impose the same order for both color curves.
In the text

	Fig. 5. Same as Fig. 3, but for 500 noise realizations for case 3. Panel a: Case 3 with one color curve u − g. Panel b: Orders of fitting selected for case 3 with one color curve u − g. Panel c: Case 3 with two color curves u − g and u − r. Panel d: Orders of fitting selected for case 3 with two color curves u − g and u − r. We impose the same order for both color curves.
In the text

	Fig. 6. Rest-frame u − g color curve and corresponding redshifted r − i color curve created from model spectra. The kink of the rest-frame color curve is slightly weakened by the K-correction.
In the text

	Fig. 7. uvw1 − b color curves of SWIFT SNe SN2017eaw (blue) (Szalai et al. 2019) and SN2020jfo (green) (Sollerman et al. 2021).
In the text

	Fig. 8. Color curves uvw1 − b (blue) and uvw1* − b (orange) created with the Dessart et al. (2013) model spectra using the SWIFT filters uvw1, b, and the modified uvw1*, where the red tail has been cut. The color curve uvw1* − b still shows a sufficiently strong kink for our time delay method to work.
In the text

	Fig. A.1. Corner plot for mock data created with spline fitting, including 100 noise realizations with high S/N. The initial input values shown in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag.
In the text

	Fig. A.2. Corner plot for mock data created with spline fitting including 500 noise realizations with low S/N. The initial input values shown in light blue are Δt = 16.3 days and AV, B − AV, A = 0.2 mag.
In the text

	Fig. A.3. Histogram of the selected orthogonal polynomial orders of fitting using the BIC for the low S/N mock data created with spline fitting.
In the text

	Fig. B.1. Corner plot combining 100 noise realizations of perfect linear mock color curves u − g. Because of the linearity of the mock data, a complete degeneracy is retrieved.
In the text

	Fig. C.1. Corner plot of 100 noise realizations of high S/N mock data with Δtsmall = 0.5 days and AV, B − AV, A = 0.2 mag
In the text

	Fig. C.2. Corner plot of 100 noise realizations of high S/N mock data with Δtlarge = 50.0 days and AV, B − AV, A = 0.2 mag
In the text

	Fig. D.1. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves setting RV = 2.0 during MCMC sampling.
In the text

	Fig. D.2. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves setting RV = 5.0 during MCMC sampling.
In the text

	Fig. D.3. Corner plot of 100 noise realizations of low S/N mock data with 35% data loss and with input time delay Δt = 16.3 days for combining u − g and u − r color curves using extinction laws by Gordon et al. (2023) with RV = 3.1 during MCMC sampling.
In the text

	Fig. E.1. Orthogonal polynomial fits to one noise realization of case 1 to demonstrate the BIC selection. Panel a: 6th order, Panel b: ${\displaystyle \underline{5^{\mathrm{th}}\text{order}}}$, Panel c: 4th order, Panel d: 3rd order, Panel e: 2nd order; The selected order is underlined and highlighted with a red frame of the plot.
In the text

Fig. E.2. Orthogonal polynomial fits to one noise realization of case 3 to demonstrate the BIC selection. For the 3rd order, the retrieved time delay and differential dust extinction deviate strongly from the input, which is visible from the different shapes of the fits in the range of the data points of image A and image B. Panel a: 6th order, Panel b: 5th order, Panel c: ${\displaystyle \underline{4^{\mathrm{th}}\text{order}}}$, Panel d: 3rd order, Panel e: 2nd order; The selected order is underlined and highlighted with a red frame of the plot.

In the text