© The Authors 2025

1. Introduction

The discovery of the accelerating expansion of the universe (Riess et al. 1998; Perlmutter et al. 1999) was undertaken using standardised type Ia supernovae (SNe Ia). The cause of this accelerating expansion remains an unsolved mystery, but the standard model of cosmology associates it to a cosmological constant Λ; or, more generically, to dark energy. The two decades between the discovery of the accelerating expansion and have now seen sample sizes of SNe Ia grow from tens of supernovae to thousands. Modern measurements of the dark energy equation-of-state parameter, w, have statistical uncertainties on the order of ∼0.02 (Brout et al. 2022; Vincenzi et al. 2024; Popovic et al. 2024) and this field of research will soon be dominated by the relative systematic uncertainties. The decreasing statistical uncertainties place a greater emphasis on understanding the systematic uncertainties that are inherent with respect to measurements of cosmology undertaken using SNe Ia.

Here, we further improve on previous works that study the correlation between the properties of SNe Ia and their host galaxies (Sullivan et al. 2010; Rigault et al. 2013; Uddin et al. 2017; Rigault et al. 2020; Popovic et al. 2021). Most cosmology analyses with SNe Ia use a variation of the SALT (Guy et al. 2010) framework to fit type Ia supernova (SNIa) light curves and standardise their brightness. SALT uses two parameters to standardise the SNIa brightness: a stretch parameter, x₁, which encompasses the luminosity dependence on light curve duration and a colour parameter, c, which describes the wavelength-dependent luminosity. Since the study of Kelly et al. (2010), Sullivan et al. (2010), these two parameters have been joined by an ad hoc correction, the ‘mass step’, that accounts for the otherwise unexplained observed correlation between the standardised brightness of the SNIa and the properties of their host galaxy.

Correlations between SNIa properties and their environment have long been studied in the context of SN cosmology. Early on, the correlation between the light curve stretch and the host parameters was demonstrated (e.g. Hamuy et al. 1996; Howell et al. 2007; Sullivan et al. 2010; Lampeitl et al. 2010). Recently, with increased statistical power, progress has been made into uncovering underlying reasons for the correlations between fitted SALT properties (e.g., x₁, c) and the properties of the host galaxies of the SNe Ia. While Popovic et al. (2021) found strong evidence that both x₁ and c change with the host galaxy mass, Nicolas et al. (2021) provided an analytical description of changing x₁ values with both host galaxy mass and the local specific star formation rate (LsSFR). However, this refined description was insufficient to explain the mass step. In the last decades, several explanations arose to find the origin of the mass-step, with the aim to correctly account for this SN astrophysical dependency while fitting cosmological parameters. Rigault et al. (2013, 2020) and Briday et al. (2022), following earlier works on SN rates (e.g. Sullivan et al. 2006) suggest that prompt and delayed SNe Ia have intrinsically different absolute magnitudes. In that model, the prompt SNe, associated with recent star formation, are fainter, and the mass-step originates from the fact that massive galaxies host few prompt SNe Ia (see also: e.g. Childress et al. 2013, 2014).

Recently, following on prior works (from e.g. Mandel et al. 2017) Brout & Scolnic (2021) provided an explanation that is disconnected from the prompt and delayed-rate model. They show that the amplitude of the mass-step depends on the SN colour and were able explain the variation (and, consequently, the mass-step) by allowing the colour-magnitude relation to change as a function of host-mass. This model, further developed in Popovic et al. (2023), Kelsey et al. (2023), and Wiseman et al. (2022), follows the assumption that the redder end of the SN Ia color distribution is caused by dust, whose properties varies as a function of host stellar mass. This dust model describes the data well (Popovic et al. 2023).

While both approaches (dust or age) have their benefits and their limitations with respect to unveiling the true origin of the observed astrophysical biases is of paramount importance for SN cosmology as the resulting modelling for deriving distances may impact the measurement of cosmological parameters. In this analysis, we use the second data release of the Zwicky Transient Facility (ZTF; Bellm et al. 2019; Graham et al. 2019; Masci et al. 2019) SNIa sample (ZTF Data Release Paper), with spectra from Blagorodnova et al. (2018), Rigault et al. (2019). This release (ZTF-Cosmo-DR2) contains ∼3000 cosmological-quality SNe Ia at z < 0.1. Combined with recent improvements in the reliability of non-spectroscopically confirmed SNIa samples (Vincenzi et al. 2024; Popovic et al. 2024), we present the largest study of the evolution of the SNIa colour and colour-magnitude relation with redshift.

This study is important for the future of studies of the equation of state, w. Combining constraints from SNe Ia with external probes provides exquisite measurements of the properties of our universe, as in Vincenzi et al. (2024) and DESI Collaboration (2024). Systematics that evolve with redshift are degenerate with different cosmologies. With recent results from DESI Collaboration (2024) that indicate a non-static w may be preferred, even-further increased control of systematic uncertainties will be needed for the next generation of measurements of cosmology with SNe Ia.

The overview of the paper is as follows: Section 2 provides an overview of the data: ZTF, SDSS, PS1, and DES, alongside quality cuts and light curve fitting parameters. Section 3 explains the methodology used to define and fit the sample to explore evolution of c with redshift and mass. Section 4 presents the results and any potential evolution of c with other SNIa parameters, while Sects. 5 and 6, respectively, present the discussion and conclusions.

2. Data

High- and low-redshift surveys require different cadences, sky-coverage, and overall survey strategies. In order to maximise the statistics and redshift range, samples of SNIa are often combined, as in Amanullah et al. (2010), Betoule et al. (2014), Brout et al. (2022), Popovic et al. (2024).

We make use of the SNIa samples from the Zwicky Transient Facility (ZTF; ZTF Data Release Paper), the Sloan Digital Sky Survey (SDSS; Frieman et al. 2008; Sako et al. 2018; Popovic et al. 2020), Pan-STARRS (PS1; Jones et al. 2018), and Dark Energy Survey 5-year (DES; Vincenzi et al. 2024). Table 1 provides a summary of the dataset used in this analysis, showing the magnitude limit of the telescope, the resulting redshift cut, and the number of SNe in the final sample. We include an additional conservative sample, with a stricter redshift cut, that is shown in parentheses. This is explained further in Sect. 3.

The second ZTF data release covers the years of 2018 to 2020. The survey overview is presented in ZTF Data Release Paper. Light curve data typical have a two-day cadence in g and r and a 5 day cadence in i. All SNe Ia have a spectroscopic classification and about half have a host-galaxy redshift with the typical precision of 10⁻⁴. For the other half, the data release relies on SN spectroscopic feature that has been shown to have an unbiased redshift precision of 10⁻³. Table 2 shows a summary of the ZTF sample compared to the largest historic low-redshift surveys of SNe Ia; while all cover a similar redshift range (z < 0.1), ZTF is the first all-sky survey with statistics an order of magnitude larger than previous samples.

The SDSS supernova programme was run for three observing seasons between 2005 to 2007. The survey overview is provided in Frieman et al. (2008). SDSS light curves have data in the ugriz bands with an average observing cadence of four days. We took light curves from Sako et al. (2018) and host galaxies and redshifts from Popovic et al. (2020). PS1 ran from 2009 to 2013, observing with the griz filters at a cadence of six observations every five days. The PS1 light curves were taken from Chambers et al. (2016), with host galaxy redshifts drawn from papers listed in Jones et al. (2018).

The DES-SN programme was run for five years, covering 23 deg² across the sky. With an average observing cadence of even days, DES took data in the griz bands across eight ’shallow’ fields and two ‘deep’ fields, covering a balance between increased maximum redshift in the deep fields, which utilise multiple exposures and an area coverage and volume in the shallow fields; the data come from Sánchez et al. (2024).

SNe Ia at peak brightness are not uniform enough to support accurate measurements of cosmological parameters on their own; they require standardisation via a light curve modelling and fitting programme. We use the SALT2 programme introduced by Guy et al. (2010) with updated training from Taylor et al. (2021). The SALT fit returns four parameters for each SNIa: t₀ the time of peak brightness; x₀ the overall light curve amplitude; c the colour parameter; and x₁, the stretch parameter related to the width of the light curve. These latter three parameters are directly used to standardise SNIa luminosity, μ, according to the Tripp estimator (Tripp 1998):

$$\begin{array}{c}\hfill \mu =m_B+\alpha x_1-\beta c-M_0,\end{array}$$(1)

where m_B = −2.5log₁₀(x₀), c and x₁ are defined above, and M₀ is the absolute magnitude of a SNIa, with c = x₁ = 0. Then, α and β are the stretch-luminosity and colour-luminosity coefficients, respectively, and are defined for a given sample of SNIa.

For the data, we apply the following conventional cuts, following papers, such as Betoule et al. (2014), Scolnic et al. (2018), Brout et al. (2022):

• σ_x1 < 1: SALT2 x₁ uncertainty < 1,

• σ_c < 0.1: SALT2 c uncertainty < 0.1,

• σ_t0 < 1: Uncertainty on fitted peak brightness epoch < 1 days,

• −4 < x₁ < 4,

• −0.3 < c < 0.8,

• T_rest,min < 5: requires at least 1 observation five days prior to peak brightness (rest frame),

• T_rest,max > 0: Requires at least 1 observation after peak brightness (rest frame).

With respect to the case for SDSS, PS1, and DES, we institute an additional cut, that the probability of a SN being a SNIa is greater than 0.9: P_Ia > 0.9. We make use of the SuperNNova classifier from Möller & de Boissière (2020), taking trained surfaces from Popovic et al. (2024) and Vincenzi et al. (2024). This cut is not applied for ZTF, which has a spectroscopic identification for the entire sample. The extra cut placed on the high redshift surveys ensures that the sample is comprised of SNe Ia and that the limiting magnitude of these surveys is determined by the survey photometry rather than the spectroscopic followup, bringing them in line with the ZTF sample.

3. Methodology

This work has been made possible thanks to the unprecedented redshift range and statistics of the ZTF sample. Combining ZTF with the higher-redshift surveys, SDSS, PS1, and DES, offers a volume-limited sample of ∼3500 SNe Ia, ranging from 0.04 < z < 0.36. This selection is larger than any previously assembled collection of SNe Ia.

Here, we outline the methodology used to assess the evolution of SNIa colour with redshift and host galaxy mass. In this section, we begin with the definition and creation of a volume-limited sample for each survey so as to remove potential biases arising from selection effects, then we split the samples into redshift and mass bins to fit the colour population within each bin.

One consequence of the need for standardisation for SNIa is that redder SNIa (those approaching c = 0.3) are dimmer than bluer supernovae. This presents a problem for investigating the colour evolution with redshift; as the redshift increases, surveys approach the magnitude limit of being able to observe objects. This Malmquist bias can impact inferred cosmologies and, indeed, our colour distributions as well.

There are two solutions to this magnitude limit. The first is the introduction of bias corrections from Kessler & Scolnic (2017), using realistic simulations to correct for limitations in observing the apparent magnitudes of SNIa. The second method is to create a ‘volume-limited’ sample. Since the purpose of this work is to find the underlying distribution of colour, we are not able to provide accurate bias corrections. Instead, we follow Nicolas et al. (2021) and define two volume limited samples for each survey: a fiducial survey and a conservative one.

To define the volume limited sample, we follow Nicolas et al. (2021). We use Equation 1 and assume that the absolute magnitude for an SNIa with c = x₁ = 0 is M₀ − 19.36. From this, we can calculate a brightness at time of maximum light M = M₀ − αx₁ + βc. Assuming a common β = 3.1 and α = 0.15 (Betoule et al. 2014; Scolnic et al. 2018; Brout et al. 2022), the faintest magnitude we would expect to see is approximately −18.31 mags, from a SNIa with x₁ = −1.65 and c = 0.25. However, we require at least one observation before the peak brightness, which requires an even fainter magnitude, namely, at −18 for our fiducial limiting magnitude. Our conservative limiting magnitude is fainter, set at −17.5. In both cases, the limiting magnitudes we place are in reference to the SALT-fitted brightness.

With the limiting magnitudes and the magnitude limit of the surveys (taken from Nicolas et al. 2021, ZTF Data Release Paper), we were able to calculate a redshift, z_lim, such that we expect the sample to be complete and free of magnitude-based selection effects below z_lim. A review of the magnitude limits, limiting redshifts, and resulting number of SNIa is shown for each survey in Table 1. The redshift distributions for the full, fiducial, and conservative are shown in Fig. 1 for each survey included in this paper. We placed our redshift cuts on the sample after calculating the limiting redshift, z_lim, from the SALT-fitted magnitude of the supernovae.

Fig. 1. Redshift distributions for ZTF, SDSS, PS1, and DES. The full sample is presented for each survey before any redshift cuts in unfilled histogram. The light fill represents the redshift range of the fiducial sample and the dark fill shows the redshift range of the conservative sample. For details, see Sect. 2.

To investigate evolution of SNIa colour c with redshift, we quantile-bin our sample to create seven evenly populated bins of redshift. For the colour evolution with host galaxy stellar masses measured relative to the mass of our Sun (i.e. Sol = 1 M_⊙, M_Star = 10^X M_⊙; whereas here we report X), we made use of evenly spaced mass bins, so as to maintain a bin edge at M_Star = 10 and to be consistent with literature involving the mass step at 10M_Star.

In both cases, we fit within the individual bins of mass or redshift separately to find the apparent colour distribution. The change of the fit results in each bin can be assessed to determine any kind of evolution.

Within the individual bins, we made use of the fitting code from Ginolin et al. (2025) to fit the apparent colour via a convolution of a Gaussian intrinsic SNIa colour, c_int, and an exponential dust component, E_dust:

$$\begin{array}{c}\hfill c_{\mathrm{obs}}=c_{\mathrm{int}}+E_{\mathrm{dust}}+{ϵ}_{\mathrm{noise}},\end{array}$$(2)

with a noise term: ϵ_noise. The Gaussian and exponential tail is drawn as:

$$\begin{array}{c}\hfill P(c)=\mathcal{N}(c|\overline{c},c_{\sigma })\ast \{\begin{array}{cc}0\hfill & \text{if}c\le 0,\hfill \\ \frac{1}{E_{\mathrm{dust}}}{e}^{-c/E_{\mathrm{dust}}}\hfill & \text{if}c>0.\hfill \end{array}\end{array}$$(3)

The minimisation process is done via MINUIT (James & Roos 1975), which also returns fitting errors that are used within the analysis. The errors recorded by MINUIT are symmetric Gaussians. We report these results with $\overline{c}$ as the mean of the Gaussian intrinsic colour distribution, c_σ as the standard deviation, and E_Dust as the τ value describing the exponential distribution.

The choice here of an intrinsic Gaussian distribution that is reddened by external dust is taken from Brout & Scolnic (2021), but has recently gained traction with SNIa cosmology, used in such works as Wiseman et al. (2022), Kelsey et al. (2023), Popovic et al. (2021). It was also applied as the default model in the cosmology analyses performed by Brout et al. (2022), Vincenzi et al. (2024), Popovic et al. (2024).

Brout & Scolnic (2021) and Popovic et al. (2023) showed that while the observed SNIa colour distribution is an acceptable probe of the columnar density of dust E_Dust (e.g. the quantity of dust through the line-of-sight), it is not able to provide strong constraints on the R_v parameter. The best indicator of the R_v of the host galaxy, from SNIa data alone, is the colour-luminosity coefficient β. To further investigate the relationship between colour and host-galaxy properties, we investigate the relationship between host galaxy mass and the β parameter.

To fit β, we modify Equation (1):

$$\begin{array}{c}\hfill {\mu }^{\ast }=m_B-{\mu }_{\mathrm{theory}}+\alpha x_1\end{array}$$(4)

where m_B and x₁ are as normally defined, and we float α alongside our β. μ_theory is the theoretical distance modulus from a Flat ΛCDM cosmology with H₀ = 74. The slope of this μ^* value, when plotted against the supernova colour c, is β. We fit this slope, simultaneously with the intercept, using MINUIT in bins of mass to determine potential β evolution.

The masses reported in this paper have been taken from their respective data releases, namely: ZTF DR2 Data Release and Vincenzi et al. (2024) for ZTF and DES, or in the case of SDSS and PS1, the updated masses from Popovic et al. (2024). All surveys use consistent PEGASE.2 code with a Kroupa (2001) initial mass function. We used the same methodology and star formation history as in Smith et al. (2020); in short, the SED of each star formation history was initialised with a metallicity of 0.004, which evolves consistently, calculated at 102 time steps from 0 to 14 Gyr. We used seven foreground dust screens, with the colour excess ranging from 0 to 0.3 magnitudes.

4. Results

We present the results of the fitting process from Sect. 3 on the samples defined in Sect. 2. In Fig. 2 shows m_B vs. redshift for the fiducial sample, alongside the demarcations for the quantile redshift bins. In the inset, we show the colour distributions of each of the redshift bins; here we investigate the cause of the apparent change in colour distribution. Significance measurements are calculated by comparing the binned results to the mean of the overall distribution and converting them to a significance measurement using the number of degrees of freedom.

Fig. 2. mB vs. redshift for the SNIa in the Fiducial Sample defined in Table 1. Each supernova is colour coded according to its colour c. The grey vertical lines demark the quantile redshift bins, and the colour histogram is provided for each redshift bin in the inset plot. See Sect. 2 for more details.

4.1. Evolution with redshift

We fit the colour distribution of the SNIa in each of the seven quantile redshift bins and present the inferred Gaussian+exponential distribution in Fig. 3.

Fig. 3. Fit results for the quantile redshift bins of the fiducial sample. The c distribution is shown for each redshift bin and colour coded to become darker with increasing redshift. For details, see Sect. 4.1.

We plotted the Gaussian and exponential parameters from Fig. 3 to more accurately assess any potential trends. We find that the Gaussian parameters that describe the SNIa colour distribution are consistent with no evolution. $\overline{c}$ changes with redshift at a confidence of < 1σ, and c_σ changes with redshift at a confidence of 3.4σ. These low signals are in contrast with the E_Dust parameter, for which we find evidence of evolution with redshift at a confidence of > 6σ. For comparison, we show the χ² values from the conservative sample that was detailed in Sect. 3. The conservative sample trends are qualitatively similar to the full sample, though the larger uncertainties that arise from the smaller sample decrease the evidence of E_Dust evolving from 6σ to 3σ (Fig. 4).

Fig. 4. Fit parameters for the Gaussian and exponential distributions plotted with increasing redshift. Alongside each parameter, the mean fit value is plotted in dash-dotted line to provide context. The fit results from the conservative sample are overlaid in triangles for comparison. See more in Sect. 4.1.

Table 3 shows the χ² breakdown for each of the parameters and their evolution with redshift. We see the trends from the fiducial sample (i.e. data that appear to prefer an evolving dust-based E_Dust parameter over a change in the intrinsic SNIa colour distribution) hold for the conservative sample, although we do find a confidence of > 5σ for E_Dust evolution in the conservative sample, instead of > 6σ.

4.2. Evolution with mass

Figure 5 shows the evolution of the host galaxy mass with redshift. There is an under-density of SNe between z = 0.06 and z = 0.1; this under-density arises from the redshift cut on ZTF at z = 0.06, leaving only SDSS SNe in the range of 0.06 < z0.1. Even with this under-density we find no evidence of the median host galaxy mass changing with redshift (< 1σ from the median of the full distribution).

Fig. 5. Host galaxy mass plotted vs redshift for the Fiducial sample. Supernovae are shown in grey, the coloured circles represent the median host galaxy mass in each quantile redshift bin. The dash-dotted black line shows the median host galaxy mass of the entire sample. The mass distribution, marginalised over redshift, is shown on the right. See Sect. 4.2 for more details.

Bolstered by the confidence that there is no correlation between the host galaxy mass and redshift, we show the evolution of the SNIa colour distribution with host galaxy mass. Figure 6 shows the change of colour distribution with host galaxy mass, binned in steps of 0.5M_Star, starting at 8M_Star and going to 11.5M_Star. As a reminder, in contrast to the quantile binning for redshift, we do not require equal statistics in each bin, so that we can place a bin edge at the historical mass step value of 10M_Star.

Fig. 6. Fit results for the host galaxy mass bins of the fiducial sample. The c distribution is shown for each mass bin, colour coded with increasing host galaxy mass. (Sect. 4.2).

Figure 7 displays the individual parameters $\overline{c}$, c_σ, and E_Dust from the fitted results presented in Fig. 6. We show the χ² values in Table 4 and we see results similar to the redshift distribution results. There is little evidence of the SNIa colour parameters changing, as $\overline{c}$ exhibits a change of < 1σ and c_σ exhibits a 3.3σ signal of evolution with host galaxy mass. Again, the conservative sample, which is plotted alongside the fiducial, does not show significantly different trends, shown in Table 4.

Fig. 7. Fit parameters for the Gaussian and exponential distributions plotted with increasing host galaxy mass. Alongside each parameter, the mean fit values is plotted in dash-dotted line to provide context. The fit results from the conservative sample are overlaid with triangles for comparison. See Sect. 4.2 for more details.

On the other hand, the E_Dust parameter displays a clear evolution, namely, > 6σ, with the mass of the host galaxy. The E_Dust values peak around the host galaxy mass of 10 /M_Star, decreasing beyond this point for both higher and lower mass galaxies. The first mass bin, however, does not follow this trend. This does not appear to be a binning effect, as it persists across different bin edges and sizes. However, the discrepancy here may be caused by the mis-association of host galaxies, or biases arising from too-faint hosts (Paulino-Afonso et al. 2022). This trend is equally strong with the conservative sample, at > 6σ.

For both the z and host galaxy mass evolution tests, we attempted to account for potential non-Ia contamination in the signal by making use of simulations from Popovic et al. (2024). While these simulations are only available for SDSS and PS1, they are a good test for any potential signal contaminant by representing the two surveys with the highest potential non-Ia contamination. We performed the redshift and host galaxy mass analysis twice: once with a simulated contamination and once with a pure-Ia sample. We applied the same P_Ia > 0.9 cut to the contaminated simulated sample as we did in the case of the real data to test any bias due to this additional cut. This is compared with the pure SNIa simulated sample. We find no significant difference ( < 0.5σ) in recovered results between the contaminated and un-contaminated samples.

4.3. Evolution of β with host galaxy mass

While the observed c distribution is a tracer of E_Dust, it does not provide strong constraints on the type or size of the dust particulates, conventionally denoted as R_v. Instead, R_v has been shown to be better tracked with supernova data by the colour-luminosity relationship, β (Brout & Scolnic 2021; Popovic et al. 2023; Johansson et al. 2021). Here, we investigate the evolution of β with the host galaxy mass to see if whether a change is present in the dust properties, alongside the observed changes in the columnar density.

Figure 8 shows the colour-luminosity relationship for seven mass bins across the fiducial sample. Again, we have fit α alongside our β and assumed a Flat ΛCDM cosmology in making these plots. We focus on the inset that shows the fitted β values for each mass bin. With the exception of the first bin, there exists a clear trend of decreasing β with increasing host galaxy mass. Considering the bins that are populated for all three redshift categories (‘full’, ‘low-z’, and ‘high-z’), we find evidence of a change in β with mass at the 3σ level, likely driven by the lowest mass bins. These results are in line with those of Popovic et al. (2023), Chen et al. (2022) and provide further hints of changing host-galaxy properties. We discuss the inclusion of the first bin in Sect. 5.

Fig. 8. Comparison of c vs. mB − μ + αx1, offset in bins of host galaxy mass, shown in the main plot. The SNe are colour-coded by their colour, c, for visual guidance. The slope (β value) for each mass bin is shown in red, alongside the fiducial beta in black dash-dotted line. The light grey dash-dotted vertical line shows the c = 0 point for each mass bin. Top: Colour histograms for each of the mass bins is shown, colour-coded according to c. Inset: β vs mass relationship for the fitted mass bins. The β values for the full sample are presented in brown points, the ZTF-only β values in teal triangle, and the high-redshift sample in olive triangle, with the fiducial β shown in grey dash-dotted line. See Sect. 4.3 for more details.

5. Discussion

5.1. Redshift and mass evolution

We find no strong evidence for an evolving dust-free SNIa colour distribution with redshift. In our fiducial case, there are hints of the standard deviation of the dust-free distribution changing (3.4σ), but < 1σ evidence for the change in the mean of the dust-free distribution. Instead, we find very strong evidence of an evolving quantity of dust with redshift (> 6σ). Our conservative case finds no evidence for dust-free evolution and a weaker trend in dust evolution (3σ), although the weaker dust evolution trend is due to larger statistical uncertainties. However, the question of E_Dust evolution with the host galaxy mass becomes trickier. We have taken into consideration the exclusion of the first bin in our mass results; this is motivated by the potential of host galaxy mis-association at lower redshift. Our lowest mass bin at M < 8M_Star is primarily dominated by low redshift galaxies; those galaxies at high redshift come with significant errors (on or above the order of 1M_Star). If a SNIa was incorrectly associated to a lower mass host, it would effectively raise the observed E_Dust as a result of contamination from the true higher-mass hosts.

Rather than a linear model across the host galaxy stellar mass for E_Dust and β, it is possible that the 10¹⁰ mass bin represents a break-point for these parameters. For stellar masses above 10¹⁰, E_Dust is consistent with no evolution within 1.7σ; below this mass, the E_Dust parameter evolves with mass at a significance of 3.8σ. The opposite occurs with β: below a mass of 10¹⁰β does not change (< 1σ change from a flat line), compared to the familiar 3σ above 10¹⁰. In either case, these results hint at mass-dependent dust properties of SNe Ia.

This marks the first time that dust information with this fidelity has been extracted from SNIa light curves alone. With the exception of the aforementioned first bin, we find a continuous increase in the columnar density, E_Dust, to a peak and inflection point around 10M_Star, followed by a decrease. This trend has been found in the literature, notably appearing in Zahid et al. (2013) and Calura et al. (2017). Zahid et al. (2013) finds an increasing extinction, A_V, until an inflection point of 10M_Star, whereby it decreases again. Interestingly, this behaviour holds for all star formation rates. Calura et al. (2017) shows increasing dust quantities with host galaxy mass, but only for large elliptical galaxies; spiral galaxies exhibit a flat relationship between mass and dust density. Rigault et al. (2020) proposes that this unique shape of the dust property evolution with mass is the byproduct of competing production modes: star formation rate (SFR) driven processes are responsible for increasing the amount of dust up until ≈10¹⁰ stellar masses, when star formation stops. Subsequently, explosion-processes driven by the large galactic stellar mass are the primary driver of dust relations for those galaxies heavier than ≈10¹⁰ stellar masses.

5.2. Evolution of β with redshift

The colour-luminosity coefficient β exhibits a strong linear relationship with the host galaxy stellar mass and this relationship holds regardless of redshift range. This follows the trends seen in Brout & Scolnic (2021) and Popovic et al. (2023), although those papers only fit in two bins of mass due to computational and statistical constraints, respectively. Previous analyses, such as those of Sullivan et al. (2010) and González-Gaitán et al. (2021), also reported evidence of multiple β values, but there was no evidence to support the notion that multiple β values alone could explain the ‘mass step’. Unlike González-Gaitán et al. (2021), we did not find any evidence of a sharp change in β at 10M_Star. Instead, a continuous β relationship as a function of the host galaxy, R_v, remains the most probable explanation for this behaviour.

As a test to ensure that this behaviour of changing β with mass is not a sampling issue, we re-fit several downsampled instances of the full data set to check potential biases on β. We did not find significant biases on β (< 2σ) for samples as small as ×20 smaller than our full data.

In addition, we tested the evolution of β with redshift, defining

$$\begin{array}{c}\hfill \beta (z)={\beta }_0+{\beta }_1\times z,\end{array}$$(5)

where β₀ is the familiar SALT β value and β₁ is the redshift-dependent term. We find β₀ = 3.25, but no strong evidence for β evolution: β₁ = 0.3 ± 0.4.

We note that earlier eras of cosmological analysis, such as JLA (Betoule et al. 2014), were dominated by systematic uncertainties arising from calibration. Today, analyses such as Popovic et al. (2024), Vincenzi et al. (2024) find that astrophysical processes not captured by light curve fitting methods are the driving force of systematic uncertainties. The results of this paper show strong evidence that the density of dust, either around the site of SNIa explosions, or within the galaxy as a whole, changes with redshift. This effect is degenerate with cosmology. The improper treatment of these effects in any cosmological analysis will preclude an accurate assessment of cosmological parameters.

6. Conclusion

In this paper, we have assembled a sample of ≈3000 SNIas that are free from magnitude-based selection effects, compiled from ZTF, SDSS, PS1, and DES. This selection marks the largest compilation of SNIa light curves to date, demonstrating the power of the ZTF sample.

With our ‘volume-limited’ sample and its more conservative counterpart, we have shown strong evidence for evolving host-galaxy properties. The amount of columnar dust density, E_Dust, changes significantly with both redshift z and the host galaxy mass. This contrasts with the lack of strong evidence for a change in the dust-free SNIa colour distribution, which is found to be well characterised by a Gaussian distribution in this work.

These results stand independent of potential non-Ia contamination, which we find to be mitigated by the placement of a Ia-probability cut on the data and our results persist for the conservative sample, providing strong evidence that these results are not a byproduct of magnitude-limited selection effects.

Additionally, we found evidence of a continuous and linear relationship between the colour-luminosity coefficient β and the host galaxy stellar mass, M_Star. Taken altogether, these results point towards further proof that future standardisation modes for SNIa light curves must account for the properties of the host galaxy and that SNe Ia are becoming useful tools for extracting the dust properties of their host galaxies.

Acknowledgments

Based on observations obtained with the Samuel Oschin Telescope 48-inch and the 60-inch Telescope at the Palomar Observatory as part of the Zwicky Transient Facility project. ZTF is supported by the National Science Foundation under Grants No. AST-1440341 and AST-2034437 and a collaboration including current partners Caltech, IPAC, the Weizmann Institute of Science, the Oskar Klein Center at Stockholm University, the University of Maryland, Deutsches Elektronen-Synchrotron and Humboldt University, the TANGO Consortium of Taiwan, the University of Wisconsin at Milwaukee, Trinity College Dublin, Lawrence Livermore National Laboratories, IN2P3, University of Warwick, Ruhr University Bochum, Northwestern University and former partners the University of Washington, Los Alamos National Laboratories, and Lawrence Berkeley National Laboratories. Operations are conducted by COO, IPAC, and UW. SED Machine is based upon work supported by the National Science Foundation under Grant No. 1106171 The ZTF forced-photometry service was funded under the Heising-Simons Foundation grant #12540303 (PI: Graham). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement n°759194 – USNAC). This work has been supported by the Agence Nationale de la Recherche of the French government through the programme ANR-21-CE31-0016-03. Y.-L.K. has received funding from the Science and Technology Facilities Council [grant number ST/V000713/1]. This work has been supported by the research project grant “Understanding the Dynamic Universe” funded by the Knut and Alice Wallenberg Foundation under Dnr KAW 2018.0067, Vetenskapsrådet, the Swedish Research Council, project 2020-03444. L.G., T.E.M.B acknowledges financial support from the Spanish Ministerio de Ciencia e Innovación (MCIN) and the Agencia Estatal de Investigación (AEI) 10.13039/501100011033 under the PID2020-115253GA-I00 HOSTFLOWS project, from Centro Superior de Investigaciones Científicas (CSIC) under the PIE project 20215AT016 and the programme Unidad de Excelencia María de Maeztu CEX2020-001058-M, and from the Departament de Recerca i Universitats de la Generalitat de Catalunya through the 2021-SGR-01270 grant. Y.-L.K. has received funding from the Science and Technology Facilities Council [grant number ST/V000713/1]. UB, GD, JHT are supported by the H2020 European Research Council grant no. 758638.

References

All Tables

All Figures

	Fig. 1. Redshift distributions for ZTF, SDSS, PS1, and DES. The full sample is presented for each survey before any redshift cuts in unfilled histogram. The light fill represents the redshift range of the fiducial sample and the dark fill shows the redshift range of the conservative sample. For details, see Sect. 2.
In the text

	Fig. 2. mB vs. redshift for the SNIa in the Fiducial Sample defined in Table 1. Each supernova is colour coded according to its colour c. The grey vertical lines demark the quantile redshift bins, and the colour histogram is provided for each redshift bin in the inset plot. See Sect. 2 for more details.
In the text

	Fig. 3. Fit results for the quantile redshift bins of the fiducial sample. The c distribution is shown for each redshift bin and colour coded to become darker with increasing redshift. For details, see Sect. 4.1.
In the text

	Fig. 4. Fit parameters for the Gaussian and exponential distributions plotted with increasing redshift. Alongside each parameter, the mean fit value is plotted in dash-dotted line to provide context. The fit results from the conservative sample are overlaid in triangles for comparison. See more in Sect. 4.1.
In the text

	Fig. 5. Host galaxy mass plotted vs redshift for the Fiducial sample. Supernovae are shown in grey, the coloured circles represent the median host galaxy mass in each quantile redshift bin. The dash-dotted black line shows the median host galaxy mass of the entire sample. The mass distribution, marginalised over redshift, is shown on the right. See Sect. 4.2 for more details.
In the text

	Fig. 6. Fit results for the host galaxy mass bins of the fiducial sample. The c distribution is shown for each mass bin, colour coded with increasing host galaxy mass. (Sect. 4.2).
In the text

	Fig. 7. Fit parameters for the Gaussian and exponential distributions plotted with increasing host galaxy mass. Alongside each parameter, the mean fit values is plotted in dash-dotted line to provide context. The fit results from the conservative sample are overlaid with triangles for comparison. See Sect. 4.2 for more details.
In the text

Fig. 8. Comparison of c vs. mB − μ + αx1, offset in bins of host galaxy mass, shown in the main plot. The SNe are colour-coded by their colour, c, for visual guidance. The slope (β value) for each mass bin is shown in red, alongside the fiducial beta in black dash-dotted line. The light grey dash-dotted vertical line shows the c = 0 point for each mass bin. Top: Colour histograms for each of the mass bins is shown, colour-coded according to c. Inset: β vs mass relationship for the fitted mass bins. The β values for the full sample are presented in brown points, the ZTF-only β values in teal triangle, and the high-redshift sample in olive triangle, with the fiducial β shown in grey dash-dotted line. See Sect. 4.3 for more details.

In the text