Characterising the Standardisation Properties of Type Ia Supernovae in the z band with Hierarchical Bayesian Modelling

Introduction

Type Ia supernovae (SNe Ia) are important cosmological probes: as “standardizable” candles, their distances can be precisely inferred from their light curves. Because of this property, SNe Ia provided the first evidence of an accelerating universe in 1998 (Riess et al. 1998; Perlmutter et al. 1999). Today, they are being used to place some of the tightest constraints on the expansion rate of the universe, $H_0$ (Riess et al. 2022; Freedman 2021), and the properties of the elusive “dark energy” which is causing the accelerated expansion of the universe (DES Collaboration et al. 2024).

With next generation cosmological surveys, such as the Vera C. Rubin Observatory’s Legacy Survey of Space and Time (Rubin-LSST; Ivezić et al. 2019) and the Nancy Grace Roman Space Telescope (Roman; Hounsell et al. 2018), systematic errors will soon dominate over statistics in the cosmological parameter inference error budget. The majority of the systematic uncertainties in the estimates of cosmological distance moduli from SNe Ia come from the standardisation of SN Ia light curves. One approach to improving models of SNe Ia light curves for cosmology is by studying SNe Ia light curves as a function

© 2025 The Authors of wavelength. Firstly, understanding how standard SNe Ia are in different filters may reveal that SNe Ia in certain wavelength regimes are relatively better standard candles (i.e. exhibit a narrower distribution of peak absolute magnitudes). Characterising these properties can inform how to dedicate observational follow-up of discovered events. Secondly, the properties of SNe Ia across wavelength can improve our knowledge of the astrophysics underlying these systems, such as properties of the dust in the host galaxy of a SN Ia. In this work, we present the first comprehensive analysis of the standardisation properties of SNe Ia in the $z$ band.

The process of standardising of SN Ia light curves involves correcting for a number of well-documented correlations between SN Ia luminosities and light curve shapes, colours, and host galaxy properties. These correlations include the broader-brighter relation (Phillips 1993; Guy et al. 2007; Mandel et al. 2011; Burns et al. 2014), which accounts for the more luminous nature of SNe Ia which have slower declining light curves, and effective apparent colour-luminosity relation (Tripp 1998), which can be decomposed (Mandel et al. 2017) into a correction for the intrinsically redder colour of less luminous SNe Ia and interstellar reddening due to dust in the host galaxy (Burns et al. 2014, 2018; Thorp et al. 2021; Mandel et al. 2022). Many cosmological analyses also correct for the so-called “mass step” – the tendency of SNe Ia in high-mass galaxies to be systematically brighter post-standardisation than those in low-mass galaxies (Kelly et al. 2010; Sullivan et al. 2010). Although the mass step is well-documented in the literature, its physical origin remains an open question.

To date, the samples of SNe Ia used in cosmological analyses have mainly, though not exclusively, been from optical telescopes (Jones et al. 2018; Riess et al. 2022; Rubin et al. 2023; DES Collaboration et al. 2024). Because of the variation observed in SNe Ia in the optical, standardisation corrections are necessary for precise cosmological analyses. However, SNe Ia observed in the near-infrared (NIR) wavelength regime have been suggested to be relatively better standard candles compared to optical wavelengths at peak pre-standardisation, i.e. prior to corrections for known light curve decline-rate and colour correlations with luminosity (Krisciunas et al. 2004; Nobili et al. 2005; Wood-Vasey et al. 2008; Mandel et al. 2009; Folatelli et al. 2010; Kattner et al. 2012; Barone-Nugent et al. 2012; Weyant et al. 2014; Stanishev et al. 2018; Burns et al. 2018; Dhawan et al. 2018; Avelino et al. 2019; Jones et al. 2022; Peterson et al. 2024). Specifically, the peak NIR magnitudes of SNe Ia pre-standardisation tend to show less residual scatter and a weaker dependence on the decline-rate parameter than in the optical. These studies may be sensitive to the sample and model used in the analysis, though ([31]).

Light in the NIR also benefits from being 5 to 11 times less sensitive to dust extinction than the optical (e.g. Fitzpatrick 1999). Because of the difficulty of disentangling the extrinsic host galaxy dust effects from variations in the intrinsic colour populations of SNe Ia (Thorp et al. 2021; Mandel et al. 2022; Thorp & Mandel 2022), dust in the host galaxy is a significant source of uncertainty in the standardisation of SNe Ia (Freedman et al. 2009). Therefore, if indeed SNe Ia exhibit a smaller intrinsic dispersion in their NIR peak magnitudes, this property and the increased transparency of the NIR to dust would make distance moduli estimates to SNe Ia from NIR observations more robust to these systematic uncertainties.

The weaker sensitivity to dust also makes the NIR well-suited to study the mass step, which has been suggested to arise from different dust law shapes, parameterised by $R_V$, in high and low mass host galaxies (Brout & Scolnic 2021; Popovic et al. 2023). If dust is the dominant driver of the mass step, a reduced mass step is expected in the Hubble residuals from distance estimates based on NIR data because of the reduced sensitivity of this wavelength regime to dust.² Indeed, Peterson et al. (2024) recently used simulations to examine the expected size of the mass step under the assumption of the dust explanation first put forward by Brout & Scolnic (2021) and found a reduced, though still non-zero, mass step in the NIR. Therefore, there has recently been elevated interest in studying SN Ia in the NIR to probe the origin of the mass step (Uddin et al. 2020; Ponder et al. 2021; Johansson et al. 2021; Thorp et al. 2021; Jones et al. 2022; Thorp & Mandel 2022; Uddin et al. 2024; Grayling et al. 2024; Peterson et al. 2024; Thorp et al. 2024).

Several works have examined the contribution of dust to the inferred mass step by estimating the mass step in individual bands from the optical to the NIR. The results from these analyses have varied from detections of non-zero mass steps that are consistent with the optical mass step in one or more NIR bands (Uddin et al. 2020; Ponder et al. 2021; Jones et al. 2022; Uddin et al. 2024) to a mass step in the NIR that is consistent with both zero and with the optical mass step (Johansson et al. 2021; Thorp & Mandel 2022; Peterson et al. 2024). The lack of evidence for a decrease in the size of the NIR mass step relative to the optical broadly suggests that the mass step cannot be explained in full by dust (see e.g., Childress et al. (2013a), for alternative explanations of the mass step). However, the hypothesis that differing properties of dust in high and low mass host galaxies explains the mass step cannot be ruled out due to the limited statistical precision of these studies.

The $z$ band straddles the optical and NIR wavelength regimes (rest-frame $\sim8000 - 9500$Å). Light in the $z$ band is 2.6x less sensitive to dust than optical bands. Yet, the standardisation properties and mass step of SNe Ia in the $z$ band have not been explored in the literature. If SNe Ia prove to be relatively better standard candles in the $z$ band, the advantage of this wavelength regime would be compounded by the fact that $z$-band data is more easily obtainable than data in the NIR. In particular, $z$-band data can be collected alongside the optical using charge coupled devices (CCDs). The brightness of the atmosphere in the NIR makes it difficult to obtain NIR observations from the ground to a comparable depth to the optical with CCDs. Therefore, NIR observations of SNe Ia must be collected as follow-up of SNe Ia that have already been discovered, subjecting these samples of SNe Ia to complex selection effects as a result. The exact consequences of these selection effects on cosmological inference are not well understood, and Rigault et al. (2020) cautions that the unknown impact of selection effects may already be biasing our inference of cosmological parameters.

Given these strengths of the $z$ band, we set out to determine the standardisation properties of SNe Ia in the $z$ band. For this study, we use a spectroscopically confirmed sample of SN Ia observed in the griz filters from the Foundation Supernova Survey Data Release 1 (Foley et al. 2018) and the Young Supernova Experiment Data Release 1 (Aleo et al. 2023). We fit the griz light curves from this sample with BAYESN, a hierarchical Bayesian model for the spectral energy distributions (SEDs) of SNe Ia in the $gBVrizYJH$ filters (Mandel et al. 2022; Grayling et al. 2024). This work represents the first investigation of the standardisation properties of SNe Ia in an individual band using BayeSN, including an individual-band mass step estimate, and the first $z$-band Hubble diagram.

Combining optical data with data in redder bands leads to the tightest constraints on distance estimates for SNe Ia (Mandel et al. 2022; Dhawan et al. 2023). Yet, the potential of the $z$ band to increase the precision of SNe Ia standardisation has not been explored in the literature to date. Improving our physical understanding of the empirical corrections applied to SN Ia light curves is of the utmost importance, as demonstrated by the current $5\sigma$ tension between early- and late-time estimates of $H_0$ [87] (Riess et al. 2022). This tension puts pressure on $\Lambda$CDM and, if determined to be physical, may be suggestive of new physics (e.g., Mörtsell & Dhawan 2018; Di Valentino et al. 2021, for a review). Furthermore, measurements of the present-day dark energy equation of state, $w_0$, and its evolution with time, $w_a$ from Adame et al. (2025); DES Collaboration et al. (2024); DESI Collaboration et al. [25] recently suggested deviations from $\Lambda$CDM. To sufficiently rule out the possibility that these tensions arise from an empirical error (Dhawan et al. 2025), uncertainties in SN Ia estimates of cosmological parameters must be reduced from the current 1.5% level to percent level precision (Riess et al. 2022; DES Collaboration et al. 2024).

The paper is structured as follows. §2 describes the data and sample selection cuts. In §3, we describe the methodology for fitting light curves and modelling the shape-luminosity relation and mass step. We present our findings in §4 and discuss these results in §5. In §6, we conclude. Throughout this analysis, we assume a fiducial flat cosmology with $H_0 = 73.24$ and $\Omega_m = 0.28$, as was adopted for the low-$z$ training of BayeSN (Mandel et al. 2022; Grayling et al. 2024). The code used in this work can be found on GitHub at https://github.com/erinhay/snia-zband.

Data

Data for this project comes from two surveys which were both carried out on the Panoramic Survey Telescope and Rapid Response System (Pan-STARRS) telescope (Chambers et al. 2016). The Foundation Supernova Survey (Foundation), which ran from 2015–2017, produced a sample of 225 SN Ia (Foley et al. 2018). The first data release from the ongoing Young Supernova Experiment (YSE), which had its first observations in 2019, produced an additional 292 spectroscopically confirmed SN Ia (Aleo et al. 2023).

As this work is the first Hubble diagram analysis from YSE DR1, we take additional steps to clean the YSE DR1 data to prepare it for cosmology. In particular, we remove observations that have been flagged as having masked central pixels, as a mask through the central pixel severely impacts the flux measurement for that observation. We also require that the measurement uncertainty for each observation be less than 0.3 mag. These cuts remove 7.93% of the total number of observations within the range of (-10, 40) phase days that is fit by BayeSN.

In addition, we update the redshifts and redshift uncertainties from the estimates reported in YSE DR1. In YSE DR1, a conservative redshift uncertainty of 0.01 was assumed for all spectroscopically-classified objects, though the redshift measurements themselves as reported on the NASA/IPAC Extragalactic Database³ and Transient Name Server⁴ were often of higher precision. If the SN Ia is matched to a host galaxy with a spectroscopic redshift, we proceed with the redshift and redshift uncertainty from the host. For the remaining 150 SNe Ia, we use the redshift estimate and uncertainty from the transient’s spectra. For the 66 SNe Ia that do not have reported redshift measurement uncertainties, we assume a conservative uncertainty of 0.01. The final versions of the data used in this analysis can be found on the GitHub repository associated with this paper. In Appendix A, we provide an abridged table that includes the information necessary to reproduce the analysis in this paper. The full version of this table is provided online.

Host Galaxy Mass Estimates

The host galaxy masses are estimated using Blast (Jones et al. 2024). Blast is a web application that matches a specified transient to its host galaxy, retrieves all available archival data for the host, and uses that data to estimate host properties via SED fitting. The archival host photometry that is available for the transients comes from one or some of the Two Micron All-Sky Survey (2MASS; Skrutskie et al. 2006), the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), the Galaxy Evolution Explorer (GALEX; Martin et al. 2005), Pan-STARRS (Chambers et al. 2016; Flewelling et al. 2020; Magnier et al. 2020), SDSS (York et al. 2000; Blanton et al. 2017; Ahumada et al. 2020), and/or the DESI Legacy Imaging Surveys (Dey et al. 2019). The photometry is fit using prospector, a tool for estimating stellar population properties via SED fitting with a flexible suite of models (Johnson et al. 2021). Blast uses the Prospector-$\alpha$ model within prospector, which uses a non-parametric star-formation history to model the host galaxy SED (Leja et al. 2017). Priors on all parameters are as listed in Jones et al. (2024). The posterior is sampled using simulation-based inference (SBI) via the SBI++ neural posterior estimator (Wang et al. 2023), which is trained on large simulated realizations of low-$z$ galaxies ($0 < z < 0.2$) and uses normalizing flows to efficiently sample the posterior distribution. The prospector estimates of the mass of the host galaxies for the SNe Ia are then made available through Blast.⁵

Sample Selection

We define our cosmological sample using cuts on spectroscopic classification, data quality, and BayeSN parameter estimates. We first remove the peculiar SNe Ia from our sample, which includes SNe Ia-91bg-like, SNe Ia-91T-like, SNe Ia-CSM, and SNe Iax. There is ongoing research into the cosmological utility of SNe Ia-91bg-like and SNe Ia-91T-like, which both have peculiar spectra relative to “normal” SNe Ia. Though some analyses have suggested that SNe Ia-91bg-like may be standardisable with an appropriate choice of the shape parametrisation (Burns et al. 2018; Graur 2024), we follow the precedent of past cosmological analyses (e.g. Rest et al. 2014; Betoule et al. 2014; Foley et al. 2018) and remove the SNe Ia-91bg-like from our sample. While SNe Ia-91T-like have been used in cosmological analyses, recent work has suggested this sub-class of SNe Ia are over-luminous post-standardisation (Phillips et al. 2022), so we also remove these objects from the sample.

Of the remaining 446 Foundation and YSE objects, we remove 102 additional SNe Ia outside the redshift range of $0.015 < z < 0.08$ for Foundation and $0.015 < z < 0.1$ for YSE to minimise Malmquist bias. The redshift limits for Foundation are adopted from previous cosmological analyses of this sample (Foley et al. 2018; Jones et al. 2019; Thorp et al. 2021). To determine the redshift limits

for YSE, we apply a KS test to the distributions of inferred light curve parameters – $\theta$ and $A_V$ (see § 3.1 for more details on the BayeSN model parameters) – in redshift bins of $d_z = 0.01$ relative to reference distributions from the SNe Ia in the range $0.015 < z < 0.05$. The distributions of $\theta$ and $A_V$ within each bin remain consistent with being drawn from the same distribution as the reference distribution up to redshift $z = 0.1$. Above this redshift, the distributions of these inferred parameters are found to be inconsistent with the reference distribution, indicating a change in population that likely originates from Malmquist bias. We find that our results are insensitive to the choice of redshift bin size.

We next impose a requirement on the quality of the light curve data. We require at least one observation in each of the optical $gri$ bands. We also require that in at least one of the optical bands, there be an observation before the time of maximum flux in the $B$ band, $T_0$.¹⁰ Regarding the $z$-band data, we require there be 3 observations in the $z$ band within (-10, 40) phase days to ensure we can get a reliable estimate of the distance modulus from a BayeSN fit to the $z$-band only light curve. For further details on the $z$-band data quality cut, see Appendix § 9. To ensure high quality light curve fits, we include only those objects with a reduced $\chi^2_{griz} < 3$.¹¹ Finally, we require uncertainties on $T_0$ of less than a day for the BayeSN fit to the full

$griz$ light curve data, as well as the fit to the optical only $gri$ data. We also enforce that the estimates of $T_0$ agree across these two fits to within 0.5 days.

To ensure the sample lies within the bounds of the original training set of the model, we impose cuts on the BayeSN light curve parameters. In particular, we make cuts on the light curve shape parameter, $-1.5 < \theta< 3$, and the host dust extinction, $A_V < 1$. We describe the BayeSN model in more detail in § 3.1. The final sample is left at 150 objects, with 41 from YSE and 105 from Foundation. The full set of cuts is summarised in Table 1.

	YSE		Foundation
	$N_{\mathrm{cut}}$	$N_{\mathrm{rem}}$	$N_{\mathrm{cut}}$	$N_{\mathrm{rem}}$
Total Spec Sample	–	292	–	225
Normal Ia	26	266	45	180
Redshift Cut6	81	185	18	162
$gri$ Data Requirement7	66	119	47	115
$z$-band Data Requirement8	62	57	1	114
Reduced $\chi^2 < 3$	3	54	3	111
$T_0$ Cut9	8	46	6	105
$A_V < 1$	0	46	0	105
$-1.5 < \theta< 3$	1	45	0	105
Final Sample		45		105

The distribution of redshifts and of host galaxy masses in our sample is shown in Figure 1. In the final sample, 132 SNe Ia have redshifts from a host galaxy spectra and 17 have redshifts from a SN spectra. Of those 17 SNe Ia, 14 have redshift measurement uncertainties of 0.01 and 3 have redshift measurement uncertainties of 0.001. In this figure, we additionally compare the sample from this work to NIR samples used in similar studies from the Harvard-Smithsonian Center for Astrophysics (CfAIR2; Wood-Vasey et al. 2008; Friedman et al. 2015), the intermediate Palomar Transient Facility (iPTF; Johansson et al. 2021), and the Carnegie Supernova Project I (CSP-I; Contreras et al. 2010; Stritzinger et al. 2011; Krisciunas et al. 2007). In Figure 2, we show the distribution of $\theta$ and $A_V$ for the final sample of 150 objects and in Figure 3, we show the distribution of stellar mass vs. star formation rate (SFR) for the host galaxies of the final sample of SNe Ia.

Methodology

Light Curve Fitting with BayeSN

BayeSN is a hierarchical Bayesian model for SN Ia spectral energy distributions (SEDs) in the optical to infrared from 3500 Å to 18000 Å (Thorp et al. 2021; Mandel et al. 2022; Grayling et al. 2024). BayeSN provides robust and continuous fits of SN Ia light curves with physically distinct parameters for intrinsic properties of the SN Ia and host galaxy dust effects. These parameters include a light curve shape parameter, $\theta$, which encodes the broader-brighter relationship. Further intrinsic colour variations among SNe Ia are captured by a set of fitted time- and wavelength-varying components, $\epsilon(\lambda,t)$. The amount of reddening and extinction from host galaxy dust are quantified by $A_V$, the host dust extinction, and $R_V$, the dust law shape parameter, according to the Fitzpatrick (1999) dust law. Accounting for these effects, BayeSN estimates a distance modulus, $\mu$, which has been corrected for shape and colour effects, from the SN Ia photometry. The BayeSN model and priors on the model parameters are described in more detail in Mandel et al. (2022).

For this analysis, we fit with the “T21” BayeSN model, trained within the range of 3500–9500 Å, described in Thorp et al. (2021). This model was trained on Foundation data in the $griz$ filters, so it is ideal for use on this data set, which is composed predominantly of Foundation data. Unless a different host galaxy dust treatment is specified, we assume a fixed $R_V = 2.61$ for all SNe Ia – the population mean $R_V$ learned in the training of the T21 model. We fit for residual scatter within our hierarchical Bayesian models, described in § 3.2, for the population distributions of peak magnitudes, the shape-luminosity relation, and the mass step after fitting with BayeSN, so we do not include the default achromatic scatter variance term, $\sigma_0^2$, in the T21 model when fitting the light curves, instead setting it to zero.

We use the BayeSN fits to the data in two ways: 1) to estimate the peak apparent magnitudes in an individual filter¹², and 2) to estimate the “corrected” distance modulus from the SNe Ia light curves. Here, we use “corrected” to denote that both intrinsic (i.e. shape and colour relations) and extrinsic (i.e. dust) corrections, as learned in the training of the model, have been applied to give the most precise distance modulus estimate for each SN Ia.

We carry out six different BayeSN fits of each light curve:

• individual-band: the fit to an individual-band light curve with all intrinsic parameters ($\theta$, $\epsilon$) being fit for, $A_V$ fixed to zero, and $T_0$ fixed to the value from an initial fit to the $gri$ data.

• $griz$: a fit to the full $griz$ light curve, with all parameters being fit for, and with $T_0$ fixed to the value from an initial fit to the $griz$ data.

• $gri$: a fit to the optical $gri$ only data, with all parameters being fit for, and with $T_0$ fixed to the value from an initial fit to the $gri$ data.

For the $griz$, $gri$ fits, we opt to use $T_0$ as estimated from the $griz$ and $gri$ data, respectively, to maintain consistency in the data that is used in these fits. For the fit to the individual-band data, we use $T_0$ from an initial fit to the $gri$ data. This choice maintains the use of optical only data for the individual optical filter fits. For the $z$-band fit, $T_0$ is anchored with completely independent $gri$ data. An example light curve with the $griz$, $gri$, and $z$-band only fits is shown in Figure 4.

The individual-band fits as described do not include a correction for host galaxy dust extinction. In an individual band, there is no colour information to constrain the effects of dust extinction. Without multi-band information, then, the posterior on $A_V$ will revert to the prior, which may not constrain well the dust effects for an individual SN Ia. For this reason, we fix $A_V$ to zero in the individual-band BayeSN fits to obtain dust-extinguished estimates. Additionally, to compute the individual-band magnitudes corrected for host galaxy extinction, we use $A_V$ from the BayeSN fit to the $griz$ data for each SN Ia. We assume dust extinction follows a Fitzpatrick (1999) dust law evaluated at the effective wavelength of the relevant Pan-STARRS band as implemented in the extinction package (Barbary 2017). By default, we use $R_V = 2.61$, which was determined in training of the T21 model (Thorp et al. 2021), although in § 4.3, we also explore the impact of host mass-dependent $R_V$ values on our mass step inferences.

There are alternative choices that can be made when fitting individual-band data with BayeSN, such as fixing $\theta$. To ensure that the model best fits the data in an individual band, we proceed with a fit which only uses data from that band. We test the effects of alternative fitting choices on the analysis in Appendix § 10 for completeness.

Population Inference

In this work, we sample from the posteriors of interest using a Hamiltonian Monte Carlo (HMC) algorithm implemented with NumPyro (Bingham et al. 2019; Phan et al. 2019). In particular, we use the No U-turn Sampler (NUTS; Hoffman & Gelman 2014) variant of HMC from NumPyro. We use the Gelman-Rubin statistic, $\hat{R}$, to assess the convergence of the Markov Chains (Gelman & Rubin 1992). The sampling algorithms used in this work have been made available on GitHub at https://github.com/erinhay/snia-zband.

We retrieve the cosmological luminosity distance, $\mu_{\Lambda\mathrm{CDM}}(z)$, from the redshift using functionality from astropy (Astropy Collaboration et al. [4, 5, 6]). The redshift of each SN Ia is corrected for peculiar velocities using the Cosmic Flows tools from Carrick et al. (2015). Throughout this work, we assume a peculiar velocity uncertainty of $\sigma_{\mathrm{pec}} = 150\ \mathrm{km\ s}^{-1}$ (Carrick et al. 2015). We additionally run the analysis with $\sigma_{\mathrm{pec}} = 250\ \mathrm{km\ s}^{-1}$ and find that the the conclusions of the work are unchanged.

Modelling the Absolute Magnitude Distribution

To understand how standard SNe Ia are in the $z$ band relative to the $g$, $r$, and $i$ bands, we are interested in modelling the distribution of absolute magnitudes in each of these filters. Below, we describe our hierarchical Bayesian model for the population mean magnitude at peak in some generic band, as well as the population residual scatter around this mean value, for a sample of SNe Ia.

We will first be interested in modelling the population mean “dust-extinguished” or “extinguished,” peak absolute magnitude, $M^{0}_{\mathrm{ext}}$. Here, we use “extinguished” to indicate that the magnitudes have not been corrected for the effects of dust in the host galaxy. Second, we model the population mean peak absolute magnitude after correction for host galaxy dust, which we refer to as the peak “intrinsic” absolute magnitude, $M^{0}_{\mathrm{int}}$. We use the phrase “at peak” to denote “at $T_{0}$.”

In our hierarchical Bayesian model, we assume that the latent peak extinguished/intrinsic absolute magnitude of SN Ia $s$ is independently sampled from the population distribution:

where $M^{0}_{\mathrm{ext/int}}$ is the population mean peak extinguished/intrinsic absolute magnitude in some band, and $\sigma^{2}_{\mathrm{res, ext/int}}$ is the variance around this population mean peak extinguished/intrinsic absolute magnitude. We use $\mathcal{N}(\mu,\sigma^{2})$ to denote a normal distribution with mean, $\mu$, and variance, $\sigma^{2}$.

We estimate the extinguished/intrinsic magnitude at peak for SN $s$ according to:

where $\hat{m}^{s}_{\mathrm{ext/int}}$ is the extinguished/intrinsic apparent magnitude at peak inferred from BayeSN and $\mu_{\Lambda\mathrm{CDM}}(\hat{z}^{s})$ is the cosmological luminosity distance evaluated at the CMB-frame redshift of the SN that has been corrected for bulk flows, $\hat{z}^{s}$.

The inferred extinguished/intrinsic magnitude at peak, $\hat{M}^{s}_{\mathrm{ext/int}}$, is an uncertain estimate of the latent peak absolute magnitude, $M^{s}_{\mathrm{ext/int}}$, for SN $s$. Therefore, the estimator of the peak magnitude can be related to the latent parameter according to:

where the uncertainty on the inferred peak magnitude is given by:

with $\sigma_{\mathrm{fit},s}$ denoting the uncertainty from the BayeSN fit for the apparent magnitude at peak of SN Ia $s$, and $\sigma_{\mathrm{cosmo},s}$ denoting the uncertainty due to redshift-based cosmological distances.

The uncertainty due to redshift-based cosmological distances arises from uncertainty in observed redshift, $\hat{z}^{s}$, according to:

where $\sigma_{z,s}$ is the uncertainty on the measured redshift and $\sigma_{\mathrm{pec}}$ is the uncertainty due to peculiar velocities. We assume the errors in fitting peak magnitudes, measuring redshifts, and the peculiar velocity corrections are all independent.

For the observed peak extinguished/intrinsic absolute magnitude of a single SN Ia, $\hat{M}^{s}_{\mathrm{ext/int}}$, the likelihood is given by,

where $\boldsymbol{\Theta}=(M^{0}_{\mathrm{ext/int}},\sigma_{\mathrm{res, ext/int}})$. The joint likelihood is thus given by:

for a sample of $N$ SNe Ia.

The joint posterior is then given by:

where we use the following hyperpriors on these hyperparameters:

in the analysis presented in § 4.1. We use $\mathcal{U}(a,b)$ to denote a uniform distribution with lower bound $a$ and upper bound $b$.

Modelling the Shape-Luminosity Relation

Next, we explore the shape-luminosity relation – the dependence of the peak intrinsic absolute magnitudes on the shape parameter $\theta$. We model the shape-luminosity relation with a simple linear relation, including a residual scatter term around this relation. We perform this analysis for the peak intrinsic absolute magnitudes in each band, which we compare to the shape parameter inferred from the $gri$ only BayeSN fit.

For a SN $s$ with a shape parameter, $\theta^{s}$, we model the dependence of the latent peak intrinsic absolute $z$-band magnitude on the shape parameter as:

where $a$ and $b$ are the slope and intercept of the linear model and $\sigma^{2}_{\mathrm{res,S-L}}$ is the residual scatter around the linear shape-luminosity model. We assume the latent peak intrinsic absolute magnitudes for a sample of SNe Ia are conditionally independent given the shape parameters.

The observed peak magnitude is an uncertain estimate of the true peak magnitude:

where the measurement uncertainty is given by Equation (4). We neglect the measurement uncertainties on the shape parameter for each SN Ia, as justified below.

The likelihood given the data for a single SN Ia $s$ is thus:

where $\boldsymbol{\Theta}=(a,b,\sigma_{\mathrm{res,S-L}})$ here. The joint likelihood for a sample of $N$ SNe Ia is then:

again assuming conditional independence given the hyperparameters of the model.

We can now sample the joint posterior according to:

where we will use the following hyperpriors on the shape-luminosity relation hyperparameters:

when presenting the results of this analysis in § 4.1.

Note that we are are regressing against the measured shape parameters, rather than the latent shape parameters. Neglecting measurement uncertainties on shape, the independent variable, can lead to regression dilution (Kelly 2007, 2012), whereby the slope in a linear model becomes biased towards zero due to mistreatment of uncertainties on the independent variable. We use forward simulations to test the bias on the estimated slope of the shape-luminosity relation, $a$, is small in this case. For this test, we regress the peak $z$-band magnitudes against $\theta_{gri}$. In our forward simulations, we assume that the shape parameter is drawn from a population unit normal distribution. Each latent shape parameter is assumed to be measured with a Gaussian error, with measurement uncertainty based on the distribution of $\theta$ uncertainties from BayeSN for the real data.

To determine realistic values for the shape-luminosity hyperparameters, we first fit the real data with a model that neglects measurement uncertainties on the shape parameter, as shown in Figure 5. Informed by this result, we assume $a=0.15$, $b=-18.3$ mag, and $\sigma_{\mathrm{res,S-L}}=0.15$. We then simulate 100 realisations of the data and fit each realisation with our model, which neglects measurement error on $\theta$. For each realisation of the data, we retrieve the fitted shape-luminosity slope by taking the posterior mean of $a$. We then examine the distribution of fitted slopes for the 100 realisations of simulated data. The sample mean of the distribution of fitted slopes is only slightly offset by $-0.006$ relative to the true value of $a=0.15$. The conclusion of this work are insensitive to a $\sim4\%$ level of bias.

Modelling the Mass Step

We are also interested in estimating the size of the mass step in the sample. The mass step is a difference in the population mean Hubble residual of SNe which explode in high mass (HM) host galaxies and those which explode in low mass (LM) host galaxies.

When estimating the mass step, we use the fitted distance modulus from BayeSN, rather than the peak magnitudes. The reason for this choice is: 1) the distance modulus from BayeSN is corrected for the shape-luminosity relation comprehensively learned from $gri$ training data, rather than the individual-band linear model adopted in § 3.3, and 2) it allows for a like-for-like comparison between a fit to data in multiple bands (e.g., a fit to the optical $gri$ data) and a fit to data in an individual band.

Recall that the Hubble residual of SN Ia $s$ is given by $\delta^{s}=\hat{\mu}^{s}-\mu_{\mathrm{\Lambda CDM}}(\hat{z}^{s})$, where $\hat{\mu}^{s}$ is the photometric distance modulus inferred from the BayeSN fit and $\mu_{\mathrm{\Lambda CDM}}(\hat{z}^{s})$ is the cosmological distance modulus from the measured redshift of the SN, $\hat{z}^{s}$. Assuming that the mass step is a step function which transitions at $\mathcal{M}_{\mathrm{step}}$, the Hubble residual of an SN $s$ will depend on the mass of the SN’s host galaxy, $\mathcal{M}^{s}_{*}$, according to:

where $\delta_{\mathrm{LM/HM}}$ is the population mean Hubble residual in the LM/HM host galaxy bin and $\sigma^{2}_{\mathrm{res,LM/HM}}$ is the residual population scatter of SNe absolute magnitudes in the LM/HM host galaxy bins. We assume that in each mass bin, the latent Hubble residual for each SN is independently and identically drawn according to Equation (19).

The inferred Hubble residual, which is an uncertain estimate of the true Hubble residual, is given by:

with some measurement uncertainty on the Hubble residual, given by Equation (4) where $\sigma^{2}_{\mathrm{fit},s}$ now comes from the uncertainty on the inferred distance modulus, $\hat{\mu}^{s}$, from BayeSN. Again we assume independence of all measurement errors.

For each mass bin, there are two population-level parameters to infer: $\boldsymbol{\Theta}=(\delta_{\mathrm{LM/HM}},\sigma_{\mathrm{res,LM/HM}})$. The likelihood in each mass bin is therefore given by:

for some SN Ia $s$.

We define the joint likelihood for a sample of $N$ SNe Ia to be:

which we use to define the joint posterior:

where we adopt the following uninformative hyperpriors on the hyperparameters:

for the analysis in § 4.3. To compute a posterior on the mass step, $\gamma= \delta_{\mathrm{HM}} - \delta_{\mathrm{LM}}$, we take the difference of the population mean Hubble residual in each mass bin for each step along the chain.

To account for measurement uncertainties on the host galaxy mass estimates, we take 100 samples from the host mass posteriors from Blast for each object and re-run the inference for each of these realisations of the data. We take the variance of the 100 resulting mass steps that we infer as an additional uncertainty that we add in quadrature with the uncertainty on the mass step inferred when using the posterior median host mass for each object. The contribution to the uncertainty is 0.008 mag on average.

Analysis

Standardisation of SNe Ia Peak z-band Magnitudes

To assess the utility of SNe Ia before applying empirical light curve corrections, we can estimate the population mean peak extinguished absolute magnitude and residual population scatter using the apparent magnitudes as extracted from BayeSN. We will carry out this analysis for each of the individual bands according to the methodology presented in § 3.2.1.

In Table 2, we report the residual scatter around the peak extinguished absolute magnitude by band found using § 3.3. For the $g$, $r$, and $i$ bands, we find $\sigma_{\mathrm{res, ext}}(g) = 0.37 \pm0.02$ mag, $\sigma_{\mathrm{res, ext}}(r) = 0.29 \pm0.02$ mag, $\sigma_{\mathrm{res, ext}}(i) = 0.22 \pm0.01$ mag, respectively. In the $z$ band, we find $\sigma_{\mathrm{res, ext}}(z) = 0.25 \pm0.02$ mag. These results demonstrate that the redder $i$ and $z$ bands benefit from being relatively better standard candles than the optical $g$ and $r$ bands. This finding is likely contributed to by the greater transparency of the redder wavelength regime to dust, in addition to a decreased intrinsic scatter in the $i$ and $z$ bands.

Filter	$\sigma_{\mathrm{res, ext}}$ [mag]	$M_{\mathrm{int}}^{0}$ [mag]	$\sigma_{\mathrm{res, int}}$ [mag]	$a$	$\sigma_{\mathrm{res, S-L}}$
$g$	$0.37 \pm0.02$	$-18.99 \pm0.03$	$0.30 \pm0.02$	$0.174 \pm0.021$	$0.239 \pm0.015$
$r$	$0.29 \pm0.02$	$-18.92 \pm0.02$	$0.26 \pm0.02$	$0.146 \pm0.019$	$0.208 \pm0.014$
$i$	$0.22 \pm0.01$	$-18.29 \pm0.02$	$0.20 \pm0.01$	$0.076 \pm0.017$	$0.190 \pm0.013$
$z$	$0.25 \pm0.02$	$-18.30 \pm0.02$	$0.24 \pm0.02$	$0.143 \pm0.017$	$0.184 \pm0.013$

Next, we correct the peak extinguished absolute magnitudes for

dust extinction to yield the peak intrinsic absolute magnitudes. We report the population mean intrinsic absolute magnitude and residual scatter around the mean for each band in Table 2. In the $g$-band, we find $\sigma_{\mathrm{res, int}}(g) = 0.30 \pm0.02$ mag, in the $r$-band, we find $\sigma_{\mathrm{res, int}}(r) = 0.26 \pm0.02$ mag, and in the $i$-band, we find $\sigma_{\mathrm{res, int}}(i) = 0.20 \pm0.01$ mag. We find the residual scatter in the $z$ band to be $\sigma_{\mathrm{res, int}}(z) = 0.24 \pm0.02$ mag. The residual scatter around the peak intrinsic absolute magnitude remains higher for the $g$ and $r$ bands than for the $i$ and $z$ bands. Therefore, SNe Ia appear to be relatively better standard candles in the $i$ and $z$ band than in the $g$- and $r$-band. Because of the $z$ band's greater transparency to dust, it shows only a small reduction in the residual scatter from peak magnitudes pre-dust correction to post-dust correction, as anticipated.

We next investigate the dependence of the $z$-band peak intrinsic magnitudes to the inferred shape parameter, $\theta$ from the $gri$ fit to the data, according to the methodology in § 3.3. We find the slope of the shape-luminosity relation to be $a = 0.143 \pm0.017$ with a residual scatter around the fitted linear relation of $\sigma_{\mathrm{res, S-L}} = 0.184 \pm0.013$ mag. Figure 5 shows $M_{z,\mathrm{int}}$ versus $\theta_{gri}$ with the fitted linear shape-luminosity relation. There is a statistically significant correlation between the peak intrinsic $z$-band magnitude and the shape of the light curve, as the estimated shape-luminosity slope is inconsistent with zero at the 8.4-$\sigma$ level.

The Hubble Diagram

We infer the scatter in the Hubble residuals from the corrected distance moduli to the SNe Ia. The Hubble diagrams and residuals from the individual-band BayeSN fits are shown in Figure 6. In Figure

7, we additionally show the $griz$ and $gri$ Hubble diagrams. In all Hubble diagrams, the mean Hubble residuals are consistent with zero, indicating good calibration of the model. We use the root mean square (RMS) of the Hubble residuals to estimate the total scatter in the Hubble residuals, as this metric is insensitive to modelling assumptions (e.g., assumed peculiar velocity uncertainty). To give a measure of the residual scatter that is unexplained by measurement uncertainties and peculiar velocities, we estimate the residual scatter using the methodology of § 3.2.1 applied to Hubble residuals.

The individual-band Hubble diagrams exhibit a total scatter of RMS = 0.200 mag in the $g$ band, RMS = 0.181 mag in the $r$ band, RMS = 0.175 mag in the $i$ band, and RMS = 0.195 mag in the $z$ band. Interestingly, the RMS decreases with increasing wavelength from the $g$ to $i$ bands, and then increases in the $z$ band. The increase in scatter in the $z$-band Hubble diagram appears to be driven by the YSE subset of the sample. Therefore, we attribute some of this increase to the lower median number of $z$-band observations (4 observations per object) compared to in the other filters (between 6–7 observations per object) in the sample. The Foundation subset of the sample, which has a more balanced number of observations in each bands, shows trends similar to those seen in the peak magnitudes in § 4.1 where the $z$-band peak magnitudes show a level of scatter that falls between the $r$- and $i$-band peak magnitudes. Based on the Hubble diagrams in this work, we can conclude that SNe Ia in the $z$ band are relatively worse standard candles than those in the $i$ band, and their relative standardisation compared to the $g$ and $r$ bands cannot be definitely determined from this sample of SNe Ia.

Analysing the residual scatter, $\sigma_{\mathrm{res}}$, of the Hubble residuals from the combined sample and individual surveys leads to similar conclusions. The residual scatter estimates for the combined sample are $\sigma_{\mathrm{res}}(g) = 0.145$ mag, $\sigma_{\mathrm{res}}(r) = 0.152$ mag, $\sigma_{\mathrm{res}}(i) = 0.146$ mag, and $\sigma_{\mathrm{res}}(z) = 0.157$ mag. When considering the Foundation sample on its own, the $z$ band again exhibits a level of scatter that is slightly reduced compared to the $r$ band, but slightly increased relative to the $i$ band. The YSE sample does not show as clear a trend, with the $g$-band Hubble diagram having the lowest residual scatter and the $z$-band Hubble diagram having the highest residual scatter.

We also compute an inverse variance weighted RMS (wRMS), where the variance includes uncertainties from cosmological based distances, redshift measurement, and the BayeSN fit, as well as residual scatter. The residual scatter for each Hubble diagram is taken from the $\sigma_{\mathrm{res}}$ values reported above. We find wRMS = 0.186 mag in the $g$ band, wRMS = 0.175 mag in the $r$ band, wRMS = 0.165 mag in the $i$ band, and wRMS = 0.183 mag in the $z$ band. These results lead to consistent conclusions with those drawn from the unweighted RMS.

When comparing the $gri z$ and $gri$ Hubble diagrams, we find that the RMS decreases from RMS = 0.173 mag in the $gri$ Hubble diagram to RMS = 0.168 mag in the $gri z$ Hubble diagram. The scatter in the optical $gri$ Hubble diagram is lower than the $z$-band Hubble diagram, though this result is not unexpected given differences in the amount of data going into each fit. Importantly, we find that combining the $z$-band data with the optical does improve distances estimates.

Again, the conclusions from the analysis of the residual scatter are not as clear cut, with both the $gri z$ and $gri$ Hubble diagrams having $\sigma_{\mathrm{res}} = 0.142$ mag. This finding is also driven by the YSE subset of the sample, as Foundation does show a decrease in residual scatter when including the $z$-band data; we comment on the differences between YSE and Foundation which may lead to this result in the following paragraph. Finally, the weighted RMS of the $gri z$ Hubble diagram is wRMS = 0.159 mag compared to the relatively higher weighted RMS from the $gri$ Hubble diagram of wRMS = 0.164 mag, which is a consistent result as from the unweighted RMS.

Further to the differences between the two surveys, the total and residual scatter of the YSE subset is consistently higher than that of the Foundation subset. We attribute this discrepancy to the small number of YSE objects which pass the requisite data cuts, which makes the total scatter more sensitive to outliers, the larger number of objects with redshifts derived from the SN Ia spectra, for which the assumed measurement uncertainty of $\sigma_z = 0.01$ dominates the error budget, and the lower cadence of the survey relative to Foundation. Indeed, the median number of $z$-band observations is five for objects from Foundation and two for objects from YSE, before data cuts. The reason for this difference can largely be attributed to YSE being an untargeted survey, unlike Foundation, so the survey scans the sky according to a predetermined schedule. The selection effects are better defined in untargeted surveys, which is beneficial in cosmological analyses, but there are trade-offs in the cadence of the light curves from YSE DR1. While we mitigate some of these differences with the requirement that there be four $z$-band observations, there are still differences in the resulting $z$-band light curve data of the two surveys.

The Mass Step

Finally, we estimate the mass step according to the methodology presented in § 3.3.1 using the Hubble residuals from each of the $gri z$, $gri$, and individual-band fits.

Notably, we include separate residual scatter terms for the Foundation and YSE subsets of the sample. When fitting the mass step separately for each survey, we find that the SNe Ia from the two surveys exhibit significantly different residual scatter properties. This finding is consistent with the trends in Hubble residual scatter as presented in § 4.2. We modify Equation (19) so that $\sigma_{\mathrm{res,LM/HM}}^{2}$ is dependent on whether SN $s$ is from Foundation or YSE. We still fit for one mass step for both samples, as we find that the mass step remains consistent when fit separately for each survey.

For our analysis, we compare two choices of the mass step location. Firstly, we use $\mathcal{M}_{\mathrm{step}} = 10^{10}\mathcal{M}_{\odot}$, as this value has been used in numerous previous analyses of the mass step (Sullivan et al. 2010; Jones et al. 2019; Johansson et al. 2021; Thorp et al. 2021; Thorp & Mandel 2022; Brout et al. 2022; Grayling et al. 2024). This choice places 49 objects in the low mass bin and 101 objects in the high mass bin. Secondly, we use $\mathcal{M}_{\mathrm{step}} = 10^{10.329}\mathcal{M}_{\odot}$, the median host galaxy mass in the sample. With this choice of mass step location, we get a more balanced split of the sample with 75 objects in each of the low and high mass bins.

The results of the analysis on the combined YSE and Foundation sample, and for each survey individually, are presented for both mass step locations under the fiducial assumption of a single $R_V$ for all SNe Ia in Table 3. With $\mathcal{M}_{\mathrm{step}} = 10^{10}\mathcal{M}_{\odot}$, we find $\gamma_{griz} = -0.126 \pm0.027$ mag for the Hubble residuals from the distance moduli from the $gri z$ data. The reported uncertainty includes a sub-dominant contribution arising from uncertainty in the host galaxy mass estimates. The residual scatter in the LM bin, $\sigma_{\mathrm{res,LM}} = 0.086 \pm0.048$ mag for YSE and $\sigma_{\mathrm{res,LM}} = 0.138 \pm0.022$ mag for Foundation. In the HM bin, we find $\sigma_{\mathrm{res,HM}} = 0.217 \pm0.036$ mag for YSE and $\sigma_{\mathrm{res,HM}} = 0.103 \pm0.013$ mag for Foundation. As in § 4.2, the residual scatter in the YSE subset of the sample is consistently larger than the residual scatter in the Foundation subset of the sample.

Fit	Sample	$\sigma_{\mathrm{res,LM}}$ [mag]	$\sigma_{\mathrm{res,HM}}$ [mag]	$\gamma$ [mag]	Significance
		Foundation	YSE	Foundation	YSE
$M_{\mathrm{step}}=10^{10}\,M_\odot$
$\mu_{griz}$	Foundation+YSE	$0.138 \pm0.022$	$0.086 \pm0.048$	$0.103 \pm0.013$	$0.217 \pm0.036$	$-0.126 \pm0.027$; 4.58
	Foundation	$0.138 \pm0.021$	–	$0.103 \pm0.013$	–	$-0.131 \pm0.031$; 4.16
	YSE	–	$0.094 \pm0.052$	–	$0.215 \pm0.035$	$-0.079 \pm0.073$; 1.08
$\mu_{gri}$	Foundation+YSE	$0.138 \pm0.023$	$0.080 \pm0.047$	$0.107 \pm0.014$	$0.223 \pm0.038$	$-0.117 \pm0.029$; 4.06
	Foundation	$0.137 \pm0.023$	–	$0.107 \pm0.014$	–	$-0.131 \pm0.032$; 4.04
	YSE	–	$0.090 \pm0.051$	–	$0.220 \pm0.038$	$-0.046 \pm0.074$; 0.63
$\mu_g$	Foundation+YSE	$0.159 \pm0.029$	$0.141 \pm0.072$	$0.117 \pm0.019$	$0.208 \pm0.046$	$-0.099 \pm0.037$; 2.72
	Foundation	$0.159 \pm0.028$	–	$0.116 \pm0.019$	–	$-0.124 \pm0.038$; 3.21
	YSE	–	$0.145 \pm0.079$	–	$0.192 \pm0.046$	$0.034 \pm0.095$; 0.36
$\mu_r$	Foundation+YSE	$0.124 \pm0.022$	$0.138 \pm0.065$	$0.131 \pm0.016$	$0.233 \pm0.041$	$-0.102 \pm0.031$; 3.33
	Foundation	$0.124 \pm0.022$	–	$0.131 \pm0.016$	–	$-0.120 \pm0.033$; 3.63
	YSE	–	$0.152 \pm0.074$	–	$0.217 \pm0.039$	$0.009 \pm0.085$; 0.11
$\mu_i$	Foundation+YSE	$0.150 \pm0.023$	$0.073 \pm0.045$	$0.102 \pm0.013$	$0.233 \pm0.039$	$-0.113 \pm0.028$; 4.00
	Foundation	$0.150 \pm0.022$	–	$0.101 \pm0.013$	–	$-0.125 \pm0.032$; 3.85
	YSE	–	$0.079 \pm0.050$	–	$0.222 \pm0.038$	$-0.036 \pm0.071$; 0.50
$\mu_z$	Foundation+YSE	$0.147 \pm0.024$	$0.113 \pm0.061$	$0.115 \pm0.015$	$0.275 \pm0.047$	$-0.105 \pm0.031$; 3.37
	Foundation	$0.147 \pm0.025$	–	$0.115 \pm0.015$	–	$-0.118 \pm0.034$; 3.45
	YSE	–	$0.129 \pm0.068$	–	$0.249 \pm0.041$	$0.026 \pm0.095$; 0.27
$M_{\mathrm{step}}=10^{10.329}\,M_\odot$
$\mu_{griz}$	Foundation+YSE	$0.135 \pm0.017$	$0.191 \pm0.043$	$0.115 \pm0.016$	$0.196 \pm0.041$	$-0.082 \pm0.027$; 3.04
	Foundation	$0.134 \pm0.017$	–	$0.115 \pm0.016$	–	$-0.080 \pm0.029$; 2.73
	YSE	–	$0.191 \pm0.045$	–	$0.197 \pm0.039$	$-0.101 \pm0.081$; 1.25
$\mu_{gri}$	Foundation+YSE	$0.136 \pm0.018$	$0.172 \pm0.043$	$0.120 \pm0.017$	$0.211 \pm0.045$	$-0.079 \pm0.028$; 2.81
	Foundation	$0.134 \pm0.018$	–	$0.121 \pm0.018$	–	$-0.078 \pm0.031$; 2.51
	YSE	–	$0.174 \pm0.047$	–	$0.216 \pm0.045$	$-0.087 \pm0.078$; 1.13
$\mu_g$	Foundation+YSE	$0.143 \pm0.022$	$0.137 \pm0.049$	$0.136 \pm0.023$	$0.223 \pm0.059$	$-0.065 \pm0.033$; 1.97
	Foundation	$0.142 \pm0.022$	–	$0.136 \pm0.022$	–	$-0.073 \pm0.036$; 2.06
	YSE	–	$0.145 \pm0.049$	–	$0.214 \pm0.056$	$-0.014 \pm0.081$; 0.18
$\mu_r$	Foundation+YSE	$0.122 \pm0.018$	$0.180 \pm0.046$	$0.147 \pm0.019$	$0.237 \pm0.051$	$-0.065 \pm0.030$; 2.14
	Foundation	$0.121 \pm0.018$	–	$0.145 \pm0.018$	–	$-0.076 \pm0.032$; 2.38
	YSE	–	$0.185 \pm0.047$	–	$0.222 \pm0.044$	$0.001 \pm0.080$; 0.01
$\mu_i$	Foundation+YSE	$0.144 \pm0.018$	$0.171 \pm0.042$	$0.111 \pm0.016$	$0.226 \pm0.047$	$-0.080 \pm0.028$; 2.89
	Foundation	$0.143 \pm0.019$	–	$0.110 \pm0.015$	–	$-0.076 \pm0.030$; 2.56
	YSE	–	$0.170 \pm0.043$	–	$0.228 \pm0.048$	$-0.076 \pm0.078$; 0.97
$\mu_z$	Foundation+YSE	$0.137 \pm0.019$	$0.195 \pm0.050$	$0.128 \pm0.018$	$0.275 \pm0.057$	$-0.074 \pm0.030$; 2.44
	Foundation	$0.136 \pm0.019$	–	$0.128 \pm0.018$	–	$-0.077 \pm0.032$; 2.41
	YSE	–	$0.189 \pm0.049$	–	$0.259 \pm0.058$	$-0.028 \pm0.094$; 0.29

In addition, we estimate the mass step in individual filters and find $\gamma_g = -0.099 \pm0.037$ mag, $\gamma_r = -0.102 \pm0.031$ mag, $\gamma_i = -0.131 \pm0.028$ mag, and $\gamma_z = -0.105 \pm0.031$ mag. The mass steps estimated in each wavelength regime are all consistent within $1\sigma$ of each other, and of the $gri z$ mass step. The residual scatter for each of these mass step estimates is also consistent within $1\sigma$ of that measured for the $gri z$ mass step in each mass bin. Therefore, we estimate a mass step at 3–4.5$\sigma$ significance that is consistent from the optical $gri$ bands to the $z$ band.

When considering a step at the median host galaxy mass of the sample, we find $\gamma_{griz} = -0.082 \pm0.027$ mag. Individually in each band we find $\gamma_g = -0.065 \pm0.033$ mag, $\gamma_r = -0.065 \pm0.030$ mag, $\gamma_i = -0.080 \pm0.028$ mag, and $\gamma_z = -0.074 \pm0.030$ mag. Again, the residual scatter terms are broadly consistent in trend with those estimated for $\mathcal{M}_{\mathrm{step}} = 10^{10}\mathcal{M}_{\odot}$. Our results are slightly less significant for this choice of mass step location. Still, there is evidence for a non-zero mass step from the optical $gri$ to the $z$ band that persists for both choices of mass step location. The results are also consistent when considering the Foundation and YSE subsets of the sample individually.

Figure 8 visualises the results by individual filter assuming a mass step location of $10^{10}\mathcal{M}_{\odot}$. Shown is the host galaxy mass versus Hubble residual for the Hubble residuals from the four Hubble diagrams constructed from $\mu_g$, $\mu_r$, $\mu_i$, and $\mu_z$. Recall that these four sets of Hubble residuals have been corrected for dust in post-processing assuming $A_V$ from the $gri z$ fit and $R_V = 2.61$. The fitted population mean Hubble residual in the low and high mass host bins, $\delta_{\mathrm{LM/HM}}$ is also plotted with $1\sigma$ uncertainties. In addition, the residual scatter in each mass bin and for each survey are reported in the figure.

We additionally fit for a mass step assuming different values of $R_V$ for SNe Ia in low mass host galaxies and for those in high mass host galaxies, rather than the fiducial assumption of a single $R_V$ for all SNe Ia. For the split $R_V$ case, we assume $R_{V,\mathrm{LM}} = 3.14$ for SNe Ia in low mass host galaxies and $R_{V,\mathrm{HM}} = 2.39$ for SNe Ia in high mass host galaxies, according to the results from Grayling et al. (2024). We use the results from the analysis in Grayling et al. (2024) that rely only on dust to explain differences in SNe Ia in high mass and low mass host galaxies. The results of this analysis are reported in Table 4. Figure 9 compares the posterior on the mass step parameters assuming a single $R_V$ to that assuming a split on $R_V$ based on host galaxy mass for a step location of $\mathcal{M}_{\mathrm{step}} = 10^{10}\mathcal{M}_{\odot}$. For this mass step location, we estimate $\gamma_z = -0.086 \pm0.032$ mag for the split $R_V$ case. For $\mathcal{M}_{\mathrm{step}} = 10^{10.329}\mathcal{M}_{\odot}$, we estimate $\gamma_z = -0.052 \pm0.030$ mag for the split $R_V$ case.

Fit	Sample	$\sigma_{\mathrm{res,LM}}$ [mag] Foundation	$\sigma_{\mathrm{res,LM}}$ [mag] YSE	$\sigma_{\mathrm{res,HM}}$ [mag] Foundation	$\sigma_{\mathrm{res,HM}}$ [mag] YSE	$\gamma$ [mag]	Significance
$M_{\mathrm{step}} = 10^{10}\,M_{\odot}$
$\mu_{griz}$	Foundation+YSE	$0.132 \pm0.022$	$0.092 \pm0.050$	$0.104 \pm0.013$	$0.217 \pm0.036$	$-0.102 \pm0.027$	$3.78$
	Foundation	$0.132 \pm0.022$	—	$0.105 \pm0.013$	—	$-0.106 \pm0.031$	$3.41$
	YSE	—	$0.105 \pm0.058$	—	$0.211 \pm0.036$	$-0.055 \pm0.074$	$0.74$
$\mu_{gri}$	Foundation+YSE	$0.138 \pm0.023$	$0.078 \pm0.048$	$0.109 \pm0.014$	$0.220 \pm0.038$	$-0.088 \pm0.029$	$3.03$
	Foundation	$0.137 \pm0.024$	—	$0.109 \pm0.014$	—	$-0.100 \pm0.033$	$3.06$
	YSE	—	$0.083 \pm0.051$	—	$0.215 \pm0.037$	$-0.015 \pm0.073$	$0.20$
$\mu_g$	Foundation+YSE	$0.158 \pm0.029$	$0.139 \pm0.074$	$0.120 \pm0.019$	$0.214 \pm0.047$	$-0.079 \pm0.036$	$2.16$
	Foundation	$0.158 \pm0.030$	—	$0.119 \pm0.018$	—	$-0.107 \pm0.038$	$2.80$
	YSE	—	$0.141 \pm0.075$	—	$0.196 \pm0.045$	$0.054 \pm0.093$	$0.58$
$\mu_r$	Foundation+YSE	$0.124 \pm0.022$	$0.138 \pm0.066$	$0.134 \pm0.016$	$0.235 \pm0.041$	$-0.072 \pm0.031$	$2.35$
	Foundation	$0.124 \pm0.021$	—	$0.134 \pm0.016$	—	$-0.088 \pm0.033$	$2.71$
	YSE	—	$0.156 \pm0.070$	—	$0.217 \pm0.040$	$0.046 \pm0.087$	$0.53$
$\mu_i$	Foundation+YSE	$0.151 \pm0.023$	$0.072 \pm0.045$	$0.103 \pm0.013$	$0.233 \pm0.039$	$-0.089 \pm0.028$	$3.16$
	Foundation	$0.152 \pm0.023$	—	$0.103 \pm0.013$	—	$-0.100 \pm0.033$	$3.04$
	YSE	—	$0.080 \pm0.050$	—	$0.224 \pm0.036$	$-0.006 \pm0.073$	$0.08$
$\mu_z$	Foundation+YSE	$0.148 \pm0.024$	$0.113 \pm0.061$	$0.116 \pm0.015$	$0.276 \pm0.048$	$-0.086 \pm0.032$	$2.73$
	Foundation	$0.148 \pm0.024$	—	$0.115 \pm0.015$	—	$-0.099 \pm0.033$	$2.96$
	YSE	—	$0.124 \pm0.069$	—	$0.246 \pm0.041$	$0.045 \pm0.093$	$0.49$
$M_{\mathrm{step}} = 10^{10.329}\,M_{\odot}$
$\mu_{griz}$	Foundation+YSE	$0.132 \pm0.018$	$0.202 \pm0.046$	$0.116 \pm0.016$	$0.201 \pm0.041$	$-0.055 \pm0.027$	$2.01$
	Foundation	$0.132 \pm0.017$	—	$0.116 \pm0.015$	—	$-0.054 \pm0.030$	$1.79$
	YSE	—	$0.202 \pm0.048$	—	$0.200 \pm0.041$	$-0.060 \pm0.082$	$0.74$
$\mu_{gri}$	Foundation+YSE	$0.141 \pm0.019$	$0.184 \pm0.046$	$0.120 \pm0.017$	$0.213 \pm0.045$	$-0.044 \pm0.029$	$1.54$
	Foundation	$0.140 \pm0.018$	—	$0.120 \pm0.017$	—	$-0.044 \pm0.032$	$1.40$
	YSE	—	$0.185 \pm0.047$	—	$0.215 \pm0.044$	$-0.047 \pm0.078$	$0.61$
$\mu_g$	Foundation+YSE	$0.143 \pm0.022$	$0.142 \pm0.050$	$0.138 \pm0.023$	$0.233 \pm0.059$	$-0.039 \pm0.033$	$1.18$
	Foundation	$0.142 \pm0.023$	—	$0.137 \pm0.021$	—	$-0.049 \pm0.037$	$1.32$
	YSE	—	$0.148 \pm0.051$	—	$0.221 \pm0.058$	$0.020 \pm0.080$	$0.25$
$\mu_r$	Foundation+YSE	$0.124 \pm0.018$	$0.196 \pm0.048$	$0.148 \pm0.019$	$0.243 \pm0.052$	$-0.030 \pm0.031$	$0.99$
	Foundation	$0.124 \pm0.019$	—	$0.147 \pm0.020$	—	$-0.041 \pm0.033$	$1.22$
	YSE	—	$0.200 \pm0.049$	—	$0.231 \pm0.049$	$0.044 \pm0.083$	$0.53$
$\mu_i$	Foundation+YSE	$0.147 \pm0.018$	$0.185 \pm0.044$	$0.111 \pm0.016$	$0.228 \pm0.048$	$-0.050 \pm0.028$	$1.81$
	Foundation	$0.147 \pm0.019$	—	$0.112 \pm0.016$	—	$-0.050 \pm0.031$	$1.59$
	YSE	—	$0.186 \pm0.044$	—	$0.223 \pm0.044$	$-0.034 \pm0.080$	$0.43$
$\mu_z$	Foundation+YSE	$0.139 \pm0.019$	$0.201 \pm0.050$	$0.128 \pm0.018$	$0.278 \pm0.058$	$-0.052 \pm0.030$	$1.71$
	Foundation	$0.138 \pm0.019$	—	$0.128 \pm0.018$	—	$-0.055 \pm0.033$	$1.69$
	YSE	—	$0.199 \pm0.050$	—	$0.260 \pm0.055$	$-0.000 \pm0.094$	$0.00$

Notably, the mass step in the $z$ band is reduced when assuming a split $R_V$ compared to the single $R_V$ case. The mass step in the optical is also reduced in this case compared to when assuming a single $R_V$. To contextualise this result, we highlight that the sample average extinction in the $z$ band when assuming $R_V=2.61$, $\langle A_{z,R_V=2.61}\rangle=0.096$ mag. For this choice of $R_V$, the $z$ band extinction is 2.60× less than the extinction in the $g$ band. When assuming $R_V=3.14$ in the low mass bin, the average $z$-band extinction of the objects in the low mass bin is larger by 0.014 mag than the average extinction in that bin when assuming $R_V=2.61$. For the high mass bin, the average $z$-band extinction when assuming $R_V=2.39$ decreases by 0.009 mag from the average in the high mass bin in the single $R_V$ case. These magnitude offsets explain the reduction in the size of the mass step when changing one’s assumption from the single $R_V$ to the split $R_V$ case, and suggests dust plays some role in the mass step. That we still estimate a consistent, non-zero mass step across wavelength in the split $R_V$ case, though, suggests that dust alone cannot explain the mass step.

Finally, we test the effect of not applying a dust correction to the individual-band distance estimates on the mass step inference. In this analysis, we are agnostic to any assumptions about the dust law shape in different host galaxies. Because of the large scatter in the dust-extinguished distances, the optical mass steps are not well constrained in this analysis variant. For the same reason, the $z$-band mass step is found to be of lower significance than under the single $R_V$ and split $R_V$ assumptions about dust. We compare the “No Dust Correction” posterior on the $z$-band mass step parameters to that of the single $R_V$ and split $R_V$ cases in Figure 9.

We visualise the inferred mass step as a function of wavelength for both the single and split $R_V$ cases in Figure 10. This figure illustrates the non-zero mass step that persists across wavelength regime, regardless of assumptions about $R_V$, which suggests that dust is not the dominant driver of the mass step. The reduced mass step estimated assuming different values of $R_V$ for SNe Ia in high mass vs low mass host galaxies relative to that estimated assuming a single $R_V$ for all objects, indicates that dust may play some role in explaining the mass step.

At present, the conclusions that can be drawn from this study are limited by the statistical power of the sample. In the future, larger samples of SNe Ia that are well-observed in the $z$ band will enable better constraints on the mass step in the $z$ band. These observational results should be compared to robust forward model simulations of a dust-driven mass step and of a mass step driven by intrinsic differences among SNe Ia to quantify the contribution of dust to the mass step. How to simulate SNe Ia light curves under different assumption about the driver(s) of the mass step is a non-trivial question that is left to future work.

Discussion

We now discuss the results of this analysis in relation to previous work investigating properties of SNe Ia across wavelength space. In § 5.1, we discuss the pre-correction residual scatter in the intrinsic magnitudes and their dependence on shape. In § 5.2, we compare the total scatter in the Hubble diagrams across wavelength space. Finally in § 5.3, we discuss previous estimates of the mass step in comparison to the mass step inferred in this work.

Standardisation of SNe Ia Peak $z$-band Magnitudes

We find that the residual scatter in the peak extinguished $z$-band magnitudes is smaller than in the $g$- and $r$-bands and larger than in the $i$-band. However, studies of the standardisation of SNe Ia in the NIR consistently find lower levels of residual scatter in the peak extinguished absolute magnitude in the NIR $YJH$ filters than this work reports in the $z$-band, regardless of differences in sample and methods (Dhawan et al. 2018; Stanishev et al. 2018; Avelino et al. 2019). While the $i$- and $z$-bands therefore do not appear to be more standard pre-correction for dust and light curve shape/colour effects compared to the NIR, they do represent a improvement from the optical $g$- and $r$-bands.

Additionally, we find a significant correlation between $M_{z,\mathrm{int}}$ and $\theta_{g r i}$, with a fitted linear model having a slope of $a = 0.143 \pm0.017$. Because different methods for light curve fitting use different definitions of a shape parameter, and this analysis is the first to use BayeSN to estimate the shape-luminosity relation, we can only compare qualitatively to the results from previous analyses. There is broad support for relatively weak shape-luminosity relations in the NIR (YJH) compared to in the optical (Kattner et al. 2012; Stanishev et al. 2018; Burns et al. 2018; Uddin et al. 2024; Peterson et al. 2024). This work affirms the presence of a shape-luminosity relation in the $z$-band, though understanding the strength of this correlation relative to optical and redder NIR filters will require data from SNe Ia observed from the optical through the NIR, including the $z$-band.

The Hubble Diagram

We find that the $z$-band Hubble diagram has a total scatter of 0.195 mag. Comparing to the total scatter in the other individual-band Hubble diagrams, we find that SNe Ia in the $z$ band are not relatively better standard candles than in the $i$ band. There are hints that SNe Ia in the $z$ band may be relatively better standard candles than in the $g$ and $r$ bands, but this result cannot be decisively claimed with the sample in this work. Notably, trends in the total and residual scatter across wavelength regimes differ for the Foundation and YSE subsets of the data, which we attribute to differences in the two surveys (e.g., being targeted vs untargeted).

We also find that adding the $z$-band data to construct a Hubble diagram from the full $g r i z$ data results in a decrease of 0.005 mag in the RMS compared to a Hubble diagram constructed from $g r i$ data only. This result is broadly consistent with previous work demonstrating the benefit of combining optical and NIR data to reduce scatter in the Hubble diagram (Avelino et al. 2019; Thorp et al. 2021; Mandel et al. 2022; Dhawan et al. 2023; Uddin et al. 2024).

Because of differences in samples and methods used to standardise SNe Ia across analyses, it is difficult to compare quantitatively for a determination of how standard SNe Ia in the $z$-band are relative to in the NIR. For such an analysis, a large uniform sample of SNe Ia observed from the optical to the NIR, including the $z$ band, is necessary. The upcoming Rubin-LSST, expected to begin survey operations in October 2025, will provide the necessary wavelength coverage and statistical leverage to draw stronger conclusions regarding the relative standardisation of SNe Ia across wavelength (Ivezić et al. 2019). Still, this work provides further evidence for the benefit of having data in redder filters (from the $z$-band to the NIR) in reducing the scatter in the Hubble residuals, as suggested by previous studies.

The Mass Step

We constrain the mass step to be $\gamma_{g r i z} = -0.126 \pm0.027$ mag for a mass step location of $10^{10}M_\odot$ in the fiducial analysis. Across choices of mass step location and differential treatments of dust for high and low mass hosts, this analysis finds evidence of a consistent mass step from the optical to the $z$-band at 1–4$\sigma$ significance. Previous analyses of the Foundation sample have yielded mass steps of similar size and significance as is reported in this analysis (Jones et al. 2019; Thorp et al. 2021).

Recently, DES, a five-year survey which observed the sky in the $g r i z$ filters with the Dark Energy Camera (DECam), produced a sample 1499 photometrically classified SNe Ia (Vincenzi et al. 2024). In their cosmological analysis, they correct for a mass step of size $-0.046 \pm0.009$ at $10^{10}M_\odot$. Also, Ginolin et al. (2025) recently analysed the mass step in a sample of 927 spectroscopically-classified SNe Ia observed in $g r i$ from the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Rigault et al. 2025). They find a much larger mass step of $-0.162 \pm0.020$ at $10^{10}M_\odot$ in a similar wavelength regime, though for a much lower redshift sample.

We largely attribute the difference seen across these analyses to the different samples. The selection function for each of these surveys will be slightly different based on whether the sample is targeted (i.e. Foundation, DES) vs. untargeted (i.e. YSE, ZTF), or spectroscopically (i.e. Foundation, YSE, ZTF) vs. photometrically classified (i.e. DES). These differences may affect the samples enough to result in the difference in the reported mass steps. In addition, Vincenzi et al. (2024) and Ginolin et al. (2025) use SALT3 and SALT2.4 models, respectively, to fit for the distance moduli (Kenworthy et al. 2021; Taylor et al. 2023), while this analysis uses BayeSN. These differences should be investigated in more detail to determine how they impact mass step inference (e.g., Grayling & Popovic 2024, for a comparison of SALT and BayeSN). Furthermore, Murakami & Scolnic (2025) recently raised that differences in analyses could lead to different results even in the same sample, so comparisons across works must be carefully considered.

In the NIR $Y$-, $J$-, and $H$-bands, mass steps of $-0.02$ to $-0.15$ mag have been reported for mass step locations of $10^{10-10.6}M_\odot$ (Uddin et al. 2020; Ponder et al. 2021; Johansson et al. 2021; Thorp & Mandel 2022; Jones et al. 2022; Uddin et al. 2024; Peterson et al. 2024). The approaches for estimating a NIR-only mass step have varied in their choices (e.g., whether they correct for a shape-luminosity relation, colour effects, or host galaxy dust; whether NIR filters are fit individually or jointly; choice of light curve fitting method).

Despite their differences, the majority of these works agree that there is evidence for a non-zero NIR mass step which is of a similar size as the optical mass step (Uddin et al. 2020; Ponder et al. 2021; Jones et al. 2022; Thorp & Mandel 2022). Based on this result, differences in dust properties in high- and low-mass host galaxies are not likely to be the physical cause for the mass step, as had been suggested by Brout & Scolnic (2021). The low statistical significance of the step size reported in these studies limits the strength of this claim. The NIR mass step found in these analyses are largely consistent within $1\sigma$ with the mass step reported in this work and with the size of the optical mass step found in each of their samples.

However, several works report mass steps which are consistent with zero in the NIR (Johansson et al. 2021; Uddin et al. 2024; Peterson et al. 2024). Johansson et al. (2021) and Uddin et al. (2024) report evidence for a mass step that is non-zero in the optical, but which is consistent with zero and the optical mass step in the NIR. Peterson et al. (2024) report a mass step that is consistent with zero in both the optical and NIR wavelength regimes. The size of the mass steps reported in these works are consistent with the size of the optical and $z$-band mass steps reported in this work despite differences in sample, light curve fitting models, and methods for correcting for dust.

As in this analysis of the optical and $z$ band mass steps, these previous studies of the mass step from the optical to the NIR largely do not constrain a change in the mass step as a function of wavelength, regardless of dust treatment. This result disfavors the hypothesis that dust explains the observed mass step in full given the current statistical precision. Though dust may play some role in explaining the mass step, as indicated by the lessened optical and $z$-band mass step for the split $R_V$ analysis in this work, dust alone does not appear to explain the mass step.

As previously mentioned, these results should be compared to those from robust forward simulations of the mass step under different assumptions about its origin. Recently, Peterson et al. (2024) used the SNANA package for cosmological analyses with SNe Ia (Kessler et al. 2009, 2019) to simulate the DESHVILS survey used in their analysis for two models of the mass step. Firstly, they assumed a dust-based mass step, à la Brout & Scolnic (2021), and secondly, a "achromatic" mass step, based on artificial variations in the SN Ia SED that result in a wavelength-independent scatter. Peterson et al. (2024) found their sample did not have sufficient precision to affirm or rule out dust as the origin of the mass step. The creation of realistic forward simulations and the statistical precision allowed by a sample of SNe Ia should be considered in such tests in the future.

The main limitation of current studies, though, is the availability of data in redder wavelength regimes. Studies of the mass step with some of the largest uniform samples of SNe Ia in the NIR have not been able to find a significant change in the size of the mass step from the optical to redder wavelength regimes. While this result provides some support for an intrinsic explanation of the mass step, it is possible that different $R_V$ values in high and low mass hosts produce a non-zero mass step in the NIR (Uddin et al. 2024; Peterson et al. 2024). Still, the NIR mass step would be lessened compared to the optical mass step in this case. Larger samples of SNe Ia observed across optical and NIR will be necessary to constrain small deviations of the mass step across wavelength space to definitively confirm or rule out a dust-based explanation for the mass step.

With larger samples of SNe Ia observed across a broad wavelength regime, the assumption of a mass “step” rather than another model for the Hubble residual dependence on host galaxy mass (e.g., a linear trend) should also be challenged. In this work, we consider only a step function because of the clear predictions made by Brout & Scolnic (2021) regarding the behaviour of a mass step as a function of wavelength, assuming the step is caused by differences in host galaxy dust properties. An alternative description of the relationship between supernova magnitude and host mass may better describe the physical process producing this observed residual dependence.

We proceed with a step in this work because a two-population model might suggest a step if mass traces the two populations cleanly (Briday et al. 2022). That said, a more complex model may be necessary if global host mass does not cleanly trace the two populations, which may arise if e.g., mass is determined to be a poor tracer or there exist more than two populations within the data. Indeed, Murakami & Scolnic (2025) recently presented a method for modelling a continuous evolution in the Hubble residuals with host properties, which they applied to the ZTF DR2 sample. In the future, more data across the optical and NIR will be necessary (Phillips et al. 2019) to determine if a different model is preferred over the fiducial step function model. Rubin-LSST is expected to provide an ideal sample for such further tests in its ten-year lifetime.

Finally, we note that global host galaxy mass is not the strongest tracer of the mass step (Briday et al. 2022; Wiseman et al. 2023). Indeed, other tracers – such as local mass or global/local star formation rate – may provide additional insight into how best to standardise SNe Ia (Sullivan et al. 2003; Gallagher et al. 2008; Childress et al. 2013b; Rigault et al. 2013; Moreno-Raya et al. 2016; Rose et al. 2019; Hakobyan et al. 2020; Rose et al. 2021; Kelsey et al. 2021; Millán-Irigoyen et al. 2022; Briday et al. 2022; Wiseman et al. 2023; Kelsey et al. 2023.

Conclusion

In this paper, we explore the cosmological utility of SNe Ia the $z$-band. Observations in the $z$-band benefit from a reduced sensitivity to dust compared to the optical, yet can be obtained to greater depth than NIR filters using CCDs. Therefore, the $z$-band sits at a sweet spot between the optical and NIR. Despite this, it has yet to be extensively studied. We present a fully Bayesian framework for analysing the properties of SNe Ia in the $z$-band, which can be applied to any cosmological analysis of SNe Ia and may be particularly useful to future studies of SNe Ia in an individual filter¹³.

We find that both the peak extinguished magnitudes and intrinsic magnitudes in the $z$-band are relatively better standard compared to the $g$-band and $r$-band, and exhibit a similar level of variance compared to the $i$-band. The residual scatter in the peak $z$-band extinguished magnitudes is $\sigma_{\mathrm{res,ext}}(z) = 0.25 \pm0.02$ mag, compared to $\sigma_{\mathrm{res,ext}}(g) = 0.37 \pm0.02$ mag, $\sigma_{\mathrm{res,ext}}(r) = 0.29 \pm0.02$ mag, and $\sigma_{\mathrm{res,ext}}(i) = 0.22 \pm0.01$ mag for the $g$-band, $r$-band and $i$-band, respectively. We also find that the $z$-band magnitudes at peak that have been corrected for extinction from dust in the host galaxy are correlated with the inferred light curve shape parameter. The linear shape-luminosity relation fit has a slope of $0.143 \pm0.017$ with a residual scatter around the relation of $0.184 \pm0.013$ mag. The $z$-band peak intrinsic magnitudes do, therefore, exhibit some dependence on the shape parameter.

We then construct the first $z$-band only Hubble diagram and find a total scatter of RMS = 0.195. In addition, the Hubble diagram shows less total scatter with optical and $z$-band data than when using only optical data. When combined with host galaxy mass information, we are able to constrain the size of the mass step in the $z$-band to be $\gamma_z = -0.105 \pm0.031$ mag for a step location of $M_{\mathrm{step}} = 10^{10}M_\odot$ when assuming $R_V = 2.61$ for all SNe Ia. This mass step is consistent with the mass steps in each of the $g$, $r$, and $i$ bands. When assuming different $R_V$ values for SNe Ia in high mass and low mass host galaxies, we estimate $\gamma_z = -0.086 \pm0.032$ for the same choice of step location. While the split $R_V$ case leads to a reduced mass step compared to the single $R_V$ case, the optical mass step is consistent with the $z$-band mass step to within $1\sigma$ regardless of assumptions about $R_V$. Therefore, this work suggests dust alone cannot explain the mass step.

While SNe Ia in the $z$-band do not appear to be standard candles, based on their dependence on the shape parameter, this analysis demonstrates that the $z$-band wavelength regime exhibits complementary properties to the optical. Observations in the $z$-band appear to improve distances constraints to SNe Ia, as the $griz$ Hubble diagram exhibits a smaller amount of total scatter than the optical only $gri$ Hubble diagram. Still, SNe Ia in the NIR do appear to be relatively better standard candles than SNe Ia in the $z$-band, as demonstrated by a higher total scatter in the $z$-band Hubble residuals compared to that in NIR Hubble residuals from other works. In the future, observations which span the optical, the $z$-band, and the NIR for the same sample of SNe Ia will be necessary to understand the relative standardisation properties of SNe Ia across these wavelength regimes.

As YSE is an ongoing survey, it will continue to discover and make

follow-up observations of SNe Ia, so the sample of SNe Ia observed in the $z$-band is growing by the day. A larger sample size will greatly improve the precision of studies of the standardisation properties of SNe Ia in the $z$-band and potentially allow for a more definitive determination of the contribution of dust to the post-standardisation dependence of the Hubble residuals on host galaxy mass. Furthermore, there is a sample of over 1,000 photometrically confirmed SNe Ia included in DR1 (Aleo et al. 2023). Although subject to contamination from non-IaS, the additional objects may prove beneficial to a study of the mass step in the $z$-band regardless. We leave the analysis of this sample and future data releases from YSE to future work.

More broadly, the $z$-band is increasingly important to understand for the modelling of the SNe Ia that will be observed by upcoming cosmological surveys. For example, Rubin-LSST will observe SNe Ia in the rest-frame $z$-band out to $z \sim0.1$ (Ivezić et al. 2019). This low-redshift sample of SNe Ia will serve as the anchor for crucial cosmological analyses in the future. Additionally, the Roman Space Telescope will observe SNe Ia out to $z \sim1.8$ (Hounsell et al. 2018), enabling precision studies of SNe Ia in redder wavelength regimes. Therefore, it is important that our models for the light curves, such as BayeSN, are trained on sufficient $z$-band data so that the behaviour of SNe Ia in this wavelength regime is well-constrained to take full advantage of the sample of SNe Ia from these upcoming instruments. Accurate and comprehensive models of SNe Ia trained across a wide wavelength regime and large samples of SNe Ia from next-generation cosmological surveys will enable estimates of the Hubble parameter and the dark energy equation of state to unprecedented precision in the next decade.

Acknowledgements

We thank J. Johansson for helpful discussions.

Supernova and astrostatisics research at Cambridge University is supported by the European Union’s Horizon 2020 research and innovation programme under European Research Council Grant Agreement No. 101002652 (PI K. Mandel) and Marie Skłodowska-Curie Grant Agreement No. 873089.

E.E.H. is supported by a Gates Cambridge Scholarship (#OPP1144). S.D. is supported by a Marie Skłodowska Curie Grant Agreement No. 890695, a Kavli Fellowship and a Junior Research Fellowship at Lucy Cavendish College. This work was funded by UK Research and Innovation (UKRI) under the UK government’s Horizon Europe funding Guarantee EP/Z000475/1. D.O.J. acknowledges support from NSF grants AST-2407632 and AST-2429450, NASA grant 80NSSC24M0023, and HST/JWST grants HST-GO-17128, HST-GO-16269, and JWST-GO-05324, awarded by the Space Telescope Science Institute (STScI), which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS5-26555. The UCSC team is supported in part by NASA grants 80NSSC18K0303, 80NSSC18K0303, 80NSSC19K1386, 80NSSC19K0113, 80NSSC20K0953, 80NSSC21K2076, 80NSSC20K0953, 80NSSC22K1513, 80NSSC23K0301, 80NSSC23K1151, and 80NSSC24K1411; NSF grant AST-1815935, the Gordon & Betty Moore Foundation, the Heising-Simons Foundation, and by a fellowship from the David and Lucile Packard Foundation to R.J.F. S.T. was supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programmes (grant agreement no. 101018897 CosmicExplorer). S.M.W. was supported by the UK Science and Technology Facilities Council (STFC). D.L. was supported by research grants (VIL16599, VIL54489) from VILLUM FONDEN. K.M.dS. acknowledges support by the NSF through grant AST-2108676 and from LSSTC through grant 2023-SFF-LFI-02-Villar. KMdS thanks the LSST-DA Data Science Fellowship Program, which is funded by LSST-DA, the Brinson Foundation, and the Moore Foundation. NE acknowledges support from NSF grant AST-2206164. G.N. gratefully acknowledges NSF support from AST-2206195, and CAREER grant AST-2239364, supported in-part by funding from Charles Simonyi, and OAC-2311355, DOE support through the Department of Physics at the University of Illinois, Urbana-Champaign (#13771275), and support from the HST Guest Observer Program through HST-GO-16764 and HST-GO-17128 (PI: R. Foley). C.G. is supported by research grants (00025501, VIL69896) from VILLUM FONDEN.

The Young Supernova Experiment (YSE) and its research infrastructure is supported by the European Research Council under the European Union’s Horizon 2020 research and innovation programme (ERC Grant Agreement 101002652, PI K. Mandel), the Heising-Simons Foundation (2018-0913, PI R. Foley; 2018-0911, PI R. Margutti), NASA (NNG17PX03C, PI R. Foley), NSF (AST-1720756, AST-1815935, PI R. Foley; AST-1909796, AST-1944985, PI R. Margutti), the David & Lucille Packard Foundation (PI R. Foley), VILLUM FONDEN (project 16599, PI J. Hjorth), and the Center for AstroPhysical Surveys (CAPS) at the National Center for Supercomputing Applications (NCSA) and the University of Illinois Urbana-Champaign.

Pan-STARRS is a project of the Institute for Astronomy of the University of Hawaii, and is supported by the NASA SSO Near Earth Observation Program under grants 80NSSC18K0971, NNX14AM74G, NNX12AR65G, NNX13AQ47G, NNX08AR22G, 80NSSC21K1572 and by the State of Hawaii. The PS1 Surveys and the PS1 public science archive have been made possible through contributions by the Institute for Astronomy, the University of Hawaii, the Pan-STARRS Project Office, the Max-Planck Society and its participating institutes, the Max Planck Institute for Astronomy, Heidelberg and the Max Planck Institute for Extraterrestrial Physics, Garching, The Johns Hopkins University, Durham University, the University of Edinburgh, the Queen’s University Belfast, the Harvard-Smithsonian Center for Astrophysics, the Las Cumbres Observatory Global Telescope Network Incorporated, the National Central University of Taiwan, STScI, NASA under grant NNX08AR22G issued through the Planetary Science Division of the NASA Science Mission Directorate, NSF grant AST-1238877, the University of Maryland, Eötvös Lorand University (ELTE), the Los Alamos National Laboratory, and the Gordon and Betty Moore Foundation.

A major upgrade of the Kast spectrograph on the Shane 3 m telescope at Lick Observatory was made possible through generous gifts from the Heising-Simons Foundation as well as William and Marina Kast. Research at Lick Observatory is partially supported by a generous gift from Google.

The data presented here were obtained (in part) with ALFOSC, which is provided by the Instituto de Astrofísica de Andalucía (IAA) under a joint agreement with the University of Copenhagen and NOT.

YSE-PZ was developed by the UC Santa Cruz Transients Team. The UCSC team is supported in part by NASA grants NNG17PX03C, 80NSSC18K0303, 80NSSC19K0113, 80NSSC19K1386, 80NSSC20K0953, 80NSSC21K2076, 80NSSC22K1513, 80NSSC22K1518, and 80NSSC23K0301; NSF grants AST–1720756, AST–1815935, and AST–1911206; grants associated with Hubble Space Telescope programs DD–14925, DD–15600, GO–15876, GO–16238, SNAP–16239, GO–16690, SNAP–16691, and GO–17128; the Gordon & Betty Moore Foundation; the Heising-Simons Foundation; fellowships from the Alfred P. Sloan Foundation and the David and Lucile Packard Foundation to RJF; Gordon and Betty Moore Foundation postdoctoral fellowships and a NASA Einstein fellowship, as administered through the NASA Hubble Fellowship program and grant HST-HF2-51462.001, to D.O.J.; and a National Science Foundation Graduate Research Fellowship, administered through grant No. DGE-1339067, to D.A.C.

Data Availability

The code and data needed to reproduce the analysis presented in this work is made publicly available on Github at https://github.com/erinhay/sn ia-zband/tree/main. The Foundation Data Release 1 photometry from Foley et al. (2018) are publicly available from the Monthly Notices of the Royal Astronomical Society at https://doi.org/10.1093/mnras/stx3136. The Young Supernova Experiment Data Release 1 from Aleo et al. (2023) are publicly available on Zenodo at https://doi.org/10.5281/zenodo.7314746.

Data Tables

In Table A1, we present an abridged version of the data used in this analysis. The table includes metadata for each SN Ia, including redshift and host galaxy mass, as well as fitted parameters, including the BayeSN parameters from each of the $griz$, $gri$, and individual-band fits and the peak magnitudes in each filter. The full table can be found online.

Sample Selection

Here, we elaborate on the motivation and justification for the $z$-band data requirement, as described in § 2.2. Because we fit the $z$-band only light curves, we need to ensure the SNe Ia are well enough sampled in the $z$ band to obtain a reliable distance estimate. Therefore, we set up a simple test to determine the amount of $z$-band data necessary to get a reliable distance estimate, which we describe below.

We simulate 80-100 light curves from the BayeSN model with $N$ $z$-band observations, where $N \in1, 2, ..., 7$. The times of observations are randomly sampled from between -5 and 15 days. We then fit these 661 unique light curves with BayeSN. Binning by number of $z$-band observations, we examine how many $2\sigma$ and $3\sigma$ outliers there are in recovered distance modulus. Regardless of number of $z$-band observations, the percent of $2\sigma$ outliers lies between 5% and 9% and does not show trends with number of $z$-band observations. Similarly, the $3\sigma$ outliers is low (between 0-1.5%) and consistent across number of observations in the light curve.

As the quantity of outliers is well calibrated across any number of observations in the light curve, we examine the distribution of uncertainty on the distance modulus for the 80-100 light curves per bin. The mean uncertainty consistently decreases from 0.22 for 1 $z$-band observation to below 0.18 for more than 5 $z$-band observations. From this metric, there is not an obvious choice of number of $z$-band observations to require, though the mean uncertainty drops below 0.2 for 3 $z$-band observations. The standard deviation of the distribution, however, drops from 0.026 for 2 $z$-band observations to 0.017 for 3 $z$-band observations – the biggest reduction when moving from 1 $z$-band observation to 7 $z$-band observations. Therefore, we proceed with the requirement of 3 $z$-band observations. Additionally, we test the results of this analysis when making a more strict cut on the amount of $z$-band data to confirm the robustness of the results. When we require 4 $z$-band observations, the sample is reduced to 124 SNe Ia, with 25 from YSE and 99 from Foundation. We find that the conclusions of the analysis remain unchanged. In particular, we find the $griz$ mass step to be $-0.064 \pm0.028$ mag and the $z$-band mass step to be $-0.061 \pm0.031$ mag at $10^{10}M_{\odot}$. We do not include an uncertainty arising from host galaxy mass estimate uncertainties here, though we do not expect it to be significant given it is a sub-dominant source of uncertainty in the main analysis. These results are of reduced significance than the fiducial analysis, but broadly align with the range of mass step estimates inferred in this work. The conclusions of the analysis remain unchanged.

In the GitHub repository associated with this paper, at https://github.com/erinhay/snia-zband, we include example notebooks for running the analysis which allow the user to change the sample according to their own data cuts. To test additional iterations of this analysis, we encourage the use of this functionality.

Alternative Fitting Methods

Fixed Theta Fits

Because we use only a subset of the data in the $z$-band only fit, the constraint on $\theta$ is not going to be the most accurate. Here, we explore the impact of fixing $\theta$ to the value determined from the fit to the optical $gri$ data, which we will denote $\theta_{gri}$, when fitting with $z$-band only data. We use $\theta$ from the $gri$ fit, rather than from the $griz$ fit, to avoid double-counting the $z$-band data.

As expected, the $z$-band Hubble residuals from the fixed $\theta$ BayeSN fits show less total scatter than those from the fit where $\theta$ is allowed to vary freely. A reduction of the Hubble residuals RMS from 0.195 to 0.190 is seen overall, with the RMS for the Foundation subset being 0.156 and for the YSE subset being 0.252. We also see a more significant detection of the mass step in the $z$ band with the fixed $\theta$ fits – at $\gamma_{z,\mathrm{fixed}\ \theta} = -0.145 \pm0.029$ mag for a step at $10^{10}M_{\odot}$. We do not include uncertainty arising from host galaxy mass estimate uncertainties here, though we do not expect it to be significant. This result leaves the conclusions of this work unchanged.

Peak Magnitudes from $griz$ Data

Alternatively to retrieving the $z$-band peak apparent magnitude from the $z$-band BayeSN fit, we can extract the peak apparent magnitude in the $z$ band from the fit to the full $griz$ light curve, as described in § 4.2. This approach relies more heavily on the assumed light curve model because the trained relation between the shape of the light curve in the $gri$ filters and the shape in the $z$ band will influence the resulting peak magnitudes. However, it allows us to leverage the multi-band data to get the best constraints on $\theta$ and $A_V$.

When using the $griz$ fit, the fitted mean peak intrinsic absolute magnitude in the $z$ band is found to be $M^{0}_{z,\mathrm{int}} = -18.31 \pm0.02$ mag, with a residual scatter of $\sigma_{\mathrm{res,int}}(z) = 0.24 \pm0.02$. Under this alternative peak $z$-band magnitude extraction routine, the results of the analysis remain consistent with the fiducial analysis (c.f. Table 2).

SNID	RA	DEC	Survey	$z_{\mathrm{final}}$	$\sigma_z$	$\log_{10}(M_\ast/M_\odot)$	$\mu_{g r i z}$	$A_{V,g r i z}$	$\theta_{g r i z}$	$T_{0,g r i z}$	Reduced $\chi^2_{g r i z}$
2016W	37.6654873	42.2359685	Foundation	0.019	0.00006700	11.008	34.477 ± 0.047	0.223 ± 0.056	1.275 ± 0.144	57417.927 ± 0.198	1.200
2016afk	155.6245508	15.0966103	Foundation	0.048	0.00000670	9.782	36.603 ± 0.063	0.236 ± 0.072	−0.075 ± 0.189	57439.853 ± 0.491	1.136
2020adkc	154.081739	−1.516317	YSE	0.092	0.00021300	10.548	37.680 ± 0.076	0.115 ± 0.076	−0.886 ± 0.238	59218.175 ± 0.569	0.761
2020evu	238.12045	24.494119	YSE	0.059	0.00000659	10.261	36.891 ± 0.077	0.450 ± 0.076	−0.646 ± 0.183	58941.498 ± 0.389	0.424
2020ewx	229.327104	13.02275	YSE	0.061	0.01000000	9.556	37.053 ± 0.059	0.143 ± 0.086	−0.830 ± 0.256	58944.728 ± 0.908	0.331
2020fhs	240.276175	19.450186	YSE	0.015	0.00000678	10.364	34.315 ± 0.065	0.510 ± 0.076	2.100 ± 0.174	58949.045 ± 0.361	0.214
2020ilb	145.266392	34.73362	YSE	0.052	0.00001000	11.073	36.435 ± 0.065	0.557 ± 0.081	−0.115 ± 0.357	58980.837 ± 0.563	0.203
2020ioz	195.94825	20.902494	YSE	0.061	0.01000000	9.474	37.361 ± 0.045	0.067 ± 0.047	−0.981 ± 0.190	58976.940 ± 0.352	0.575
2020rii	249.016723	33.209832	YSE	0.086	0.00000905	10.591	37.687 ± 0.061	0.176 ± 0.066	0.574 ± 0.258	59077.994 ± 0.345	2.353
2020rmg	17.280127	−15.997169	YSE	0.038	0.00011300	11.214	36.133 ± 0.062	0.079 ± 0.052	−0.009 ± 0.174	59091.215 ± 0.369	1.140

SNID

$\mu_{g r i}$

$A_{V,g r i}$

$\theta_{g r i}$

$T_{0,g r i}$

$\mu_z$

$\theta_z$

Number of Obs. $g$

Number of Obs. $r$

Number of Obs. $i$

Number of Obs. $z$

Time of First Obs. $g$ [MJD]

Time of First Obs. $r$ [MJD]

Time of First Obs. $i$ [MJD]

Time of First Obs. $z$ [MJD]

2016W

34.449 ± 0.065

0.240 ± 0.069

1.349 ± 0.156

57417.906 ± 0.206

34.623 ± 0.065

1.000 ± 0.250

7

7

7

7

57410.29

57410.29

57410.29

57410.29

2016afk

36.554 ± 0.079

0.297 ± 0.089

−0.221 ± 0.224

57439.637 ± 0.571

36.779 ± 0.060

−0.169 ± 0.328

7

11

7

6

57437.40

57437.40

57437.40

57437.40

2020adkc

37.677 ± 0.079

0.101 ± 0.074

−0.829 ± 0.243

59218.212 ± 0.588

37.753 ± 0.129

−0.630 ± 0.789

4

7

4

3

59230.60

59208.55

59212.64

59208.54

2020evu

36.992 ± 0.141

0.367 ± 0.122

−0.704 ± 0.212

58941.395 ± 0.382

37.082 ± 0.088

−0.448 ± 0.463

8

7

4

6

58931.52

58931.52

58958.59

58941.54

2020ewx

36.998 ± 0.079

0.207 ± 0.113

−0.867 ± 0.276

58944.965 ± 0.917

37.135 ± 0.136

−0.086 ± 0.927

5

5

4

3

58955.46

58945.50

58942.49

58951.45

2020fhs

34.367 ± 0.098

0.401 ± …

…

…

…

…

…

…

…

…

…

…

…

…

SNID	Peak App. Magnitude [mag] (Dust Corrected) $g$	Peak App. Magnitude [mag] (Dust Corrected) $r$	Peak App. Magnitude [mag] (Dust Corrected) $i$	Peak App. Magnitude [mag] (Dust Corrected) $z$	Dust Extinction $g$	Dust Extinction $r$	Dust Extinction $i$	Dust Extinction $z$
2016W	15.701 ± 0.009	15.739 ± 0.019	15.422 ± 0.066	16.181 ± 0.042	0.271	0.179	0.129	0.104
2016afk	17.652 ± 0.022	17.663 ± 0.025	17.573 ± 0.063	18.934 ± 0.057	0.287	0.189	0.137	0.110
2020adkc	18.329 ± 0.092	17.995 ± 0.034	18.213 ± 0.042	19.247 ± 0.042	0.042	0.092	0.060	0.050
2020evu	17.711 ± 0.021	17.799 ± 0.083	17.826 ± 0.030	18.540 ± 0.029	0.287	0.362	0.261	0.210
2020ewx	17.857 ± 0.098	17.932 ± 0.029	18.001 ± 0.188	18.750 ± 0.078	0.174	0.115	0.083	0.067
2020fhs	17.369 ± 0.052	17.428 ± 0.050	17.496 ± 0.119	17.929 ± 0.046	0.608	0.447	0.323	0.260
2020ilb	17.853 ± 0.098	17.932 ± 0.029	17.826 ± 0.030	17.929 ± 0.078	0.678	0.474	0.323	0.260
2020ioz	18.048 ± 0.052	18.222 ± 0.030	18.173 ± 0.176	19.066 ± 0.074	0.081	0.054	0.039	0.031