Binning is Sinning: Redemption for Hubble Diagram Using Photometrically Classified Type Ia Supernovae

Following the discovery of cosmic acceleration using a few dozen Type Ia supernovae (SNe Ia) (Riess et al. 1998; Perlmutter et al. 1999), increasingly large SN Ia samples have been used to improve measurements of the dark energy equation of state parameter, w. While the most precise w measurements are based on spectroscopically confirmed samples, photometric imaging surveys have been discovering far more supernovae than spectroscopic resources can observe. Existing SN surveys include the Sloan Digital Sky Survey-II (SDSS), ⁵ Supernova Legacy Survey (SNLS), Panoramic Survey Telescope and Rapid Response System-1 (PS1), ⁶ and Dark Energy Survey (DES); ⁷ future wide-area surveys that will overwhelm spectroscopic resources include the Legacy Survey of Space and Time (LSST) ⁸ and the Nancy Grace Roman Space Telescope. ⁹

To make full use of these large SN Ia samples, photometric identification using broadband filters has been developed over the past decade. Photometric methods include template matching (Sako et al. 2011) and machine learning (Lochner et al. 2016; Möller & de Boissière 2020; Qu et al. 2021), and they determine the probability (P_Ia) for each event to be an SN Ia. A framework to incorporate the resulting P_Ia was developed to measure cosmological parameters. This framework, called "Bayesian Estimation Applied to Multiple Species" (BEAMS; Kunz et al. 2007; Hlozek et al. 2012), was first used in an SN-cosmology analysis for the PS1 photometric sample (Jones et al. 2018).

As part of the DES SN-cosmology analysis, BEAMS was extended to "BEAMS with Bias Corrections" (BBC; Kessler & Scolnic 2017; hereafter KS17), a fitting procedure designed to produce a Hubble diagram (HD) that is corrected for selection biases and for contamination from core-collapse SNe (SNCC) and peculiar SNe Ia. BBC has been used in the SN Ia cosmology analysis for spectroscopic samples from Pantheon (Scolnic et al. 2018), DES (DES Collaboration 2019), and Pantheon+ (Brout et al. 2022a). BBC has also been used on a photometric sample from PS1 (Jones et al. 2019) and to examine contamination biases for the photometric DES sample (Vincenzi et al. 2023).

The BBC fit is performed in redshift bins to determine nuisance parameters and SN Ia distances that are independent of cosmological parameters, which enables more flexible use of cosmology-fitting programs. BBC therefore produces both a redshift binned and unbinned HD. Previous analyses using BBC took advantage of the binned HD to reduce computation time in cosmology-fitting programs. However, Brout et al. (2021, hereafter BHS21) showed that while the statistical uncertainty is the same using a binned or unbinned HD, the systematic uncertainty is ∼1.5 times smaller using an unbinned HD and unbinned covariance matrix. The uncertainty reduction is from an effect known as self-calibration (Faccioli et al. 2011).

BHS21 demonstrated the uncertainty reduction using BBC with a spectroscopically confirmed sample. Here we expand the use of unbinned HDs to photometric samples where BEAMS is used. In anticipation of very large samples in future analyses, we also explore the possibility of reducing computation time with a smaller HD and covariance matrix and still benefit from self-calibration: a rebinned HD in the space of redshift, color, and stretch. This choice of variables is motivated by the color-dependent systematic explored in BHS21. While the unbinned approach is optimal, the rebinned approach may be useful for the many intermediate simulation tests prior to unblinding.

We validate the unbinning and rebinning methods using simulations of DES that include SNe Ia, SNCC, and peculiar SNe Ia. The simulation and analysis presented here are similar to those in Vincenzi et al. (2023), and all analysis software used in this analysis is publicly available. The software for simulations, light-curve fitting, BBC, and cosmology fitting is from the SuperNova ANAlysis package (SNANA; Kessler et al. 2009). ¹⁰ The photometric classification software is from SuperNNova (SNN; Möller & de Boissière 2020). ¹¹ For workflow orchestration we used Pippin (Hinton & Brout 2020). ¹²

The outline of this Letter is as follows. The SALT2 and BBC formalism is reviewed in Section 2. The unbinning and rebinning procedures are presented in Section 3. The validation analysis is described in Section 4, and the validation results are given in Section 5.

Using the SALT2 framework (Guy et al. 2010), BBC is a fitting procedure that delivers an HD corrected for selection effects and corrected for contamination. BBC incorporates three main features: (1) BEAMS (Kunz et al. 2007), (2) fitting in redshift bins to avoid dependence on cosmological parameters (Marriner et al. 2011), and (3) detailed simulation to correct distance biases (Kessler et al. 2019a).

For each event, a SALT2 light-curve fit determines the time of peak brightness (t₀), stretch (x₁), color parameter (c), and amplitude (x₀) with ${m}_{B}\equiv -2.5{\mathrm{log}}_{10}({x}_{0})$. Within BBC, the measured distance modulus is defined by the Tripp equation (Tripp 1998):

$$\begin{eqnarray}&&\mu ={m}_{B}+\alpha {x}_{1}-\beta c+{ \mathcal M }-{\rm{\Delta }}{\mu }_{\mathrm{bias}}\,,\end{eqnarray} \tag{ 1 }$$

where α and β are the stretch- and color-luminosity parameters, ${ \mathcal M }$ is a global offset, and Δμ_bias is a distance-bias correction for each event, μ − μ_true, determined from a large simulation. The Δμ_bias value for each event is evaluated by interpolating in a five-dimensional space of {z, x₁, c, α, β}. In Equation (1), there is an implicit SN index i for μ, m_B , x₁, c, and Δμ_bias; this index is suppressed for readability. To simplify this study, host–SN correlations have been ignored in Equation (1), and also in the simulations used for validation.

Following Section 5 of KS17 and making a few simplifications for this review, the BBC fit maximizes a likelihood of the form ${ \mathcal L }={\prod }_{i=1}^{N}{{ \mathcal L }}_{i}$, where ${{ \mathcal L }}_{i}$ for event i is

$$\begin{eqnarray}&&{{ \mathcal L }}_{i}={P}_{\mathrm{Ia},i}{D}_{\mathrm{Ia},i}+(1-{P}_{\mathrm{Ia},i}){D}_{\mathrm{CC},i}\,,\end{eqnarray} \tag{ 2 }$$

where P_Ia,i is the photometric classification probability for event i to be an SN Ia. The SN Ia component of ${{ \mathcal L }}_{i}$ is ${D}_{\mathrm{Ia}}\sim \exp [-{\chi }_{\mathrm{HR}}^{2}/2]$, where ${\chi }_{\mathrm{HR}}^{2}={\mathrm{HR}}^{2}/{\sigma }_{\mu }^{2}$, HR is a Hubble residual described below, and σ_μ is the uncertainty on μ in Equation (1) as shown in Equation (3) of KS17. The non–SN Ia (contamination) component, D_CC, is evaluated from a simulation.

To remove the dependence on cosmological parameters in the BBC fit, we follow Marriner et al. (2011) and define the Hubble residual for the ith SN (HR_i ) as

$$\begin{eqnarray}&&{\mathrm{HR}}_{i}\equiv {\mu }_{i}-[{\mu }_{\mathrm{ref}}({z}_{i},{{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}})+{M}_{\zeta }],\end{eqnarray} \tag{ 3 }$$

where ${\mu }_{\mathrm{ref}}={\mu }_{\mathrm{ref}}({z}_{i},{{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}})$ is a reference distance computed from redshift z_i and an arbitrary choice of reference cosmology parameters denoted by ${{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}$. Our choice for ${{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}$ is flatness and

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}\equiv \{{{\rm{\Omega }}}_{{\rm{M}}}=0.3,w=-1\}.\end{eqnarray} \tag{ 4 }$$

M_ζ are fitted distance offsets in redshift bins denoted by ζ. An important concept is that using a cosmological model for μ_ref in Equation (3) is a convenience, not a necessity. For example, μ_ref could be replaced with a polynomial function of redshift or any function that approximates the distance-redshift relation within each redshift bin.

The BBC fit determines α, β, γ, M_ζ , and an intrinsic scatter term (σ_int) added to the distance uncertainties (Section 3) that results in a reduced χ² of one. The final binned HD is obtained as follows. First, each binned redshift (z_ζ ) is computed from

$$\begin{eqnarray}&&{z}_{\zeta }={\mu }^{-1}[{\overline{{\mu }_{\mathrm{ref}}}}_{\zeta }],\end{eqnarray} \tag{ 5 }$$

where μ⁻¹ is an inverse-distance function that numerically determines redshift from the weighted average of μ_ref in redshift bin ζ. The weight for each event is ${\sigma }_{\mu }^{-2}$. Next, the BBC-fitted distance in each redshift bin (μ_ζ ) is

$$\begin{eqnarray}&&{\mu }_{\zeta }={\overline{{\mu }_{\mathrm{ref}}}}_{\zeta }+{M}_{\zeta }.\end{eqnarray} \tag{ 6 }$$

and the collection of {z_ζ , μ_ζ } is the binned HD corrected for selection effects and contamination. The uncertainty on μ_ζ is the BBC-fitted uncertainty for M_ζ . If a different choice of ${\boldsymbol{ \mathcal C }}$ is used for μ_ref, the fitted M_ζ will change but the μ_ζ remain the same. For spectroscopically confirmed samples with all ${{P}_{\mathrm{Ia}}}_{,i}=1$, the unbinned HD is the collection of {z_i , μ_i } where the μ_i are computed from Equation (1) using the BBC-fitted parameters and each distance uncertainty (σ_μ,i ) is computed from Equation (3) in KS17. This procedure is an approximation that we rigorously test (Section 5.1) with high-statistics simulations.

For an unbinned HD, we use the BBC-fitted parameters and compute the distances defined in Equation (1). The unbinned distance uncertainties (σ_μ,unbin,i ), however, are not the naively computed distance uncertainties (σ_μ,i ) for a spectroscopically confirmed sample. To determine σ_μ,unbin,i , we require that the weighted average uncertainty in each redshift bin is equal to σ_M,ζ , the BBC-fitted uncertainty on M_ζ ;

$$\begin{eqnarray}&&1/{\sigma }_{M,\zeta }^{2}=\displaystyle \sum _{i\in \zeta }1/{\sigma }_{\mu ,\mathrm{unbin},i}^{2}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&=\displaystyle \sum _{i\in \zeta }{{ \mathcal P }}_{{\rm{B}}(\mathrm{Ia}),i}/{\left[{S}_{\zeta }{\sigma }_{\mu ,i}\right]}^{2},\end{eqnarray} \tag{ 8 }$$

where i is the SN index within redshift bin ζ,

$$\begin{eqnarray}&&{{ \mathcal P }}_{{\rm{B}}(\mathrm{Ia}),i}=\displaystyle \frac{{P}_{\mathrm{Ia},i}{D}_{\mathrm{Ia},i}}{{P}_{\mathrm{Ia},i}{D}_{\mathrm{Ia},i}+(1-{P}_{\mathrm{Ia},i}){D}_{\mathrm{CC},i}}\end{eqnarray} \tag{ 9 }$$

is the BEAMS probability for event i to be an SN Ia, and S_ζ is a ζ-dependent uncertainty scale that is computed to satisfy Equation (8). We find that S_ζ is a few percent greater than 1 because of small correlations between the fitted parameters (α, β, M_ζ ). Equation (8) is an ad hoc assumption and does not have a rigorous derivation. From Equations (7)–(8), the unbinned distance uncertainty is

$$\begin{eqnarray}&&{\sigma }_{\mu ,\mathrm{unbin},i}={S}_{\zeta }{\sigma }_{\mu ,i}/\sqrt{{{ \mathcal P }}_{{\rm{B}}(\mathrm{Ia}),i}}.\end{eqnarray} \tag{ 10 }$$

As a crosscheck, the weighted average of the Hubble residuals (〈HR_ζ 〉) should be zero for each redshift bin:

$$\begin{eqnarray}&&\langle {\mathrm{HR}}_{\zeta }\rangle =\left[\displaystyle \sum _{i\in \zeta }{\mathrm{HR}}_{i}{W}_{i}\right]/\left[\displaystyle \sum _{i\in \zeta }{W}_{i}\right],\end{eqnarray} \tag{ 11 }$$

where ${W}_{i}={\sigma }_{\mu ,\mathrm{unbin},i}^{-2}$ and HR_i is defined in Equation (3).

To limit the size of the HD and still benefit from reduced systematics, we propose rebinning in the space of redshift, stretch, and color, denoted by ζ = {z, x₁, c}. The distance modulus in each 3D ζ cell is a weighted average of distances in the cell,

$$\begin{eqnarray}&&{\mu }_{{\boldsymbol{\zeta }}}=\displaystyle \sum _{i\in {\boldsymbol{\zeta }}}{\mu }_{i}{W}_{i}/\displaystyle \sum _{i\in {\boldsymbol{\zeta }}}{W}_{i},\end{eqnarray} \tag{ 12 }$$

and following Equation (7) the uncertainty on μ _ζ is

$$\begin{eqnarray}&&1/{\sigma }_{\mu ,{\boldsymbol{\zeta }}}^{2}=\displaystyle \sum _{i\in {\boldsymbol{\zeta }}}1/{\sigma }_{\mu ,\mathrm{unbin},i}^{2}.\end{eqnarray} \tag{ 13 }$$

We test the unbinning and rebinning procedure (Section 3) by analyzing 50 simulated data-sized samples that closely follow Vincenzi et al. (2023). Each simulation corresponds to the 5 yr DES photometric sample for events with an accurate spectroscopic redshift of the host galaxy, combined with a spectroscopically confirmed low-redshift (LOWZ) sample (z < 0.1). The LOWZ sample uses the cadence and signal-to-noise ratio (S/N) properties for the Carnegie Supernova Project (CSP), ¹³ Center for Astrophysics (CFA3, CFA4; Hicken et al. 2009, 2012), and Foundation Supernova Survey (Foley et al. 2018).

The simulated models include:

• 1.  
  SNe Ia generated from the SALT2 model in Guy et al. (2010) using trained model parameters from Betoule et al. (2014);
• 2.  
  SNCC generated from spectral energy distribution (SED) templates in Vincenzi et al. (2019);
• 3.  
  Peculiar SN Iax using the SED model from Kessler et al. (2019b) ¹⁴ and extinction correction from Vincenzi et al. (2021);
• 4.  
  Peculiar 91bg-like SNe Ia using the SED model from Kessler et al. (2019b); and
• 5.  
  ΛCDM with Ω_M = 0.311, Ω_Λ = 0.689, and w = −1.

All simulated events are analyzed as follows:

• 1.  
  Use the SALT2 light-curve fit for each event to determine {t₀, m_B , x₁, c}.
• 2.  
  Apply selection requirements (cuts):

  • (a)  
    At least two passbands with maximum S/N >5;
  • (b)  
    At least one observation before t₀;
  • (c)  
    At least one observation >10 days after t₀ (rest frame);
  • (d)  
    ∣x₁∣ < 3 and ∣c∣ < 0.3;
  • (e)  
    Fitted uncertainties σ_x1 < 1.0 and σ_t0 < 2 days;
  • (f)  
    SALT2 light-curve fit probability >0.001; and
  • (g)  
    Valid bias correction in the BBC fit (see below).

  For the 50 samples, the average number of events passing cuts is 1897 (1622, 172, and 103 for DES, Foundation, and LOWZ, respectively).
• 3.  
  Determine P_Ia using the "SuperNNova" photometric classification (Möller & de Boissière 2020) ¹⁵ based on recurrent neural networks.
• 4.  
  Use the BBC fit to determine redshift-binned HD corrected for selection effects and non–SN Ia contamination. We use 20 z-bins, with bin size proportional to (1 + z)³ so that there is finer z-binning at lower redshift.
• 5.  
  Create statistical+systematic covariance matrix as in Conley et al. (2011).
• 6.  
  Use methods from Section 3 to produce an unbinned HD and two rebinned HDs. The first rebinned HD has 2 stretch and 4 color bins (Rebin2x4), and the total number of HD bins is 20 × 2 × 4 = 160. The second rebinned HD has 4 stretch and 8 color bins (Rebin4x8) and a total of 640 HD bins.
• 7.  
  Perform cosmology fit using a fast minimization program that combines an SN Ia HD with a cosmic microwave background prior that uses an R-shift parameter computed from the same cosmology parameters as in the simulated samples. To match the constraining power from Planck Collaboration et al. (2020), the R-uncertainty is σ_R = 0.006. We fit for Ω_M and w using the wCDM model, and we also fit for Ω_M, w₀, and w_a using the w₀ w_a CDM model.
• 8.  
  For the w₀ w_a CDM model, the figure of merit (FoM) is computed based on the dark energy task force (DETF) definition in Albrecht et al. (2006),

  $$\begin{eqnarray}&&\mathrm{FoM}={\left[\sigma ({w}_{0})\times \sigma ({w}_{a})\times \sqrt{1-{\rho }^{2}}\right]}^{-1},\end{eqnarray} \tag{ 14 }$$

  where ρ is the reduced covariance between w₀ and w_a .

Here we consider 70 systematic uncertainties that include the following:

• 1.  
  For each of the 34 passbands, shift the zero-point using the uncertainty from Brout et al. (2022b).
• 2.  
  For each of the 34 passbands, shift the filter transmission wavelength using the uncertainty from Brout et al. (2022b).
• 3.  
  Correlated zero-point shift, 0.00714λ/micron, corresponding to the Hubble Space Telescope calibration uncertainty for primary reference C26202.
• 4.  
  Galactic extinction uncertainty is 4%.

For each of the 68 zero-point and wavelength systematics, the SALT2 model is retrained and the shift is propagated in the simulated data. M. Vincenzi et al. (2023, in preparation) present the complete set of systematic uncertainties that includes calibration covariances.

5.1. BBC Sensitivity to Reference Cosmology

We begin by evaluating the sensitivity of fitted BBC distances to the reference cosmology ${{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}$ defined in Equation (3). Using the wCDM model to vary ${{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}$, we vary Ω_M up to ±0.1 with fixed w = −1, and we vary w by up to ±0.2 with fixed Ω_M = 0.3. We define a sensitivity metric to be

$$\begin{eqnarray}&&{\mathrm{STD}}_{\mu }=\mathrm{STD}({\mu }_{i,{\boldsymbol{ \mathcal C }}}-{\mu }_{i,{{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}}),\end{eqnarray} \tag{ 15 }$$

where STD is the standard deviation, ${\mu }_{i,{{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}}$ are unbinned distances (Equation (1)) using ${{\boldsymbol{ \mathcal C }}}_{\mathrm{ref}}$ (Equation (4)) in the BBC fit, and ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$ are unbinned distances from using a different ${\boldsymbol{ \mathcal C }}$ in the BBC fit. Results for eight wCDM model variants are shown in Table 1, and we find STD_μ ∼ 10⁻⁴ mag, which is about 1000 times smaller than the intrinsic scatter.

Table 1. STD_μ for Different ${\boldsymbol{ \mathcal C }}$ Choices

${\boldsymbol{ \mathcal C }}$ Variant	STDμ × 104
ΩM = 0.20	1.0
ΩM = 0.25	1.7
ΩM = 0.35	0.4
ΩM = 0.40	0.7
w = −1.2	0.7
w = −1.1	0.5
w = −0.9	0.3
w = −0.8	0.6
p0 (${\mu }_{i,{\boldsymbol{ \mathcal C }}}=40$)	59
p3 a	6.8
p6 b	3.1

Notes.

^a Third-order polynomial fit to the ΛCDM model. ^b Sixth-order polynomial fit to the ΛCDM model.

Download table as:  ASCIITypeset image

The last three rows of Table 1 are based on a polynomial function of redshift for ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$ to illustrate the BBC performance with poorer ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$ estimates. A constant ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$ (p0) results in STD_μ that is more than 1 order of magnitude larger than for the wCDM models, but is still well below 0.01 mag and thus works surprisingly well for such a poor ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$ estimate. Using third- and sixth-order polynomial fits to the baseline ΛCDM model (p3 and p6) works much better than constant ${\mu }_{i,{\boldsymbol{ \mathcal C }}}$, but still not quite as well as for the wCDM models.

5.2. Uncertainty Scale and HR Check

For the Ia+CC samples, the uncertainty scale (S_ζ in Equations (8) and (10)) versus redshift bin is shown in the upper panel of Figure 1, averaged over the 50 data samples. The average S_ζ value is ∼1.01 at all redshifts, and is thus a small correction. The Hubble residual crosscheck (〈HR_ζ 〉 defined in Equation (11)) versus redshift is shown in the lower panel of Figure 1. The values are within 0.001 mag of the expected value of zero.

Zoom In Zoom Out Reset image size

5.3. Bias Results for wCDM

5.4. Bias Results for w₀ w_a CDM

5.5. Binning Option Consistency

5.6. Impact of ${{ \mathcal P }}_{{\rm{B}}(\mathrm{Ia})}$ Term

As a follow-up to the original binned Hubble diagram (HD) from BBC, we have developed methods to derive cosmological results from an unbinned HD and from a rebinned HD in the space of redshift, stretch, and color. Averaging analysis results from 50 simulated data samples, we find biases consistent with zero and bias constraints almost 1 order of magnitude smaller than the single-sample uncertainty. We also find that using an unbinned HD results in a reduced total uncertainty consistent with BHS21. This conclusion holds for both the wCDM and w₀ w_a CDM models, and we find the same results with or without photometric contamination.

Using a rebinned HD with two stretch and four color bins (Rebin2x4), we recover unbiased cosmology results and also benefit from the reduced uncertainty in the unbinned HD. Using more bins (four stretch and eight color bins), there is still no bias and the total uncertainty is similar to the unbinned case. With ∼2000 events in the DES unbinned HD, the rebinned cosmology-fitting speed is only a factor of few faster compared to the unbinned case. With anticipated future samples in the 10⁴–10⁵ range, the rebinned HD size does not increase and therefore the cosmology-fitting speed improvement will be much more significant.

While our unbiased results are encouraging, we note that Mitra et al. (2023) reported a significant cosmology bias using an unbinned HD from a simulated LSST data sample of pure SNe Ia. We therefore recommend repeating our bias tests on simulated data for future analyses.

R.K. is supported by DOE grant DE-SC0009924. P.A acknowledges parts of this research were carried out on the traditional lands of the Ngunnawal people. We pay our respects to their elders past, present, and emerging. P.A. was supported by an Australian Government Research Training Program (RTP) Scholarship. We acknowledge the University of Chicago's Research Computing Center for their support of this work.