[table]capposition=top

APS/123-QED

# Improved Treatment of Host-Galaxy Correlations in Cosmological Analyses With
Type Ia Supernovae

Brodie Popovic111Email<EMAIL_ADDRESS>Dillon
Brout2,32,3footnotemark: 2,3, Richard Kessler4,54,5footnotemark: 4,5, Dan
Scolnic111Email<EMAIL_ADDRESS>Lisa Lu 66footnotemark: 6
1Department of Physics, Duke University, Durham, NC, 27708, USA.  2 Center for
Astrophysics, Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  3 NASA Einstein Fellow  4Department of Astronomy and Astrophysics, The
University of Chicago, Chicago, IL 60637, USA.  5Kavli Institute for
Cosmological Physics, University of Chicago, Chicago, IL 60637, USA.
6Independent Scholar

###### Abstract

Improving the use of Type Ia supernovae (SNIa) as standard candles requires a
better approach to incorporate the relationship between SNIa and the
properties of their host galaxies. Using a spectroscopically-confirmed sample
of $\sim$1600 SNIa, we develop the first empirical model of underlying
populations for SNIa light-curve properties that includes their dependence on
host-galaxy stellar mass. These populations are important inputs to
simulations that are used to model selection effects and correct distance
biases within the BEAMS with Bias Correction (BBC) framework. Here we improve
BBC to also account for SNIa-host correlations, and we validate this technique
on simulated data samples. We recover the input relationship between SNIa
luminosity and host-galaxy stellar mass (the mass step, $\gamma$) to within
$0.004~{}$mags, which is a factor of 5 improvement over the previous method
that results in a $\gamma$-bias of ${\sim}0.02$. We adapt BBC for a novel
dust-based model of intrinsic brightness variations, which results in a
greatly reduced mass step for data ($\gamma=0.017\pm 0.008$), and for
simulations ($\gamma=0.006\pm 0.007$). Analysing simulated SNIa, the biases on
the dark energy equation-of-state, $w$, vary from $\Delta w=0.006(5)$ to
$0.010(5)$ with our new BBC method; these biases are significantly smaller
than the $0.02(5)$ $w$-bias using previous BBC methods that ignore SNIa-host
correlations.

##  1\. Introduction

The standardisation of measurements of Type Ia Supernovae (SNIa) led to the
discovery that the universe is expanding at an increasing rate (Riess et al.,
1998; Perlmutter et al., 1999). A possible cause of this expansion, called
‘dark energy’, remains an unsolved mystery in cosmology to this day. In the
intervening decades since the discovery of dark energy, SNIa cosmology has
grown from sample sizes of tens to thousands. When combined with constraints
from the Cosmic Microwave Background, the statistical uncertainties in
measurements of the dark energy equation-of-state parameter, $w$, are on the
order of $\sim 0.04$, (Scolnic et al., 2018; Jones et al., 2018; Brout et al.,
2019b). Therefore, understanding systematic uncertainties on the level of
$\sim 0.01$ is needed to better measure the nature of dark energy. To further
reduce systematic uncertainties, here we improve the treatment of the
correlation between supernova properties and host-galaxy properties.

Most SNIa cosmology analyses use the SALT2 (Guy et al., 2010) framework that
relies on two parameters to standardise the SNIa brightness: a colour ($c$)
describing the wavelength-dependent luminosity and a light-curve stretch
($x_{1}$) describing the luminosity dependence on light-curve duration. In
addition to these SALT2 luminosity correlations, studies have shown a
correlation between the standardised brightness of SNIa and properties of
their host-galaxy, such as mass (e.g. Sullivan et al., 2010), Star Formation
Rate (e.g. Uddin et al., 2017), as well as other properties (Rose et al.,
2020). While $c$ and $x_{1}$ are rigorously included in the SALT2 model, the
correlation with host properties is typically included as an ad-hoc correction
to the supernova brightness. Additionally, correlations between $c$ and host-
properties as well as $x_{1}$ and host-properties have been ignored when this
ad-hoc correction is applied. In this analysis, we focus on host-galaxy
stellar mass ($M_{\textrm{stellar}}$) as a representative host-galaxy
property.

To predict and correct for biases in the measurement of distance modulus
values from standardisation methods, analyses have increasingly relied on
simulations (Kessler et al., 2009a; Betoule et al., 2014; Scolnic et al.,
2018; Brout et al., 2019b). These simulations used the publicly available
SuperNova ANAlysis package (SNANA: Kessler et al., 2009b); a detailed
description of the simulations is given in Kessler et al. (2019). To simulate
SNIa samples, simulations rely on measurements of underlying populations of
$c$ and $x_{1}$. Early SNIa cosmology analyses tuned these populations based
on visual comparisons of fitted $c$ and $x_{1}$ distributions between data and
simulations and lacked rigorous procedure. This motivated Scolnic & Kessler,
2016 (SK16) to formalise a method for determining the underlying populations.
By utilizing large simulations, SK16 determined the underlying populations by
tracking how selection effects, noise and intrinsic scatter cause populations
to ‘migrate’ from underlying distributions to observed distributions. Here, we
improve on SK16 by determining underlying $c$ and $x_{1}$ populations as a
function of $M_{\textrm{stellar}}$. Our population fitting code is publicly
available111https://github.com/bap37/ParentPops, along with sample
inputs.222T.B.D.

Simulated samples can be used to determine redshift-dependent bias
corrections, as done in the analysis for SDSS-II (Kessler et al., 2009a) and
for the Joint Lightcurve Analysis (JLA; Betoule et al., 2014). Subtle $c$ and
$x_{1}$ biases were shown in SK16 for simulated distance moduli, suggesting
that 1-dimensional redshift-dependent bias corrections are not sufficient.
This motivated Kessler & Scolnic (2017) to develop the BEAMS with Bias
Corrections (BBC) method, which corrects for observed light-curve fit
parameters based on rigorous simulations with corrections computed in
5-dimensional parameter space (redshift, colour, stretch, colour-luminosity
relationship, and stretch-luminosity relationship). Hereafter, we refer to
this method as BBC5D, which has been used in Dark Energy Survey 3 year
analysis (DES3YR; Abbott et al., 2019; Brout et al., 2019b), analyses of the
Pan-STARRs sample (Scolnic et al., 2018; Jones et al., 2018), and in a re-
analysis of the Sloan Digital Sky Survey-II (SDSS) (Popovic et al., 2019). We
refer to the redshift-only correction as BBC1D. A summary of these methods is
presented in Table 1.

Using SNANA to simulate DES3YR with a known input value of $\gamma$, Smith et
al. (2020) showed that the BBC5D-recovery of $\gamma$ from an input value is
significantly biased when correlations between $c$/$x_{1}$ and
$M_{\textrm{stellar}}$ are included in the simulation. To address the $\gamma$
bias found in Smith et al. (2020), we improve the BBC formalism from BBC5D to
BBC7D by introducing two new dimensions in the bias corrections.

Method | Bias-Correction | Comment | Citation
---|---|---|---
| Dependence |  |
BBC1D | $z$ | Biascor grid of $z$: correct $\mu(z).$ | Marriner et al. (2011)
|  |  | Kessler & Scolnic (2017)
BBC5D | $\\{z,x_{1},c,\alpha,\beta\\}$ | Biascor grid of $\alpha,\beta,c,x_{1},z$: correct $m_{B}(z),c$ and $x_{1}$ | Kessler & Scolnic (2017)
BBC7D | $\\{z,x_{1},c,\alpha,\beta,\theta,M_{\textrm{stellar}}\\}$ | Biascor grid of $\alpha,\beta,\c{,}x_{1},z,\gamma,M_{\rm stellar}$: correct $m_{B}(z),c$ and $x_{1}$. | This Work
BBC-BS20 | $\\{z,x_{1},c,M_{\textrm{stellar}}\\}$ | Biascor based on BS20 model, correct $\mu(z)$ in grid of $c,x_{1},z,M_{\rm stellar}$ | This Work

Table 1Breakdown of different bias correcting methods.

Brout & Scolnic, 2020, hereafter BS20, present a new explanation for the mass
step $\gamma$ as a consequence of varying extinction ratios for high and low
mass host galaxies. BS20 extended the SALT2 formalism to include two sources
of colour variation: intrinsic and dust. This formalism, however, is
incompatible with BBC5D and BBC7D; therefore we adapt BBC to work with BS20
(BBC-BS20).

Here, we address the concerns about corrections for the mass step raised by
Smith et al. (2020), and confirm the origin of the mass step presented in
BS20. Section 2 provides an overview of the data used in this work. Section 3
is an overview of population-fitting techniques and selection effects. In
Section 4 we review the methodology introduced in previous studies and expand
on them to track evolution with host-galaxy properties. Section 5 details
changes to the BBC formalism, the efficacy of which is presented in Section 6.
Finally, discussion and conclusions are presented in Section 7.

##  2\. Data

In this analysis, we compile light-curves from 6 SNIa samples that have been
spectroscopically confirmed. The samples used here are from the Sloan Digital
Sky Survey-II (SDSS; York et al., 2000; Sako et al., 2018), the Supernova
Legacy Survey (SNLS; Sullivan et al., 2010), Pan-STARRs (PS1; Jones et al.,
2018; Scolnic et al., 2018), the Foundation Supernova Survey (Foley et al.,
2018) and the Dark Energy Survey (DES; Brout et al., 2019a). The low-redshift
(Low-z) supernovae include samples from the Carnegie Supernova Project (CSP)
(Stritzinger et al., 2011) and Harvard-Smithsonian Center for Astrophysics
(CfA3-4) (Hicken et al., 2009b, a, 2012).

Figure 1.— The redshift distribution of the combined Foundation, Low-z, SDSS,
DES, SNLS, and PS1 data sample.

These samples are flux-calibrated to the SuperCal system (Scolnic et al.,
2015). Host-galaxy masses are taken from past analyses, notably we use the
masses for SDSS, PS1, and SNLS provided in the Pantheon sample (Scolnic et
al., 2018). For Foundation, we use the masses provided by Jones et al. (2018).
For DES, we use updated masses provided by Smith et al. (2020) and Wiseman et
al. (2020). These analyses derive host-galaxy masses from SED fitting to
broadband photometry. To measure parent populations, we apply the following
selection requirements (cuts):

  * •

light-curve fit probability $P_{\textrm{fit}}>0.01$ that varies for each
survey

  * •

fitted SNIa colour $|c|<.3$

  * •

fitted SNIa colour uncertainty $\sigma_{c}<.2$

  * •

fitted SNIa stretch $|x_{1}|<3$

  * •

fitted SNIa stretch uncertainty $\sigma_{x1}<1$

  * •

fitted SNIa time-of-peak-brightness error $\sigma_{t0}<2$

  * •

at least five observations

  * •

at least one observation before peak brightness

  * •

at least one observation after peak brightness

However, for Low-z, we drop the $c$/$x_{1}$ uncertainty cuts to increase
statistics. A summary of the total SNIa for each sample is given in Table 2.
The resulting redshift distribution is shown in Figure 1. The $c$, $x_{1}$,
and $M_{\textrm{stellar}}$ distributions for our combined data sample are
presented in Figure 2, along with the mean values for $c$ and $x_{1}$ as a
function of $M_{\textrm{stellar}}$.

Sample | # SNIa
---|---
Low-z | 270
Found | 120
SDSS | 385
SNLS | 237
DES | 227
PS1 | 308
Total | 1547
Table 2Summary of Data Statistics

##  3\.  Simulations and Analysis

###  3.1. Simulations

Simulations of supernovae are needed to correct for biases arising from
inefficiencies and from the light-curve fitting process. Bias corrections are
needed to measure populations for stretch and colour, and to measure
distances. Therefore we generate simulations to assess the impact of
systematics on cosmological measurements. To this end, we use SNANA to
simulate realistic samples of SNIa with appropriate noise, observing
conditions, cadence, and detection efficiencies. We generate these simulations
with the Flat $\Lambda$-CDM cosmology from Planck Collaboration et al. (2016):
$\Omega_{\textrm{M}}=0.315$ (matter density at $z=0$) and $w=-1$.

Figure 2.— The distributions of $c$ and $x_{1}$ for the combined Low-z, DES,
PS1, Foundation, SNLS, and SDSS spectroscopic samples as a function of host-
galaxy mass. Circles are the data; purple represent individual data points and
black are the average value for that bin. The dash-dotted line is 0.

The SNANA simulation is presented in Kessler et al. (2019), and here we
provide a brief overview. The simulation begins by generating a Spectral
Energy Distribution and applying dimming from cosmic expansion, lensing,
galactic extinction, and redshift effects. Next, telescope-specific
information (PSFs, zero point, sky noise, cadence) is used to simulate
measured fluxes and uncertainties in each passband. Finally, selection
efficiencies from the detection pipeline and spectroscopic targeting are
applied to select events for analysis. We modify the spectroscopic efficiency
function, Eff${}_{\textrm{spec}}$, for the SNLS survey to better fit the data
(Appendix A1).

These simulations rely on previously derived frameworks that describe
brightness variations among the SNIa population, which we call intrinsic
scatter. Here, we use the G10 (Guy et al., 2010) and C11 (Chotard et al.,
2011) scatter models, translated into SED variants in Kessler et al. (2013),
along with the BS20 model. G10 and C11 are spectral variation models that
differ in the amount of variation ascribed to chromatic versus achromatic
scatter. For G10, roughly 70% of scatter is achromatic and 30% is chromatic,
whereas for C11, roughly 25% of scatter is achromatic and 75% is chromatic.
The BS20 model does not have an explicit SED variation model; instead it
includes dust parameter populations that are not part of the SALT2 framework.

To accurately characterise our systematic and statistical uncertainties, we
simulate 100 data-sized samples.

###  3.2. Analysis

We fit the light-curves using SNANA with the SALT2 model as developed in Guy
et al. (2010) and updated in Betoule et al. (2014). For each SNIa, the SALT2
fit gives four fitted parameters: $m_{B}$, the log of the fitted light-curve
amplitude $x_{0}$; $x_{1}$, the stretch parameter corresponding to light-curve
width,; $c$, the light-curve colour; and $t_{0}$, the time of peak brightness.
From the fitted SALT2 parameters, we infer distance modulus values with a
modified version of the Tripp distance estimator (Tripp, 1998). Following the
BBC formalism in Kessler & Scolnic (2017), the distance modulus ($\mu)$ is
defined as:

$\mu=m_{B}+\alpha x_{1}-\beta
c-M_{z_{i}}+\delta\mu_{\textrm{host}}+\delta\mu_{\textrm{bias}}$ (1)

where $\alpha$ and $\beta$ are global nuisance parameters relating to stretch
and colour respectively. $M_{z_{i}}$ is the distance offset in discrete
redshift bins denoted by $z_{i}$. $\delta\mu_{\textrm{host}}$ is the
luminosity correction for the mass step, and is defined as:

$\delta\mu_{\textrm{host}}=\gamma\times(1+e^{(X_{\textrm{stellar}}-S)/\tau_{X}})^{-1}-\frac{\gamma}{2}~{},$
(2)

where $\gamma$ is the magnitude of the SNIa luminosity difference between SNIa
in high and low mass galaxies, $X_{\textrm{stellar}}$ is
$\log_{10}$$M_{\textrm{stellar}}$, $S\sim 10$ is the step location, and
$\tau_{X}$ is the width of the step. Finally, $\delta\mu_{\textrm{bias}}$ is
the distance bias correction. As summarised in Table 1, for BBC1D,
$\delta\mu_{\textrm{bias}}$ is determined in bins of redshift using
simulations with a fixed $\alpha$ and $\beta$. For higher-dimensional
corrections (BBC5D/BBC7D/BBC-BS20), the corrections are applied in bins of
redshift along with other SNIa and host galaxy properties.

Furthermore, the total distance modulus error $\sigma^{2}$ is described as

$\sigma^{2}=\sigma^{2}_{N}+\sigma^{2}_{\textrm{lens}}+\sigma^{2}_{\sigma
z}+\sigma^{2}_{\textrm{int}}$ (3)

where $\sigma_{N}$ includes the uncertainties from SALT2 parameters and their
covariances (see Equation 3 in Kessler & Scolnic (2017)),
$\sigma_{\textrm{lens}}=0.055z$ is the uncertainty from weak gravitational
lensing (Scolnic et al., 2018), $\sigma_{\sigma z}$ is the uncertainty from
peculiar velocity (Carrick et al., 2015), and $\sigma_{\textrm{int}}$ is the
intrinsic scatter contribution determined in the BBC fitting process.

The BBC fit maximises a likelihood that includes separate terms for SNIa and
core collapse (CC) SNe. Because our sample includes only spectroscopically
confirmed events, the CC term is excluded. The BBC fit determines $\alpha$,
$\beta$, $\gamma$, $\sigma_{\textrm{int}}$ over the entire redshift range, and
a distance modulus in each redshift bin ($z_{i}$). Following Brout, Hinton &
Scolnic (2020), because we are not assessing the covariances of systematics, a
$z$-binned Hubble diagram results in equivalent cosmological constraints as
that of an unbinned set. These binned distances are combined with priors in
cosmological fits. Here, $\Omega_{\textrm{M}}$ and $w$ are obtained with
“wfit”, a $\chi^{2}$ minimisation program using MINUIT (James & Roos, 1975).
We use a Gaussian $\Omega_{\textrm{M}}$ prior with a mean of 0.315 and width
of $0.01$.

##  4\. Modeling SNIa Populations With Host Galaxy Correlations

Here, we determine the underlying correlations between SNIa parameters and
host-galaxy properties given different models of intrinsic scatter (G10, C11:
Section 2; BS20: Section 4.4). In Section 4.1, we review the methodology from
SK16, who developed an underlying parent population method to ensure that
resulting simulated distributions of $c$ and $x_{1}$ match the data. Naively,
we would extend SK16 to a three dimensional population that includes
$M_{\textrm{stellar}}$, however, the SK16 method relies on a robust model of
simulated uncertainties that has been validated for $c$ and $x_{1}$, but not
for $M_{\textrm{stellar}}$. Therefore, we use the SK16 method for $c$ and
$x_{1}$, and introduce a different method for $M_{\textrm{stellar}}$ which
does not rely on modeling $M_{\textrm{stellar}}$ uncertainties.

We fit our underlying parameter populations separately for Low-z, Foundation,
and the combined high-$z$ collection of DES, SDSS, SNLS, and PS1. As with
SK16, we find notable differences between low and high redshift surveys. The
difference in $M_{\textrm{stellar}}$ distributions and selection effects
between Low-z and Foundation (Figure 3 in Jones et al., 2019) motivate the
separate fits.

We derive the dependence of underlying parameter populations on host-galaxy
mass for the G10 and C11 models (4.1 and 4.2). In Section 4.3 and 4.4, we
review the BS20 framework and describe improvements to modeling the parent
populations.

###  4.1. Review of Method to Determine Uncorrelated Parent Populations

Following SK16, we extend their stretch population description for high-
redshift surveys (PS1, SDSS, DES, and SNLS) by replacing the asymmetric
Gaussian with an asymmetric generalised normal distribution that depends on
four parameters,

$P(x_{1})=\begin{cases}e^{(-|x_{1}-\overline{x}_{1}|^{n}/n\sigma_{-}^{n})}&\textrm{if}~{}x_{1}\leq\overline{x}_{1}\\\
e^{|-|x_{1}-\overline{x}_{1}|^{n}/n\sigma_{+}^{n})}&\textrm{if}~{}x_{1}>\overline{x}_{1}\end{cases}$
(4)

where $\overline{x}_{1}$ is the value at peak probability of the asymmetric
generalised normal distribution, $n$ is the shape, and $\sigma_{-}$ and
$\sigma_{+}$ are the width parameters for negative and positive $x_{1}$
values, respectively. The motivation for this distribution is discussed in
Section 6.1. For $n=2$, Equation 4 reduces to the three parameter asymmetric
Gaussian distribution used in SK16. This parametrisation is also used to
describe the SNIa colour $c$.

For Low-z and Foundation, the stretch distribution is double-peaked. We
therefore take an alternative approach following Scolnic et al. (2018) and use
a double Gaussian model,

$P(x_{1})=A_{1}\times
e^{(-|x_{1}-\overline{x_{1}}_{1}|^{2}/2\sigma_{1}^{2})}+A_{2}\times
e^{(-|x_{1}-\overline{x_{1}}_{2}|^{2}/2\sigma_{2}^{2})}$ (5)

where $A_{i}$ is the weight, $\overline{x_{1}}_{i}$ is the mean value, and
$\sigma_{i}$ is the standard deviation of the respective Gaussian.

A simulation with a flat distribution of true colour and true stretch is
generated, and SALT2 light-curve fits determine the measured colour and
measured stretch for the SNIa that pass light-curve quality cuts. The true and
measured values from these flat simulations are used to compute a migration
matrix ($X$ for stretch and $C$ for colour) that captures the migration from
the input distribution to the observed distribution through selection effects,
measurement noise and intrinsic scatter. Each component of the matrix,
$X_{ij}$, describes the likelihood that a true value in an input stretch bin
$i$ $(x_{1_{i}})$ migrates to a measured value in stretch bin $j$
$(x_{1_{j}})$, and similarly for colour.

For an underlying stretch population, $P(x_{1})$, we define a binned
distribution $\vec{P}_{x}$ with components $P_{xi}$. SK16 defines
$\vec{\Delta}_{x}$, the data-simulation difference vector, as

$\displaystyle\vec{\Delta}_{x_{1}}=\begin{bmatrix}o_{x1}\\\ o_{x2}\\\
\vdots\\\ o_{xn}\end{bmatrix}-\begin{bmatrix}X_{1,1}&X_{1,2}&\cdots&X_{1,n}\\\
X_{2,1}&X_{2,2}&\cdots&X_{2,n}\\\ \vdots&\vdots&\ddots&\vdots\\\
X_{d,1}&X_{d,2}&\cdots&X_{d,n}\end{bmatrix}\times\begin{bmatrix}P_{x1}\\\
P_{x2}\\\ \vdots\\\ P_{xn}\end{bmatrix}$ (6)

where the measured distribution vector $\vec{o}_{x}$ has the same binning as
$X$ and $\vec{P}_{x}$.

For the $\chi^{2}$ calculation, the associated data error vector is
$\vec{e}_{x}=[e_{x1},e_{x2},...,e_{xn}]$, where $e_{x_{i}}=\sqrt{o_{x_{i}}}$
for ${o_{x_{i}}}>0$ and $e_{x_{i}}=1$ for ${o_{x_{i}}}=0$; while not
technically correct for a Poisson distribution (Baker & Cousins, 1984), SK16
has shown this this error approximation is sufficient. The four parameters
that describe $P(x_{1})$ are determined by minimising the $\chi_{x}^{2}$
defined as

$\chi_{x}^{2}=\sum_{i=1}^{n}\left(\frac{\Delta_{xi}}{e_{xi}}\right)^{2}~{}.$
(7)

For the colour distribution, $\chi^{2}_{c}$ is defined similarly using
$\vec{\Delta}_{c}$ and $\vec{e}_{c}$.

SK16 performed a grid search for their parameters, where here we use a Monte
Carlo minimisation procedure using the emcee python package (Foreman-Mackey et
al., 2013). We compare our population parameters to those of SK16 and find
that we replicate their results to within $1\sigma$ for parameters describing
the colour and stretch distributions (Equation 4 with $n=2$). We find that
this assumption works well for the $c$ population, and therefore we fix $n=2$
(see Appendix A2). For $x_{1}$, however, we find that fixing $n=3$ works
better and the resulting $\chi_{x1}^{2}$ is smaller by $\sim 3$ compared to
$n=2$.

###  4.2. Extending Parent Populations to Include Mass Dependence

The ideal approach for including mass-dependent correlations is to replace
$\vec{o}_{x1}$ in Equation 6 with a 2-dimensional array of stretch and host-
galaxy mass, $\vec{o}_{x1,M}$, where the subscript is $M=$
$M_{\textrm{stellar}}$. Similarly, the migration matrix $X$ and probability
vector $\vec{P}_{x}$ would also be extended to include $M_{\textrm{stellar}}$.

However, the $M_{\textrm{stellar}}$ measurements lack a well defined
uncertainty and therefore the migration matrix is not as well determined in
the $M_{\textrm{stellar}}$ dimension as it is for $c$ or $x_{1}$. We therefore
assume that the measured $M_{\textrm{stellar}}$ is the true
$M_{\textrm{stellar}}$. To reduce the dependence on $M_{\textrm{stellar}}$
uncertainties, we implement a 2 step process. First, we fit the parent
populations by minimising $\chi_{x}^{2}$ and $\chi_{c}^{2}$ (Eq. 7) in
$M_{\textrm{stellar}}$ bins. Ideally we would perform this minimisation in
small $M_{\textrm{stellar}}$ bins, however, the statistics of SNIa per bin is
insufficient. As a compromise, we fit in relatively large
$M_{\textrm{stellar}}$ bins of $1.2\times 10^{10}~{}M_{\odot}$ but use a small
step size of $0.2\times 10^{10}~{}M_{\odot}$. Although the
$M_{\textrm{stellar}}$ bins are strongly correlated, we have used simulated
data samples to validate our method. In the second step, we re-weight the
simulated $M_{\textrm{stellar}}$ distribution to match the data.

###  4.3. REVIEW OF BS20

BS20 observed a significant ($>10\sigma$) colour-dependent Hubble scatter, and
a $5\sigma$ colour-dependence on the mass step ($\gamma$). They model this
effect with simulations using the SALT2 model combined with host galaxy dust.
They define 3 contributions to the observed colour,

$c_{\textrm{obs}}=c_{\textrm{int}}+E_{\textrm{dust}}+\epsilon_{\rm noise}$ (8)

where $c_{\textrm{int}}$ is the intrinsic SNIa colour, $E_{\textrm{dust}}$ is
the dust extinction, and $\epsilon_{\rm noise}$ is measurement noise. These
parameters contribute to a total change in true SNIa brightness of

$\Delta m_{B,\textrm{obs}}=\beta_{\textrm{SN}}\times
c_{\textrm{int}}+R_{V}\times E_{\textrm{dust}}$ (9)

where $\beta_{\textrm{SN}}$ is introduced as the correlation coefficient
between the intrinsic colour and SNIa luminosity, and $R_{V}$ is the dust
extinction ratio. Components of the BS20 model that are intrinsic to the SNIa,
i.e. $c_{\rm int}$ and $\beta_{\rm SN}$, are assumed to be independent of
$M_{\textrm{stellar}}$. Thus, BS20 determine the following parameters
describing the distributions of dust and colour: Gaussian distribution of
intrinsic color ($\bar{c}_{\rm int}$, $\sigma_{c\textrm{int}}$), Gaussian
distribution of intrinsic color-luminosity coefficient ($\bar{\beta}_{\rm SN}$
& $\sigma_{\beta_{\rm SN}}$), and Gaussian distribution of dust extinction
ratios ($\bar{R}_{V}$ & $\sigma_{Rv}$) for low and high mass, and exponential
distribution of the dust extinction ($\tau_{E}$) for low and high mass.

The BS20 parameters are determined in a forward-modeling fitting process by
minimising the $\chi_{\rm TOT}^{2}$ in Eq. 8 of BS20. The $\chi_{\rm TOT}^{2}$
includes constraints enforcing consistency between data and simulations for:
the colour distribution, the Hubble diagram scatter versus colour, the BBC-
fitted Hubble residual versus colour, and BBC-fitted $\beta$. BS20 fixes the
$\alpha$ in the simulations to be that obtained from the dataset using BBC1D,
and fix $\gamma=0$.

###  4.4. Upgrades to Parent Populations for Dust Based Scatter Models

For the BS20 model, we derive the $x_{1}$ populations following the method
discussed in Section 4.4. However, we keep the same colour distribution
presented in BS20 because their dust model parameterisation is incompatible
with the SK16 method. While this upgrade to the BS20 population parameters is
minimal, there are significant upgrades to the BBC formalism presented in
Section 5.3. These upgrades to the BBC formalism are essential because BS20 is
incompatible with higher dimensional BBC procedures. A new and significantly
improved formalism for determining the underlying population parameters for
BS20 is under development and will be presented in a future work.

##  5\. Correcting Distance Biases

Here we review the BBC methodology for correcting distance biases (Section
5.1) and describe improvements. In Section 5.2, we detail the improvements to
account for correlations between SNIa brightness and host-galaxy properties.
Finally, in Section 5.3, we describe improvements to the BBC formalism
required for bias corrections using the BS20 model.

###  5.1. Review of 5D Bias Corrections

A brief overview of the BBC process is described in Section 3.2 and here we
elaborate on the bias corrections component for BBC5D. Kessler & Scolnic
(2017) expanded on the 1D method by defining a bias-corrected distance where
the individual Tripp components of $m_{B}$, $c$, and $x_{1}$ are corrected.
They define a bias-corrected distance,

$\displaystyle\mu^{*}$ $\displaystyle=m_{B}^{*}+\alpha x_{1}^{*}-\beta
c^{*}-M_{z_{i}}$ (10)
$\displaystyle=(m_{B}-{\delta}_{m_{B}})+\alpha(x_{1}-{\delta}_{x_{1}})-\beta(c-{\delta}_{c})-M_{z_{i}}$
$\displaystyle=m_{B}+\alpha x_{1}-\beta
c-M_{z_{i}}-\delta\mu_{\textrm{bias}}(z,x_{1},c,\alpha,\beta)$

where bias-corrected quantities are denoted with a star superscript and

$\delta\mu_{\textrm{bias}}\equiv({\delta}_{m_{B}}+\alpha{\delta}_{x_{1}}-\beta{\delta}_{c})$
(11)

Measurement noise and intrinsic scatter preclude calculating the exact bias
correction for each event; therefore, $\delta\mu_{\textrm{bias}}$ is
interpolated in 5D cells of $\\{z,x_{1},c,\alpha,\beta\\}$. The
$\delta\mu_{\textrm{bias}}$ term is calculated with a large simulation
designated as a ‘BiasCor’, by comparing the observed values to the simulated
ones for each SNIa parameter. The first three dimensions ($z$,$x_{1}$,$c$) are
interpolated using the simulated populations. However, the luminosity
correlation coefficients ($\alpha$, $\beta$) are single valued and not
described with a population model akin to $z,x_{1},c$. Therefore, the bias
corrections for $\alpha$ and $\beta$ are determined on a $2\times 2$ grid that
brackets the values of $\alpha$ and $\beta$ found in the data. This grid
enables for an interpolation of the BiasCor sample at the value of the
proposed $\alpha$ and $\beta$ in each iteration of the BBC fit.

###  5.2. Improving Bias Corrections to Account for SNIa-Host Correlations

Smith et al. (2020) find that the fitted $\gamma$ from BBC5D is biased when
analysing simulated samples that include correlations between $c/x_{1}$ and
host-galaxy mass. To account for these biases in the mass step, we introduce
two new dimensions to the $\delta\mu_{\textrm{bias}}$ term in Equation 10:
$\theta$, which is a magnitude shift, and $M_{\textrm{stellar}}$.

The new $\theta$ dimension is incorporated into the BiasCor by adding a
magnitude shift of $+\theta$ to a random half of the simulation, and $-\theta$
to the other half; there is no correlation between $\theta$ and host
properties, and thus $\theta\neq\gamma$. Using BiasCor with $\theta$, the BBC
fit allows for a mag-shift ($\delta\mu_{\textrm{host}}$) as an arbitrary
function of $X=$ $\log_{10}$$M_{\textrm{stellar}}$ and other parameters. While
previous cosmology analyses have used $\delta\mu_{\textrm{host}}=-\gamma/2$
for $X>10^{10}$ and $\delta\mu_{\textrm{host}}=+\gamma/2$ for $X<10^{10}$,
here we adopt the more general $\delta\mu_{\textrm{host}}$ function with step
location $S$ and step width $\tau_{M}$ (Equation 2). At each step of the BBC
fit, the value of $\delta\mu_{\textrm{host}}$ for each event is used to
interpolate the BiasCor between $\pm\theta$ such that

$\displaystyle\delta\mu_{\textrm{bias}}={\delta\mu_{\textrm{bias}}}(\vec{x}_{5},X,-\theta)$
(12)
$\displaystyle+f\times[{\delta\mu_{\textrm{bias}}}(\vec{x}_{5},X,+\theta)-{\delta\mu_{\textrm{bias}}}(\vec{x}_{5},X,-\theta)]$

where $f=(\delta\mu_{\textrm{host}}+\theta)/2\theta$ and
$\vec{x}_{5}=\\{z,x_{1},c,\alpha,\beta\\}$. It is important that
$\theta>\delta\mu_{\textrm{bias}}$ to ensure a valid interpolation.

The addition of $\theta$ and $M_{\textrm{stellar}}$ results in changing the 5D
$\delta\mu_{\textrm{bias}}$ term in Equation 10 to a 7D
$\delta\mu_{\textrm{bias}}$ as follows,

$\displaystyle\delta\mu_{\textrm{bias}}(z,x_{1},c,\alpha,\beta)\to$ (13)
$\displaystyle\delta\mu_{\textrm{bias}}(z,x_{1},c,\alpha,\beta,\theta,M_{\textrm{stellar}})~{}.$

For the first six BiasCor dimensions, $\delta\mu_{\textrm{bias}}$ is
interpolated. For the $M_{\textrm{stellar}}$ dimension,
$\delta\mu_{\textrm{bias}}$ is evaluated in discrete bins to avoid
interpolating across the luminosity step at $10^{10}M_{\odot}$.

The effect of $\theta$ on the bias corrections is shown in Figure 3 as a
function of redshift. Here we set $\theta=\pm 0.06$, the effects of which is
negligible for $z<0.8$ and increases to $\sim 0.02$ mag at high-$z$. Because
$\theta$ is independent of supernova parameters, it can be used for
investigating correlations between any host galaxy property and SNIa
luminosity, not just host-galaxy stellar mass.

Figure 3.— Bias correction vs. redshift for $\theta=+0.06$ (blue) and
$\theta=-0.06$ (orange). The biased quantity is indicated on each panel.

###  5.3. Changes to the BBC formalism for BS20

While the SALT2 model is an accurate description of SNIa light curves, it is
nonetheless an approximation that ignores the difference between intrinsic
color variation and dust. Previous SNIa cosmology analyses have used the same
SALT2 model in both light curve fitting and BiasCor; this means that the
BiasCor corrects for selection effects, but does not correct for biases in the
SALT2 model. BS20 attempts to provide a more accurate light-curve model that
can be used for the BiasCor, and here we update BBC to be compatible with
BS20.

The BBC5D formalism is not compatible with BS20 for three reasons: 1) BBC
intrinsic scatter is characterised by a single colour-independent
$\sigma_{\textrm{int}}$ while BS20 uses a dust-dependent scatter that is
poorly characterised by a single $\sigma_{\textrm{int}}$, 2) BBC assumes the
single SALT2 colour-luminosity relation $\beta$, whereas BS20 uses
distributions for intrinsic $\beta_{\rm SN}$ and dust $R_{V}$ (Eq. 9), and 3)
$c_{\rm int}$ and $\beta_{\rm SN}$ in BS20 refer to intrinsic colour and
intrinsic colour-luminosity relationship, while SALT2 $c$ and $\beta$ include
both dust and intrinsic properties.

Ideally, to incorporate BS20, we would decompose the fitted colour $c$ into
$c_{\rm int}$ and $c_{\rm dust}$ and compute a true $\beta$ in the simulation.
While we plan to address these issues in a future work, they require
significant updates to both training and fitting code; instead, we make
several approximations. Here we describe BBC updates to be compatible with the
BS20 model:

  1. 1.

Replace the SALT2 model in the BiasCor with the BS20 model.

  2. 2.

Include $M_{\textrm{stellar}}$ as in Section 5.2, but not $\theta$ because the
mass step is predicted by the BS20 model without the $\gamma$ parameter.

  3. 3.

Remove grid-interpolation of $\alpha$ and $\beta$ because the BS20 simulations
are forward-modeled and result in BBC-fitted $\alpha$ and $\beta$ that are in
agreement with that of the data.

  4. 4.

$\delta\mu_{\textrm{bias}}$ is computed for distances instead of bias-
correcting each SALT2 parameter.

These changes are referred to as BBC-BS20 and result in a
$\delta\mu_{\textrm{bias}}$ dimensionality of
$\\{z,x_{1},c,M_{\textrm{stellar}}\\}$. For BBC-BS20, Equation 10 becomes:

$\mu^{*}=m_{B}+\alpha x_{1}-\beta c-M_{z_{i}}-\delta\mu_{\textrm{bias}}$ (14)

and Equation 11 becomes

$\delta\mu_{\textrm{bias}}=m_{B}^{\textrm{B}C}+\alpha
x_{1}^{\textrm{B}C}-\beta c^{\textrm{B}C}-M_{z_{i}}-\mu_{\textrm{true}}~{}.$
(15)

The true distance modulus for the BiasCor is $\mu_{\rm true}$, and the
superscript BC denotes that these values are from SALT2 light-curve fits to
the BiasCor events. We note that Equations 14 and 15 can be used with BBC5D
and the SALT2 model, and yield consistent results compared to using Equations
10 and 11.

##  6\. Results

Here we evaluate the accuracy of our population modelling and bias correction
approaches. Section 6.1 discusses the results of the parent population
modelling and how well our simulations match the data. Section 6.2 presents a
comparison of BBC7D to previous approaches for a variety of metrics. Section
6.3 contains the results of the new BBC-BS20 method for dust-scatter models. A
discussion of these results is in Section 7.

Figure 4.— The distributions of $c$ and $x_{1}$ for the combined DES, PS1,
SNLS, and SDSS spectroscopic samples as a function of host-galaxy mass.
Circles are the data; purple represent individual data points and black are
the average value for that bin. Solid lines are the G10 parent population
parameters: green is $\sigma_{+}$, blue is the peak value, and orange is
$\sigma_{-}$. The dash-dotted line is 0. The bin averages for the data are
presented in red-lined black dots, for the simulations, red-lined yellow
triangles. The histograms on the side show data in circles and simulated
results in dashed line. For both $c$ and $x_{1}$ we find good agreement
between the mean value in data and sims, with a $\chi^{2}/N=1.5$ and
$\chi^{2}/N=0.65$ respectively.  Figure 5.— The redshift evolution of $c$,
$x_{1}$, and host-galaxy stellar mass for the combined sample. The data is
presented in purple circles, and averaged bin values for the data are shown in
red-lined black circles. The averaged bin values for the SK16 simulations are
shown in red-lined blue squares, red-lined yellow triangles represent the
average bin values for this work. Individual simulated events are not shown
for clarity. We present the $\chi^{2}/N$ values for each parameter. SK16
assumes a flat distribution of $M_{\textrm{stellar}}$with redshift.

###  6.1. Parent Populations

Using the simulations described in Section 3.1, we follow the process
presented in Section 4.2 to determine the underlying populations of stretch
and colour when assuming the G10 or C11 scatter models. For both scatter
models, the population parameters for the generalised normal distribution,
$\sigma_{+}$, $\overline{x}_{1}$, $\sigma_{-}$, are presented in Appendix A2
for each individual survey. These population parameters are provided in steps
of $0.2\times 10^{10}M_{\textrm{stellar}}$.

In Figure 4, we show the comparison between sim and data for $c$ vs.
$M_{\textrm{stellar}}$ and for $x_{1}$ vs. $M_{\textrm{stellar}}$. We also
show the $M_{\textrm{stellar}}$ dependence of the underlying combined SDSS,
DES, SNLS, and PS1 parent population. The Foundation and Low-z parent
populations are shown separately in Figures 10 and 11 for colour and stretch
respectively. There are three notable results. First, there is excellent
agreement between the observed mean values for data and simulation. In the
combined SDSS, DES, SNLS, and PS1 sample, this parent population also
characterises the constituent individual surveys. Second, we find that both
the intrinsic $x_{1}$ and $c$ populations depend on $M_{\textrm{stellar}}$.

We evaluate the significance of $M_{\textrm{stellar}}$ dependence by comparing
each curve in Figure 4 to a null model that has no $M_{\textrm{stellar}}$
dependence. For the $x_{1}$ distribution, $\sigma_{+}$ is consistent with no
$M_{\textrm{stellar}}$ dependence with a confidence of 99.9%. The probability
that $\sigma_{-}$ has no $M_{\textrm{stellar}}$ dependence, however, is 1E-8.
This means that the observed $x_{1}$ dependence on mass is driven by
$\sigma_{-}$ increasing with increasing masses. Similarly, $\overline{x_{1}}$
has a low probability of being independent of $M_{\textrm{stellar}}$ at only
0.01%, but does not correlate with the observed distribution and is therefore
unlikely to be the primary driver. For the colour distribution, $\sigma_{-}$
and $\overline{c}$ are consistent with no $M_{\textrm{stellar}}$ dependence
with a confidence of 90%. Similarly to $x_{1}$, the probability that the
faint-side $\sigma_{+}$ has no $M_{\textrm{stellar}}$ dependence is 99.9%.
Taken together, this constitutes the third notable result: the
$M_{\textrm{stellar}}$ dependence of the observed $c$ and $x_{1}$
distributions is driven by increasing faint-side widths ($\sigma_{-}$ for
$x_{1}$ and $\sigma_{+}$ for $c$) rather than a shift in the mean values.

We also investigate the observed redshift dependence of our $c$ and $x_{1}$
distributions. Figure 5 shows that the observed colour and stretch
distributions depend on redshift for two reasons: 1) separate populations for
low and high $z$ (Section 4), and 2) selection effects. The data-simulation
agreement is reasonable, with $\chi^{2}/N=$ 4.8, 1.5, and 3.5 for stretch,
colour, and $M_{\textrm{stellar}}$, respectively.

Figure 6.— The $x_{1}$ distributions for Foundation and Low-z surveys. Data is
presented in black circles and simulation is presented in blue histogram.

To model the double-peaked $x_{1}$ distribution in Low-z and Foundation, we
found it necessary to implement a prior requiring that
$\overline{x_{1}}_{1}<0$ and $\overline{x_{1}}_{2}>0$ (Eq. 5). While we are
able to capture the double-peak behaviour for both surveys, shown in Figure 6,
there is a slight data-sim discrepancy near $x_{1}\sim-2$ for the Low-z
sample.

###  6.2. Impact on Cosmology Using Bias Corrections with Host Properties

With our simulations from Section 3.1 and the parent populations shown in
Figure 4, we create 100 simulated data sets and a large BiasCor to test the
accuracy of our bias correction methods. We compare four different bias
correction methods for the G10 and C11 models. Two methods use BBC1D, with and
without $M_{\textrm{stellar}}$ information included in the BiasCor (designated
‘No Mass’, similar to that of Betoule et al. (2014)). The other two methods
use the BBC5D and new BBC7D method. Overall, we find our new method
significantly improves upon BBC5D. In Table 1, we show the fitted nuisance
parameters, Hubble scatter (RMS of $\mu_{\rm fit}-\mu_{\rm true}$), $w$ and
$w$-bias defined as $\Delta w=(w_{\rm fit}-w_{\rm true})$, averaged over the
100 simulations.

For the G10 scatter model, we find an $\alpha$ bias below $1\%$ for all BBC
methods. This broad agreement holds well for the recovery of $\beta$ as well.
The BBC5D approach results in a significant $\gamma$ bias of
$\Delta\gamma=0.019\pm 0.001$ mag. BBC7D reduces this bias by a factor of
$\sim 5$, although a 2 millimag $\gamma$-bias remains with $2\sigma$
significance. These trends in $\alpha$, $\beta$, and $\gamma$ are
qualitatively similar in the C11 results. However, it is worth noting that the
BBC1D recovery of $\beta$ and $\gamma$ for the C11 scatter model is
significantly biased.

For both G10 and C11, we see a general decrease in Hubble scatter with
increasing dimensionality of the BBC methodology. BBC5D and BBC7D have
comparable Hubble scatter with each other for both G10 and C11, and both are
$\sim 10\%$ smaller compared to using BBC1D.

The presence of host-galaxy correlations in the BBC1D BiasCor does not make a
significant impact on recovered $w$ values for G10 or C11. In the case with
and without host-galaxy correlations, the $w$-bias is $\sim 0.02\pm 0.005$ for
G10 and $\sim 0.01\pm 0.005$ for C11. BBC5D does not significantly effect this
BBC1D $w$-bias, recovering $-0.0201\pm 0.0048$ and $-0.0237\pm 0.0046$ for G10
and C11 respectively. The BBC7D $w$-bias for both G10 and C11 is a factor of 2
smaller than their 5D counterparts: $w$-bias$=0.0086\pm 0.0046$ for G10 and
($0.0105\pm 0.0046$ for C11. Both BBC7D biases have $\sim 2\sigma$
significance. Overall, BBC7D has the smallest $w$-bias as well as the smallest
biases for other nuisance parameters and scatter.

In Figure 7, we show the binned distance modulus residuals as a function of
redshift for BBC1D, BBC5D, and BBC7D. BBC7D has smaller $\mu$ residuals
(mostly $\sim 0.005$ across the $z$ range) than the BBC5D approach. All BBC
$\mu$ residuals have a significant excess, $\mu_{\rm BBC}-\mu_{\rm
true}=0.01\pm 0.004$ ($2.4\sigma$), around $z=0.1$.

The accuracy of the reported distance modulus error, $\sigma_{\mu}$, in
comparison to Hubble Residual (HR) scatter is another important metric in
determining the effectiveness of BBC. Figure 8 shows BBC1D and BBC7D
$\sigma_{\mu}$ and HR scatter as a function of SNIa colour, $c$. We see better
agreement for BBC7D compared to BBC1D. Figure 8 and Table 1 show that for G10
and C11, the BBC7D methods have the smallest scatter.

Method | Model | input $\alpha=0.145$ | input $\beta=3.1$ | input $\gamma=0.05$ |  | input $w=-1$ |
---|---|---|---|---|---|---|---
|  | Fitted $\alpha$ | Fitted $\beta$ | Fitted $\gamma$ | Hubble scatter 222The error in the Hubble scatter is the R.M.S divided by $\sqrt{100}$. | Fitted $w$ | $\Delta_{w}~{}(N_{\sigma_{w}})$
BBC1D (No Mass) 333Does not include $M_{\textrm{stellar}}$ dependent parent populations in the BiasCor. | G10 | $0.1442\pm 0.0005$ | $3.1154\pm 0.0062$ | $0.0475\pm 0.0010$ | $0.149\pm 0.005$ | $-0.9828\pm 0.0047$ | $~{}0.017~{}(3.7\sigma)$
BBC1D | G10 | $0.1443\pm 0.0005$ | $3.1160\pm 0.0062$ | $0.0474\pm 0.0010$ | $0.149\pm 0.005$ | $-0.9777\pm 0.0047$ | $~{}0.022~{}(4.7\sigma)$
BBC5D | G10 | $0.1440\pm 0.0005$ | $3.0901\pm 0.0061$ | $0.0313\pm 0.0010$ | $0.138\pm 0.005$ | $-1.0201\pm 0.0046$ | $-0.0201~{}(4.4\sigma)$
BBC7D | G10 | $0.1451\pm 0.0005$ | $3.0891\pm 0.0062$ | $0.0464\pm 0.0010$ | $0.137\pm 0.005$ | $-1.0086\pm 0.0046$ | $-0.0086~{}(1.9\sigma)$
Method | Model | input $\alpha=0.145$ | input $\beta=3.8$ | input $\gamma=0.05$ |  | input $w=-1$ |
|  | Fitted $\alpha$ | Fitted $\beta$ | Fitted $\gamma$ | Hubble scatter | Fitted $w$ | $\Delta_{w}~{}(N_{\sigma_{w}})$
BBC1D (No Mass) | C11 | $0.1458\pm 0.0005$ | $3.2285\pm 0.0068$ | $0.0445\pm 0.0011$ | $0.151\pm 0.005$ | $-0.9896\pm 0.0048$ | $~{}0.0104~{}(2.6\sigma)$
BBC1D | C11 | $0.1459\pm 0.0005$ | $3.2283\pm 0.0068$ | $0.0432\pm 0.0011$ | $0.151\pm 0.005$ | $-0.9961\pm 0.0048$ | $~{}0.0039~{}(0.8\sigma)$
BBC5D | C11 | $0.1444\pm 0.0005$ | $3.7740\pm 0.0073$ | $0.0287\pm 0.0009$ | $0.132\pm 0.005$ | $-1.0237\pm 0.0046$ | $-0.0237~{}(5.2\sigma)$
BBC7D | C11 | $0.1460\pm 0.0005$ | $3.7813\pm 0.0071$ | $0.0466\pm 0.0010$ | $0.131\pm 0.005$ | $-1.0105\pm 0.0046$ | $-0.0105~{}(2.3\sigma)$
Method | Model | input $\alpha=0.145$ | input $\beta=\textrm{ND}$ 444Not Determined: For the BS20 model, input SALT2 $\beta$ $(\beta\neq\beta_{\rm SN})$ is not defined. | input $\gamma=0$ |  | input $w=-1$ |
|  | Fitted $\alpha$ | Fitted $\beta$ | Fitted $\gamma$ | Hubble scatter | Fitted $w$ | $\Delta_{w}~{}(N_{\sigma_{w}})$
BBC1D (No Mass) | BS20 | $0.1458\pm 0.0006$ | $2.7526\pm 0.0075$ | $0.0533\pm 0.0009$ | $0.121\pm 0.005$ | $-0.9545\pm 0.0045$ | $~{}0.0455~{}(10.1\sigma)$
BBC1D | BS20 | $0.1461\pm 0.0006$ | $2.7565\pm 0.0070$ | $0.0467\pm 0.0009$ | $0.116\pm 0.005$ | $-0.9572\pm 0.0045$ | $~{}0.0428~{}(9.5\sigma)$
BBC-BS20 | BS20 | $0.1417\pm 0.0004$ | $3.0770\pm 0.0050$ | $0.0006\pm 0.0007$ | $0.104\pm 0.005$ | $-0.9943\pm 0.0040$ | $~{}0.0057~{}(1.4\sigma)$
space | space | space | space | space | space | space |
|  |  |  |  |  |  |

Table 3Fitted values and uncertainties 11footnotemark: 1 averaged over 100
simulations  Figure 7.— The binned $\mu$ residuals for the G10 scatter model
using three different BBC methods on the same set of simulations.
Uncertainties are shown for BBC7D only but are representative for the other
methods. Figure 8.— The Hubble Scatter (line) and the BBC-fitted distance
modulus uncertainties, $\sigma_{\mu}$, (triangle) are plotted for three
different scatter models. We expect good agreement between the two for an
effective BBC method.

###  6.3. Impact on Cosmology Using Bias Corrections for BS20

The bottom tier of Table 1 shows how well BBC1D and BBC-BS20 perform on
simulations generated with the BS20 model. While simulations of the BS20 model
do not have an input $\gamma$ like the aforementioned simulations G10 and C11,
the simulation includes host-mass vs colour-luminosity dependence that matches
that of the data. As reported in BS20, we find that using BBC1D results in
$\gamma\sim 0.06$, consistent with the value found in the data. BBC-BS20
results in $\gamma=0.0006\pm 0.0007$, consistent with no mass step. Applying
BBC1D to BS20 simulations performs well in recovering $\alpha$
($\Delta{\alpha}{\sim}0.001\pm 0.0006$), while BBC-BS20 has a slightly larger
$\alpha$-bias ($\Delta{\alpha}\sim 0.003\pm 0.0004$). To validate the
predicted $\gamma$-reduction, we apply the BBC-BS20 method on real data and
find

$\gamma=0.0166\pm 0.0076$ (16)

which is more than $\times 3$ smaller than previous measurements.

The BS20 model is characterised by a distribution of $\beta_{\rm SN}$ and
$R_{V}$, and thus we cannot determine an input $\beta$ for the simulation. For
this reason, the BS20 section of Table 1 refers to the input $\beta$ as Not
Defined (ND), and we are unable to compute $\beta$ biases. BS20 found their
model parameters such that their BBC1D $\beta$ from simulations matched that
of the data ($\beta\sim 3.06\pm 0.06$). Since BS20 used an overly simplistic
assumption about $M_{\textrm{stellar}}$ distributions, their fitted model
parameters are approximate. Here we update the $M_{\textrm{stellar}}$
distribution to match the data (Figure 2) but do not update the BS20 model
parameters; we find BBC1D correction on BS20 simulations results in $\beta\sim
2.76\pm 0.01$, which no longer matches that of the data. However, when running
BBC-BS20 on our simulations and real data we find good agreement:
$\beta=3.077\pm 0.005$ (sim) and $\beta=3.165\pm 0.045$ (data) respectively.
This BBC1D $\beta$ discrepancy warrants the need for further improvements to
the BS20 modeling.

BBC1D results in a $w$-bias of $\sim 0.04\pm 0.004$, while the BBC-BS20
$w$-bias is significantly reduced and consistent with 0, $\sim 0.006\pm
0.004$. We find that the Hubble Scatter in the BS20 simulations is $\sim 10\%$
smaller than that of G10 or C11, which is illustrated in Figure 8.

##  7\. Discussion and Conclusion

The utilisation of host galaxy information in SNIa standardisation analyses
has become common in recent analyses. Typically, host galaxy properties have
been included in the one of the final stages of the analysis to make an
additive correction to SNIa luminosities before measuring cosmological
parameters from the set of SNIa distances. However, this approach is
conceptually flawed because it does not account for subtle biases due to the
correlations of SNIa properties and host galaxies. Here, we have determined
the underlying populations of SNIa properties and their correlations with host
galaxies so that we can trace and correct for these biases.

Our approach for these corrections is to modify the BBC method to allow for
higher dimensional bias corrections. If the intrinsic scatter model is known,
BBC determines the input mass step ($\gamma$) with a bias of 0.004 mag for all
models, and determines $w$ with biases consistent with 0 for G10 and BS20. For
$\gamma$, the BBC7D bias is a factor of 5 smaller than for BBC5D; for $w$, the
BBC7D bias is a factor of 2 smaller.

We do not address potential biases if the intrinsic scatter model is not
known. Using BBC5D, Scolnic et al. (2018) and Brout et al. (2019b) showed that
the difference in $w$ between assuming the G10 and C11 is $\sim 0.03$. Using
BBC1D, Brout & Scolnic (2020) showed that the $w$-difference when doing 1D
corrections between BS20, G10, and C11 was as much as much as 4%. However, in
that study, they did not simultaneously refit for a mass step, and it is
unclear how this approximation would affect the cosmological bias. This
intrinsic scatter uncertainty will be quantified in an upcoming analysis
(Popovic et al. in prep.), which will determine optimal parameters of the BS20
model and address the discrepancy in the fitted 1D $\beta$ values between data
and simulation.

In determining the underlying populations, we follow a purely empirical
approach that does not rely on theoretical models to relate SNIa properties to
host galaxy properties. Nicolas et al. (2020) use theoretical modeling to
predict an evolution of the relation between the $x_{1}$ distribution with
redshift due to the different sampling of progenitor systems at high-$z$
compared to low-$z$. While we are unable to show whether this evolution is due
to progenitor systematics or survey selection effects, our analysis
illustrates how the model from Nicolas et al. (2020) can be incorporated in a
cosmological analysis.

From our empirical approach, we find that peak probabilities of the parent
distributions of color and stretch do not evolve significantly with host-
galaxy stellar mass, but rather that the observed dependence on
$M_{\textrm{stellar}}$ is explained by increasing asymmetries of the parent
distributions. For $c$, this finding agrees with the model proposed in BS20,
who suggested that dust is responsible for one-sided SNIa $c$ scatter towards
redder colors. Similarly, Rigault et al. (2018) propose that local specific
star formation rate (lsSFR) is a tracer of the progenitor age which itself
traces the SNIa $x_{1}$ distribution. They find that the parent distributions
of $x_{1}$ for low lsSFR and high lsSFR are inconsistent. We find that the
observed correlation between $x_{1}$ and mass is driven by an increasing faint
side standard deviation ($\sigma_{-}$) of the parent distribution, which
suggests a more subtle relationship between $x_{1}$ and mass that can be
described by a non-evolving peak probability and an evolving asymmetry.

Finally, while we focused here on using host-galaxy mass to correlate with
SNIa properties, our methods can be applied to other host galaxy properties.
Furthermore, while we used a step function for the mass dependence, an
arbitrary functional form can be used with BBC7D. Ultimately, our improved BBC
method will be evaluated by varying the SNIa model such that simulated
distributions match those of the data, and propagating these variations to the
BiasCor used by BBC.

###  7.1. Acknowledgements

Dillon Brout acknowledges support for this work was provided by NASA through
the NASA Hubble Fellowship grant HST-HF2-51430.001 awarded by the Space
Telescope Science Institute, which is operated by Association of Universities
for Research in Astronomy, Inc., for NASA, under contract NAS5-26555.

Richard Kessler acknowledges that this work was supported in part by the Kavli
Institute for Cosmological Physics at the University of Chicago through an
endowment from the Kavli Foundation and its founder Fred Kavli. R.K. is
supported by DOE grant DE-SC0009924. Dan Scolnic is supported by DOE grant DE-
SC0010007 and the David and Lucile Packard Foundation.

###  .2. A1. EFFICIENCIES AND ERRATA

For the SNLS sample, we improve the spectroscopic efficiency function
described as a function of peak $i$-band mag, Eff${}_{\textrm{spec}(i)}$.
Figure 9 shows that the JLA Eff${}_{\textrm{spec}}$, also used in BS20,
underestimates the number of supernovae at redshifts greater than $0.75$. Our
update is to replace Eff${}_{\textrm{spec}(i)}$ with
Eff${}_{\textrm{spec}(i-0.3)}$.

Figure 9.— The SNLS redshift distribution for the data (solid black line with
points), simulations with JLA spectroscopic efficiency (orange dash-dot line),
and updated spectroscopic efficiency presented here (blue solid line). At
higher redshifts, the new efficiency better matches the data.

###  .3. A2. PARENT POPULATION PARAMETERS

Here we show the parent populations as a function of mass for each survey,
characterised by an asymmetric generalised normal distribution with three
parameters: a peak probability and two standard deviations. See Equation 4 for
rigorous definition. For $c$, the $n$ value is fixed at $n=2$. For $x_{1}$,
$n=3$. For the Low-z and Foundation surveys, we present a double-Gaussian
approach characterised in Equation 5. Population fit results for $c$ are in
Tables 4, 6, 9 for Low-z, Foundation, and the combined DES, SDSS, SNLS, and
PS1 sample, respectively. The population fit results for $x_{1}$ are in Tables
5 7, 8 for Low-z, Foundation, and the combined DES, SDSS, SNLS, and PS1
sample, respectively. These populations are shown in Figures 10 and 11 for $c$
and $x_{1}$ for the G10 scatter model.

Figure 10.— The G10 $c$ parent populations as a function of host-galaxy mass for several spectroscopic surveys. The mean of the asymmetric distribution is presented in blue, the peak in orange, the bright-side standard deviation in green, and the faint-side standard deviation in red. Errors are included for each parameter and the 68% confidence interval for each is shown in grey fill.  Figure 11.— The G10 $x_{1}$ parent populations as a function of host-galaxy mass for several spectroscopic surveys. The mean of the asymmetric distribution is presented in blue, the peak in orange, the faint-side standard deviation in green, and the bright-side standard deviation in red. Errors are included for each parameter and the 68% confidence interval for each is shown in grey fill. For Foundation and Low-z, two peak values are shown by purple lines, and the line thickness denotes the relative weight of each peak. Mass | mean (G10) |  $\sigma_{-}$ (G10) |  $\sigma_{+}$ (G10) | mean (C11) |  $\sigma_{-}$ (C11) |  $\sigma_{+}$ (C11)
---|---|---|---|---|---|---
9.8 | $-0.043\pm 0.042$ | $0.046\pm 0.037$ | $0.164\pm 0.039$ | $-0.042\pm 0.044$ | $0.043\pm 0.038$ | $0.163\pm 0.044$
10.0 | $-0.072\pm 0.02$ | $0.018\pm 0.016$ | $0.151\pm 0.035$ | $-0.069\pm 0.023$ | $0.02\pm 0.019$ | $0.157\pm 0.041$
10.2 | $-0.075\pm 0.017$ | $0.014\pm 0.012$ | $0.134\pm 0.022$ | $-0.067\pm 0.018$ | $0.015\pm 0.013$ | $0.128\pm 0.025$
10.4 | $-0.074\pm 0.015$ | $0.013\pm 0.01$ | $0.154\pm 0.018$ | $-0.065\pm 0.016$ | $0.013\pm 0.011$ | $0.147\pm 0.018$
10.6 | $-0.073\pm 0.014$ | $0.012\pm 0.009$ | $0.136\pm 0.015$ | $-0.062\pm 0.015$ | $0.012\pm 0.01$ | $0.126\pm 0.016$
10.8 | $-0.072\pm 0.02$ | $0.017\pm 0.013$ | $0.132\pm 0.017$ | $-0.061\pm 0.018$ | $0.016\pm 0.012$ | $0.118\pm 0.017$
11.0 | $-0.067\pm 0.023$ | $0.021\pm 0.015$ | $0.138\pm 0.022$ | $-0.061\pm 0.021$ | $0.017\pm 0.014$ | $0.126\pm 0.022$
11.2 | $-0.044\pm 0.034$ | $0.036\pm 0.022$ | $0.12\pm 0.029$ | $-0.044\pm 0.031$ | $0.028\pm 0.021$ | $0.111\pm 0.028$
11.4 | $-0.019\pm 0.039$ | $0.043\pm 0.027$ | $0.094\pm 0.037$ | $-0.024\pm 0.039$ | $0.033\pm 0.025$ | $0.09\pm 0.037$
11.6 | $0.029\pm 0.035$ | $0.064\pm 0.03$ | $0.036\pm 0.029$ | $0.019\pm 0.033$ | $0.047\pm 0.03$ | $0.035\pm 0.031$
Table 4LOWZ parent colour populations Mass | Weight (G10) |  $\overline{x}_{1}$ (G10) |  $\sigma$ (G10) | Weight (C11) |  $\overline{x}_{1}$ (C11) |  $\sigma$ (C11)
---|---|---|---|---|---|---
9.8 | $1.874\pm 1.023$ | $-1.485\pm 0.644$ | $1.438\pm 0.454$ | $3.239\pm 1.204$ | $0.402\pm 0.181$ | $0.595\pm 0.153$
10.0 | $2.402\pm 1.125$ | $-1.917\pm 0.342$ | $0.835\pm 0.52$ | $3.254\pm 1.169$ | $0.273\pm 0.142$ | $0.693\pm 0.123$
10.2 | $3.271\pm 1.181$ | $-1.884\pm 0.102$ | $0.502\pm 0.107$ | $2.265\pm 1.144$ | $0.418\pm 0.635$ | $0.805\pm 0.282$
10.4 | $3.371\pm 1.172$ | $-1.836\pm 0.163$ | $0.458\pm 0.169$ | $1.46\pm 0.579$ | $0.24\pm 0.116$ | $0.816\pm 0.118$
10.6 | $3.374\pm 1.179$ | $-1.816\pm 0.058$ | $0.432\pm 0.043$ | $1.367\pm 0.665$ | $0.495\pm 0.603$ | $0.866\pm 0.333$
10.8 | $3.365\pm 1.186$ | $-1.78\pm 0.055$ | $0.433\pm 0.041$ | $1.366\pm 0.539$ | $0.394\pm 0.108$ | $0.73\pm 0.097$
11.0 | $3.333\pm 1.161$ | $-1.788\pm 0.059$ | $0.452\pm 0.046$ | $1.471\pm 0.592$ | $0.418\pm 0.101$ | $0.691\pm 0.111$
11.2 | $3.362\pm 1.165$ | $-1.678\pm 0.064$ | $0.451\pm 0.06$ | $1.04\pm 0.544$ | $0.536\pm 0.127$ | $0.578\pm 0.12$
11.4 | $3.326\pm 1.181$ | $-1.613\pm 0.085$ | $0.408\pm 0.111$ | $0.917\pm 0.481$ | $0.577\pm 0.163$ | $0.729\pm 0.202$
11.6 | $3.42\pm 1.143$ | $-1.719\pm 0.377$ | $1.341\pm 0.624$ | $0.696\pm 0.61$ | $1.023\pm 0.827$ | $1.662\pm 0.688$
Table 5LOWZ parent stretch populations Mass | mean (G10) |  $\sigma_{-}$ (G10) |  $\sigma_{+}$ (G10) | mean (C11) |  $\sigma_{-}$ (C11) |  $\sigma_{+}$ (C11)
---|---|---|---|---|---|---
10.0 | $0.01\pm 0.049$ | $0.071\pm 0.03$ | $0.218\pm 0.066$ | $-0.046\pm 0.05$ | $0.048\pm 0.034$ | $0.212\pm 0.06$
10.2 | $0.004\pm 0.041$ | $0.06\pm 0.031$ | $0.192\pm 0.082$ | $-0.069\pm 0.041$ | $0.037\pm 0.029$ | $0.243\pm 0.041$
10.4 | $-0.004\pm 0.04$ | $0.064\pm 0.025$ | $0.2\pm 0.055$ | $-0.07\pm 0.034$ | $0.033\pm 0.023$ | $0.211\pm 0.054$
10.6 | $0.01\pm 0.047$ | $0.067\pm 0.027$ | $0.181\pm 0.057$ | $-0.055\pm 0.039$ | $0.042\pm 0.026$ | $0.206\pm 0.057$
10.8 | $0.011\pm 0.051$ | $0.069\pm 0.028$ | $0.168\pm 0.062$ | $-0.028\pm 0.048$ | $0.061\pm 0.031$ | $0.198\pm 0.067$
11.0 | $0.024\pm 0.05$ | $0.072\pm 0.029$ | $0.085\pm 0.043$ | $0.026\pm 0.049$ | $0.086\pm 0.031$ | $0.115\pm 0.08$
Table 6Foundation parent colour populations Mass | Weight (G10) |  $\overline{x}_{1}$ (G10) |  $\sigma$ (G10) | Weight (C11) |  $\overline{x}_{1}$ (C11) |  $\sigma$ (C11)
---|---|---|---|---|---|---
10.0 | $2.943\pm 1.3$ | $-0.972\pm 0.506$ | $0.97\pm 0.554$ | $2.217\pm 1.367$ | $0.645\pm 0.543$ | $0.955\pm 0.487$
10.2 | $3.36\pm 1.158$ | $-1.167\pm 0.153$ | $0.493\pm 0.192$ | $1.819\pm 0.99$ | $0.682\pm 0.28$ | $0.654\pm 0.334$
10.4 | $3.317\pm 1.179$ | $-1.152\pm 0.097$ | $0.471\pm 0.101$ | $1.285\pm 0.743$ | $0.731\pm 0.278$ | $0.631\pm 0.319$
10.6 | $3.342\pm 1.162$ | $-1.173\pm 0.092$ | $0.559\pm 0.103$ | $1.694\pm 1.011$ | $0.942\pm 0.189$ | $0.423\pm 0.308$
10.8 | $3.32\pm 1.179$ | $-1.181\pm 0.097$ | $0.522\pm 0.1$ | $1.048\pm 0.732$ | $0.896\pm 0.276$ | $0.574\pm 0.399$
11.0 | $3.257\pm 1.215$ | $-1.264\pm 0.152$ | $0.618\pm 0.256$ | $0.904\pm 0.644$ | $0.845\pm 0.451$ | $0.914\pm 0.581$
Table 7Foundation parent stretch populations Mass | mean (G10) |  $\sigma_{-}$ (G10) |  $\sigma_{+}$ (G10) | mean (C11) |  $\sigma_{-}$ (C11) |  $\sigma_{+}$ (C11)
---|---|---|---|---|---|---
8.2 | $-0.079\pm 0.026$ | $0.021\pm 0.016$ | $0.098\pm 0.025$ | $-0.07\pm 0.029$ | $0.023\pm 0.019$ | $0.085\pm 0.027$
8.4 | $-0.082\pm 0.02$ | $0.016\pm 0.013$ | $0.098\pm 0.018$ | $-0.048\pm 0.03$ | $0.024\pm 0.018$ | $0.062\pm 0.023$
8.6 | $-0.085\pm 0.015$ | $0.014\pm 0.01$ | $0.088\pm 0.014$ | $-0.057\pm 0.021$ | $0.016\pm 0.012$ | $0.047\pm 0.02$
8.8 | $-0.077\pm 0.023$ | $0.027\pm 0.015$ | $0.1\pm 0.017$ | $-0.072\pm 0.022$ | $0.018\pm 0.014$ | $0.068\pm 0.018$
9.0 | $-0.084\pm 0.019$ | $0.022\pm 0.013$ | $0.108\pm 0.015$ | $-0.077\pm 0.018$ | $0.014\pm 0.011$ | $0.076\pm 0.015$
9.2 | $-0.081\pm 0.019$ | $0.024\pm 0.013$ | $0.11\pm 0.015$ | $-0.054\pm 0.025$ | $0.024\pm 0.015$ | $0.049\pm 0.018$
9.4 | $-0.074\pm 0.02$ | $0.028\pm 0.013$ | $0.111\pm 0.015$ | $-0.065\pm 0.021$ | $0.02\pm 0.013$ | $0.063\pm 0.015$
9.6 | $-0.063\pm 0.019$ | $0.034\pm 0.013$ | $0.114\pm 0.016$ | $-0.06\pm 0.02$ | $0.02\pm 0.013$ | $0.06\pm 0.015$
9.8 | $-0.062\pm 0.02$ | $0.034\pm 0.014$ | $0.117\pm 0.016$ | $-0.034\pm 0.026$ | $0.033\pm 0.017$ | $0.046\pm 0.018$
10.0 | $-0.086\pm 0.012$ | $0.012\pm 0.009$ | $0.138\pm 0.012$ | $-0.063\pm 0.017$ | $0.017\pm 0.012$ | $0.075\pm 0.012$
10.2 | $-0.08\pm 0.015$ | $0.018\pm 0.011$ | $0.138\pm 0.013$ | $-0.073\pm 0.014$ | $0.014\pm 0.009$ | $0.096\pm 0.011$
10.4 | $-0.079\pm 0.016$ | $0.019\pm 0.011$ | $0.139\pm 0.013$ | $-0.079\pm 0.015$ | $0.014\pm 0.01$ | $0.104\pm 0.011$
10.6 | $-0.063\pm 0.018$ | $0.03\pm 0.013$ | $0.125\pm 0.014$ | $-0.073\pm 0.016$ | $0.017\pm 0.011$ | $0.098\pm 0.012$
10.8 | $-0.064\pm 0.018$ | $0.033\pm 0.012$ | $0.125\pm 0.014$ | $-0.079\pm 0.017$ | $0.018\pm 0.013$ | $0.106\pm 0.011$
11.0 | $-0.068\pm 0.018$ | $0.031\pm 0.013$ | $0.128\pm 0.015$ | $-0.085\pm 0.014$ | $0.014\pm 0.011$ | $0.11\pm 0.01$
11.2 | $-0.062\pm 0.019$ | $0.036\pm 0.013$ | $0.123\pm 0.015$ | $-0.085\pm 0.015$ | $0.015\pm 0.011$ | $0.108\pm 0.011$
11.4 | $-0.065\pm 0.022$ | $0.039\pm 0.016$ | $0.125\pm 0.019$ | $-0.086\pm 0.019$ | $0.018\pm 0.013$ | $0.109\pm 0.016$
11.6 | $-0.074\pm 0.025$ | $0.034\pm 0.017$ | $0.118\pm 0.023$ | $-0.084\pm 0.021$ | $0.02\pm 0.014$ | $0.097\pm 0.018$
11.8 | $-0.098\pm 0.023$ | $0.022\pm 0.016$ | $0.138\pm 0.029$ | $-0.097\pm 0.026$ | $0.022\pm 0.017$ | $0.111\pm 0.026$
12.0 | $-0.073\pm 0.031$ | $0.03\pm 0.023$ | $0.118\pm 0.049$ | $-0.06\pm 0.04$ | $0.036\pm 0.028$ | $0.091\pm 0.047$
Table 8Combined SDSS, DES, SNLS, PS1 parent colour populations Mass | mean (G10) |  $\sigma_{-}$ (G10) |  $\sigma_{+}$ (G10) | mean (C11) |  $\sigma_{-}$ (C11) |  $\sigma_{+}$ (C11)
---|---|---|---|---|---|---
8.2 | $0.908\pm 0.288$ | $0.8\pm 0.241$ | $0.249\pm 0.201$ | $0.559\pm 0.227$ | $0.726\pm 0.751$ | $0.729\pm 0.687$
8.4 | $0.39\pm 0.349$ | $0.377\pm 0.245$ | $0.631\pm 0.244$ | $0.473\pm 0.312$ | $0.584\pm 0.611$ | $0.725\pm 0.573$
8.6 | $0.511\pm 0.416$ | $0.633\pm 0.305$ | $0.575\pm 0.287$ | $0.832\pm 0.281$ | $0.835\pm 0.212$ | $0.307\pm 0.19$
8.8 | $0.537\pm 0.365$ | $0.571\pm 0.265$ | $0.464\pm 0.244$ | $0.791\pm 0.253$ | $0.796\pm 0.182$ | $0.288\pm 0.175$
9.0 | $0.232\pm 0.277$ | $0.374\pm 0.202$ | $0.688\pm 0.188$ | $0.458\pm 0.273$ | $0.575\pm 0.196$ | $0.532\pm 0.19$
9.2 | $0.242\pm 0.253$ | $0.421\pm 0.176$ | $0.685\pm 0.171$ | $0.375\pm 0.287$ | $0.543\pm 0.201$ | $0.625\pm 0.198$
9.4 | $0.173\pm 0.273$ | $0.43\pm 0.206$ | $0.761\pm 0.172$ | $0.59\pm 0.315$ | $0.766\pm 0.229$ | $0.517\pm 0.214$
9.6 | $0.43\pm 0.31$ | $0.681\pm 0.228$ | $0.591\pm 0.21$ | $0.52\pm 0.295$ | $0.756\pm 0.21$ | $0.54\pm 0.201$
9.8 | $0.761\pm 0.257$ | $1.023\pm 0.177$ | $0.336\pm 0.192$ | $0.608\pm 0.246$ | $0.892\pm 0.171$ | $0.48\pm 0.178$
10.0 | $0.871\pm 0.225$ | $1.198\pm 0.155$ | $0.274\pm 0.169$ | $0.763\pm 0.239$ | $1.088\pm 0.166$ | $0.352\pm 0.174$
10.2 | $0.938\pm 0.202$ | $1.422\pm 0.142$ | $0.218\pm 0.15$ | $0.757\pm 0.262$ | $1.272\pm 0.179$ | $0.352\pm 0.193$
10.4 | $0.793\pm 0.26$ | $1.553\pm 0.187$ | $0.323\pm 0.189$ | $0.441\pm 0.307$ | $1.138\pm 0.212$ | $0.594\pm 0.22$
10.6 | $0.868\pm 0.221$ | $1.741\pm 0.162$ | $0.222\pm 0.149$ | $0.67\pm 0.309$ | $1.429\pm 0.227$ | $0.389\pm 0.218$
10.8 | $0.697\pm 0.26$ | $1.763\pm 0.186$ | $0.33\pm 0.176$ | $0.597\pm 0.338$ | $1.517\pm 0.246$ | $0.433\pm 0.237$
11.0 | $0.655\pm 0.285$ | $1.836\pm 0.207$ | $0.353\pm 0.191$ | $0.585\pm 0.35$ | $1.579\pm 0.258$ | $0.431\pm 0.243$
11.2 | $0.27\pm 0.334$ | $1.614\pm 0.25$ | $0.612\pm 0.227$ | $0.513\pm 0.344$ | $1.616\pm 0.257$ | $0.463\pm 0.241$
11.4 | $0.062\pm 0.41$ | $1.571\pm 0.301$ | $0.745\pm 0.281$ | $0.73\pm 0.243$ | $1.79\pm 0.192$ | $0.253\pm 0.17$
11.6 | $0.079\pm 0.512$ | $1.433\pm 0.398$ | $0.681\pm 0.325$ | $0.637\pm 0.31$ | $1.671\pm 0.258$ | $0.323\pm 0.204$
11.8 | $0.309\pm 0.512$ | $1.852\pm 0.428$ | $0.535\pm 0.363$ | $0.608\pm 0.392$ | $1.828\pm 0.547$ | $0.352\pm 0.274$
12.0 | $-0.14\pm 0.845$ | $1.947\pm 0.69$ | $1.137\pm 0.672$ | $0.379\pm 0.605$ | $2.12\pm 0.538$ | $0.793\pm 0.657$
Table 9Combined SDSS, DES, SNLS, PS1 parent stretch populations

## References

  * Abbott et al. (2019) Abbott, T. M. C., Allam, S., Andersen, P., et al., 2019, ApJ, 872, 2, L30, arXiv:1811.02374
  * Baker & Cousins (1984) Baker, S., Cousins, R., 1984, Nucl. Instrum. Methods Phys. Res., 221, A437
  * Betoule et al. (2014) Betoule, M., Kessler, R., Guy, J., et al., 2014, A&A, 568, A22, arXiv:1401.4064
  * Brout et al. (2020) Brout, D., Hinton, S., Scolnic, D., 2020, arXiv:2012.05900
  * Brout et al. (2019a) Brout, D., Sako, M., Scolnic, D., et al., 2019a, The Astrophysical Journal, 874, 1, 106, ISSN 1538-4357
  * Brout & Scolnic (2020) Brout, D., Scolnic, D., 2020, arXiv:2004.10206
  * Brout et al. (2019b) Brout, D., Scolnic, D., Kessler, R., et al., 2019b, The Astrophysical Journal, 874, 2, 150, ISSN 1538-4357
  * Carrick et al. (2015) Carrick, J., Turnbull, S. J., Lavaux, G., Hudson, M. J., 2015, MNRAS, 450, 317, arXiv:1504.04627
  * Chotard et al. (2011) Chotard, N., Gangler, E., Aldering, G., et al., 2011, A&A, 529, L4, arXiv:1103.5300
  * Foley et al. (2018) Foley, R. J., Scolnic, D., Rest, A., et al., 2018, MNRAS, 475, 193, arXiv:1711.02474
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., Goodman, J., 2013, PASP, 125, 925, 306, arXiv:1202.3665
  * Guy et al. (2010) Guy, J., Sullivan, M., Conley, A., et al., 2010, A&A, 523, A7, arXiv:1010.4743
  * Hicken et al. (2009a) Hicken, M., Challis, P., Jha, S., et al., 2009a, ApJ, 700, 331, arXiv:0901.4787
  * Hicken et al. (2012) Hicken, M., Challis, P., Kirshner, R. P., et al., 2012, ApJS, 200, 12, arXiv:1205.4493
  * Hicken et al. (2009b) Hicken, M., Wood-Vasey, W. M., Blondin, S., et al., 2009b, ApJ, 700, 1097, arXiv:0901.4804
  * James & Roos (1975) James, F., Roos, M., 1975, Computer Physics Communications, 10, 6, 343
  * Jones et al. (2018) Jones, D. O., Riess, A. G., Scolnic, D. M., et al., 2018, ApJ, 867, 108, arXiv:1805.05911
  * Jones et al. (2019) Jones, D. O., Scolnic, D. M., Foley, R. J., et al., 2019, The Astrophysical Journal, 881, 1, 19, ISSN 1538-4357
  * Kessler et al. (2009a) Kessler, R., Becker, A. C., Cinabro, D., et al., 2009a, ApJS, 185, 32, arXiv:0908.4274
  * Kessler et al. (2009b) Kessler, R., Bernstein, J. P., Cinabro, D., et al., 2009b, Publications of the Astronomical Society of the Pacific, 121, 1028, arXiv:0908.4280
  * Kessler et al. (2013) Kessler, R., Guy, J., Marriner, J., et al., 2013, ApJ, 764, 48, arXiv:1209.2482
  * Kessler et al. (2019) Kessler, R., Narayan, G., Avelino, A., et al., 2019, Publications of the Astronomical Society of the Pacific, 131, 1003, 094501, arXiv:1903.11756
  * Kessler & Scolnic (2017) Kessler, R., Scolnic, D., 2017, ApJ, 836, 56, arXiv:1610.04677
  * Marriner et al. (2011) Marriner, J., Bernstein, J. P., Kessler, R., et al., 2011, ApJ, 740, 72, arXiv:1107.4631
  * Nicolas et al. (2020) Nicolas, N., Rigault, M., Copin, Y., et al., 2020, arXiv e-prints, arXiv:2005.09441, arXiv:2005.09441
  * Perlmutter et al. (1999) Perlmutter, S., Aldering, G., Goldhaber, G., et al., 1999, ApJ, 517, 565, astro-ph/9812133
  * Planck Collaboration et al. (2016) Planck Collaboration, Ade, P. A. R., Aghanim, N., et al., 2016, A&A, 594, A14, arXiv:1502.01590
  * Popovic et al. (2019) Popovic, B., Scolnic, D., Kessler, R., 2019, arXiv e-prints, arXiv:1910.05228, arXiv:1910.05228
  * Riess et al. (1998) Riess, A. G., Filippenko, A. V., Challis, P., et al., 1998, AJ, 116, 1009, astro-ph/9805201
  * Rigault et al. (2018) Rigault, M., Brinnel, V., Aldering, G., et al., 2018, arXiv e-prints, arXiv:1806.03849
  * Rose et al. (2020) Rose, B. M., Rubin, D., Cikota, A., et al., 2020, ApJ, 896, 1, L4, arXiv:2002.12382
  * Sako et al. (2018) Sako, M., Bassett, B., Becker, A. C., et al., 2018, Publications of the Astronomical Society of the Pacific, 130, 064002, arXiv:1401.3317
  * Scolnic et al. (2015) Scolnic, D., Casertano, S., Riess, A., et al., 2015, ApJ, 815, 2, 117, arXiv:1508.05361
  * Scolnic & Kessler (2016) Scolnic, D., Kessler, R., 2016, ApJ, 822, L35, arXiv:1603.01559
  * Scolnic et al. (2018) Scolnic, D., Kessler, R., Brout, D., et al., 2018, ApJ, 852, L3, arXiv:1710.05845
  * Smith et al. (2020) Smith, M., Sullivan, M., Wiseman, P., et al., 2020, arXiv e-prints, arXiv:2001.11294, arXiv:2001.11294
  * Stritzinger et al. (2011) Stritzinger, M. D., Phillips, M. M., Boldt, L. N., et al., 2011, AJ, 142, 156, arXiv:1108.3108
  * Sullivan et al. (2010) Sullivan, M., Conley, A., Howell, D. A., et al., 2010, MNRAS, 406, 782, arXiv:1003.5119
  * Tripp (1998) Tripp, R., 1998, A&A, 331, 815
  * Uddin et al. (2017) Uddin, S. A., Mould, J., Lidman, C., Ruhlmann-Kleider, V., Zhang, B. R., 2017, ApJ, 848, 1, 56, arXiv:1709.05830
  * Wiseman et al. (2020) Wiseman, P., Smith, M., Childress, M., et al., 2020, MNRAS, 495, 4, 4040, arXiv:2001.02640
  * York et al. (2000) York, D. G., Adelman, J., Anderson, J., John E., et al., 2000, AJ, 120, 3, 1579, arXiv:astro-ph/0006396