# A galaxy-driven model of type Ia supernova luminosity variations

P. Wiseman, M. Vincenzi, M. Sullivan, L. Kelsey, B. Popovic, B. Rose, D.

Brout, T.M. Davis, C. Frohmaier, L. Galbany, et al.

► To cite this version:

P. Wiseman, M. Vincenzi, M. Sullivan, L. Kelsey, B. Popovic, et al. A galaxy-driven model of type Ia supernova luminosity variations. Monthly Notices of the Royal Astronomical Society, 2022, 515 (3), pp.4587-4605. 10.1093/mnras/stac1984. hal-03728180

HAL Id: hal-03728180

https://hal.science/hal-03728180v1

Submitted on 8 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A galaxy-driven model of type Ia supernova luminosity variations

P. Wiseman<sup>1,★</sup> M. Vincenzi<sup>1,2</sup> M. Sullivan<sup>1</sup> L. Kelsey<sup>1,3</sup> B. Popovic,<sup>2</sup> B. Rose,<sup>2</sup> D. Brout,<sup>4,†</sup>

T. M. Davis<sup>5</sup> C. Frohmaier<sup>1,3</sup> L. Galbany<sup>6,7</sup> C. Lidman,<sup>8,9</sup> A. Möller,<sup>10</sup> D. Scolnic<sup>2</sup>

M. Smith<sup>11</sup> M. Aguena,<sup>12</sup> S. Allam,<sup>13</sup> F. Andrade-Oliveira,<sup>14</sup> J. Annis,<sup>13</sup> E. Bertin,<sup>15,16</sup> S. Bocquet,<sup>17</sup>

D. Brooks,<sup>18</sup> D. L. Burke,<sup>19,20</sup> A. Carnero Rosell,<sup>12,21,22</sup> M. Carrasco Kind,<sup>23,24</sup> J. Carretero,<sup>25</sup>

F. J. Castander,<sup>6,7</sup> M. Costanzi,<sup>26,27,28</sup> M. E. S. Pereira,<sup>29</sup> S. Desai,<sup>30</sup> H. T. Diehl,<sup>13</sup> P. Doel,<sup>18</sup> S. Everett,<sup>31</sup>

I. Ferrero,<sup>32</sup> D. Friedel,<sup>23</sup> J. Frieman,<sup>13,33</sup> J. García-Bellido,<sup>34</sup> M. Gatti,<sup>35</sup> E. Gaztanaga,<sup>6,7</sup> D. Gruen,<sup>17</sup>

J. Gschwend,<sup>12,36</sup> G. Gutierrez,<sup>13</sup> S. R. Hinton,<sup>5</sup> D. L. Hollowood,<sup>37</sup> K. Honscheid,<sup>38,39</sup> D. J. James,<sup>4</sup>

M. March,<sup>35</sup> F. Menanteau,<sup>23,24</sup> R. Miquel,<sup>25,40</sup> R. Morgan,<sup>41</sup> A. Palmese,<sup>42</sup> F. Paz-Chinchón,<sup>23,43</sup>

A. Pieres,<sup>12,36</sup> A. A. Plazas Malagón,<sup>44</sup> A. K. Romer,<sup>45</sup> E. Sanchez,<sup>46</sup> V. Scarpine,<sup>13</sup> I. Sevilla-Noarbe,<sup>46</sup>

M. Soares-Santos,<sup>14</sup> E. Suchyta,<sup>47</sup> G. Tarle,<sup>14</sup> C. To,<sup>38</sup> and T. N. Varga<sup>48,49,50</sup> (DES Collaboration)

*Affiliations are listed at the end of the paper*

Accepted 2022 July 11. Received 2022 July 6; in original form 2022 May 3

## ABSTRACT

Type Ia supernovae (SNe Ia) are used as standardizable candles to measure cosmological distances, but differences remain in their corrected luminosities which display a magnitude step as a function of host galaxy properties such as stellar mass and rest-frame *U*-*R* colour. Identifying the cause of these steps is key to cosmological analyses and provides insight into SN physics. Here we investigate the effects of SN progenitor ages on their light-curve properties using a galaxy-based forward model that we compare to the Dark Energy Survey 5-yr SN Ia sample. We trace SN Ia progenitors through time and draw their light-curve width parameters from a bimodal distribution according to their age. We find that an intrinsic luminosity difference between SNe of different ages cannot explain the observed trend between step size and SN colour. The data split by stellar mass are better reproduced by following recent work implementing a step in total-to-selective dust extinction ratio (*R*<sub>V</sub>) between low- and high-mass hosts, although an additional intrinsic luminosity step is still required to explain the data split by host galaxy *U*-*R*. Modelling the *R*<sub>V</sub> step as a function of galaxy age provides a better match overall. Additional age versus luminosity steps marginally improve the match to the data, although most of the step is absorbed by the width versus luminosity coefficient α. Furthermore, we find no evidence that α varies with SN age.

**Key words:** supernovae: general – dust, extinction – galaxies: evolution – cosmology: observations.

## 1 INTRODUCTION

Empirical relationships among light-curve properties of type Ia supernovae (SNe Ia) are routinely used to standardize their peak brightness (Pskovskii 1977; Phillips 1993; Riess, Press & Kirshner 1996; Tripp 1998; Mandel et al. 2017) and facilitate their use as distance indicators in cosmological measurements (Riess et al. 1998; Perlmutter et al. 1999; Sullivan et al. 2011; Betoule et al. 2014; DES Collaboration 2018; Scolnic et al. 2018; Brout et al. 2022). This standardization typically exploits the ‘faster–fainter’ relation between peak brightness and SN light-curve width, and the ‘bluer–brighter’ relation between peak brightness and SN optical colour. A detailed astrophysical understanding of these relationships is not yet in place, but they are assumed to derive from physical

processes intrinsic to the SN progenitors and their explosions (e.g. SN explosion mechanism and physics, or progenitor white dwarf mass, age or metallicity), or in the circumstellar medium surrounding the explosion, or in the interstellar medium between the SN and the observer (e.g. dust). Established correlations between the observed SN light-curve width and the properties of the SN host galaxy, with faster-declining SNe more readily located in more massive, less actively star-forming galaxies, presumably reflect some of this underlying astrophysics.

After standardizing SN luminosities for SN light-curve width and SN colour, any remaining difference between inferred distances and those predicted by the cosmological model are termed ‘Hubble residuals’. These Hubble residuals are also observed to correlate with host galaxy properties, both those local to the SN and the integrated global properties. The most commonly used of these correlations is the between Hubble residual and host galaxy stellar mass (Kelly et al. 2010; Lampeitl et al. 2010; Sullivan et al. 2010; Childress et al. 2013), which is typically represented by a step function in

<footnote>* E-mail: P.S.Wiseman@soton.ac.uk</footnote>

<footnote>† Einstein Fellow</footnote>
SN peak brightness of order 0.1 mag at a host galaxy stellar mass of around $log(M_*/M_\odot) = 10$. A similar step in peak brightness is also observed when considering other host galaxy properties, such as specific star-formation rate (sSFR; Sullivan et al. 2010; Rigault et al. 2020), gas phase metallicity (Gallagher et al. 2005; Childress et al. 2013), rest-frame galaxy colour (Roman et al. 2018; Kelsey et al. 2021), and emission line equivalent widths (Dixon et al. 2022). By assuming that the step is caused by two populations of SNe Ia that can be (imperfectly) probed by one of these environmental tracers, Briday et al. (2022) showed that the size of the step is correlated with the level of population contamination across the step location of the given host galaxy tracer, with spectroscopically observed local sSFR and global stellar mass displaying the least contamination and largest steps. One interpretation of local sSFR being the best step tracer is that SNe Ia of different ages have different standardized luminosities, perhaps representing two (or more) SN Ia populations. Indeed, the case for multiple populations is strengthened by Nicolas et al. (2021) who show that the distribution of light-curve widths can be modelled as a combination of a young, slow-declining population and an old, fast-declining population.

Recent analyses have revealed a dichotomy in sizes of host galaxy–Hubble residual steps when the SN samples are divided based on the colour of the SNe. In the 3-yr spectroscopically-confirmed sample of SNe Ia from the Dark Energy Survey (DES3YR; Brout et al. 2019a), redder SNe have large steps as a function of host stellar mass, while bluer SNe have negligible steps (Kelsey et al. 2021). This effect is also present in the larger, photometrically-classified DES5YR sample Kelsey et al. (2022, hereafter K22), who also showed that the effect is present when considering host galaxy rest-frame $U-R$ colour in place of stellar mass.

To explain the complex relationships between host galaxy Hubble residual steps and SN colour Brout et al. 2021 (hereafter BS21) introduced a framework where the slope of the dust extinction law along the SN line-of-sight, as governed by the total-to-selective extinction parameter $R_V$, is different in low mass (mean $R_V = 2.75$) and high mass (mean $R_V = 1.5$) SN host galaxies. The effect of this model is that dustier SNe are redder and fainter, but the level of dimming is different for a SN with the same colour depending on the stellar mass of its host galaxy. Recent analysis by Popovic et al. (2021a) of the larger Pantheon + SN Ia compilation (Scolnic et al. 2021) confirms the results from BS21, and shows that correcting SN Ia luminosities using this model may improve their use as cosmological probes. However, although the BS21 model reproduces the divergence of the stellar mass step as a function of SN colour, K22 show that after correcting for the trend of Hubble residual versus SN colour when split by stellar mass, there remains a small residual step in host galaxy $U-R$, with SNe in redder galaxies appearing brighter after correction. As this $U-R$ colour is related to stellar population age, this suggests there may also be an age-based effect in addition to dust.

In this work, we build upon the model of BS21 to qualitatively investigate the relationships between host galaxy mass, $R_V$, galaxy colour, and age, with the aim of explaining the residual $U-R$ step observed by K22 by expanding upon the multiple SN Ia age populations modelled by Nicolas et al. (2021).

Most modern simulations of SN Ia data sets use the SuperNova Analysis (SNANA; Kessler et al. 2009) software which constructs an accurate representation of how a large sample of SNe would appear in a survey, at the image level, by passing synthesized light curves through telescope responses and survey cadences. In this analysis we build a simplified toy model of the DES-SN survey by approximating many of the complex steps involved in SNANA, such as the light

curve generation and fitting. This streamlining reduces computation time, allowing us to explore a broad parameter space while tracing individual stellar populations through a large simulation of galaxies. Results of this qualitative analysis will then be used to inform full scale simulations using SNANA in the future.

We base our simulation on a forward model of host galaxy properties, and test the effects of galaxy and SN age on the light curves, dust extinctions and HRs of SNe Ia. While most SN simulations use empirical relationships to assign SNe to host galaxies that have minimal information (usually stellar mass and SFR), here we use the stellar mass assembly models of Wiseman et al. (2021, hereafter W21) in order to create a library of host galaxies for which the stellar age distribution of the galaxy is known at any given look-back time. We then derive a SN progenitor age distribution for each galaxy in the library, allowing us to trace the age of each simulated SN. In Section 2 we describe the galaxy simulation process and how we use stellar population synthesis models to translate star-formation histories into observables used in K22: stellar mass and rest-frame colour. In Section 3 we describe various models that we use to sample light-curve parameters for the SNe, paying notable attention to the SN progenitor age and host stellar ages. Section 4 presents a validation that the model approximates the observed SN parameter distributions and light-curve parameter–host galaxy relations. In Section 5 we evaluate the models based on their ability to reproduce the Hubble residual trends with SN colour, and in Section 6 we discuss the interpretations and implications of the results.

Where appropriate we assume a flat $\Lambda$CDM cosmological model with $H_0 = 70 \text{ km s}^{-1} \text{ Mpc}^{-1}$ and $\Omega_M = 0.3$. Uncertainties are given as $1\sigma$ confidence intervals.

## 2 GALAXY SIMULATIONS

The main aim of this work is to create a realistic simulation of a population of SNe Ia with an accurate representation of the links between the light-curve properties of the SNe and the global properties of their host galaxies. This library is then used to link the host galaxy properties to Hubble residuals. In this section, we outline the method used to simulate the SN Ia host galaxies and their properties. In Section 3 we describe how we use the host galaxy library to simulate the SNe themselves.

We base our host galaxy simulation on the stellar mass assembly models of W21, which themselves are based on those of Childress, Wolf & Zahid (2014). These simple analytical models of galaxy evolution are derived empirically from observed relationships between stellar mass, SFR and redshift, and take into account mass loss and quenching. The ability of the models to reproduce observed stellar mass and SFR distributions is presented in Childress et al. (2014) and W21.

### 2.1 Updated galaxy models

In the simulations of W21, galaxies evolve following a deterministic set of equations, such that a simulated galaxy represents the average galaxy of a given stellar mass and redshift. While such a prescription is adequate for an overall modelling of SN rates as a function of mass and redshift, it does not result in an accurate distribution of host galaxy star formation rates, which is important when considering host galaxy colours. To mitigate this issue, we update our prescription of galaxy quenching to be more stochastic. We also add bursts of star formation to the models. These two additions are detailed in turn below.
### 2.1.1 Quenching

Each galaxy track in the simulation is represented by 100 implementations. At each time step of the simulation, each galaxy implementation is designated as either ‘non-quenched’ or ‘quenched’. Quenched galaxies remain quenched, while non-quenched galaxies are subject to being quenched according to a quenching probability. This probability is a function of the stellar mass, similar to the average quenching penalty in the Childress et al. (2014) and W21 implementations of the model,

$$p_Q(M_*, z) = \frac{1}{2} \left[ 1 - \operatorname{erf} \left( \frac{\log(M_*) - \log(M_Q(z))}{\sigma_Q} \right) \right] - f_Q, \quad (1)$$

where $f_Q$ is the fraction of the 100 implementations that have already quenched. In order to accurately reproduce the star formation main sequence and quenched fraction, we slightly adjust the reference quenching mass from W21 by 0.2 dex to

$$\log(M_Q(z)/M_\odot) = 10.377 + 0.636z, \quad (2)$$

and the transition scale $\sigma_Q$ from 1.1 to 1.5.

### 2.1.2 Bursts of star formation

A second update to the models of W21 is the addition of random bursts of star formation to the models. We add instantaneous, delta-function bursts that form stars at a rate of half of the main sequence value for the mass and redshift of the galaxy. Bursts occur at a constant probability of 0.05 for every time step in the simulation.

## 2.2 Simulating observable galaxy properties

We build upon the model SFHs by ascribing observable properties to each of the simulated galaxies. In particular, we are interested in the observed and rest frame spectral energy distributions, which are the observed data used to infer global galaxy properties such as stellar mass, SFR, and stellar population age. For the latter, we focus on the rest frame absolute magnitudes in Bessell following Roman et al. (2018), Kelsey et al. (2021), and K22. In the simulation, we link the model SFHs (with known stellar mass, SFR, and age distributions) to synthetic broad-band photometry using a similar method to that widely used in the inference of those observables: we convolve the SFHs with a grid of spectral templates of simple stellar populations (SSPs, Section 2.2.1) of different ages, summing up the individual SSP spectra according to the relative weights given by the SFH to produce a final composite spectrum $S_{gal}$. The luminosities of the combined spectra are scaled by the stellar mass of the galaxy, resulting in an accurate representation of the intrinsic rest-frame galaxy spectra. We add nebular emission from ionized gas (Section 2.2.2), and apply a reddening due to interstellar dust via attenuation (Section 2.2.3). This grid of host galaxies acts as a host galaxy library from which we sample a realistic SN Ia host population in Section 3.

### 2.2.1 Template libraries

The choice of spectral templates affects the relationship between galaxy properties and observables. We use the common Bruzual & Charlot (2003, hereafter BC03) library, and for consistency with the SN host galaxy literature we generate SSPs with a Chabrier (2003) initial mass function. The BC03 code draws upon the Padova 1994 evolutionary tracks (Bertelli et al. 1994). For simplicity, all SSPs are generated at solar metallicity ($Z_\odot$).

### 2.2.2 Nebular emission

In addition to the stellar continuum, our spectral templates incorporate nebular emission lines from H II regions. The BC03 models do not incorporate nebular emission, and thus we add it separately following Boquien et al. (2019). The line strengths are set relative to the ionizing photon flux, and are determined by the ratios at solar metallicity presented in Inoue (2011). The ionization parameter $U$ is set at $\log(U) = -2$.

### 2.2.3 Dust attenuation

To arrive at realistic simulation of galaxy spectra, we apply the effects of dust attenuation to the final scaled spectral templates. We add dust according to the attenuation law of Cardelli, Clayton & Mathis (1989). For each galaxy there are 15 realizations with dust added in steps of equal size in V-band extinction $A_V$ in the range $0 \le A_V \le 1.5$ mag; the final value of $A_V$ for each host is chosen at a later stage along with the SN parameters, such that host has the same reddening $E(B-V)$ as the SN. To calculate the host extinction $A_V$ we follow the results of Salim, Boquien & Lee (2018) and use a total-to-selective extinction $R_V$ of 2.61 for star-forming galaxies with $\log(M/M_\odot) < 9.5$, 2.99 for $9.5 \le \log(M/M_\odot) < 10.5$, and 3.47 for $\log(M/M_\odot) \ge 10.5$. For passive galaxies with $\log(sSFR) < -11$ we follow Salim et al. (2018) and assume the quiescent galaxy $R_V$ of 2.61.

### 2.2.4 Synthetic galaxy photometry

We convert our simulated galaxy spectra into photometric observables by multiplying them with the transmission functions of desired passbands. This is possible with any filter at any wavelength – the galaxy templates extend from the far-ultraviolet to the far-infrared – but here we focus on optical wavelengths. For rest-frame galaxy photometry we calculate absolute magnitudes in the UVBVI filters of Bessell (1990) following K22, and in the observer frame we calculate apparent magnitudes in the DES griz filters (e.g. $m_g^{host}$, $m_f^{host}$; Flaugher et al. 2015) for use in SN survey simulations (e.g. Vincenzi et al. 2021).

The resulting library of simulated galaxy stellar masses and rest-frame colours are shown in Fig. 1, with DES SN Ia hosts shown to illustrate how the simulation covers the necessary stellar-mass-colour parameter space from which to accurately draw SN hosts.

# 3 SUPERNOVAE

We next describe how we simulate a sample of SNe Ia. This includes ensuring a realistic association with their host galaxies, accurate distributions of their light-curve parameters, and the principle effects introduced by observational noise. We discuss each in turn.

## 3.1 SN host galaxy association

We start with our simulated galaxy models from Section 2, and take a flat distribution in stellar mass over the range $7 < \log(M_*/M_\odot) < 12$. We use a SN Ia delay-time distribution (DTD), which describes the rate at which SNe Ia occur as a function of the delay time $\tau_A$ since an episode of star formation, and convolve this DTD with each galaxy’s SFH. In our simulations we use the DTD of W21, which takes a power-law form with index $-1.13$ and takes effect 40 Myr after star formation. The result is a SN Ia rate per year in each model galaxy ($R_G$). We then calculate a volumetric rate of SNe by multiplying $R_G$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Rest frame galaxy U-R colour versus stellar mass for the full simulated galaxy catalogue (blue contours) using the star-formation history models from W21, updated to include better quenching prescription. Magenta points are the observed stellar masses and colours from the DES5YR sample of SN Ia host galaxies (K22), and indicate where SN Ia hosts lie in this parameter space. The vertical dotted line indicates the log (M*/M⊙) = 10 location assumed for mass-related transitions throughout this work, while the horizontal dashed line is located at the corresponding point used to separate galaxy colours, U-R = 1.</figcaption></figure>

by a stellar-mass function ($\phi(M_*, z)$) at the redshift of the simulated galaxy. The combination $R_G \times \phi(M_*, z)$ is then used as a weighting to select a representative sample of SN Ia host galaxies in which SNe can be simulated. We use the redshift-dependent stellar-mass function from ZFOURGE (Tomczak et al. 2014).

The relationship between observed galaxy colour and SN and galaxy ages is shown in Fig. 2. There is a distinction between the host U-R populations for young and old SNe, although there is a significant fraction of old SNe that exist in blue hosts. The relationship between galaxy colour and galaxy age is even clearer, with very little crossover between the old-and-red, and young-and-blue populations.

## 3.2 Light-curve properties

For each simulated SN in a given host galaxy, we then assign ‘observed’ light-curve properties by drawing randomly from statistical distributions. We use the standard SALT2 light-curve fitter (Guy et al. 2007) framework of $x_1$ (light-curve width, or ‘stretch’) and c (SN peak rest-frame optical colour), together with a Tripp (1998)-like standardization relation that relates the inferred distance modulus $\mu_{obs}$ to SN Ia light-curve parameters. In most current SN Ia analyses, this takes the form

$$ \mu_{\text{obs}} = m_B - M_B + \alpha x_1 - \beta c + \mu_{\text{bias}}, \quad (3) $$

where $M_B$ is the peak SN absolute magnitude in the rest-frame B band, $m_B$ is the peak rest-frame apparent magnitude, the $\alpha$ and $\beta$ coefficients parametrize the stretch–luminosity and colour–luminosity relationships, and the $\mu_{bias}$ term is a correction for biases introduced by survey selection effects. Hubble residuals, $\mu_{res}$, can then be calculated using

$$ \mu_{\text{res}} = \mu_{\text{obs}} - \mu_{\text{model}}, \quad (4) $$

where $\mu_{model}$ depends on the cosmological parameters, which in this analysis we hold fixed at our fiducial model.

This simple framework assumes all colour variation in SNe Ia can be captured in a simple ‘$\beta c$’ term. However, we are explicitly interested in this analysis in the effect of dust extinction and reddening on SNe Ia, and thus we also assign each simulated SN a colour excess (i.e. $E(B-V)$) due to dust reddening that we denote $E_{SN}$, and a total-to-selective extinction parameter $R_V$. We next outline the methods used to draw the samples for each light-curve and dust parameter. A summary of the main model variations used in this work is presented in Table 1.

### 3.2.1 SN Ia absolute magnitude

We simulate SNe using a fixed peak B-band absolute magnitude of $M_B = -19.325$ mag. We allow for the possibility of host galaxy ‘steps’ of size $\gamma$, where SNe on either side of a threshold in a given host galaxy property X (such as stellar mass) have a factor of $\Gamma \times \gamma$ added to their simulated absolute magnitude, where $\Gamma = +0.5$ for hosts below the threshold and $\Gamma = -0.5$ for hosts above the threshold $X_{split}$. Thus,

$$ M_B = -19.325 + \Gamma_X \times \gamma_X \text{ mag}, \quad (5) $$

with

$$ \Gamma_X = \begin{cases} 0.5, & X < X_{\text{split}} \\ -0.5, & X \ge X_{\text{split}} \end{cases}, \quad (6) $$

where $X_{step}$ is the location of the step in the given host parameter.

### 3.2.2 SN $x_1$

We include two prescriptions for the SALT2 $x_1$ parameter. The first parametrizes the $x_1$ distribution as an asymmetric Gaussian, while the second assumes two Gaussians dependent on the age of the SN progenitor’s stellar population.

The first $x_1$ population model is used in BS21, which in turn follows Scolnic & Kessler (2016), where $x_1$ is drawn from an asymmetric Gaussian model defined by parameters $\sigma_-$, $\sigma_+$, $\bar{x}_1$,<fnref n="1" /> giving a probability distribution

$$ P(x_1) = \begin{cases} e^{\left[-(x_1 - \bar{x}_1)^2 / 2\sigma_+^2\right]}, & x_1 < \bar{x}_1 \\ e^{\left[-(x_1 - \bar{x}_1)^2 / 2\sigma_-^2\right]}, & x_1 > \bar{x}_1 \end{cases}. \quad (7) $$

The second $x_1$ model builds upon the evolving double-Gaussian of Nicolas et al. (2021) with the $x_1$ probability distribution represented as a Gaussian mixture model – a combination of two Gaussian distributions, one for a ‘young’ SN population and one an ‘old’ population. In the Nicolas et al. (2021) model, environmental age is traced by the local specific SFR (lsSFR). SNe in young (high lsSFR) environments are assumed to belong solely to the old $x_1$ population, while SNe in old (low lsSFR) environments can be drawn from either the old or young $x_1$ modes, with the relative weight of each determined by a mixing probability $a_{x_1}$. Since the fraction of old and young SNe varies with redshift, so does the $x_1$ distribution (see fig. 8 in Nicolas et al. 2021).

Instead of using lsSFR as an observational indicator of progenitor age, our simulations allow us to track the progenitors directly. For

<footnote n="1">¹We use the notation $\bar{x}_1$ instead of the conventional $\mu_{x_1}$ to avoid confusion with the distance modulus $\mu$, and maintain this notation for all parameter means.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Left: Galaxy rest frame U-R colour versus SN progenitor age at the time of explosion for SNe in the simulation. Contours represent lines of equal relative density. The horizontal line represents the colour U-R = 1 used to split observed populations, while the vertical line is the age used to split the young and old SN populations. Right: As left, but galaxy rest frame U-R colour versus galaxy mass weighted mean stellar age. The vertical line shows the 3 Gyr split used for distinguishing young and old galaxies.</figcaption></figure>

each galaxy, the likelihood that a SN exploding at time $t_0$ comes from a progenitor of age $\tau_A$ is

$$P(t_0, \tau_A) = \psi(t_0 - \tau_A)\Phi(\tau_A), \qquad (8)$$

where $\psi(t_0 - \tau_A)$ is the SFH and $\Phi(\tau_A)$ the DTD. We sample SN age $\tau_A$ from this distribution, and subsequently derive $x_1$ from the age by adapting the Nicolas et al. (2021) prescription

$$x_1 \sim \begin{cases} \mathcal{N}(\bar{x}_{1,1}, \sigma_{x_1,1}), & \tau_A < \tau_{A,\text{thresh}} \\ a_{x_1} \times \mathcal{N}(\bar{x}_{1,1}, \sigma_{x_1,1}) + (1 - a_{x_1}) \times \mathcal{N}(\bar{x}_{1,2}, \sigma_{x_1,2}), & \tau_A > \tau_{A,\text{thresh}} \end{cases} \qquad (9)$$

where $\bar{x}_{1,[1,2]}$ and $\sigma_{x_1,[1,2]}$ are the mean and standard deviation of the young and old Gaussians, respectively, $a_{x_1}$ is the mixing coefficient which corresponds to the probability that the progenitor belongs to the young $x_1$ mode, and $\tau_{A, thresh}$ is the age which separates young and old timescales.

### 3.2.3 SN colour

We follow the method of BS21, which itself builds on Jha, Riess & Kirshner (2007), Mandel, Narayan & Kirshner (2011), and Mandel et al. (2017). These model the observed SN Ia colour distribution as a combination of an intrinsic colour $c_{int}$ and the colour excess due to dust reddening $E_{SN}$. $c_{int}$ is normally distributed

$$c_{\text{int}} \sim \mathcal{N}(\bar{c}_{\text{int}}, \sigma_{c_{\text{int}}}), \qquad (10)$$

and $E_{SN}$ is drawn from an exponential distribution following BS21

$$P(E_{\text{SN}}) = \frac{1}{\tau_E} e^{-E_{\text{SN}}/\tau_E}, \qquad (11)$$

where $\tau_E$ is the mean reddening of the SN population.

### 3.2.4 Extinction law $R_V$

Related to the SN dust reddening $E_{SN}$ is the slope of the reddening law $R_V$, such that the absolute V-band extinction $A_V$ is

$$A_V = R_V E_{\text{SN}}, \qquad (12)$$

and for the B-band magnitudes typically used in SN Ia analyses,

$$A_B = R_B E_{\text{SN}} \equiv (R_V + 1) E_{\text{SN}}. \qquad (13)$$

We investigate three models for how $R_V$ varies among SNe. These are inspired by the model of BS21, where $R_V$ follows a step with stellar mass, and motivated by the evidence that the $x_1$ distribution is related to SN or host galaxy age. Our $R_V$ that depend on some host galaxy property Y are:

(i) BS21 ($Y = M_*$): $R_V$ is drawn from a Gaussian distribution,

$$R_V \sim \mathcal{N}(\bar{R}_V, \sigma_{R_V}), \qquad (14)$$

with $\bar{R}_V$ and $\sigma_{R_V}$ taking on different values for host galaxies on either side of $Y_{split}$, a split point in stellar mass;

(ii) Age $R_V$ ($Y = \tau_G$): The same model as BS21, but using mass-weighted mean stellar age of the host galaxy, $\tau_G$, instead of host stellar mass, as the parameter to split $R_V$. $\tau_G$ is linked to both the observed stellar mass and the U-R colour;

(iii) Age $R_V$ Linear ($Y = \tau_A$): $\bar{R}_V$ follows a linear relationship with progenitor age $\tau_A$. We parametrize the relationship by fixing the value of $\bar{R}_V$ for very young ($t_1$) and very old ($t_2$) progenitors as

$$\bar{R}_V = R_V^{t_1} + (\tau_A - t_1) \frac{R_V^{t_2} - R_V^{t_1}}{t_2 - t_1}, \qquad (15)$$

and the standard deviation $\sigma_{R_V}$ following a similar relationship in order to encapsulate having more diverse dust in younger systems:

$$\sigma_{R_V} = \sigma_{R_V}^{t_1} + (\tau_A - t_1) \frac{\sigma_{R_V}^{t_2} - \sigma_{R_V}^{t_1}}{t_2 - t_1}. \qquad (16)$$

### 3.2.5 Mean reddening $\tau_E$

For setting the mean reddening $\tau_E$ we follow BS21, fixing the value separately in low- and high-mass galaxies. We also implement a slight variation where the split is defined on the mass-weighted mean stellar age of the galaxy.
### 3.2.6 Colour-luminosity coefficient $\beta$

As per BS21 we do not assume a universal relationship between SN peak brightness and $c_{int}$. Instead we allow the coefficient relating them, denoted $\beta_{sim}$, to vary between SNe, and draw it from a Gaussian distribution,

$$ \beta_{sim} \sim \mathcal{N}(\bar{\beta}, \sigma_{\beta}), \qquad (17) $$

with mean $\bar{\beta}$ and standard deviation $\sigma_{\beta}$. In this model $\beta_{sim}$ is independent of any host galaxy, progenitor, or dust property.

### 3.2.7 Stretch-luminosity coefficient $\alpha$

For the baseline model we adopt a value of 0.15 for the stretch-luminosity coefficient $\alpha_{sim}$. In cases where a mass or age luminosity step are introduced, the measured value of $\alpha$ can be different to the input value $\alpha_{sim}$. In such cases we vary $\alpha_{sim}$ such that the measured value matches that found by K22. We further investigate the effects of $\alpha$ and its importance in age-based models in Section 6.3.1.

## 3.3 Observables

Given a set of simulated SN Ia light-curve parameters, the synthetic $m_B$ is then calculated by analogy to equation (3), modified to explicitly distinguish dust from SN colour variations, i.e.

$$ m_B = M_B + \mu(z)_{sim} - \alpha_{sim}x_1 + \beta_{sim}c_{int} + E_{SN}(R_V + 1) + \Delta m_B, \qquad (18) $$

where $\mu(z)_{sim}$ is the distance modulus at the simulated redshift $z$ given by the reference cosmology.

Even following the standardization of SNe Ia there remains additional scatter in their observed $M_B$. We describe this in our simulated SNe using the parameter $\sigma_{int}$. We draw an offset ($\Delta M_B$) to $M_B$ from this scatter according to the distribution,

$$ \Delta M_B \sim \mathcal{N}(0, \sigma_{int}). \qquad (19) $$

The value of $\sigma_{int}$ is tuned so that in combination with the rest of the model the overall unexplained scatter in the Hubble residuals matches that of the data, which is of the order ~0.1 mag.

## 3.4 Application to DES5YR data

The models introduced above are generic and intended to approximate the intrinsic universal properties of SNe Ia. Here we describe the data set against which we compare the simulations, and the survey-specific additions to the simulation that are necessary in order to replicate the selection effects of the data set.

We use an identical sample to that presented in K22: namely a photometrically classified SN Ia sample using 5 yr of DES observations (DES5YR), as described in Vincenzi et al. (2021) and Möller et al. (2022), with a redshift limit of $z \le 0.7$. Despite the redshift cut, the sample is magnitude-limited, and thus there are selection effects which we must approximate, which we discuss below.

### 3.4.1 Host galaxy spectroscopic efficiency

To be included in the sample, DES SNe Ia require a spectroscopic redshift of the host galaxy, and the efficiency of obtaining one in the spectroscopic follow-up campaigns is a strong function of host galaxy apparent magnitude, $m_p^{\text{host}}$. Vincenzi et al. (2021) presented the DES

SN spectroscopic-redshift efficiency function $\epsilon_{z_{spec}}(m_p^{\text{host}})$, which we apply to our simulated SN hosts. In the simulation, we assign a probability of successfully obtaining a spectroscopic redshift for each galaxy at random from a Bernoulli distribution with success probability $p = \epsilon_{z_{spec}}(m_p^{\text{host}})$. In this work we do not model further survey efficiencies (e.g. SN detection efficiency; Brout et al. 2019b; Kessler et al. 2019). Therefore, instead of simulating a population of SNe that follows the true intrinsic volumetric rate as a function of redshift, we sample redshift from an ad hoc probability distribution that matches the observed DES5YR redshift distribution. We find that sampling with a probability that increases directly proportional to redshift (up until $z = 0.7$) matches the data well. We validate this $z$ distribution in Section 4.

### 3.4.2 Supernova uncertainties

Without fully simulating SN light curves and their subsequent observation by the survey, we must approximate the uncertainties on our SN light-curve parameters. We simulate uncertainties $\sigma_{m_B}$ on $m_B$ by approximating the relationship between $\sigma_{m_B, \text{obs}}$ and $m_B, \text{obs}$ in the DES5YR data. The uncertainty on $x_1$ ($\sigma_{x_1}$) is then simulated according to $\sigma_{m_B}$ following a linear least-squares fit to the observed DES5YR $\sigma_{m_B}$ and $\sigma_{x_1}$. The same procedure is followed for the colour uncertainty, $\sigma_c$. A full description of the uncertainties is presented in Appendix A.

### 3.4.3 Bias corrections

With simulated brightnesses and light-curve parameters, their respective uncertainties, as well as a redshift distribution that all match the data, the simulation mimics the selection effects and observational biases inherent in the data. We therefore calculate the distance moduli to each simulated SN Ia using an identical method to K22 using the BEAMS with Bias Corrections (BBC; Kessler & Scolnic 2017) framework. BBC follows the method of Marriner et al. (2011), a $\chi^2$ minimization of Hubble residuals binned in redshift. The BEAMS part of BBC is not used, as we do not deal with contamination from core-collapse SNe and all simulated SNe are assigned as SNe Ia with probability 1.0. In BBC, we fit for equation (3). We do not fit for the additional host step terms $\gamma_X$, as we investigate how these terms evolve as a function of SN colour in this work.

The simplest form of bias correction ($\mu_{bias}$; equation (3)) implemented in BBC is dependent on redshift only: a one-dimensional correction ('BBC-1D'). However, Scolnic & Kessler (2016) showed that the asymmetric nature of the $c$ and $x_1$ distributions leads to distance biases as a function of those parameters, which must be corrected as a function of $\alpha$, $x_1$, $\beta$, and $c$ in addition to $z$ (a 'BBC-5D' bias correction). In implementations of such 5D corrections, the treatment of relationships between $x_1$ and host stellar mass must also be carefully accounted for to avoid introducing subtle biases (Smith et al. 2020). Popovic et al. (2021b) introduced a further two dimensions to account for such correlations between light-curves and host galaxies ('BBC-7D'), as well as a method that assumes the BS21 colour model ('BBC-BS20').

In this work we use the 1D (redshift only) bias correction to validate and compare models, and also compare our best models to the data with the BBC-BS20 corrections. We include a $x_1$-host mass correlation in the bias correction simulation as per Vincenzi et al. (2021), and for BBC-BS20 we use the best-fitting $R_V$ values from BS21 in the bias correction simulation. Both approaches require that the simulation adequately includes the same observational biases
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Average magnitude of a one-dimensional (BBC-1D; blue circles) and four-dimensional (BBC-BS20; magenta diamonds) bias correction as a function of redshift for a representative simulation. Points show the median in redshift bins; point sizes are proportional to the number of SNe in each redshift bin, while shaded regions show the 16th–84th percentiles.</figcaption></figure>

<table><caption>**Table 1.** Outline of the different SN Ia models used in this analysis.</caption><thead><tr><th>Model name</th><th>SN x<sub>1</sub><br>distribution</th><th>SN c<br>distribution</th><th>R<sub>V</sub><br>distribution</th></tr></thead><tbody><tr><td>BS21</td><td>SK16</td><td>Gauss + exp (mass)</td><td>Step (stellar mass)</td></tr><tr><td>Age R<sub>V</sub></td><td>N21</td><td>Gauss + exp (age)</td><td>Step (galaxy age)</td></tr><tr><td>Age R<sub>V</sub> Linear</td><td>N21</td><td>Gauss + exp (age)</td><td>Linear (WD age)</td></tr></tbody></table>

and selection effects that have been introduced on to the data, which although not guaranteed has been well validated in the following section. The average and range of bias correction magnitudes for both cases are shown in Fig. 3 and show that, in particular for the BBC-1D case, the corrections are small compared to the step size in the redshift range considered.

### 3.5 Baseline model

For the validation of the model (Section 4), we introduce baseline simulations on top of which the dust and luminosity step parameters are subsequently varied (Section 5). For the BS21 model we fix the $x_1$ and c population parameters similarly to those used in BS21. We use the $x_{1, sim}$ parameters from Scolnic & Kessler (2016), while for $c_{sim}$ we use the values in Table 2. The parameters for both $x_{1, sim}$ and $c_{sim}$ used in our baseline age-based model are also presented in Table 2. Dust parameters $R_V$ are held constant at 2.75, although $R_V$ is not important for the validation presented below. The width-luminosity coefficient is fixed at $\alpha_{sim} = 0.15$ and the colour-luminosity coefficient has mean $\bar{\beta} = 1.98$ and standard deviation $\sigma_{\beta} = 0.35$.

## 4 VALIDATION

Before presenting the results, we describe the metrics used to validate the simulations introduced in Sections 2 and 3. We validate our model by comparing it to the DES5YR SN Ia data presented in K22. A set of qualitative and quantitative validations are outlined as follows. In the following sections, we define the chi-squared statistic on a histogram with both the data and simulation in equivalent bins,

$$ \chi^2 = \sum_i \frac{(N_i^{\text{data}} - N_i^{\text{sim}})^2}{e_i^2}, \qquad (20) $$

where $N_i^{\text{data}}$ is the count of the observed data in the $i$th bin, $N_i^{\text{sim}}$ is the equivalent for the simulation, and $e_i^2$ the error on the data

approximated by $\sqrt{N_i^{\text{data}}}$.

$$ \chi_v^2 = \chi^2 / N_{\text{bins}}. \qquad (21) $$

### 4.1 Host galaxies

We validate the host galaxy simulations by investigating the one-dimensional distributions of redshift, stellar mass ($M_*$), and rest-frame $U-R$ (Fig. 4). The overall trend of the redshift distribution is recovered by construction, as outlined in Section 3.4.1. The distribution of observed $M_*$ is well-replicated. The observed distribution of $U-R$ is bimodal, and this general behaviour is recovered by the model, although the peak of red galaxies in the simulation lies around $U-R \sim 1.5$, compared to the data which peak closer to 1.8. This difference could be due to a number of reasons, such as the simplistic treatment of metallicity and an incomplete implementation of galaxy quenching. While not perfect, we proceed satisfied that both the observed and simulated distributions can be separated at $U-R = 1$. Since the Hubble residual steps tested using this split point, differing shape of the red peak of the distribution should not affect the results significantly.

We further validate this argument by comparing the observed and simulated correlations between $M_*$ and $U-R$ (Fig. 1): the simulations accurately replicate the positive trend between $U-R$ and $M_*$, as well as the characteristic step from majority blue to majority red galaxies as their mass increases beyond $\log(M_*/M_\odot) = 10$. This qualitative validation reinforces trust that our galaxy evolution model reproduces, on a general level, the complex relationships between galaxy stellar mass, age, and colour, and that the DTD-SFH convolution selects a representative sample of SN Ia hosts from the global galaxy population.

### 4.2 Light-curve parameters

We next validate the simulated distributions of SN $x_1$ and c, and their relationships with host galaxy attributes $M_*$ and $U-R$. The simulated distributions are shown in Fig. 5, and the parameters used for the baseline age-based model are presented in Table 2.

The c distribution is well-matched by the ‘Gaussian plus dust reddening’ model. We find the data are described well ($\chi_v^2 = 0.95$) using similar parameters to BS21 for both the intrinsic Gaussian and exponential tail, with the mean of the colour distribution shifting marginally bluer. We find a single value of $\tau_E = 0.135$ describes the data well, which implies there is no strong difference in the mean reddening for SNe in different galaxies, whether they are split by their stellar mass or their age.

Both the asymmetric Gaussian model and Gaussian mixture model provide good fits to the $x_1$ distribution, with the mixture model offering a slightly improved fit around the peak. The best-matching parameters for the mixture model are similar to those from N21, with the only significant difference being the mixture coefficient $a_{x_1}$, for which N21 report $\sim 0.5$ while we find 0.38. In our model, the value of $a_{x_1}$ is degenerate with the choice of age threshold $\tau_A$, an investigation into which is left for future work. There is a degree of overfitting, with the best-fitting model having $\chi_v^2 = 0.65$, due to the large uncertainties on measured $x_1$ in the SN Ia data.

The relationship between $x_1$ and host galaxy properties is shown in Fig. 6. The N21 model implemented in our simulations reproduces
<table><caption>**Table 2.** SN population parameters for the age-based models. Explanation of parameters can be found in Section 3.2.</caption><thead><tr><th>x&#773;<sub>1,1</sub></th><th>&sigma;<sub>x<sub>1,1</sub></sub></th><th>x&#773;<sub>1,2</sub></th><th>&sigma;<sub>x<sub>1,2</sub></sub></th><th>a<sub>x<sub>1</sub></sub></th><th>&tau;<sub>A, thresh</sub></th><th>c&#773;<sub>int</sub></th><th>&sigma;<sub>cint</sub></th><th>&tau;<sub>E, 1</sub></th><th>&tau;<sub>E, 2</sub></th></tr></thead><tbody><tr><td>0.22</td><td>0.61</td><td>-1.22</td><td>0.56</td><td>0.38</td><td>0.75 Gyr</td><td>-0.075</td><td>0.042</td><td>0.135</td><td>0.135</td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Population distributions of observed DES5YR SNe Ia compared to SNe Ia samples simulated in this work, in redshift, host galaxy stellar mass, and host galaxy rest-frame *U*-*R* colour. Counts from the simulation have been scaled to match the total number of DES5YR SNe.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left: distribution of SN SALT2 colour *c* for DES5YR and simulations drawing *c* from an intrinsic Gaussian and an exponential dust distribution. Right: distribution of SN SALT2 stretch *x*<sub>1</sub> for DES5YR and simulations drawing from an asymmetric Gaussian (yellow dashed) as per Scolnic & Kessler (2016) and a double Gaussian based on SN progenitor age (red solid) as per N21.</figcaption></figure>

the general trends in the data, but does not achieve the same strength of relation, despite the one-dimensional $x_1$ distribution being adequately replicated (Fig. 5). The likely cause is the inconsistency between the simulated and observed *U*-*R* distribution, due to the under-representation of passive/red galaxies in the simulation. The underprediction of the slope of the $x_1$ versus $M_*$ relation has consequences in the recovery of the $\alpha$ parameter (Dixon 2021; Rose et al. 2021).

## 5 RESULTS

The principle aim of this work is to explain the observed trends between Hubble residual and SN *c*, with the sample divided by both

host galaxy $M_*$ and host galaxy *U*-*R* colour. In this section, we present the results of simulating a sample of SN Ia host galaxies (Section 2), light-curve properties (Section 3.2), and fitting their Hubble residuals (Section 3.4.3) using different assumptions of the generative models for $x_1$ and $R_V$ and various forms of intrinsic luminosity step.

For each of our model variations, all of the parameters defining the population distributions are held fixed apart from three parameters of interest, for which we simulate a course grid. These parameters are the two extinction law $R_V$ parameters and the luminosity step size. We vary $R_{V,1}$ between [1.5 – 2.75] in steps of 0.25, $R_{V,2}$ in the range [2.5 – 3.75] with steps of 0.25, and $\gamma_X$ between [0 – 0.2] with steps of 0.05. Future analyses will employ inference techniques to precisely
<figure><img src="image_7.png" /><figcaption>Figure 6. Left: Distribution of SN $x_1$ as a function of host galaxy stellar mass for DES5YR data and simulations using the N21 double Gaussian $x_1$ model based on SN progenitor age; Right: as left, but for host galaxy U-R colour.</figcaption></figure>

constrain the best-fitting parameters of these models for their use in cosmological analyses, as has been performed for the BS21 model in Popovic et al. (2021a). The relationship between simulated SN $R_V$ as a function of host galaxy stellar mass and U-R colour is shown in Fig. 7, for a representative simulation with $R_{V,1} = 1.75$, $R_{V,2} = 3.0$. It is clear that an $R_V$ step in age is better recovered by U-R colour than stellar mass, and that stellar mass steps are not particularly well recovered by splitting colour at U-R = 1, highlighting the need for both diagnostics.

Each model is compared to the Hubble residuals measured in the DES5YR data of K22 split by $M_*$ and U-R. A summary of the results of our simulations and comparison to the data is presented in Tables 3 and 4.

Since SNe are simulated stochastically by sampling from a number of interlinked distributions and a numerical host galaxy library, an analytical evaluation of the relationship between Hubble residual and $c$ in the model is not possible. We thus compare the simulated and observed trends in a similar way to the one-dimensional parameter distributions (Section 4.2).

For each bin $i$ in a host galaxy property $X$, we group simulated SNe in the same bins $c_j$ as used for the data. The reduced chi-squared is then

$$ \chi_{i, \mu_{\text{res}}}^2 = \frac{1}{N_{\text{bins}}} \sum_{i}^{X} \sum_{j}^{N_{\text{bins}}} \frac{(\mu_{\text{res, data, i, j}} - \mu_{\text{res, sim, i, j}})^2}{e_{i, j}^2}, \quad (22) $$

where $\mu_{\text{res, data, i, j}}$ and $\mu_{\text{res, sim, i, j}}$ are the weighted mean Hubble residuals in each host bin $i$ and colour bin $j$ for data and simulations, respectively, and $e_{i, j}$ is the standard error on the means of the data.

Each of the following sections present the results for one model of how $R_V$ varies with host or SN properties. For each model we present the data and simulations (Figs 8–10) split by their stellar mass (left-hand panels) and U-R (right-hand panels). The values from our grid search that result in the smallest $\chi_{i, \mu_{\text{res}}}^2$ for each model when split by stellar mass and U-R are presented in Tables 3 and 4, respectively. Unless otherwise stated, parameters reported below are determined using a 1D (BBC-1D) bias correction. Corresponding figures and tables for BBC-BS20 can be found in Appendix B. We stress that the $R_V$ values here are indicative only, since the generative models for the simulated data and bias corrections are not the same. A thorough fitting via the method of Popovic et al. (2021a) is left for future work.

## 5.1 Fixed $R_V$, no luminosity step

We begin by simulating SNe from a baseline model with no luminosity step and $R_V$ fixed at 2.5 across all host galaxies. The resulting Hubble residuals shown in Fig. 8 (panels a and b) display a trend with SN colour $c$, as expected due to the effects of parameter migration (Scolnic & Kessler 2016), but they do not reproduce the differences between low and high stellar mass, or between blue and red host galaxies. To test whether a simple luminosity step as a function of SN age can explain the diverging Hubble residuals, we add a luminosity age step $\gamma_{\tau_A} = 0.2$ mag to the simulation (Fig. 8, panels c and d). Again, the diverging Hubble residuals are not reproduced. Notably the simulated 0.2 mag luminosity step is not recovered in the Hubble residuals, with a step of only around 0.05 mag appearing. The reasons behind this are discussed in Section 6.3.

## 5.2 $R_V$ split on stellar mass (BS21)

Motivated by the clear signal in the data of diverging Hubble residuals as a function of SN $c$, we continue by implementing the model of BS21. Our BS21 simulation is shown in Fig. 9 (upper panels), and qualitatively replicates the trend in the data when split by stellar mass. Simulations with $R_V = 3.0$ in low mass galaxies and 1.75 in high mass hosts match the data well ($\chi_v^2 = 1.63$), values that are consistent with those found by BS21 and Popovic et al. (2021a).

However, this model that describes the mass-split data adequately performs less well when the data are split by host U-R colour ($\chi_v^2 = 2.81$), and neither the difference between Hubble residuals in blue nor red SNe are matched as well. The best match ($\chi_v^2 = 1.70$), shown in the lower panels of Fig. 9, has a steeper extinction law in high mass galaxies with $R_V = 1.5$, as well as the addition of an intrinsic 0.1 mag luminosity step on stellar mass. However, the addition of this step has the effect of drastically reducing the quality of the match when the data are split by stellar mass ($\chi_v^2 = 4.93$).

## 5.3 $R_V$ split on galaxy age

The BS21 model is able to reproduce the observed trends between Hubble residual and SN $c$ when split by stellar mass but the model requires different parameters and an intrinsic luminosity step when matching the data split by host U-R colour. Motivated by a hypothesis that dust parameters could be driven by stellar age rather than galaxy mass, we implement the 'Age $R_V$' model, which is identical to the BS21 model but with $R_V$ described by a step function at a mean
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left: Relationship between total-to-selective dust extinction $R_V$ and host galaxy stellar mass for the three models used in this work. Right: As left, but shown against host galaxy U-R colour.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** SN Ia Hubble residual versus SN colour, with the DES5YR SN Ia sample shown as the data points and the simulations shown as the solid lines. Left-hand panels show SNe split by their host stellar mass at split at log (M/M\*\*) = 10, and right-hand panels show SNe split by their host colour at U-R = 1. Panel (a) shows simulations from the baseline model with $R_V$ fixed across all hosts and no intrinsic luminosity step. Hubble residuals for data and simulations have been measured with BBC-1D. b): as the left panel, but for the data and simulations split at a host galaxy rest-frame U-R colour of U-R = 1. The $\chi^2$ statistic for both comparisons is shown. c): as (a) but for simulations with an additional 0.2 mag intrinsic luminosity step at a SN progenitor age of 0.75 Gyr. d): as (c) but with SNe split by their host galaxy colour.</figcaption></figure>

stellar age of 3 Gyr rather than log (M/M☉) = 10. This model also incorporates the age-based Nicolas et al. (2021) model of $\chi_1$.

As with the mass-based model, the better matching simulations are those with differing $\bar{R}_V$ in young and old galaxies, with larger differences in $\bar{R}_V$ required compared to the BS21 model. When splitting Hubble residuals on stellar mass (Fig. 10 upper panels, Table 3), the favoured models ($\chi^2$ = 1.94) include a similar difference in $R_V$

between old (1.5) and young (2.75) hosts and no intrinsic luminosity step is required. When splitting the data on U-R this model provides a significantly better match than the mass-$R_V$-step model of BS21 (when implemented with no further luminosity step, i.e. Fig. 10(b) has $\chi^2$ = 1.70 compared to Fig. 9(b) with $\chi^2$ = 2.44). The best match to the U-R split data has a smaller difference between $R_V$ values (1.5 and 2.5) plus additional 0.1 mag intrinsic luminosity step
<figure><img src="image_10.png" /><figcaption>**Figure 9.** SN Ia Hubble residual versus SN colour, as in Fig. 8, but for models with $R_V$ that varies with host galaxy stellar mass. (a): Simulations from the BS21 model with $\bar{R}_V = 2.5$ in low mass host galaxies and $\bar{R}_V = 1.5$ in high mass host galaxies. Hubble residuals for data and simulations have been measured with BBC-1D; (b): as (a), but for the data and simulations split host galaxy $U-R = 1$. The model that describes the mass-split data well does not reproduce the full difference between Hubble residuals split on $U-R$. (c) and (d): as (a) and (b) but for the model parameters that best describe the $U-R$ data (right-hand panel): $\bar{R}_V = 3.0$ in low mass host galaxies and $\bar{R}_V = 1.5$ in high mass host galaxies, and an intrinsic luminosity step at $\log(M/M_\odot) = 10$ of size 0.1 mag.</figcaption></figure>

<table><caption>**Table 3.** Dust population parameters for the models presented in this work when compared to DES5YR data split at host stellar mass $\log(M_*/M_\odot) = 10$ (left-hand panels of Figs 9 and 10.) Y represents the host galaxy or SN property used to determine $R_V$, while X is the property on which an intrinsic luminosity step (of magnitude $\gamma_X$; equation 6) is placed.</caption><thead><tr><th>Model Name</th><th>BiasCor</th><th>R&#773;<sub>V,1</sub></th><th>&sigma;<sub>R<sub>V</sub>,1</sub></th><th>R&#773;<sub>V,2</sub></th><th>&sigma;<sub>R<sub>V</sub>,2</sub></th><th>Y</th><th>Y<sub>split</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>X</th><th>X<sub>split</sub></th><th>&gamma;<sub>X</sub></th><th>&chi;<sup>2</sup><sub>V,&mu;<sub>res</sub></sub>(M<sub>*</sub>)</th></tr></thead><tbody><tr><td>BS21</td><td>1D</td><td>1.75</td><td>1.0</td><td>3.0</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>-</td><td>-</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.0</td><td>1.63</td></tr><tr><td>BS21</td><td>1D</td><td>1.75</td><td>1.0</td><td>3.0</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>-</td><td>-</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.0</td><td>2.03</td></tr><tr><td>Age R<sub>V</sub></td><td>1D</td><td>1.5</td><td>1.0</td><td>3.0</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>-</td><td>-</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.0</td><td>1.94</td></tr><tr><td>Age R<sub>V</sub></td><td>1D</td><td>1.5</td><td>1.0</td><td>2.75</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>-</td><td>-</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.0</td><td>1.86</td></tr><tr><td>Age R<sub>V</sub> Linear</td><td>1D</td><td>1.5</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>A</sub></td><td>-</td><td>0.1 Gyr</td><td>10 Gyr</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.15</td><td>2.27</td></tr></tbody></table>

<table><caption>**Table 4.** Dust population parameters for the models presented in this work when compared to DES5YR data split at host $U-R = 1$ (right-hand panels of Figs 9 and 10.).</caption><thead><tr><th>Model Name</th><th>BiasCor</th><th>R&#773;<sub>V,1</sub></th><th>&sigma;<sub>R<sub>V</sub>,1</sub></th><th>R&#773;<sub>V,2</sub></th><th>&sigma;<sub>R<sub>V</sub>,2</sub></th><th>Y</th><th>Y<sub>split</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>X</th><th>X<sub>split</sub></th><th>&gamma;<sub>X</sub></th><th>&chi;<sup>2</sup><sub>V,&mu;<sub>res</sub>(U-R)</sub></th></tr></thead><tbody><tr><td>BS21</td><td>1D</td><td>1.5</td><td>1.0</td><td>3.0</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>-</td><td>-</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.05</td><td>1.70</td></tr><tr><td>BS21</td><td>1D</td><td>1.5</td><td>1.0</td><td>3.0</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>-</td><td>-</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.10</td><td>1.59</td></tr><tr><td>Age R<sub>V</sub></td><td>1D</td><td>1.75</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>-</td><td>-</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.15</td><td>1.70</td></tr><tr><td>Age R<sub>V</sub></td><td>1D</td><td>1.5</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>-</td><td>-</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.1</td><td>1.65</td></tr><tr><td>Age R<sub>V</sub> Linear</td><td>1D</td><td>1.5</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>A</sub></td><td>-</td><td>0.1 Gyr</td><td>10 Gyr</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.20</td><td>1.99</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** SN Ia Hubble residual versus SN colour, as in Fig. 9 but for the model with $R_V$ changing with galaxy age. (a): The model that best fits the SN data when split by stellar mass includes $\bar{R}_V = 3.0$ in young hosts, $\bar{R}_V = 1.5$ in old hosts and no intrinsic luminosity step. (b): as (a), but showing SNe split at host galaxy $U-R = 1$. (c) and (d): as for (a) and (b) but for the best-fitting model to the data split by host $U-R$ colour: $\bar{R}_V = 2.5$ in young hosts, $\bar{R}_V = 1.75$ in old hosts, and a 0.15 mag intrinsic luminosity step at a progenitor age of 0.75 Gyr.</figcaption></figure>

on stellar mass, or an even smaller $R_V$ difference (1.75 and 2.5) and a 0.15 mag step on stellar age (Fig. 10, lower panels). As described in the results of the previous models, a much smaller step is evident in the Hubble residuals than the age step input into the SNe. We address this effect in Section 6.3.

## 5.4 $R_V$ linear with SN age

Modelling $R_V$ as a linear function of SN progenitor age is not as successful as the previous step-based models, with consistently higher $\chi^2$ (Tables 3 and 4). This $R_V$ model is the only one that also requires an age-luminosity step to explain both the mass-split and $U-R$ split data, with $\gamma_X \ge 0.15$ mag.

# 6 DISCUSSION

The results of our comparisons of simulations to data support the notion that the SN colour–luminosity relation in SNe Ia is not linear with SN colour, consistent with the model of BS21 where the dust extinction law along the line-of-sight to SNe Ia correlates with the global properties of host galaxies. The reason for, and best tracer of, the difference in $R_V$ is less clear and is discussed in the following sections. As with the previous section, the discussion here focuses on BBC-1D bias corrections, but the results are consistent when using BBC-BS20 bias corrections.

## 6.1 $R_V$ in young, low mass galaxies

We find no significant difference in the fit quality when the $R_V$ changes with stellar mass or stellar age, the two of which are themselves strongly correlated. When implementing the BS21 model we find $\bar{R}_{V,2}$ is broadly consistent with the results from BS21 and Popovic et al. (2021a) when splitting on stellar mass. When $R_V$ changes with stellar age the simulations are more consistent with the data across both stellar mass and $U-R$, but $2.5 \le \bar{R}_{V,2} \le 3.0$ are generally smaller than the 3.0 in Popovic et al. (2021a). These values for $R_V$ consistent with those typically measured in SN Ia light curves and spectra for similar hosts, with Cikota, Deustua & Marleau (2016) finding $R_V$ values of $2.71 \pm 1.58$ in spiral galaxies.

## 6.2 $R_V$ in old, high mass galaxies

In all of the models trialled the $R_V$ value in high mass/old galaxies ($\bar{R}_{V,1}$) is required to be significantly lower than that for low mass/young galaxies in order to replicate the trends of Hubble residual against SN colour. Generally, lower values are preferred when the data are split on $U-R$. Values in the range $1.5 \le \bar{R}_{V,1} \le 1.75$ indicate the dust in older, more massive SN hosts is composed mainly of small grains leading to a steep extinction curve. How this observation relates to studies of the general galaxy population is complicated by different relationships seen between dust and different galaxy descriptors. While Reddy et al. (2018) and Salim et al. (2018) show that $R_V$ increases with stellar mass for star-forming galaxies, Salim
et al. (2018) also show that $R_V$ is lower in passive (and starburst) galaxies than those on the star forming main sequence, and that the $R_V$ of star-forming galaxies is more closely related to its overall extinction $A_V$ than its stellar mass. However, the quiescent galaxies in Salim et al. (2018) still have an average $R_V$ of 2.61, far greater than that found here and indicating that the SNe do not trace the average dust composition in the galaxy. We note however that the SN $R_V$ measurements are somewhat dependent on the assumed SALT2 colour law, rendering direct comparisons less informative.

## 6.3 The effect of intrinsic luminosity steps

While models of an $R_V$ step in stellar mass or galaxy age provide reasonable matches to the Hubble residuals split by either stellar mass or host $U-R$ colour, they do not fully match both host properties simultaneously – different $R_V$ steps are required to match mass and to match $U-R$. Furthermore, K22 show that there are residual steps in the data in one parameter (e.g. stellar mass) if the effect of another (e.g. $U-R$) is subtracted. We find that in general the step seen in our simulated HRs is smaller in $U-R$ than in $M_*$, but that step is rarely fully removed using $R_V$ only and that adding an intrinsic luminosity step is required. However, such a step introduces a number of its own issues.

Firstly we consider adding a step based on the SN age. We find that adding a 0.2 mag step has only a marginal effect on the HRs – instead, the best-fitting value $\alpha$ is much smaller than the simulated value $\alpha_{sim}$. This effect is caused in this model by the dependence on age of both the intrinsic step and $x_1$, which act on the SN luminosity in the opposite direction to each other. Such an affect has been observed before (e.g. Rose et al. 2021), who suggest that steps $\gamma_X$ should be fit simultaneously with the other nuisance parameters. However, with the step size evolving with SN colour, this is not trivial.

An intrinsic mass step is affected much less absorbed by $\alpha$, with a 0.1 mag step able to explain the $U-R$ data for most $R_V$ model and bias correction combinations (e.g. Fig. 9). However, this means that the discrepancy between models with an intrinsic mass step and the data split by mass is much larger than when the intrinsic step is on age and no particularly good combinations are found.

### 6.3.1 Varying the width-luminosity coefficient $\alpha$

Neither an intrinsic luminosity step based on SN age or host mass are able to simultaneously model the data split by stellar mass and by $U-R$. One logical explanation for this is that the relationships between SN age, $x_1$, and $\alpha$ are linked with different strengths to the host parameters $M_*$ and $U-R$. The $U-R$ colour traces age more directly than $M_*$, and with the bimodal N21 model $x_1$ is strongly linked to SN age, meaning intrinsic steps in luminosity are washed out by being absorbed into the alpha parameter, and this occurs more strongly in the age-like $U-R$ than $M_*$ which is less related to SN age. To counter this effect, we hypothesize that the two $x_1$ populations could follow different intrinsic values of $\alpha$ while simultaneously having different mean absolute magnitudes.

To test the two-$\alpha$ hypothesis we run a further set of simulations. We fix $R_V$ to those from the best Age-$R_V$ model as measured with BBC-1D and as compared to data split by mass: ($R_{V,1} = 1.5$, $R_{V,2} = 3.0$). We simulate over a new grid, varying $\alpha_{sim}$ from 0.05 to 0.30 for the young and old populations, as well as the intrinsic luminosity step $\gamma_X$. The combinations of $\alpha_{sim}$ and $\gamma_X$ that we investigate are constrained by the fitted $\alpha$. For $X$, we trial both mass and age steps between zero and 0.25 mag.

The results of varying $\alpha_{sim}$ are summarized in Fig. 11, for which the models included a 0.1 mag step on $M_*$. Varying $\alpha_{sim}$ for the young and old SNe makes no difference to the inferred step in either $M_*$ or $U-R$ and does not solve the discrepancy between the size of recovered steps – whatever the values of $\alpha_{sim}$, the resulting HRs are split much more strongly by $M_*$ than by $U-R$. The same results were found when including a 0.2 mag luminosity step on SN age. These simulations thus indicate a universal value for $\alpha$ even if the SNe are divided into two populations of $x_1$. Such a lack of bimodality for $\alpha$ is consistent with previous results, e.g. Sullivan et al. (2011).

### 6.3.2 Reconciling mass and $U-R$ steps

With multiple $\alpha$ values improbable, we investigate whether the discrepancy is caused by limitations in our model. For example, the smooth SFHs lead to a single Gaussian distribution of $U-R$ while the data show bimodality, while the simulation also slightly underestimates the number of hosts in the stellar mass range $9 \le \log(M/M_{\odot})$. Although Fig. 1 shows that the models do reproduce the $M_*$-$U-R$ relation in general, there may be subtle second order effects introducing inconsistencies for individual simulated SNe. Such effects include the treatment of metallicity which is fixed at solar in the simulation but evolves over both the mass and redshift range of the data (e.g. Tremonti et al. 2004; Zahid et al. 2014), the lack of bursts of star formation in the simulation, the choice of BC03 templates for the simulation, and incorrect modelling of survey selection effects. These effects will be investigated in subsequent work.

A second possibility is that the inconsistencies between the model and data are caused by physical effects that have not been fully incorporated or accurately implemented in the model. The relationships between the $R_V$ along a SN sight line and its host galaxy properties are evidently more complicated that a simple step function of galaxy age or stellar mass – when dealing with integrated galaxy attenuation, $R_V$ is inversely correlated with $A_V$, which we have not included in the model for SN $R_V$. Meanwhile, $A_V$ correlates strongly with stellar mass (Zahid et al. 2013), meaning the $R_V$ of galaxies also increases with stellar mass (Salim et al. 2018). On the contrary, the results of BS21, Popovic et al. (2021a), and this work, indicate that $R_V$ on SN sightlines decreases with stellar mass, suggesting that SN sightlines are systematically different than those of the integrated galaxy light. Meanwhile, Salim et al. (2018) also show that galaxies with old stellar populations have lower $R_V$ than those with young populations, which is in line with the age – $R_V$ steps found here. Because our the models track individual stellar populations, it is possible to use our method to compare the predictions to measurements from integral field spectroscopy (e.g. Galbany et al. 2014, 2016, 2018) to measure dust properties in the very local environments of SNe, which we defer to future work.

Our model has no implementation of a relation between galaxy mass and metallicity. Metallicity is strongly correlated with stellar mass and more weakly with SFR (e.g. Yates, Kauffmann & Guo 2012), and is tied to the SFH (Bellstedt et al. 2020) as well as affecting the strengths of nebular emission lines (e.g. Kewley & Ellison 2008), themselves affecting the $U-R$ colour. If there is an intrinsic correlation between metallicity and SN luminosity (e.g. Höflich, Wheeler & Thielemann 1998; Kasen, Röpke & Woosley 2009; Moreno-Raya et al. 2016) it may affect the observed mass and age steps in different amounts, but be hidden amongst the more dominant $R_V$ effects.
MATPLOTLIB (Hunter 2007), SCIPY (Virtanen et al. 2020), and PANDAS (McKinney 2010).

## ACKNOWLEDGEMENTS

We thank the referee Mickael Rigault for comments which significantly improved this paper. P.W. acknowledges support from the Science and Technology Facilities Council (STFC) grant ST/R000506/1. MS acknowledges support from EU/FP7-ERC grant 615929. LK thanks the UKRI Future Leaders Fellowship for support through the grant MR/T01881X/1. LG acknowledges financial support from the Spanish Ministerio de Ciencia e Innovación (MCIN), the Agencia Estatal de Investigación (AEI) 10.13039/501100011033, and the European Social Fund (ESF) 'Investing in your future' under the 2019 Ramón y Cajal program RYC2019-027683-I and the PID2020-115253GA-100 HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC) under the PIE project 20215AT016, and the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement number 759194 - USNAC).

Funding for the DES Projects has been provided by the U.S. Department of Energy, the U.S. National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute of Cosmological Physics at the University of Chicago, the Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft, and the Collaborating Institutions in the Dark Energy Survey.

The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of Cambridge, Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NSF's NOIRLab, the University of Nottingham, The Ohio State University, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, Texas A&M University, and the OzDES Membership Consortium.

Based in part on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2012B-0001; PI: J. Frieman), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.

The DES data management system is supported by the National Science Foundation under Grant Numbers AST-1138766 and AST-1536171. The DES participants from Spanish institutions are partially supported by MICINN under grants ESP2017-89838,

PGC2018-094773, PGC2018-102021, SEV-2016-0588, SEV-2016-0597, and MDM-2015-0509, some of which include ERDF funds from the European Union. IFAE is partially funded by the CERCA program of the Generalitat de Catalunya. Research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Program (FP7/2007-2013) including ERC grant agreements 240672, 291329, and 306478. We acknowledge support from the Brazilian Instituto Nacional de Ciência e Tecnologia (INCT) do e-Universo (CNPq grant 465376/2014-2).

This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.

## DATA AVAILABILITY

The light curves for the DES-SN photometric SN Ia catalogue and associated host galaxy data will be made available as part of the DES5YR SN cosmology analysis at https://des.ncsa.illinois.edu/releases/sn.

## REFERENCES

Astropy Collaboration et al., 2018, AJ, 156, 123

Bellstedt S. et al., 2020, MNRAS, 498, 5581

Bertelli G. et al., 1994, A&AS, 106, 275

Bessell M. S., 1990, PASP, 102, 1181

Betoule M. et al., 2014, A&A, 568, A22

Boquien M., Burgarella D., Roehlly Y., Buat V., Ciesla L., Corre D., Inoue A. K., Salas H., 2019, A&A, 622, A103

Briday M. et al., 2022, A&A, 657, A22

Brout D., Scolnic D., 2021, ApJ, 909, 26

Brout D. et al., 2019a, ApJ, 874, 106

Brout D. et al., 2019b, ApJ, 874, 150

Brout D. et al., 2022, preprint (arXiv:2202.04077)

Bruzual G., Charlot S., 2003, MNRAS, 344, 1000 (BC03)

Cardelli J. A., Clayton G. C., Mathis J. S., 1989, ApJ, 345, 245

Chabrier G., 2003, PASP, 115, 763

Childress M. et al., 2013, ApJ, 770, 108

Childress M. J., Wolf C., Zahid H. J., 2014, MNRAS, 445, 1898

Cikota A., Deustua S., Marleau F., 2016, ApJ, 819, 152

DES Collaboration, 2018, Phys. Rev. Lett., 122, 171301

Dixon S., 2021, PASP, 133, 054501

Dixon M. et al., 2022, preprint (arXiv:2206.12085)

Flaugher B. et al., 2015, AJ, 150, 150

Galbany L. et al., 2014, A&A, 572, 22

Galbany L. et al., 2016, MNRAS, 455, 4087

Galbany L. et al., 2018, ApJ, 855, 107

Gallagher J. S. et al., 2005, ApJ, 634, 210

Gandhi P. J., Wetzel A., Hopkins P. F., Shappee B. J., Wheeler C., Faucher-Giguère C.-A., 2022, preprint (arXiv:2202.10477)

Guy J. et al., 2007, A&A, 466, 11

Harris C. R. et al., 2020, Nature, 585, 357

Höflich P., Wheeler J., Thielemann F., 1998, ApJ, 495, 617

Hunter J. D., 2007, Comp. Sci. Eng., 9, 99

Inoue A. K., 2011, MNRAS, 415, 2920

Jha S., Riess A. G., Kirshner R. P., 2007, ApJ, 659, 122

Kasen D., Röpke F. K., Woosley S. E., 2009, Nature, 460, 869

Kelly P. L., Hicken M., Burke D. L., Mandel K. S., Kirshner R. P., 2010, ApJ, 715, 743

Kelsey L. et al., 2021, MNRAS, 501, 4861

Kelsey L., Sullivan M., Wiseman P., Smith M., Vincenzi M., Scolnic D., Brout D., 2022, preprint (arxiv:2208.01357; K22)

Kessler R., Scolnic D., 2017, ApJ, 836, 56
Kessler R. et al., 2009, ApJS, 185, 32

Kessler R. et al., 2019, MNRAS, 485, 1171

Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

Lampeitl H. et al., 2010, MNRAS, 401, 2331

McKinney W., 2010, in van der Walt S., Millman J., eds, Proc. of the 9th Python in Science Conf, p. 56

Mandel K. S., Narayan G., Kirshner R. P., 2011, ApJ, 731, 120

Mandel K. S., Scolnic D. M., Shariff H., Foley R. J., Kirshner R. P., 2017, ApJ, 842, 93

Marriner J. et al., 2011, ApJ, 740, 72

Möller A. et al., 2022, MNRAS, 514, 5159

Moreno-Raya M. E., Mollá M., López-Sánchez Á. R., Galbany L., Vílchez J. M., Rosell A. C., Domínguez I., 2016, ApJ, 818, L19

Nicolas N. et al., 2021, A&A, 649, A74

Perlmutter S. et al., 1999, ApJ, 517, 565

Phillips M. M., 1993, ApJ, 413, L105

Popovic B., Brout D., Kessler R., Scolnic D., 2021a, preprint (arXiv:2112.04456)

Popovic B., Brout D., Kessler R., Scolnic D., Lu L., 2021b, ApJ, 913, 49

Pskovskii I. P., 1977, SvA, 21, 675

Reddy N. A. et al., 2018, ApJ, 853, 56

Riess A. G., Press W. H., Kirshner R. P., 1996, ApJ, 473, 88

Riess A. G. et al., 1998, AJ, 116, 1009

Rigault M. et al., 2020, A&A, 644, A176

Roman M. et al., 2018, A&A, 615, A68

Rose B. M., Rubin D., Strolger L., Garnavich P. M., 2021, ApJ, 909, 28

Salim S., Boquien M., Lee J. C., 2018, ApJ, 859, 11

Scolnic D., Kessler R., 2016, ApJ, 822, L35

Scolnic D. M. et al., 2018, ApJ, 859, 101

Scolnic D. et al., 2021, preprint (arXiv:2112.03863)

Smith M. et al., 2020, MNRAS, 494, 4426

Sullivan M. et al., 2010, MNRAS, 406

Sullivan M. et al., 2011, ApJ, 737, 102

Tomczak A. R. et al., 2014, ApJ, 783, 85

Tremonti C. A. et al., 2004, AJ, 613, 898

Tripp R., 1998, A&A, 331, 815

Vincenzi M. et al., 2021, MNRAS, 505, 2819

Virtanen P. et al., 2020, Nature Methods, 17, 261

Wiseman P. et al., 2021, MNRAS, 506, 3330 (W21)

Yates R. M., Kauffmann G., Guo Q., 2012, MNRAS, 422, 215

Zahid H. J., Yates R. M., Kewley L. J., Kudritzki R. P., 2013, ApJ, 763, 92

Zahid H. J., Dima G. I., Kudritzki R.-P., Kewley L. J., Geller M. J., Hwang H. S., Silverman J. D., Kashino D., 2014, ApJ, 791, 130

<figure><img src="image_12.png" /><figcaption>**Figure A1.** The relationship between peak B-band brightness $m_B$ and its uncertainty, $\sigma_{m_B}$. Points are from the DES5YR data and the relationship used in equation (A2).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure A2.** Relationships between $\sigma_{m_B}$ and $\sigma_{x_1}$ (upper), and $\sigma_c$ (lower). Points are from the DES5YR data and the relationship used in equation (A2).</figcaption></figure>

## APPENDIX A: SIMULATING UNCERTAINTIES

We simulate uncertainties $\sigma_{m_B}$ on $m_B$ by approximating the relationship between $\sigma_{m_B,obs}$ and $m_B,_{obs}$ in the observed data set of K22.

$$ \sigma_{m_B} \sim \text{Max} \left\{ \frac{N(\bar{\sigma}_{m_B}, \sigma_{\sigma_{m_B}})}{0.025} \right\}, \qquad (A1) $$

where

$$ \bar{\sigma}_{m_B} = 10^{(0.395(m_B - 1.5) - 10)} + 0.03, \qquad (A2) $$

where $\bar{\sigma}_{m_B}$ is represented by the dashed line in Fig. A1. The scatter on the uncertainties $\sigma_{\sigma_{m_B}}$ also increases as a function of $m_B$:

$$ \sigma_{\sigma_{m_B}} = \text{Max} \left\{ \begin{array}{l} 0.003(m_B - 20) \\ 0.003 \end{array} \right\}. \qquad (A3) $$

The uncertainty $\sigma_{m_B}$ is calculated after the intrinsic $m_B$ has been adjusted for $x_1$ and $c$ via equation (18).

Uncertainties on $x_1$ and $c$ are estimated by linear least squares fits to the DES5YR $\sigma_{x_1} - \sigma_{m_B}$ and $\sigma_c - \sigma_{m_B}$ relations, respectively. The observed and fitted relationships between the uncertainties are shown in Figs A1 and A2. We use these uncertainties to add noise to $m_B, x_1$, and $c$ by drawing them randomly from Gaussian distributions centred at 0 as per equation (19).

## APPENDIX B: BBC-BS20 BIAS CORRECTIONS

Here we present results in the same way as Section 5 but with Hubble residuals measured using a BBC-BS20 bias correction, which inherently assumes $R_V$ changing with stellar mass. The models shown in Figs B1–B3 correspond to the best match parameters when splitting on stellar mass (blue) and $U-R$ (green), which are presented in Tables B1 and B2, respectively.
<figure><img src="image_14.png" /><figcaption>**Figure B1.** BS21 models as Fig. 9 but where Hubble residuals for data and simulations have been measured with BBC-BS20. (a): Simulations from the BS21 model with $\bar{R}_V = 3.25$ in low mass host galaxies and $\bar{R}_V = 2.5$ in high mass host galaxies; (b): as a, but for the data and simulations split host galaxy $U-R=1$. (c) and (d): as (a) and (b) but for the model parameters that best describe the $U-R$ data (right-hand panel): $\bar{R}_V = 3.25$ in low mass host galaxies and $\bar{R}_V = 2.25$ in low mass host galaxies, and an intrinsic luminosity step at $\log (M/M_{\odot}) = 10$ of size 0.1 mag.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure B2.** Age-$R_V$ split models as Fig. 10, but where Hubble residuals for data and simulations have been measured with BBC-BS20.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure B3. As Fig. B3 but for the models best-fitting when introducing an intrinsic stellar mass step instead of a SN age step.</figcaption></figure>

<table><caption>Table B1. Dust population parameters for the models presented in this work when compared to DES5YR data split at host log (M*/M⊙) = 10 and fit with BBC-BS20 (left-hand panels of Figs B1–B3.).</caption><thead><tr><th>Model Name</th><th>BiasCor</th><th>R&#773;<sub>V,1</sub></th><th>&sigma;<sub>R<sub>V</sub>,1</sub></th><th>R&#773;<sub>V,2</sub></th><th>&sigma;<sub>R<sub>V</sub>,2</sub></th><th>Y</th><th>Y<sub>split</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>X</th><th>X<sub>split</sub></th><th>&gamma;X</th><th>&chi;<sup>2</sup><sub>V,&mu;<sub>res</sub></sub>(U &minus; R)</th></tr></thead><tbody><tr><td>BS21</td><td>BS20</td><td>2.5</td><td>1.0</td><td>3.25</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>&mdash;</td><td>&mdash;</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.00</td><td>1.16</td></tr><tr><td>BS21</td><td>BS20</td><td>2.25</td><td>1.0</td><td>3.0</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>&mdash;</td><td>&mdash;</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.05</td><td>1.23</td></tr><tr><td>Age R<sub>V</sub></td><td>BS20</td><td>1.75</td><td>1.0</td><td>3.25</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>&mdash;</td><td>&mdash;</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.0</td><td>1.11</td></tr><tr><td>Age R<sub>V</sub></td><td>BS20</td><td>1.75</td><td>1.0</td><td>3.25</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>&mdash;</td><td>&mdash;</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.0</td><td>1.11</td></tr><tr><td>Age R<sub>V</sub> Linear</td><td>BS20</td><td>1.5</td><td>1.0</td><td>3.75</td><td>1.0</td><td>&tau;<sub>A</sub></td><td>&mdash;</td><td>0.1 Gyr</td><td>10 Gyr</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.20</td><td>1.40</td></tr></tbody></table>

<table><caption>Table B2. Dust population parameters for the models presented in this work when compared to DES5YR data split at host U-R = 1 and fit with BBC-BS20 (right-hand panels of Figs B1–B3.).</caption><thead><tr><th>Model Name</th><th>BiasCor</th><th>R&#773;<sub>V,1</sub></th><th>&sigma;<sub>R<sub>V</sub>,1</sub></th><th>R&#773;<sub>V,2</sub></th><th>&sigma;<sub>R<sub>V</sub>,2</sub></th><th>Y</th><th>Y<sub>split</sub></th><th>Y<sub>1</sub></th><th>Y<sub>2</sub></th><th>X</th><th>X<sub>split</sub></th><th>&gamma;X</th><th>&chi;<sup>2</sup><sub>V,&mu;<sub>res</sub></sub>(U &minus; R)</th></tr></thead><tbody><tr><td>BS21</td><td>BS20</td><td>2.25</td><td>1.0</td><td>3.25</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>&mdash;</td><td>&mdash;</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.1</td><td>0.80</td></tr><tr><td>BS21</td><td>BS20</td><td>1.75</td><td>1.0</td><td>3.75</td><td>1.0</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>&mdash;</td><td>&mdash;</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.05</td><td>0.76</td></tr><tr><td>Age R<sub>V</sub></td><td>BS20</td><td>2.0</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>&mdash;</td><td>&mdash;</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.15</td><td>0.89</td></tr><tr><td>Age R<sub>V</sub></td><td>BS20</td><td>1.5</td><td>1.0</td><td>2.5</td><td>1.0</td><td>&tau;<sub>G</sub></td><td>3 Gyr</td><td>&mdash;</td><td>&mdash;</td><td>M<sub>*</sub></td><td>10<sup>10</sup> M<sub>&odot;</sub></td><td>0.05</td><td>0.82</td></tr><tr><td>Age R<sub>V</sub> Linear</td><td>BS20</td><td>1.5</td><td>1.0</td><td>3.75</td><td>1.0</td><td>&tau;<sub>A</sub></td><td>&mdash;</td><td>0.1 Gyr</td><td>10 Gyr</td><td>&tau;<sub>A</sub></td><td>0.75 Gyr</td><td>0.20</td><td>1.22</td></tr></tbody></table>

<footnote n="1">¹School of Physics and Astronomy, University of Southampton, Southampton, SO17 1BJ, UK</footnote>

<footnote n="2">²Department of Physics, Duke University, Durham, NC 27708, USA</footnote>

<footnote n="3">³Institute of Cosmology and Gravitation, University of Portsmouth, Portsmouth PO1 3FX, UK</footnote>

<footnote n="4">⁴Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="5">⁵School of Mathematics and Physics, University of Queensland, Brisbane, QLD 4072, Australia</footnote>

<footnote n="6">⁶Institut d'Estudis Espacials de Catalunya (IEEC), E-08034 Barcelona, Spain</footnote>

<footnote n="7">⁷Institute of Space Sciences (ICE, CSIC), Campus UAB, Carrer de Can Magrans, s/n, E-08193 Barcelona, Spain</footnote>

<footnote n="8">⁸Centre for Gravitational Astrophysics, College of Science, The Australian National University, ACT 2601, Australia</footnote>

<footnote n="9">⁹The Research School of Astronomy and Astrophysics, Australian National University, ACT 2601, Australia</footnote>

<footnote n="10">¹⁰Centre for Astrophysics & Supercomputing, Swinburne University of Technology, Victoria 3122, Australia</footnote>

<footnote n="11">¹¹Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, IP2I Lyon / IN2P3, IMR 5822, F-69622 Villeurbanne, France</footnote>

<footnote n="12">¹²Laboratório Interinstitucional de e-Astronomia - LIneA, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil</footnote>

<footnote n="13">¹³Fermi National Accelerator Laboratory, P. O. Box 500, Batavia, IL 60510, USA</footnote>

<footnote n="14">¹⁴Department of Physics, University of Michigan, Ann Arbor, MI 48109, USA</footnote>

<footnote n="15">¹⁵CNRS, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="16">¹⁶Sorbonne Universités, UPMC Univ Paris 06, UMR 7095, Institut d'Astrophysique de Paris, F-75014, Paris, France</footnote>

<footnote n="17">¹⁷University Observatory, Faculty of Physics, Ludwig-Maximilians-Universität, Scheinerstr. 1, D-81679 Munich, Germany</footnote>

<footnote n="18">¹⁸Department of Physics & Astronomy, University College London, Gower Street, London WC1E 6BT, UK</footnote>
<footnote n="19">¹⁹Kavli Institute for Particle Astrophysics & Cosmology, P. O. Box 2450, Stanford University, Stanford, CA 94305, USA</footnote>

<footnote n="20">²⁰SLAC National Accelerator Laboratory, Menlo Park, CA 94025, USA</footnote>

<footnote n="21">²¹Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain</footnote>

<footnote n="22">²²Universidad de La Laguna, Dpto. Astrofísica, E-38206 La Laguna, Tenerife, Spain</footnote>

<footnote n="23">²³Center for Astrophysical Surveys, National Center for Supercomputing Applications, 1205 West Clark St., Urbana, IL 61801, USA</footnote>

<footnote n="24">²⁴Department of Astronomy, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA</footnote>

<footnote n="25">²⁵Institut de Física d'Altes Energies (IFAE), The Barcelona Institute of Science and Technology, Campus UAB, 08193 Bellaterra (Barcelona) Spain</footnote>

<footnote n="26">²⁶Astronomy Unit, Department of Physics, University of Trieste, via Tiepolo 11, I-34131 Trieste, Italy</footnote>

<footnote n="27">²⁷INAF-Osservatorio Astronomico di Trieste, via G. B. Tiepolo 11, I-34143 Trieste, Italy</footnote>

<footnote n="28">²⁸Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy</footnote>

<footnote n="29">²⁹Hamburger Sternwarte, Universität Hamburg, Gojenbergsweg 112, D-21029 Hamburg, Germany</footnote>

<footnote n="30">³⁰Department of Physics, IIT Hyderabad, Kandi, Telangana 502285, India</footnote>

<footnote n="31">³¹Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Dr., Pasadena, CA 91109, USA</footnote>

<footnote n="32">³²Institute of Theoretical Astrophysics, University of Oslo, P.O. Box 1029 Blindern, NO-0315 Oslo, Norway</footnote>

<footnote n="33">³³Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL 60637, USA</footnote>

<footnote n="34">³⁴Instituto de Física Teórica UAM/CSIC, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="35">³⁵Department of Physics and Astronomy, University of Pennsylvania, Philadelphia, PA 19104, USA</footnote>

<footnote n="36">³⁶Observatório Nacional, Rua Gal. José Cristino 77, Rio de Janeiro, RJ - 20921-400, Brazil</footnote>

<footnote n="37">³⁷Santa Cruz Institute for Particle Physics, Santa Cruz, CA 95064, USA</footnote>

<footnote n="38">³⁸Center for Cosmology and Astro-Particle Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="39">³⁹Department of Physics, The Ohio State University, Columbus, OH 43210, USA</footnote>

<footnote n="40">⁴⁰Institució Catalana de Recerca i Estudis Avançats, E-08010 Barcelona, Spain</footnote>

<footnote n="41">⁴¹Physics Department, 2320 Chamberlin Hall, University of Wisconsin-Madison, 1150 University Avenue Madison, WI 53706-1390, USA</footnote>

<footnote n="42">⁴²Department of Astronomy, University of California, Berkeley, 501 Campbell Hall, Berkeley, CA 94720, USA</footnote>

<footnote n="43">⁴³Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="44">⁴⁴Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544, USA</footnote>

<footnote n="45">⁴⁵Department of Physics and Astronomy, Pevensey Building, University of Sussex, Brighton, BN1 9QH, UK</footnote>

<footnote n="46">⁴⁶Centro de Investigaciones Energéticas, Medioambientales y Tecnológicas (CIEMAT), Madrid, 28040, Spain</footnote>

<footnote n="47">⁴⁷Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA</footnote>

<footnote n="48">⁴⁸Excellence Cluster Origins, Boltzmannstr. 2, D-85748 Garching, Germany</footnote>

<footnote n="49">⁴⁹Max Planck Institute for Extraterrestrial Physics, Giessenbachstrasse, D-85748 Garching, Germany</footnote>

<footnote n="50">⁵⁰Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians Universität München, Scheinerstr. 1, D-81679 München, Germany</footnote>

<footnote>This paper has been typeset from a LaTeX file prepared by the author.</footnote>