Cosmological Results from the RAISIN Survey: Using Type Ia Supernovae in the Near Infrared as a Novel Path to Measure the Dark Energy Equation of State

The RAISIN program was carried out in cycle 20 through HST-GO 13046 (hereafter RAISIN1; PI: Kirshner) and cycle 23 through HST-GO 14216 (hereafter RAISIN2; PI: Kirshner). RAISIN1 followed 23 spectroscopically classified SNe from the Pan-STARRS MDS (Chambers et al. 2016), and RAISIN2 followed 23 spectroscopically classified SNe from DES (The Dark Energy Survey Collaboration 2005). RAISIN1 observed SNe in the redshift range of 0.22 ≤ z ≤ 0.50 and RAISIN2 observed SNe at 0.35 ≤ z ≤ 0.61 to observe at redshifts where the available HST filters overlap with the rest-frame YJH filters to minimize K-correction uncertainties. The median redshifts of PS1 and DES SNe aligned well with these targeted redshift ranges.

Due to occasional poor weather during PS1 and DES observing seasons and the need for HST template imaging after each SN had faded, RAISIN observations for each program extended over a period of 1.5−2 yr; RAISIN1 observations were taken from 2012 October 29 to 2014 June 17, and RAISIN2 observations occurred from 2015 September 28 to 2017 November 21.

To select RAISIN candidates, we identified Pan-STARRS and DES SNe that were likely to be SNe Ia discovered before maximum light. These candidates were then classified spectroscopically; classifications for MDS and DES SNe were carried out either by the RAISIN team in collaboration with the survey teams or by the MDS or DES teams themselves. For RAISIN1, we selected candidates that had been discovered by the Pan-STARRS team and were consistent with having a light-curve phase approximately 5−10 days before maximum light. For RAISIN2, we similarly selected candidates based on their apparent light-curve phase and additionally required that the photometric redshifts of their host galaxies were consistent with the target redshift range of ∼0.5 ± 0.1.

Most classifications were from Magellan (16 SNe) and Gemini South (13 SNe), with additional classifications from Gemini North (four SNe), the MMT (six SNe), the AAT (three SNe), and Keck (three SNe). Each classification was determined using the SN IDentification software (SNID; Blondin & Tonry 2007), which uses cross-correlation matching to template SNe to yield SN types, light-curve phases, and redshifts. Representative spectra for RAISIN targets are shown in Figure 1. Spectroscopically classified SNe Ia with rising light curves within the target redshift range were submitted to STScI for scheduling as nondisruptive targets of opportunity on HST. To improve the distance measurement, two additional epochs were obtained for each SN after the initial HST observation, spaced by approximately 5 rest-frame days. Template images for each SN were taken at least 6 months after the initial observations, with the exception of a single SN for which no template was needed because it was located in a region with no apparent galaxy light.

Zoom In Zoom Out Reset image size

2.2.1. The Pan-STARRS Medium Deep Survey

2.2.2. The Dark Energy Survey

2.2.3. The Carnegie Supernova Project

2.2.4. Spectroscopic Classifications and Redshifts

The HST data for RAISIN1 SNe were taken with three epochs each of F125W (approximately the J band) and F160W (approximately the H band). For RAISIN2 SNe, only the F160W band was used for most SNe; five RAISIN2 SNe included data in both filters. Figure 2 illustrates the rest-frame wavelengths probed by the HST data as a function of redshift.

Zoom In Zoom Out Reset image size

The first HST observations occurred at a median of 12.2 observer-frame days relative to B-band maximum light for RAISIN1 (+9.2 rest-frame days) and 13.9 observer-frame days (+9.5 rest-frame days) for RAISIN2 owing to the time required for SN discovery, classification, and the latency in the HST scheduling. The median S/N is 12.5 for RAISIN1 and 12.9 for RAISIN2, with template observations taken at an average of 189 days after the first SN detection. See Appendix A for a discussion of how we correct for a small amount of late-time SN flux in some template images.

Examples of RAISIN observations are shown in Figure 3. Optical photometry and spectra of these SNe are provided in the online data release accompanying this paper to allow these measurements to be of use in future NIR analyses such as SIRAH (HST-GO 15889).

Zoom In Zoom Out Reset image size

Our team measured photometry from the RAISIN images according to the following steps:

• 1.  
  For each epoch, FLT images, which are images that have been bias-subtracted, dark-subtracted, and flat-fielded, are drizzled ²⁸ together, and the final (template) epoch is subtracted.
• 2.  
  We then measured SN centroids from the difference images and performed 04 aperture photometry on the difference images with aperture corrections from HST. We verified these aperture corrections using drizzled images of the CALSPEC standard star P330E.
• 3.  
  Using artificial stars injected into the data frames, we corrected the resulting measurements and uncertainties for host galaxy noise.
• 4.  
  We calibrated the data using publicly available zero-points that were measured from observations of four white dwarf standards (GD153, G191B2B, GD71, GRW+70D5824) and the G-type star P330E. ²⁹

We describe these steps in additional detail in Appendix A. Coordinates of the RAISIN SNe and the final RAISIN photometry are given in Appendix B.

SN Ia standardization models rarely make use of the NIR owing to the small number of well-calibrated, high-cadence NIR light curves. In the past decade, only the SALT2 (Guy et al. 2007, 2010) and SiFTO (Conley et al. 2008) light-curve fitters have been used for measurements of w, with SALT2 being by far the most common model. SALT2, however, extends only to a central filter wavelength of approximately 7000 Å, with the recently developed SALT3 model (Kenworthy et al. 2021) now extending this wavelength to 8700 Å to include the rest-frame z band.

However, there are now two well-vetted light-curve models extending to the NIR: the SNooPy model (Burns et al. 2011, 2014) and the BayeSN model (Thorp et al. 2021; Mandel et al. 2022). SNooPy benefits from over a decade of use by the community and includes light-curve templates for the uBVgriYJH bands that have been trained using well-calibrated CSP light curves. BayeSN is a hierarchical Bayesian spectral energy distribution (SED) model for SNe Ia (BayeSN) that models time-dependent optical−NIR (to 1.8 μm) SN Ia SEDs as a combination of physically distinct intrinsic spectral components and host galaxy dust effects. By having an SED-based model continuous in both time and wavelength, BayeSN removes the need to compare high-redshift photometric data to an approximate rest-frame filter and effectively allows a K-correction that varies with light-curve stretch. At the same time, BayeSN has yet to be integrated in publicly available tools like SNANA (Kessler et al. 2010), so the simulation of distance-dependent biases for a cosmology analysis is not yet feasible.

Though we incorporate both SNooPy and BayeSN in our analysis, we use SNooPy as our baseline method owing to the greater ease of determining distance biases and because BayeSN was under development during the time when much of this analysis was carried out. The advantages of BayeSN’s statistical framework are somewhat mitigated owing to the limited number of NIR observational epochs in the RAISIN data. Still, we find good consistency between the distance measurements from these two models, as demonstrated in Section 5.2.1; we plan to use BayeSN in a future optical+NIR analysis of these data.

3.1.1. SNooPy Model Philosophy

The SNooPy model philosophy assumes that the observed variation in SN colors is caused by extinction in the SN host galaxy but allows the selective-to-total extinction ratio R_V to be specified by the user—the nominal method in this analysis—or fit by the data. We also allow SNooPy to fit negative extinctions for a more agnostic treatment of SN color; SNooPy assumes the luminosity versus color relationship to be linear with no minimum value (the same philosophy as the Tripp relation; Tripp 1998). This effectively means that R_V should be interpreted as a luminosity versus color trend in this work, rather than a physical dust law.

SNooPy also includes options to correct for the dependence of SN luminosity on light-curve shape. We adopted the “EBV_model2,” which uses the s_BV parameter to determine the relation between SN light-curve shape and luminosity (Burns et al. 2018). The s_BV parameter is a stretch parameter that is insensitive to extinction because it takes advantage of the fact that the delay between time of maximum light and the time of reddest color is insensitive to reddening (Lira et al. 1998). To enable easier determinations of distance biases and systematic uncertainties in later stages of the analysis, we incorporated this SNooPy model into the SNANA software (Section 3.5; Kessler et al. 2010) and used the SNooPy iYJH templates generated from low-z CSP data to fit our NIR data and measure distances. Example RAISIN light curves with SNooPy fits are shown in Figure 4.

Zoom In Zoom Out Reset image size

Because RAISIN NIR light curves include just three post-maximum light epochs per SN, we are unable to make robust shape and color corrections or determine the time of maximum light from NIR data alone; we find a Hubble residual rms of 0.24 mag when fitting A_V and 0.28 mag when fitting s_BV from NIR-only RAISIN data (in Appendix C, we will use simulations to confirm that larger rms values are expected given the small number of post-maximum light epochs per SN).

Therefore, our nominal distance measurement approach is to use the SNooPy model to fit only the NIR SN amplitude; we use optical data to determine the time of maximum light and fix the SN shape and color to a nominal value. ³⁰ This approach yields a Hubble residual rms of 0.18 mag.

We relax either the NIR-only aspect of this analysis or the nominal approach of fitting only the amplitude of the NIR data—i.e., not fitting s_BV or A_V —in several places in this analysis as described below:

• 1.  
  Selection cuts. Optical data are used to make sample selection cuts to remove SNe with high values of E(B − V), which cannot be estimated from NIR data alone. This is discussed in Section 3.4.
• 2.  
  Time of maximum light. Because we do not have NIR data near maximum light in the high-z sample, we must use optical data to constrain the phase of the SN observations.
• 3.  
  Scatter model. SNe Ia were discovered and spectroscopically classified at optical wavelengths, so we determine the consequence of these selection effects on our sample. To estimate distance-dependent selection effects for the SNooPy model, we use optical+NIR CSP data to estimate the covariance matrix of the distance residuals after estimating and correcting for s_BV and A_V . This is described below in Section 3.3.
• 4.  
  Bias corrections: To determine and correct for the redshift dependence of the intrinsic populations of s_BV , we fit s_BV using the NIR data. Though the resulting Hubble residual scatter increases, as discussed above, this is a useful tool to compare data distributions of s_BV to simulated distributions of s_BV . This is described in Appendix C.
• 5.  
  Wavelength dependence of the mass step. In Section 5.1, we compare estimates of the host galaxy mass step between optical and NIR data in a self-consistent manner that accounts for the correlation between s_BV and galaxy mass. Therefore, in this section we will correct the NIR Hubble residuals for their dependence on the optical+NIR-measured s_BV parameter. We use both the SNooPy derived s_BV −NIR luminosity correlation and the directly measured dependence of our NIR Hubble residuals on s_BV . An s_BV -corrected mass step is not used for cosmological parameter estimation but only to investigate the physics and wavelength dependence of the mass step.
• 6.  
  Comparing to optical and optical+NIR measurements of w. In Section 4.1, we compare the baseline results to results that fit the SNooPy model to optical and optical+NIR data.

Because we compare SNe Ia across a large span of redshift, we must determine distance-dependent biases, which are caused by the increased likelihood of discovering higher-luminosity SNe at larger redshifts in magnitude-limited surveys; these can be up to ∼5% in distance for the PS1 and DES samples (Scolnic et al. 2018; Brout et al. 2019b). To determine these biases, we build a dispersion covariance matrix from the optical to NIR to understand the ways in which optical SN Ia selection effects—which are caused by PS1 and DES selection criteria—affect distances measured in the NIR.

Because there has been no measurement of the covariance in magnitude between different SNooPy bands, we estimated this model from the low-z CSP data in our sample. To generate this model, we use low-z SNe at 0.015 < z < 0.1 to ensure that the results are minimally dependent on cosmological parameters while also mitigating the effect of peculiar velocity uncertainties at z < 0.015 (Peterson et al. 2021).

We then simultaneously fit the SNooPy model to all the optical and NIR bands in the low-z CSP sample to estimate the time of maximum light, s_BV , and A_V for each SN. Keeping these parameters fixed, we then fit the low-z data with the SNooPy model for each of the BVgriYJH bands individually. We use the resulting Hubble residuals in each band to generate a band-to-band covariance matrix that can be used to generate large Monte Carlo simulations and determine wavelength-dependent distance biases. We note that the additional distance scatter predicted from a peculiar velocity uncertainty of 250 km s⁻¹ (Scolnic et al. 2018) is smaller than the scatter observed in each band at the median redshift z = 0.023, and we find that the observed optical-to-NIR covariances are very similar when restricting the sample to z > 0.025 or z > 0.03.

The resulting covariance matrix is shown in Figure 5 and is now included in the public SNANA software. ³¹ We find similar dispersion measurements in most bands, with values ranging from ∼0.116 to 0.136 mag; the maximum dispersion is in the H band, and the minimum is in Y, but the band-to-band differences are not highly significant in this sample. Covariance is on the order of ∼0.015 mag² between neighboring optical bands (correlations of ≳0.95) and ∼0.010 mag² between NIR-to-NIR and optical-to-NIR bands (correlations of ∼0.78−0.86). These correlations can be understood as the correlation between Hubble residuals measured using one band at a time. Slightly higher correlations are found when we examine the Uddin et al. (2020) data, which might be due to the exclusion of fast-declining, red, and low-z SNe in the present analysis. NIR-derived distances in our sample are well correlated with optical distances—implying a significant wavelength-independent component of SN Hubble residual scatter—but are less correlated with optical bands than optical bands are with each other, implying that NIR data contain additional information useful for constraining SN Ia distances.

Zoom In Zoom Out Reset image size

Next, we apply selection cuts to our sample to ensure that each SN has well-measured photometry and distances. First, as discussed in Appendix A.1, we remove five SNe with host galaxy noise contributions that add greater than 0.1 mag to the photometric errors, based on an injected star analysis, as these are indications of subtraction residuals. Due to PSF variation caused by the breathing of HST, it is difficult to have clean subtractions on galaxies that are bright relative to the SN; this is because even few-percent variations in the PSF when multiplied by a bright galaxy core can cause large residuals relative to the brightness of the SN. This PSF variation will unavoidably bias our sample against SNe in high surface brightness environments, making it even more important to better characterize the SN–host galaxy environment correlation. We visually inspect each difference image to ensure that the remaining SNe appear to have high-quality subtractions. We also remove SN PS1-520107, which has an unclear spectroscopic classification. We note that SN 2006bt is peculiar (Foley et al. 2010), but we chose to include it for consistency across the redshift range, as similar events would be impossible to identify and exclude in our high-z sample given the noisier spectra at those redshifts. However, SN 2006bt subsequently fails Chauvenet’s criterion (see below). We further remove SNe that lack either optical or NIR data before maximum light, as this makes it impossible to accurately determine the phase of a given observation. We also ensure that pre-maximum light data exist both before and after removing >10σ SNooPy light-curve fit outliers from the data; if an SN had a misidentified time of maximum light, we found that an SN could appear to have pre-maximum light data when in reality it did not, but this issue was fixed once we removed outlying data points from the fit.

Next, we remove SNe measured to have E(B − V) >0.3 mag from optical data. Though only a single SN in our sample fails this cut, these SNe could in principle cause large outliers on our NIR-only Hubble diagram given that we perform light-curve fits with fixed A_V . High-A_V SNe could also have some sensitivity to differences in R_V that are difficult to measure owing to the extrinsic/intrinsic degeneracy in SN colors (e.g., Mandel et al. 2022). We also remove SNe with extreme values of the s_BV parameter that are not well represented in the SNooPy training sample. These include s_BV > 1.18 and s_BV < 0.75, values for which SNooPy is not well trained. Cutting at s_BV > 0.75 also excludes 91bg-like and transitional SNe Ia (see Burns et al. 2014; Gall et al. 2018), which removes SNe that would be poorly standardized by our baseline method that does not fit for s_BV . At s_BV > 1.18 in particular, only a single high-stretch SN is included in the SNooPy model training. SNe with s_BV measurement uncertainties >0.2 are also removed, as we cannot reliably estimate whether they pass cuts. Finally, two SNe—SN 2006bt and DES15C1nhv—fail Chauvenet’s criterion as 2.8σ and 3.1σ outliers on the Hubble diagram. A sample of 38 objects (the high-z sample) is expected to have just 0.07 3.1σ (or greater) outliers on average, and a sample with 44 objects (the CSP sample) is expected to have just 0.2 >2.8σ outliers; both numbers are below the Chauvenet threshold of 0.5. We treat the low- and high-z distributions differently in this calculation, as they have substantially different photometric data coverage and scatter (Section 5.2.2).

The full set of distances, light-curve parameter measurements (using optical data), and cuts on those data for both the low-z and RAISIN samples are given in Appendix E. For each subsample, selection cuts are summarized in Table 1.

Using the dispersion model from Section 3.3, we use the SNANA (Kessler et al. 2010) simulation framework to generate realizations of each subsample in this analysis (CSP, RAISIN1, and RAISIN2). SNANA is the “industry standard” tool for accurately simulating large SN Ia samples for cosmological parameter measurements, and for this reason both the MDS and DES teams have already developed sophisticated SNANA simulations of their surveys that yield accurate realizations of their survey data. SNANA includes sky noise, survey zero-points, survey cadences, host galaxy properties, realistic shape and color distributions for the SN population, and the quantitative criteria that were used to trigger follow-up observations.

These simulations are presented in Appendix C. In brief, we based our simulations on preexisting CSP (Kessler et al. 2019), PS1 (Scolnic et al. 2018), and DES (Kessler et al. 2019) simulations, but we applied them to a SNooPy model that we implemented within the SNANA framework (Figure 6). We then estimated intrinsic population distributions of the s_BV values using the aggregate NIR data, and although NIR data are only weakly dependent on A_V , we simulate a nominal exponential A_V scale of τ = 0.2 mag and vary this scale length in our systematic error budget. This exponential scale is well matched to our measurements from optical observations, from which we estimate τ ≃ 0.15−0.18 mag; preliminary optical+NIR analysis of the RAISIN sample from BayeSN also gives τ = 0.21 ± 0.04.

Zoom In Zoom Out Reset image size

The simulations then give the distance bias: the average difference between simulated and measured distance as a function of redshift, shown in Figure 7. We correct the measured distances for the average bias at each SN redshift. The predicted optical distance biases have similar sizes to bias corrections generated for the Pantheon analysis (Scolnic et al. 2018) using the Guy et al. (2010) optical scatter model and the SALT2 standardization model. The NIR bias corrections are dominated by the difference in mean stretch and A_V between the low- and high-z samples. The uncertainties on the differences in these distributions are used to determine systematic uncertainties on the bias corrections.

Zoom In Zoom Out Reset image size

The step-like dependence of SN Hubble residuals on host galaxy mass has been well established (Kelly et al. 2010; Lampeitl et al. 2010; Sullivan et al. 2010), and recently there have been reported—but somewhat conflicting—detections of the host mass step in the NIR (Uddin et al. 2020; Johansson et al. 2021; Ponder et al. 2021). Constraining the wavelength dependence of the mass step could help to distinguish among theories for the origin of the step, e.g., dust properties (Brout & Scolnic 2021) versus progenitor metallicity (Rose et al. 2021b). We attempt to measure and correct for a host galaxy mass step in the RAISIN data. We measure host galaxy masses ourselves to ensure consistency across the sample using a procedure described in Appendix D. We simultaneously estimate the mass step and the average Hubble residual in three redshift bins to ensure that our measurement is independent of cosmological parameters. The resulting masses are shown in Figure 8.

Zoom In Zoom Out Reset image size

We treat systematic uncertainties in much the same way as previous cosmological analyses (e.g., Scolnic et al. 2018; Brout et al. 2019b; Jones et al. 2019): we determine 1σ uncertainties in data and analysis parameters (e.g., peculiar velocities, Milky Way (MW) reddening, calibration, and bias corrections), apply each systematic to the data or analysis, and construct a covariance matrix of the distances resulting from each analysis choice or data modification. The systematic uncertainty covariance matrix is then added to the (diagonal-only) statistical covariance matrix. The systematic uncertainty covariance matrix is given by

for a given set of systematics S_n applied to each light curve or distance measurement. For the nth systematic, is the change in distance after applying systematic S_n to an SN at redshift z_j . The size of each systematic uncertainty is given by σ(S_n ).

Details about individual systematic uncertainties are given below. Compared to previous analyses, our choice of the SNooPy model and the limited sample size allows us to omit a few second-order systematic uncertainties, including redshift evolution in nuisance parameters and redshift evolution of the host galaxy mass step. We also do not include the choice of SN Ia dispersion model as a systematic uncertainty because we have constrained this directly from the low-z Hubble diagram itself. We are therefore not subject to the G10/C11 scatter model degeneracy that causes significant uncertainties in optical-only analyses (e.g., Brout et al. 2019b). The redshift dependence of a subset of our systematic uncertainties is shown in Figure 9.

Zoom In Zoom Out Reset image size

3.7.1. Calibration

3.7.2. Bias Corrections

3.7.3. NIR SN Standardization Model

3.7.4. Mass Step

3.7.5.  K-corrections

3.7.6. Peculiar Velocities and MW E(B−V)

We use a combination of CosmoMC (Lewis & Bridle 2002) and CosmoSIS ³⁴ (Zuntz et al. 2015) with likelihood chains from Planck Collaboration et al. (2020), including temperature, polarization, and lensing, to measure cosmological parameters from the vector of SN distances and the systematic uncertainty covariance matrix. To measure SN-only constraints, we use CosmoMC to sample the likelihood, while for combining SNe with Planck, we use the CosmoSIS importance sampler to combine SN constraints with results from the Planck chains. Luminosity distances d_L from the wCDM model, in Mpc, are given by

Here Ω_m , Ω_Λ, and Ω_k are the cosmic matter density, dark energy density, and spatial curvature, respectively. We use the following function to estimate cosmological parameters:

where C is the covariance matrix from the combined statistical and systematic uncertainties, and is the distance modulus measured from the data as described in Section 3.2 and including the bias corrections and host mass step as discussed in Sections 3.5 and 3.6. The model distance modulus . is the offset between the assumed SNooPy absolute magnitude of an SN Ia, which is also degenerate with the assumed value of H₀, and the best-fit global value.

By using CosmoMC to constrain cosmological parameters from SNe alone, Figure 11 shows the constraints on nonflat CDM relative to Pantheon and the cosmic acceleration discovery sample from Riess et al. (1998). These data thereby confirm cosmic acceleration with high significance, using a new method and wavelength range that has less sensitivity to dust extinction uncertainties. In a flat universe, from RAISIN SNe alone, we find Ω_m = 0.258 ± 0.090.

Zoom In Zoom Out Reset image size

With CMB constraints from Planck Collaboration et al. (2020) we find a dark energy equation of state of 1 + w = −0.17 ± 0.12 (Figure 12). This measurement is consistent with the previous NIR measurement from Freedman et al. (2009) of 1 + w = −0.05 ± 0.13 (stat) ± 0.09 (sys) but with smaller uncertainties.

Zoom In Zoom Out Reset image size

Statistical and systematic uncertainties in this analysis are approximately equal. The dominant systematic uncertainty in this measurement is the bias correction—predominantly due to uncertainty in the intrinsic s_BV distribution—which will improve with larger sample sizes and improved light-curve fitting methods. The second-largest systematic uncertainty is the size of the NIR mass step (Section 5.1), which will similarly improve with larger sample sizes. The full list of systematic uncertainties on the wCDM model is shown in Table 2 and visually in Figure 13. Though we have attempted to be comprehensive in our accounting of systematic uncertainties, we note that there may be substantial additional uncertainties in this novel measurement that will become more apparent only with larger sample sizes.

Zoom In Zoom Out Reset image size

In addition to measuring w, we measure H₀ with NIR SN data alone. We use Cepheid distances with CSP and RAISIN SN data to measure the Hubble constant from the distance ladder approach and use SN+CMB data to constrain the H₀ parameter with an “inverse distance ladder” approach (e.g., Cuesta et al. 2015) that combines the angular size of the CMB sound horizon with wCDM constraints from CSP, RAISIN, and Planck.

To measure H₀ from the distance ladder, we use publicly available Cepheid or tip of the red giant branch (TRGB) distances for all 10 SNe with CSP photometry and with independent TRGB or Cepheid measurements to calibrate the NIR luminosity of SNe Ia. Nine of these SNe—SN 2006D, SN 2007af, SN 2007on, SN 2007 sr, SN 2009Y, SN 2011iv, SN 2012fr, SN 2012ht, and SN 2015F ³⁵ —have redshifts of z < 0.01 and therefore were not included in our baseline sample for measuring w. One SN, SN 2007A, was already included in our sample, as it has a redshift of z = 0.017.

We use a combination of available Cepheid and TRGB distance calibrations (Freedman et al. 2019; Anand et al. 2022; Riess et al. 2021) and take the agnostic approach of giving all available distances as originally provided equal weight. We use Cepheid distances from Riess et al. (2021) for five SNe: SN 2006D, SN 2007A, SN 2009Y, SN 2012ht, and SN 2015F. We use the mean of available Cepheid and TRGB distances for three SNe: SN 2007af (Freedman et al. 2019; Riess et al. 2021), SN 2007sr (Freedman et al. 2019; Anand et al. 2022; Riess et al. 2021), and SN 2012fr (Freedman et al. 2019; Anand et al. 2022; Riess et al. 2021). Finally, we use the mean TRGB distances for two SNe from Freedman et al. (2019) and Anand et al. (2022), SN 2007on and SN 2011iv, which are in the same host galaxy of NGC 1404. These two SNe are fast decliners (Gall et al. 2018 find s_BV = 0.57 ± 0.01 and 0.64 ± 0.01 for SN 2007on and SN 2011iv, respectively), with decline rates below our nominal s_BV cuts, so we provide results with and without them.

We use these Cepheid and TRGB distances to calibrate the SN luminosity of the CSP and RAISIN samples following the method of Riess et al. (2021). If we include SN 2007on and SN 2011iv in our H₀ measurement, we measure H₀ = 77.5 ±2.1 km s⁻¹ Mpc⁻¹ from Cepheid+TRGB and H₀ = 76.6 ±2.6 km s⁻¹ Mpc⁻¹ from TRGB alone. Excluding the two fast-declining SNe in NGC 1404, however, we measure H₀ =75.9 ± 2.2 km s⁻¹ Mpc⁻¹ from Cepheid+TRGB and H₀ =72.9 ± 3.0 km s⁻¹ Mpc⁻¹ from TRGB distances alone. We note that the statistical significance of this TRGB measurement with a reduced sample is limited, as we have just three TRGB calibrator galaxies and the measurement is statistically consistent with other recent TRGB-based measurements (Freedman et al. 2019; Anand et al. 2022). The majority of the uncertainty derives from the error in the mean of the modest samples of NIR SN calibrators.

Because our distance measurement method does not correct for any dependence of SN luminosity on s_BV and our Hubble flow sample excludes such fast-declining objects (we require s_BV > 0.75), it is possible that these fast decliners, which make up a large fraction of the calibrator sample, could bias our results; therefore, our baseline result excludes the fast decliners.

From the inverse distance ladder in a wCDM model, we find a model parameter H₀ = 71.2 ± 3.8 km s⁻¹ Mpc⁻¹ (Figure 14). The value of H₀ is consistent with the conventional result of 67.4 ± 0.5 assuming ΛCDM with Planck data alone. Independent approaches to measuring SN distances at cosmological redshifts such as this one demonstrate the ways in which future measurements can help determine what role, if any, the present dark energy could play in the tension between local (e.g., Riess et al. 2016; Huang et al. 2020; Pesce et al. 2020; Wong et al. 2020; Foreman-Mackey 2016; Riess et al. 2021) and CMB H₀ measurements (e.g., Dhawan et al. 2020).

Zoom In Zoom Out Reset image size

Though the uncertainties on these measurements are large, we find that moving to NIR wavelengths does not appear to change the size of the difference between CMB and local H₀ measurements (see also Burns et al. 2018; Dhawan et al. 2018).

Finally, Table 3 compares our NIR-only measurement of w to measurements using optical data alone and optical+NIR data from SNooPy. In computing the optical and optical+NIR measurements, we apply the same base set of systematic uncertainties but add calibration uncertainties at the level of 3 mmag for each PS1 band, 5 mmag for each DES band, and 1% for each optical CSP band (Krisciunas et al. 2017; Burke et al. 2018; Scolnic et al. 2018).

Interestingly, we find substantial shifts when considering optical or optical+NIR data versus NIR data alone. In the NIR-only case versus the optical+NIR case in which NIR data primarily have the effect of constraining the SN color, w shifts by up to 0.11. Though the statistical significance of this shift is not entirely clear given the correlated data and some correlation in the systematic uncertainties, this shift could point to the need for revisions to the NIR or optical SN standardization model in future work, particularly given forthcoming surveys such as the Roman Space Telescope SN survey that will observe in the NIR (Rose et al. 2021a).

Although the sensitivity of the NIR-only measurement to bias correction uncertainties yields a result with significantly higher systematic errors in our baseline analysis, larger sample sizes would result in more precise bias estimations from a better constraint on the redshift-dependent stretch distribution of SNe. However, given that NIR-only approaches will always have difficulty constraining the potential evolution of dust extinction with redshift, we believe that the optimal approach in future work will be an optical+NIR measurement; Table 3 shows that this measurement has the lowest total systematic uncertainty budget of our three wavelength regimes. Though the overall improvement in the precision of w is small, we note that alternative distance measurement methods such as BayeSN (Mandel et al. 2022) appear to more optimally weight the NIR versus optical data and yield significantly reduced distance uncertainties when NIR data are included.

Using our nominal distance measurement method, we find a maximum likelihood size of the mass step of Δ_M =0.022 ± 0.047 mag at a step location of 10 dex and 0.026 ±0.040 mag at the Ponder et al. (2021) location of 10.44 dex (the step location is fixed during each fit). These results are used to generate the distance measurements in Section 4 above. However, the differences between the NIR and optical measurements are affected by the fact that stretch is correlated with host galaxy mass; our nominal distance method does not fit for stretch, and there is a median difference of Δs_BV =−0.099 between the high- and low-mass subsets of our sample. Correcting for this stretch differential, as done in an optical-only analysis, would have the effect of making the mass step larger.

To measure the mass step in a way that is more comparable to the way optical mass steps are measured, we test the impact of applying an s_BV correction to the NIR distances using the best-fit s_BV from optical+NIR data. We see a strong trend between Hubble residual and optical+NIR-derived s_BV as shown in Figure 15 and find that the scatter in our Hubble residuals is reduced by 10% by making this correction. We do not see a significant change in the dependence of Hubble residual on s_BV in the low-z data (slope −0.622 ± 0.189) compared to the high-z data (slope −0.702 ± 0.244).

Zoom In Zoom Out Reset image size

After s_BV fitting, we measure a larger mass step, as we would predict, of 0.072 ± 0.041. This result is shown in Figure 16. If, instead, we use the location of the mass step from Ponder et al. (2021) of 10.44 dex, we find a mass step of 0.057 ± 0.035 mag. If we adopt the approximate linear, empirical correction between s_BV and Hubble residual shown in Figure 15, instead of using the default SNooPy s_BV −luminosity relation, we find consistent steps of 0.086 ± 0.043 mag and 0.094 ± 0.037 mag for 10 and 10.44 dex step locations, respectively. We conclude that there is ∼2σ evidence for a mass step in these data even if we do not a priori assume that the SNooPy luminosity−s_BV relation is correct.

Zoom In Zoom Out Reset image size

Using the high-z data alone, we measure a mass step of 0.114 ± 0.081 mag at a step location of 10 dex and a step of 0.122 ± 0.091 mag at 10.44 dex. Unfortunately, the sample is too small to constrain the size of the step with these high-z data alone, which yield a statistically insignificant 1σ evidence for the NIR mass step.

If we use the optical data to measure the mass step with the same SNe, we find a step of 0.11 ± 0.03 mag from SNooPy and 0.08 ± 0.03 mag from SALT3 (the same step size is measured from SALT2). These are consistent with optical mass steps in the literature, which can range from 0.04 ± 0.019 mag from DES (Smith et al. 2020) to 0.119 ± 0.032 from the Nearby Supernova Factory (Aldering et al. 2002; Rigault et al. 2020); the largest data set (Pantheon; Scolnic et al. 2018) yields a step of 0.072 ± 0.01. ³⁶

These mass step results qualitatively support—albeit with limited significance—the findings of Ponder et al. (2021) and Uddin et al. (2020), who also found evidence that the NIR mass step was of order the same size as the optical mass step.

5.2.1. NIR Distances from BayeSN

Although we have not yet simulated distance biases with the BayeSN SED model, we did compare our SNooPy distances before bias correction to raw BayeSN distances that assume a fixed stretch and zero host dust extinction A_V , analogous to the approach we adopted for SNooPy. We fix the BayeSN primary intrinsic component θ = −1, approximately equivalent to our chosen SNooPy s_BV = 1. We use the same time of B-band maximum light for both SNooPy and BayeSN and fit only the NIR light-curve data to obtain photometric distance estimates.

As demonstrated in Figure 17, these distances are consistent across the redshift range of our sample, with the average difference between BayeSN and SNooPy distances changing by just 2 mmag from the low-z to high-z samples. As SNooPy has rarely been used at high redshift, this consistency between the two models gives us confidence that the raw distances that we measure with SNooPy are robust.

Zoom In Zoom Out Reset image size

5.2.2. Optical Distances from SALT3 and SNooPy

In Table 4, we explore the differences between several methods of measuring NIR, optical, and optical+NIR distances. We compare the baseline NIR-only distances to optical+NIR distances with both the baseline R_V = 1.5 and an MW-like R_V = 3.1. We also compare to optical distances from SALT3, and finally, we use an approach where optical+NIR data are first used to determine s_BV and then s_BV is fixed to that value in an NIR-only fit that measures the distance. This essentially allows the NIR distances to include the anticipated s_BV −luminosity relation.

Table 4. Optical+NIR Distance Comparisons

	Δμ (mag)
Optical+NIR (RV = 1.5)	0.011 ± 0.030
Optical+NIR (RV = 3.1)	−0.065 ± 0.032
Optical with SALT3	0.005 ± 0.038
Optical+NIR sBV , NIR dist.	−0.016 ± 0.027

Note. For optical+NIR analysis methods without bias corrections, the difference between low- and high-z distances relative to the baseline NIR-only method.

Download table as:  ASCIITypeset image

We find Hubble residual differences of ∼1%−2% between our baseline NIR method and methods using optical data and the default SNooPy value of R_V = 1.5 (Folatelli et al. 2010). With R_V = 3.1, the difference in distances is much larger, at 0.065 mag between low and high redshift, though we note that using this value of R_V increases the Hubble residual scatter by 28%. In an NIR-only analysis, however, using R_V = 3.1 changes the distances by just ∼0.02−0.04 mag owing to the lower sensitivity of the NIR to the value of R_V . Because we allow negative extinction in this analysis and because SNooPy does not model the intrinsic correlation between color and luminosity, the R_V parameter is more analogous to the nuisance parameter β in the Tripp relation than to a physical dust law.

The difference between our baseline distances and SALT3 distances is just 0.005 ± 0.038 mag. Although this difference in distance is small, the bias corrections shown in Figure 7 are ∼0.05 mag smaller in the NIR-only distances, increasing the potential discrepancy between SALT3 and the baseline NIR measurement. However, given the large uncertainties on the distance in measurement, a SALT3-based (or SALT2-based) analysis would still yield a consistent measurement of w.

Finally, the procedure of adjusting the NIR distances by a value of s_BV measured in the optical results in a statistically insignificant −0.016 ± 0.027 mag change. The lack of a small s_BV correction in our baseline analysis therefore does not appear to have a significant impact on the resulting cosmological parameter estimation. We also explore altering the semiarbitrary s_BV cuts in our baseline distance measurements, particularly given that the s_BV values tend to cluster near ∼1.18, perhaps due to a limited high-s_BV training sample in SNooPy. Fortunately, we find that our analysis is relatively insensitive to the minimum and maximum s_BV ; both a “loose” cut of 0.6 < s_BV < 1.3 and a “tight” cut of 0.8 < s_BV < 1.15 change the average Hubble residual by just 7 mmag when comparing SNe at z < 0.1 to those at z > 0.1. The mass step measurement is changed by less than 0.01 mag.

In Figure 18, we examine the dispersion of Hubble residuals for different methods that use optical data. We find that the Hubble residual dispersion is lowest when using the optical+NIR with R_V = 1.5 (0.140 mag) or the SALT3 model (0.130 mag). SALT3 appears to standardize SNe Ia more precisely than SNooPy at high redshift, while SNooPy distances are more precise than SALT3 distances for the CSP data on which SNooPy was trained. The SALT3 results are nearly the same as those from SALT2, which was used for most previous optical-only analyses.

Zoom In Zoom Out Reset image size

If an optical versus NIR difference becomes statistically significant in future work, it could point to unrecognized systematic uncertainties in analyses with optical data, NIR data, or both. In optical analyses, the G10 scatter model is a source of significant uncertainty in correcting for distance-dependent biases (Scolnic et al. 2018; Brout et al. 2019b), and proposed alternative scatter models could shift w by 4% (Brout & Scolnic 2021). In the NIR, differences in stretch−luminosity correlations between Burns et al. (2014) and Burns et al. (2018), for example, may point to systematic uncertainties in SNooPy NIR standardizations. These results show that NIR observations constitute an extremely useful consistency check on optical analyses.

Our Hubble diagram compares SNe observed near maximum light in the low-z CSP sample to SNe observed only at later phases in the higher-z RAISIN sample; the initial NIR epoch in RAISIN occurred at +9 rest-frame days on average. This could introduce systematic uncertainty if the SN standardization model is not as accurate at late times as it is near maximum light, where most previous NIR data exist. Although we account for this systematic uncertainty by using two SN Ia standardization models, both use much of the same training data.

To investigate potential biases due to the phase at which each SN was observed in the NIR, we investigated whether there is a trend in Hubble residuals as a function of the median phase at which each SN was observed. By comparing the Hubble residuals of SNe observed at a median phase of <+ 15 days to those observed at a median phase of >+ 15 days (+15 days is approximately the sample median), we detected a “step” in the NIR Hubble residuals at 3σ significance. This step appears to show that RAISIN SNe Ia observed at a phase of >+ 15 days are fainter by 0.164 ± 0.051 mag than RAISIN SNe Ia observed at a median phase of <+ 15 days. Though this trend is potentially concerning, we note that the difference in stretch between the early- and late-observed samples is 0.12, which could explain ∼0.05 mag of the ∼0.16 mag difference (Figure 15). This selection bias could be caused by the fact that low-stretch SNe are intrinsically fainter, so it may have been harder to classify them and subsequently trigger HST observations prior to maximum light.

To inspect potential biases in the SNooPy model as a function of phase, in Figure 19 we show the SNooPy Y-band template alongside the CSP and RAISIN observations after the observed data have been converted to the rest frame. Both the high- and low-z data agree relatively well with the SNooPy model, albeit with slight offsets at +10 days < phase <+20 days of ∼0.05 mag. These offsets suggest that the late-phase NIR SNooPy model should be adjusted as more data become available but would only shift the mean high-z distance measured here by approximately 0.02 mag. The variance of the NIR data also appears to increase slightly near the phase at which most of the later RAISIN observations occur, while observations near maximum light have the lowest scatter. The slightly higher scatter of the RAISIN magnitudes compared to CSP might also be an indication of K-correction uncertainties in the high-z distances. We observe similar trends when comparing the J-band SNooPy model to our data.

Zoom In Zoom Out Reset image size

The apparent dependence of Hubble residual on the phase at which an SN was observed could also be due to subtle errors in the relationship between SN luminosity and s_BV or the fact that the NIR SN color−luminosity dependence differs significantly from the relation that would be expected from dust extinction alone. Although retraining the SNooPy model is beyond the scope of this work and will likely require significant additional rest-frame NIR data to be robustly redetermined, additional low-z data and model retraining are important ways in which future NIR analyses could be improved.

The NIR dispersion model measured both in this work and in the BayeSN SED model shows that NIR distances offer significant independent leverage to constrain distance measurements compared to those obtained only with optical data. Further, NIR data, in combination with optical, provide additional leverage to improve constraints on both R_V and A_V , as demonstrated in recent analyses that use NIR data to break the degeneracy between host galaxy dust and intrinsic SN physics (Thorp et al. 2021; Mandel et al. 2022).

In this study, we see that, even without shape and color corrections, just three epochs of NIR data at ∼10−20 days after maximum light are sufficient to constrain dark energy properties. It is possible that even fewer epochs or filters would still yield competitive cosmological constraints, though additional epochs and filters give the valuable ability to avoid data outliers and to constrain dust properties. Much of the systematic uncertainty in this analysis stems from the fact that these SNe were selected for spectroscopic follow-up in the optical. In large, untargeted NIR samples like those from Roman, this source of systematic uncertainty will vanish.

SN Ia analyses and standardization models in the NIR are currently limited by the sample sizes of well-calibrated data sets, both photometric and spectroscopic, which introduce substantial calibration systematics and uncertainties in the precision of the standardization model training itself. New, well-calibrated NIR samples from CSP-II, SIRAH, VISTA, UKIRT (Konchady et al. 2022), and eventually Roman will be needed for improved model training. Simultaneous NIR spectroscopy will also be important for a better understanding of K-corrections. The SIRAH sample in particular, with extremely well-calibrated photometry and grism spectroscopy tied to ground-based NIR observations, will be a powerful link between precise NIR cosmological analyses from space and the ground.

This analysis shows that the NIR+optical measurement of w has similar precision to the optical measurement alone. Although we find that systematic uncertainties on w are lowest when we use NIR data in combination with optical, there is no improvement in the precision of cosmological parameter measurements when adding NIR data in this analysis. This is likely because our light-curve fitting approaches are not yet optimally weighting the SN data at different wavelengths, a problem that previous studies have found can be addressed with improved uncertainty treatments in SN standardization models (Mandel et al. 2011, 2022).

The underlying cause of the host galaxy mass step is a subject of debate, with numerous explorations of alternative global host galaxy dependencies (Hayden et al. 2013; Pan et al. 2014) or relationships between Hubble residuals and local host galaxy properties near the SN location (Rigault et al. 2013; Jones et al. 2015a; Rigault et al. 2015; Jones et al. 2018b; Kim et al. 2018; Rigault et al. 2020; Roman et al. 2018; Kelsey et al. 2021). The degeneracy between physical effects, potential data reduction artifacts (Solak et al. 2021), and dust has made measuring the wavelength dependence of the host galaxy mass step, as well as other host galaxy dependencies, extremely important for next-generation cosmological measurements. Measurement of whether a step-like dependence of NIR Hubble residuals on global or local galaxy color (e.g., Roman et al. 2018; Jones et al. 2018b; Kelsey et al. 2021) would be particularly interesting to understand the effects of dust or progenitor ages on SN distance measurements. Previous measurements of the mass step in the NIR have been limited by sparse data (Ponder et al. 2021) and the need to include optical data with the NIR observations to constrain the best-fit parameters such as s_BV and A_V (Uddin et al. 2020). Flexible simulation-based approaches (e.g., Pierel et al. 2021) will be necessary to fully understand the impact of different SN–host relationships on cosmological parameter measurements.

Although the galaxy-targeted nature of the CSP sample has made measuring the mass step in this work difficult (nearly all CSP galaxies are in the “high-mass” regime), we do find evidence suggesting that the mass step in the NIR is about the size of the mass step in the optical. We note that our results are not in tension with other recent NIR mass step results, including those of Johansson et al. (2021), who saw that the mass step was insignificant in the JH bands and less significant (∼2σ) in the Y band. Our distances primarily probe the rest-frame Y band, albeit with some i- and J-band overlap, where Johansson et al. (2021) find a mass step of approximately 0.07 ± 0.03 mag. However, our results do not favor the conclusion of Brout & Scolnic (2021) that the mass step is caused by variations in dust properties, as the effect of variation in R_V should decrease by a factor of ∼3 between, e.g., rest-frame V and rest-frame Y. Further study of the origin of the mass step is clearly needed, and larger infrared samples may prove useful.

Future measurements of this and other host galaxy dependencies as a function of wavelength may be necessary to separate intrinsic color, dust, other SN Ia physics, or even multiple SN Ia progenitor models. For example, when measuring the mass step, we did not fit for individual R_V values, as was done in Johansson et al. (2021), or simultaneously constrain the population distribution of R_V values, as was done in Thorp et al. (2021); these methods can help to disentangle the role of dust in the mass step’s size. Although SNe Ia are used as an empirical tool for measuring distance, understanding the objects and their connection to stellar populations is also valuable in itself. The evolution of stellar populations and of SN Ia progenitors across the long span of cosmic time that Roman will provide gives an opportunity to understand the nature of SN progenitors and their environments and to use them for more precise cosmological measurements.

SNe Ia on top of bright host galaxies can have biased measurements or underestimated errors. The additional noise introduced by a bright host galaxy tends to create significant subtraction residuals by increasing the effect of small PSF changes between the template and SN images. We occasionally see significant dipole-like subtraction artifacts.

To model this effect, we used TinyTim to create up to 25 fake stars near the same host galaxy as each SN. We require that the fake stars are at least 07 from the SN position and from each other, which is >4 times the FWHM of the F160W PSF. We cannot measure the noise at the exact location of the SN, but we can gain a statistical sense of the noise as a function of host galaxy surface brightness from these fake stars.

The results of this test are shown in Figure 20. Small error bars give the uncertainty on the average magnitude in each surface brightness bin, and large error bars give the dispersion of the fake magnitudes in each bin. To compute more realistic uncertainties, we then subtract the measured magnitude uncertainties in quadrature from the total magnitude dispersion of the fake stars. These quantities give the extra uncertainty that should be added in quadrature to each SN magnitude as a function of the brightness of the host galaxy underneath the SN. SNe with estimated host galaxy uncertainties of greater than 0.1 mag are removed from the sample; we chose this cut based on visual inspection, which shows that such SNe typically have subtraction artifacts likely due to misalignments or “breathing” of HST and are unlikely to provide reliable, unbiased photometry. Our team attempted applying convolutional subtraction algorithms to the data including HOTPANTS (Becker 2015) and ZOGY (Zackay et al. 2016) but did not achieve improved results for these bright hosts. Five SNe were removed from the final analysis based on excessive host galaxy noise.

Zoom In Zoom Out Reset image size

The RAISIN programs took HST template imaging between 129 and 285 days after the estimated time of maximum light for subtracting from the HST images with SN light. To test that our time interval between SN images and template images was sufficient, we used early- and late-time observations of SN 2012fr, which has YJ imaging at ∼150 days after maximum light from Contreras et al. (2018). With these data, we can extrapolate to later times to estimate the residual flux of RAISIN SNe at the time the template was taken.

For all RAISIN SNe, we use this estimated template flux to correct the magnitudes of RAISIN SN observations for which the closest rest-frame band is Y. Due to the faster decline of SNe in the J band, the J-band correction is expected to be negligible. To account for variability in the typical decline rate of SNe at late times in the NIR, we created a version of our analysis in which we conservatively increase the predicted late-time template flux by 0.5 mag for each RAISIN SN and correct the photometry for this revised template flux prediction. We include this analysis variant in our systematic error budget (Section 3.7). Unfortunately, SN 2012fr is the only SN for which high-quality NIR data exist at ∼6 months after maximum light. We predict that 36 of 47 RAISIN SNe would have declined by <5 mag, giving expected magnitude corrections that are more than 1% for these SNe, with a maximum predicted correction of ∼3.7%. We note that this source of uncertainty could be eliminated by revisiting the sites of RAISIN SNe to obtain templates with zero flux from the SNe.

Finally, to verify both our photometric measurement technique and our drizzling procedure, we used TinyTim to simulate fake, geometrically distorted stars in the FLT images. Then, we drizzled those images, subtracted them by a template, and measured output photometry on the sources detected in those images.

The results of these tests are shown in Figure 21. The median F125W and F160W biases are ≤1 mmag for bright sources and ∼1% for all sources. However, because these tests use only a single epoch of data, this 1% bias may be due to inaccurate fake star centroids. In the real RAISIN data this bias is likely not present because we use multiple epochs to determine the precise PSF centroid.

Zoom In Zoom Out Reset image size

Using the RAISIN simulations, we test the assumption that NIR-only distance measurements are most precise when s_BV and A_V are fixed to a constant value. From simulations, we measure a Hubble residual rms of 0.16 mag, slightly lower than the 0.19 mag measured from the real RAISIN data, when s_BV and A_V are fixed to a constant. When we attempt to use the NIR data alone to fit for A_V , the Hubble residual rms increases to 0.28 mag, and when we fit for both s_BV and A_V , we measure an rms of 0.25 mag. However, when we use the NIR data alone to fit for s_BV , we find that the rms is unchanged from the s_BV fixed case; this may be because the simulated SNooPy model does not include the increase in scatter at later phases that we observe and discuss in Section 5.3. However, despite this modest discrepancy, both the simulations and real RAISIN data appear to show no benefit from fitting to A_V and/or s_BV with NIR data alone.