SN H0pe: The First Measurement of $H_{0}$ from a Multiply-Imaged Type Ia Supernova, Discovered by JWST

Abstract

The first James Webb Space Telescope (JWST) Near InfraRed Camera (NIRCam) imaging in the field of the galaxy cluster PLCK G165.7+67.0 ($z=0.35$) uncovered a Type Ia supernova (SN Ia) at $z=1.78$, called “SN H0pe.” Three different images of this one SN were detected as a result of strong gravitational lensing, each one traversing a different path in spacetime, thereby inducing a relative delay in the arrival of each image. Follow-up JWST observations of all three SN images enabled photometric and rare spectroscopic measurements of the two relative time delays. Following strict blinding protocols which oversaw a live unblinding and regulated post-unblinding changes, these two measured time delays were compared to the predictions of seven independently constructed cluster lens models to measure a value for the Hubble constant, $H_{0}=$ $71.8^{+9.8}_{-7.6}$ km s^-1 Mpc^-1. The range of admissible $H_{0}$ values predicted across the lens models limits further precision, reflecting the well-known degeneracies between lens model constraints and time delays. It has long been theorized that a way forward is to leverage a standard candle, however this has not been realized until now. For the first time, the lens models are evaluated by their agreement with the SN absolute magnification, breaking these degeneracies and producing our best estimate, $H_{0}=$ $75.4^{+8.1}_{-5.5}$ km s^-1 Mpc^-1 . This is the first precision measurement of $H_{0}$ from a multiply-imaged SN Ia, and provides a measurement in a rarely utilized redshift regime. This result agrees with other local universe measurements, yet exceeds the value of $H_{0}$ derived from the early Universe with $\gtrsim 90\%$ confidence, increasing evidence of the Hubble tension. With the precision provided by only four more events, this approach could solidify this disagreement to $>3\sigma$.

1 Introduction

Determination of the current expansion rate of the Universe (the Hubble Constant, $H_{0}$) across a wide range of redshifts provides a fundamental test of the standard cosmological model, $\Lambda$CDM. While this model succeeds in reproducing an abundance of observed cosmological phenomena, disagreement between independent measurements of $H_{0}$ at early and late times in the Universe raises doubts about the model’s reliability. This tension primarily results from two independent precision methods for measuring $H_{0}$: SNe Ia using the local-distance-ladder approach and CMB observation in conjunction with the $\Lambda$CDM model. Confirming or resolving the $H_{0}$ discrepancy is crucial for fundamental physics including neutrinos’ effects, the nature of dark energy, and the spatial curvature of the Universe.

For the local distance-ladder, the ‘Supernova $H_{0}$ for the Equation of State’ (SH0ES) team’s (Riess et al., 2022) imaging of Cepheid variable stars with the Hubble Space Telescope (HST) yields $H_{0}=73.04\pm 1.04$ km s^-1 Mpc^-1. In contrast, CMB measurements (Planck Collaboration et al., 2020) from the Planck satellite favor $H_{0}=67.4\pm 0.6$ km s^-1 Mpc^-1. This apparent ${>}5\sigma$ disagreement between two premier methods is the basis of what is known as the ‘Hubble tension’. Examining other independent $H_{0}$ measurements can delineate the underlying discrepancy. Local-distance-ladder measurements using HST imaging of the ‘tip of the red-giant branch’ (TRGB) stars to calibrate SNe Ia, rather than Cepheids, from the Carnegie-Chicago Hubble Program (CCHP, Freedman et al., 2019), the Extragalactic Distance Database (EDD, Anand et al., 2022), and Comparative Analysis of TRGBs (CATs, Scolnic et al., 2023) yield a range in $H_{0}$ of 70–73 km s^-1 Mpc^-1. Through comparison of measured colors to measured surface-brightness fluctuations (SBF) calibrated on TRGB and Cepheids, Blakeslee et al. (2021) found $H_{0}=73.3\pm 0.7\pm 2.4$ km s^-1 Mpc^-1, while Garnavich et al. (2023) found $74.6\pm 0.9\pm 2.7$ km s^-1 Mpc^-1 by adding SNe Ia as the top rung. In a ‘single-rung’ distance ladder, the Maser Cosmology Project (MCP, Pesce et al., 2020) combined six maser distances and peculiar velocities to find $H_{0}=73.9\pm 3.0$ km s^-1 Mpc^-1. Another single-rung method uses time variability between the multiple images of strongly gravitationally lensed quasars. The Time Delay COSMOgraphy (TDCOSMO) collaboration leveraged a joint measurement across seven lensed quasars to find $H_{0}=74.2^{+1.6}_{-1.6}$ km s^-1 Mpc^-1 assuming power-law mass profiles for the lens and a flat $\Lambda$CDM cosmology, and $H_{0}=73.3^{+5.8}_{-5.8}$ km s^-1 Mpc^-1 using free-form mass profiles (Birrer et al. 2020; see also Millon et al. 2020, Wong et al. 2020, and Shajib et al. 2023). This growing catalog of independent measurements has already provided insights into the Hubble tension (Verde et al., 2023), and new measurements at intermediate redshifts $1<z<10$ are imperative in elucidating when in cosmic history this discrepancy occurs.

Refsdal (1964) first recognized that strong gravitational lensing of SNe enables another independent method for measuring $H_{0}$. The ‘strong’ regime of gravitational lensing refers to the deflection of the light of a distant background object (star, galaxy, or QSO) by a foreground galaxy or cluster of galaxies such that the image of a single object is observed at multiple positions on the sky. For each observed image of an “image system,” the background object’s light traverses a different path length, resulting in different travel times for photons to each image. Transient events within the lensed object (e.g., a supernova or quasar variability) enable measurement of the time delay between images. The time delay, together with a model of the lensing mass distribution, yields the geometric distance to the source object, a quantity inversely proportional to $H_{0}$.

Although initially envisioned for SNe, inference of $H_{0}$ through strong-lensing time delays was first applied to multiply-imaged quasars lensed by single foreground galaxies (Suyu et al., 2010). The technique, known as time-delay cosmography, was not realized for SNe until the discovery of ‘SN Refsdal’, a multiply-imaged Type II SN found in the MACS J1149.6+2223 cluster field (Kelly et al., 2015). Crucially, strongly-lensed SNe in galaxy-cluster lenses provide an independent route for measuring $H_{0}$ whose systematics are distinct from the distance-ladder, quasar-time-delay cosmography, and CMB approaches (Suyu et al., 2023). SN Refsdal was observed as four images in an Einstein cross around a single elliptical-galaxy lens (a ‘galaxy–galaxy’ lens). A fifth image with an $\sim$1 year relative time delay was predicted (Sharon & Johnson, 2015; Diego et al., 2016; Oguri, 2015) and later observed (Kelly et al., 2016). The relatively long time delay enabled SN Refsdal to be the first lensed SN for which meaningful constraints of $H_{0}$ could be inferred, and the nature of the galaxy-cluster-scale lens gave these constraints distinct systematics from previous galaxy–galaxy lens methods (Grillo et al., 2018; Vega-Ferrero et al., 2018). The most robust constraints to date from SN Refsdal were provided by Kelly et al. (2023), who inferred $H_{0}=64.8^{+4.4}_{-4.3}$ km s^-1 Mpc^-1 or $H_{0}=66.6^{+4.1}_{-3.3}$ km s^-1 Mpc^-1 depending on the set of lens models used.

On the heels of the discovery of SN Refsdal, multiple other lensed supernova have been observed. These include other Type II SN, such as the $z\sim 3$ SN in the Abell 370 cluster field which allowed detailed studied of shock cooling in a red supergiant SN (Chen et al., 2022), as well as SN Ia. ‘iPTF16geu’ at $z=0.409$ (Goobar et al., 2017) was the first SN Ia to have multiple images observed. It was followed by ‘SN Requiem’ at $z=1.95$ (Rodney et al., 2021), and ‘SN Zwicky’ at $z=0.3554$ (Goobar et al., 2023; Pierel et al., 2023). Both iPTF16geu and Zwicky were galaxy-scale lenses, whose $\sim$day-long time delays limited inference of $H_{0}$ to $\gtrsim$40% precision (Dhawan et al., 2020; Pierel et al., 2023). By contrast, the galaxy-cluster-scale lens of SN Requiem has a nearly decade-long time delay, which will yield a precision measurement of $H_{0}$ in 2037 when the SN counterimage is predicted to appear (Rodney et al., 2021). Across all of these SNe, SN Refsdal is the only to provide precision $H_{0}$ constraints to date.

SN H0pe is a triply-imaged SN Ia discovered in the galaxy-cluster field PLCK G165.7+67.0 (‘G165’, $z=0.348$). It is a massive lensing galaxy cluster ($M_{\rm tot}=(2.6\pm 0.3)\times 10^{14}$ M_⊙, Frye et al., 2019; Pascale et al., 2022a), that induces a relatively long $\sim$100 day time delay between the first and third image, offering only the second-ever opportunity for precision SN time-delay cosmography and the first with a SN Ia. The three supernova images were discovered in JWST/NIRCam imaging in the lensed, NIR-bright host galaxy ‘Arc 2’, which itself is triply imaged (Frye et al., 2023). Polletta et al. (2023) first measured the Arc 2 redshift as $z=1.783\pm 0.002$. The SN H0pe discovery paper (Frye et al., 2024) used followup JWST/NIRSpec spectroscopy to identify SN H0pe as a SN Ia, to confirm the SN host galaxy redshift to higher precision ($z=1.7834\pm 0.0005$), and to spectroscopically-confirm a coincident galaxy overdensity at the redshift of SN H0pe. The known spectroscopic evolution of SNe Ia allowed Chen et al. (2024) to determine the SN phases in the three images of SN H0pe and thereby measure relative time delays of $-122.3^{+43.7}_{-43.8}$ days and $-49.3^{+12.2}_{-14.7}$ days between images a–b and c–b respectively. Image b is the last to arrive and consequently shows the earliest phase of the SN, which was near peak brightness at the time of discovery. This is the first time spectroscopy has provided useful time-delay constraints. Chen et al. (2024) also classified SN H0pe as a normal SN Ia. Pierel et al. (2024) used three epochs of JWST NIRCam six-band imaging to photometrically measure relative image time delays of $-116^{+10.8}_{-9.3}$ days and $-48.3^{+3.6}_{-4.0}$ days for images a–b and c–b respectively.

Independently through the photometric and spectroscopic datasets respectively, Pierel et al. (2024) and Chen et al. (2024) use the standard-candle nature of the SN to constrain the absolute magnifications induced by strong lensing. Absolute magnification measurements have long been thought to break degeneracies among mass models inherent in $H_{0}$ measurements from strong-lens time delays (Falco et al., 1985; Kolatt & Bartelmann, 1998; Oguri & Kawano, 2003; Schneider & Sluse, 2013), yet have previously been unavailable for time delay cosmography. For example in the case of the cluster-lensed SN Refsdal, the inferred $H_{0}$ was found to vary across twenty-three different lens mass models despite all satisfying the same lensing constraints. However, since $H_{0}$ was found to strongly correlate with the predicted magnifications, magnification measurements could have broken this degeneracy (Liu & Oguri, 2024). SN H0pe provides the first opportunity to harness the standard candle for time delay cosmography.

This paper presents the first precision measurement of $H_{0}$ from time-delay cosmography of a strongly-lensed SN Ia. The outline is as follows. Section 2 summarizes the underlying theory of time-delay cosmography. Section 3 describes the HST and JWST observations that led to SN H0pe’s discovery and subsequent analysis. Sections 4 and 5 briefly summarize the photometric light-curve fitting and spectroscopic age-dating analyses. Section 6 details the construction of the seven contributing lens models. Section 7 outlines the process through which $H_{0}$ was derived from measured observables. Section 8 outlines the guidelines for the blinded inference. Section 9 presents the results of the $H_{0}$ inference. Section 10 discusses the unique advantages provided by the SN Ia, the robustness of the various inference approaches, and the overall error budget. Finally, section 11 discusses future work and how SN H0pe fits into the broader picture of time-delay cosmography.

Throughout this analysis, we enforced a blinding between the seven independent lens models, the photometrically-derived time delays, and the spectroscopically-derived time delays from each other. Changes made after the unblinding and the impacts of each change are explicitly discussed in Appendix C.

2 Time Delay Cosmography

Here we briefly outline the how time delay, $\Delta t_{i,j}$, between any two lensed images $i$ and $j$ of a multiply-imaged source depends on the positions of the observed images with respect to the source, the local lensing potential, and the assumed cosmological model. For an individual source at angular position $\beta$ with a corresponding lensed image observed at angular position $\theta$, the time delay can be expressed as:

$$\Delta t(\theta)=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}[\frac{1}{2}(\theta-\beta)^{2}-\psi(\theta)]\quad,$$ (1)

where $z_{l}$ is the redshift of the lensing cluster, $\psi(\theta)$ is the lensing potential at the observed image position $\theta$, and $D_{l}$, $D_{s}$, and $D_{ls}$ are the angular diameter distances to the lens, the source, and between the lens and source respectively. If instead we measure the time delay between any two images of a set of multiple images:

$$\Delta t_{i,j}(\theta_{i},\theta_{j})=D_{\Delta t}\Delta\tau_{i,j}(\theta_{i},\theta_{j},\beta)\quad,$$ (2)

where $\tau$ is the Fermat potential, $\tau(\theta,\beta)=\frac{1}{2}(\theta-\beta)^{2}-\psi(\theta)$ (Schneider, 1985; Blandford & Narayan, 1986). $D_{\Delta t}$ is the so-called ‘time delay distance’, $D_{\Delta t}=\frac{1+z_{l}}{c}\frac{D_{l}D_{s}}{D_{ls}}$, and is dependent on the angular diameter distances and hence the cosmological model (Refsdal, 1964; Schneider et al., 1992; Suyu et al., 2010). $D_{\Delta t}$ is inversely proportional to $H_{0}$, while being only weakly dependent on other cosmological parameters (Bonvin et al., 2017). Hence, if the time delay between any two images is known, as well as the redshifts of the lens and source, then a model of the lensing potential is necessary to constrain $H_{0}$.

In practice a lens model can be constructed from lensing observables, such as the observed positions of multiply-imaged galaxies or information of the foreground lensing galaxies, and assumes a fiducial value of $H_{0}$, $H_{0}^{\rm fid}=70$ km s^-1 Mpc^-1, to predict a fiducial time delay $\Delta t_{i,j}^{\rm fid}$. Following equation 1, the true measured time delay for any image pair can be related to this lens model fiducial time delay by rescaling $H_{0}$, in turn inferring $H_{0}$:

$$\Delta t_{i,j}^{\rm pred}(H_{0})=\Delta t_{i,j}^{\rm fid}\times\frac{70~{}{\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}}{H_{0}}\quad.$$ (3)

Galaxy-cluster lenses can present multiple sets of image systems corresponding to sources at different redshifts, allowing for precision constraints on the lens model. Crucially, the availability of spectroscopic redshifts for two or more image systems at different redshifts helps to break the mass-sheet degeneracy, a necessity for measurement of $H_{0}$ through time delays (Falco et al., 1985; Suyu et al., 2020; Schneider & Sluse, 2014). With additional image systems, the cluster potential can be constrained to high accuracy (Johnson & Sharon, 2016). In the case of G165, five systems at three unique redshifts are confirmed spectroscopically in addition to 16 image systems with photometric constraints, breaking the mass-sheet degeneracy and making the cluster well-suited for precision lens modeling.

3 Observations

The discovery paper for SN H0pe (Frye et al., 2024) supplied details regarding the JWST imaging and spectroscopic data sets. In brief, SN H0pe was initially discovered in NIRCam images from the PEARLS GTO program (PI: Windhorst, PID: 1176; Windhorst et al., 2023). The discovery images were taken on 31 March 2023 in eight filters. The SN was recognized by comparing the F150W PEARLS NIRCam imaging to WFC3/IR F160W HST images taken on 30 May 2016 (Figure 1). The discovery prompted a JWST disruptive Director’s Discretionary Time (DDT) program which acquired two additional epochs of NIRCam imaging in six filters (PI: Frye, PID: 4446) on 22 April 2023 and 9 May 2023. Photometric analysis of the galaxy cluster field was presented by Frye et al. (2024), and photometric light curve fitting of the SN was presented by Pierel et al. (2024).

The DDT program also acquired JWST/NIRSpec Mico-shutter assembly (MSA) spectra (Frye et al., 2024; Chen et al., 2024) in the PRISM, G140M, and G235M gratings targeting the SN host galaxy and the two brightest SN images. The spectra, taken on 22 April, confirmed the Ia type from the identification of the requisite blueshifted [Si II]$\lambda$6355 feature, Ca H&K, and other absorption-line features commonly found in Type Ia SNe (Figure 2).

Figure 1: JWST/NIRCam color image in the central region of G165. Insets show closeup of the boxed region depicting the three images of the host galaxy Arc 2 prior to SN H0pe in HST WFC3/IR F160W imaging from 2016 (lower left) and during its appearances in JWST/NIRCam F150W imaging from 2023 (upper right). The three images of the Type Ia SN H0pe appear only in the 2023 images, and their positions are marked by the blue boxes, where the ‘A’, ‘B’, and ‘C’ refer to images $a$, $b$, and $c$. North is up and east is to the left.

Figure 2: High signal-to-noise NIRSpec spectrum of SN H0pe. The G140M (green) and G235M (blue) grating spectra are shown. Some prominent line features are identified, such as the [Si II] line blueshifted to rest-frame $\sim 6150$ Å that is a requiste line feature of a Type Ia supernova. A couple of features are also traced to the SN host galaxy, such as NaD, H$\alpha$, [N II]$\lambda\lambda$6548,6584 (magenta dotted-line). The observed wavelength is given on the lower abscissa and the rest wavelength in the frame of the host galaxy is provided on the upper abscissa.

4 Photometric Time Delays and Magnifications

The light curve evolution of SNe Ia is well-studied (e.g., Hsiao et al., 2007; Pierel et al., 2021), such that precision measurements of the time delay can be made by fitting measured SN colors across multi-epoch imaging. Given that there are three epochs of NIRCam imaging, and each epoch contains three SN images, the SN evolution is photometrically sampled at nine points in time. Here we briefly summarize the light curve fitting analysis of Pierel et al. (2024). The halo of the bright galaxy host, a significant source of noise, is reconstructed and subtracted through simultaneous fitting of its surface brightness profile and cluster lens model with the GLEE software package (Suyu et al., 2010; Suyu, 2012, details in Caminha et al., in prep.). The photometry of the SN images is then fit simultaneously to BayeSN SN Ia models (Mandel et al., 2022; Grayling et al., 2024). Time delays are inferred by fitting a common set of light curve parameters to each of the SN images. The magnification ratios are estimated by normalizing the light curve to image $b$, and then computing the ratios of images $a$ and $c$ relative to $b$. Finally, the absolute magnifications are estimated by comparing the peak F277W flux ($\sim$rest-frame Y band) predicted by the BayeSN models, to the population of field SNe Ia.

Statistical uncertainties of the inferred time delays and magnifications are estimated by bootstrapping the photometric uncertainties assuming Gaussian errors. Systematic uncertainties resulting from the BayeSN model itself are produced in a similar manner. Systematic uncertainties from chromatic microlensing are estimated by examining the range of fits to a set of 1000 simulated SN H0pe light curves with variable realizations of the foreground microlensing. Systematics from millilensing by dark matter subhalos in the foreground cluster lens are achromatic, and applied directly as a correction to the distribution of inferred magnifications. Finally, the $\sim 0.1$ mag intrinsic scatter of SN Ia absolute magnitudes (Scolnic et al., 2018) is applied to the magnification distribution.

5 Spectroscopic time delays and Magnifications

Spectra from multiple images of a strongly-lensed SN obtained during a single observing epoch enable the measurement of time delays via spectroscopic-based phase analysis. Previously, the lack of quality simultaneously-obtained spectra have precluded such an analysis. Since SN H0pe was discovered from the PEARLS NIRCam Epoch 1 observation when its multiple images were near peak brightnesses, the rapid acquisition of spectra in a “disruptive” JWST DDT program using the JWST NIRSpec’s MSA enabled the identification of the SN as Type Ia through the detection of the requisite blueshifted [Si II]$\lambda$6155 absorption-line feature.

The identification of SN H0pe as a SN Ia allowed for more straightforward template fitting to determine the phases between images using spectroscopic methods. Below, we describe the process of measuring time delays through spectroscopic age-dating from NIRSpec data. We refer to Chen et al. (2024) for details.

Extraction of the SN spectrum from the NIRSpec data is performed using custom software designed to cope with background subtraction, host galaxy contamination, and bad pixels (e.g., cosmic rays). The error contribution from these potential systematics is included in the overall uncertainties, along with other sources of error such as slit loss and flux calibration. The relative delays between the images’ emitting times are estimated by determining their SN phases (days post-B-band peak brightness) through a joint MCMC fit of the SN spectra against various spectroscopic templates (e.g. Hsiao et al., 2007). The magnification is estimated by comparing the rest-frame $Y$-band magnitude of the best fit spectral models with that of an unmagnified SN Ia at $z=1.783$. Systematics from microlensing, millilensing, and the $\sim 0.1$ mag intrinsic scatter in SN Ia absolute magnitudes are addressed precisely as described in the photometric light curve fitting section outlined previously.

6 Lens Modeling

To obtain the cluster lens model predicted time delays and magnifications, seven different teams agreed to construct independent lens models for this analysis. All lens modelers were supplied with the same input constraints. These consist of 21 image systems, five of them with spectroscopic redshifts, and 161 cluster member galaxies with the same positions, F200W brightnesses, and morphological parameters. The image system and cluster member catalogs can be found in Frye et al. (2024). We emphasize that each lens model was created without knowledge of the photometric or spectroscopic time delay and magnification measurements, and the results of any given model are also blinded from the other modeling teams. All lens models assume a fiducial flat cosmology $\Omega_{m}=0.3$, $\Omega_{\lambda}=0.7$, and $H_{0}^{\rm fid}=70$ km s^-1 Mpc^-1 . The time delays and magnification predictions from the lens models are measured at the lens model-predicted positions of the images rather than the observed positions.

We briefly summarize the general methods of lens models included in this analysis, the assumptions made by each, and the expected systematics. A full discussion of each lens model can be found in Appendix A. The seven lens models in this analysis are broadly split into two categories. Parametric models (models 1, 2, 3) assume the mass distribution is entirely characterized as the superposition of analytic profiles. These models typically impose large-scale profiles to describe the cluster-scale dark matter halo, and smaller halos to describe each galaxy, which include priors based on the galaxy brightness or morphology (Zitrin et al., 2015; Oguri, 2010; Jullo et al., 2007). Non-parametric models (models 4 and 6), on the other hand, make fewer assumptions on the profile of the cluster-scale dark matter, and allow a flexible grid to describe the mass distribution (Diego et al., 2005; Sendra et al., 2014). Some models (models 5 and 7) are so-called ‘semi-parametric’, and may make use of both analytic profiles and more flexible non-parametric distributions. Parametric models benefit from a small number of free parameters and mass distributions which more closely mirror the galaxy configuration of the cluster, while non-parametric methods are given greater freedom in the shape of the underlying mass distribution and include a number of regularization parameters and stopping criteria to avoid overfitting and unphysical solutions (Diego et al., 2005; Ponente & Diego, 2011; Cha & Jee, 2022).

One of the contributing models, model 4, leverages weak gravitational lensing to construct their lens model (Cha et al., 2023), applies shape distortions of background galaxies in the galaxy cluster field as constraints on the cluster mass distribution, and extends far outside the locus of strong-lensing observables in the cluster center (e.g., Kaiser & Squires, 1993; Hoekstra et al., 2000; Jee et al., 2006; Bradač et al., 2005). By incorporating weak-lensing into strong-lensing models, it is possible to further constrain the steepness of the mass profile which directly impacts the predicted time delay (Bradač et al., 2005; Jee et al., 2006; Zitrin et al., 2015).

Lens model diversity broadly accounts for systematics regarding how the lensing mass and hence potential is parametrized, however a number of systematics are expected to affect each lens model similarly. These include the thin-sheet mass assumption (Falco et al., 1985; Birrer et al., 2016), line-of-sight structure (Dalal et al., 2005), multiplane lensing (McCully et al., 2014), millilensing from dark matter subhalos (Dai et al., 2020; Gilman et al., 2020), and the effects of large-scale structure not captured within the field of view (Wong et al., 2011).

7 Inferring $H_{0}$ Through Bayesian Analysis

We followed a Bayesian approach to constrain $H_{0}$ similar to the case of SN Refsdal in Kelly et al. (2023). For any individual lens model, $H_{0}$ can be inferred by rescaling the lens-model predicted time delays to match the photometrically and/or spectroscopically measured time delays. Their are two measured time delays, $\Delta t_{a,b}$ and $\Delta t_{c,b}$ for the image a–b and image c–b relative time delays respectively. We can infer $H_{0}$ from both time delays simultaneously, which is advantageous compared to doing so separately because the time delays may be covariant with one another.

7.1 From Time Delays Alone

We begin with a minimal set of observables made up of only the two time delays $\mathcal{O}:\{\Delta t_{a,b},\Delta t_{c,b}\}$, however this may be expanded to include other observables covariant with the time delays. For simplicity, we first consider only the photometric results. The goal is to determine the probability of $H_{0}$ given the observed light curve (LC) $P(H_{0}|LC)$, and marginalize over each lens model $M_{l}$. The light curve yields measurements of each observable in $\mathcal{O}$, while each lens model $M_{l}$ makes predictions for the same observables assuming the fiducial $H_{0}$. For each model, the fiducial $H_{0}$ is rescaled such that the model predictions best match the measurements from the light curve. Using Bayes’s theorem, we can write the probability for value of $H_{0}$ from a single lens model as:

$$p_{m}(H_{0}|{\rm LC})\propto P(H_{0})P(M_{l})\int P(\mathcal{O}|M_{l};H_{0})\\ \times P(\mathcal{O}|{\rm LC})d\mathcal{O}_{1}...d\mathcal{O}_{n}\quad,$$ (4)

where $n$ is the number of observables. Assuming the systematics of the lens models are not fully independent from one another, the overall likelihood is a weighted sum of the individual model likelihoods:

$$P(H_{0}|\rm{LC})=P(H_{0})\sum^{N}_{l}\int P(\mathcal{O}|M_{l};H_{0})\\ \times P(\mathcal{O}|{\rm LC})d\mathcal{O}_{1}...d\mathcal{O}_{n}\quad,$$ (5)

where $N$ is the number of lens models. This formalism can be understood as enforcing two priors: 1) following Eq. 4 each individual lens model will prefer values of $H_{0}$ that allows the lens model to best reproduces $\mathcal{O}$ and 2) following Eq. 5 the lens models which best reproduce the observables will have the greatest weight in the sum and hence hold the most influence over the $H_{0}$ inferred across all lens models. For a set of observables containing only the two time delays, this implies that 1) each lens model will favor values of $H_{0}$ for which the model-predicted time delays (rescaled following Eq, 3) best match the two measured time delays and 2) the lens models which best match the ratio of the two time delays are the models which effectively determine the inferred $H_{0}$.

7.2 Including Absolute Magnifications

Constraints on any observables which are correlated with the time delays effectively translate into constraints on the time delays themselves. The SN absolute magnifications are found to correlate with the time delays in the lens models, and, by including them into $\mathcal{O}$, it is possible to infer $H_{0}$ to greater accuracy and precision than time delays alone. The addition of the image a, b, and c absolute magnifications ($\mu$) increases the total number of observables to five:

$$\mathcal{O}:\{\Delta t_{a,b},\Delta t_{c,b},\mu_{a},\mu_{b},\mu_{c}\}\quad.$$ (6)

The measurements for these observables are shown in Fig. 3 and Table 1. Once again following the formalism of Eq. 4 and Eq. 5, the lens model time delay–magnification covariance implies that inducting magnification into $\mathcal{O}$ affects the $H_{0}$ inferred by each individual lens model, as well as the $H_{0}$ inferred across all lens models Because the absolute magnification was unavailable for the Type II SN Refsdal, this marks the first case where it can be meaningfully leveraged for time delay cosmography. We refer to § 10.2 for further discussion).

7.3 Including Spectroscopic Constraints

The spectroscopic data set yields measurements of all five observables independently from the photometric data set. Though the spectroscopic measurements are somewhat less precise, it is possible to constrain the observables to even greater precision through a joint-fitting approach. Assuming the spectroscopic and photometric constraints are fully independent, the likelihood of Eq. 5 can be modified to include the spectroscopic measurements, such that we determine the probability of $H_{0}$ given both the observed light curve (LC), and the observed spectroscopic evolution (SP):

$$P(H_{0}|\rm{LC,SP})=P(H_{0})\sum^{N}_{l}\int P(\mathcal{O}|M_{l};H_{0})P(\mathcal{O}|{\rm LC})\\ \times P(\mathcal{O}|{\rm Sp})d\mathcal{O}_{1}...d\mathcal{O}_{n}$$ (7)

Here $P(\mathcal{O}|{\rm LC})$ and $P(\mathcal{O}|{\rm Sp})$ refer to the posteriors for the photometric and spectroscopic constraints respectively.

Altogether, we employ three approaches to the inference of $H_{0}$ following this formalism. First, we infer $H_{0,{\rm TD-only}}$, which uses only the photometric constraints following Eq. 5, and uses only the two time delays in the set of observables, $\mathcal{O}:\{\Delta t_{a,b},\Delta t_{c,b}\}$. For the second approach, we similarly follow Eq. 5 with only photometric constraints, but include absolute magnifications, following the full set of observables given in Eq. 6 to infer $H_{0,{\rm phot-only}}$. Finally in the third approach, we follow Eq. 7 to infer $H_{0,{\rm phot+spec}}$, which includes all available constraints from photometry and spectroscopy. We note that $H_{0,{\rm TD-only}}$ inference provided identical results when spectroscopic constraints are added, and hence we do not include it as a separate approach. Across all approaches, we assign a top-hat prior for the Hubble constant: $H_{0}\in(30,130)~{}{\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$. The parametrization of $P(\mathcal{O}|M_{l};H_{0})$, $P(\mathcal{O}|{\rm LC})$, and $P(\mathcal{O}|{\rm SP})$ from the lens model, light curve fitting, and spectroscopic age-dating posteriors respectively, as well as the integration of the likelihood, are described in Appendix B.

Table 1: Magnifications and time delays From Cluster Lens Models, Photometry, and Spectroscopy

#	Reference	$\Delta t_{a,b}$	$\Delta t_{c,b}$	$|\mu_{a}|$	$|\mu_{b}|$	$|\mu_{c}|$
1	GLAFIC (Oguri, 2010, 2021)	$-105.18^{+5.16}_{-7.89}$	$-50.71^{+3.37}_{-5.09}$	$8.02^{+0.64}_{-0.57}$	$12.23^{+1.30}_{-1.51}$	$9.32^{+0.80}_{-0.87}$
2	Zitrin-analytic (Zitrin et al., 2015)	$-105.52^{+5.09}_{-6.32}$	$-41.06^{+14.27}_{-13.19}$	$11.25^{+1.05}_{-0.90}$	$16.03^{+1.77}_{-1.81}$	$14.48^{+1.91}_{-1.65}$
3	LENSTOOL (Kneib et al., 2011)	$-105.68^{+1.51}_{-1.48}$	$-52.27^{+2.06}_{-2.44}$	$6.67^{+0.11}_{-0.14}$	$9.82^{+0.23}_{-0.31}$	$8.94^{+0.21}_{-0.29}$
4	MARS (Cha & Jee, 2022)	$-136.27^{+20.23}_{-20.32}$	$-63.72^{+29.24}_{-28.46}$	$6.82^{+0.50}_{-0.44}$	$9.17^{+0.80}_{-0.72}$	$7.55^{+0.86}_{-0.70}$
5	Chen et al. (2020)	$-112.25^{+6.43}_{-6.57}$	$-53.35^{+2.71}_{-2.99}$	$6.52^{+0.24}_{-0.22}$	$10.35^{+0.46}_{-0.41}$	$6.68^{+0.24}_{-0.22}$
6	WSLAP+ (Diego et al., 2005)	$-228.33^{+19.74}_{-19.92}$	$-124.42^{+23.75}_{-23.07}$	$28.37^{+6.34}_{-7.15}$	$63.70^{+16.70}_{-17.34}$	$36.04^{+9.57}_{-9.52}$
7	Zitrin-LTM (Zitrin et al., 2009)	$-96.45^{+5.41}_{-5.85}$	$-27.64^{+6.52}_{-4.80}$	$5.66^{+0.15}_{-0.14}$	$9.77^{+0.47}_{-0.51}$	$8.80^{+0.67}_{-0.46}$
Photometry	Pierel et al. (2024)a	$-116.60^{+10.83}_{-9.29}$	$-48.30^{+3.58}_{-3.97}$	$4.43^{+1.52}_{-1.60}$	$8.00^{+3.42}_{-2.34}$	$6.43^{+1.25}_{-1.13}$
Spectroscopy	Chen et al. (2024)a	$-122^{+43.7}_{-43.8}$	$-49.3^{+12.2}_{-14.7}$	$10.93^{+8.63}_{-5.16}$	$13.22^{+7.49}_{-2.33}$	$7.14^{+1.53}_{-1.66}$

Note. — Top: The seven contributing cluster lens models, their primary references, and predicted time delays (in days) and magnifications given as the 68% confidence interval. Bottom: Measurements of the corresponding quantities from photometry and spectroscopy.

Figure 3: The absolute magnifications and relative time delays predicted by each of the seven contributing lens models (markers) and the measured quantities from photometric light curve fitting (vertical black lines and gray swaths represent the median and 68% confidence interval respectively) and spectroscopic age-dating (orange lines and swaths). The middle row is zoomed in for visual clarity. The empty purple markers of ‘Model 2*’ represent the blinded results from Model 2 which contained an error in the time delays that was amended after unblinding (see § 8). The lens model predicted time delays assume a fiducial $H_{0}=70~{}{\rm km}~{}{\rm s}^{-1}~{}{\rm Mpc}^{-1}$. Following Eq. 3, the time delays are inversely proportional to $H_{0}$, and rescaled for a range of values of $H_{0}$ to infer the total probability distribution.

8 Blinded Analysis

Blinding protocols were installed to reduce the impact of human bias. They established the information that could shared between subgroups of our team, the live unblinding process, and the rules for which changes were allowable post-unblinding. The values of the time delays and the implied magnifications from both light curve analysis and spectroscopic analysis were not revealed to one another, nor to any of the lens modeling teams. The modelers assembled regularly to agree upon the set lensing constraints, such as the catalogs of lensed image positions, spectroscopic redshifts, and photometry of the cluster galaxies as stated explicitly in Frye et al. (2024), and then connstructed the individual lens models independently. More specifically, for the blinded inference of $H_{0}$, the relative time delays and absolute magnifications from photometry and spectroscopy were blinded via a median subtraction of the probability distribution function (PDF) of each measurement. Once exception is the spectroscopic time delay and magnification measurements were unblinded to the lens model 5 predictions. This is because the spectroscopic measurements and the lens model construction were performed by the same team, however the lens model was completed and locked prior to obtaining the spectroscopic results, and hence was blinded from the spectroscopic time delays and the magnifications.

The photometric and spectroscopic analysis timelines were different. To accommodate the staggered timelines, the initial live unblinding was performed only for $H_{0,{\rm phot-only}}$, which the spectroscopic analysis remained blinded to until the $H_{0,{\rm phot+spec}}$ value was unblinded at a later date. Post-unblinding changes were only allowed to address errors. We defined an error to be a genuine mistake for which there is an objectively correct solution rather than a subjective choice. Following the $H_{0,{\rm phot-only}}$ unblinding, five errors were identified by this definition and corrected. Amending these errors ultimately resulted in a $\sim 1\sigma$ decrease in $H_{0,{\rm phot-only}}$. While the resolution of these errors impacts the inference of $H_{0}$, the methodology was unchanged, and the details of each error and its impact on both the observables and the inferred $H_{0}$ are given in Appendix C.

9 Results

Following the most straightforward approach, we inferred $H_{0}$ using only the photometric time delays, yielding $H_{0,{\rm TD-only}}=$ $71.8^{+9.8}_{-7.6}$ km s^-1 Mpc^-1. In this approach, the lens models are weighted based on their ability to reproduce the ratio between the two photometrically measured time delays. This inference is dominated by five models: model 5 (27%), model 1 (21%), model 3 (16%), model 2 (16%), and model 4 (15%). We found that including the spectroscopic time delays did not change the result due to their agreement with the photometric time delays, but $\sim 4$ times larger uncertainties.Model 6 in particular is truncated by the $130$ km s^-1 Mpc^-1 upper bound of the top-hat $H_{0}$ prior, resulting in only 3% of the weighting. Removal of this prior results in a dominant weighting of $\sim 20\%$ for $H_{0,{\rm TD-only}}$ , while the magnification-weighted approaches remain unaffected.

We next included magnification into the suite of observables, which yields $H_{0,{\rm phot-only}}=$ $74.4^{+13.0}_{-6.0}$km s^-1 Mpc^-1 from photometric constraints. Because the lens model-predicted time delays and magnifications are correlated, the magnifications influence both the $H_{0}$ inferred by an individual lens model, as well as the relative weighting each model receives in the $H_{0}$ inferred across lens models. The weighting is such that the models which best reproduce the two measured time delays and three measured absolute magnifications simultaneously are given greater contribution to the PDF of $H_{0,{\rm phot-only}}$. The result is primarily influenced by four lens models: model 5 is given 33% of the weighting, model 3 is given 28%, model 4 is given 23%, and model 1 is given 16%. No single model dominates the inference, instead the four subdominant contributions individually predict different $H_{0}$’s which ultimately broaden the resulting PDF.

Finally, we obtained the highest precision estimate by making use of all available constraints, inferring $H_{0,{\rm phot+spec}}=$ $75.4^{+8.1}_{-5.5}$ km s^-1 Mpc^-1  from joint photometric-spectroscopic constraints. This inference is dominated by model 5 ($63\%$ of the weighting), with subdominant contributions by model 1 (19%), model 4 (14%), and model 3 (4%). For each approach, the $H_{0}$ predictions of each individual model as well as their assigned weights are shown in Table 2. These results are also illustrated in Figure 4, which shows the PDF of each lens model as well as the total PDF of $H_{0}$.

We compare the inferred $H_{0}$ from SN H0pe to the results of other methods for measuring $H_{0}$ as well as SN Refsdal in Figure 5. The combined photometric-spectroscopic inference, $H_{0,{\rm phot+spec}}$, provides the strongest constraints, and is in relative agreement with most other late-time measurements, a less than $<1\sigma$ tension with measurements from SH0ES, CCHP, and TDCOSMO (Murakami et al., 2023; Freedman et al., 2019; Wong et al., 2020). By contrast, it shows marginal disagreement, a $\sim 1.4-1.5\sigma$ tension, with early-time measurements, with $P(H_{0}<67.5{\rm~{}km~{}s}^{-1}{\rm~{}Mpc}^{-1})<7\%$. Notably this is also in $\sim 1.3-1.5\sigma$ tension with the results of SN Refsdal from Kelly et al. (2023). However follow-up analysis from Liu & Oguri (2024) used the same lens modeling technique which dominated the Kelly et al. (2023) result, but varied it across range of systematics to measure $H_{0}=70.0^{+4.7}_{-4.9}$ km s^-1 Mpc^-1, in relative agreement with this work.

Table 2: $H_{0}$ Inferences

Lens Model	$H_{0,{\rm phot+spec}}$	$H_{0,{\rm phot-only}}$	$H_{0,{\rm TD-only}}$	Weightphot+spec	Weightphot-only	WeightTD-only
1	$76.05^{+6.41}_{-6.01}$	$75.55^{+6.71}_{-6.31}$	$69.74^{+7.11}_{-5.81}$	$0.19$	$0.16$	$0.21$
2	$83.55^{+6.61}_{-5.91}$	$82.05^{+7.01}_{-6.11}$	$64.83^{+7.21}_{-6.41}$	$0.00$	$0.00$	$0.16$
3	$69.94^{+4.40}_{-3.90}$	$70.54^{+4.60}_{-4.10}$	$70.94^{+5.11}_{-4.50}$	$0.04$	$0.28$	$0.16$
4	$93.96^{+15.02}_{-14.31}$	$93.66^{+14.91}_{-14.21}$	$85.96^{+16.02}_{-14.21}$	$0.14$	$0.22$	$0.15$
5	$74.54^{+5.11}_{-4.70}$	$73.84^{+5.31}_{-4.80}$	$73.84^{+5.51}_{-5.01}$	$0.63$	$0.33$	$0.27$
6	$>130$	$>130$	$>130$	$0.00$	$0.00$	$0.03$
7	$53.92^{+4.90}_{-4.40}$	$57.03^{+5.01}_{-4.50}$	$54.82^{+5.31}_{-4.80}$	$0.00$	$0.00$	$0.03$
Total	$75.4^{+8.1}_{-5.5}$	$74.4^{+13.0}_{-6.0}$	$71.8^{+9.8}_{-7.6}$	$1.00$	$1.00$	$1.00$

Note. — The seven contributing cluster lens models with their individual predictions for the three $H_{0}$ (km s^-1 Mpc^-1) inferences. Also listed are the relative weightings for each model in the likelihood across the three inferences. We also present the results for the total inference resulting from Eq. 5 and Eq. 7.

Figure 4: Inferred $H_{0}$ from comparing lens model predicted time delays and absolute magnifications to those measured from (in order of increasing constraints) only time delays (top), only photometry (middle), and combined photometry and spectroscopy (bottom). The summed PDF (dashed curves) represents the PDF of $H_{0}$ inferred by the total likelihoods of Eq. 5 and Eq. 7. Individual lens model inferences (solid curves) are shown with normalization equal to their weighting, and the median prediction of each model is plotted as markers below the curves; both the weightings and median $H_{0}$’s are given in Table 2. The median of model 6 pushes up to the $H_{0}\leq 130$ km s^-1 Mpc^-1 prior in all cases, and is instead shown with an arrow.

Figure 5: $H_{0}$ measurements from SN H0pe compared to other independent methods, broadly divided into late- and early-time measurements. SN time delay cosmography measurements include those from SN Refsdal (Kelly et al., 2023) using all contributing lens models and using a subset of two models (Grillo-g and Oguri-a*). Quasar time delay cosmography measurements include those using fully analytic mass profiles (stars+NFW; Millon et al., 2020) and free-form mass profiles (”free” mass profile; Birrer et al., 2020) from lenses within the time delay COSMOgraphy (TDCOSMO) collaboration as well as the Sloan Lens ACS (SLACS) sample (Bolton et al., 2006). Distance-ladder measurements include those from SH0ES (Murakami et al., 2023), the CCHP (Freedman et al., 2019), the EDD (Anand et al., 2022), CATs (Scolnic et al., 2023), Miras (Huang et al., 2020), SBF from SNe (Garnavich et al., 2023) and TRGBs+Cepheids (Blakeslee et al., 2021). Measurements using megamasers from the MCP (Pesce et al., 2020) and using GW170817 as a standard siren (Wang et al., 2023) are also included. Early-time measurements include those from CMB anisotropies from the Planck collaboration (Planck Collaboration et al., 2020) and from large-scale-structure clustering using SNe, baryon acoustic oscillation (BAO), and big-bang nucleosynthesis (BBN) data (Schöneberg et al., 2022). Figure is adapted from Bonvin & Millon (2020).

Figure 6: The effect of different interpolation methods for estimating the unlensed SNe Ia apparent magnitude at $z=1.783$ on $H_{0}$. We compare the inferred $H_{0}$ under three interpolation approaches: polynomial, Gaussian process, and kinematic expansion. The interpolation method affects the inferred magnification of each SN image, which in turn influences the $H_{0}$ inference. The impact is $\sim 3\%$ for both photometric-only constraints and joint photometric-spectroscopic constraints, well within the $1\sigma$ uncertainties.

Figure 7: 2D posterior distributions for time delays and absolute magnifications predicted by lens model 2 (black) and its corresponding unblinded model where only the absolute magnification measurements from light curve fitting are used as constraints in the modeling (purple). The contours represent the 68, 95, and 99% confidence intervals, and the measured values from light curve fitting are overplotted (blue lines). As seen in the black contours, the blinded lens model predicted magnifications and time delays are inversely correlated. When the magnifications are instead used as constraints in the model, the predicted time delays (purple contours) become larger, roughly following the covariance of the blinded model.

10 Discussion

10.1 Time Delay Cosmography With a SN Ia

SNe Ia are standard candles that are advantageous for time delay cosmography compared to SNe II which only provide magnification ratios. In addition, the spectroscopic evolution of SNe Ia is also standardizable unlike a SNe II. To assess the impacts of SN type, we repeat the inference of $H_{0}$, but modify the set of observables to include only constraints available to SNe II: the photometric measurements of the time delays and the magnification ratios. This exercise returns $H_{\rm 0,Type~{}II}=71.7^{+10.7}_{-6.6}$ km s^-1 Mpc^-1, providing no significant improvement over $H_{\rm 0,TD-only}$. In more detail, all lens models which satisfy the time delay ratios also satisfy the measured magnification ratios to within $1\sigma$, suggesting the magnification ratio has little leverage to break the degeneracy between the lens models and $H_{0}$. By contrast, the photometrically and spectroscopically measured absolute magnifications favor only a single lens model (Figure 4), breaking this degeneracy and improving the precision of the $H_{0}$ inference.

10.2 Lens Model Covariance

The contributing lens models exhibit varying degrees of covariance between the predicted absolute magnifications and the time delays. The likelihoods of Eq. 5 and Eq. 7 depend on both the time delays and magnifications, implying this covariance affects both the $H_{0}$ inference made by each individual lens model as well as their relative weightings. We find the lens model predicted absolute magnifications and time delays are inversely correlated, such that lower magnifications are associated with larger time delays (see 2D posteriors in Fig. 7) and hence larger values of $H_{0}$. This covariance affirms that absolute magnification is a powerful observable for inferring $H_{0}$, as is corroborated by Liu & Oguri (2024), who in the case of SN Refsdal found a similar inverse correlation.

Across the full set of models, the parametric models tended to exhibit a stronger correlation between magnification and time delays than the non-parametric models, perhaps owing to the inflexibility of fully-parametrized mass distributions enforcing a stricter parameter space. We also found that using magnification as a constraint in the lens modeling produces results consistent with the expectations of the covariance (see purple contours of Fig. 7). The specifics of these ‘unblinded’ models are detailed in Appendix D.

10.3 SN Ia Apparent Magnitude Interpolation

To ensure the magnifications for SN H0pe are measured independently of cosmology, both the spectroscopic and photometric methods compared the peak apparent magnitude of each SN image to the distribution of field SNe Ia. Due to the relatively small number of spectroscopically confirmed SNe Ia at $z>1$ (Scolnic et al., 2018), the field SN Ia peak magnitudes had to be fit with a smooth model to enable an interpolated inference at $z=1.783$. Before the un-blinding step in this analysis a choice was made to perform this fit with a Gaussian process, but after the un-blinding we explored two additional options for the fitting to understand the impact for our measurement. We fit the field SN Ia peak magnitudes with a second-order polynomial in $log(z)$-space (allowing all coefficients to vary), as well as with a second-order cosmology-independent kinematic expansion (see Riess et al., 2022, Equation 4).

We find that the range predicted by these three methods is fairly large, $\sim 0.3$mag, but that the Gaussian process result is situated directly in the center of the polynomial and kinematic expansion fits (see Figure 6 of Pierel et al., 2024). We cannot not change this method after un-blinding, as there is no evidence to suggest any of the three approaches are in error. Nevertheless we evaluate the impact of the interpolation method on the $H_{0}$ inference, as this value directly translates to the magnification measurements and therefore affects the $H_{0}$ prediction. The inferences across the three interpolation methods are shown in Figure 6; the median $H_{0}$ varies at most $\sim 3\%$ in both the photometric and the joint photometric-spectroscopic inferences. As more SNe Ia at $z>1.5$ are spectroscopically confirmed, the three methods are expected to converge while still remaining independent of cosmology.

10.4 Error Budget

The lens model predictions and the photometric-spectroscopic measurements both contribute significantly to the error budget, resulting in the $\sim 9\%$ uncertainty in the joint photometric-spectroscopic inference. This inference is dominated by lens model 5 which has both time delay uncertainties at the $\sim 5\%$ level. However, as is apparent in Table 2, model 5 alone produces $\sim 6.5\%$ uncertainty on $H_{0,{\rm phot+spec}}$, implying that the subdominant contributions of other lens models must be driving up the uncertainty to $\sim 9\%$. This could be interpreted as the systematic uncertainty between lens models, and could be reduced by the inclusion of new lensing constraints (e.g., Liu & Oguri, 2024), or by stricter magnification constraints from photometry and spectroscopy to more strongly discriminate between models.

The photometrically inferred time delays each carry $\sim 8\%$ uncertainties. Notably, the inclusion of two high precision time delays is more constraining than a single time delay of similar precision, and for the purpose of estimating the error budget can be approximated as a single effective time delay uncertainty with an up to a $\sqrt{2}$ reduction in uncertainty. Indeed if we perform the $H_{0}$ inference assuming zero uncertainty from the lens models, then the uncertainty on $H_{0}$ decreases to $\sim 5-6\%$, consistent with the single time delay with $\frac{8\%}{\sqrt{2}}\sim 5.6\%$ uncertainty (following error propagation through Eq. 3). This implies that, the total $\sim 9\%$ $H_{0}$ uncertainty is composed of a two-fifths contribution from photometry and spectroscopy, and a three-fifths contribution from lens modeling. Furthermore, the lens model uncertainty can be decomposed into one-fourth ‘measurement’ uncertainty (originating from a single model) and three-fourths ‘systematic’ uncertainty (the spread between contributing models).

Ignoring the model 5 uncertainty results in only a minor reduction in the $H_{0}$ uncertainty, and hence further reduction must come from improvements in photometry or spectroscopy. Follow-up JWST imaging after SN H0pe has faded will provide a template image for precision difference image photometry, and is expected to significantly improve magnification constraints and nearly halve photometric time delay uncertainties, resulting in as small as a $\sim 3.5\%$ effective time delay uncertainty. Given an effective lens model time delay uncertainty $\sim 7.5\%$ (which approximately accounts for the model 5 uncertainties and subdominant contributions of other lens models), this would translate to $\sim 7\%$ uncertainty on $H_{0}$. Should the improved magnification constraints prove to entirely favor a single lens model, then this could push to $\lesssim 5\%$.

11 Conclusions

SN H0pe is only the second opportunity for precision SN time delay cosmography. We measured the probability distribution of $H_{0}$ by comparing the photometrically and spectroscopically measured time delays and magnifications to those predicted by seven independently created lens models, inferring $H_{0}=$ $75.4^{+8.1}_{-5.5}$ km s^-1 Mpc^-1 . We found the lens model predicted time delays and magnifications were correlated, such that including magnifications could break degeneracies between the lens models and $H_{0}$, influencing both the $H_{0}$ prediction of each individual lens model as well as weighing its contribution to the total inference. We found this constraint was a powerful tool to select for lens model accuracy, disfavoring lens model predictions that otherwise would have been admitted had the SN been a Type II. The lens model predictions and combined photometric-spectroscopic measurements both contribute significantly to the error budget, resulting in the $\sim 9\%$ uncertainty in the joint photometric-spectroscopic inference. Follow-up epochs of JWST NIRCam imaging after SN H0pe has faded will improve the precision of the SN photometry, and hence also both the time delay and magnification estimates. This combined with the additional constraining power provided by two precision time delays may allow for constraints at the $\lesssim 5\%$ level, comparable to or exceeding that of SN Refsdal.

SN H0pe is the first cluster-lensed SN where the methodology for the $H_{0}$ inference could be developed in advance of the lens model construction and photometric and spectroscopic measurements. All lens modelers received identical lensing constraints, providing a novel opportunity to ensure consistency between lens models and test lens systematics. We found it was crucial to establish rules for blinding due to the number of independent teams contributing to the analysis, which included a formalized process to address errors. We distinguished objective errors from subjective choices, and only admitted changes to the former following unblinding. These protocols were successfully implemented for SN H0pe and provide a timely point of reference for growing the catalog of $H_{0}$ measurements from cluster-lensed SNe.

Upcoming surveys from the Vera C. Rubin and Nancy Grace Roman observatories as well as continuing JWST programs will increase the sample of cluster-lensed SNe (e.g., SN Encore), which promises to ultimately drive $H_{0}$ constraints down to percent-level precision.

12 Acknowledgements

We would like to highlight the teams who contributed the lens models for this work, which include Masamune Oguri (model 1), Ashish K. Meena, Nick Foo, and Adi Zitrin (models 2 and 7), Patrick S. Kamieneski (model 3), Sangjun Cha and M. James Jee (model 4), Wenlei Chen (model 5), and Jose M. Diego (model 6). We thank Chris McKee, Aliza Beverage, James Sullivan, and Emma Turtelboom for useful conversations. We thank the JWST Project at NASA GSFC and JWST Program at NASA HQ for their monumental effort in ensuring the success of the JWST mission. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes (MAST) at the STScI, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with JWST programs 1176 and 4446. This work is also based on observations made with the NASA/ESA Hubble Space Telescope (HST). The data were obtained from the Barbara A. Mikulski Archive for Space Telescopes (MAST) at the STScI, which is operated by the Association of Universities for Research in Astronomy (AURA) Inc., under NASA contract NAS 5-26555 for HST. This work was supported by JSPS KAKENHI Grant Numbers JP22H01260, JP20H05856, JP22K21349. M.J.J. acknowledges support for the current research from the National Research Foundation (NRF) of Korea under the programs 2022R1A2C1003130 and RS-2023-00219959. AZ acknowledges support by Grant No. 2020750 from the United States-Israel Binational Science Foundation (BSF) and Grant No. 2109066 from the United States National Science Foundation (NSF); by the Ministry of Science & Technology, Israel; and by the Israel Science Foundation Grant No. 864/23. P.L.K. acknowledges funding from NSF grants AST-1908823 and AST-2308051. AG acknowledges support from the Swedish National Space Agency, Dnr 2023-00226.

Appendix A The Seven Cluster Lens Models

All models were constructed using the same set of lensing constraints and were blind from one another as well as from the photometrically and spectroscopically measured time delays and magnifications. The exception is model 5, for which the modelers knew the spectroscopic results but did not include them in the lens-model constraints. Model 7 was constructed prior to the unblinding but was not present in the unblinding due to logistical difficulties, however the model later included into the analysis.

Model 1 implements the well-tested GLAFIC software (Oguri, 2010, 2021), which utilizes a parametric approach to solve for the analytic lens profiles. GLAFIC adopts an adaptive-mesh grid to increase computation efficiency, essentially increasing resolution in regions with larger magnification gradients (e.g., near critical curves). The model used five Navarro–Frenk–White (NFW) profiles to model the dark-matter halos of the cluster, corresponding to the two merging mass components. The cluster galaxies were modeled with pseudo-Jaffe ellipsoid profiles scaled according to each galaxy’s observed F200W flux density. To allow additional flexibility to obtain a better fit, an external shear component was included. Higher-order multipole perturbations $m=3$ and 4 (where $m=2$ represents the external shear) were also included.

Model 2 was created using a revised version of the parametric approach detailed by Zitrin et al. (2015). The new version (sometimes referred to as Zitrin-Analytic; see also Furtak et al. 2023) is no longer limited to a grid resolution. This approach operates similarly to other parametric lens-modeling techniques and has been successfully applied to JWST imaging of other clusters (e.g., Pascale et al., 2022b; Meena et al., 2023). The model employs two primary components, each described by an analytic profile. The cluster galaxies were parametrized as double pseudo-isothermal elliptical mass-density distributions (dPIE, Elíasdóttir et al., 2007) scaled by their luminosities following common relations. The cluster-scale dark matter halo was parametrized with a pseudo-isothermal elliptical mass-density distribution (PIEMD; Kassiola & Kovner, 1993). Constituting a minor component, the masses of four central galaxies were scaled independently to allow further flexibility. While the model can incorporate galaxy ellipticities, for the current model all cluster galaxy profiles were assumed to be circular. The model optimization was performed in the source plane via a Markov Chain Monte Carlo, which was sampled to derive the uncertainties. For more details on the methodology, see Furtak et al. (2023).

Model 3 was constructed using lenstool, which uses parametric mass profiles constrained simply by the locations of the members of multiply-imaged systems. In this case, the only mass profiles employed were the $z\approx 0.35$ cluster members identified by Frye et al. (2024) plus two cluster-scale halo profiles for the dominant components of the merging cluster, all having a PIEMD profile. Positions and shape parameters (ellipticity and orientation) of the cluster members were held fixed to their F200W morphologies. The mass of each cluster member was scaled in proportion to its F200W flux, following Limousin et al. (2005). A characteristic $L_{\star}$ galaxy at the cluster redshift ($m_{\rm F200W}=17.0$ AB mag) was set to a velocity dispersion $\sigma_{0}^{\star}=120~{}{\rm km~{}s}^{-1}$, a core radius $r_{\rm core}^{\star}=0.15$ kpc, and a cut radius $r_{\rm cut}^{\star}=30$ kpc. These choices are consistent with other work (e.g., Limousin et al., 2008; Johnson et al., 2014). As the impact of cluster galaxies is generally modest and poorly constrained in the model, their masses were fixed according to their scaling relative to $L_{\star}$ with the exception of the BCG, whose scaling was optimized independently. For the two cluster halos, their positions, ellipticities, orientations, velocity dispersions (mass), and core radii were left free to vary. The halos’ cut radii are poorly constrained and therefore were held fixed at 1 Mpc.

Model 4 was constructed by the MARS algorithm (Cha & Jee, 2022), which is a free-form, grid-based approach. MARS consists of two terms, one a chi-square minimization of SL multiple images in the source plane and the other a regularization term. The MARS algorithm adopts maximum cross-entropy for the regularization. This approach can suppress fluctuations and generate a quasi-unique solution for a given SL multiple images. For lens modeling, MARS can consider not only SL multiple images but also WL galaxies (Cha et al., 2023). To account for WL, MARS utilizes an additional term which aims to minimize differences between observed and predicted reduced shear. For this study, MARS provided lens models from both SL-only and SL+WL datasets to predict magnifications and time delays. The SL-only model consists of a $100\times 100$ pixel grid covering $90\arcsec\times 90\arcsec$ around the two BCGs. The SL+WL model consists of a $400\times 400$ pixel grid encompassing the entire $360\arcsec\times 360\arcsec$ FoV. One of the interesting differences between the SL-only and the SL+WL reconstructed mass distributions is that the SL-only model shows an offset of the mass peak from the northeastern core because of the deficiency of SL images, but the SL+WL mass map exhibits good agreement with the luminosity distributions of the galaxy-cluster members.

Model 5 uses the lens-modeling software originally developed for a blind prediction of the SN Refsdal reappearance and time delays (Chen et al., 2020). This approach makes minimal assumptions about the mass distribution at the galaxy-cluster scale, similar to non-parametric methods, but uses a parametric approach at the individual-galaxy scale. This semi-parametric method parameterizes the dark-matter halos of the cluster galaxies as symmetric analytic Navarro–Frenk-White profiles (Navarro et al., 1997) with the same scale radius and assumes the halo mass scales in proportion to stellar flux. For the cluster-scale dark-matter distribution, however, the model makes no assumptions about symmetry. Instead, the global lensing potential is deformed by a set of curving functions that apply numerical perturbations to the potential. These curving functions are inherently smooth on the cluster scale and are applied to minimize the separation between source-plane positions of lensed images of image families. While the only previous application of this model was to the MACS J1149.5+2223 cluster, its predictions of the reappearance position of SN Refsdal’s image SX as well as Refsdal’s four relative time delays and magnification ratios are consistent with other well-tested parametric models (Kelly et al., 2023).

Model 6 is a hybrid model, implementing the WSLAP+ (Diego et al., 2005) software, where the large-scale component of the lens model is described by a grid of Gaussians, and the small-scale component traces the light distribution of member galaxies. The model is optimized by solving a system of two linear equations per lensing constraint (one for the $x$ coordinate and one for the $y$ coordinate) of the form $\theta=\beta-\alpha(\theta,M)$ with $\theta$ the observed position of the lensed galaxy, $\beta$ the unknown position of the galaxy in the source plane, and $\alpha(\theta,M)$ the deflection angle at the position $\theta$ produced by the unknown mass distribution $M$.

For G165, the distribution of grid points for the smooth component was built recursively with a first solution obtained from a uniform grid that was then used to derive an adaptive grid that increased the spatial resolution in regions with higher mass concentrations. Using the adaptive grid, WSLAP+ derived 20 solutions varying the initial guess from the optimization and the redshifts of the systems with photometric redshifts. These models were used to derive predicted magnifications and time delays at the three observed positions of SN H0pe. The lack of lensing constraints in the east portion of the cluster and between the two groups resulted in large variations in the predicted magnification. In particular, the critical curve at $z=1.78$ exhibits a dual behaviour with some models predicting a single critical curve that connects both groups and other models predicting two disjoint critical curves, one around each group.

Model 7 uses an updated version of the Zitrin-LTM software (Zitrin et al. 2009, Zitrin et al. 2015; see also Broadhurst et al. 2005), which was modified to accommodate time delay cosmography of galaxy cluster lenses following SN Refsdal (see Zitrin, 2021). In this semi-parametric approach, it is assumed that the light traces the mass (LTM), such that observed cluster-member light is assumed to trace the underlying dark matter. The stellar mass of each cluster galaxy is represented analytically as power-law mass surface-density profile with relative weights determined by the galaxy F200W brightness. The profiles are then simultaneously smoothed with a Gaussian kernel to approximate the projected cluster-scale dark matter distribution, which receives a free weighting relative to the parametric galaxy component. An external shear component with free amplitude and position angle is also applied to broadly account for systematics like large scale structure. Finally, additional flexibility is included for individual galaxies of interest (e.g., the BCG or nearby perturbers), where the position angle, ellipticity, and mass-to-light ratio can be fit freely. As only five image systems have spectroscopic redshifts, all other systems have redshifts left free in the model, using the photometric redshift estimate as an initial guess. The minimization of the model is done in the image plane, using a MCMC.

Appendix B Probability Distribution Function Parametrization and Likelihood Integration

The probability distribution functions (PDFs) from the photometric light curve fitting, spectroscopic fitting, and the lens-model fitting are parametrized as Gaussian-mixture models (GMMs). Each fit (photometry, spectroscopy, or lens model), $F$, generated a few hundred to thousands of samples from the posterior distribution, $P(\mathcal{O}|F)$, which was then fit by a linear combination of $N$ multivariate Gaussians such that:

$$P(\mathcal{O}|F)\approx\sum^{N}_{i=1}w_{i}\mathcal{N}_{F}(\mathcal{O}|\mu_{i},\mathcal{C}_{i})\quad,$$ (B1)

where $w_{i}$ is the weight assigned to each normalized multivariate Gaussian, $\mu_{i}$ and $\mathcal{C}_{i}$ are the associated mean vectors and covariance matrices respectively, and $N$ was chosen for each model based on the Bayesian information criterion (BIC) to avoid overfitting to artifacts which may result from undersampling. Thus, the integrand in Eq. 5 becomes the product of two GMMs or equivalently the sum of products of multivariate Gaussians:

$$\int P(\mathcal{O}|M_{l};H_{0})P(\mathcal{O}|{\rm LC})d\mathcal{O}_{1}...d\mathcal{O}_{n}=\int\sum^{N}_{i=1}\sum^{N}_{j=1}w_{i}w_{j}\mathcal{N}_{M_{l}}(\mathcal{O}|\mu_{i},\mathcal{C}_{i})\mathcal{N}_{LC}(\mathcal{O}|\mu_{j},\mathcal{C}_{j})d\mathcal{O}_{1}...d\mathcal{O}_{n}\quad,$$ (B2)

where $i$ indices refer to the Gaussians which make up the GMM representing the lens model posterior, and the $j$ indices refer to the Gaussians which make up the GMM representing the light-curve fitting posterior. This formalism is advantageous in that this integral can solved analytically as the multiplication of two multivariate Gaussians is itself a multivariate Gaussian with a known renormalization.

The combined photometric–spectroscopic likelihood of Eq. 7 follows the same formalism, where the posterior from spectroscopic fitting is also parametrized as a GMM, and the integrand is the distributive product of three GMMs. Because the product of any two multivariate Gaussians is itself a Gaussian, the product of 3 multivariate Gaussians is also a Gaussian.

Appendix C Errors Corrected Post-Unblinding

Here we discuss the initial unblinded value of $H_{0}$, the five errors found subsequent to the unblinding, and how amending the errors affected the inference. These errors result from genuine mistakes which had objective fixes. The following events are stated in chronological order beginning with the unblinding event. The initial unblinding was conducted using only photometric constraints, and all errors were found while blinded to the spectroscopic results. Hence, all values given pertain to $H_{0,{\rm phot-only}}$.

The inference of $H_{0,{\rm phot-only}}$ was unblinded live by revealing the median values of the distribution of measured magnifications and time delays from light curve fitting, which yielded $H_{0,{\rm phot-only}}=81.3^{+15.9}_{-8.5}$ km s^-1 Mpc^-1.

Post-unblinding, the first error was a miscommunication between the light-curve fitting team and the team performing the $H_{0}$ inference that resulted in the exclusion of both the ${\sim}0.1$ mag intrinsic scatter of SN Ia magnitudes and the scatter due to millilensing from the inference. Applying these effects broadened the measured magnification uncertainties by ${\sim}25\%$, which resulted in an updated value of $H_{0,{\rm phot-only}}=79.6^{+15.0}_{-8.7}$ km s^-1 Mpc^-1.

Second, the light-curve fitting team found an error in the simulations which estimate the systematic uncertainty due to microlensing. The uncertainties were based on the bias and standard deviation for distributions of measured magnifications relative to the simulated magnifications. However, the measured magnifications were recorded as simply the ratio of measured amplitudes, which ignores the additional uncertainty from other light curve parameters. This was updated to match the final method used to measure the true magnifications (see Pierel et al., 2024, for details), which takes the full uncertainty into account. This further increased magnification uncertainties on images $a$ and $b$ by $\sim$10–20% while image $c$ remained mostly unaffected because of its phase at the second infrared maximum. The same simulations were found to be using a different version of the BayeSN model (Mandel et al., 2022) for simulations and fitting, which led to an artificial bias in fitted parameters. Fixing this error decreased the uncertainty of both time delays by $\sim$20%, altogether resulting in an updated value of $H_{0,{\rm phot-only}}=76.9^{+14.2}_{-8.2}$ km s^-1 Mpc^-1.

Third, the light-curve fitting team discovered a numerical error in the removal of bias induced by systematics on the magnification and time delays. The median of each retrieved magnification and time delay distribution, created from the simulations discussed above (and see Pierel et al., 2024, for more details) was removed from the measured values to account for biases identified by fitting the simulations (following Kelly et al., 2023). After unblinding, it was discovered that this bias was being removed with a sign error. This shifted the median magnifications by $\lesssim$5%, the median $\Delta t_{c,b}$ by $\sim$4%, and $\Delta t_{a,b}$ by $\lesssim$1%. These resulted in an updated value of $H_{0,{\rm phot-only}}=75.9^{+15}_{-8.4}$ km s^-1 Mpc^-1.

Fourth, the lens-modeling team for model 3 corrected an error in the sampling of the lens model results which had been used for estimation of the posterior distribution in the inference. This was caught serendipitously when the lens modeling team was requested to generate samples for other lens-model parameters to gauge the dominant source of lens-model uncertainty. The improper sampling resulted in a loss of covariance information between the magnification and time delays and significantly overestimated the time delay uncertainties by a factor of $\sim$5. Re-generating the samples correctly did not require any changes in the lens model itself and did not change the median predictions of the magnifications or time delays. The lens model was given similar weight in the $H_{0,{\rm phot-only}}$ inference both before and after making this adjustment, resulting in an updated value of $H_{0,{\rm phot-only}}=$ $74.4^{+13.0}_{-6.0}$ km s^-1 Mpc^-1.

Fifth, the lens modeling team for model 2 found an error in the calculation of the time delay from the lens model, where the source position sampled was offset by 1 pixel (0$\farcs$03) in the source plane. This did not require any modification of the underlying lens model, only a change in how time delays were measured from it. As this significantly changed the time delay distribution for this model, both the blinded values and the error-corrected values are shown Table 1 and Figure 3. Both $H_{0,{\rm phot-only}}$ and $H_{0,{\rm phot+spec}}$ were unaffected by this change as a result of negligible weighting both before and after the correction. $H_{0,{\rm TD-only}}$, however, the model weighting rose from $\sim$0 to $\sim$16%, changing the inference from $H_{0,{\rm TD-only}}=73.2^{+10.5}_{-6.4}$ km s^-1 Mpc^-1 to $H_{0,{\rm TD-only}}=$ $71.8^{+9.8}_{-7.6}$ km s^-1 Mpc^-1.

Appendix D Unblinded Lens Models

After the blinded models were done, the contributing lens-model teams also had the option to construct unblinded lens models, where absolute magnifications measured from light-curve fitting were used as constraints within the lens models. Only the teams for models 2 and 4 did so, but these two models represent both parametric and non-parametric types. Figure 7 illustrates parameter changes and their covariance for model 2. Reassuringly, the unblinded model 2 posterior follows the expectations from the correlations seen in the blinded model posterior, where the median predicted time delay increases due to the now-lower magnifications. The non-parametric model 4, however, saw very little change in its time delay predictions. The blinded version of model 4 predicts fairly accurate magnifications with the exception of $\mu_{a}$, but $\mu_{a}$ is the only magnification which lacks significant covariance between the time delays, resulting in only minor changes to the time delays in the unblinded model.

For both unblinded models, we infer only $H_{0,{\rm TD-only}}$ because magnifications became constraints in the lens models rather than predicted observables. While model 4 gives similar results in both the blinded and unblinded cases, the unblinded model 2 gives a much larger $H_{0}$ than the blinded model. Interestingly, both unblinded models predict the time delay ratio more accurately than their blinded counterparts, implying they would receive generally higher weighting in the $H_{0,{\rm TD-only}}$ inference. While it is difficult to make strong conclusions based on a sample of only two models, it is interesting to ask whether unblinding the magnifications may better address systematics between lens models, allowing for more modeling approaches to accurately reproduce the observables and hence contribute to the $H_{0}$ inference. As the sample of cluster-lensed SNe Ia grows, functionality for magnification constraints will likely become more ubiquitous in lens-modeling techniques, creating an opportunity for a more comprehensive view of lens-model systematics.

Appendix E Testing Lens Model Systematics

All seven lens models were fully blinded from the time delays, magnification ratios, and absolute magnifications inferred from the photometric light curve fitting of Pierel et al. (2024). While these measurements provide strong constraining power for weighting model contributions to $H_{0}$, they also introduced an opportunity to test lens model systematics using the absolute magnification at the galaxy cluster scale with a multiply-imaged source. Note that this is a more stringent test than lower magnification singly-imaged cases such as SN HFF14Tom in the Abell 2744 cluster field (Rodney et al., 2016). Here we broadly follow Rodney et al. (2016) to provide an in-depth discussion of the performance of the contributing models, which in turn determines the weighted contributions to the $H_{0}$ measurement.

Parametric versus Non-Parametric: While the contributing models are broadly divided into two categories, parametric and non-parametric, each model is distinct from one another. Here we compare the most free-form models, models 4 and 6, to the three classically parametric models, models 1, 2, and 3.

We find the non-parametric models have longer relative time delays (assuming the fiducial $H_{0}$ value), and hence higher $H_{0}$ inferences than the parametric models. The inferred $H_{0}$’s also have higher uncertainties. This may be owing to the greater flexibility of non-parametric models, which can explore solutions for the lens potential inaccessible to rigid fully-parametric profiles. Owing to the lack of sufficient lensing constraints about the northeastern core of the cluster, parametric methods do effectively reduce the parameter space by assuming the analytic profiles are centered at the locations of the bright cluster galaxies. As shown in Table 2, both parametric and non-parametric models receive significant weighting under all three inferences, where the weighting reflects the agreement between model predictions of observables and their corresponding measurements from photometry and spectroscopy. As the non-parametric models tended to predict longer time delays, this emphasizes the need for a wide variety of lens model approaches to fully capture the range of possible solutions.

Strong versus Strong+Weak: Model 4 makes use of both weak lensing and strong lensing constraints. While only this model was used for the $H_{0}$ inference, an additional version of the model using only strong lensing was also constructed for testing purposes only. This comparison offers a rare opportunity to evaluate the benefits weak lensing. The addition of weak lensing constraints are found to generally increase the time delay uncertainties without greatly impacting their median, while decreasing uncertainties on the magnifications and generally shifting their median lower. This significantly improves the agreement of both the lens model-predicted time delay ratios and the predicted magnifications with those measured by photometry and spectroscopy. As a result, the model receives a much greater weighting in the $H_{0,{\rm phot+spec}}$ inference, increasing the weighting to 14% for strong+weak constraints rather than $<1\%$ for strong-lensing alone. This suggests that weak lensing constraints are crucial for non-parametric modeling methods, especially when the available number of strong-lensing constraints is low.